1/* Finds Mersenne primes using the Lucas-Lehmer test
2 *
3 * Tom St Denis, tomstdenis@gmail.com
4 */
5#include <time.h>
6#include <tommath.h>
7
8int
9is_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;
65LBL_MU:mp_clear (&u);
66LBL_N:mp_clear (&n);
67  return res;
68}
69
70/* square root of a long < 65536 */
71long
72i_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 */
90int
91isprime (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
104int
105main (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}
141
142/* $Source: /cvs/libtom/libtommath/etc/mersenne.c,v $ */
143/* $Revision: 1.3 $ */
144/* $Date: 2006/03/31 14:18:47 $ */
145