lib/silcmath/silcprimegen.c

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