1 /**
2 * Number Theory Functions
3 *
4 * Copyright:
5 * (C) 1999-2007 Jack Lloyd
6 * (C) 2014-2015 Etienne Cimon
7 *
8 * License:
9 * Botan is released under the Simplified BSD License (see LICENSE.md)
10 */
11 module botan.math.numbertheory.numthry;
12
13 import botan.constants;
14 static if (BOTAN_HAS_PUBLIC_KEY_CRYPTO):
15
16 public import botan.math.bigint.bigint;
17 public import botan.math.numbertheory.pow_mod;
18 public import botan.math.numbertheory.primes;
19 public import botan.utils.types;
20 import botan.rng.rng;
21 import botan.hash.hash;
22 import botan.utils.parsing;
23 import std.algorithm;
24 import botan.math.numbertheory.reducer;
25 import botan.utils.bit_ops;
26 import botan.math.mp.mp_core;
27 import botan.utils.rounding;
28 import botan.algo_factory.algo_factory : AlgorithmFactory;
29 import botan.constants;
30 import std.conv : to;
31
32 /**
33 * Fused multiply-add
34 * Params:
35 * a = an integer
36 * b = an integer
37 * c = an integer
38 * Returns: (a*b)+c
39 */
40 /*
41 * Multiply-Add Operation
42 */
43 BigInt mulAdd()(auto const ref BigInt a, auto const ref BigInt b, auto const ref BigInt c)
44 {
45 if (c.isNegative() || c.isZero())
46 throw new InvalidArgument("mulAdd: Third argument must be > 0");
47
48 BigInt.Sign sign = BigInt.Positive;
49 if (a.sign() != b.sign())
50 sign = BigInt.Negative;
51
52 const size_t a_sw = a.sigWords();
53 const size_t b_sw = b.sigWords();
54 const size_t c_sw = c.sigWords();
55
56 BigInt r = BigInt(sign, std.algorithm.max(a.length + b.length, c_sw) + 1);
57 SecureVector!word workspace = SecureVector!word(r.length);
58
59 bigint_mul(r.mutablePtr(), r.length, workspace.ptr, a.ptr, a.length, a_sw, b.ptr, b.length, b_sw);
60
61 const size_t r_size = std.algorithm.max(r.sigWords(), c_sw);
62 bigint_add2(r.mutablePtr(), r_size, c.ptr, c_sw);
63 return r;
64 }
65
66
67 /**
68 * Fused subtract-multiply
69 * Params:
70 * a = an integer
71 * b = an integer
72 * c = an integer
73 * Returns: (a-b)*c
74 */
75 BigInt subMul()(auto const ref BigInt a, auto const ref BigInt b, auto const ref BigInt c)
76 {
77 if (a.isNegative() || b.isNegative())
78 throw new InvalidArgument("subMul: First two arguments must be >= 0");
79
80 BigInt r = a.dup;
81 r -= b;
82 r *= c;
83 return r.move();
84 }
85
86 /**
87 * Return the absolute value
88 * Params:
89 * n = an integer
90 * Returns: absolute value of n
91 */
92 BigInt abs()(auto const ref BigInt n) { return n.abs(); }
93
94 /**
95 * Compute the greatest common divisor
96 * Params:
97 * a = positive integer x
98 * b = positive integer y
99 * Returns: gcd(x,y)
100 */
101 BigInt gcd()(auto const ref BigInt a, auto const ref BigInt b)
102 {
103 if (a.isZero() || b.isZero()) return BigInt(0);
104 if (a == 1 || b == 1) return BigInt(1);
105
106 BigInt x = a.dup, y = b.dup;
107 x.setSign(BigInt.Positive);
108 y.setSign(BigInt.Positive);
109 size_t shift = std.algorithm.min(lowZeroBits(x), lowZeroBits(y));
110
111 x >>= shift;
112 y >>= shift;
113
114 while (x.isNonzero())
115 {
116 x >>= lowZeroBits(x);
117 y >>= lowZeroBits(y);
118 if (x >= y) { x -= y; x >>= 1; }
119 else { y -= x; y >>= 1; }
120 }
121
122 return (y << shift);
123 }
124
125 /**
126 * Least common multiple
127 * Params:
128 * a = a positive integer x
129 * b = a positive integer y
130 * Returns: z, smallest integer such that z % x == 0 and z % y == 0
131 */
132 BigInt lcm()(auto const ref BigInt a, auto const ref BigInt b)
133 {
134 auto a_b = (a * b);
135 auto gcd_a_b = gcd(a, b);
136 return a_b/gcd_a_b;
137 }
138
139
140 /**
141 * Params:
142 * x = an integer
143 * Returns: (x*x)
144 */
145 BigInt square()(auto const ref BigInt x)
146 {
147 const size_t x_sw = x.sigWords();
148
149 BigInt z = BigInt(BigInt.Positive, roundUp!size_t(2*x_sw, 16));
150 SecureVector!word workspace = SecureVector!word(z.length);
151
152 bigint_sqr(z.mutablePtr(), z.length, workspace.ptr, x.ptr, x.length, x_sw);
153 return z;
154 }
155
156 /**
157 * Modular inversion
158 * Params:
159 * n = a positive integer x
160 * mod = a positive integer modulus
161 * Returns: y st (x*y) % modulus == 1
162 */
163 BigInt inverseMod()(auto const ref BigInt n, auto const ref BigInt mod)
164 {
165 if (mod.isZero())
166 throw new BigInt.DivideByZero();
167 if (mod.isNegative() || n.isNegative())
168 throw new InvalidArgument("inverseMod: arguments must be non-negative");
169
170 if (n.isZero() || (n.isEven() && mod.isEven()))
171 return BigInt(0); // fast fail checks
172
173 if (mod.isOdd())
174 return inverseModOddModulus(n, mod);
175
176 BigInt u = mod.dup, v = n.dup;
177 BigInt A = BigInt(1);
178 BigInt B = BigInt(0);
179 BigInt C = BigInt(0);
180 BigInt D = BigInt(1);
181
182 while (u.isNonzero())
183 {
184 const size_t u_zero_bits = lowZeroBits(u);
185 u >>= u_zero_bits;
186 foreach (size_t i; 0 .. u_zero_bits)
187 {
188 if (A.isOdd() || B.isOdd())
189 { A += n; B -= mod; }
190 A >>= 1; B >>= 1;
191 }
192
193 const size_t v_zero_bits = lowZeroBits(v);
194 v >>= v_zero_bits;
195 foreach (size_t i; 0 .. v_zero_bits)
196 {
197 if (C.isOdd() || D.isOdd())
198 { C += n; D -= mod; }
199 C >>= 1; D >>= 1;
200 }
201
202 if (u >= v) { u -= v; A -= C; B -= D; }
203 else { v -= u; C -= A; D -= B; }
204 }
205
206 if (v != 1)
207 return BigInt(0); // no modular inverse
208
209 while (D.isNegative()) D += mod;
210 while (D >= mod) D -= mod;
211
212 return D.move();
213 }
214
215
216 /**
217 * Compute the Jacobi symbol. If n is prime, this is equivalent
218 * to the Legendre symbol.
219 * @see http://mathworld.wolfram.com/JacobiSymbol.html
220 *
221 * Params:
222 * a = is a non-negative integer
223 * n = is an odd integer > 1
224 * Returns: (n / m)
225 */
226 int jacobi()(auto const ref BigInt a, auto const ref BigInt n)
227 {
228 if (a.isNegative())
229 throw new InvalidArgument("jacobi: first argument must be non-negative");
230 if (n.isEven() || n < 2)
231 throw new InvalidArgument("jacobi: second argument must be odd and > 1");
232
233 BigInt x = a.dup, y = n.dup;
234 int J = 1;
235
236 while (y > 1)
237 {
238 x %= y;
239 if (x > y / 2)
240 {
241 x = y - x;
242 if (y % 4 == 3)
243 J = -J;
244 }
245 if (x.isZero())
246 return 0;
247
248 size_t shifts = lowZeroBits(x);
249 x >>= shifts;
250 if (shifts % 2)
251 {
252 word y_mod_8 = y % 8;
253 if (y_mod_8 == 3 || y_mod_8 == 5)
254 J = -J;
255 }
256
257 if (x % 4 == 3 && y % 4 == 3)
258 J = -J;
259 std.algorithm.swap(x, y);
260 }
261 return J;
262 }
263
264 /**
265 * Modular exponentation
266 * Params:
267 * base = an integer base b
268 * exp = a positive exponent x
269 * mod = a positive modulus m
270 * Returns: (b^x) % m
271 */
272 BigInt powerMod()(auto const ref BigInt base, auto const ref BigInt exp, auto const ref BigInt mod)
273 {
274 auto pow_mod = scoped!PowerMod(mod);
275
276 /*
277 * Calling setBase before setExponent means we end up using a
278 * minimal window. This makes sense given that here we know that any
279 * precomputation is wasted.
280 */
281 pow_mod.setBase(base);
282 pow_mod.setExponent(exp);
283 return pow_mod.execute();
284 }
285
286
287 /**
288 * Compute the square root of x modulo a prime using the
289 * Shanks-Tonnelli algorithm
290 *
291 * Params:
292 * a = the input x
293 * p = the prime p
294 * Returns: y such that (y*y)%p == x, or -1 if no such integer
295 */
296
297 /*
298 * Shanks-Tonnelli algorithm
299 */
300 BigInt ressol()(auto const ref BigInt a, auto const ref BigInt p)
301 {
302
303 if (a == 0)
304 return BigInt(0);
305 else if (a < 0)
306 throw new InvalidArgument("ressol(): a to solve for must be positive");
307
308 if (p == 2)
309 return a.dup;
310 else if (p <= 1)
311 throw new InvalidArgument("ressol(): prime must be > 1");
312 else if(p.isEven())
313 throw new InvalidArgument("ressol(): invalid prime");
314 if (jacobi(a, p) != 1) { // not a quadratic residue
315 auto bi = -BigInt(1);
316 return bi.move();
317 }
318
319 if (p % 4 == 3)
320 return powerMod(a, ((p+1) >> 2), p);
321
322 size_t s = lowZeroBits(p - 1);
323 BigInt q = p >> s;
324
325 q -= 1;
326 q >>= 1;
327
328 ModularReducer mod_p = ModularReducer(p);
329
330 BigInt r = powerMod(a, q, p);
331 BigInt n = mod_p.multiply(a, mod_p.square(r));
332 r = mod_p.multiply(r, a);
333
334 if (n == 1)
335 return r.move();
336
337 // find random non quadratic residue z
338 BigInt z = 2;
339 while (jacobi(z, p) == 1) // while z quadratic residue
340 ++z;
341
342 BigInt c = powerMod(z, (q << 1) + 1, p);
343
344 while (n > 1)
345 {
346 q = n.dup();
347
348 size_t i = 0;
349 while (q != 1)
350 {
351 q = mod_p.square(q);
352 ++i;
353 if (i >= s) {
354 auto bi = -BigInt(1);
355 return bi.move();
356 }
357 }
358
359 c = powerMod(c, BigInt.powerOf2(s-i-1), p);
360 r = mod_p.multiply(r, c);
361 c = mod_p.square(c);
362 n = mod_p.multiply(n, c);
363 s = i;
364 }
365
366 return r.move();
367 }
368
369 /*
370 * Compute -input^-1 mod 2^MP_WORD_BITS. Returns zero if input
371 * is even. If input is odd, input and 2^n are relatively prime
372 * and an inverse exists.
373 */
374 word montyInverse(word input)
375 {
376 word b = input;
377 word x2 = 1, x1 = 0, y2 = 0, y1 = 1;
378
379 // First iteration, a = n+1
380 word q = bigint_divop(1, 0, b);
381 word r = (MP_WORD_MAX - q*b) + 1;
382 word x = x2 - q*x1;
383 word y = y2 - q*y1;
384
385 word a = b;
386 b = r;
387 x2 = x1;
388 x1 = x;
389 y2 = y1;
390 y1 = y;
391
392 while (b > 0)
393 {
394 q = a / b;
395 r = a - q*b;
396 x = x2 - q*x1;
397 y = y2 - q*y1;
398
399 a = b;
400 b = r;
401 x2 = x1;
402 x1 = x;
403 y2 = y1;
404 y1 = y;
405 }
406
407 // Now invert in addition space
408 y2 = (MP_WORD_MAX - y2) + 1;
409
410 return y2;
411 }
412
413 /**
414 * Params:
415 * n = a positive integer x
416 * Returns: count of the zero bits in x, or, equivalently, the largest
417 * value of n such that 2^n divides x evenly. Returns zero if
418 * n is less than or equal to zero.
419 */
420 size_t lowZeroBits()(auto const ref BigInt n)
421 {
422 size_t low_zero = 0;
423
424 if (n.isPositive() && n.isNonzero())
425 {
426 for (size_t i = 0; i != n.length; ++i)
427 {
428 const word x = n.wordAt(i);
429
430 if (x)
431 {
432 low_zero += ctz(x);
433 break;
434 }
435 else
436 low_zero += BOTAN_MP_WORD_BITS;
437 }
438 }
439
440 return low_zero;
441 }
442
443 /**
444 * Check for primality using Miller-Rabin
445 * Params:
446 * n = a positive integer to test for primality
447 * rng = a random number generator
448 * prob = chance of false positive is bounded by 1/2**prob
449 * is_random = true if n was randomly chosen by us
450 * Returns: true if all primality tests passed, otherwise false
451 */
452 bool isPrime()(auto const ref BigInt n, RandomNumberGenerator rng, size_t prob = 56, bool is_random = false)
453 {
454 import std.range : assumeSorted, SortedRange, empty;
455 if (n == 2)
456 return true;
457 if (n <= 1 || n.isEven())
458 return false;
459
460 // Fast path testing for small numbers (<= 65521)
461 if (n <= PRIMES[PRIME_TABLE_SIZE-1])
462 {
463 const ushort num = cast(ushort) n.wordAt(0);
464 auto r = assumeSorted(PRIMES[0..$]);
465 return !r.equalRange(num).empty;
466 }
467
468 const size_t test_iterations = mrTestIterations(n.bits(), prob, is_random);
469 const BigInt n_minus_1 = n - 1;
470 const size_t s = lowZeroBits(n_minus_1);
471 FixedExponentPowerModImpl pow_mod = ThreadMem.alloc!FixedExponentPowerModImpl(n_minus_1 >> s, n);
472 scope(exit) ThreadMem.free(pow_mod);
473 ModularReducer reducer = ModularReducer(n);
474
475 foreach (size_t i; 0 .. test_iterations)
476 {
477 auto bi = BigInt(2);
478 const BigInt a = BigInt.randomInteger(rng, bi, n_minus_1);
479 BigInt y = pow_mod(a);
480 if (mrWitness(y, reducer, n_minus_1, s))
481 return false;
482 }
483
484 return true;
485 }
486
487 bool quickCheckPrime(const ref BigInt n, RandomNumberGenerator rng)
488 { return isPrime(n, rng, 32); }
489
490 bool checkPrime(const ref BigInt n, RandomNumberGenerator rng)
491 { return isPrime(n, rng, 56); }
492
493 bool verifyPrime(const ref BigInt n, RandomNumberGenerator rng)
494 { return isPrime(n, rng, 80); }
495
496 /**
497 * Randomly generate a prime
498 * Params:
499 * rng = a random number generator
500 * bits = how large the resulting prime should be in bits
501 * coprime = a positive integer the result should be coprime to
502 * equiv = a non-negative number that the result should be
503 equivalent to modulo equiv_mod
504 * modulo = the modulus equiv should be checked against
505 * Returns: random prime with the specified criteria
506 */
507 BigInt randomPrime()(RandomNumberGenerator rng,
508 size_t bits, const ref BigInt coprime,
509 size_t equiv = 1, size_t modulo = 2)
510 {
511 if (bits <= 1)
512 throw new InvalidArgument("randomPrime: Can't make a prime of " ~ to!string(bits) ~ " bits");
513 else if (bits == 2)
514 return ((rng.nextByte() % 2) ? BigInt(2) : BigInt(3));
515 else if (bits == 3)
516 return ((rng.nextByte() % 2) ? BigInt(5) : BigInt(7));
517 else if (bits == 4)
518 return ((rng.nextByte() % 2) ? BigInt(11) : BigInt(13));
519
520 if (coprime <= 0)
521 throw new InvalidArgument("randomPrime: coprime must be > 0");
522 if (modulo % 2 == 1 || modulo == 0)
523 throw new InvalidArgument("randomPrime: Invalid modulo value");
524 if (equiv >= modulo || equiv % 2 == 0)
525 throw new InvalidArgument("randomPrime: equiv must be < modulo, and odd");
526
527 while (true)
528 {
529 BigInt p = BigInt(rng, bits);
530
531 // Force lowest and two top bits on
532 p.setBit(bits - 1);
533 p.setBit(bits - 2);
534 p.setBit(0);
535
536 if (p % modulo != equiv)
537 p += (modulo - p % modulo) + equiv;
538
539 const size_t sieve_size = std.algorithm.min(bits / 2, PRIME_TABLE_SIZE);
540 SecureVector!ushort sieve = SecureVector!ushort(sieve_size);
541
542 for (size_t j = 0; j != sieve.length; ++j)
543 sieve[j] = cast(ushort)( p % PRIMES[j]);
544
545 size_t counter = 0;
546 while (true)
547 {
548 if (counter == 4096 || p.bits() > bits)
549 break;
550
551 bool passes_sieve = true;
552 ++counter;
553 p += modulo;
554
555 if (p.bits() > bits)
556 break;
557
558 for (size_t j = 0; j != sieve.length; ++j)
559 {
560 sieve[j] = cast(ushort)((sieve[j] + modulo) % PRIMES[j]);
561 if (sieve[j] == 0)
562 passes_sieve = false;
563 }
564 auto one = BigInt(1);
565 auto p_1 = p - one;
566 auto gcd_coprime = gcd(p_1, coprime);
567 if (!passes_sieve || gcd_coprime != one)
568 continue;
569 if (isPrime(p, rng, 64, true))
570 return p.move;
571 }
572 }
573 }
574
575 /// ditto
576 BigInt randomPrime()(RandomNumberGenerator rng,
577 size_t bits, const BigInt coprime = BigInt(1),
578 size_t equiv = 1, size_t modulo = 2)
579 {
580 return randomPrime(rng, bits, coprime, equiv, modulo);
581 }
582
583 /**
584 * Return a random 'safe' prime, of the form p=2*q+1 with q prime
585 * Params:
586 * rng = a random number generator
587 * bits = is how long the resulting prime should be
588 * Returns: prime randomly chosen from safe primes of length bits
589 */
590 BigInt randomSafePrime(RandomNumberGenerator rng, size_t bits)
591 {
592 if (bits <= 64)
593 throw new InvalidArgument("randomSafePrime: Can't make a prime of " ~ to!string(bits) ~ " bits");
594
595 BigInt p;
596 do
597 p = (randomPrime(rng, bits - 1) << 1) + 1;
598 while (!isPrime(p, rng, 64, true));
599 return p;
600 }
601
602 /**
603 * Generate DSA parameters using the FIPS 186 kosherizer
604 * Params:
605 * rng = a random number generator
606 * af = an algorithm factory
607 * p_out = where the prime p will be stored
608 * q_out = where the prime q will be stored
609 * pbits = how long p will be in bits
610 * qbits = how long q will be in bits
611 * Returns: random seed used to generate this parameter set
612 */
613 Vector!ubyte generateDsaPrimes(RandomNumberGenerator rng,
614 AlgorithmFactory af,
615 ref BigInt p_out, ref BigInt q_out,
616 size_t pbits, size_t qbits)
617 {
618 while (true)
619 {
620 Vector!ubyte seed = Vector!ubyte(qbits / 8);
621 rng.randomize(seed.ptr, seed.length);
622
623 if (generateDsaPrimes(rng, af, p_out, q_out, pbits, qbits, seed))
624 return seed;
625 }
626 }
627
628
629 /**
630 * Generate DSA parameters using the FIPS 186 kosherizer
631 * Params:
632 * rng = a random number generator
633 * af = an algorithm factory
634 * p_out = where the prime p will be stored
635 * q_out = where the prime q will be stored
636 * pbits = how long p will be in bits
637 * qbits = how long q will be in bits
638 * seed_c = the seed used to generate the parameters
639 * Returns: true if seed generated a valid DSA parameter set, otherwise
640 false. p_out and q_out are only valid if true was returned.
641 */
642 bool generateDsaPrimes()(RandomNumberGenerator rng,
643 AlgorithmFactory af,
644 ref BigInt p_out, ref BigInt q_out,
645 size_t pbits, size_t qbits,
646 auto const ref Vector!ubyte seed_c)
647 {
648 if (!fips1863ValidSize(pbits, qbits))
649 throw new InvalidArgument(
650 "FIPS 186-3 does not allow DSA domain parameters of " ~ to!string(pbits) ~ "/" ~ to!string(qbits) ~ " bits long");
651
652 if (seed_c.length * 8 < qbits)
653 throw new InvalidArgument("Generating a DSA parameter set with a " ~ to!string(qbits) ~
654 "long q requires a seed at least as many bits long");
655
656 Unique!HashFunction hash = af.makeHashFunction("SHA-" ~ to!string(qbits));
657
658 const size_t HASH_SIZE = hash.outputLength;
659
660 struct Seed
661 {
662 public:
663 this()(auto const ref Vector!ubyte s) { m_seed = s.dup(); }
664
665 ref T opCast(T : Vector!ubyte)() { return m_seed; }
666
667 alias m_seed this;
668
669 ref Seed opUnary(string op)()
670 if (op == "++")
671 {
672 for (size_t j = m_seed.length; j > 0; --j)
673 if (++m_seed[j-1])
674 break;
675 return this;
676 }
677 private:
678 Vector!ubyte m_seed;
679 }
680
681 Seed seed = Seed(seed_c);
682
683 q_out.binaryDecode(hash.process(seed));
684 q_out.setBit(qbits-1);
685 q_out.setBit(0);
686
687 if (!isPrime(q_out, rng))
688 return false;
689
690 const size_t n = (pbits-1) / (HASH_SIZE * 8), b = (pbits-1) % (HASH_SIZE * 8);
691
692 BigInt X;
693 Vector!ubyte V = Vector!ubyte(HASH_SIZE * (n+1));
694
695 foreach (size_t j; 0 .. 4096)
696 {
697 for (size_t k = 0; k <= n; ++k)
698 {
699 ++seed;
700 hash.update(seed);
701 hash.flushInto(&V[HASH_SIZE * (n-k)]);
702 }
703
704 X.binaryDecode(&V[HASH_SIZE - 1 - b/8], V.length - (HASH_SIZE - 1 - b/8));
705 X.setBit(pbits-1);
706
707 p_out = X - (X % (q_out*2) - 1);
708
709 if (p_out.bits() == pbits && isPrime(p_out, rng))
710 return true;
711 }
712 return false;
713 }
714
715 /*
716 * Check if this size is allowed by FIPS 186-3
717 */
718 bool fips1863ValidSize(size_t pbits, size_t qbits)
719 {
720 if (qbits == 160)
721 return (pbits == 512 || pbits == 768 || pbits == 1024);
722
723 if (qbits == 224)
724 return (pbits == 2048);
725
726 if (qbits == 256)
727 return (pbits == 2048 || pbits == 3072);
728
729 return false;
730 }
731
732 /*
733 * If the modulus is odd, then we can avoid computing A and C. This is
734 * a critical path algorithm in some instances and an odd modulus is
735 * the common case for crypto, so worth special casing. See note 14.64
736 * in Handbook of Applied Cryptography for more details.
737 */
738 BigInt inverseModOddModulus(const ref BigInt n, const ref BigInt mod)
739 {
740 BigInt u = mod.dup;
741 BigInt v = n.dup;
742 BigInt B = BigInt(0);
743 BigInt D = BigInt(1);
744
745 while (u.isNonzero())
746 {
747 const size_t u_zero_bits = lowZeroBits(u);
748 u >>= u_zero_bits;
749 foreach (size_t i; 0 .. u_zero_bits)
750 {
751 if (B.isOdd())
752 { B -= mod; }
753 B >>= 1;
754 }
755
756 const size_t v_zero_bits = lowZeroBits(v);
757 v >>= v_zero_bits;
758 foreach (size_t i; 0 .. v_zero_bits)
759 {
760 if (D.isOdd())
761 { D -= mod; }
762 D >>= 1;
763 }
764
765 if (u >= v) { u -= v; B -= D; }
766 else { v -= u; D -= B; }
767 }
768
769 if (v != 1)
770 return BigInt(0); // no modular inverse
771
772 while (D.isNegative()) D += mod;
773 while (D >= mod) D -= mod;
774
775 return D.move();
776 }
777
778 bool mrWitness(T : ModularReducer)(ref BigInt y,
779 auto ref T reducer_n,
780 const ref BigInt n_minus_1, size_t s)
781 {
782 if (y == 1 || y == n_minus_1)
783 return false;
784
785 foreach (size_t i; 1 .. s)
786 {
787 y = reducer_n.square(y);
788
789 if (y == 1) // found a non-trivial square root
790 return true;
791 if (y == n_minus_1) // -1, trivial square root, so give up
792 return false;
793 }
794
795 return true; // fails Fermat test
796 }
797
798 size_t mrTestIterations(size_t n_bits, size_t prob, bool random)
799 {
800 const size_t base = (prob + 2) / 2; // worst case 4^-t error rate
801
802 /*
803 * For randomly chosen numbers we can use the estimates from
804 * http://www.math.dartmouth.edu/~carlp/PDF/paper88.pdf
805 *
806 * These values are derived from the inequality for p(k,t) given on
807 * the second page.
808 */
809 if (random && prob <= 80)
810 {
811 if (n_bits >= 1536)
812 return 2; // < 2^-89
813 if (n_bits >= 1024)
814 return 4; // < 2^-89
815 if (n_bits >= 512)
816 return 5; // < 2^-80
817 if (n_bits >= 256)
818 return 11; // < 2^-80
819 }
820
821 return base;
822 }