monotone

monotone Mtn Source Tree

Root/botan/numthry.cpp

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

Archive Download this file

Branches

Tags

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