2^11213-1 prime in 10 seconds

1 files changed, 90 insertions, 0 deletions

diff --git a/mersenne_tricks_nogcd_wtest.c b/mersenne_tricks_nogcd_wtest.c
new file mode 100644
index 0000000..1f13949
--- /dev/null
+++ b/mersenne_tricks_nogcd_wtest.c
@@ -0,0 +1,90 @@
+#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;
+}