#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);
   }



   
main()

{

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

 long primes[50002];
 long sqrtpr[50002];   /* see if you can have larger sizes */


 long biglen;


 biglen =50001;   /*ESTABLISHES HOW HIGH WE SHALL GO*/

 /*!!!!!!!!!!!!!!!!!!!!!!!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 <= biglen; ++j)
   {
    primes[j] = 1;
    sqrtpr[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 */

 biglen = 50001;

 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;

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

 x = 21;
 for (j = 73; j <= biglen; j = j + 2)
   {
    j = nextprime(j,y);
    primes[x] = j;
    x = x + 1;
   }

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

 numprimesless = x - 1;


 /*!!!!!!!!!!!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 */



 sqrtpr[1] = 1;
 sqrtpr[2] = 1;
 sqrtpr[3] = 2;
 sqrtpr[4] = 2;
 sqrtpr[5] = 3;
 sqrtpr[6] = 3;
 sqrtpr[7] = 4;
 sqrtpr[8] = 4;
 sqrtpr[9] = 4;
 sqrtpr[10]= 4;
 sqrtpr[11]= 5;
 sqrtpr[12]= 5;
 sqrtpr[13]= 6;
 sqrtpr[14]= 6;
 sqrtpr[15]= 6;
 sqrtpr[16]= 6;
 sqrtpr[17]= 7;
 sqrtpr[18]= 7;
 sqrtpr[19]= 8;
 sqrtpr[20]= 8;

 x = 8;
 for (j = 21; j <= biglen; j=j+2)
  {
   y = 0;
   y = primecheck(j,y);
   if (y == 1) { x = x + 1; }
   sqrtpr[j]  = x;
   sqrtpr[j+1]= x;
  }


 /*!!!!!!!!!!e x p l a n a t i o n   o f   a b o v e !!!!!!!!!!!
    this is a little more tricky than the process used to store the
    primes; in checking the primality of a number N, it is sufficient
    to check wrt all the primes <= Sqrt(N); for instance, if N = 400,
    then Sqrt(N) = 20; all primes <= 20 are {2,3,5,7,11,13,17,19}, so
    it is sufficient to check N mod j for j in {2,3,5,...,19}, although
    for convienience, we will choose N odd st it is not necc to check
    N mod 2.

    explanation:
    for (j=21,...,+2)  => this moves us forward two units at a time
			  we will then determine what value goes in
			  sqrtpr[j]; this integer will be the number of
			  primes we need to check (again, ignore 2)
    y = 0              => initialization step
    y = primecheck(j,y)=> if j is prime, y = 1 else y = 0
    if (y==1) {++x}    => if y = 1 this is because j was prime; we now
			  have a new prime factor <= sqrt(N), so we
			  must increase the number of prime factors to
			  check by one, hence x -> x+1;
			  if y = 0, j was not prime, hence the only
			  possible prime factors are still just those
			  x ones we had before, and no need to increase
			  x.
    sqrtpr[j]  = x     => for a number N, we will find the integer
			  j = [Sqrt(N)], sqrtpr[j] now tells us how
			  many primes we must check.
    sqrtpr[j+1]= x     => NOTE: care must be taken with the even values
			  we note, however, that the number of primes
			  <= 2m is the same as those <= 2m - 1  */





 /*!!!!!!!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:      " P R I M E 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!!!!!!!!!*/



   fp = fopen("primedat.fil", "w");    /* opens for writing */
   fprintf(fp, "%ld ", numprimesless); /* writes numprimesless*/
   for (j = 1; j <= numprimesless; ++j)
      {
	fprintf(fp, "%ld ", primes[j]);
      }
   printf("\n");
   fprintf(fp, "%ld ", biglen);
   for (j = 1; j <= biglen; ++j)
      {
	fprintf(fp, "%ld ", sqrtpr[j]);
      }

   fclose(fp);                        /*closes file */



}