Use Montgomery representation by default for Finite Field

- Fix montyMagic, modular inversion mode 2^2k was missing an iteration
- Fix test for buffer size in BigInt serialization
- Add UINT/Hex serialization for finite fields
- Montgomery conversion and redc

constantine.nimble

@@ -21,3 +21,4 @@ task test, "Run all tests":
   test "",                  "tests/test_bigints.nim"
   test "",                  "tests/test_bigints_multimod.nim"
   test "",                  "tests/test_bigints_vs_gmp.nim"
+  test "",                  "tests/test_finite_fields.nim"

constantine/config/curves.nim

@@ -18,57 +18,51 @@ import
 #
 # ############################################################
 
-func montyMagic(M: static BigInt): static Word {.inline.} =
+# Note: The declareCurve macro builds an exported montyMagic*(C: static Curve) on top of the one below
+func montyMagic(M: BigInt): BaseType =
   ## Returns the Montgomery domain magic constant for the input modulus:
   ##   -1/M[0] mod LimbSize
   ## M[0] is the least significant limb of M
   ## M must be odd and greater than 2.
+  # We use BaseType for return value because static distinct type
+  # confuses Nim semchecks [UPSTREAM BUG]
+  # We don't enforce compile-time evaluation here
+  # because static BigInt[bits] also causes semcheck troubles [UPSTREAM BUG]
 
-  # Test vectors: https://www.researchgate.net/publication/4107322_Montgomery_modular_multiplication_architecture_for_public_key_cryptosystems
+  # Modular inverse algorithm:
-  # on p354
+  # Explanation p11 "Dumas iterations" based on Newton-Raphson:
-  # Reference C impl: http://www.hackersdelight.org/hdcodetxt/mont64.c.txt
+  # - Cetin Kaya Koc (2017), https://eprint.iacr.org/2017/411
-  # Explanation: https://eprint.iacr.org/2017/1057.pdf
+  # - Jean-Guillaume Dumas (2012), https://arxiv.org/pdf/1209.6626v2.pdf
-
+  # - Colin Plumb (1994), http://groups.google.com/groups?selm=1994Apr6.093116.27805%40mnemosyne.cs.du.edu
-  # ######################################################################
+  # Other sources:
-  # Implementation of modular multiplicative inverse
+  # - https://crypto.stackexchange.com/questions/47493/how-to-determine-the-multiplicative-inverse-modulo-64-or-other-power-of-two
-  # Assuming 2 positive integers a and m the modulo
+  # - https://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two
-  #
+  # - http://marc-b-reynolds.github.io/math/2017/09/18/ModInverse.html
-  # 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` (2^LimbSize) 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/M0 mod LimbSize
   # we can either negate the resulting x of `ax(2 - ax) ≡ 1 (mod 2^(2k))`
   # or do ax(2 + ax) ≡ 1 (mod 2^(2k))
+  #
+  # To get the the modular inverse of 2^k' with arbitrary k' (like k=63 in our case)
+  # we can do modInv(a, 2^64) mod 2^63 as mentionned in Koc paper.
 
-  const
+  let
-    M0 = M.limbs[0]
+    M0 = BaseType(M.limbs[0])
-    k = log2(WordBitSize)
+    k = log2(WordPhysBitSize)
 
-  result = M0                # Start from an inverse of M0 modulo 2, M0 is odd and it's own inverse
+  result = M0                 # Start from an inverse of M0 modulo 2, M0 is odd and it's own inverse
-  for _ in 0 ..< k:
+  for _ in 0 ..< k:           # at each iteration we get the inverse mod(2^2k)
-    result *= 2 + M * result # x' = x(2 + ax) (`+` to avoid negating at the end)
+    result *= 2 + M0 * result # x' = x(2 + ax) (`+` to avoid negating at the end)
+
+  # Our actual word size is 2^63 not 2^64
+  result = result and BaseType(MaxWord)
 
 # ############################################################
 #

constantine/config/curves_parser.nim

