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

ltm 0.30 orig import
author Matt Johnston <matt@ucc.asn.au>
date Mon, 31 May 2004 18:25:22 +0000
parents
children d8254fc979e9
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
1 /* Finds Mersenne primes using the Lucas-Lehmer test
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
2 *
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
3 * Tom St Denis, [email protected]
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
4 */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
5 #include <time.h>
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
6 #include <tommath.h>
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
7
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
8 int
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
9 is_mersenne (long s, int *pp)
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
10 {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
11 mp_int n, u;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
12 int res, k;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
13
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
14 *pp = 0;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
15
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
16 if ((res = mp_init (&n)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
17 return res;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
18 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
19
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
20 if ((res = mp_init (&u)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
21 goto __N;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
22 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
23
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
24 /* n = 2^s - 1 */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
25 if ((res = mp_2expt(&n, s)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
26 goto __MU;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
27 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
28 if ((res = mp_sub_d (&n, 1, &n)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
29 goto __MU;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
30 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
31
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
32 /* set u=4 */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
33 mp_set (&u, 4);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
34
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
35 /* for k=1 to s-2 do */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
36 for (k = 1; k <= s - 2; k++) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
37 /* u = u^2 - 2 mod n */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
38 if ((res = mp_sqr (&u, &u)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
39 goto __MU;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
40 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
41 if ((res = mp_sub_d (&u, 2, &u)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
42 goto __MU;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
43 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
44
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
45 /* make sure u is positive */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
46 while (u.sign == MP_NEG) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
47 if ((res = mp_add (&u, &n, &u)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
48 goto __MU;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
49 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
50 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
51
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
52 /* reduce */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
53 if ((res = mp_reduce_2k (&u, &n, 1)) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
54 goto __MU;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
55 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
56 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
57
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
58 /* if u == 0 then its prime */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
59 if (mp_iszero (&u) == 1) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
60 mp_prime_is_prime(&n, 8, pp);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
61 if (*pp != 1) printf("FAILURE\n");
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
62 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
63
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
64 res = MP_OKAY;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
65 __MU:mp_clear (&u);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
66 __N:mp_clear (&n);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
67 return res;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
68 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
69
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
70 /* square root of a long < 65536 */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
71 long
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
72 i_sqrt (long x)
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
73 {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
74 long x1, x2;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
75
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
76 x2 = 16;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
77 do {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
78 x1 = x2;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
79 x2 = x1 - ((x1 * x1) - x) / (2 * x1);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
80 } while (x1 != x2);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
81
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
82 if (x1 * x1 > x) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
83 --x1;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
84 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
85
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
86 return x1;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
87 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
88
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
89 /* is the long prime by brute force */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
90 int
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
91 isprime (long k)
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
92 {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
93 long y, z;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
94
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
95 y = i_sqrt (k);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
96 for (z = 2; z <= y; z++) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
97 if ((k % z) == 0)
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
98 return 0;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
99 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
100 return 1;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
101 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
102
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
103
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
104 int
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
105 main (void)
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
106 {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
107 int pp;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
108 long k;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
109 clock_t tt;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
110
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
111 k = 3;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
112
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
113 for (;;) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
114 /* start time */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
115 tt = clock ();
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
116
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
117 /* test if 2^k - 1 is prime */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
118 if (is_mersenne (k, &pp) != MP_OKAY) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
119 printf ("Whoa error\n");
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
120 return -1;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
121 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
122
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
123 if (pp == 1) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
124 /* count time */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
125 tt = clock () - tt;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
126
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
127 /* display if prime */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
128 printf ("2^%-5ld - 1 is prime, test took %ld ticks\n", k, tt);
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
129 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
130
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
131 /* goto next odd exponent */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
132 k += 2;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
133
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
134 /* but make sure its prime */
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
135 while (isprime (k) == 0) {
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
136 k += 2;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
137 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
138 }
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
139 return 0;
86e0b50a9b58 ltm 0.30 orig import
Matt Johnston <matt@ucc.asn.au>
parents:
diff changeset
140 }