bignum-random-prime.c

   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 */
  67 struct trial_div_info {
  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
  73 primes[NPRIMES] = {
  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
  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
 122 };
 123
 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},
 169 };
 170
 171 /* Element j gives the index of the first prime of size 3+j bits */
 172 static uint8_t
 173 prime_by_size[9] = {
 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
 181 miller_rabin_pocklington(mpz_t n, mpz_t nm1, mpz_t nm1dq, mpz_t a)
 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
 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,
 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 }