Mercurial > dropbear

comparison libtommath/etc/pprime.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 /* Generates provable primes
2 *
3 * See http://iahu.ca:8080/papers/pp.pdf for more info.
4 *
5 * Tom St Denis, [email protected], http://tom.iahu.ca
6 */
7 #include <time.h>
8 #include "tommath.h"
9
10 int n_prime;
11 FILE *primes;
12
13 /* fast square root */
14 static mp_digit
15 i_sqrt (mp_word x)
16 {
17 mp_word x1, x2;
18
19 x2 = x;
20 do {
21 x1 = x2;
22 x2 = x1 - ((x1 * x1) - x) / (2 * x1);
23 } while (x1 != x2);
24
25 if (x1 * x1 > x) {
26 --x1;
27 }
28
29 return x1;
30 }
31
32
33 /* generates a prime digit */
34 static void gen_prime (void)
35 {
36 mp_digit r, x, y, next;
37 FILE *out;
38
39 out = fopen("pprime.dat", "wb");
40
41 /* write first set of primes */
42 r = 3; fwrite(&r, 1, sizeof(mp_digit), out);
43 r = 5; fwrite(&r, 1, sizeof(mp_digit), out);
44 r = 7; fwrite(&r, 1, sizeof(mp_digit), out);
45 r = 11; fwrite(&r, 1, sizeof(mp_digit), out);
46 r = 13; fwrite(&r, 1, sizeof(mp_digit), out);
47 r = 17; fwrite(&r, 1, sizeof(mp_digit), out);
48 r = 19; fwrite(&r, 1, sizeof(mp_digit), out);
49 r = 23; fwrite(&r, 1, sizeof(mp_digit), out);
50 r = 29; fwrite(&r, 1, sizeof(mp_digit), out);
51 r = 31; fwrite(&r, 1, sizeof(mp_digit), out);
52
53 /* get square root, since if 'r' is composite its factors must be < than this */
54 y = i_sqrt (r);
55 next = (y + 1) * (y + 1);
56
57 for (;;) {
58 do {
59 r += 2; /* next candidate */
60 r &= MP_MASK;
61 if (r < 31) break;
62
63 /* update sqrt ? */
64 if (next <= r) {
65 ++y;
66 next = (y + 1) * (y + 1);
67 }
68
69 /* loop if divisible by 3,5,7,11,13,17,19,23,29 */
70 if ((r % 3) == 0) {
71 x = 0;
72 continue;
73 }
74 if ((r % 5) == 0) {
75 x = 0;
76 continue;
77 }
78 if ((r % 7) == 0) {
79 x = 0;
80 continue;
81 }
82 if ((r % 11) == 0) {
83 x = 0;
84 continue;
85 }
86 if ((r % 13) == 0) {
87 x = 0;
88 continue;
89 }
90 if ((r % 17) == 0) {
91 x = 0;
92 continue;
93 }
94 if ((r % 19) == 0) {
95 x = 0;
96 continue;
97 }
98 if ((r % 23) == 0) {
99 x = 0;
100 continue;
101 }
102 if ((r % 29) == 0) {
103 x = 0;
104 continue;
105 }
106
107 /* now check if r is divisible by x + k={1,7,11,13,17,19,23,29} */
108 for (x = 30; x <= y; x += 30) {
109 if ((r % (x + 1)) == 0) {
110 x = 0;
111 break;
112 }
113 if ((r % (x + 7)) == 0) {
114 x = 0;
115 break;
116 }
117 if ((r % (x + 11)) == 0) {
118 x = 0;
119 break;
120 }
121 if ((r % (x + 13)) == 0) {
122 x = 0;
123 break;
124 }
125 if ((r % (x + 17)) == 0) {
126 x = 0;
127 break;
128 }
129 if ((r % (x + 19)) == 0) {
130 x = 0;
131 break;
132 }
133 if ((r % (x + 23)) == 0) {
134 x = 0;
135 break;
136 }
137 if ((r % (x + 29)) == 0) {
138 x = 0;
139 break;
140 }
141 }
142 } while (x == 0);
143 if (r > 31) { fwrite(&r, 1, sizeof(mp_digit), out); printf("%9d\r", r); fflush(stdout); }
144 if (r < 31) break;
145 }
146
147 fclose(out);
148 }
149
150 void load_tab(void)
151 {
152 primes = fopen("pprime.dat", "rb");
153 if (primes == NULL) {
154 gen_prime();
155 primes = fopen("pprime.dat", "rb");
156 }
157 fseek(primes, 0, SEEK_END);
158 n_prime = ftell(primes) / sizeof(mp_digit);
159 }
160
161 mp_digit prime_digit(void)
162 {
163 int n;
164 mp_digit d;
165
166 n = abs(rand()) % n_prime;
167 fseek(primes, n * sizeof(mp_digit), SEEK_SET);
168 fread(&d, 1, sizeof(mp_digit), primes);
169 return d;
170 }
171
172
173 /* makes a prime of at least k bits */
174 int
175 pprime (int k, int li, mp_int * p, mp_int * q)
176 {
177 mp_int a, b, c, n, x, y, z, v;
178 int res, ii;
179 static const mp_digit bases[] = { 2, 3, 5, 7, 11, 13, 17, 19 };
180
181 /* single digit ? */
182 if (k <= (int) DIGIT_BIT) {
183 mp_set (p, prime_digit ());
184 return MP_OKAY;
185 }
186
187 if ((res = mp_init (&c)) != MP_OKAY) {
188 return res;
189 }
190
191 if ((res = mp_init (&v)) != MP_OKAY) {
192 goto LBL_C;
193 }
194
195 /* product of first 50 primes */
196 if ((res =
197 mp_read_radix (&v,
198 "19078266889580195013601891820992757757219839668357012055907516904309700014933909014729740190",
199 10)) != MP_OKAY) {
200 goto LBL_V;
201 }
202
203 if ((res = mp_init (&a)) != MP_OKAY) {
204 goto LBL_V;
205 }
206
207 /* set the prime */
208 mp_set (&a, prime_digit ());
209
210 if ((res = mp_init (&b)) != MP_OKAY) {
211 goto LBL_A;
212 }
213
214 if ((res = mp_init (&n)) != MP_OKAY) {
215 goto LBL_B;
216 }
217
218 if ((res = mp_init (&x)) != MP_OKAY) {
219 goto LBL_N;
220 }
221
222 if ((res = mp_init (&y)) != MP_OKAY) {
223 goto LBL_X;
224 }
225
226 if ((res = mp_init (&z)) != MP_OKAY) {
227 goto LBL_Y;
228 }
229
230 /* now loop making the single digit */
231 while (mp_count_bits (&a) < k) {
232 fprintf (stderr, "prime has %4d bits left\r", k - mp_count_bits (&a));
233 fflush (stderr);
234 top:
235 mp_set (&b, prime_digit ());
236
237 /* now compute z = a * b * 2 */
238 if ((res = mp_mul (&a, &b, &z)) != MP_OKAY) { /* z = a * b */
239 goto LBL_Z;
240 }
241
242 if ((res = mp_copy (&z, &c)) != MP_OKAY) { /* c = a * b */
243 goto LBL_Z;
244 }
245
246 if ((res = mp_mul_2 (&z, &z)) != MP_OKAY) { /* z = 2 * a * b */
247 goto LBL_Z;
248 }
249
250 /* n = z + 1 */
251 if ((res = mp_add_d (&z, 1, &n)) != MP_OKAY) { /* n = z + 1 */
252 goto LBL_Z;
253 }
254
255 /* check (n, v) == 1 */
256 if ((res = mp_gcd (&n, &v, &y)) != MP_OKAY) { /* y = (n, v) */
257 goto LBL_Z;
258 }
259
260 if (mp_cmp_d (&y, 1) != MP_EQ)
261 goto top;
262
263 /* now try base x=bases[ii] */
264 for (ii = 0; ii < li; ii++) {
265 mp_set (&x, bases[ii]);
266
267 /* compute x^a mod n */
268 if ((res = mp_exptmod (&x, &a, &n, &y)) != MP_OKAY) { /* y = x^a mod n */
269 goto LBL_Z;
270 }
271
272 /* if y == 1 loop */
273 if (mp_cmp_d (&y, 1) == MP_EQ)
274 continue;
275
276 /* now x^2a mod n */
277 if ((res = mp_sqrmod (&y, &n, &y)) != MP_OKAY) { /* y = x^2a mod n */
278 goto LBL_Z;
279 }
280
281 if (mp_cmp_d (&y, 1) == MP_EQ)
282 continue;
283
284 /* compute x^b mod n */
285 if ((res = mp_exptmod (&x, &b, &n, &y)) != MP_OKAY) { /* y = x^b mod n */
286 goto LBL_Z;
287 }
288
289 /* if y == 1 loop */
290 if (mp_cmp_d (&y, 1) == MP_EQ)
291 continue;
292
293 /* now x^2b mod n */
294 if ((res = mp_sqrmod (&y, &n, &y)) != MP_OKAY) { /* y = x^2b mod n */
295 goto LBL_Z;
296 }
297
298 if (mp_cmp_d (&y, 1) == MP_EQ)
299 continue;
300
301 /* compute x^c mod n == x^ab mod n */
302 if ((res = mp_exptmod (&x, &c, &n, &y)) != MP_OKAY) { /* y = x^ab mod n */
303 goto LBL_Z;
304 }
305
306 /* if y == 1 loop */
307 if (mp_cmp_d (&y, 1) == MP_EQ)
308 continue;
309
310 /* now compute (x^c mod n)^2 */
311 if ((res = mp_sqrmod (&y, &n, &y)) != MP_OKAY) { /* y = x^2ab mod n */
312 goto LBL_Z;
313 }
314
315 /* y should be 1 */
316 if (mp_cmp_d (&y, 1) != MP_EQ)
317 continue;
318 break;
319 }
320
321 /* no bases worked? */
322 if (ii == li)
323 goto top;
324
325 {
326 char buf[4096];
327
328 mp_toradix(&n, buf, 10);
329 printf("Certificate of primality for:\n%s\n\n", buf);
330 mp_toradix(&a, buf, 10);
331 printf("A == \n%s\n\n", buf);
332 mp_toradix(&b, buf, 10);
333 printf("B == \n%s\n\nG == %d\n", buf, bases[ii]);
334 printf("----------------------------------------------------------------\n");
335 }
336
337 /* a = n */
338 mp_copy (&n, &a);
339 }
340
341 /* get q to be the order of the large prime subgroup */
342 mp_sub_d (&n, 1, q);
343 mp_div_2 (q, q);
344 mp_div (q, &b, q, NULL);
345
346 mp_exch (&n, p);
347
348 res = MP_OKAY;
349 LBL_Z:mp_clear (&z);
350 LBL_Y:mp_clear (&y);
351 LBL_X:mp_clear (&x);
352 LBL_N:mp_clear (&n);
353 LBL_B:mp_clear (&b);
354 LBL_A:mp_clear (&a);
355 LBL_V:mp_clear (&v);
356 LBL_C:mp_clear (&c);
357 return res;
358 }
359
360
361 int
362 main (void)
363 {
364 mp_int p, q;
365 char buf[4096];
366 int k, li;
367 clock_t t1;
368
369 srand (time (NULL));
370 load_tab();
371
372 printf ("Enter # of bits: \n");
373 fgets (buf, sizeof (buf), stdin);
374 sscanf (buf, "%d", &k);
375
376 printf ("Enter number of bases to try (1 to 8):\n");
377 fgets (buf, sizeof (buf), stdin);
378 sscanf (buf, "%d", &li);
379
380
381 mp_init (&p);
382 mp_init (&q);
383
384 t1 = clock ();
385 pprime (k, li, &p, &q);
386 t1 = clock () - t1;
387
388 printf ("\n\nTook %ld ticks, %d bits\n", t1, mp_count_bits (&p));
389
390 mp_toradix (&p, buf, 10);
391 printf ("P == %s\n", buf);
392 mp_toradix (&q, buf, 10);
393 printf ("Q == %s\n", buf);
394
395 return 0;
396 }