From 3d5b522d9e790f4b8718d9a0b3808ab6920c07fe Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Wed, 5 Jan 2022 12:53:39 +0000
Subject: [PATCH 01/28] Add modular exponentiation

---
 src/bigints.nim    | 25 +++++++++++++++++++++++++
 tests/tbigints.nim | 36 ++++++++++++++++++++++++++++++++++++
 2 files changed, 61 insertions(+)

diff --git a/src/bigints.nim b/src/bigints.nim
index 794bdb3..ef14fba 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -959,3 +959,28 @@ iterator `..<`*(a, b: BigInt): BigInt =
   while res < b:
     yield res
     inc res
+
+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]`.
+  if modulus < 0:
+    raise newException(ValueError, "modulus must be strictly positive")
+  if modulus == 0:
+    raise newException(DivByZeroDefect, "modulus must be nonzero")
+  if modulus == 1:
+    return zero
+  if exponent < 0:
+    # let baseInv = inverse_mod(base, modulus)
+    # return powmod(baseInv, -exponent, modulus)
+    raise newException(ValueError, "Not yet Implemented!")
+
+  var
+    exp = exponent
+    basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
+  result = one
+  while exp != zero:
+    if (exp and one) != zero:
+      result = (result * basePow) mod modulus
+    basePow = (basePow * basePow) mod modulus
+    exp = exp shr 1
+
diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index cd80bf6..63cd19f 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -379,5 +379,41 @@ template main() =
     doAssert pow(zero, 0) == one
     doAssert pow(zero, 1) == zero
 
+  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
+
+  block: # powmod with negative base
+    let a = "1986599".initBigInt
+    let p = "10230581".initBigInt
+    doAssert powmod(-a, 2.initBigInt, p) == "6199079".initBigInt
+
 static: main()
 main()

From daec81c7de12ef99baca7f2c8944932a612bec98 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 6 Jan 2022 22:05:26 +0000
Subject: [PATCH 02/28] Add modular inverse to complete modular exponentiation

---
 src/bigints.nim    | 29 ++++++++++++++++++++++++++---
 tests/tbigints.nim | 18 ++++++++++++++++++
 2 files changed, 44 insertions(+), 3 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index ef14fba..0e12ad4 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -960,6 +960,27 @@ iterator `..<`*(a, b: BigInt): BigInt =
     yield res
     inc res
 
+func invmod*(a, modulus: BigInt): BigInt =
+    ## Compute the modular inverse of `a` by Euclide's method. 
+    ## The return value is always in the range `[1, modulus-1]`
+    if a == 0:
+      raise newException(DivByZeroDefect, "0 has no modular inverse")
+    var
+      r0 = a
+      r1 = modulus
+      u = one
+      u1 = zero
+      q, rt, ut : BigInt
+    while r1 > 0:
+      q = r0 div r1
+      rt = r0
+      ut = u
+      r0 = r1
+      u = u1
+      r1 = rt - q * r1
+      u1 = ut - q * u1
+    result = ((u 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]`.
@@ -970,9 +991,11 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
   if modulus == 1:
     return zero
   if exponent < 0:
-    # let baseInv = inverse_mod(base, modulus)
-    # return powmod(baseInv, -exponent, modulus)
-    raise newException(ValueError, "Not yet Implemented!")
+    if base == zero:
+      return zero
+    else:
+      let baseInv = invmod(base, modulus)
+      return powmod(baseInv, -exponent, modulus)
 
   var
     exp = exponent
diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index 63cd19f..b609a25 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -379,6 +379,17 @@ template main() =
     doAssert pow(zero, 0) == one
     doAssert pow(zero, 1) == zero
 
+  block: # invmod
+    # with prime modulus
+    let a = "30292868".initBigInt
+    let p = "60449131".initBigInt # p is prime
+    doAssert invmod(a, p) == "51713091".initBigInt
+    # with composite modulus
+    let b = "2472018".initBigInt
+    let n = "3917515".initBigInt # 5 * 7 * 19 * 43 * 137
+    doAssert invmod(b, n) == "2622632".initBigInt
+
+
   block: # powmod
     let a = "30292868".initBigInt
     let p = "60449131".initBigInt # p is prime
@@ -415,5 +426,12 @@ template main() =
     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
+
+
 static: main()
 main()

From 7235a40caeaecf4c6a657f19ac704c12ee3c6808 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Sat, 8 Jan 2022 18:45:42 +0000
Subject: [PATCH 03/28] Fix indentation and docstring for modular inverse

