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