/*
silcprimegen.c
Author: Pekka Riikonen <priikone@silcnet.org>
Copyright (C) 1997 - 2005 Pekka Riikonen
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
*/
/* Created: Mon Dec 8 16:35:37 GMT+0200 1997 */
/* $Id: silcprimegen.c,v 1.10.4.1 2005/04/24 15:24:12 priikone Exp $ */
#include "silcincludes.h"
/*
Fixed primetable for small prime division. We use this primetable to
test if possible prime is divisible any of these. Primetable is NULL
terminated. This same primetable can be found in SSH and PGP except that
SSH don't have prime 2 in the table. This sort of prime table is easily
created, for example, using sieve of Erastosthenes.
Just to include; if somebody is interested about Erastosthenes. Creating a
small prime table using sieve of Erastosthenes would go like this: Write
down a list of integers in range of 2 to n. Here in example n = 20.
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Then mark all multiples of 2 (Note: 2 is a prime number):
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
x x x x x x x x x
Move to the next unmarked number, 3, then mark all multiples of 3.
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
x x x x x x x x x x x
And continue doing this, marking all multiples of the next unmarked
number until there are now new unmarked numbers. And there you have a
prime table. So in this example, primes between 2 - 20 are:
2 3 5 7 11 13 17 19
*/
static SilcUInt32 primetable[] =
{
2, 3, 5, 7, 11, 13, 17, 19,
23, 29, 31, 37, 41, 43, 47, 53,
59, 61, 67, 71, 73, 79, 83, 89,
97, 101, 103, 107, 109, 113, 127, 131,
137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223,
227, 229, 233, 239, 241, 251, 257, 263,
269, 271, 277, 281, 283, 293, 307, 311,
313, 317, 331, 337, 347, 349, 353, 359,
367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457,
461, 463, 467, 479, 487, 491, 499, 503,
509, 521, 523, 541, 547, 557, 563, 569,
571, 577, 587, 593, 599, 601, 607, 613,
617, 619, 631, 641, 643, 647, 653, 659,
661, 673, 677, 683, 691, 701, 709, 719,
727, 733, 739, 743, 751, 757, 761, 769,
773, 787, 797, 809, 811, 821, 823, 827,
829, 839, 853, 857, 859, 863, 877, 881,
883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997,
1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049,
1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163,
1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283,
1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571,
1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619,
1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949,
1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069,
2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311,
2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377,
2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657,
2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801,
2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939,
2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167,
3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301,
3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347,
3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491,
3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541,
3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671,
3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797,
3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863,
3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003,
4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,
4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211,
4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259,
4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337,
4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481,
4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547,
4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621,
4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673,
4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813,
4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909,
4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967,
4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011,
5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167,
5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233,
5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309,
5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399,
5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507,
5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573,
5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653,
5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711,
5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849,
5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897,
5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007,
6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073,
6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211,
6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271,
6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329,
6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379,
6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563,
6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637,
6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701,
6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779,
6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907,
6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971,
6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027,
7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121,
7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253,
7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349,
7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457,
7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517,
7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621,
7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691,
7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757,
7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853,
7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009,
8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087,
8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161,
8167, 8171, 8179, 8191, 0
};
/* Find appropriate prime. It generates a number by taking random bytes.
It then tests the number that it's not divisible by any of the small
primes and then it performs Fermat's prime test. I thank Rieks Joosten
(r.joosten@pijnenburg.nl) for such a good help with prime tests.
If argument verbose is TRUE this will display some status information
about the progress of generation. */
bool silc_math_gen_prime(SilcMPInt *prime, SilcUInt32 bits, bool verbose,
SilcRng rng)
{
unsigned char *numbuf = NULL;
SilcUInt32 i, b, k;
SilcUInt32 *spmods;
SilcMPInt r, base, tmp, tmp2, oprime;
bool valid = FALSE;
silc_mp_init(&r);
silc_mp_init(&base);
silc_mp_init(&tmp);
silc_mp_init(&tmp2);
silc_mp_init(&oprime);
silc_mp_set_ui(&base, 2);
SILC_LOG_DEBUG(("Generating new prime"));
while (valid == FALSE) {
/* Get random number */
if (rng)
numbuf = silc_rng_get_rn_data(rng, (bits / 8));
else
numbuf = silc_rng_global_get_rn_data((bits / 8));
if (!numbuf)
return FALSE;
/* Convert into MP and set the size */
silc_mp_bin2mp(numbuf, (bits / 8), prime);
silc_mp_mod_2exp(prime, prime, bits);
/* Empty buffer */
memset(numbuf, 0, (bits / 8));
silc_free(numbuf);
/* Set highest bit */
silc_mp_set_ui(&tmp, 1);
silc_mp_mul_2exp(&tmp, &tmp, bits - 1);
silc_mp_or(prime, prime, &tmp);
/* Number could be even number, so we'll make it odd. */
silc_mp_set_ui(&tmp, 1);
silc_mp_or(prime, prime, &tmp);
/* Init modulo table with the prime candidate and the primes
in the primetable. */
spmods = silc_calloc(1, sizeof(primetable) * sizeof(SilcUInt32));
for (i = 0; primetable[i] != 0; i++) {
silc_mp_mod_ui(&tmp, prime, primetable[i]);
spmods[i] = silc_mp_get_ui(&tmp);
}
/* k is added by 2, this way we skip all even numbers (prime is odd). */
silc_mp_set(&oprime, prime);
for (k = 0;; k += 2) {
silc_mp_add_ui(&oprime, prime, k);
/* See if the prime has a divisor in primetable[].
* If it has then it's a composite. We add k to the
* original modulo value, k is added by 2 on every roll.
* This way we don't have to re-init the whole table if
* the number is composite.
*/
for (b = 0; b < i; b++) {
silc_mp_set_ui(&tmp2, spmods[b]);
silc_mp_add_ui(&tmp2, &tmp2, k);
silc_mp_mod_ui(&tmp, &tmp2, primetable[b]);
/* If mod is 0, the number is composite */
if (silc_mp_cmp_ui(&tmp, 0) == 0)
break;
}
if (b < i)
continue;
/* Passed the quick test. */
/* Does the prime pass the Fermat's prime test.
* r = 2 ^ p mod p, if r == 2, then p is probably a prime.
*/
silc_mp_pow_mod(&r, &base, &oprime, &oprime);
if (silc_mp_cmp_ui(&r, 2) != 0) {
if (verbose) {
printf(".");
fflush(stdout);
}
continue;
}
/* Passed the Fermat's test. And I don't think
* we have to do more tests. If anyone wants to do more
* tests, MP library has probability of prime test:
*
* if (silc_mp_probab_prime_p(prime, 25))
* break; found a probable prime
*
* However, this is very slow.
*/
/* XXX: this could be done if some argument, say strict_primes, is
TRUE when we are willing to spend more time on the prime test and
to get, perhaps, better primes. */
silc_mp_set(prime, &oprime);
break;
}
/* Check highest bit */
silc_mp_div_2exp(&tmp, prime, bits - 1);
if (silc_mp_get_ui(&tmp) == 1) {
valid = TRUE;
break;
}
}
silc_free(spmods);
silc_mp_uninit(&r);
silc_mp_uninit(&base);
silc_mp_uninit(&tmp);
silc_mp_uninit(&tmp2);
silc_mp_uninit(&oprime);
return valid;
}
/* Performs primality testings for given number. Returns TRUE if the
number is probably a prime. */
bool silc_math_prime_test(SilcMPInt *p)
{
SilcMPInt r, base, tmp;
int i, ret = 0;
silc_mp_init(&r);
silc_mp_init(&tmp);
silc_mp_init(&base);
silc_mp_set_ui(&base, 2);
SILC_LOG_DEBUG(("Testing probability of prime"));
/* See if the number is divisible by any of the
small primes in primetable[]. */
for (i = 0; primetable[i] != 0; i++) {
silc_mp_mod_ui(&tmp, p, primetable[i]);
/* If mod is 0, the number is composite */
if (silc_mp_cmp_ui(&tmp, 0) == 0)
ret = -1;
}
/* Does the prime pass the Fermat's prime test.
* r = 2 ^ p mod p, if r == 2, then p is probably a prime.
*/
silc_mp_pow_mod(&r, &base, p, p);
if (silc_mp_cmp_ui(&r, 2) != 0)
ret = -1;
silc_mp_uninit(&r);
silc_mp_uninit(&tmp);
silc_mp_uninit(&base);
if (ret)
return FALSE;
/* Number is probably a prime */
return TRUE;
}
syntax highlighted by Code2HTML, v. 0.9.1