Mercurial > dropbear

comparison libtommath/etc/mersenne.c @ 284:eed26cff980b

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 6c790cad5a7fa866ad062cb3a0c279f7ba788583) to branch 'au.asn.ucc.matt.dropbear' (head fff0894a0399405a9410ea1c6d118f342cf2aa64)
author Matt Johnston <matt@ucc.asn.au>
date Wed, 08 Mar 2006 13:23:49 +0000
parents
children 5ff8218bcee9
comparison
equal deleted inserted replaced
283:bd240aa12ba7 284:eed26cff980b
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 }