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