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