/* primegen.0.4.c
 * j. vehent - july 2009
 * 
 * this code test potential prime number with Rabin-Miller
 * and if it finds one, it tries to divide it by every single 
 * prime number before it (from 3 to square root of number)
*/

#include <stdio.h>
#include <gmp.h>

#define MAXPRIMELENGTH 1048576
#define PRIMELISTFILE "prime_list.txt"
#define TRICKEDLIST "tricked_list.txt"

int main()
{
	/* init a big number to value "3"
	*/
	mpz_t bign;
	mpz_init_set_str(bign,"3",10);

	mpz_t rootsquare_of_bign;
	mpz_init_set_str(rootsquare_of_bign,"1",10);

	int rootsquare_reached = 0;

	mpz_t two;
	mpz_init_set_str(two,"2",10);

	int isnotprime = 0;

	FILE *fd;
   	FILE *fdtricked;
	char p[MAXPRIMELENGTH];
	mpz_t bp;

	// if file already exist, open it and get the last value
	if(NULL!=(fd=(fopen(PRIMELISTFILE,"r"))))
	{
		//go to the end of the file
		while(fgets(p,MAXPRIMELENGTH,fd));
		// init bign at the last value of p
		mpz_init_set_str(bign,p,0);
	}
	else
	{
		// if it doesn't, create it
		if(NULL==(fd=(fopen(PRIMELISTFILE,"w"))))
		{
			fprintf(stderr, "UNABLE TO CREATE %s\n", PRIMELISTFILE);
			return -98;
		}
		fprintf(fd,"2\n");
	}
	fclose(fd);

	/* 2 tests are processed on each numbers :
		1. rabin miller primality test
		2. div with all the other prime from 2 to square root of the number
	*/
	while(1)
	{
		//go to next odd number
		mpz_add(bign,bign,two);

		// rabin miller and some trivial div tests
		if(mpz_probab_prime_p(bign, 10) != 0)
		{

			mpz_sqrt(rootsquare_of_bign, bign);
			rootsquare_reached = 0;
			isnotprime = 0;

		         // div bign with all previous primes
			if(NULL == (fd = fopen(PRIMELISTFILE, "r")))
			{
				fprintf(stderr, "ERROR OPENING %s !!!!",PRIMELISTFILE);
				return -99;
			}

			while(fgets(p, MAXPRIMELENGTH, fd) && rootsquare_reached == 0 && isnotprime == 0);
			{
				// convert p from char to bigint
				if (0 != (mpz_init_set_str(bp, p, 10)))
				{
					fprintf(stderr, "ERROR READING NUMBER FROM %s\n",PRIMELISTFILE);
					return -97;
				}

				// if prime from list is higher than square root of bign, stop here
				if (mpz_cmp(bp, rootsquare_of_bign) > 0)
				{
					rootsquare_reached = 1;
				}
				else
				{
					if(mpz_divisible_p(bign,bp) != 0)
						isnotprime = 1;
				}
			};
			fclose(fd);
			// add new prime at the end of the list
			if(isnotprime == 0)
			{
				// print new prime to output file
				if (NULL == (fd = fopen(PRIMELISTFILE,"a+")))
				{
					fprintf(stderr,"ERROR OPENING %s !!!!",PRIMELISTFILE);
					return -99;
				}
				mpz_out_str(fd,10,bign);
				fprintf(fd,"\n");
				fclose(fd);
			}

        	 	// keep prime number that tricked the rabin-miller test
         		else
		         {
				// print bign and bigp to the tricked list
				if(NULL==(fdtricked=(fopen(TRICKEDLIST,"a+"))))
				{
	        	       		// if list doesn't exist, create it and write the number
        	       			if(NULL==(fdtricked=(fopen(TRICKEDLIST,"w"))))
               				{
	                	  		fprintf(stderr, "UNABLE TO CREATE %s\n", PRIMELISTFILE);
        	          			return -98;
	               			}
		            	}
        		    	mpz_out_str(fdtricked,10,bign);
            			fprintf(fdtricked," passed Rabin-Miller and is divisible by ");
		            	mpz_out_str(fdtricked,10,bp);
				fclose(fdtricked);
         		}
		}
	}
	return -1;
}