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(orig) + s1*n(orig) == m(current)
23 * r2*m(orig) + s2*n(orig) == n(current) */
26 r
= r1
; s
= s1
; g
= m
;
29 // Subtract q times the second invariant from the first invariant.
30 m
.divideWithRemainder(n
, q
);
31 r1
-= q
*r2
; s1
-= q
*s2
;
34 r
= r2
; s
= s2
; g
= n
;
37 // Subtract q times the first invariant from the second invariant.
38 n
.divideWithRemainder(m
, q
);
39 r2
-= q
*r1
; s2
-= q
*s1
;
43 BigUnsigned
modinv(const BigInteger
&x
, const BigUnsigned
&n
) {
45 extendedEuclidean(x
, n
, g
, r
, s
);
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
50 throw "BigInteger modinv: x and n have a common factor";
53 BigUnsigned
modexp(const BigInteger
&base
, const BigUnsigned
&exponent
,
54 const BigUnsigned
&modulus
) {
55 BigUnsigned ans
= 1, base2
= (base
% modulus
).getMagnitude();
56 BigUnsigned::Index i
= exponent
.bitLength();
57 // For each bit of the exponent, most to least significant...
63 // And multiply if the bit is a 1.
64 if (exponent
.getBit(i
)) {