Mercurial > dropbear

comparison etc/mersenne.c @ 2:86e0b50a9b58 libtommath-orig ltm-0.30-orig

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ltm 0.30 orig import
author Matt Johnston <matt@ucc.asn.au>
date Mon, 31 May 2004 18:25:22 +0000
parents
children d8254fc979e9
comparison
equal deleted inserted replaced
-1:000000000000 2:86e0b50a9b58
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 __N;
22 }
23
24 /* n = 2^s - 1 */
25 if ((res = mp_2expt(&n, s)) != MP_OKAY) {
26 goto __MU;
27 }
28 if ((res = mp_sub_d (&n, 1, &n)) != MP_OKAY) {
29 goto __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 __MU;
40 }
41 if ((res = mp_sub_d (&u, 2, &u)) != MP_OKAY) {
42 goto __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 __MU;
49 }
50 }
51
52 /* reduce */
53 if ((res = mp_reduce_2k (&u, &n, 1)) != MP_OKAY) {
54 goto __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 __MU:mp_clear (&u);
66 __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 }