1 #include "BigIntegerAlgorithms.hh"
3 BigUnsigned gcd(BigUnsigned a, BigUnsigned b) {
5 // Neat in-place alternating technique.
9 a.divideWithRemainder(b, trash);
12 b.divideWithRemainder(a, trash);
16 void extendedEuclidean(BigInteger m, BigInteger n,
17 BigInteger &g, BigInteger &r, BigInteger &s) {
18 if (&g == &r || &g == &s || &r == &s)
19 throw "BigInteger extendedEuclidean: Outputs are aliased";
20 BigInteger r1(1), s1(0), r2(0), s2(1), q;
22 * r1*m + s1*n == m(orig)
23 * r2*m + s2*n == n(orig) */
26 r = r1; s = s1; g = m;
29 m.divideWithRemainder(n, q);
30 r1 -= q*r2; s1 -= q*s2;
33 r = r2; s = s2; g = n;
36 n.divideWithRemainder(m, q);
37 r2 -= q*r1; s2 -= q*s1;
41 BigUnsigned modinv(const BigInteger &x, const BigUnsigned &n) {
43 extendedEuclidean(x, n, g, r, s);
45 // r*x + s*n == 1, so r*x === 1 (mod n), so r is the answer.
46 return (r % n).getMagnitude(); // (r % n) will be nonnegative
48 throw "BigInteger modinv: x and n have a common factor";
51 BigUnsigned modexp(const BigInteger &base, const BigUnsigned &exponent,
52 const BigUnsigned &modulus) {
53 BigUnsigned ans = 1, base2 = (base % modulus).getMagnitude();
54 BigUnsigned::Index i = exponent.getLength();
55 // For each block of the exponent, most to least significant...
58 BigUnsigned::Blk eb = exponent.getBlock(i);
59 // For each bit, most to least significant...
60 for (BigUnsigned::Blk mask = ~((~BigUnsigned::Blk(0)) >> 1);
61 mask != 0; mask >>= 1) {
62 // Square and maybe multiply.