@@ -139,7 +139,7 @@ macro declareCurves*(curves: untyped): untyped =
   let C = ident"C"
   let matchingBigInt = ident"matchingBigInt"
   result.add newProc(
-    name = matchingBigInt,
+    name = nnkPostFix.newTree(ident"*", matchingBigInt),
     params = [ident"untyped", newIdentDefs(C, nnkStaticTy.newTree(Curve))],
     body = nnkBracketExpr.newTree(bindSym"BigInt", nnkBracketExpr.newTree(cbs, C)),
     procType = nnkTemplateDef
@@ -158,16 +158,13 @@ macro declareCurves*(curves: untyped): untyped =
       nnkGenericParams.newTree(newIdentDefs(
         C, nnkStaticTy.newTree(Curve), newEmptyNode()
       )),
+      # TODO: where should I put the nnkCommentStmt?
       nnkObjectTy.newTree(
         newEmptyNode(),
         newEmptyNode(),
         nnkRecList.newTree(
-          nnkCommentStmt.newTree "All operations on a field are modulo P",
-          nnkCommentStmt.newTree "P being the prime modulus of the Curve C",
-          nnkCommentStmt.newTree "Internally, data is stored in Montgomery n-residue form",
-          nnkCommentStmt.newTree "with the magic constant chosen for convenient division (a power of 2 depending on P bitsize)"
           newIdentDefs(
-            nnkPostfix.newTree(ident"*", ident"value"),
+            nnkPostfix.newTree(ident"*", ident"nres"),
             newCall(matchingBigInt, C)
           )
         )
@@ -203,7 +200,7 @@ macro declareCurves*(curves: untyped): untyped =
 
   # proc MontyMagic(curve: static Curve): static Word
   result.add newProc(
-    name = nnkPostfix.newTree(ident"*", ident"MontyMagic"),
+    name = nnkPostfix.newTree(ident"*", ident"montyMagic"),
     params = [
       ident"auto",
       newIdentDefs(
@@ -213,10 +210,14 @@ macro declareCurves*(curves: untyped): untyped =
     ],
     body = newCall(
       ident"montyMagic",
-      newCall(ident"Mod", ident"curve")
+      # curve.Mod().nres
+      nnkDotExpr.newTree(
+        newCall(ident"Mod", ident"curve"),
+        ident"nres"
+      )
+
     ),
-    procType = nnkFuncDef,
+    procType = nnkFuncDef
-    pragmas = nnkPragma.newTree(ident"compileTime")
   )
 
   # echo result.toStrLit()

constantine/io/io_bigints.nim

@@ -181,13 +181,14 @@ func serializeRawUint*(
         dstEndianness: static Endianness) =
   ## Serialize a bigint into its canonical big-endian or little endian
   ## representation.
-  ## A destination buffer of size "BigInt.bits div 8" at minimum is needed.
+  ## A destination buffer of size "(BigInt.bits + 7) div 8" at minimum is needed,
+  ## i.e. bits -> byte conversion rounded up
   ##
   ## If the buffer is bigger, output will be zero-padded left for big-endian
   ## or zero-padded right for little-endian.
   ## I.e least significant bit is aligned to buffer boundary
 
-  assert dst.len >= static(BigInt.bits div 8), "BigInt -> Raw int conversion: destination buffer is too small"
+  assert dst.len >= (BigInt.bits + 7) div 8, "BigInt -> Raw int conversion: destination buffer is too small"
 
   when BigInt.bits == 0:
     zeroMem(dst, dst.len)
@@ -339,16 +340,7 @@ func toHex*(big: BigInt, order: static Endianness = bigEndian): string =
   ## Note. Leading zeros are not removed.
   ## Result is prefixed with 0x
   ##
-  ## This is a raw memory dump. Output will be padded with 0
+  ## Output will be padded with 0s to maintain constant-time.
-  ## if the big int does not use the full memory allocated for it.
-  ##
-  ## Regardless of the machine endianness the output will be big-endian hex.
-  ##
-  ## For example a BigInt representing 10 will be
-  ##   - 0x0a                for BigInt[8]
-  ##   - 0x000a              for BigInt[16]
-  ##   - 0x00000000_0000000a for BigInt[64]
-  ## (underscore added for docuentation readability only)
   ##
   ## CT:
   ##   - no leaks

constantine/io/io_fields.nim

@@ -8,7 +8,8 @@
 
 import
   ./io_bigints,
-  ../math/finite_fields
+  ../config/curves,
+  ../math/[bigints_checked, finite_fields]
 
 # ############################################################
 #
@@ -16,8 +17,34 @@ import
 #
 # ############################################################
 
-func fromUint*(dst: var Fp,
+func fromUint*(dst: var Fq,
                src: SomeUnsignedInt) =
   ## Parse a regular unsigned integer
   ## and store it into a BigInt of size `bits`
-  dst.value.fromRawUint(cast[array[sizeof(src), byte]](src), cpuEndian)
+  dst.nres.fromRawUint(cast[array[sizeof(src), byte]](src), cpuEndian)
+  dst.nres.unsafeMontgomeryResidue(Fq.C.Mod.nres)
+
+func serializeRawUint*(dst: var openarray[byte],
+                       src: Fq,
+                       dstEndianness: static Endianness) =
+  ## Serialize a finite field element to its canonical big-endian or little-endian
+  ## representation
+  ## With `bits` the number of bits of the field modulus
+  ## a buffer of size "(bits + 7) div 8" at minimum is needed
+  ## i.e. bits -> byte conversion rounded up
+  ##
+  ## If the buffer is bigger, output will be zero-padded left for big-endian
+  ## or zero-padded right for little-endian.
+  ## I.e least significant bit is aligned to buffer boundary
+  serializeRawUint(dst, src.toBig(), dstEndianness)
+
+func toHex*(f: Fq, order: static Endianness = bigEndian): string =
+  ## Stringify a finite field element to hex.
+  ## Note. Leading zeros are not removed.
+  ## Result is prefixed with 0x
+  ##
+  ## Output will be padded with 0s to maintain constant-time.
+  ##
+  ## CT:
+  ##   - no leaks
+  result = f.toBig().toHex()

constantine/math/bigints_checked.nim

@@ -110,3 +110,24 @@ func reduce*[aBits, mBits](r: var BigInt[mBits], a: BigInt[aBits], M: BigInt[mBi
   # but we don't want to inline it as it would increase codesize, better have Nim
   # pass a pointer+length to a fixed session of the BSS.
   reduce(r.view, a.view, M.view)
+
+func unsafeMontgomeryResidue*[mBits](nres: var BigInt[mBits], N: BigInt[mBits]) =
+  ## Convert a BigInt from its natural representation
+  ## to the Montgomery n-residue form
+  ##
+  ## `nres` is modified in-place
+  ##
+  ## Caller must take care of properly switching between
+  ## the natural and montgomery domain.
+  ## Nesting Montgomery form is possible by applying this function twice.
+  montgomeryResidue(nres.view, N.view)
+
+func unsafeRedc*[mBits](nres: var BigInt[mBits], N: BigInt[mBits], montyMagic: static BaseType) =
+  ## Convert a BigInt from its Montgomery n-residue form
+  ## to the natural representation
+  ##
+  ## `nres` is modified in-place
+  ##
+  ## Caller must take care of properly switching between
+  ## the natural and montgomery domain.
+  redc(nres.view, N.view, Word(montyMagic))

constantine/math/bigints_raw.nim

@@ -181,6 +181,12 @@ template checkValidModulus(m: BigIntViewConst) =
   debug:
     assert not m[^1].isZero.bool, "Internal Error: the modulus must use all declared bits"
 
+template checkOddModulus(m: BigIntViewConst) =
+  ## CHeck that the modulus is odd
+  ## and valid for use in the Montgomery n-residue representation
+  debug:
+    assert bool(BaseType(m[0]) and 1), "Internal Error: the modulus must be odd to use the Montgomery representation."
+
 debug:
   func `$`*(a: BigIntViewAny): string =
     let len = a.numLimbs()
@@ -375,3 +381,56 @@ func reduce*(r: BigIntViewMut, a: BigIntViewAny, M: BigIntViewConst) =
     # Now shift-left the copied words while adding the new word modulo M
     for i in countdown(aOffset, 0):
       r.shlAddMod(a[i], M)
+
+func montgomeryResidue*(a: BigIntViewMut, N: BigIntViewConst) =
+  ## Transform a bigint ``a`` from it's natural representation (mod N)
+  ## to a the Montgomery n-residue representation
+  ## i.e. Does "a * (2^LimbSize)^W (mod N), where W is the number
+  ## of words needed to represent n in base 2^LimbSize
+  ##
+  ## `a`: The source BigInt in the natural representation. `a` in [0, N) range
+  ## `N`: The field modulus. N must be odd.
+  ##
+  ## Important: `a` is overwritten
+  # Reference: https://eprint.iacr.org/2017/1057.pdf
+  checkValidModulus(N)
+  checkOddModulus(N)
+  checkMatchingBitlengths(a, N)
+
+  let nLen = N.numLimbs()
+  for i in countdown(nLen, 1):
+    a.shlAddMod(Zero, N)
+
+func redc*(a: BigIntViewMut, N: BigIntViewConst, montyMagic: Word) =
+  ## Transform a bigint ``a`` from it's Montgomery N-residue representation (mod N)
+  ## to the regular natural representation (mod N)
+  ##
+  ## i.e. Does "a * ((2^LimbSize)^W)^-1 (mod N), where W is the number
+  ## of words needed to represent n in base 2^LimbSize
+  ##
+  ## This is called a Montgomery Reduction
+  ## The Montgomery Magic Constant is µ = -1/N mod N
+  ## is used internally and can be precomputed with montyMagic(Curve)
+  # References:
+  #   - https://eprint.iacr.org/2017/1057.pdf (Montgomery)
+  #     page: Radix-r interleaved multiplication algorithm
+  #   - https://en.wikipedia.org/wiki/Montgomery_modular_multiplication#Montgomery_arithmetic_on_multiprecision_(variable-radix)_integers
+  #   - http://langevin.univ-tln.fr/cours/MLC/extra/montgomery.pdf
+  #     Montgomery original paper
+  #
+  checkValidModulus(N)
+  checkOddModulus(N)
+  checkMatchingBitlengths(a, N)
+
+  let nLen = N.numLimbs()
+  for i in 0 ..< nLen:
+    let z0 = Word(BaseType(a[0]) * BaseType(montyMagic)) and MaxWord
+
+    var carry = DoubleWord(0)
+    for j in 0 ..< nLen:
+      let z = DoubleWord(a[i]) + unsafeExtPrecMul(z0, N[i]) + carry
+      carry = z shr WordBitSize
+      if j != 0:
+        a[j] = Word(z) and MaxWord
+
+    a[^1] = Word(carry)

constantine/math/finite_fields.nim

@@ -42,16 +42,21 @@ debug:
 # No exceptions allowed
 {.push raises: [].}
 
-func toMonty*[C: static Curve](a: Fq[C]): Montgomery[C] =
+# ############################################################
-  ## Convert a big integer over Fq to its montgomery representation
+#
-  ## over Fq.
+#                        Conversion
-  ## i.e. Does "a * (2^LimbSize)^W (mod p), where W is the number
+#
-  ## of words needed to represent p in base 2^LimbSize
+# ############################################################
 
-  result = a
+func fromBig*(T: type Fq, src: BigInt): T =
-  for i in static(countdown(C.Mod.limbs.high, 1)):
+  ## Convert a BigInt to its Montgomery form
-    shiftAdd(result, 0)
+  result.nres = src
+  result.nres.unsafeMontgomeryResidue(Fq.C.Mod)
 
+func toBig*[C: static Curve](src: Fq[C]): auto =
+  ## Convert a finite-field element to a BigInt in natral representation
+  result = src.nres
+  result.unsafeRedC(C.Mod.nres, montyMagic(C))
 
 # ############################################################
 #
@@ -66,7 +71,7 @@ template add(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] =
   ##
   ## a and b MAY be the same buffer
   ## a and b MUST have the same announced bitlength (i.e. `bits` static parameters)
-  add(a.value, b.value, ctl)
+  add(a.nres, b.nres, ctl)
 
 template sub(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] =
   ## Constant-time big integer in-place optional substraction
@@ -75,7 +80,7 @@ template sub(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] =
   ##
   ## a and b MAY be the same buffer
   ## a and b MUST have the same announced bitlength (i.e. `bits` static parameters)
-  sub(a.value, b.value, ctl)
+  sub(a.nres, b.nres, ctl)
 
 # ############################################################
 #
@@ -83,9 +88,6 @@ template sub(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] =
 #
 # ############################################################
 
-# No exceptions allowed
-{.push raises: [].}
-
 func `+=`*(a: var Fq, b: Fq) =
   ## Addition over Fq
   var ctl = add(a, b, CtTrue)

tests/test_finite_fields.nim

@@ -11,6 +11,8 @@ import  unittest, random,
         ../constantine/io/io_fields,
         ../constantine/config/curves
 
+import ../constantine/io/io_bigints
+
 static: doAssert defined(testingCurves), "This modules requires the -d:testingCurves compile option"
 
 proc main() =
@@ -24,7 +26,15 @@ proc main() =
         z.fromUint(90'u32)
 
         x += y
-        check: bool(z == x)
+
+        var x_bytes: array[8, byte]
+        x_bytes.serializeRawUint(x, cpuEndian)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          90'u64 == cast[uint64](x_bytes)
 
       block:
         var x, y, z: Fq[Fake101]
@@ -44,7 +54,15 @@ proc main() =
         z.fromUint(1'u32)
 
         x += y
-        check: bool(z == x)
+
+        var x_bytes: array[8, byte]
+        x_bytes.serializeRawUint(x, cpuEndian)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          1'u64 == cast[uint64](x_bytes)
 
     test "Substraction mod 101":
       block:
@@ -55,7 +73,15 @@ proc main() =
         z.fromUint(70'u32)
 
         x -= y
-        check: bool(z == x)
+
+        var x_bytes: array[8, byte]
+        x_bytes.serializeRawUint(x, cpuEndian)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          70'u64 == cast[uint64](x_bytes)
 
       block:
         var x, y, z: Fq[Fake101]
@@ -65,7 +91,15 @@ proc main() =
         z.fromUint(0'u32)
 
         x -= y
-        check: bool(z == x)
+
+        var x_bytes: array[8, byte]
+        x_bytes.serializeRawUint(x, cpuEndian)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          0'u64 == cast[uint64](x_bytes)
 
       block:
         var x, y, z: Fq[Fake101]
@@ -75,6 +109,14 @@ proc main() =
         z.fromUint(100'u32)
 
         x -= y
-        check: bool(z == x)
+
+        var x_bytes: array[8, byte]
+        x_bytes.serializeRawUint(x, cpuEndian)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          100'u64 == cast[uint64](x_bytes)
 
 main()