Mercurial > dropbear
annotate etc/pprime.c @ 200:c5c969ed76f3 libtommath
propagate from branch 'au.asn.ucc.matt.ltm-orig' (head 7fa10cba9535de3461cedb14b877c24858826204)
to branch 'au.asn.ucc.matt.dropbear.ltm' (head fc26f60de0370ab0a281fa41a2d13fb17c9d90a8)
| author | Matt Johnston <matt@ucc.asn.au> |
|---|---|
| date | Wed, 11 May 2005 16:15:27 +0000 |
| parents | d8254fc979e9 |
| children |
| rev | line source |
|---|---|
| 2 | 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) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
192 goto LBL_C; |
| 2 | 193 } |
| 194 | |
| 195 /* product of first 50 primes */ | |
| 196 if ((res = | |
| 197 mp_read_radix (&v, | |
| 198 "19078266889580195013601891820992757757219839668357012055907516904309700014933909014729740190", | |
| 199 10)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
200 goto LBL_V; |
| 2 | 201 } |
| 202 | |
| 203 if ((res = mp_init (&a)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
204 goto LBL_V; |
| 2 | 205 } |
| 206 | |
| 207 /* set the prime */ | |
| 208 mp_set (&a, prime_digit ()); | |
| 209 | |
| 210 if ((res = mp_init (&b)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
211 goto LBL_A; |
| 2 | 212 } |
| 213 | |
| 214 if ((res = mp_init (&n)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
215 goto LBL_B; |
| 2 | 216 } |
| 217 | |
| 218 if ((res = mp_init (&x)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
219 goto LBL_N; |
| 2 | 220 } |
| 221 | |
| 222 if ((res = mp_init (&y)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
223 goto LBL_X; |
| 2 | 224 } |
| 225 | |
| 226 if ((res = mp_init (&z)) != MP_OKAY) { | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
227 goto LBL_Y; |
| 2 | 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 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
239 goto LBL_Z; |
| 2 | 240 } |
| 241 | |
| 242 if ((res = mp_copy (&z, &c)) != MP_OKAY) { /* c = a * b */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
243 goto LBL_Z; |
| 2 | 244 } |
| 245 | |
| 246 if ((res = mp_mul_2 (&z, &z)) != MP_OKAY) { /* z = 2 * a * b */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
247 goto LBL_Z; |
| 2 | 248 } |
| 249 | |
| 250 /* n = z + 1 */ | |
| 251 if ((res = mp_add_d (&z, 1, &n)) != MP_OKAY) { /* n = z + 1 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
252 goto LBL_Z; |
| 2 | 253 } |
| 254 | |
| 255 /* check (n, v) == 1 */ | |
| 256 if ((res = mp_gcd (&n, &v, &y)) != MP_OKAY) { /* y = (n, v) */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
257 goto LBL_Z; |
| 2 | 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 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
269 goto LBL_Z; |
| 2 | 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 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
278 goto LBL_Z; |
| 2 | 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 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
286 goto LBL_Z; |
| 2 | 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 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
295 goto LBL_Z; |
| 2 | 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 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
303 goto LBL_Z; |
| 2 | 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 */ | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
312 goto LBL_Z; |
| 2 | 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; | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
349 LBL_Z:mp_clear (&z); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
350 LBL_Y:mp_clear (&y); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
351 LBL_X:mp_clear (&x); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
352 LBL_N:mp_clear (&n); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
353 LBL_B:mp_clear (&b); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
354 LBL_A:mp_clear (&a); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
355 LBL_V:mp_clear (&v); |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
2
diff
changeset
|
356 LBL_C:mp_clear (&c); |
| 2 | 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 } |