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