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