bignum-random-prime.c

   1 /* bignum-random-prime.c
   2  *
   3  * Generation of random provable primes.
   4  */
   5
   6 /* nettle, low-level cryptographics library
   7  *
   8  * Copyright (C) 2010 Niels Möller
   9  *
  10  * The nettle library is free software; you can redistribute it and/or modify
  11  * it under the terms of the GNU Lesser General Public License as published by
  12  * the Free Software Foundation; either version 2.1 of the License, or (at your
  13  * option) any later version.
  14  *
  15  * The nettle library is distributed in the hope that it will be useful, but
  16  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  17  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  18  * License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public License
  21  * along with the nettle library; see the file COPYING.LIB.  If not, write to
  22  * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  23  * MA 02111-1301, USA.
  24  */
  25
  26 #if HAVE_CONFIG_H
  27 # include "config.h"
  28 #endif
  29
  30 #ifndef RANDOM_PRIME_VERBOSE
  31 #define RANDOM_PRIME_VERBOSE 0
  32 #endif
  33
  34 #include <assert.h>
  35 #include <stdlib.h>
  36
  37 #if RANDOM_PRIME_VERBOSE
  38 #include <stdio.h>
  39 #define VERBOSE(x) (fputs((x), stderr))
  40 #else
  41 #define VERBOSE(x)
  42 #endif
  43
  44 #include "bignum.h"
  45
  46 #include "macros.h"
  47
  48 /* Use a table of p_2 = 3 to p_{172} = 1021, used for sieving numbers
  49    of up to 20 bits. */
  50
  51 #define NPRIMES 171
  52 #define TRIAL_DIV_BITS 20
  53 #define TRIAL_DIV_MASK ((1 << TRIAL_DIV_BITS) - 1)
  54
  55 /* A 20-bit number x is divisible by p iff
  56
  57      ((x * inverse) & TRIAL_DIV_MASK) <= limit
  58 */
  59 struct trial_div_info {
  60   uint32_t inverse; /* p^{-1} (mod 2^20) */
  61   uint32_t limit;   /* floor( (2^20 - 1) / p) */
  62 };
  63
  64 static const uint16_t
  65 primes[NPRIMES] = {
  66   3,5,7,11,13,17,19,23,
  67   29,31,37,41,43,47,53,59,
  68   61,67,71,73,79,83,89,97,
  69   101,103,107,109,113,127,131,137,
  70   139,149,151,157,163,167,173,179,
  71   181,191,193,197,199,211,223,227,
  72   229,233,239,241,251,257,263,269,
  73   271,277,281,283,293,307,311,313,
  74   317,331,337,347,349,353,359,367,
  75   373,379,383,389,397,401,409,419,
  76   421,431,433,439,443,449,457,461,
  77   463,467,479,487,491,499,503,509,
  78   521,523,541,547,557,563,569,571,
  79   577,587,593,599,601,607,613,617,
  80   619,631,641,643,647,653,659,661,
  81   673,677,683,691,701,709,719,727,
  82   733,739,743,751,757,761,769,773,
  83   787,797,809,811,821,823,827,829,
  84   839,853,857,859,863,877,881,883,
  85   887,907,911,919,929,937,941,947,
  86   953,967,971,977,983,991,997,1009,
  87   1013,1019,1021,
  88 };
  89
  90 static const uint32_t
  91 prime_square[NPRIMES+1] = {
  92   9,25,49,121,169,289,361,529,
  93   841,961,1369,1681,1849,2209,2809,3481,
  94   3721,4489,5041,5329,6241,6889,7921,9409,
  95   10201,10609,11449,11881,12769,16129,17161,18769,
  96   19321,22201,22801,24649,26569,27889,29929,32041,
  97   32761,36481,37249,38809,39601,44521,49729,51529,
  98   52441,54289,57121,58081,63001,66049,69169,72361,
  99   73441,76729,78961,80089,85849,94249,96721,97969,
 100   100489,109561,113569,120409,121801,124609,128881,134689,
 101   139129,143641,146689,151321,157609,160801,167281,175561,
 102   177241,185761,187489,192721,196249,201601,208849,212521,
 103   214369,218089,229441,237169,241081,249001,253009,259081,
 104   271441,273529,292681,299209,310249,316969,323761,326041,
 105   332929,344569,351649,358801,361201,368449,375769,380689,
 106   383161,398161,410881,413449,418609,426409,434281,436921,
 107   452929,458329,466489,477481,491401,502681,516961,528529,
 108   537289,546121,552049,564001,573049,579121,591361,597529,
 109   619369,635209,654481,657721,674041,677329,683929,687241,
 110   703921,727609,734449,737881,744769,769129,776161,779689,
 111   786769,822649,829921,844561,863041,877969,885481,896809,
 112   908209,935089,942841,954529,966289,982081,994009,1018081,
 113   1026169,1038361,1042441,1L<<20
 114 };
 115
 116 static const struct trial_div_info
 117 trial_div_table[NPRIMES] = {
 118   {699051,349525},{838861,209715},{748983,149796},{953251,95325},
 119   {806597,80659},{61681,61680},{772635,55188},{866215,45590},
 120   {180789,36157},{1014751,33825},{793517,28339},{1023001,25575},
 121   {48771,24385},{870095,22310},{217629,19784},{710899,17772},
 122   {825109,17189},{281707,15650},{502135,14768},{258553,14364},
 123   {464559,13273},{934875,12633},{1001449,11781},{172961,10810},
 124   {176493,10381},{203607,10180},{568387,9799},{788837,9619},
 125   {770193,9279},{1032063,8256},{544299,8004},{619961,7653},
 126   {550691,7543},{182973,7037},{229159,6944},{427445,6678},
 127   {701195,6432},{370455,6278},{90917,6061},{175739,5857},
 128   {585117,5793},{225087,5489},{298817,5433},{228877,5322},
 129   {442615,5269},{546651,4969},{244511,4702},{83147,4619},
 130   {769261,4578},{841561,4500},{732687,4387},{978961,4350},
 131   {133683,4177},{65281,4080},{629943,3986},{374213,3898},
 132   {708079,3869},{280125,3785},{641833,3731},{618771,3705},
 133   {930477,3578},{778747,3415},{623751,3371},{40201,3350},
 134   {122389,3307},{950371,3167},{1042353,3111},{18131,3021},
 135   {285429,3004},{549537,2970},{166487,2920},{294287,2857},
 136   {919261,2811},{636339,2766},{900735,2737},{118605,2695},
 137   {10565,2641},{188273,2614},{115369,2563},{735755,2502},
 138   {458285,2490},{914767,2432},{370513,2421},{1027079,2388},
 139   {629619,2366},{462401,2335},{649337,2294},{316165,2274},
 140   {484655,2264},{65115,2245},{326175,2189},{1016279,2153},
 141   {990915,2135},{556859,2101},{462791,2084},{844629,2060},
 142   {404537,2012},{457123,2004},{577589,1938},{638347,1916},
 143   {892325,1882},{182523,1862},{1002505,1842},{624371,1836},
 144   {69057,1817},{210787,1786},{558769,1768},{395623,1750},
 145   {992745,1744},{317855,1727},{384877,1710},{372185,1699},
 146   {105027,1693},{423751,1661},{408961,1635},{908331,1630},
 147   {74551,1620},{36933,1605},{617371,1591},{506045,1586},
 148   {24929,1558},{529709,1548},{1042435,1535},{31867,1517},
 149   {166037,1495},{928781,1478},{508975,1458},{4327,1442},
 150   {779637,1430},{742091,1418},{258263,1411},{879631,1396},
 151   {72029,1385},{728905,1377},{589057,1363},{348621,1356},
 152   {671515,1332},{710453,1315},{84249,1296},{959363,1292},
 153   {685853,1277},{467591,1274},{646643,1267},{683029,1264},
 154   {439927,1249},{254461,1229},{660713,1223},{554195,1220},
 155   {202911,1215},{753253,1195},{941457,1190},{776635,1187},
 156   {509511,1182},{986147,1156},{768879,1151},{699431,1140},
 157   {696417,1128},{86169,1119},{808997,1114},{25467,1107},
 158   {201353,1100},{708087,1084},{1018339,1079},{341297,1073},
 159   {434151,1066},{96287,1058},{950765,1051},{298257,1039},
 160   {675933,1035},{167731,1029},{815445,1027},
 161 };
 162
 163 /* Element j gives the index of the first prime of size 3+j bits */
 164 static uint8_t
 165 prime_by_size[9] = {
 166   1,3,5,10,17,30,53,96,171
 167 };
 168
 169 /* Combined Miller-Rabin test to the base a, and checking the
 170    conditions from Pocklington's theorem. */
 171 static int
 172 miller_rabin_pocklington(mpz_t n, mpz_t nm1, mpz_t nm1dq, mpz_t a)
 173 {
 174   mpz_t r;
 175   mpz_t y;
 176   int is_prime = 0;
 177
 178   /* Avoid the mp_bitcnt_t type for compatibility with older GMP
 179      versions. */
 180   unsigned k;
 181   unsigned j;
 182
 183   VERBOSE(".");
 184
 185   if (mpz_even_p(n) || mpz_cmp_ui(n, 3) < 0)
 186     return 0;
 187
 188   mpz_init(r);
 189   mpz_init(y);
 190
 191   k = mpz_scan1(nm1, 0);
 192   assert(k > 0);
 193
 194   mpz_fdiv_q_2exp (r, nm1, k);
 195
 196   mpz_powm(y, a, r, n);
 197
 198   if (mpz_cmp_ui(y, 1) == 0 || mpz_cmp(y, nm1) == 0)
 199     goto passed_miller_rabin;
 200
 201   for (j = 1; j < k; j++)
 202     {
 203       mpz_powm_ui (y, y, 2, n);
 204
 205       if (mpz_cmp_ui (y, 1) == 0)
 206         break;
 207
 208       if (mpz_cmp (y, nm1) == 0)
 209         {
 210         passed_miller_rabin:
 211           /* We know that a^{n-1} = 1 (mod n)
 212
 213              Remains to check that gcd(a^{(n-1)/q} - 1, n) == 1 */
 214           VERBOSE("x");
 215
 216           mpz_powm(y, a, nm1dq, n);
 217           mpz_sub_ui(y, y, 1);
 218           mpz_gcd(y, y, n);
 219           is_prime = mpz_cmp_ui (y, 1) == 0;
 220           VERBOSE(is_prime ? "\n" : "");
 221           break;
 222         }
 223
 224     }
 225
 226   mpz_clear(r);
 227   mpz_clear(y);
 228
 229   return is_prime;
 230 }
 231
 232 /* The algorithm is based on the following special case of
 233    Pocklington's theorem:
 234
 235    Assume that n = 1 + f q, where q is a prime, q > sqrt(n) - 1. If we
 236    can find an a such that
 237
 238      a^{n-1} = 1 (mod n)
 239      gcd(a^f - 1, n) = 1
 240
 241    then n is prime.
 242
 243    Proof: Assume that n is composite, with smallest prime factor p <=
 244    sqrt(n). Since q is prime, and q > sqrt(n) - 1 >= p - 1, q and p-1
 245    are coprime, so that we can define u = q^{-1} (mod (p-1)). The
 246    assumption a^{n-1} = 1 (mod n) implies that also a^{n-1} = 1 (mod
 247    p). Since p is prime, we have a^{(p-1)} = 1 (mod p). Now, r =
 248    (n-1)/q = (n-1) u (mod (p-1)), and it follows that a^r = a^{(n-1)
 249    u} = 1 (mod p). Then p is a common factor of a^r - 1 and n. This
 250    contradicts gcd(a^r - 1, n) = 1, and concludes the proof.
 251
 252    If n is specified as k bits, we need q of size ceil(k/2) + 1 bits
 253    (or more) to make the theorem apply.
 254 */
 255
 256 /* Generate a prime number p of size bits with 2 p0q dividing (p-1).
 257    p0 must be of size >= ceil(bits/2) + 1. The extra factor q can be
 258    omitted. If top_bits_set is one, the top most two bits are one,
 259    suitable for RSA primes. */
 260 void
 261 _nettle_generate_pocklington_prime (mpz_t p, mpz_t r,
 262                                     unsigned bits, int top_bits_set,
 263                                     void *ctx, nettle_random_func *random,
 264                                     const mpz_t p0,
 265                                     const mpz_t q,
 266                                     const mpz_t p0q)
 267 {
 268   mpz_t r_min, r_range, pm1,a;
 269
 270   assert (2*mpz_sizeinbase (p0, 2) > bits + 1);
 271
 272   mpz_init (r_min);
 273   mpz_init (r_range);
 274   mpz_init (pm1);
 275   mpz_init (a);
 276
 277   if (top_bits_set)
 278     {
 279       /* i = floor (2^{bits-3} / p0q), then 3I + 3 <= r <= 4I, with I
 280          - 2 possible values. */
 281       mpz_set_ui (r_min, 1);
 282       mpz_mul_2exp (r_min, r_min, bits-3);
 283       mpz_fdiv_q (r_min, r_min, p0q);
 284       mpz_sub_ui (r_range, r_min, 2);
 285       mpz_mul_ui (r_min, r_min, 3);
 286       mpz_add_ui (r_min, r_min, 3);
 287     }
 288   else
 289     {
 290       /* i = floor (2^{bits-2} / p0q), I + 1 <= r <= 2I */
 291       mpz_set_ui (r_range, 1);
 292       mpz_mul_2exp (r_range, r_range, bits-2);
 293       mpz_fdiv_q (r_range, r_range, p0q);
 294       mpz_add_ui (r_min, r_range, 1);
 295     }
 296   for (;;)
 297     {
 298       uint8_t buf[1];
 299
 300       nettle_mpz_random (r, ctx, random, r_range);
 301       mpz_add (r, r, r_min);
 302
 303       /* Set p = 2*r*p0q + 1 */
 304       mpz_mul_2exp(r, r, 1);
 305       mpz_mul (pm1, r, p0q);
 306       mpz_add_ui (p, pm1, 1);
 307
 308       assert(mpz_sizeinbase(p, 2) == bits);
 309
 310       /* Should use GMP trial division interface when that
 311          materializes, we don't need any testing beyond trial
 312          division. */
 313       if (!mpz_probab_prime_p (p, 1))
 314         continue;
 315
 316       random(ctx, sizeof(buf), buf);
 317
 318       mpz_set_ui (a, buf[0] + 2);
 319
 320       if (q)
 321         {
 322           mpz_t e;
 323           int is_prime;
 324
 325           mpz_init (e);
 326
 327           mpz_mul (e, r, q);
 328           is_prime = miller_rabin_pocklington(p, pm1, e, a);
 329           mpz_clear (e);
 330
 331           if (is_prime)
 332             break;
 333         }
 334       else if (miller_rabin_pocklington(p, pm1, r, a))
 335         break;
 336     }
 337   mpz_clear (r_min);
 338   mpz_clear (r_range);
 339   mpz_clear (pm1);
 340   mpz_clear (a);
 341 }
 342
 343 /* Generate random prime of a given size. Maurer's algorithm (Alg.
 344    6.42 Handbook of applied cryptography), but with ratio = 1/2 (like
 345    the variant in fips186-3). */
 346 void
 347 nettle_random_prime(mpz_t p, unsigned bits, int top_bits_set,
 348                     void *random_ctx, nettle_random_func *random,
 349                     void *progress_ctx, nettle_progress_func *progress)
 350 {
 351   assert (bits >= 3);
 352   if (bits <= 10)
 353     {
 354       unsigned first;
 355       unsigned choices;
 356       uint8_t buf;
 357
 358       assert (!top_bits_set);
 359
 360       random (random_ctx, sizeof(buf), &buf);
 361
 362       first = prime_by_size[bits-3];
 363       choices = prime_by_size[bits-2] - first;
 364
 365       mpz_set_ui (p, primes[first + buf % choices]);
 366     }
 367   else if (bits <= 20)
 368     {
 369       unsigned long highbit;
 370       uint8_t buf[3];
 371       unsigned long x;
 372       unsigned j;
 373
 374       assert (!top_bits_set);
 375
 376       highbit = 1L << (bits - 1);
 377
 378     again:
 379       random (random_ctx, sizeof(buf), buf);
 380       x = READ_UINT24(buf);
 381       x &= (highbit - 1);
 382       x |= highbit | 1;
 383
 384       for (j = 0; prime_square[j] <= x; j++)
 385         {
 386           unsigned q = x * trial_div_table[j].inverse & TRIAL_DIV_MASK;
 387           if (q <= trial_div_table[j].limit)
 388             goto again;
 389         }
 390       mpz_set_ui (p, x);
 391     }
 392   else
 393     {
 394       mpz_t q, r;
 395
 396       mpz_init (q);
 397       mpz_init (r);
 398
 399      /* Bit size ceil(k/2) + 1, slightly larger than used in Alg. 4.62
 400         in Handbook of Applied Cryptography (which seems to be
 401         incorrect for odd k). */
 402       nettle_random_prime (q, (bits+3)/2, 0, random_ctx, random,
 403                            progress_ctx, progress);
 404
 405       _nettle_generate_pocklington_prime (p, r, bits, top_bits_set,
 406                                           random_ctx, random,
 407                                           q, NULL, q);
 408
 409       if (progress)
 410         progress (progress_ctx, 'x');
 411
 412       mpz_clear (q);
 413       mpz_clear (r);
 414     }
 415 }