crypto/external/bsd/heimdal/dist/lib/hcrypto/libtommath/etc/mersenne.c

   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
  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
  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
  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
 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  */