From 6ce825c75b43b58febee7c70f45509fde389d9d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niels=20M=C3=B6ller?= <nisse@lysator.liu.se>
Date: Wed, 26 May 2010 16:03:33 +0200
Subject: [PATCH] (_nettle_generate_pocklington_prime): New function. Rely on
mpz_probab_prime_p (for lack of a trial division function) for trial
division. (nettle_random_prime): Rewritten. Uses the prime table for the
smallest sizes, then trial division using a new set of tables, and then
Maurer's algorithm, calling the new _nettle_generate_pocklington_prime for
the final search.
Rev: nettle/ChangeLog:1.80
Rev: nettle/bignum-random-prime.c:1.4
Rev: nettle/bignum.h:1.5
---
ChangeLog | 12 +-
bignum-random-prime.c | 361 +++++++++++++++++++++++++-----------------
bignum.h | 7 +
3 files changed, 230 insertions(+), 150 deletions(-)
diff --git a/ChangeLog b/ChangeLog
index 7009c07f..4ce046e9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,7 +1,17 @@
+2010-05-26 Niels M�ller <nisse@lysator.liu.se>
+
+ * bignum-random-prime.c (_nettle_generate_pocklington_prime): New
+ function. Rely on mpz_probab_prime_p (for lack of a trial division
+ function) for trial division.
+ (nettle_random_prime): Rewritten. Uses the prime table for the
+ smallest sizes, then trial division using a new set of tables, and
+ then Maurer's algorithm, calling the new
+ _nettle_generate_pocklington_prime for the final search.
+
2010-05-25 Niels M�ller <nisse@lysator.liu.se>
* testsuite/dsa-test.c (test_main): Updated for dsa testing
- changes.
+ changes.
* testsuite/dsa-keygen-test.c (test_main): Test dsa256.
diff --git a/bignum-random-prime.c b/bignum-random-prime.c
index 601590ae..1fa7ee4a 100644
--- a/bignum-random-prime.c
+++ b/bignum-random-prime.c
@@ -45,73 +45,126 @@
#include "macros.h"
-/* Use a table of p_2 = 3 to p_{172} = 1021, multiplied to 32-bit or
- 64-bit size. */
-
-struct sieve_element {
- /* Product of some small primes. */
- unsigned long prod;
- /* Square of the smallest one. */
- unsigned long p2;
+/* Use a table of p_2 = 3 to p_{172} = 1021, used for sieving numbers
+ of up to 20 bits. */
+
+#define NPRIMES 171
+#define TRIAL_DIV_BITS 20
+#define TRIAL_DIV_MASK ((1 << TRIAL_DIV_BITS) - 1)
+
+/* A 20-bit number x is divisible by p iff
+
+ ((x * inverse) & TRIAL_DIV_MASK) <= limit
+*/
+struct trial_div_info {
+ uint32_t inverse; /* p^{-1} (mod 2^20) */
+ uint32_t limit; /* floor( (2^20 - 1) / p) */
+};
+
+static const uint16_t
+primes[NPRIMES] = {
+ 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,
+};
+
+static const uint32_t
+prime_square[NPRIMES+1] = {
+ 9,25,49,121,169,289,361,529,
+ 841,961,1369,1681,1849,2209,2809,3481,
+ 3721,4489,5041,5329,6241,6889,7921,9409,
+ 10201,10609,11449,11881,12769,16129,17161,18769,
+ 19321,22201,22801,24649,26569,27889,29929,32041,
+ 32761,36481,37249,38809,39601,44521,49729,51529,
+ 52441,54289,57121,58081,63001,66049,69169,72361,
+ 73441,76729,78961,80089,85849,94249,96721,97969,
+ 100489,109561,113569,120409,121801,124609,128881,134689,
+ 139129,143641,146689,151321,157609,160801,167281,175561,
+ 177241,185761,187489,192721,196249,201601,208849,212521,
+ 214369,218089,229441,237169,241081,249001,253009,259081,
+ 271441,273529,292681,299209,310249,316969,323761,326041,
+ 332929,344569,351649,358801,361201,368449,375769,380689,
+ 383161,398161,410881,413449,418609,426409,434281,436921,
+ 452929,458329,466489,477481,491401,502681,516961,528529,
+ 537289,546121,552049,564001,573049,579121,591361,597529,
+ 619369,635209,654481,657721,674041,677329,683929,687241,
+ 703921,727609,734449,737881,744769,769129,776161,779689,
+ 786769,822649,829921,844561,863041,877969,885481,896809,
+ 908209,935089,942841,954529,966289,982081,994009,1018081,
+ 1026169,1038361,1042441,1L<<20
};
-static const struct sieve_element
-sieve_table[] = {
- {111546435, 9}, /* 3 -- 23 */
- {58642669, 841}, /* 29 -- 43 */
- {600662303, 2209}, /* 47 -- 67 */
- {33984931, 5041}, /* 71 -- 83 */
- {89809099, 7921}, /* 89 -- 103 */
- {167375713, 11449}, /* 107 -- 127 */
- {371700317, 17161}, /* 131 -- 149 */
- {645328247, 22801}, /* 151 -- 167 */
- {1070560157, 29929}, /* 173 -- 191 */
- {1596463769, 37249}, /* 193 -- 211 */
- {11592209, 49729}, /* 223 -- 229 */
- {13420567, 54289}, /* 233 -- 241 */
- {16965341, 63001}, /* 251 -- 263 */
- {20193023, 72361}, /* 269 -- 277 */
- {23300239, 78961}, /* 281 -- 293 */
- {29884301, 94249}, /* 307 -- 313 */
- {35360399, 100489}, /* 317 -- 337 */
- {42749359, 120409}, /* 347 -- 353 */
- {49143869, 128881}, /* 359 -- 373 */
- {56466073, 143641}, /* 379 -- 389 */
- {65111573, 157609}, /* 397 -- 409 */
- {76027969, 175561}, /* 419 -- 431 */
- {84208541, 187489}, /* 433 -- 443 */
- {94593973, 201601}, /* 449 -- 461 */
- {103569859, 214369}, /* 463 -- 479 */
- {119319383, 237169}, /* 487 -- 499 */
- {133390067, 253009}, /* 503 -- 521 */
- {154769821, 273529}, /* 523 -- 547 */
- {178433279, 310249}, /* 557 -- 569 */
- {193397129, 326041}, /* 571 -- 587 */
- {213479407, 351649}, /* 593 -- 601 */
- {229580147, 368449}, /* 607 -- 617 */
- {250367549, 383161}, /* 619 -- 641 */
- {271661713, 413449}, /* 643 -- 653 */
- {293158127, 434281}, /* 659 -- 673 */
- {319512181, 458329}, /* 677 -- 691 */
- {357349471, 491401}, /* 701 -- 719 */
- {393806449, 528529}, /* 727 -- 739 */
- {422400701, 552049}, /* 743 -- 757 */
- {452366557, 579121}, /* 761 -- 773 */
- {507436351, 619369}, /* 787 -- 809 */
- {547978913, 657721}, /* 811 -- 823 */
- {575204137, 683929}, /* 827 -- 839 */
- {627947039, 727609}, /* 853 -- 859 */
- {666785731, 744769}, /* 863 -- 881 */
- {710381447, 779689}, /* 883 -- 907 */
- {777767161, 829921}, /* 911 -- 929 */
- {834985999, 877969}, /* 937 -- 947 */
- {894826021, 908209}, /* 953 -- 971 */
- {951747481, 954529}, /* 977 -- 991 */
- {1019050649, 994009}, /* 997 -- 1013 */
- {1040399, 1038361}, /* 1019 -- 1021 */
+static const struct trial_div_info
+trial_div_table[NPRIMES] = {
+ {699051,349525},{838861,209715},{748983,149796},{953251,95325},
+ {806597,80659},{61681,61680},{772635,55188},{866215,45590},
+ {180789,36157},{1014751,33825},{793517,28339},{1023001,25575},
+ {48771,24385},{870095,22310},{217629,19784},{710899,17772},
+ {825109,17189},{281707,15650},{502135,14768},{258553,14364},
+ {464559,13273},{934875,12633},{1001449,11781},{172961,10810},
+ {176493,10381},{203607,10180},{568387,9799},{788837,9619},
+ {770193,9279},{1032063,8256},{544299,8004},{619961,7653},
+ {550691,7543},{182973,7037},{229159,6944},{427445,6678},
+ {701195,6432},{370455,6278},{90917,6061},{175739,5857},
+ {585117,5793},{225087,5489},{298817,5433},{228877,5322},
+ {442615,5269},{546651,4969},{244511,4702},{83147,4619},
+ {769261,4578},{841561,4500},{732687,4387},{978961,4350},
+ {133683,4177},{65281,4080},{629943,3986},{374213,3898},
+ {708079,3869},{280125,3785},{641833,3731},{618771,3705},
+ {930477,3578},{778747,3415},{623751,3371},{40201,3350},
+ {122389,3307},{950371,3167},{1042353,3111},{18131,3021},
+ {285429,3004},{549537,2970},{166487,2920},{294287,2857},
+ {919261,2811},{636339,2766},{900735,2737},{118605,2695},
+ {10565,2641},{188273,2614},{115369,2563},{735755,2502},
+ {458285,2490},{914767,2432},{370513,2421},{1027079,2388},
+ {629619,2366},{462401,2335},{649337,2294},{316165,2274},
+ {484655,2264},{65115,2245},{326175,2189},{1016279,2153},
+ {990915,2135},{556859,2101},{462791,2084},{844629,2060},
+ {404537,2012},{457123,2004},{577589,1938},{638347,1916},
+ {892325,1882},{182523,1862},{1002505,1842},{624371,1836},
+ {69057,1817},{210787,1786},{558769,1768},{395623,1750},
+ {992745,1744},{317855,1727},{384877,1710},{372185,1699},
+ {105027,1693},{423751,1661},{408961,1635},{908331,1630},
+ {74551,1620},{36933,1605},{617371,1591},{506045,1586},
+ {24929,1558},{529709,1548},{1042435,1535},{31867,1517},
+ {166037,1495},{928781,1478},{508975,1458},{4327,1442},
+ {779637,1430},{742091,1418},{258263,1411},{879631,1396},
+ {72029,1385},{728905,1377},{589057,1363},{348621,1356},
+ {671515,1332},{710453,1315},{84249,1296},{959363,1292},
+ {685853,1277},{467591,1274},{646643,1267},{683029,1264},
+ {439927,1249},{254461,1229},{660713,1223},{554195,1220},
+ {202911,1215},{753253,1195},{941457,1190},{776635,1187},
+ {509511,1182},{986147,1156},{768879,1151},{699431,1140},
+ {696417,1128},{86169,1119},{808997,1114},{25467,1107},
+ {201353,1100},{708087,1084},{1018339,1079},{341297,1073},
+ {434151,1066},{96287,1058},{950765,1051},{298257,1039},
+ {675933,1035},{167731,1029},{815445,1027},
};
-#define SIEVE_SIZE (sizeof(sieve_table) / sizeof(sieve_table[0]))
+/* Element j gives the index of the first prime of size 3+j bits */
+static uint8_t
+prime_by_size[9] = {
+ 1,3,5,10,17,30,53,96,171
+};
/* Combined Miller-Rabin test to the base a, and checking the
conditions from Pocklington's theorem. */
@@ -176,19 +229,14 @@ miller_rabin_pocklington(mpz_t n, mpz_t nm1, mpz_t nm1dq, mpz_t a)
return is_prime;
}
-/* Generate random prime of a given size. Maurer's algorithm (Alg.
- 6.42 Handbook of applied cryptography), but with ratio = 1/2 (like
- the variant in fips186-3). FIXME: Force primes to start with two
- one bits? */
-
/* The algorithm is based on the following special case of
Pocklington's theorem:
- Assume that n = 1 + r q, where q is a prime, q > sqrt(n) - 1. If we
+ Assume that n = 1 + f q, where q is a prime, q > sqrt(n) - 1. If we
can find an a such that
a^{n-1} = 1 (mod n)
- gcd(a^r - 1, n) = 1
+ gcd(a^f - 1, n) = 1
then n is prime.
@@ -203,42 +251,98 @@ miller_rabin_pocklington(mpz_t n, mpz_t nm1, mpz_t nm1dq, mpz_t a)
If n is specified as k bits, we need q of size ceil(k/2) + 1 bits
(or more) to make the theorem apply.
- */
+*/
+
+/* Generate a prime number p of size bits with 2 p0q dividing (p-1).
+ p0 must be of size >= ceil(bits/2) + 1. The extra factor q can be
+ omitted. */
void
-nettle_random_prime(mpz_t p, unsigned bits,
- void *ctx, nettle_random_func random)
+_nettle_generate_pocklington_prime (mpz_t p, unsigned bits,
+ void *ctx, nettle_random_func random,
+ const mpz_t p0,
+ const mpz_t q,
+ const mpz_t p0q)
{
- assert (bits >= 6);
- if (bits < 20)
+ mpz_t i, r, pm1,a;
+
+ assert (2*mpz_sizeinbase (p0, 2) > bits + 1);
+
+ mpz_init (i);
+ mpz_init (r);
+ mpz_init (pm1);
+ mpz_init (a);
+
+ /* i = floor (2^{bits-2} / p0q) */
+ mpz_init_set_ui (i, 1);
+ mpz_mul_2exp (i, i, bits-2);
+ mpz_fdiv_q (i, i, p0q);
+
+ for (;;)
{
- unsigned long highbit;
- uint8_t buf[3];
- unsigned long x;
- unsigned j;
+ uint8_t buf[1];
- /* Small cases:
+ /* Generate r in the range i + 1 <= r <= 2*i */
+ nettle_mpz_random (r, ctx, random, i);
+ mpz_add (r, r, i);
+ mpz_add_ui (r, r, 1);
- 3 bits: 5 or 7
- 4 bits: 11, 13, 15
- 5 bits: 17, 19, 23, 29, 31
+ /* Set p = 2*r*p0q + 1 */
+ mpz_mul_2exp(r, r, 1);
+ mpz_mul (pm1, r, p0q);
+ mpz_add_ui (p, pm1, 1);
- With 3 bits, no sieving is done, since candidates are smaller
- than 3^2 = 9 (and this is ok; all odd 3-bit numbers are
- prime).
+ assert(mpz_sizeinbase(p, 2) == bits);
+
+ /* Should use GMP trial division interface when that
+ materializes, we don't need any testing beyond trial
+ division. */
+ if (!mpz_probab_prime_p (p, 1))
+ continue;
+
+ random(ctx, sizeof(buf), buf);
+
+ mpz_set_ui (a, buf[0] + 2);
+
+ if (q)
+ mpz_mul (r, r, q);
+
+ if (miller_rabin_pocklington(p, pm1, r, a))
+ break;
+ }
+ mpz_clear (i);
+ mpz_clear (r);
+ mpz_clear (pm1);
+ mpz_clear (a);
+}
- With 4 bits, sieving with the first value, 3*5*...*23 doesn't
- work, since this includes the primes 11 and 13 in the
- interval. Of the odd numbers in the interval, 9, 11, 13, 15,
- only the factors of three need be discarded.
+/* Generate random prime of a given size. Maurer's algorithm (Alg.
+ 6.42 Handbook of applied cryptography), but with ratio = 1/2 (like
+ the variant in fips186-3). FIXME: Force primes to start with two
+ one bits? */
+void
+nettle_random_prime(mpz_t p, unsigned bits,
+ void *ctx, nettle_random_func random)
+{
+ assert (bits >= 3);
+ if (bits <= 10)
+ {
+ unsigned first;
+ unsigned choices;
+ uint8_t buf;
- With 5 bits, we still sieve with only the first value, which
- includes three of the primes in the interval. Of the odd
- numbers in the interval, 17, 19, (21), 23, (25), (27), 29,
- 31, we need to discard multiples of 3 and 5 only.
+ random (ctx, sizeof(buf), &buf);
- With 6 bits, we sieve with only the first value (since 63 <
- 29^2), and there's no problem.
- */
+ first = prime_by_size[bits-3];
+ choices = prime_by_size[bits-2] - first;
+
+ mpz_set_ui (p, primes[first + buf % choices]);
+ }
+ else if (bits <= 20)
+ {
+ unsigned long highbit;
+ uint8_t buf[3];
+ unsigned long x;
+ unsigned j;
highbit = 1L << (bits - 1);
@@ -248,69 +352,28 @@ nettle_random_prime(mpz_t p, unsigned bits,
x &= (highbit - 1);
x |= highbit | 1;
- mpz_set_ui (p, x);
- for (j = 0; j < SIEVE_SIZE && x >= sieve_table[j].p2; j++)
- if (mpz_gcd_ui (NULL, p, sieve_table[j].prod) != 1)
- goto again;
+ for (j = 0; prime_square[j] <= x; j++)
+ {
+ unsigned q = x * trial_div_table[j].inverse & TRIAL_DIV_MASK;
+ if (q <= trial_div_table[j].limit)
+ goto again;
+ }
+ mpz_set_ui (p, x);
}
else
{
- mpz_t q, r, nm1, t, a, i;
- unsigned j;
+ mpz_t q;
mpz_init (q);
- mpz_init (r);
- mpz_init (nm1);
- mpz_init (t);
- mpz_init (a);
- mpz_init (i);
/* Bit size ceil(k/2) + 1, slightly larger than used in Alg. 4.62
in Handbook of Applied Cryptography (which seems to be
incorrect for odd k). */
nettle_random_prime (q, (bits+3)/2, ctx, random);
- /* i = floor (2^{bits-2} / q) */
- mpz_init_set_ui (i, 1);
- mpz_mul_2exp (i, i, bits-2);
- mpz_fdiv_q (i, i, q);
-
- for (;;)
- {
- uint8_t buf[1];
-
- /* Generate r in the range i + 1 <= r <= 2*i */
- nettle_mpz_random (r, ctx, random, i);
- mpz_add (r, r, i);
- mpz_add_ui (r, r, 1);
-
- /* Set p = 2*r*q + 1 */
- mpz_mul_2exp(r, r, 1);
- mpz_mul (nm1, r, q);
- mpz_add_ui (p, nm1, 1);
-
- assert(mpz_sizeinbase(p, 2) == bits);
-
- for (j = 0; j < SIEVE_SIZE; j++)
- {
- if (mpz_gcd_ui (NULL, p, sieve_table[j].prod) != 1)
- goto composite;
- }
-
- random(ctx, sizeof(buf), buf);
-
- mpz_set_ui (a, buf[0] + 2);
-
- if (miller_rabin_pocklington(p, nm1, r, a))
- break;
- composite:
- ;
- }
+ _nettle_generate_pocklington_prime (p, bits, ctx, random,
+ q, NULL, q);
+
mpz_clear (q);
- mpz_clear (r);
- mpz_clear (nm1);
- mpz_clear (t);
- mpz_clear (a);
- mpz_clear (i);
}
}
diff --git a/bignum.h b/bignum.h
index 759e7618..9db38825 100644
--- a/bignum.h
+++ b/bignum.h
@@ -89,6 +89,13 @@ void
nettle_random_prime(mpz_t p, unsigned bits,
void *ctx, nettle_random_func random);
+void
+_nettle_generate_pocklington_prime (mpz_t p, unsigned bits,
+ void *ctx, nettle_random_func random,
+ const mpz_t p0,
+ const mpz_t q,
+ const mpz_t p0q);
+
/* sexp parsing */
struct sexp_iterator;
--
GitLab