userapps/opensource/sshd/libtommath/etc/mersenne.c

   1 /* Finds Mersenne primes using the Lucas-Lehmer test
   2  *
   3  * Tom St Denis, tomstdenis@iahu.ca
   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 }