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                          SilcRng rng)
 200 {
 201   unsigned char *numbuf = NULL;
 202   SilcUInt32 i, b, k;
 203   SilcUInt32 *spmods;
 204   SilcMPInt r, base, tmp, tmp2, oprime;
 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   /* Get random number and assure that the first digit is not zero since
 217      our conversion routines does not like the first digit being zero. */
 218   do {
 219     if (numbuf) {
 220       memset(numbuf, 0, (bits / 8));
 221       silc_free(numbuf);
 222     }
 223     if (rng)
 224       numbuf = silc_rng_get_rn_string(rng, (bits / 8));
 225     else
 226       numbuf = silc_rng_global_get_rn_string((bits / 8));
 227     if (!numbuf)
 228       return FALSE;
 229   } while (numbuf[0] == '0');
 230
 231   /* Convert into MP and set the size */
 232   silc_mp_set_str(prime, numbuf, 16);
 233   silc_mp_mod_2exp(prime, prime, bits);
 234
 235   /* Empty buffer */
 236   memset(numbuf, 0, (bits / 8));
 237   silc_free(numbuf);
 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);               /* OR operator */
 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   silc_free(spmods);
 305   silc_mp_uninit(&r);
 306   silc_mp_uninit(&base);
 307   silc_mp_uninit(&tmp);
 308   silc_mp_uninit(&tmp2);
 309   silc_mp_uninit(&oprime);
 310
 311   return TRUE;
 312 }
 313
 314 /* Performs primality testings for given number. Returns TRUE if the
 315    number is probably a prime. */
 316
 317 bool silc_math_prime_test(SilcMPInt *p)
 318 {
 319   SilcMPInt r, base, tmp;
 320   int i, ret = 0;
 321
 322   silc_mp_init(&r);
 323   silc_mp_init(&tmp);
 324   silc_mp_init(&base);
 325   silc_mp_set_ui(&base, 2);
 326
 327   SILC_LOG_DEBUG(("Testing probability of prime"));
 328
 329   /* See if the number is divisible by any of the
 330      small primes in primetable[]. */
 331   for (i = 0; primetable[i] != 0; i++) {
 332     silc_mp_mod_ui(&tmp, p, primetable[i]);
 333
 334     /* If mod is 0, the number is composite */
 335     if (silc_mp_cmp_ui(&tmp, 0) == 0)
 336       ret = -1;
 337   }
 338
 339   /* Does the prime pass the Fermat's prime test.
 340    * r = 2 ^ p mod p, if r == 2, then p is probably a prime.
 341    */
 342   silc_mp_pow_mod(&r, &base, p, p);
 343   if (silc_mp_cmp_ui(&r, 2) != 0)
 344     ret = -1;
 345
 346   silc_mp_uninit(&r);
 347   silc_mp_uninit(&tmp);
 348   silc_mp_uninit(&base);
 349
 350   if (ret)
 351     return FALSE;
 352
 353   /* Number is probably a prime */
 354   return TRUE;
 355 }