lib/hcrypto/libtommath/etc/pprime.c

   1 /* Generates provable primes
   2  *
   3  * See http://gmail.com:8080/papers/pp.pdf for more info.
   4  *
   5  * Tom St Denis, tomstdenis@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 = labs(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 $ */