1/*	$NetBSD: mersenne.c,v 1.1.1.1 2011/04/13 18:15:06 elric 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
10int
11is_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;
67LBL_MU:mp_clear (&u);
68LBL_N:mp_clear (&n);
69  return res;
70}
71
72/* square root of a long < 65536 */
73long
74i_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 */
92int
93isprime (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
106int
107main (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