Mercurial > dropbear

comparison libtommath/etc/mersenne.c @ 389:5ff8218bcee9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
propagate from branch 'au.asn.ucc.matt.ltm.dropbear' (head 2af95f00ebd5bb7a28b3817db1218442c935388e) to branch 'au.asn.ucc.matt.dropbear' (head ecd779509ef23a8cdf64888904fc9b31d78aa933)
author Matt Johnston <matt@ucc.asn.au>
date Thu, 11 Jan 2007 03:14:55 +0000
parents eed26cff980b
children 60fc6476e044
comparison
equal deleted inserted replaced
388:fb54020f78e1 389:5ff8218bcee9
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) {
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 $ */