Botan  1.10.9
numthry.cpp
Go to the documentation of this file.
1 /*
2 * Number Theory Functions
3 * (C) 1999-2011 Jack Lloyd
4 *
5 * Distributed under the terms of the Botan license
6 */
7 
8 #include <botan/numthry.h>
9 #include <botan/reducer.h>
10 #include <botan/internal/bit_ops.h>
11 #include <algorithm>
12 
13 namespace Botan {
14 
15 namespace {
16 
17 /*
18 * Miller-Rabin Primality Tester
19 */
20 class MillerRabin_Test
21  {
22  public:
23  bool is_witness(const BigInt& nonce);
24  MillerRabin_Test(const BigInt& num);
25  private:
26  BigInt n, r, n_minus_1;
27  size_t s;
28  Fixed_Exponent_Power_Mod pow_mod;
29  Modular_Reducer reducer;
30  };
31 
32 /*
33 * Miller-Rabin Test, as described in Handbook of Applied Cryptography
34 * section 4.24
35 */
36 bool MillerRabin_Test::is_witness(const BigInt& a)
37  {
38  if(a < 2 || a >= n_minus_1)
39  throw Invalid_Argument("Bad size for nonce in Miller-Rabin test");
40 
41  BigInt y = pow_mod(a);
42  if(y == 1 || y == n_minus_1)
43  return false;
44 
45  for(size_t i = 1; i != s; ++i)
46  {
47  y = reducer.square(y);
48 
49  if(y == 1) // found a non-trivial square root
50  return true;
51 
52  if(y == n_minus_1) // -1, trivial square root, so give up
53  return false;
54  }
55 
56  // If we reached here then n fails the Fermat test
57  return true;
58  }
59 
60 /*
61 * Miller-Rabin Constructor
62 */
63 MillerRabin_Test::MillerRabin_Test(const BigInt& num)
64  {
65  if(num.is_even() || num < 3)
66  throw Invalid_Argument("MillerRabin_Test: Invalid number for testing");
67 
68  n = num;
69  n_minus_1 = n - 1;
71  r = n_minus_1 >> s;
72 
73  pow_mod = Fixed_Exponent_Power_Mod(r, n);
74  reducer = Modular_Reducer(n);
75  }
76 
77 /*
78 * Miller-Rabin Iterations
79 */
80 size_t miller_rabin_test_iterations(size_t bits, size_t level)
81  {
82  struct mapping { size_t bits; size_t verify_iter; size_t check_iter; };
83 
84  static const mapping tests[] = {
85  { 50, 55, 25 },
86  { 100, 38, 22 },
87  { 160, 32, 18 },
88  { 163, 31, 17 },
89  { 168, 30, 16 },
90  { 177, 29, 16 },
91  { 181, 28, 15 },
92  { 185, 27, 15 },
93  { 190, 26, 15 },
94  { 195, 25, 14 },
95  { 201, 24, 14 },
96  { 208, 23, 14 },
97  { 215, 22, 13 },
98  { 222, 21, 13 },
99  { 231, 20, 13 },
100  { 241, 19, 12 },
101  { 252, 18, 12 },
102  { 264, 17, 12 },
103  { 278, 16, 11 },
104  { 294, 15, 10 },
105  { 313, 14, 9 },
106  { 334, 13, 8 },
107  { 360, 12, 8 },
108  { 392, 11, 7 },
109  { 430, 10, 7 },
110  { 479, 9, 6 },
111  { 542, 8, 6 },
112  { 626, 7, 5 },
113  { 746, 6, 4 },
114  { 926, 5, 3 },
115  { 1232, 4, 2 },
116  { 1853, 3, 2 },
117  { 0, 0, 0 }
118  };
119 
120  for(size_t i = 0; tests[i].bits; ++i)
121  {
122  if(bits <= tests[i].bits)
123  {
124  if(level >= 2)
125  return tests[i].verify_iter;
126  else if(level == 1)
127  return tests[i].check_iter;
128  else if(level == 0)
129  return std::max<size_t>(tests[i].check_iter / 4, 1);
130  }
131  }
132 
133  return level > 0 ? 2 : 1; // for large inputs
134  }
135 
136 }
137 
138 /*
139 * Return the number of 0 bits at the end of n
140 */
141 size_t low_zero_bits(const BigInt& n)
142  {
143  size_t low_zero = 0;
144 
145  if(n.is_positive() && n.is_nonzero())
146  {
147  for(size_t i = 0; i != n.size(); ++i)
148  {
149  word x = n[i];
150 
151  if(x)
152  {
153  low_zero += ctz(x);
154  break;
155  }
156  else
157  low_zero += BOTAN_MP_WORD_BITS;
158  }
159  }
160 
161  return low_zero;
162  }
163 
164 /*
165 * Calculate the GCD
166 */
167 BigInt gcd(const BigInt& a, const BigInt& b)
168  {
169  if(a.is_zero() || b.is_zero()) return 0;
170  if(a == 1 || b == 1) return 1;
171 
172  BigInt x = a, y = b;
174  y.set_sign(BigInt::Positive);
175  size_t shift = std::min(low_zero_bits(x), low_zero_bits(y));
176 
177  x >>= shift;
178  y >>= shift;
179 
180  while(x.is_nonzero())
181  {
182  x >>= low_zero_bits(x);
183  y >>= low_zero_bits(y);
184  if(x >= y) { x -= y; x >>= 1; }
185  else { y -= x; y >>= 1; }
186  }
187 
188  return (y << shift);
189  }
190 
191 /*
192 * Calculate the LCM
193 */
194 BigInt lcm(const BigInt& a, const BigInt& b)
195  {
196  return ((a * b) / gcd(a, b));
197  }
198 
199 /*
200 * Find the Modular Inverse
201 */
203  {
204  if(mod.is_zero())
205  throw BigInt::DivideByZero();
206  if(mod.is_negative() || n.is_negative())
207  throw Invalid_Argument("inverse_mod: arguments must be non-negative");
208 
209  if(n.is_zero() || (n.is_even() && mod.is_even()))
210  return 0;
211 
212  BigInt x = mod, y = n, u = mod, v = n;
213  BigInt A = 1, B = 0, C = 0, D = 1;
214 
215  while(u.is_nonzero())
216  {
217  size_t zero_bits = low_zero_bits(u);
218  u >>= zero_bits;
219  for(size_t i = 0; i != zero_bits; ++i)
220  {
221  if(A.is_odd() || B.is_odd())
222  { A += y; B -= x; }
223  A >>= 1; B >>= 1;
224  }
225 
226  zero_bits = low_zero_bits(v);
227  v >>= zero_bits;
228  for(size_t i = 0; i != zero_bits; ++i)
229  {
230  if(C.is_odd() || D.is_odd())
231  { C += y; D -= x; }
232  C >>= 1; D >>= 1;
233  }
234 
235  if(u >= v) { u -= v; A -= C; B -= D; }
236  else { v -= u; C -= A; D -= B; }
237  }
238 
239  if(v != 1)
240  return 0;
241 
242  while(D.is_negative()) D += mod;
243  while(D >= mod) D -= mod;
244 
245  return D;
246  }
247 
248 /*
249 * Modular Exponentiation
250 */
251 BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
252  {
253  Power_Mod pow_mod(mod);
254  pow_mod.set_base(base);
255  pow_mod.set_exponent(exp);
256  return pow_mod.execute();
257  }
258 
259 /*
260 * Test for primaility using Miller-Rabin
261 */
262 bool primality_test(const BigInt& n,
264  size_t level)
265  {
266  const size_t PREF_NONCE_BITS = 128;
267 
268  if(n == 2)
269  return true;
270  if(n <= 1 || n.is_even())
271  return false;
272 
273  // Fast path testing for small numbers (<= 65521)
274  if(n <= PRIMES[PRIME_TABLE_SIZE-1])
275  {
276  const word num = n.word_at(0);
277 
278  for(size_t i = 0; PRIMES[i]; ++i)
279  {
280  if(num == PRIMES[i])
281  return true;
282  if(num < PRIMES[i])
283  return false;
284  }
285 
286  return false;
287  }
288 
289  if(level > 2)
290  level = 2;
291 
292  const size_t NONCE_BITS = std::min(n.bits() - 2, PREF_NONCE_BITS);
293 
294  MillerRabin_Test mr(n);
295 
296  if(mr.is_witness(2))
297  return false;
298 
299  const size_t tests = miller_rabin_test_iterations(n.bits(), level);
300 
301  for(size_t i = 0; i != tests; ++i)
302  {
303  BigInt nonce;
304  while(nonce < 2 || nonce >= (n-1))
305  nonce.randomize(rng, NONCE_BITS);
306 
307  if(mr.is_witness(nonce))
308  return false;
309  }
310 
311  return true;
312  }
313 
314 }
const size_t PRIME_TABLE_SIZE
Definition: numthry.h:220
word word_at(size_t n) const
Definition: bigint.h:238
size_t ctz(T n)
Definition: bit_ops.h:93
bool is_odd() const
Definition: bigint.h:164
bool is_even() const
Definition: bigint.h:158
BigInt n
Definition: numthry.cpp:26
BigInt gcd(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:167
bool is_negative() const
Definition: bigint.h:245
std::invalid_argument Invalid_Argument
Definition: exceptn.h:20
size_t size() const
Definition: bigint.h:284
Modular_Reducer reducer
Definition: numthry.cpp:29
BigInt n_minus_1
Definition: numthry.cpp:26
bool is_nonzero() const
Definition: bigint.h:170
Fixed_Exponent_Power_Mod pow_mod
Definition: numthry.cpp:28
size_t bits() const
Definition: bigint.cpp:253
bool primality_test(const BigInt &n, RandomNumberGenerator &rng, size_t level)
Definition: numthry.cpp:262
RandomNumberGenerator * rng
Definition: global_rng.cpp:165
void randomize(RandomNumberGenerator &rng, size_t bitsize=0)
Definition: big_rand.cpp:29
GMP_MPZ exp
Definition: gmp_powm.cpp:29
GMP_MPZ mod
Definition: gmp_powm.cpp:29
void set_base(const BigInt &) const
Definition: pow_mod.cpp:83
size_t low_zero_bits(const BigInt &n)
Definition: numthry.cpp:141
BigInt inverse_mod(const BigInt &n, const BigInt &mod)
Definition: numthry.cpp:202
BigInt power_mod(const BigInt &base, const BigInt &exp, const BigInt &mod)
Definition: numthry.cpp:251
bool is_positive() const
Definition: bigint.h:251
BigInt r
Definition: numthry.cpp:26
const u16bit BOTAN_DLL PRIMES[]
Definition: primes.cpp:12
bool is_zero() const
Definition: bigint.h:176
BigInt square(const BigInt &x) const
Definition: reducer.h:39
void set_sign(Sign sign)
Definition: bigint.cpp:291
BigInt lcm(const BigInt &a, const BigInt &b)
Definition: numthry.cpp:194
BigInt execute() const
Definition: pow_mod.cpp:109
size_t s
Definition: numthry.cpp:27
void set_exponent(const BigInt &) const
Definition: pow_mod.cpp:96
GMP_MPZ base
Definition: gmp_powm.cpp:29