monotone

monotone Mtn Source Tree

Root/botan/numthry.cpp

1/*************************************************
2* Number Theory Source File *
3* (C) 1999-2006 The Botan Project *
4*************************************************/
5
6#include <botan/numthry.h>
7#include <botan/libstate.h>
8#include <algorithm>
9
10namespace Botan {
11
12namespace {
13
14/*************************************************
15* Miller-Rabin Iterations *
16*************************************************/
17u32bit miller_rabin_test_iterations(u32bit bits, bool verify)
18 {
19 struct mapping { u32bit bits; u32bit verify_iter; u32bit check_iter; };
20
21 static const mapping tests[] = {
22 { 50, 55, 25 },
23 { 100, 38, 22 },
24 { 160, 32, 18 },
25 { 163, 31, 17 },
26 { 168, 30, 16 },
27 { 177, 29, 16 },
28 { 181, 28, 15 },
29 { 185, 27, 15 },
30 { 190, 26, 15 },
31 { 195, 25, 14 },
32 { 201, 24, 14 },
33 { 208, 23, 14 },
34 { 215, 22, 13 },
35 { 222, 21, 13 },
36 { 231, 20, 13 },
37 { 241, 19, 12 },
38 { 252, 18, 12 },
39 { 264, 17, 12 },
40 { 278, 16, 11 },
41 { 294, 15, 10 },
42 { 313, 14, 9 },
43 { 334, 13, 8 },
44 { 360, 12, 8 },
45 { 392, 11, 7 },
46 { 430, 10, 7 },
47 { 479, 9, 6 },
48 { 542, 8, 6 },
49 { 626, 7, 5 },
50 { 746, 6, 4 },
51 { 926, 5, 3 },
52 { 1232, 4, 2 },
53 { 1853, 3, 2 },
54 { 0, 0, 0 }
55 };
56
57 for(u32bit j = 0; tests[j].bits; ++j)
58 {
59 if(bits <= tests[j].bits)
60 if(verify)
61 return tests[j].verify_iter;
62 else
63 return tests[j].check_iter;
64 }
65 return 2;
66 }
67
68}
69
70/*************************************************
71* Return the number of 0 bits at the end of n *
72*************************************************/
73u32bit low_zero_bits(const BigInt& n)
74 {
75 if(n.is_zero()) return 0;
76
77 u32bit bits = 0, max_bits = n.bits();
78 while((n.get_bit(bits) == 0) && bits < max_bits)
79 ++bits;
80 return bits;
81 }
82
83/*************************************************
84* Calculate the GCD *
85*************************************************/
86BigInt gcd(const BigInt& a, const BigInt& b)
87 {
88 if(a.is_zero() || b.is_zero()) return 0;
89 if(a == 1 || b == 1) return 1;
90
91 BigInt x = a, y = b;
92 x.set_sign(BigInt::Positive);
93 y.set_sign(BigInt::Positive);
94 u32bit shift = std::min(low_zero_bits(x), low_zero_bits(y));
95
96 x >>= shift;
97 y >>= shift;
98
99 while(x.is_nonzero())
100 {
101 x >>= low_zero_bits(x);
102 y >>= low_zero_bits(y);
103 if(x >= y) { x -= y; x >>= 1; }
104 else { y -= x; y >>= 1; }
105 }
106
107 return (y << shift);
108 }
109
110/*************************************************
111* Calculate the LCM *
112*************************************************/
113BigInt lcm(const BigInt& a, const BigInt& b)
114 {
115 return ((a * b) / gcd(a, b));
116 }
117
118/*************************************************
119* Find the Modular Inverse *
120*************************************************/
121BigInt inverse_mod(const BigInt& n, const BigInt& mod)
122 {
123 if(mod.is_zero())
124 throw BigInt::DivideByZero();
125 if(mod.is_negative() || n.is_negative())
126 throw Invalid_Argument("inverse_mod: arguments must be non-negative");
127
128 if(n.is_zero() || (n.is_even() && mod.is_even()))
129 return 0;
130
131 BigInt x = mod, y = n, u = mod, v = n;
132 BigInt A = 1, B = 0, C = 0, D = 1;
133
134 while(u.is_nonzero())
135 {
136 u32bit zero_bits = low_zero_bits(u);
137 u >>= zero_bits;
138 for(u32bit j = 0; j != zero_bits; ++j)
139 {
140 if(A.is_odd() || B.is_odd())
141 { A += y; B -= x; }
142 A >>= 1; B >>= 1;
143 }
144
145 zero_bits = low_zero_bits(v);
146 v >>= zero_bits;
147 for(u32bit j = 0; j != zero_bits; ++j)
148 {
149 if(C.is_odd() || D.is_odd())
150 { C += y; D -= x; }
151 C >>= 1; D >>= 1;
152 }
153
154 if(u >= v) { u -= v; A -= C; B -= D; }
155 else { v -= u; C -= A; D -= B; }
156 }
157
158 if(v != 1)
159 return 0;
160
161 while(D.is_negative()) D += mod;
162 while(D >= mod) D -= mod;
163
164 return D;
165 }
166
167/*************************************************
168* Modular Exponentiation *
169*************************************************/
170BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod)
171 {
172 Power_Mod pow_mod(mod);
173 pow_mod.set_base(base);
174 pow_mod.set_exponent(exp);
175 return pow_mod.execute();
176 }
177
178/*************************************************
179* Do simple tests of primality *
180*************************************************/
181s32bit simple_primality_tests(const BigInt& n)
182 {
183 const s32bit NOT_PRIME = -1, UNKNOWN = 0, PRIME = 1;
184
185 if(n == 2)
186 return PRIME;
187 if(n <= 1 || n.is_even())
188 return NOT_PRIME;
189
190 if(n <= PRIMES[PRIME_TABLE_SIZE-1])
191 {
192 const word num = n.word_at(0);
193 for(u32bit j = 0; PRIMES[j]; ++j)
194 {
195 if(num == PRIMES[j]) return PRIME;
196 if(num < PRIMES[j]) return NOT_PRIME;
197 }
198 return NOT_PRIME;
199 }
200
201 u32bit check_first = std::min(n.bits() / 32, PRIME_PRODUCTS_TABLE_SIZE);
202 for(u32bit j = 0; j != check_first; ++j)
203 if(gcd(n, PRIME_PRODUCTS[j]) != 1)
204 return NOT_PRIME;
205
206 return UNKNOWN;
207 }
208
209/*************************************************
210* Fast check of primality *
211*************************************************/
212bool check_prime(const BigInt& n)
213 {
214 return run_primality_tests(n, 0);
215 }
216
217/*************************************************
218* Test for primality *
219*************************************************/
220bool is_prime(const BigInt& n)
221 {
222 return run_primality_tests(n, 1);
223 }
224
225/*************************************************
226* Verify primality *
227*************************************************/
228bool verify_prime(const BigInt& n)
229 {
230 return run_primality_tests(n, 2);
231 }
232
233/*************************************************
234* Verify primality *
235*************************************************/
236bool run_primality_tests(const BigInt& n, u32bit level)
237 {
238 s32bit simple_tests = simple_primality_tests(n);
239 if(simple_tests) return (simple_tests == 1) ? true : false;
240 return passes_mr_tests(n, level);
241 }
242
243/*************************************************
244* Test for primaility using Miller-Rabin *
245*************************************************/
246bool passes_mr_tests(const BigInt& n, u32bit level)
247 {
248 const u32bit PREF_NONCE_BITS = 40;
249
250 if(level > 2)
251 level = 2;
252
253 MillerRabin_Test mr(n);
254
255 if(!mr.passes_test(2))
256 return false;
257
258 if(level == 0)
259 return true;
260
261 const u32bit NONCE_BITS = std::min(n.bits() - 1, PREF_NONCE_BITS);
262
263 const bool verify = (level == 2);
264
265 u32bit tests = miller_rabin_test_iterations(n.bits(), verify);
266
267 BigInt nonce;
268 for(u32bit j = 0; j != tests; ++j)
269 {
270 if(verify) nonce = random_integer(NONCE_BITS);
271 else nonce = PRIMES[j];
272
273 if(!mr.passes_test(nonce))
274 return false;
275 }
276 return true;
277 }
278
279/*************************************************
280* Miller-Rabin Test *
281*************************************************/
282bool MillerRabin_Test::passes_test(const BigInt& a)
283 {
284 if(a < 2 || a >= n_minus_1)
285 throw Invalid_Argument("Bad size for nonce in Miller-Rabin test");
286
287 global_state().pulse(PRIME_TESTING);
288
289 BigInt y = pow_mod(a);
290 if(y == 1 || y == n_minus_1)
291 return true;
292
293 for(u32bit j = 1; j != s; ++j)
294 {
295 global_state().pulse(PRIME_TESTING);
296 y = reducer.square(y);
297
298 if(y == 1)
299 return false;
300 if(y == n_minus_1)
301 return true;
302 }
303 return false;
304 }
305
306/*************************************************
307* Miller-Rabin Constructor *
308*************************************************/
309MillerRabin_Test::MillerRabin_Test(const BigInt& num)
310 {
311 if(num.is_even() || num < 3)
312 throw Invalid_Argument("MillerRabin_Test: Invalid number for testing");
313
314 n = num;
315 n_minus_1 = n - 1;
316 s = low_zero_bits(n_minus_1);
317 r = n_minus_1 >> s;
318
319 pow_mod = Fixed_Exponent_Power_Mod(r, n);
320 reducer = Modular_Reducer(n);
321 }
322
323}

Archive Download this file

Branches

Tags

Quick Links:     www.monotone.ca    -     Downloads    -     Documentation    -     Wiki    -     Code Forge    -     Build Status