lib/silcmath/silcprimegen.c

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