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 + 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 }