Home | History | Annotate | Download | only in etc 
      1 /* Finds Mersenne primes using the Lucas-Lehmer test
      2  *
      3  * Tom St Denis, tomstdenis (at) gmail.com
      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 }
    141 
    142 /* $Source: /cvs/libtom/libtommath/etc/mersenne.c,v $ */
    143 /* $Revision: 1.3 $ */
    144 /* $Date: 2006/03/31 14:18:47 $ */
    145