Commit | Line | Data |
---|---|---|
0afe80d5 MM |
1 | #include "BigIntegerAlgorithms.hh" |
2 | ||
3 | BigUnsigned gcd(BigUnsigned a, BigUnsigned b) { | |
4 | BigUnsigned trash; | |
5 | // Neat in-place alternating technique. | |
6 | for (;;) { | |
7 | if (b.isZero()) | |
8 | return a; | |
9 | a.divideWithRemainder(b, trash); | |
10 | if (a.isZero()) | |
11 | return b; | |
12 | b.divideWithRemainder(a, trash); | |
13 | } | |
14 | } | |
15 | ||
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; | |
21 | /* Invariants: | |
22 | * r1*m + s1*n == m(orig) | |
23 | * r2*m + s2*n == n(orig) */ | |
24 | for (;;) { | |
25 | if (n.isZero()) { | |
26 | r = r1; s = s1; g = m; | |
27 | return; | |
28 | } | |
29 | m.divideWithRemainder(n, q); | |
30 | r1 -= q*r2; s1 -= q*s2; | |
31 | ||
32 | if (m.isZero()) { | |
33 | r = r2; s = s2; g = n; | |
34 | return; | |
35 | } | |
36 | n.divideWithRemainder(m, q); | |
37 | r2 -= q*r1; s2 -= q*s1; | |
38 | } | |
39 | } | |
40 | ||
41 | BigUnsigned modinv(const BigInteger &x, const BigUnsigned &n) { | |
42 | BigInteger g, r, s; | |
43 | extendedEuclidean(x, n, g, r, s); | |
44 | if (g == 1) | |
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 | |
47 | else | |
48 | throw "BigInteger modinv: x and n have a common factor"; | |
49 | } | |
50 | ||
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... | |
56 | while (i > 0) { | |
57 | i--; | |
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. | |
63 | ans *= ans; | |
64 | ans %= modulus; | |
65 | if (eb & mask) { | |
66 | ans *= base2; | |
67 | ans %= modulus; | |
68 | } | |
69 | } | |
70 | } | |
71 | return ans; | |
72 | } |