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

int main(int argc, char *argv[]){
	unsigned long long p, i2pp1, k=1, ik2pp1;
	int isprime;
	int composite=0, cantdidate=0;

	mpz_t prime, pprime;
	mpz_t one, z2exp, candidate;
	mpz_t z2p, z2pp1, k2p, k2pp1, gcd, sqrt, temp, m, tm;
	mpz_t a, q, r, final;
	mpz_init_set_ui(prime, 2);
	mpz_init_set_ui(one, 1);
	mpz_inits(pprime, z2exp, candidate, z2p, z2pp1, k2p, k2pp1, gcd, sqrt, temp, m, tm, a, q, r, final, NULL);
	
	clock_t total_start = clock();

	while(1) {
		clock_t start = clock();
		p = mpz_get_ui(prime);
		i2pp1 = 2*p+1 % 8;
		if(i2pp1 == 1 || i2pp1 == 7) {
			mpz_mul_2exp(z2p, prime, 1);
			mpz_set(k2p, z2p);
			mpz_add_ui(z2pp1, z2p, 1);
			isprime = mpz_probab_prime_p(z2pp1, 50);
			if(isprime==1) {printf("2pp1 only probably prime for p=%d\n", p);}
			if(isprime>0) {
				if(p%4 == 3) { /*printf("%d composite\n", p); */goto next; }
				mpz_mul_2exp(z2exp, one, p);
				mpz_sub_ui(candidate, z2exp, 1);
				if(mpz_divisible_p(candidate, z2pp1)!=0){  /*printf("%d composite\n", p);*/ goto next;}
			}
		}
		else {
			mpz_mul_2exp(z2p, prime, 1);
			mpz_set(k2p,z2p);
			mpz_add_ui(z2pp1,z2p,1);

			mpz_mul_2exp(z2exp, one, p);
			mpz_sub_ui(candidate, z2exp, 1);
		}
		mpz_set_ui(m,1);
		for(int i=0;i<p*(32-__builtin_clz(64-__builtin_clzll(p)))>>1;i++) {
			mpz_add(temp, k2p, z2p);
			mpz_swap(k2p, temp);
			mpz_add_ui(k2pp1, k2p, 1);
			ik2pp1 = mpz_get_ui(k2pp1);
			ik2pp1 = ik2pp1%8;
			if(ik2pp1==1 || ik2pp1==7) {
				// divexact tested and found to be slower.
				if(mpz_divisible_p(candidate, k2pp1)!=0){ clock_t end = clock() - start; /*printf("%d composite in %f\n", p, (double) end / CLOCKS_PER_SEC);*/ composite += 1; goto next;}
			}
			k+=1;
		}
		clock_t end = clock() - start;
		//printf("%d a candidate in %f\n", p, (double) end / CLOCKS_PER_SEC);
		// test
		mpz_set_ui(final, 9);
		for(unsigned long j=1; j<p; j++) {
			mpz_mul(a, final, final);
			mpz_tdiv_q_2exp(q, a, p);
			mpz_tdiv_r_2exp(r, a, p);
			mpz_add(final, q, r);
			mpz_tdiv_q_2exp(q, final, p);
			mpz_tdiv_r_2exp(r, final, p);
			mpz_add(final, q, r);
		}
		clock_t total = clock() - total_start;
		if(mpz_cmp_ui(final, 9)==0) {
			printf("%d prime in %f\n", p, (double) total / CLOCKS_PER_SEC);
			//printf("%d\n", p*(32-__builtin_clz(64-__builtin_clzll(p)))>>1);
		}

		cantdidate += 1;

		next:
		//clock_t total = clock() - total_start;
		mpz_nextprime(pprime, prime);
		mpz_swap(prime, pprime);
	}


	mpz_clears(prime, pprime, one, z2exp, candidate, z2p, z2pp1, k2p, k2pp1, gcd, sqrt, temp, m, tm, a, q, r, final, NULL);

	return 0;
}