lib/silcmath/silcprimegen.c

   1 /*
   2
   3   silcprimegen.c
   4
   5   Author: Pekka Riikonen <priikone@silcnet.org>
   6
   7   Copyright (C) 1997 - 2007 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;
 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     numbuf = silc_malloc((((bits + 7) / 8) + 1) * sizeof(*numbuf));
 218     if (!numbuf)
 219       return FALSE;
 220
 221     /* Get random number */
 222     if (rng)
 223       silc_rng_get_rn_data(rng, (bits / 8), numbuf, (bits / 8));
 224     else
 225       silc_rng_global_get_rn_data(rng, (bits / 8), numbuf, (bits / 8));
 226
 227     /* Convert into MP and set the size */
 228     silc_mp_bin2mp(numbuf, (bits / 8), prime);
 229     silc_mp_mod_2exp(prime, prime, bits);
 230
 231     /* Empty buffer */
 232     memset(numbuf, 0, (bits / 8));
 233     silc_free(numbuf);
 234
 235     /* Set highest bit */
 236     silc_mp_set_ui(&tmp, 1);
 237     silc_mp_mul_2exp(&tmp, &tmp, bits - 1);
 238     silc_mp_or(prime, prime, &tmp);
 239
 240     /* Number could be even number, so we'll make it odd. */
 241     silc_mp_set_ui(&tmp, 1);
 242     silc_mp_or(prime, prime, &tmp);
 243
 244     /* Init modulo table with the prime candidate and the primes
 245        in the primetable. */
 246     spmods = silc_calloc(1, sizeof(primetable) * sizeof(SilcUInt32));
 247     for (i = 0; primetable[i] != 0; i++) {
 248       silc_mp_mod_ui(&tmp, prime, primetable[i]);
 249       spmods[i] = silc_mp_get_ui(&tmp);
 250     }
 251
 252     /* k is added by 2, this way we skip all even numbers (prime is odd). */
 253     silc_mp_set(&oprime, prime);
 254     for (k = 0;; k += 2) {
 255       silc_mp_add_ui(&oprime, prime, k);
 256
 257       /* See if the prime has a divisor in primetable[].
 258        * If it has then it's a composite. We add k to the
 259        * original modulo value, k is added by 2 on every roll.
 260        * This way we don't have to re-init the whole table if
 261        * the number is composite.
 262        */
 263       for (b = 0; b < i; b++) {
 264         silc_mp_set_ui(&tmp2, spmods[b]);
 265         silc_mp_add_ui(&tmp2, &tmp2, k);
 266         silc_mp_mod_ui(&tmp, &tmp2, primetable[b]);
 267
 268         /* If mod is 0, the number is composite */
 269         if (silc_mp_cmp_ui(&tmp, 0) == 0)
 270           break;
 271       }
 272       if (b < i)
 273         continue;
 274
 275       /* Passed the quick test. */
 276
 277       /* Does the prime pass the Fermat's prime test.
 278        * r = 2 ^ p mod p, if r == 2, then p is probably a prime.
 279        */
 280       silc_mp_pow_mod(&r, &base, &oprime, &oprime);
 281       if (silc_mp_cmp_ui(&r, 2) != 0) {
 282         if (verbose) {
 283           printf(".");
 284           fflush(stdout);
 285         }
 286         continue;
 287       }
 288
 289       /* Passed the Fermat's test. And I don't think
 290        * we have to do more tests. If anyone wants to do more
 291        * tests, MP library has probability of prime test:
 292        *
 293        * if (silc_mp_probab_prime_p(prime, 25))
 294        *     break;                          found a probable prime
 295        *
 296        * However, this is very slow.
 297        */
 298       /* XXX: this could be done if some argument, say strict_primes, is
 299          TRUE when we are willing to spend more time on the prime test and
 300          to get, perhaps, better primes. */
 301       silc_mp_set(prime, &oprime);
 302       break;
 303     }
 304
 305     /* Check highest bit */
 306     silc_mp_div_2exp(&tmp, prime, bits - 1);
 307     if (silc_mp_get_ui(&tmp) == 1) {
 308       valid = TRUE;
 309       break;
 310     }
 311   }
 312
 313   silc_free(spmods);
 314   silc_mp_uninit(&r);
 315   silc_mp_uninit(&base);
 316   silc_mp_uninit(&tmp);
 317   silc_mp_uninit(&tmp2);
 318   silc_mp_uninit(&oprime);
 319
 320   return valid;
 321 }
 322
 323 /* Performs primality testings for given number. Returns TRUE if the
 324    number is probably a prime. */
 325
 326 SilcBool silc_math_prime_test(SilcMPInt *p)
 327 {
 328   SilcMPInt r, base, tmp;
 329   int i, ret = 0;
 330
 331   silc_mp_init(&r);
 332   silc_mp_init(&tmp);
 333   silc_mp_init(&base);
 334   silc_mp_set_ui(&base, 2);
 335
 336   SILC_LOG_DEBUG(("Testing probability of prime"));
 337
 338   /* See if the number is divisible by any of the
 339      small primes in primetable[]. */
 340   for (i = 0; primetable[i] != 0; i++) {
 341     silc_mp_mod_ui(&tmp, p, primetable[i]);
 342
 343     /* If mod is 0, the number is composite */
 344     if (silc_mp_cmp_ui(&tmp, 0) == 0) {
 345       SILC_LOG_DEBUG(("Number is not prime"));
 346       silc_mp_uninit(&r);
 347       silc_mp_uninit(&tmp);
 348       silc_mp_uninit(&base);
 349       return FALSE;
 350     }
 351   }
 352
 353   /* Does the prime pass the Fermat's prime test.
 354    * r = 2 ^ p mod p, if r == 2, then p is probably a prime.
 355    */
 356   silc_mp_pow_mod(&r, &base, p, p);
 357   if (silc_mp_cmp_ui(&r, 2) != 0)
 358     ret = -1;
 359
 360   silc_mp_uninit(&r);
 361   silc_mp_uninit(&tmp);
 362   silc_mp_uninit(&base);
 363
 364   if (ret) {
 365     SILC_LOG_DEBUG(("Number is not prime"));
 366     return FALSE;
 367   }
 368
 369   /* Number is probably a prime */
 370   return TRUE;
 371 }