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