---
 src/bigints.nim | 39 ++++++++++++++++++++-------------------
 1 file changed, 20 insertions(+), 19 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 0e12ad4..92b322b 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -961,25 +961,26 @@ iterator `..<`*(a, b: BigInt): BigInt =
     inc res
 
 func invmod*(a, modulus: BigInt): BigInt =
-    ## Compute the modular inverse of `a` by Euclide's method. 
-    ## The return value is always in the range `[1, modulus-1]`
-    if a == 0:
-      raise newException(DivByZeroDefect, "0 has no modular inverse")
-    var
-      r0 = a
-      r1 = modulus
-      u = one
-      u1 = zero
-      q, rt, ut : BigInt
-    while r1 > 0:
-      q = r0 div r1
-      rt = r0
-      ut = u
-      r0 = r1
-      u = u1
-      r1 = rt - q * r1
-      u1 = ut - q * u1
-    result = ((u mod modulus) + modulus) mod modulus
+  ## Compute the modular inverse of `a` modulo `modulus`. 
+  ## The return value is always in the range `[1, modulus-1]`
+  # Euclide's method
+  if a == 0:
+    raise newException(DivByZeroDefect, "0 has no modular inverse")
+  var
+    r0 = a
+    r1 = modulus
+    u = one
+    u1 = zero
+    q, rt, ut : BigInt
+  while r1 > 0:
+    q = r0 div r1
+    rt = r0
+    ut = u
+    r0 = r1
+    u = u1
+    r1 = rt - q * r1
+    u1 = ut - q * u1
+  result = ((u mod modulus) + modulus) mod modulus
 
 func powmod*(base, exponent, modulus: BigInt): BigInt =
   ## Compute modular exponentation of `base` with power `exponent` modulo `modulus`.

From 63979d6a471a633726970b5ccab775d2a3861c2d Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Sat, 8 Jan 2022 19:15:31 +0000
Subject: [PATCH 04/28] Test if integer is invertible modulo composite modulus

---
 src/bigints.nim | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/bigints.nim b/src/bigints.nim
index 92b322b..57801d6 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -980,6 +980,8 @@ func invmod*(a, modulus: BigInt): BigInt =
     u = u1
     r1 = rt - q * r1
     u1 = ut - q * u1
+  if r0 != one:
+    raise newException(ValueError, $a & " has no modular inverse")
   result = ((u mod modulus) + modulus) mod modulus
 
 func powmod*(base, exponent, modulus: BigInt): BigInt =

From 11c08ea58b6cc42972a09d6598218bf225996e18 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Mon, 10 Jan 2022 13:43:51 +0000
Subject: [PATCH 05/28] Changed to if..elif..else construct and add tests for
 zero

---
 src/bigints.nim    | 28 ++++++++++++++--------------
 tests/tbigints.nim |  2 ++
 2 files changed, 16 insertions(+), 14 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 57801d6..73a634e 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -989,24 +989,24 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
   ## The return value is always in the range `[0, modulus-1]`.
   if modulus < 0:
     raise newException(ValueError, "modulus must be strictly positive")
-  if modulus == 0:
+  elif modulus.isZero:
     raise newException(DivByZeroDefect, "modulus must be nonzero")
-  if modulus == 1:
+  elif modulus == 1:
     return zero
-  if exponent < 0:
-    if base == zero:
+  elif exponent < 0:
+    if base.isZero:
       return zero
     else:
       let baseInv = invmod(base, modulus)
       return powmod(baseInv, -exponent, modulus)
-
-  var
-    exp = exponent
-    basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
-  result = one
-  while exp != zero:
-    if (exp and one) != zero:
-      result = (result * basePow) mod modulus
-    basePow = (basePow * basePow) mod modulus
-    exp = exp shr 1
+  else:
+    var
+      exp = exponent
+      basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
+    result = one
+    while exp != zero:
+      if (exp and one) != zero:
+        result = (result * basePow) mod modulus
+      basePow = (basePow * basePow) mod modulus
+      exp = exp shr 1
 
diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index b609a25..e4d41a3 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -420,6 +420,8 @@ template main() =
     # 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

From 9f8811dccc0ab37069fff1bd6fd4b215844b1725 Mon Sep 17 00:00:00 2001
From: dlesnoff <54949944+dlesnoff@users.noreply.github.com>
Date: Mon, 10 Jan 2022 23:49:46 +0000
Subject: [PATCH 06/28] Optimize loop checks in powmod

Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com>
---
 src/bigints.nim | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 133e5da..17d7a6c 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1001,8 +1001,8 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
       exp = exponent
       basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
     result = one
