/* primegen
 * jve - mai 2008
 */
#include <stdio.h>
#include <gmp.h>

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

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

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

	int squarereached = 0;

	mpz_t two;
	mpz_init_set_str(two,"2",10);
	
	int isnotprime = 0;

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

	if(NULL==(fd=(fopen(PRIMELISTFILE,"w"))))
	{
		fprintf(stderr,"UNABLE TO CREATE %s\n",PRIMELISTFILE);
		return -98;
	}
	else
	{
		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(square_bign, bign);
			squarereached = 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) && squarereached == 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, square_bign) > 0)
				{
					squarereached = 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);
			}
		}
	}

	return -1;
}