|
282
|
1 /* Finds Mersenne primes using the Lucas-Lehmer test |
|
|
2 * |
|
|
3 * Tom St Denis, [email protected] |
|
|
4 */ |
|
|
5 #include <time.h> |
|
|
6 #include <tommath.h> |
|
|
7 |
|
|
8 int |
|
|
9 is_mersenne (long s, int *pp) |
|
|
10 { |
|
|
11 mp_int n, u; |
|
|
12 int res, k; |
|
|
13 |
|
|
14 *pp = 0; |
|
|
15 |
|
|
16 if ((res = mp_init (&n)) != MP_OKAY) { |
|
|
17 return res; |
|
|
18 } |
|
|
19 |
|
|
20 if ((res = mp_init (&u)) != MP_OKAY) { |
|
|
21 goto LBL_N; |
|
|
22 } |
|
|
23 |
|
|
24 /* n = 2^s - 1 */ |
|
|
25 if ((res = mp_2expt(&n, s)) != MP_OKAY) { |
|
|
26 goto LBL_MU; |
|
|
27 } |
|
|
28 if ((res = mp_sub_d (&n, 1, &n)) != MP_OKAY) { |
|
|
29 goto LBL_MU; |
|
|
30 } |
|
|
31 |
|
|
32 /* set u=4 */ |
|
|
33 mp_set (&u, 4); |
|
|
34 |
|
|
35 /* for k=1 to s-2 do */ |
|
|
36 for (k = 1; k <= s - 2; k++) { |
|
|
37 /* u = u^2 - 2 mod n */ |
|
|
38 if ((res = mp_sqr (&u, &u)) != MP_OKAY) { |
|
|
39 goto LBL_MU; |
|
|
40 } |
|
|
41 if ((res = mp_sub_d (&u, 2, &u)) != MP_OKAY) { |
|
|
42 goto LBL_MU; |
|
|
43 } |
|
|
44 |
|
|
45 /* make sure u is positive */ |
|
|
46 while (u.sign == MP_NEG) { |
|
|
47 if ((res = mp_add (&u, &n, &u)) != MP_OKAY) { |
|
|
48 goto LBL_MU; |
|
|
49 } |
|
|
50 } |
|
|
51 |
|
|
52 /* reduce */ |
|
|
53 if ((res = mp_reduce_2k (&u, &n, 1)) != MP_OKAY) { |
|
|
54 goto LBL_MU; |
|
|
55 } |
|
|
56 } |
|
|
57 |
|
|
58 /* if u == 0 then its prime */ |
|
|
59 if (mp_iszero (&u) == 1) { |
|
|
60 mp_prime_is_prime(&n, 8, pp); |
|
|
61 if (*pp != 1) printf("FAILURE\n"); |
|
|
62 } |
|
|
63 |
|
|
64 res = MP_OKAY; |
|
|
65 LBL_MU:mp_clear (&u); |
|
|
66 LBL_N:mp_clear (&n); |
|
|
67 return res; |
|
|
68 } |
|
|
69 |
|
|
70 /* square root of a long < 65536 */ |
|
|
71 long |
|
|
72 i_sqrt (long x) |
|
|
73 { |
|
|
74 long x1, x2; |
|
|
75 |
|
|
76 x2 = 16; |
|
|
77 do { |
|
|
78 x1 = x2; |
|
|
79 x2 = x1 - ((x1 * x1) - x) / (2 * x1); |
|
|
80 } while (x1 != x2); |
|
|
81 |
|
|
82 if (x1 * x1 > x) { |
|
|
83 --x1; |
|
|
84 } |
|
|
85 |
|
|
86 return x1; |
|
|
87 } |
|
|
88 |
|
|
89 /* is the long prime by brute force */ |
|
|
90 int |
|
|
91 isprime (long k) |
|
|
92 { |
|
|
93 long y, z; |
|
|
94 |
|
|
95 y = i_sqrt (k); |
|
|
96 for (z = 2; z <= y; z++) { |
|
|
97 if ((k % z) == 0) |
|
|
98 return 0; |
|
|
99 } |
|
|
100 return 1; |
|
|
101 } |
|
|
102 |
|
|
103 |
|
|
104 int |
|
|
105 main (void) |
|
|
106 { |
|
|
107 int pp; |
|
|
108 long k; |
|
|
109 clock_t tt; |
|
|
110 |
|
|
111 k = 3; |
|
|
112 |
|
|
113 for (;;) { |
|
|
114 /* start time */ |
|
|
115 tt = clock (); |
|
|
116 |
|
|
117 /* test if 2^k - 1 is prime */ |
|
|
118 if (is_mersenne (k, &pp) != MP_OKAY) { |
|
|
119 printf ("Whoa error\n"); |
|
|
120 return -1; |
|
|
121 } |
|
|
122 |
|
|
123 if (pp == 1) { |
|
|
124 /* count time */ |
|
|
125 tt = clock () - tt; |
|
|
126 |
|
|
127 /* display if prime */ |
|
|
128 printf ("2^%-5ld - 1 is prime, test took %ld ticks\n", k, tt); |
|
|
129 } |
|
|
130 |
|
|
131 /* goto next odd exponent */ |
|
|
132 k += 2; |
|
|
133 |
|
|
134 /* but make sure its prime */ |
|
|
135 while (isprime (k) == 0) { |
|
|
136 k += 2; |
|
|
137 } |
|
|
138 } |
|
|
139 return 0; |
|
|
140 } |