-    while exp != zero:
-      if (exp and one) != zero:
+    while not exp.isZero:
+      if (exp.limbs[0] and 1) != 0:
         result = (result * basePow) mod modulus
       basePow = (basePow * basePow) mod modulus
       exp = exp shr 1

From 2eb864cdeb7354fd7573d5a1134a8cf49d64825c Mon Sep 17 00:00:00 2001
From: dlesnoff <54949944+dlesnoff@users.noreply.github.com>
Date: Mon, 10 Jan 2022 23:50:40 +0000
Subject: [PATCH 07/28] Update docstring of invmod

Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com>
---
 src/bigints.nim | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 17d7a6c..107982f 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -960,7 +960,7 @@ iterator `..<`*(a, b: BigInt): BigInt =
 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]`
-  # Euclide's method
+  # extended Euclidean algorithm
   if a == 0:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
   var

From c20bcf6a57f5a2837e4828ba025e0b636c222bcd Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Wed, 12 Jan 2022 22:29:55 +0000
Subject: [PATCH 08/28] Add runnableExamples and some recommandations

---
 src/bigints.nim | 56 +++++++++++++++++++++++++++++--------------------
 1 file changed, 33 insertions(+), 23 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 107982f..ac8d9e4 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -960,30 +960,43 @@ iterator `..<`*(a, b: BigInt): BigInt =
 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 a == 0:
+  if modulus < 0:
+    raise newException(ValueError, "modulus must be strictly positive")
+  elif modulus.isZero:
+    raise newException(DivByZeroDefect, "modulus must be nonzero")
+  elif a == 0:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
-  var
-    r0 = a
-    r1 = modulus
-    u = one
-    u1 = zero
-    q, rt, ut : BigInt
-  while r1 > 0:
-    q = r0 div r1
-    rt = r0
-    ut = u
-    r0 = r1
-    u = u1
-    r1 = rt - q * r1
-    u1 = ut - q * u1
-  if r0 != one:
-    raise newException(ValueError, $a & " has no modular inverse")
-  result = ((u mod modulus) + modulus) mod modulus
+  else:
+    var
+      r0 = a
+      r1 = modulus
+      u = one
+      u1 = zero
+      q, rt, ut : BigInt
+    while r1 > 0:
+      q = r0 div r1
+      rt = r0
+      ut = u
+      r0 = r1
+      u = u1
+      r1 = rt - q * r1
+      u1 = ut - q * u1
+    if r0 != one:
+      raise newException(ValueError, $a & " has no modular inverse")
+    result = ((u 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:
+    let two = 2.initBigInt
+    let three = 3.initBigInt
+    let seven = 7.initBigInt
+    assert powmod(two, three, seven) == one
   if modulus < 0:
     raise newException(ValueError, "modulus must be strictly positive")
   elif modulus.isZero:
@@ -991,11 +1004,8 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
   elif modulus == 1:
     return zero
   elif exponent < 0:
-    if base.isZero:
-      return zero
-    else:
-      let baseInv = invmod(base, modulus)
-      return powmod(baseInv, -exponent, modulus)
+    let baseInv = invmod(base, modulus)
+    return powmod(baseInv, -exponent, modulus)
   else:
     var
       exp = exponent

From 4428d18bc96bdb0f5a04f6a482f6ac76dcbe6423 Mon Sep 17 00:00:00 2001
From: dlesnoff <54949944+dlesnoff@users.noreply.github.com>
Date: Wed, 12 Jan 2022 22:30:25 +0000
Subject: [PATCH 09/28] Rename variables invmod

Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com>
---
 src/bigints.nim | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 107982f..cb89243 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -966,20 +966,20 @@ func invmod*(a, modulus: BigInt): BigInt =
   var
     r0 = a
     r1 = modulus
-    u = one
-    u1 = zero
-    q, rt, ut : BigInt
+    s0 = one
+    s1 = zero
   while r1 > 0:
-    q = r0 div r1
-    rt = r0
-    ut = u
+    let
+      q = r0 div r1
+      rk = r0 - q * r1
+      sk = s0 - q * s1
     r0 = r1
-    u = u1
-    r1 = rt - q * r1
-    u1 = ut - q * u1
+    r1 = rk
+    s0 = s1
+    s1 = sk
   if r0 != one:
     raise newException(ValueError, $a & " has no modular inverse")
-  result = ((u mod modulus) + modulus) mod 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`.

From 010845b04304e1e30024330ce583adb725bc5a4a Mon Sep 17 00:00:00 2001
From: dlesnoff <54949944+dlesnoff@users.noreply.github.com>
Date: Thu, 13 Jan 2022 00:02:47 +0000
Subject: [PATCH 10/28] Apply suggestions from code review

