1 /* Generates provable primes 2 * 3 * See http://gmail.com:8080/papers/pp.pdf for more info. 4 * 5 * Tom St Denis, tomstdenis (at) gmail.com, http://tom.gmail.com 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 } 397 398 /* $Source: /cvs/libtom/libtommath/etc/pprime.c,v $ */ 399 /* $Revision: 1.3 $ */ 400 /* $Date: 2006/03/31 14:18:47 $ */ 401