BigIntegerAlgorithms.cc

   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(orig) + s1*n(orig) == m(current)
  23          * r2*m(orig) + s2*n(orig) == n(current) */
  24         for (;;) {
  25                 if (n.isZero()) {
  26                         r = r1; s = s1; g = m;
  27                         return;
  28                 }
  29                 // Subtract q times the second invariant from the first invariant.
  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                 }
  37                 // Subtract q times the first invariant from the second invariant.
  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();
  56         BigUnsigned::Index i = exponent.bitLength();
  57         // For each bit of the exponent, most to least significant...
  58         while (i > 0) {
  59                 i--;
  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;
  66                         ans %= modulus;
  67                 }
  68         }
  69         return ans;
  70 }