lib/silcmath/silcprimegen.c

   1 /*
   2
   3   silcprimegen.c
   4
   5   Author: Pekka Riikonen <priikone@poseidon.pspt.fi>
   6
   7   Copyright (C) 1997 - 2000 Pekka Riikonen
   8
   9   This program is free software; you can redistribute it and/or modify
  10   it under the terms of the GNU General Public License as published by
  11   the Free Software Foundation; either version 2 of the License, or
  12   (at your option) any later version.
  13
  14   This program is distributed in the hope that it will be useful,
  15   but WITHOUT ANY WARRANTY; without even the implied warranty of
  16   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17   GNU General Public License for more details.
  18
  19 */
  20 /*
  21  * Created: Mon Dec  8 16:35:37 GMT+0200 1997
  22  */
  23 /*
  24  * $Id$
  25  * $Log$
  26  * Revision 1.1  2000/06/27 11:36:51  priikone
  27  * Initial revision
  28  *
  29  *
  30  */
  31
  32 #include "silcincludes.h"
  33
  34 /* XXX This must be temporary solution!! yucky! */
  35 /* Global random pool used for all prime generation. All primes generated
  36    in SILC uses this same pool. Before primes can be generated one must
  37    call silc_math_primegen_init. */
  38 static SilcRng primegen_rng;
  39
  40 /*
  41    Fixed primetable for small prime division. We use this primetable to
  42    test if possible prime is divisible any of these. Primetable is NULL
  43    terminated. This same primetable can be found in SSH and PGP except that
  44    SSH don't have prime 2 in the table. This sort of prime table is easily
  45    created, for example, using sieve of Erastosthenes.
  46
  47    Just to include; if somebody is interested about Erastosthenes. Creating a
  48    small prime table using sieve of Erastosthenes would go like this: Write
  49    down a list of integers in range of 2 to n. Here in example n = 20.
  50
  51    2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  52
  53    Then mark all multiples of 2 (Note: 2 is a prime number):
  54
  55    2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  56        x   x   x    x     x     x     x     x     x
  57
  58    Move to the next unmarked number, 3, then mark all multiples of 3.
  59
  60    2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  61        x   x   x x  x     x     x  x  x     x     x
  62
  63    And continue doing this, marking all multiples of the next unmarked
  64    number until there are now new unmarked numbers. And there you have a
  65    prime table. So in this example, primes between 2 - 20 are:
  66
  67    2 3 5 7 11 13 17 19
  68
  69 */
  70
  71 static unsigned int primetable[] =
  72 {
  73   2, 3, 5, 7, 11, 13, 17, 19,
  74   23, 29, 31, 37, 41, 43, 47, 53,
  75   59, 61, 67, 71, 73, 79, 83, 89,
  76   97, 101, 103, 107, 109, 113, 127, 131,
  77   137, 139, 149, 151, 157, 163, 167, 173,
  78   179, 181, 191, 193, 197, 199, 211, 223,
  79   227, 229, 233, 239, 241, 251, 257, 263,
  80   269, 271, 277, 281, 283, 293, 307, 311,
  81   313, 317, 331, 337, 347, 349, 353, 359,
  82   367, 373, 379, 383, 389, 397, 401, 409,
  83   419, 421, 431, 433, 439, 443, 449, 457,
  84   461, 463, 467, 479, 487, 491, 499, 503,
  85   509, 521, 523, 541, 547, 557, 563, 569,
  86   571, 577, 587, 593, 599, 601, 607, 613,
  87   617, 619, 631, 641, 643, 647, 653, 659,
  88   661, 673, 677, 683, 691, 701, 709, 719,
  89   727, 733, 739, 743, 751, 757, 761, 769,
  90   773, 787, 797, 809, 811, 821, 823, 827,
  91   829, 839, 853, 857, 859, 863, 877, 881,
  92   883, 887, 907, 911, 919, 929, 937, 941,
  93   947, 953, 967, 971, 977, 983, 991, 997,
  94   1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049,
  95   1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
  96   1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163,
  97   1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
  98   1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283,
  99   1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
 100   1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
 101   1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
 102   1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
 103   1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571,
 104   1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619,
 105   1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
 106   1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
 107   1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
 108   1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
 109   1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949,
 110   1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
 111   2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069,
 112   2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
 113   2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
 114   2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
 115   2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311,
 116   2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377,
 117   2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
 118   2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
 119   2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
 120   2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657,
 121   2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
 122   2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
 123   2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801,
 124   2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
 125   2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939,
 126   2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
 127   3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
 128   3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167,
 129   3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
 130   3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301,
 131   3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347,
 132   3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
 133   3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491,
 134   3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541,
 135   3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
 136   3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671,
 137   3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
 138   3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797,
 139   3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863,
 140   3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
 141   3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003,
 142   4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
 143   4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,
 144   4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211,
 145   4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259,
 146   4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337,
 147   4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
 148   4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481,
 149   4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547,
 150   4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621,
 151   4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673,
 152   4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
 153   4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813,
 154   4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909,
 155   4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967,
 156   4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011,
 157   5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
 158   5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167,
 159   5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233,
 160   5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309,
 161   5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399,
 162   5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
 163   5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507,
 164   5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573,
 165   5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653,
 166   5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711,
 167   5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
 168   5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849,
 169   5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897,
 170   5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007,
 171   6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073,
 172   6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
 173   6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211,
 174   6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271,
 175   6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329,
 176   6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379,
 177   6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
 178   6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563,
 179   6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637,
 180   6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701,
 181   6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779,
 182   6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
 183   6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907,
 184   6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971,
 185   6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027,
 186   7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121,
 187   7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
 188   7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253,
 189   7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349,
 190   7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457,
 191   7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517,
 192   7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
 193   7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621,
 194   7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691,
 195   7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757,
 196   7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853,
 197   7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
 198   7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009,
 199   8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087,
 200   8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161,
 201   8167, 8171, 8179, 8191, 0
 202 };
 203
 204
 205 /* Find appropriate prime. It generates a number by taking random bytes.
 206    It then tests the number that it's not divisible by any of the small
 207    primes and then it performs Fermat's prime test. I thank Rieks Joosten
 208    (r.joosten@pijnenburg.nl) for such a good help with prime tests.
 209
 210    If argument verbose is TRUE this will display some status information
 211    about the progress of generation. */
 212
 213 int silc_math_gen_prime(SilcInt *prime, unsigned int bits, int verbose)
 214 {
 215   unsigned char *numbuf;
 216   unsigned int i, b, k;
 217   unsigned int *spmods;
 218   SilcInt r, base, tmp, tmp2, oprime;
 219
 220   /* XXX */
 221   assert(primegen_rng != NULL);
 222
 223   silc_mp_init(&r);
 224   silc_mp_init_set_ui(&base, 2);
 225   silc_mp_init(&tmp);
 226   silc_mp_init(&tmp2);
 227   silc_mp_init(&oprime);
 228
 229   SILC_LOG_DEBUG(("Generating new prime"));
 230
 231   /* Get random number */
 232   numbuf = silc_rng_get_rn_string(primegen_rng, (bits / 8));
 233   if (!numbuf)
 234     return FALSE;
 235
 236   /* Convert into MP and set the size */
 237   silc_mp_set_str(prime, numbuf, 16);
 238   silc_mp_mod_2exp(prime, prime, bits);
 239
 240   /* Empty buffer */
 241   memset(numbuf, 0, (bits / 8));
 242   silc_free(numbuf);
 243
 244   /* Number could be even number, so we'll make it odd. */
 245   silc_mp_set_ui(&tmp, 1);
 246   silc_mp_ior(prime, prime, &tmp);              /* OR operator */
 247
 248   /* Init modulo table with the prime candidate and the primes
 249      in the primetable. */
 250   spmods = silc_malloc(sizeof(primetable) * sizeof(unsigned int));
 251   for (i = 0; primetable[i] != 0; i++) {
 252     silc_mp_mod_ui(&tmp, prime, primetable[i]);
 253     spmods[i] = silc_mp_get_ui(&tmp);
 254   }
 255
 256   /* k is added by 2, this way we skip all even numbers (prime is odd). */
 257   silc_mp_set(&oprime, prime);
 258   for (k = 0;; k += 2) {
 259     silc_mp_add_ui(&oprime, prime, k);
 260
 261     /* See if the prime has a divisor in primetable[].
 262      * If it has then it's a composite. We add k to the
 263      * original modulo value, k is added by 2 on every roll.
 264      * This way we don't have to re-init the whole table if
 265      * the number is composite.
 266      */
 267     for (b = 0; b < i; b++) {
 268       silc_mp_set_ui(&tmp2, spmods[b]);
 269       silc_mp_add_ui(&tmp2, &tmp2, k);
 270       silc_mp_mod_ui(&tmp, &tmp2, primetable[b]);
 271
 272       /* If mod is 0, the number is composite */
 273       if (silc_mp_cmp_ui(&tmp, 0) == 0)
 274         break;
 275     }
 276     if (b < i)
 277       continue;
 278
 279     /* Passed the quick test. */
 280
 281     /* Does the prime pass the Fermat's prime test.
 282      * r = 2 ^ p mod p, if r == 2, then p is probably a prime.
 283      */
 284     silc_mp_powm(&r, &base, &oprime, &oprime);
 285     if (silc_mp_cmp_ui(&r, 2) != 0) {
 286       if (verbose) {
 287         printf(".");
 288         fflush(stdout);
 289       }
 290       continue;
 291     }
 292
 293     /* Passed the Fermat's test. And I don't think
 294      * we have to do more tests. If anyone wants to do more
 295      * tests, MP library has probability of prime test:
 296      *
 297      * if (silc_mp_probab_prime_p(prime, 25))
 298      *     break;                          found a probable prime
 299      *
 300      * However, this is very slow.
 301      */
 302     /* XXX: this could be done if some argument, say strict_primes, is
 303        TRUE when we are willing to spend more time on the prime test and
 304        to get, perhaps, better primes. */
 305     silc_mp_set(prime, &oprime);
 306     break;
 307   }
 308
 309   silc_free(spmods);
 310   silc_mp_clear(&r);
 311   silc_mp_clear(&base);
 312   silc_mp_clear(&tmp);
 313   silc_mp_clear(&tmp2);
 314   silc_mp_clear(&oprime);
 315
 316   return TRUE;
 317 }
 318
 319 /* Performs primality testings for given number. Returns TRUE if the
 320    number is probably a prime. */
 321
 322 int silc_math_prime_test(SilcInt *p)
 323 {
 324   SilcInt r, base, tmp;
 325   int i, ret = 0;
 326
 327   silc_mp_init(&r);
 328   silc_mp_init(&tmp);
 329   silc_mp_init_set_ui(&base, 2);
 330
 331   SILC_LOG_DEBUG(("Testing probability of prime"));
 332
 333   /* See if the number is divisible by any of the
 334      small primes in primetable[]. */
 335   for (i = 0; primetable[i] != 0; i++) {
 336     silc_mp_mod_ui(&tmp, p, primetable[i]);
 337
 338     /* If mod is 0, the number is composite */
 339     if (silc_mp_cmp_ui(&tmp, 0) == 0)
 340       ret = -1;
 341   }
 342
 343   /* Does the prime pass the Fermat's prime test.
 344    * r = 2 ^ p mod p, if r == 2, then p is probably a prime.
 345    */
 346   silc_mp_powm(&r, &base, p, p);
 347   if (silc_mp_cmp_ui(&r, 2) != 0)
 348     ret = -1;
 349
 350   silc_mp_clear(&r);
 351   silc_mp_clear(&tmp);
 352   silc_mp_clear(&base);
 353
 354   if (ret)
 355     return FALSE;
 356
 357   /* Number is probably a prime */
 358   return TRUE;
 359 }
 360
 361 /* XXX This must temporary solution!! */
 362 /* Initializes the random pool used to generated primes */
 363
 364 void silc_math_primegen_init()
 365 {
 366   SILC_LOG_DEBUG(("Start"));
 367
 368   if (primegen_rng == NULL) {
 369     primegen_rng = silc_rng_alloc();
 370     silc_rng_init(primegen_rng);
 371   }
 372 }
 373
 374 /* XXX This must temporary solution!! */
 375 /* Uninitializes random pool */
 376
 377 void silc_math_primegen_uninit()
 378 {
 379   silc_rng_free(primegen_rng);
 380 }