I see no opposition to these changes.

Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com>
---
 src/bigints.nim | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index d2b5df5..06f209d 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -964,11 +964,11 @@ func invmod*(a, modulus: BigInt): BigInt =
     invmod(3.initBigInt, 7.initBigInt) = 5.initBigInt
 
   # extended Euclidean algorithm
-  if modulus < 0:
-    raise newException(ValueError, "modulus must be strictly positive")
-  elif modulus.isZero:
+  if modulus.isZero:
     raise newException(DivByZeroDefect, "modulus must be nonzero")
-  elif a == 0:
+  elif modulus < 0:
+    raise newException(ValueError, "modulus must be strictly positive")
+  elif a.isZero:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
   else:
     var
@@ -993,14 +993,11 @@ 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:
-    let two = 2.initBigInt
-    let three = 3.initBigInt
-    let seven = 7.initBigInt
-    assert powmod(two, three, seven) == one
-  if modulus < 0:
-    raise newException(ValueError, "modulus must be strictly positive")
-  elif modulus.isZero:
+    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
   elif exponent < 0:

From 1a0dc03b2da3dc92d6a609d7d155ee5486c83c1f Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 13 Jan 2022 00:11:51 +0000
Subject: [PATCH 11/28] Avoid a function call if exponent is negative

---
 src/bigints.nim | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 06f209d..be0a154 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -966,7 +966,7 @@ func invmod*(a, modulus: BigInt): BigInt =
   # extended Euclidean algorithm
   if modulus.isZero:
     raise newException(DivByZeroDefect, "modulus must be nonzero")
-  elif modulus < 0:
+  elif modulus.isNegative:
     raise newException(ValueError, "modulus must be strictly positive")
   elif a.isZero:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
@@ -1000,13 +1000,11 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
     raise newException(ValueError, "modulus must be strictly positive")
   elif modulus == 1:
     return zero
-  elif exponent < 0:
-    let baseInv = invmod(base, modulus)
-    return powmod(baseInv, -exponent, modulus)
   else:
     var
-      exp = exponent
+      base = (if exponent.isNegative == false: base else: invmod(base, modulus))
       basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
+      exp = (if exponent.isNegative == false: exponent else: -exponent)
     result = one
     while not exp.isZero:
       if (exp.limbs[0] and 1) != 0:

From ceb2f35a91a90c9d13b3e3fdf71ebe5112e15220 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 13 Jan 2022 00:17:44 +0000
Subject: [PATCH 12/28] Remove ternary operators

---
 src/bigints.nim | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index be0a154..bf85655 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1002,13 +1002,17 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
     return zero
   else:
     var
-      base = (if exponent.isNegative == false: base else: invmod(base, modulus))
+      base = base
+      exponent = exponent
+    if exponent.isNegative == true:
+      base = invmod(base, modulus)
+      exponent = -exponent
+    var
       basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
-      exp = (if exponent.isNegative == false: exponent else: -exponent)
     result = one
-    while not exp.isZero:
-      if (exp.limbs[0] and 1) != 0:
+    while not exponent.isZero:
+      if (exponent.limbs[0] and 1) != 0:
         result = (result * basePow) mod modulus
       basePow = (basePow * basePow) mod modulus
-      exp = exp shr 1
+      exponent = exponent shr 1
 

From 60c6a6d0705ef56ccfaf1539f222f8e75b9c6aee Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 13 Jan 2022 11:32:15 +0000
Subject: [PATCH 13/28] Fix <0 comparison test

---
 src/bigints.nim | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index bf85655..a364961 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1004,7 +1004,7 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
     var
       base = base
       exponent = exponent
-    if exponent.isNegative == true:
+    if exponent < 0:
       base = invmod(base, modulus)
       exponent = -exponent
     var

From 2412bcee809570cad7f31d19c113878d27728e75 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Wed, 5 Jan 2022 12:53:39 +0000
Subject: [PATCH 14/28] Add modular exponentiation Merge with latests changes

---
 src/bigints.nim    | 25 +++++++++++++++++++++++++
 tests/tbigints.nim | 35 +++++++++++++++++++++++++++++++++++
 2 files changed, 60 insertions(+)

diff --git a/src/bigints.nim b/src/bigints.nim
index ec6aaf9..ae6f30b 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1005,3 +1005,28 @@ iterator `..<`*(a, b: BigInt): BigInt =
   while res < b:
     yield res
     inc res
