userapps/opensource/sshd/libtommath/etc/pprime.c

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