nettle  3.7.3
About: Nettle is a low-level cryptographic library.
  Fossies Dox: nettle-3.7.3.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

bignum-random-prime.c
Go to the documentation of this file.
1 /* bignum-random-prime.c
2 
3  Generation of random provable primes.
4 
5  Copyright (C) 2010, 2013 Niels Möller
6 
7  This file is part of GNU Nettle.
8 
9  GNU Nettle is free software: you can redistribute it and/or
10  modify it under the terms of either:
11 
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.
15 
16  or
17 
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.
21 
22  or both in parallel, as here.
23 
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.
28 
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/.
32 */
33 
34 #if HAVE_CONFIG_H
35 # include "config.h"
36 #endif
37 
38 #ifndef RANDOM_PRIME_VERBOSE
39 #define RANDOM_PRIME_VERBOSE 0
40 #endif
41 
42 #include <assert.h>
43 #include <stdlib.h>
44 
45 #if RANDOM_PRIME_VERBOSE
46 #include <stdio.h>
47 #define VERBOSE(x) (fputs((x), stderr))
48 #else
49 #define VERBOSE(x)
50 #endif
51 
52 #include "bignum.h"
53 #include "hogweed-internal.h"
54 #include "macros.h"
55 
56 /* Use a table of p_2 = 3 to p_{172} = 1021, used for sieving numbers
57  of up to 20 bits. */
58 
59 #define NPRIMES 171
60 #define TRIAL_DIV_BITS 20
61 #define TRIAL_DIV_MASK ((1 << TRIAL_DIV_BITS) - 1)
62 
63 /* A 20-bit number x is divisible by p iff
64 
65  ((x * inverse) & TRIAL_DIV_MASK) <= limit
66 */
68  uint32_t inverse; /* p^{-1} (mod 2^20) */
69  uint32_t limit; /* floor( (2^20 - 1) / p) */
70 };
71 
72 static const uint16_t
74  3,5,7,11,13,17,19,23,
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,
95  1013,1019,1021,
96 };
97 
98 static const uint32_t
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
122 };
123 
124 static const struct trial_div_info
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},
169 };
170 
171 /* Element j gives the index of the first prime of size 3+j bits */
172 static uint8_t
174  1,3,5,10,17,30,53,96,171
175 };
176 
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
179  prime. */
180 static int
182 {
183  mpz_t r;
184  mpz_t y;
185  int is_prime = 0;
186 
187  /* Avoid the mp_bitcnt_t type for compatibility with older GMP
188  versions. */
189  unsigned k;
190  unsigned j;
191 
192  VERBOSE(".");
193 
194  if (mpz_even_p(n) || mpz_cmp_ui(n, 3) < 0)
195  return 0;
196 
197  mpz_init(r);
198  mpz_init(y);
199 
200  k = mpz_scan1(nm1, 0);
201  assert(k > 0);
202 
203  mpz_fdiv_q_2exp (r, nm1, k);
204 
205  mpz_powm(y, a, r, n);
206 
207  if (mpz_cmp_ui(y, 1) == 0 || mpz_cmp(y, nm1) == 0)
208  goto passed_miller_rabin;
209 
210  for (j = 1; j < k; j++)
211  {
212  mpz_powm_ui (y, y, 2, n);
213 
214  if (mpz_cmp_ui (y, 1) == 0)
215  break;
216 
217  if (mpz_cmp (y, nm1) == 0)
218  {
219  passed_miller_rabin:
220  /* We know that a^{n-1} = 1 (mod n)
221 
222  Remains to check that gcd(a^{(n-1)/q} - 1, n) == 1 */
223  VERBOSE("x");
224 
225  mpz_powm(y, a, nm1dq, n);
226  mpz_sub_ui(y, y, 1);
227  mpz_gcd(y, y, n);
228  is_prime = mpz_cmp_ui (y, 1) == 0;
229  VERBOSE(is_prime ? "\n" : "");
230  break;
231  }
232 
233  }
234 
235  mpz_clear(r);
236  mpz_clear(y);
237 
238  return is_prime;
239 }
240 
241 /* The most basic variant of Pocklingtons theorem:
242 
243  Assume that q^e | (n-1), with q prime. If we can find an a such that
244 
245  a^{n-1} = 1 (mod n)
246  gcd(a^{(n-1)/q} - 1, n) = 1
247 
248  then any prime divisor p of n satisfies p = 1 (mod q^e).
249 
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
252  this by d.
253 
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.
258 
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).
261 
262 
263  * Variant, slightly stronger than Fact 4.59, HAC:
264 
265  Assume n = 1 + 2rq, q an odd prime, r <= 2q, and
266 
267  a^{n-1} = 1 (mod n)
268  gcd(a^{(n-1)/q} - 1, n) = 1
269 
270  Then n is prime.
271 
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.
276 
277  In bits, the requirement is that #n <= 2 #q, then
278 
279  r = (n-1)/2q < 2^{#n - #q} <= 2^#q = 2 2^{#q-1}< 2 q
280 
281 
282  * Another variant with an extra test (Variant of Fact 4.42, HAC):
283 
284  Assume n = 1 + 2rq, n odd, q an odd prime, 8 q^3 >= n
285 
286  a^{n-1} = 1 (mod n)
287  gcd(a^{(n-1)/q} - 1, n) = 1
288 
289  Also let x = floor(r / 2q), y = r mod 2q,
290 
291  If y^2 - 4x is not a square, then n is prime.
292 
293  Proof (adapted from Maurer, Journal of Cryptology, 8 (1995)):
294 
295  Assume n is composite. There are at most two factors, both odd,
296 
297  n = (1+2m_1 q)(1+2m_2 q) = 1 + 4 m_1 m_2 q^2 + 2 (m_1 + m_2) q
298 
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.
302 
303  We have the bound
304 
305  m_1 + m_2 < 2q / m_2 + m_2 <= 2q + 1 (maximum value for m_2 = 1)
306 
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.
309 
310  Next, write r = (n-1)/2q = 2 m_1 m_2 q + m_1 + m_2.
311 
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
314 
315  m^2 - y m + x = 0
316 
317  which has integer roots iff y^2 - 4 x is the square of an integer.
318 
319  In bits, the requirement is that #n <= 3 #q, then
320 
321  n < 2^#n <= 2^{3 #q} = 8 2^{3 (#q-1)} < 8 q^3
322 */
323 
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. */
329 void
331  unsigned bits, int top_bits_set,
332  void *ctx, nettle_random_func *random,
333  const mpz_t p0,
334  const mpz_t q,
335  const mpz_t p0q)
336 {
337  mpz_t r_min, r_range, pm1, a, e;
338  int need_square_test;
339  unsigned p0_bits;
340  mpz_t x, y, p04;
341 
342  p0_bits = mpz_sizeinbase (p0, 2);
343 
344  assert (bits <= 3*p0_bits);
345  assert (bits > p0_bits);
346 
347  need_square_test = (bits > 2 * p0_bits);
348 
349  mpz_init (r_min);
350  mpz_init (r_range);
351  mpz_init (pm1);
352  mpz_init (a);
353 
354  if (need_square_test)
355  {
356  mpz_init (x);
357  mpz_init (y);
358  mpz_init (p04);
359  mpz_mul_2exp (p04, p0, 2);
360  }
361 
362  if (q)
363  mpz_init (e);
364 
365  if (top_bits_set)
366  {
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);
375  }
376  else
377  {
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);
383  }
384 
385  for (;;)
386  {
387  uint8_t buf[1];
388 
389  nettle_mpz_random (r, ctx, random, r_range);
390  mpz_add (r, r, r_min);
391 
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);
396 
397  assert(mpz_sizeinbase(p, 2) == bits);
398 
399  /* Should use GMP trial division interface when that
400  materializes, we don't need any testing beyond trial
401  division. */
402  if (!mpz_probab_prime_p (p, 1))
403  continue;
404 
405  random(ctx, sizeof(buf), buf);
406 
407  mpz_set_ui (a, buf[0] + 2);
408 
409  if (q)
410  {
411  mpz_mul (e, r, q);
412  if (!miller_rabin_pocklington(p, pm1, e, a))
413  continue;
414 
415  if (need_square_test)
416  {
417  /* Our e corresponds to 2r in the theorem */
418  mpz_tdiv_qr (x, y, e, p04);
419  goto square_test;
420  }
421  }
422  else
423  {
424  if (!miller_rabin_pocklington(p, pm1, r, a))
425  continue;
426  if (need_square_test)
427  {
428  mpz_tdiv_qr (x, y, r, p04);
429  square_test:
430  /* We have r' = 2r, x = floor (r/2q) = floor(r'/2q),
431  and y' = r' - x 4q = 2 (r - x 2q) = 2y.
432 
433  Then y^2 - 4x is a square iff y'^2 - 16 x is a
434  square. */
435 
436  mpz_mul (y, y, y);
437  mpz_submul_ui (y, x, 16);
438  if (mpz_perfect_square_p (y))
439  continue;
440  }
441  }
442 
443  /* If we passed all the tests, we have found a prime. */
444  break;
445  }
446  mpz_clear (r_min);
447  mpz_clear (r_range);
448  mpz_clear (pm1);
449  mpz_clear (a);
450 
451  if (need_square_test)
452  {
453  mpz_clear (x);
454  mpz_clear (y);
455  mpz_clear (p04);
456  }
457  if (q)
458  mpz_clear (e);
459 }
460 
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). */
464 void
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)
468 {
469  assert (bits >= 3);
470  if (bits <= 10)
471  {
472  unsigned first;
473  unsigned choices;
474  uint8_t buf;
475 
476  assert (!top_bits_set);
477 
478  random (random_ctx, sizeof(buf), &buf);
479 
480  first = prime_by_size[bits-3];
481  choices = prime_by_size[bits-2] - first;
482 
483  mpz_set_ui (p, primes[first + buf % choices]);
484  }
485  else if (bits <= 20)
486  {
487  unsigned long highbit;
488  uint8_t buf[3];
489  unsigned long x;
490  unsigned j;
491 
492  assert (!top_bits_set);
493 
494  highbit = 1L << (bits - 1);
495 
496  again:
497  random (random_ctx, sizeof(buf), buf);
498  x = READ_UINT24(buf);
499  x &= (highbit - 1);
500  x |= highbit | 1;
501 
502  for (j = 0; prime_square[j] <= x; j++)
503  {
504  unsigned q = x * trial_div_table[j].inverse & TRIAL_DIV_MASK;
505  if (q <= trial_div_table[j].limit)
506  goto again;
507  }
508  mpz_set_ui (p, x);
509  }
510  else
511  {
512  mpz_t q, r;
513 
514  mpz_init (q);
515  mpz_init (r);
516 
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);
522 
523  _nettle_generate_pocklington_prime (p, r, bits, top_bits_set,
524  random_ctx, random,
525  q, NULL, q);
526 
527  if (progress)
528  progress (progress_ctx, 'x');
529 
530  mpz_clear (q);
531  mpz_clear (r);
532  }
533 }
static uint8_t prime_by_size[9]
static const uint32_t prime_square[171+1]
static const struct trial_div_info trial_div_table[171]
#define TRIAL_DIV_MASK
static const uint16_t primes[171]
void _nettle_generate_pocklington_prime(mpz_t p, mpz_t r, unsigned bits, int top_bits_set, void *ctx, nettle_random_func *random, const mpz_t p0, const mpz_t q, const mpz_t p0q)
void nettle_random_prime(mpz_t p, unsigned bits, int top_bits_set, void *random_ctx, nettle_random_func *random, void *progress_ctx, nettle_progress_func *progress)
static int miller_rabin_pocklington(mpz_t n, mpz_t nm1, mpz_t nm1dq, mpz_t a)
#define VERBOSE(x)
#define NPRIMES
void nettle_mpz_random(mpz_t x, void *ctx, nettle_random_func *random, const mpz_t n)
Definition: bignum-random.c:64
#define x
#define j
#define k
#define NULL
Definition: getopt1.c:58
#define READ_UINT24(p)
Definition: macros.h:74
int mpz_cmp_ui(const mpz_t u, unsigned long v)
Definition: mini-gmp.c:1866
int mpz_perfect_square_p(const mpz_t u)
Definition: mini-gmp.c:3271
void mpz_sub_ui(mpz_t r, const mpz_t a, unsigned long b)
Definition: mini-gmp.c:1947
void mpz_mul_ui(mpz_t r, const mpz_t u, unsigned long int v)
Definition: mini-gmp.c:2048
void mpz_add(mpz_t r, const mpz_t a, const mpz_t b)
Definition: mini-gmp.c:2008
void mpz_submul_ui(mpz_t r, const mpz_t u, unsigned long int v)
Definition: mini-gmp.c:2140
void mpz_fdiv_q_2exp(mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
Definition: mini-gmp.c:2477
void mpz_set_ui(mpz_t r, unsigned long int x)
Definition: mini-gmp.c:1483
size_t mpz_sizeinbase(const mpz_t u, int base)
Definition: mini-gmp.c:4167
void mpz_gcd(mpz_t g, const mpz_t u, const mpz_t v)
Definition: mini-gmp.c:2720
mp_bitcnt_t mpz_scan1(const mpz_t u, mp_bitcnt_t starting_bit)
Definition: mini-gmp.c:4099
void mpz_powm_ui(mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
Definition: mini-gmp.c:3173
void mpz_init(mpz_t r)
Definition: mini-gmp.c:1410
void mpz_fdiv_q(mpz_t q, const mpz_t n, const mpz_t d)
Definition: mini-gmp.c:2311
void mpz_powm(mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
Definition: mini-gmp.c:3068
void mpz_add_ui(mpz_t r, const mpz_t a, unsigned long b)
Definition: mini-gmp.c:1938
void mpz_mul(mpz_t r, const mpz_t u, const mpz_t v)
Definition: mini-gmp.c:2058
void mpz_tdiv_qr(mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
Definition: mini-gmp.c:2299
void mpz_clear(mpz_t r)
Definition: mini-gmp.c:1435
int mpz_probab_prime_p(const mpz_t n, int reps)
Definition: mini-gmp.c:3561
void mpz_mul_2exp(mpz_t r, const mpz_t u, mp_bitcnt_t bits)
Definition: mini-gmp.c:2096
int mpz_cmp(const mpz_t a, const mpz_t b)
Definition: mini-gmp.c:1877
__mpz_struct mpz_t[1]
Definition: mini-gmp.h:77
#define mpz_even_p(z)
Definition: mini-gmp.h:131
void nettle_random_func(void *ctx, size_t length, uint8_t *dst)
Definition: nettle-types.h:75
void nettle_progress_func(void *ctx, int c)
Definition: nettle-types.h:79