+
+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]`.
+  if modulus < 0:
+    raise newException(ValueError, "modulus must be strictly positive")
+  if modulus == 0:
+    raise newException(DivByZeroDefect, "modulus must be nonzero")
+  if modulus == 1:
+    return zero
+  if exponent < 0:
+    # let baseInv = inverse_mod(base, modulus)
+    # return powmod(baseInv, -exponent, modulus)
+    raise newException(ValueError, "Not yet Implemented!")
+
+  var
+    exp = exponent
+    basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
+  result = one
+  while exp != zero:
+    if (exp and one) != zero:
+      result = (result * basePow) mod modulus
+    basePow = (basePow * basePow) mod modulus
+    exp = exp shr 1
+
diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index 279a58f..757d16c 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -448,6 +448,41 @@ proc main() =
     doAssert toSignedInt[int32](n) == none(int32)
     doAssert toSignedInt[int](n) == none(int)
 
+  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
+
+  block: # powmod with negative base
+    let a = "1986599".initBigInt
+    let p = "10230581".initBigInt
+    doAssert powmod(-a, 2.initBigInt, p) == "6199079".initBigInt
 
 static: main()
 main()

From 189b061e3d1f7b0bc367c377dbadb32c22ae512b Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Sun, 16 Jan 2022 10:33:46 +0000
Subject: [PATCH 15/28] Rebase tests

---
 src/bigints.nim    | 29 ++++++++++++++++++++++++++---
 tests/tbigints.nim | 17 +++++++++++++++++
 2 files changed, 43 insertions(+), 3 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index ae6f30b..05a3320 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1006,6 +1006,27 @@ iterator `..<`*(a, b: BigInt): BigInt =
     yield res
     inc res
 
+func invmod*(a, modulus: BigInt): BigInt =
+    ## Compute the modular inverse of `a` by Euclide's method. 
+    ## The return value is always in the range `[1, modulus-1]`
+    if a == 0:
+      raise newException(DivByZeroDefect, "0 has no modular inverse")
+    var
+      r0 = a
+      r1 = modulus
+      u = one
+      u1 = zero
+      q, rt, ut : BigInt
+    while r1 > 0:
+      q = r0 div r1
+      rt = r0
+      ut = u
+      r0 = r1
+      u = u1
+      r1 = rt - q * r1
+      u1 = ut - q * u1
+    result = ((u 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]`.
@@ -1016,9 +1037,11 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
   if modulus == 1:
     return zero
   if exponent < 0:
-    # let baseInv = inverse_mod(base, modulus)
-    # return powmod(baseInv, -exponent, modulus)
-    raise newException(ValueError, "Not yet Implemented!")
+    if base == zero:
+      return zero
+    else:
+      let baseInv = invmod(base, modulus)
+      return powmod(baseInv, -exponent, modulus)
 
   var
     exp = exponent
diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index 757d16c..b36e347 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -447,6 +447,16 @@ proc main() =
     doAssert toSignedInt[int8](n) == none(int8)
     doAssert toSignedInt[int32](n) == none(int32)
     doAssert toSignedInt[int](n) == none(int)
+=======
+  block: # invmod
+    # with prime modulus
+    let a = "30292868".initBigInt
+    let p = "60449131".initBigInt # p is prime
+    doAssert invmod(a, p) == "51713091".initBigInt
+    # with composite modulus
+    let b = "2472018".initBigInt
+    let n = "3917515".initBigInt # 5 * 7 * 19 * 43 * 137
+    doAssert invmod(b, n) == "2622632".initBigInt
 
   block: # powmod
     let a = "30292868".initBigInt
@@ -484,5 +494,12 @@ proc main() =
     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
+
+
 static: main()
 main()

From e3a38237796b43171b40daaa28063870327f48c6 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Sat, 8 Jan 2022 18:45:42 +0000
Subject: [PATCH 16/28] Fix indentation and docstring for modular inverse

---
 src/bigints.nim | 39 ++++++++++++++++++++-------------------
 1 file changed, 20 insertions(+), 19 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 05a3320..2ab0552 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1007,25 +1007,26 @@ iterator `..<`*(a, b: BigInt): BigInt =
     inc res
 
 func invmod*(a, modulus: BigInt): BigInt =
