crypto/external/bsd/heimdal/dist/lib/hcrypto/libtommath/etc/pprime.c

   1 /*      $NetBSD: pprime.c,v 1.1.1.2 2014/04/24 12:45:39 pettai 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  */