Mercurial > dropbear
annotate etc/mersenne.c @ 190:d8254fc979e9 libtommath-orig LTM_0.35
Initial import of libtommath 0.35
| author | Matt Johnston <matt@ucc.asn.au> |
|---|---|
| date | Fri, 06 May 2005 08:59:30 +0000 |
| parents | 86e0b50a9b58 |
| children |
| rev | line source |
|---|---|
| 2 | 1 /* Finds Mersenne primes using the Lucas-Lehmer test |
| 2 * | |
| 3 * Tom St Denis, [email protected] | |
| 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
21 goto LBL_N; |
| 2 | 22 } |
| 23 | |
| 24 /* n = 2^s - 1 */ | |
| 25 if ((res = mp_2expt(&n, s)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
26 goto LBL_MU; |
| 2 | 27 } |
| 28 if ((res = mp_sub_d (&n, 1, &n)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
29 goto LBL_MU; |
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
39 goto LBL_MU; |
| 2 | 40 } |
| 41 if ((res = mp_sub_d (&u, 2, &u)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
42 goto LBL_MU; |
| 2 | 43 } |
| 44 | |
| 45 /* make sure u is positive */ | |
| 46 while (u.sign == MP_NEG) { | |
| 47 if ((res = mp_add (&u, &n, &u)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
48 goto LBL_MU; |
| 2 | 49 } |
| 50 } | |
| 51 | |
| 52 /* reduce */ | |
| 53 if ((res = mp_reduce_2k (&u, &n, 1)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
54 goto LBL_MU; |
| 2 | 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; | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
65 LBL_MU:mp_clear (&u); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
66 LBL_N:mp_clear (&n); |
| 2 | 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 } |