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: | |
b2a47e1a MM |
22 | * r1*m(orig) + s1*n(orig) == m(current) |
23 | * r2*m(orig) + s2*n(orig) == n(current) */ | |
0afe80d5 MM |
24 | for (;;) { |
25 | if (n.isZero()) { | |
26 | r = r1; s = s1; g = m; | |
27 | return; | |
28 | } | |
b2a47e1a | 29 | // Subtract q times the second invariant from the first invariant. |
0afe80d5 MM |
30 | m.divideWithRemainder(n, q); |
31 | r1 -= q*r2; s1 -= q*s2; | |
32 | ||
33 | if (m.isZero()) { | |
34 | r = r2; s = s2; g = n; | |
35 | return; | |
36 | } | |
b2a47e1a | 37 | // Subtract q times the first invariant from the second invariant. |
0afe80d5 MM |
38 | n.divideWithRemainder(m, q); |
39 | r2 -= q*r1; s2 -= q*s1; | |
40 | } | |
41 | } | |
42 | ||
43 | BigUnsigned modinv(const BigInteger &x, const BigUnsigned &n) { | |
44 | BigInteger g, r, s; | |
45 | extendedEuclidean(x, n, g, r, s); | |
46 | if (g == 1) | |
47 | // r*x + s*n == 1, so r*x === 1 (mod n), so r is the answer. | |
48 | return (r % n).getMagnitude(); // (r % n) will be nonnegative | |
49 | else | |
50 | throw "BigInteger modinv: x and n have a common factor"; | |
51 | } | |
52 | ||
53 | BigUnsigned modexp(const BigInteger &base, const BigUnsigned &exponent, | |
54 | const BigUnsigned &modulus) { | |
55 | BigUnsigned ans = 1, base2 = (base % modulus).getMagnitude(); | |
b2a47e1a MM |
56 | BigUnsigned::Index i = exponent.bitLength(); |
57 | // For each bit of the exponent, most to least significant... | |
0afe80d5 MM |
58 | while (i > 0) { |
59 | i--; | |
b2a47e1a MM |
60 | // Square. |
61 | ans *= ans; | |
62 | ans %= modulus; | |
63 | // And multiply if the bit is a 1. | |
64 | if (exponent.getBit(i)) { | |
65 | ans *= base2; | |
0afe80d5 | 66 | ans %= modulus; |
0afe80d5 MM |
67 | } |
68 | } | |
69 | return ans; | |
70 | } |