-    ## Compute the modular inverse of `a` by Euclide's method. 
-    ## The return value is always in the range `[1, modulus-1]`
-    if a == 0:
-      raise newException(DivByZeroDefect, "0 has no modular inverse")
-    var
-      r0 = a
-      r1 = modulus
-      u = one
-      u1 = zero
-      q, rt, ut : BigInt
-    while r1 > 0:
-      q = r0 div r1
-      rt = r0
-      ut = u
-      r0 = r1
-      u = u1
-      r1 = rt - q * r1
-      u1 = ut - q * u1
-    result = ((u mod modulus) + modulus) mod modulus
+  ## Compute the modular inverse of `a` modulo `modulus`. 
+  ## The return value is always in the range `[1, modulus-1]`
+  # Euclide's method
+  if a == 0:
+    raise newException(DivByZeroDefect, "0 has no modular inverse")
+  var
+    r0 = a
+    r1 = modulus
+    u = one
+    u1 = zero
+    q, rt, ut : BigInt
+  while r1 > 0:
+    q = r0 div r1
+    rt = r0
+    ut = u
+    r0 = r1
+    u = u1
+    r1 = rt - q * r1
+    u1 = ut - q * u1
+  result = ((u mod modulus) + modulus) mod modulus
 
 func powmod*(base, exponent, modulus: BigInt): BigInt =
   ## Compute modular exponentation of `base` with power `exponent` modulo `modulus`.

From e52effc5a05c1b817381fe3251effb5ff3e3aa3a Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Sat, 8 Jan 2022 19:15:31 +0000
Subject: [PATCH 17/28] Test if integer is invertible modulo composite modulus

---
 src/bigints.nim | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/bigints.nim b/src/bigints.nim
index 2ab0552..0544383 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1026,6 +1026,8 @@ func invmod*(a, modulus: BigInt): BigInt =
     u = u1
     r1 = rt - q * r1
     u1 = ut - q * u1
+  if r0 != one:
+    raise newException(ValueError, $a & " has no modular inverse")
   result = ((u mod modulus) + modulus) mod modulus
 
 func powmod*(base, exponent, modulus: BigInt): BigInt =

From 8cced8ae8bfb9b9dc0786f96d05b1afad5e8cd6e Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Mon, 10 Jan 2022 13:43:51 +0000
Subject: [PATCH 18/28] Changed to if..elif..else construct and add tests for
 zero

---
 src/bigints.nim    | 28 ++++++++++++++--------------
 tests/tbigints.nim |  2 ++
 2 files changed, 16 insertions(+), 14 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 0544383..574a07f 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1035,24 +1035,24 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
   ## The return value is always in the range `[0, modulus-1]`.
   if modulus < 0:
     raise newException(ValueError, "modulus must be strictly positive")
-  if modulus == 0:
+  elif modulus.isZero:
     raise newException(DivByZeroDefect, "modulus must be nonzero")
-  if modulus == 1:
+  elif modulus == 1:
     return zero
-  if exponent < 0:
-    if base == zero:
+  elif exponent < 0:
+    if base.isZero:
       return zero
     else:
       let baseInv = invmod(base, modulus)
       return powmod(baseInv, -exponent, modulus)
-
-  var
-    exp = exponent
-    basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
-  result = one
-  while exp != zero:
-    if (exp and one) != zero:
-      result = (result * basePow) mod modulus
-    basePow = (basePow * basePow) mod modulus
-    exp = exp shr 1
+  else:
+    var
+      exp = exponent
+      basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
+    result = one
+    while exp != zero:
+      if (exp and one) != zero:
+        result = (result * basePow) mod modulus
+      basePow = (basePow * basePow) mod modulus
+      exp = exp shr 1
 
diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index b36e347..96e20bc 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -488,6 +488,8 @@ proc main() =
     # 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

From 4bbfa5b8aa6476f0a98b52d0cee1c348ac13e605 Mon Sep 17 00:00:00 2001
From: dlesnoff <54949944+dlesnoff@users.noreply.github.com>
Date: Mon, 10 Jan 2022 23:49:46 +0000
Subject: [PATCH 19/28] Optimize loop checks in powmod

Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com>
---
 src/bigints.nim | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 574a07f..3989420 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1050,8 +1050,8 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
       exp = exponent
       basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
     result = one
-    while exp != zero:
-      if (exp and one) != zero:
+    while not exp.isZero:
+      if (exp.limbs[0] and 1) != 0:
         result = (result * basePow) mod modulus
       basePow = (basePow * basePow) mod modulus
       exp = exp shr 1

