1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
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;
}
|
