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