Add modular exponentiation and modular inverse, closes #45 (#65)

Author: dlesnoff

examples/rc_modexp.nim

@@ -1,6 +1,5 @@
 ## Arbitrary precision integers.
 
-
 import std/[algorithm, bitops, math, options]
 
 type
@@ -1019,3 +1018,63 @@ iterator `..<`*(a, b: BigInt): BigInt =
   while res < b:
     yield res
     inc res
+
+func invmod*(a, modulus: BigInt): BigInt =
+  ## Compute the modular inverse of `a` modulo `modulus`. 
+  ## The return value is always in the range `[1, modulus-1]`
+  runnableExamples:
+    invmod(3.initBigInt, 7.initBigInt) = 5.initBigInt
+
+  # extended Euclidean algorithm
+  if modulus.isZero:
+    raise newException(DivByZeroDefect, "modulus must be nonzero")
+  elif modulus.isNegative:
+    raise newException(ValueError, "modulus must be strictly positive")
+  elif a.isZero:
+    raise newException(DivByZeroDefect, "0 has no modular inverse")
+  else:
+    var
+      r0 = ((a mod modulus) + modulus) mod modulus
+      r1 = modulus
+      s0 = one
+      s1 = zero
+    while r1 > 0:
+      let
+        q = r0 div r1
+        rk = r0 - q * r1
+        sk = s0 - q * s1
+      r0 = r1
+      r1 = rk
+      s0 = s1
+      s1 = sk
+    if r0 != one:
+      raise newException(ValueError, $a & " has no modular inverse modulo " & $modulus)
+    result = ((s0 mod modulus) + modulus) mod modulus
+
+func powmod*(base, exponent, modulus: BigInt): BigInt =
+  ## Compute modular exponentation of `base` with power `exponent` modulo `modulus`.
+  ## The return value is always in the range `[0, modulus-1]`.
+  runnableExamples:
+    assert powmod(2.initBigInt, 3.initBigInt, 7.initBigInt) == 1.initBigInt
+  if modulus.isZero:
+    raise newException(DivByZeroDefect, "modulus must be nonzero")
+  elif modulus.isNegative:
+    raise newException(ValueError, "modulus must be strictly positive")
+  elif modulus == 1:
+    return zero
+  else:
+    var
+      base = base
+      exponent = exponent
+    if exponent < 0:
+      base = invmod(base, modulus)
+      exponent = -exponent
+    var
+      basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
+    result = one
+    while not exponent.isZero:
+      if (exponent.limbs[0] and 1) != 0:
+        result = (result * basePow) mod modulus
+      basePow = (basePow * basePow) mod modulus
+      exponent = exponent shr 1
+

src/bigints.nim

@@ -380,6 +380,82 @@ proc main() =
     doAssert pow(zero, 0) == one
     doAssert pow(zero, 1) == zero
 
+  block: # invmod
+    # with prime modulus
+    let a = "30292868".initBigInt
+    let b = "48810860".initBigInt
+    let p = "60449131".initBigInt # p is prime
+    doAssert invmod(a, p) == "51713091".initBigInt
+    doAssert invmod(-b, p) == "31975542".initBigInt
+    # with composite modulus
+    let c = "2472018".initBigInt
+    let n = "3917515".initBigInt # 5 * 7 * 19 * 43 * 137
+    let d = "1831482".initBigInt
+    let e = "2502552".initBigInt
+    let f = "2086033".initBigInt
+    let h = "1414963".initBigInt
+    doAssert invmod(c, n) == "2622632".initBigInt
+    doAssert invmod(one, n) == one
+    doAssert invmod(n-one, n) == n-one
+
+    doAssert invmod( d, n) == h
+    doAssert invmod(-d, n) == e
+    doAssert invmod( f, n) == e
+    doAssert invmod(-f, n) == h
+    doAssert invmod( e, n) == f
+    doAssert invmod(-e, n) == d
+    doAssert invmod( h, n) == d
+    doAssert invmod(-h, n) == f
+
+    doAssertRaises(DivByZeroDefect): discard invmod(zero, n)
+    doAssertRaises(DivByZeroDefect): discard invmod(one, zero)
+    doAssertRaises(ValueError): discard invmod(one, -7.initBigInt)
+    doAssertRaises(ValueError): discard invmod(3.initBigInt, 18.initBigInt) # 3 is not invertible since gcd(3, 18) = 3 != 1
+
+  block: # powmod
+    let a = "30292868".initBigInt
+    let p = "60449131".initBigInt # p is prime
+    let two = 2.initBigInt
+    doAssert powmod(a, two, p) == "25760702".initBigInt
+    # Fermat's little theorem: a^p ≡ a mod p
+    doAssert powmod(a, p, p) == a
+    # Euler's identity a^(p-1) \equiv 1 \bmod p
+    doAssert powmod(a, p - one, p) == one
+    # We can invert a using Euler's identity / Fermat's little theorem
+    doAssert powmod(a, p - two, p) == "51713091".initBigInt
+    # We can reduce the exponent modulo phi(p) = p - 1, since p is prime
+    doAssert powmod(a, 2.initBigInt*p, p) == (a * a mod p)
+
+    let p2 = 761.initBigInt
+    var a2 = 1.initBigInt
+    # Fermat's little theorem: a^p ≡ a mod p
+    while a2 < p2:
+      doAssert powmod(a2, p2, p2) == a2
+      a2.inc
+
+  block: # Composite modulus
+    let a = "2472018".initBigInt
+    let n = "3917515".initBigInt # 5 * 7 * 19 * 43 * 137
+    let euler_phi = "2467584".initBigInt
+    doAssert powmod(a, 52.initBigInt, n) == "2305846".initBigInt
+    doAssert powmod(a, euler_phi, n) == one
+    # Edge cases
+    doAssert powmod(a, one, n) == a
+    doAssert powmod(a, zero, n) == one
+    doAssert powmod(zero, zero, n) == one
+    doAssert powmod(zero, one, n) == zero
+
+  block: # powmod with negative base
+    let a = "1986599".initBigInt
+    let p = "10230581".initBigInt
+    doAssert powmod(-a, 2.initBigInt, p) == "6199079".initBigInt
+
+  block: # powmod with negative exponent
+    let a = "1912".initBigInt
+    let p = "5297".initBigInt
+    doAssert powmod(a, -1.initBigInt, p) == "1460".initBigInt
+    doAssert powmod(a, one-p, p) == one
+
   block: # div/mod
     doAssertRaises(DivByZeroDefect): discard one div zero
     doAssertRaises(DivByZeroDefect): discard one mod zero
@@ -455,6 +531,5 @@ proc main() =
     doAssert pred(a, 3) == initBigInt(4)
     doAssert succ(a, 3) == initBigInt(10)
 
-
 static: main()
 main()