1 /* bignum-random-prime.c
3 Generation of random provable primes.
5 Copyright (C) 2010, 2013 Niels Möller
7 This file is part of GNU Nettle.
9 GNU Nettle is free software: you can redistribute it and/or
10 modify it under the terms of either:
12 * the GNU Lesser General Public License as published by the Free
13 Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
18 * the GNU General Public License as published by the Free
19 Software Foundation; either version 2 of the License, or (at your
20 option) any later version.
22 or both in parallel, as here.
24 GNU Nettle is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 General Public License for more details.
29 You should have received copies of the GNU General Public License and
30 the GNU Lesser General Public License along with this program. If
31 not, see http://www.gnu.org/licenses/.
38 #ifndef RANDOM_PRIME_VERBOSE
39 #define RANDOM_PRIME_VERBOSE 0
45 #if RANDOM_PRIME_VERBOSE
47 #define VERBOSE(x) (fputs((x), stderr))
53 #include "hogweed-internal.h"
56 /* Use a table of p_2 = 3 to p_{172} = 1021, used for sieving numbers
60 #define TRIAL_DIV_BITS 20
61 #define TRIAL_DIV_MASK ((1 << TRIAL_DIV_BITS) - 1)
63 /* A 20-bit number x is divisible by p iff
65 ((x * inverse) & TRIAL_DIV_MASK) <= limit
67 struct trial_div_info {
68 uint32_t inverse; /* p^{-1} (mod 2^20) */
69 uint32_t limit; /* floor( (2^20 - 1) / p) */
75 29,31,37,41,43,47,53,59,
76 61,67,71,73,79,83,89,97,
77 101,103,107,109,113,127,131,137,
78 139,149,151,157,163,167,173,179,
79 181,191,193,197,199,211,223,227,
80 229,233,239,241,251,257,263,269,
81 271,277,281,283,293,307,311,313,
82 317,331,337,347,349,353,359,367,
83 373,379,383,389,397,401,409,419,
84 421,431,433,439,443,449,457,461,
85 463,467,479,487,491,499,503,509,
86 521,523,541,547,557,563,569,571,
87 577,587,593,599,601,607,613,617,
88 619,631,641,643,647,653,659,661,
89 673,677,683,691,701,709,719,727,
90 733,739,743,751,757,761,769,773,
91 787,797,809,811,821,823,827,829,
92 839,853,857,859,863,877,881,883,
93 887,907,911,919,929,937,941,947,
94 953,967,971,977,983,991,997,1009,
99 prime_square[NPRIMES+1] = {
100 9,25,49,121,169,289,361,529,
101 841,961,1369,1681,1849,2209,2809,3481,
102 3721,4489,5041,5329,6241,6889,7921,9409,
103 10201,10609,11449,11881,12769,16129,17161,18769,
104 19321,22201,22801,24649,26569,27889,29929,32041,
105 32761,36481,37249,38809,39601,44521,49729,51529,
106 52441,54289,57121,58081,63001,66049,69169,72361,
107 73441,76729,78961,80089,85849,94249,96721,97969,
108 100489,109561,113569,120409,121801,124609,128881,134689,
109 139129,143641,146689,151321,157609,160801,167281,175561,
110 177241,185761,187489,192721,196249,201601,208849,212521,
111 214369,218089,229441,237169,241081,249001,253009,259081,
112 271441,273529,292681,299209,310249,316969,323761,326041,
113 332929,344569,351649,358801,361201,368449,375769,380689,
114 383161,398161,410881,413449,418609,426409,434281,436921,
115 452929,458329,466489,477481,491401,502681,516961,528529,
116 537289,546121,552049,564001,573049,579121,591361,597529,
117 619369,635209,654481,657721,674041,677329,683929,687241,
118 703921,727609,734449,737881,744769,769129,776161,779689,
119 786769,822649,829921,844561,863041,877969,885481,896809,
120 908209,935089,942841,954529,966289,982081,994009,1018081,
121 1026169,1038361,1042441,1L<<20
124 static const struct trial_div_info
125 trial_div_table[NPRIMES] = {
126 {699051,349525},{838861,209715},{748983,149796},{953251,95325},
127 {806597,80659},{61681,61680},{772635,55188},{866215,45590},
128 {180789,36157},{1014751,33825},{793517,28339},{1023001,25575},
129 {48771,24385},{870095,22310},{217629,19784},{710899,17772},
130 {825109,17189},{281707,15650},{502135,14768},{258553,14364},
131 {464559,13273},{934875,12633},{1001449,11781},{172961,10810},
132 {176493,10381},{203607,10180},{568387,9799},{788837,9619},
133 {770193,9279},{1032063,8256},{544299,8004},{619961,7653},
134 {550691,7543},{182973,7037},{229159,6944},{427445,6678},
135 {701195,6432},{370455,6278},{90917,6061},{175739,5857},
136 {585117,5793},{225087,5489},{298817,5433},{228877,5322},
137 {442615,5269},{546651,4969},{244511,4702},{83147,4619},
138 {769261,4578},{841561,4500},{732687,4387},{978961,4350},
139 {133683,4177},{65281,4080},{629943,3986},{374213,3898},
140 {708079,3869},{280125,3785},{641833,3731},{618771,3705},
141 {930477,3578},{778747,3415},{623751,3371},{40201,3350},
142 {122389,3307},{950371,3167},{1042353,3111},{18131,3021},
143 {285429,3004},{549537,2970},{166487,2920},{294287,2857},
144 {919261,2811},{636339,2766},{900735,2737},{118605,2695},
145 {10565,2641},{188273,2614},{115369,2563},{735755,2502},
146 {458285,2490},{914767,2432},{370513,2421},{1027079,2388},
147 {629619,2366},{462401,2335},{649337,2294},{316165,2274},
148 {484655,2264},{65115,2245},{326175,2189},{1016279,2153},
149 {990915,2135},{556859,2101},{462791,2084},{844629,2060},
150 {404537,2012},{457123,2004},{577589,1938},{638347,1916},
151 {892325,1882},{182523,1862},{1002505,1842},{624371,1836},
152 {69057,1817},{210787,1786},{558769,1768},{395623,1750},
153 {992745,1744},{317855,1727},{384877,1710},{372185,1699},
154 {105027,1693},{423751,1661},{408961,1635},{908331,1630},
155 {74551,1620},{36933,1605},{617371,1591},{506045,1586},
156 {24929,1558},{529709,1548},{1042435,1535},{31867,1517},
157 {166037,1495},{928781,1478},{508975,1458},{4327,1442},
158 {779637,1430},{742091,1418},{258263,1411},{879631,1396},
159 {72029,1385},{728905,1377},{589057,1363},{348621,1356},
160 {671515,1332},{710453,1315},{84249,1296},{959363,1292},
161 {685853,1277},{467591,1274},{646643,1267},{683029,1264},
162 {439927,1249},{254461,1229},{660713,1223},{554195,1220},
163 {202911,1215},{753253,1195},{941457,1190},{776635,1187},
164 {509511,1182},{986147,1156},{768879,1151},{699431,1140},
165 {696417,1128},{86169,1119},{808997,1114},{25467,1107},
166 {201353,1100},{708087,1084},{1018339,1079},{341297,1073},
167 {434151,1066},{96287,1058},{950765,1051},{298257,1039},
168 {675933,1035},{167731,1029},{815445,1027},
171 /* Element j gives the index of the first prime of size 3+j bits */
174 1,3,5,10,17,30,53,96,171
177 /* Combined Miller-Rabin test to the base a, and checking the
178 conditions from Pocklington's theorem, nm1dq holds (n-1)/q, with q
181 miller_rabin_pocklington(mpz_t n, mpz_t nm1, mpz_t nm1dq, mpz_t a)
187 /* Avoid the mp_bitcnt_t type for compatibility with older GMP
194 if (mpz_even_p(n) || mpz_cmp_ui(n, 3) < 0)
200 k = mpz_scan1(nm1, 0);
203 mpz_fdiv_q_2exp (r, nm1, k);
205 mpz_powm(y, a, r, n);
207 if (mpz_cmp_ui(y, 1) == 0 || mpz_cmp(y, nm1) == 0)
208 goto passed_miller_rabin;
210 for (j = 1; j < k; j++)
212 mpz_powm_ui (y, y, 2, n);
214 if (mpz_cmp_ui (y, 1) == 0)
217 if (mpz_cmp (y, nm1) == 0)
220 /* We know that a^{n-1} = 1 (mod n)
222 Remains to check that gcd(a^{(n-1)/q} - 1, n) == 1 */
225 mpz_powm(y, a, nm1dq, n);
228 is_prime = mpz_cmp_ui (y, 1) == 0;
229 VERBOSE(is_prime ? "\n" : "");
241 /* The most basic variant of Pocklingtons theorem:
243 Assume that q^e | (n-1), with q prime. If we can find an a such that
246 gcd(a^{(n-1)/q} - 1, n) = 1
248 then any prime divisor p of n satisfies p = 1 (mod q^e).
250 Proof (Cohen, 8.3.2): Assume p is a prime factor of n. The central
251 idea of the proof is to consider the order, modulo p, of a. Denote
254 a^{n-1} = 1 (mod n) implies a^{n-1} = 1 (mod p), hence d | (n-1).
255 Next, the condition gcd(a^{(n-1)/q} - 1, n) = 1 implies that
256 a^{(n-1)/q} != 1, hence d does not divide (n-1)/q. Since q is
257 prime, this means that q^e | d.
259 Finally, we have a^{p-1} = 1 (mod p), hence d | (p-1). So q^e | d |
260 (p-1), which gives the desired result: p = 1 (mod q^e).
263 * Variant, slightly stronger than Fact 4.59, HAC:
265 Assume n = 1 + 2rq, q an odd prime, r <= 2q, and
268 gcd(a^{(n-1)/q} - 1, n) = 1
272 Proof: By Pocklington's theorem, any prime factor p satisfies p = 1
273 (mod q). Neither 1 or q+1 are primes, hence p >= 1 + 2q. If n is
274 composite, we have n >= (1+2q)^2. But the assumption r <= 2q
275 implies n <= 1 + 4q^2, a contradiction.
277 In bits, the requirement is that #n <= 2 #q, then
279 r = (n-1)/2q < 2^{#n - #q} <= 2^#q = 2 2^{#q-1}< 2 q
282 * Another variant with an extra test (Variant of Fact 4.42, HAC):
284 Assume n = 1 + 2rq, n odd, q an odd prime, 8 q^3 >= n
287 gcd(a^{(n-1)/q} - 1, n) = 1
289 Also let x = floor(r / 2q), y = r mod 2q,
291 If y^2 - 4x is not a square, then n is prime.
293 Proof (adapted from Maurer, Journal of Cryptology, 8 (1995)):
295 Assume n is composite. There are at most two factors, both odd,
297 n = (1+2m_1 q)(1+2m_2 q) = 1 + 4 m_1 m_2 q^2 + 2 (m_1 + m_2) q
299 where we can assume m_1 >= m_2. Then the bound n <= 8 q^3 implies m_1
300 m_2 < 2q, restricting (m_1, m_2) to the domain 0 < m_2 <
301 sqrt(2q), 0 < m_1 < 2q / m_2.
305 m_1 + m_2 < 2q / m_2 + m_2 <= 2q + 1 (maximum value for m_2 = 1)
307 And the case m_1 = 2q, m_2 = 1 can be excluded, because it gives n
308 > 8q^3. So in fact, m_1 + m_2 < 2q.
310 Next, write r = (n-1)/2q = 2 m_1 m_2 q + m_1 + m_2.
312 If follows that m_1 + m_2 = y and m_1 m_2 = x. m_1 and m_2 are
313 thus the roots of the equation
317 which has integer roots iff y^2 - 4 x is the square of an integer.
319 In bits, the requirement is that #n <= 3 #q, then
321 n < 2^#n <= 2^{3 #q} = 8 2^{3 (#q-1)} < 8 q^3
324 /* Generate a prime number p of size bits with 2 p0q dividing (p-1).
325 p0 must be of size >= ceil(bits/3). The extra factor q can be
326 omitted (then p0 and p0q should be equal). If top_bits_set is one,
327 the topmost two bits are set to one, suitable for RSA primes. Also
328 returns r = (p-1)/p0q. */
330 _nettle_generate_pocklington_prime (mpz_t p, mpz_t r,
331 unsigned bits, int top_bits_set,
332 void *ctx, nettle_random_func *random,
337 mpz_t r_min, r_range, pm1, a, e;
338 int need_square_test;
342 p0_bits = mpz_sizeinbase (p0, 2);
344 assert (bits <= 3*p0_bits);
345 assert (bits > p0_bits);
347 need_square_test = (bits > 2 * p0_bits);
354 if (need_square_test)
359 mpz_mul_2exp (p04, p0, 2);
367 /* i = floor (2^{bits-3} / p0q), then 3I + 3 <= r <= 4I, with I
368 - 2 possible values. */
369 mpz_set_ui (r_min, 1);
370 mpz_mul_2exp (r_min, r_min, bits-3);
371 mpz_fdiv_q (r_min, r_min, p0q);
372 mpz_sub_ui (r_range, r_min, 2);
373 mpz_mul_ui (r_min, r_min, 3);
374 mpz_add_ui (r_min, r_min, 3);
378 /* i = floor (2^{bits-2} / p0q), I + 1 <= r <= 2I */
379 mpz_set_ui (r_range, 1);
380 mpz_mul_2exp (r_range, r_range, bits-2);
381 mpz_fdiv_q (r_range, r_range, p0q);
382 mpz_add_ui (r_min, r_range, 1);
389 nettle_mpz_random (r, ctx, random, r_range);
390 mpz_add (r, r, r_min);
392 /* Set p = 2*r*p0q + 1 */
393 mpz_mul_2exp(r, r, 1);
394 mpz_mul (pm1, r, p0q);
395 mpz_add_ui (p, pm1, 1);
397 assert(mpz_sizeinbase(p, 2) == bits);
399 /* Should use GMP trial division interface when that
400 materializes, we don't need any testing beyond trial
402 if (!mpz_probab_prime_p (p, 1))
405 random(ctx, sizeof(buf), buf);
407 mpz_set_ui (a, buf[0] + 2);
412 if (!miller_rabin_pocklington(p, pm1, e, a))
415 if (need_square_test)
417 /* Our e corresponds to 2r in the theorem */
418 mpz_tdiv_qr (x, y, e, p04);
424 if (!miller_rabin_pocklington(p, pm1, r, a))
426 if (need_square_test)
428 mpz_tdiv_qr (x, y, r, p04);
430 /* We have r' = 2r, x = floor (r/2q) = floor(r'/2q),
431 and y' = r' - x 4q = 2 (r - x 2q) = 2y.
433 Then y^2 - 4x is a square iff y'^2 - 16 x is a
437 mpz_submul_ui (y, x, 16);
438 if (mpz_perfect_square_p (y))
443 /* If we passed all the tests, we have found a prime. */
451 if (need_square_test)
461 /* Generate random prime of a given size. Maurer's algorithm (Alg.
462 6.42 Handbook of applied cryptography), but with ratio = 1/2 (like
463 the variant in fips186-3). */
465 nettle_random_prime(mpz_t p, unsigned bits, int top_bits_set,
466 void *random_ctx, nettle_random_func *random,
467 void *progress_ctx, nettle_progress_func *progress)
476 assert (!top_bits_set);
478 random (random_ctx, sizeof(buf), &buf);
480 first = prime_by_size[bits-3];
481 choices = prime_by_size[bits-2] - first;
483 mpz_set_ui (p, primes[first + buf % choices]);
487 unsigned long highbit;
492 assert (!top_bits_set);
494 highbit = 1L << (bits - 1);
497 random (random_ctx, sizeof(buf), buf);
498 x = READ_UINT24(buf);
502 for (j = 0; prime_square[j] <= x; j++)
504 unsigned q = x * trial_div_table[j].inverse & TRIAL_DIV_MASK;
505 if (q <= trial_div_table[j].limit)
517 /* Bit size ceil(k/2) + 1, slightly larger than used in Alg. 4.62
518 in Handbook of Applied Cryptography (which seems to be
519 incorrect for odd k). */
520 nettle_random_prime (q, (bits+3)/2, 0, random_ctx, random,
521 progress_ctx, progress);
523 _nettle_generate_pocklington_prime (p, r, bits, top_bits_set,
528 progress (progress_ctx, 'x');