|
1
|
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 __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 __V; |
|
|
201 } |
|
|
202 |
|
|
203 if ((res = mp_init (&a)) != MP_OKAY) { |
|
|
204 goto __V; |
|
|
205 } |
|
|
206 |
|
|
207 /* set the prime */ |
|
|
208 mp_set (&a, prime_digit ()); |
|
|
209 |
|
|
210 if ((res = mp_init (&b)) != MP_OKAY) { |
|
|
211 goto __A; |
|
|
212 } |
|
|
213 |
|
|
214 if ((res = mp_init (&n)) != MP_OKAY) { |
|
|
215 goto __B; |
|
|
216 } |
|
|
217 |
|
|
218 if ((res = mp_init (&x)) != MP_OKAY) { |
|
|
219 goto __N; |
|
|
220 } |
|
|
221 |
|
|
222 if ((res = mp_init (&y)) != MP_OKAY) { |
|
|
223 goto __X; |
|
|
224 } |
|
|
225 |
|
|
226 if ((res = mp_init (&z)) != MP_OKAY) { |
|
|
227 goto __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 __Z; |
|
|
240 } |
|
|
241 |
|
|
242 if ((res = mp_copy (&z, &c)) != MP_OKAY) { /* c = a * b */ |
|
|
243 goto __Z; |
|
|
244 } |
|
|
245 |
|
|
246 if ((res = mp_mul_2 (&z, &z)) != MP_OKAY) { /* z = 2 * a * b */ |
|
|
247 goto __Z; |
|
|
248 } |
|
|
249 |
|
|
250 /* n = z + 1 */ |
|
|
251 if ((res = mp_add_d (&z, 1, &n)) != MP_OKAY) { /* n = z + 1 */ |
|
|
252 goto __Z; |
|
|
253 } |
|
|
254 |
|
|
255 /* check (n, v) == 1 */ |
|
|
256 if ((res = mp_gcd (&n, &v, &y)) != MP_OKAY) { /* y = (n, v) */ |
|
|
257 goto __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 __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 __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 __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 __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 __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 __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 __Z:mp_clear (&z); |
|
|
350 __Y:mp_clear (&y); |
|
|
351 __X:mp_clear (&x); |
|
|
352 __N:mp_clear (&n); |
|
|
353 __B:mp_clear (&b); |
|
|
354 __A:mp_clear (&a); |
|
|
355 __V:mp_clear (&v); |
|
|
356 __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 } |