From 6bff7166987d82f22f95386d86b2cadf26509753 Mon Sep 17 00:00:00 2001
From: dlesnoff <54949944+dlesnoff@users.noreply.github.com>
Date: Mon, 10 Jan 2022 23:50:40 +0000
Subject: [PATCH 20/28] Update docstring of invmod

Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com>
---
 src/bigints.nim | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 3989420..db4350f 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1009,7 +1009,7 @@ iterator `..<`*(a, b: BigInt): BigInt =
 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]`
-  # Euclide's method
+  # extended Euclidean algorithm
   if a == 0:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
   var

From 961466f742f7abcf922705c23312608ea0e98324 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Wed, 12 Jan 2022 22:29:55 +0000
Subject: [PATCH 21/28] Add runnableExamples and some recommandations

---
 src/bigints.nim | 56 +++++++++++++++++++++++++++++--------------------
 1 file changed, 33 insertions(+), 23 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index db4350f..908f5a2 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1009,30 +1009,43 @@ iterator `..<`*(a, b: BigInt): BigInt =
 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 a == 0:
+  if modulus < 0:
+    raise newException(ValueError, "modulus must be strictly positive")
+  elif modulus.isZero:
+    raise newException(DivByZeroDefect, "modulus must be nonzero")
+  elif a == 0:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
-  var
-    r0 = a
-    r1 = modulus
-    u = one
-    u1 = zero
-    q, rt, ut : BigInt
-  while r1 > 0:
-    q = r0 div r1
-    rt = r0
-    ut = u
-    r0 = r1
-    u = u1
-    r1 = rt - q * r1
-    u1 = ut - q * u1
-  if r0 != one:
-    raise newException(ValueError, $a & " has no modular inverse")
-  result = ((u mod modulus) + modulus) mod modulus
+  else:
+    var
+      r0 = a
+      r1 = modulus
+      u = one
+      u1 = zero
+      q, rt, ut : BigInt
+    while r1 > 0:
+      q = r0 div r1
+      rt = r0
+      ut = u
+      r0 = r1
+      u = u1
+      r1 = rt - q * r1
+      u1 = ut - q * u1
+    if r0 != one:
+      raise newException(ValueError, $a & " has no modular inverse")
+    result = ((u 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:
+    let two = 2.initBigInt
+    let three = 3.initBigInt
+    let seven = 7.initBigInt
+    assert powmod(two, three, seven) == one
   if modulus < 0:
     raise newException(ValueError, "modulus must be strictly positive")
   elif modulus.isZero:
@@ -1040,11 +1053,8 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
   elif modulus == 1:
     return zero
   elif exponent < 0:
-    if base.isZero:
-      return zero
-    else:
-      let baseInv = invmod(base, modulus)
-      return powmod(baseInv, -exponent, modulus)
+    let baseInv = invmod(base, modulus)
+    return powmod(baseInv, -exponent, modulus)
   else:
     var
       exp = exponent

From a07a52aaaf1c50fb29f4a498c9479e4e4121eed5 Mon Sep 17 00:00:00 2001
From: dlesnoff <54949944+dlesnoff@users.noreply.github.com>
Date: Thu, 13 Jan 2022 00:02:47 +0000
Subject: [PATCH 22/28] Apply suggestions from code review

I see no opposition to these changes.

Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com>
---
 src/bigints.nim | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 908f5a2..d7eccc4 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1013,11 +1013,11 @@ func invmod*(a, modulus: BigInt): BigInt =
     invmod(3.initBigInt, 7.initBigInt) = 5.initBigInt
 
   # extended Euclidean algorithm
-  if modulus < 0:
-    raise newException(ValueError, "modulus must be strictly positive")
-  elif modulus.isZero:
+  if modulus.isZero:
     raise newException(DivByZeroDefect, "modulus must be nonzero")
-  elif a == 0:
+  elif modulus < 0:
+    raise newException(ValueError, "modulus must be strictly positive")
+  elif a.isZero:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
   else:
     var
@@ -1042,14 +1042,11 @@ 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:
-    let two = 2.initBigInt
-    let three = 3.initBigInt
-    let seven = 7.initBigInt
-    assert powmod(two, three, seven) == one
-  if modulus < 0:
-    raise newException(ValueError, "modulus must be strictly positive")
-  elif modulus.isZero:
+    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
   elif exponent < 0:

From 8371690e5fe6be1334ad9edb986d0e4cc0a1f1b0 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 13 Jan 2022 00:11:51 +0000
Subject: [PATCH 23/28] Avoid a function call if exponent is negative

---
 src/bigints.nim | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index d7eccc4..0c0f235 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1015,7 +1015,7 @@ func invmod*(a, modulus: BigInt): BigInt =
   # extended Euclidean algorithm
   if modulus.isZero:
     raise newException(DivByZeroDefect, "modulus must be nonzero")
-  elif modulus < 0:
+  elif modulus.isNegative:
     raise newException(ValueError, "modulus must be strictly positive")
   elif a.isZero:
     raise newException(DivByZeroDefect, "0 has no modular inverse")
@@ -1049,13 +1049,11 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
     raise newException(ValueError, "modulus must be strictly positive")
   elif modulus == 1:
     return zero
-  elif exponent < 0:
-    let baseInv = invmod(base, modulus)
-    return powmod(baseInv, -exponent, modulus)
   else:
     var
-      exp = exponent
+      base = (if exponent.isNegative == false: base else: invmod(base, modulus))
       basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
+      exp = (if exponent.isNegative == false: exponent else: -exponent)
     result = one
     while not exp.isZero:
       if (exp.limbs[0] and 1) != 0:

From cc1ba71e32f8141d0d2271852267b8fdba9c22e8 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 13 Jan 2022 00:17:44 +0000
Subject: [PATCH 24/28] Remove ternary operators

---
 src/bigints.nim | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 0c0f235..45bfc2f 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1051,13 +1051,17 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
     return zero
   else:
     var
-      base = (if exponent.isNegative == false: base else: invmod(base, modulus))
+      base = base
+      exponent = exponent
+    if exponent.isNegative == true:
+      base = invmod(base, modulus)
+      exponent = -exponent
+    var
       basePow = ((base mod modulus) + modulus) mod modulus # Base stays in [0, m-1]
-      exp = (if exponent.isNegative == false: exponent else: -exponent)
     result = one
-    while not exp.isZero:
-      if (exp.limbs[0] and 1) != 0:
+    while not exponent.isZero:
+      if (exponent.limbs[0] and 1) != 0:
         result = (result * basePow) mod modulus
       basePow = (basePow * basePow) mod modulus
-      exp = exp shr 1
+      exponent = exponent shr 1
 

From c4a4710a1d6afe49787033a8ddc78ef46033477e Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 13 Jan 2022 11:32:15 +0000
Subject: [PATCH 25/28] Fix <0 comparison test

---
 src/bigints.nim | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/bigints.nim b/src/bigints.nim
index 45bfc2f..84b57c2 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -1053,7 +1053,7 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
     var
       base = base
       exponent = exponent
-    if exponent.isNegative == true:
+    if exponent < 0:
       base = invmod(base, modulus)
       exponent = -exponent
     var

From da464b1eab6d82cd70e3647cd7bf61d3f5805d46 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Sun, 16 Jan 2022 10:39:46 +0000
Subject: [PATCH 26/28] remove a merge conflict

---
 tests/tbigints.nim | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index 96e20bc..78ace9f 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -447,7 +447,7 @@ proc main() =
     doAssert toSignedInt[int8](n) == none(int8)
     doAssert toSignedInt[int32](n) == none(int32)
     doAssert toSignedInt[int](n) == none(int)
-=======
+
   block: # invmod
     # with prime modulus
     let a = "30292868".initBigInt

From c057c6e1ebb4e4960b021f6bf46b67aec89f3fd3 Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Wed, 19 Jan 2022 18:41:58 +0000
Subject: [PATCH 27/28] Remove modexp example

---
 examples/rc_modexp.nim | 20 --------------------
 1 file changed, 20 deletions(-)
 delete mode 100644 examples/rc_modexp.nim

diff --git a/examples/rc_modexp.nim b/examples/rc_modexp.nim
deleted file mode 100644
index a4c7fec..0000000
--- a/examples/rc_modexp.nim
+++ /dev/null
@@ -1,20 +0,0 @@
-# Solution for http://rosettacode.org/wiki/Modular_exponentiation
-
-import bigints
-
-proc powmod(b, e, m: BigInt): BigInt =
-  assert e >= 0
-  var e = e
-  var b = b
-  result = initBigInt(1)
-  while e > 0:
-    if e mod 2 == 1:
-      result = (result * b) mod m
-    e = e div 2
-    b = (b.pow 2) mod m
-
-var
-  a = initBigInt("2988348162058574136915891421498819466320163312926952423791023078876139")
-  b = initBigInt("2351399303373464486466122544523690094744975233415544072992656881240319")
-
-echo powmod(a, b, 10.pow 40)

From 65c5766bd1f45df53f1a44ec9cbf695ee62b988f Mon Sep 17 00:00:00 2001
From: Dimitri Lesnoff <dimitri.lesnoff@gmail.com>
Date: Thu, 20 Jan 2022 22:25:45 +0000
Subject: [PATCH 28/28] Add more tests for invmod

---
 tests/tbigints.nim | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/tests/tbigints.nim b/tests/tbigints.nim
index 9ddb20c..a1ce0a7 100644
--- a/tests/tbigints.nim
+++ b/tests/tbigints.nim
@@ -391,10 +391,22 @@ proc main() =
     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) == "2502552".initBigInt
+
+    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)