#include <stdlib.h>
#include <stdio.h>
#include <math.h>


primecheck(long pr, long check)
   {
    long m,mm,m1,m2;

    /* works as follows: this tests to see if pr is prime.
       we only need to test primes < sqrt(pr). we won't be
       efficient, as this runs so quickly as pr < 10^5 --
       we are just using this to initialize the array of
       startup primes -- if we wanted to check primality at
       larger numbers, we'd use a different algorithm */

    /* method as follows: see if our number is divisible by
       3. If yes, stop, else check divisibility by numbers of
       the form 6n + 1, 6n - 1. The default is check = 1, the
       number is prime. If it is divisible by an odd number <
       its sqrt, we change check to 0.  */

    check = 1;                          /* check = 1 prime, =0 not */
    mm =  (long) floor( sqrt( (double) pr) );  /* largest int < Sqrt(pr)  */
    mm =  (long) floor( mm / 6 );
    if (pr % 3 == 0) check = 0;
    else
      for (m = 1; m <= mm + 1; ++m)
       {
       	m1 = 6*m - 1;
       	m2 = 6*m + 1;
       	if (pr % m1 == 0) { check = 0; m = mm + 1; }
       	else
       	if (pr % m2 == 0) { check = 0; m = mm + 1; }
       }

    return(check);
   }


   /*BIG WARNING: THE ABOVE PROCEDURE, THOUGH ADEQUATE FOR OUR PURPOSES,
     HAS A FLAW: will not recognize primes < 10; this has to do with the
     algorithm used--we take the Sqrt(of our prime); we then divide by 6
     and check the odd factors <= Sqrt(of our prime) mod 6 + 1;
     the problem is for small primes, we check pr mod pr, which == 0;
     while this can be remedied, since we will not be using this for small
     primes, we shall leave this as is in the sake of expediency
   */



nextprime(long pr, long check)
   {

   /* this algorithm calculates the next prime after pr */

    long m,mm;
    check = 0;
    while( check == 0)
      {
       check = primecheck(pr, check);
       pr = pr + 2;
      }
    return(pr - 2);
   }



   
void main()

{

 FILE *fp;
 long numprimesless,p,pp,x,y,z,j,k,l;
 int  xint, yint;

 long primes[450];

 long counter, biglen;


 biglen = 450001;   /*ESTABLISHES HOW HIGH WE SHALL GO*/
 biglen = 10000000;           

 /*!!!!!!!!!!!!!!!!!!!!!!!I m p o r t a n t!!!!!!!!!!!!!!!!!!!
	we shall use biglen to represent the size of the array used
	to store the prime values; as the memory capacities of the
	computers at WNSL and at my house differ by orders of
	magnitude, this number will be different at each place */
 /* Wright Nuclear Structures Laboratory */
   
 for (j = 0; j <= 400; ++j)
   {
    primes[j] = 1;
   }

   /* the above initializes every number to be prime. We'll save
      the information in primes[j]. sqrtpr[j] will be a very
      convenient array: for a number n, if you want to see if it's
      prime you just have to check up to the largest prime < sqrt(n).
      sqrtpr tells you where that last prime is */

 primes[1] = 2;
 primes[2] = 3;
 primes[3] = 5;
 primes[4] = 7;
 primes[5] = 11;
 primes[6] = 13;
 primes[7] = 17;
 primes[8] = 19;
 primes[9] = 23;
 primes[10]= 29;
 primes[11]= 31;
 primes[12]= 37;
 primes[13]= 41;
 primes[14]= 43;
 primes[15]= 47;
 primes[16]= 53;
 primes[17]= 59;
 primes[18]= 61;
 primes[19]= 67;
 primes[20]= 71;

 fp = fopen("eprdat.fil", "w");    /* opens for writing */
 for (j = 1; j <= 20; ++j)
   {
     fprintf(fp, "%ld ", primes[j]);
   }

 /* initializes the first few primes, as our algorithm
    is not designed to handle small primes */

 x = 21; counter = 20; j = 73;
 while (counter <= biglen)
   {
    j = nextprime(j,y);
    fprintf(fp, "%ld ", j);
    x = x + 1;
    counter = counter + 1;
    j = j + 2;
   }

   /* calculates the primes up to biglen  */
   /* x is a counter for number of primes */

 /*!!!!!!!!!!!s u m m a r y   o f   a b o v e!!!!!!!!!!!!!!!
     first we input the first 20 prime numbers; the procedure/
	function nextprime(j,y) gives the first prime >= j.
	we put this value into the block primes[x], increment x, and
	we continue cycling through j's until we reach the next prime
	which is now stored in primes[x+1]; we continue this process
	until we have stored the first biglen primes
     numprimesless is, as the name implies, simply the number of primes
	less than or equal to biglen; we will need this number for
	writing and reading our information to the files */



 /*!!!!!!!F I L E   W R I T I N G   I N F O R M A T I O N!!!!!!!!!!!
    (1) first we write numprimesless to our file; this will be needed
	  when we recall the list of primes <= biglen
    (2) next we will write the actual list of primes <=biglen
    (3) then, for consistency, we will write biglen
    (4) we then write the sqrtpr[] numbers


     writing to:      " E P R D A T . F I L "




   !!!!!!F I L E    W R I T I N G   I N F O R M A T I O N!!!!!!!!!*/



   printf("\n");
   fclose(fp);                        /*closes file */



}