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