From 67d038c650a1634e08c1bcaf5dda15a65fc3d433 Mon Sep 17 00:00:00 2001
From: mratsim <mamy_github@numforge.co>
Date: Sun, 2 Dec 2018 03:45:13 +0100
Subject: [PATCH] Implement modular inverse mod 2^k for Montgomery
 multiplication

---
 constantine/field_fp.nim | 66 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 65 insertions(+), 1 deletion(-)

diff --git a/constantine/field_fp.nim b/constantine/field_fp.nim
index e77aa7a..753115e 100644
--- a/constantine/field_fp.nim
+++ b/constantine/field_fp.nim
@@ -17,11 +17,16 @@ import
   ./word_types, ./bigints
 
 type
-  Fp[P: static BigInt] = object
+  Fp*[P: static BigInt] = object
     ## P is a prime number
     ## All operations on a field are modulo P
     value: type(P)
 
+  Montgomery*[M: static BigInt] = object
+    ## All operations in the Montgomery domain
+    ## are modulo M. M **must** be odd
+    value: type(M)
+
 # ############################################################
 #
 #                         Aliases
@@ -58,3 +63,62 @@ func `+`*(a, b: Fp): Fp =
   var ctl = add(result, b, True)
   ctl = ctl or not sub(result, Fp.P, False)
   sub(result, Fp.P, ctl)
+
+# ############################################################
+#
+#                Montgomery domain primitives
+#
+# ############################################################
+
+from bitops import fastLog2
+  # This will only be used at compile-time
+  # so no constant-time worries (it is constant-time if using the De Bruijn multiplication)
+
+func montyInv(M: static BigInt): static Limb =
+  ## Returns the Montgomery domain
+  ## magic number: -1/M[0] mod LimbSize
+  ## M[0] is the least significant limb of M
+  ## M must be odd.
+
+  # ######################################################################
+  # Implementation of modular multiplication inverse
+  # Assuming 2 positive integers a and m the modulo
+  #
+  # We are looking for z that solves `az ≡ 1 mod m`
+  #
+  # References:
+  #   - Knuth, The Art of Computer Programming, Vol2 p342
+  #   - Menezes, Handbook of Applied Cryptography (HAC), p610
+  #     http://cacr.uwaterloo.ca/hac/about/chap14.pdf
+
+  # Starting from the extended GCD formula (Bezout identity),
+  # `ax + by = gcd(x,y)` with input x,y and outputs a, b, gcd
+  # We assume a and m are coprimes, i.e. gcd is 1, otherwise no inverse
+  # `ax + my = 1` <=> `ax + my ≡ 1 mod m` <=> `ax ≡ 1 mod m`
+
+  # For Montgomery magic number, we are in a special case
+  # where a = M and m = 2^LimbSize.
+  # For a and m to be coprimes, a must be odd.
+  #
+  # M being a power of 2 greatly simplifies computation:
+  #  - https://crypto.stackexchange.com/questions/47493/how-to-determine-the-multiplicative-inverse-modulo-64-or-other-power-of-two
+  #  - http://groups.google.com/groups?selm=1994Apr6.093116.27805%40mnemosyne.cs.du.edu
+  #  - https://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two
+  #  - https://eprint.iacr.org/2017/411
+
+  # We have the following relation
+  # ax ≡ 1 (mod 2^k) <=> ax(2 - ax) ≡ 1 (mod 2^(2k))
+  #
+  # To get -1/M[0] mod LimbSize
+  #   <=> -1/M0 mod LS
+  #   <=> M0 x ≡ -1 (mod LS)
+  # we can either negate the resulting x of `ax(2 - ax) ≡ 1 (mod 2^(2k))`
+  # or do ax(2 + ax) ≡ 1 (mod 2^(2k))
+
+  const
+    M0 = M.limbs[0]
+    log2Limb = fastLog2(Limb.sizeof * 8)
+
+  result = M                 # Start from an inverse of M0 modulo 2, M0 is odd and it's own inverse
+  for _ in 1 ..< log2Limb:
+    result *= 2 + M * result # x' = x(2 + ax) (`+` to avoid negating at the end)