Decimal conversion (#139)

* Add constant-time fromDecimal conversion. Add warnings on intended purposes of hex/decimals * introduce setuint + cosmetic fixes Wordbitsize -> Wordbitwidth in comments * Add decimal conversion (non-constant-time) * fix comments [skip ci]
2025-02-13 11:16:20 +00:00 · 2021-01-29 20:42:36 +01:00 · 2021-01-29 20:42:36 +01:00 · 95e23339b2
commit 95e23339b2
parent 47daefde1f
18 changed files with 631 additions and 202 deletions
--- a/constantine/arithmetic.nim
+++ b/constantine/arithmetic.nim
@ -7,7 +7,7 @@
 # at your option. This file may not be copied, modified, or distributed except according to those terms.

 import
-  arithmetic/bigints,
+  arithmetic/[bigints, bigints_montgomery],
  arithmetic/[
    finite_fields,
    finite_fields_inversion,
@ -17,6 +17,7 @@ import

 export
  bigints,
+  bigints_montgomery,
  finite_fields,
  finite_fields_inversion,
  finite_fields_square_root,
--- a/constantine/arithmetic/README.md
+++ b/constantine/arithmetic/README.md
@ -15,7 +15,9 @@ for example secp256k1 and BN254.
 It also enables compiler unrolling, inlining and register optimization,
 where code size is not an issue for example for multi-precision addition.

-## References
+## Algorithms
+
+### Finite field multiplication

 - Analyzing and Comparing Montgomery Multiplication Algorithms
  Cetin Kaya Koc and Tolga Acar and Burton S. Kaliski Jr.
--- a/constantine/arithmetic/bigints.nim
+++ b/constantine/arithmetic/bigints.nim
@ -10,8 +10,7 @@ import
  ../config/[common, type_bigint],
  ../primitives,
  ./limbs,
-  ./limbs_modular,
-  ./limbs_montgomery
+  ./limbs_modular

 export BigInt

@ -74,6 +73,10 @@ func setOne*(a: var BigInt) =
  ## Set a BigInt to 1
  a.limbs.setOne()

+func setUint*(a: var BigInt, n: SomeUnsignedInt) =
+  ## Set a BigInt to a machine-sized integer ``n``
+  a.limbs.setUint(n)
+
 func czero*(a: var BigInt, ctl: SecretBool) =
  ## Set ``a`` to 0 if ``ctl`` is true
  a.limbs.czero(ctl)
@ -254,10 +257,6 @@ func double*(r: var BigInt, a: BigInt): SecretBool =
  ## Returns the carry
  (SecretBool) sum(r.limbs, a.limbs, a.limbs)

-func div2*(a: var BigInt) =
-  ## In-place divide ``a`` by 2
-  a.limbs.shiftRight(1)
-
 func cneg*(a: var BigInt, ctl: CTBool) =
  ## Conditional negation.
  ## Negate if ``ctl`` is true
@ -407,6 +406,18 @@ func `*`*(b: static int, a: BigInt): BigInt {.noinit, inline.} =
  result = a
  result *= b

+# Division by constants
+# ------------------------------------------------------------
+
+func div2*(a: var BigInt) =
+  ## In-place divide ``a`` by 2
+  a.limbs.shiftRight(1)
+
+func div10*(a: var BigInt): SecretWord =
+  ## In-place divide ``a`` by 10
+  ## and return the remainder
+  a.limbs.div10()
+
 # ############################################################
 #
 #                   Modular BigInt
@ -455,149 +466,5 @@ func invmod*[bits](r: var BigInt[bits], a, M, mp1div2: BigInt[bits]) =
  one.setOne()
  r.steinsGCD(a, one, M, mp1div2)

-# ############################################################
-#
-#                 Montgomery Arithmetic
-#
-# ############################################################
-
-func montyResidue*(mres: var BigInt, a, N, r2modM: BigInt, m0ninv: static BaseType, canUseNoCarryMontyMul: static bool) =
-  ## Convert a BigInt from its natural representation
-  ## to the Montgomery n-residue form
-  ##
-  ## `mres` is overwritten. It's bitlength must be properly set before calling this procedure.
-  ##
-  ## Caller must take care of properly switching between
-  ## the natural and montgomery domain.
-  ## Nesting Montgomery form is possible by applying this function twice.
-  ##
-  ## The Montgomery Magic Constants:
-  ## - `m0ninv` is µ = -1/N (mod M)
-  ## - `r2modM` is R² (mod M)
-  ## with W = M.len
-  ## and R = (2^WordBitSize)^W
-  montyResidue(mres.limbs, a.limbs, N.limbs, r2modM.limbs, m0ninv, canUseNoCarryMontyMul)
-
-func redc*[mBits](r: var BigInt[mBits], a, M: BigInt[mBits], m0ninv: static BaseType, canUseNoCarryMontyMul: static bool) =
-  ## Convert a BigInt from its Montgomery n-residue form
-  ## to the natural representation
-  ##
-  ## `mres` is modified in-place
-  ##
-  ## Caller must take care of properly switching between
-  ## the natural and montgomery domain.
-  let one = block:
-    var one {.noInit.}: BigInt[mBits]
-    one.setOne()
-    one
-  redc(r.limbs, a.limbs, one.limbs, M.limbs, m0ninv, canUseNoCarryMontyMul)
-
-func montyMul*(r: var BigInt, a, b, M: BigInt, negInvModWord: static BaseType, canUseNoCarryMontyMul: static bool) =
-  ## Compute r <- a*b (mod M) in the Montgomery domain
-  ##
-  ## This resets r to zero before processing. Use {.noInit.}
-  ## to avoid duplicating with Nim zero-init policy
-  montyMul(r.limbs, a.limbs, b.limbs, M.limbs, negInvModWord, canUseNoCarryMontyMul)
-
-func montySquare*(r: var BigInt, a, M: BigInt, negInvModWord: static BaseType, canUseNoCarryMontyMul: static bool) =
-  ## Compute r <- a^2 (mod M) in the Montgomery domain
-  ##
-  ## This resets r to zero before processing. Use {.noInit.}
-  ## to avoid duplicating with Nim zero-init policy
-  montySquare(r.limbs, a.limbs, M.limbs, negInvModWord, canUseNoCarryMontyMul)
-
-func montyPow*[mBits: static int](
-       a: var BigInt[mBits], exponent: openarray[byte],
-       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
-       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
-      ) =
-  ## Compute a <- a^exponent (mod M)
-  ## ``a`` in the Montgomery domain
-  ## ``exponent`` is a BigInt in canonical big-endian representation
-  ##
-  ## Warning ⚠️ :
-  ## This is an optimization for public exponent
-  ## Otherwise bits of the exponent can be retrieved with:
-  ## - memory access analysis
-  ## - power analysis
-  ## - timing analysis
-  ##
-  ## This uses fixed window optimization
-  ## A window size in the range [1, 5] must be chosen
-
-  const scratchLen = if windowSize == 1: 2
-                     else: (1 shl windowSize) + 1
-  var scratchSpace {.noInit.}: array[scratchLen, Limbs[mBits.wordsRequired]]
-  montyPow(a.limbs, exponent, M.limbs, one.limbs, negInvModWord, scratchSpace, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
-
-func montyPowUnsafeExponent*[mBits: static int](
-       a: var BigInt[mBits], exponent: openarray[byte],
-       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
-       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
-      ) =
-  ## Compute a <- a^exponent (mod M)
-  ## ``a`` in the Montgomery domain
-  ## ``exponent`` is a BigInt in canonical big-endian representation
-  ##
-  ## Warning ⚠️ :
-  ## This is an optimization for public exponent
-  ## Otherwise bits of the exponent can be retrieved with:
-  ## - memory access analysis
-  ## - power analysis
-  ## - timing analysis
-  ##
-  ## This uses fixed window optimization
-  ## A window size in the range [1, 5] must be chosen
-
-  const scratchLen = if windowSize == 1: 2
-                     else: (1 shl windowSize) + 1
-  var scratchSpace {.noInit.}: array[scratchLen, Limbs[mBits.wordsRequired]]
-  montyPowUnsafeExponent(a.limbs, exponent, M.limbs, one.limbs, negInvModWord, scratchSpace, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
-
-from ../io/io_bigints import exportRawUint
-# Workaround recursive dependencies
-
-func montyPow*[mBits, eBits: static int](
-       a: var BigInt[mBits], exponent: BigInt[eBits],
-       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
-       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
-      ) =
-  ## Compute a <- a^exponent (mod M)
-  ## ``a`` in the Montgomery domain
-  ## ``exponent`` is any BigInt, in the canonical domain
-  ##
-  ## This uses fixed window optimization
-  ## A window size in the range [1, 5] must be chosen
-  ##
-  ## This is constant-time: the window optimization does
-  ## not reveal the exponent bits or hamming weight
-  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
-  expBE.exportRawUint(exponent, bigEndian)
-
-  montyPow(a, expBE, M, one, negInvModWord, windowSize, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
-
-func montyPowUnsafeExponent*[mBits, eBits: static int](
-       a: var BigInt[mBits], exponent: BigInt[eBits],
-       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
-       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
-      ) =
-  ## Compute a <- a^exponent (mod M)
-  ## ``a`` in the Montgomery domain
-  ## ``exponent`` is any BigInt, in the canonical domain
-  ##
-  ## Warning ⚠️ :
-  ## This is an optimization for public exponent
-  ## Otherwise bits of the exponent can be retrieved with:
-  ## - memory access analysis
-  ## - power analysis
-  ## - timing analysis
-  ##
-  ## This uses fixed window optimization
-  ## A window size in the range [1, 5] must be chosen
-  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
-  expBE.exportRawUint(exponent, bigEndian)
-
-  montyPowUnsafeExponent(a, expBE, M, one, negInvModWord, windowSize, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
-
 {.pop.} # inline
 {.pop.} # raises no exceptions
--- a/constantine/arithmetic/bigints_montgomery.nim
+++ b/constantine/arithmetic/bigints_montgomery.nim
@ -0,0 +1,167 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+import
+  ../config/[common, type_bigint],
+  ../primitives,
+  ../io/io_bigints,
+  ./limbs,
+  ./limbs_modular,
+  ./limbs_montgomery,
+  ./bigints
+
+# No exceptions allowed
+{.push raises: [].}
+{.push inline.}
+
+# ############################################################
+#
+#                 Montgomery Arithmetic
+#
+# ############################################################
+
+func montyResidue*(mres: var BigInt, a, N, r2modM: BigInt, m0ninv: static BaseType, canUseNoCarryMontyMul: static bool) =
+  ## Convert a BigInt from its natural representation
+  ## to the Montgomery n-residue form
+  ##
+  ## `mres` is overwritten. It's bitlength must be properly set before calling this procedure.
+  ##
+  ## Caller must take care of properly switching between
+  ## the natural and montgomery domain.
+  ## Nesting Montgomery form is possible by applying this function twice.
+  ##
+  ## The Montgomery Magic Constants:
+  ## - `m0ninv` is µ = -1/N (mod M)
+  ## - `r2modM` is R² (mod M)
+  ## with W = M.len
+  ## and R = (2^WordBitWidth)^W
+  montyResidue(mres.limbs, a.limbs, N.limbs, r2modM.limbs, m0ninv, canUseNoCarryMontyMul)
+
+func redc*[mBits](r: var BigInt[mBits], a, M: BigInt[mBits], m0ninv: static BaseType, canUseNoCarryMontyMul: static bool) =
+  ## Convert a BigInt from its Montgomery n-residue form
+  ## to the natural representation
+  ##
+  ## `mres` is modified in-place
+  ##
+  ## Caller must take care of properly switching between
+  ## the natural and montgomery domain.
+  let one = block:
+    var one {.noInit.}: BigInt[mBits]
+    one.setOne()
+    one
+  redc(r.limbs, a.limbs, one.limbs, M.limbs, m0ninv, canUseNoCarryMontyMul)
+
+func montyMul*(r: var BigInt, a, b, M: BigInt, negInvModWord: static BaseType, canUseNoCarryMontyMul: static bool) =
+  ## Compute r <- a*b (mod M) in the Montgomery domain
+  ##
+  ## This resets r to zero before processing. Use {.noInit.}
+  ## to avoid duplicating with Nim zero-init policy
+  montyMul(r.limbs, a.limbs, b.limbs, M.limbs, negInvModWord, canUseNoCarryMontyMul)
+
+func montySquare*(r: var BigInt, a, M: BigInt, negInvModWord: static BaseType, canUseNoCarryMontyMul: static bool) =
+  ## Compute r <- a^2 (mod M) in the Montgomery domain
+  ##
+  ## This resets r to zero before processing. Use {.noInit.}
+  ## to avoid duplicating with Nim zero-init policy
+  montySquare(r.limbs, a.limbs, M.limbs, negInvModWord, canUseNoCarryMontyMul)
+
+func montyPow*[mBits: static int](
+       a: var BigInt[mBits], exponent: openarray[byte],
+       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
+       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
+      ) =
+  ## Compute a <- a^exponent (mod M)
+  ## ``a`` in the Montgomery domain
+  ## ``exponent`` is a BigInt in canonical big-endian representation
+  ##
+  ## Warning ⚠️ :
+  ## This is an optimization for public exponent
+  ## Otherwise bits of the exponent can be retrieved with:
+  ## - memory access analysis
+  ## - power analysis
+  ## - timing analysis
+  ##
+  ## This uses fixed window optimization
+  ## A window size in the range [1, 5] must be chosen
+
+  const scratchLen = if windowSize == 1: 2
+                     else: (1 shl windowSize) + 1
+  var scratchSpace {.noInit.}: array[scratchLen, Limbs[mBits.wordsRequired]]
+  montyPow(a.limbs, exponent, M.limbs, one.limbs, negInvModWord, scratchSpace, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
+
+func montyPowUnsafeExponent*[mBits: static int](
+       a: var BigInt[mBits], exponent: openarray[byte],
+       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
+       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
+      ) =
+  ## Compute a <- a^exponent (mod M)
+  ## ``a`` in the Montgomery domain
+  ## ``exponent`` is a BigInt in canonical big-endian representation
+  ##
+  ## Warning ⚠️ :
+  ## This is an optimization for public exponent
+  ## Otherwise bits of the exponent can be retrieved with:
+  ## - memory access analysis
+  ## - power analysis
+  ## - timing analysis
+  ##
+  ## This uses fixed window optimization
+  ## A window size in the range [1, 5] must be chosen
+
+  const scratchLen = if windowSize == 1: 2
+                     else: (1 shl windowSize) + 1
+  var scratchSpace {.noInit.}: array[scratchLen, Limbs[mBits.wordsRequired]]
+  montyPowUnsafeExponent(a.limbs, exponent, M.limbs, one.limbs, negInvModWord, scratchSpace, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
+
+from ../io/io_bigints import exportRawUint
+# Workaround recursive dependencies
+
+func montyPow*[mBits, eBits: static int](
+       a: var BigInt[mBits], exponent: BigInt[eBits],
+       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
+       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
+      ) =
+  ## Compute a <- a^exponent (mod M)
+  ## ``a`` in the Montgomery domain
+  ## ``exponent`` is any BigInt, in the canonical domain
+  ##
+  ## This uses fixed window optimization
+  ## A window size in the range [1, 5] must be chosen
+  ##
+  ## This is constant-time: the window optimization does
+  ## not reveal the exponent bits or hamming weight
+  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
+  expBE.exportRawUint(exponent, bigEndian)
+
+  montyPow(a, expBE, M, one, negInvModWord, windowSize, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
+
+func montyPowUnsafeExponent*[mBits, eBits: static int](
+       a: var BigInt[mBits], exponent: BigInt[eBits],
+       M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int,
+       canUseNoCarryMontyMul, canUseNoCarryMontySquare: static bool
+      ) =
+  ## Compute a <- a^exponent (mod M)
+  ## ``a`` in the Montgomery domain
+  ## ``exponent`` is any BigInt, in the canonical domain
+  ##
+  ## Warning ⚠️ :
+  ## This is an optimization for public exponent
+  ## Otherwise bits of the exponent can be retrieved with:
+  ## - memory access analysis
+  ## - power analysis
+  ## - timing analysis
+  ##
+  ## This uses fixed window optimization
+  ## A window size in the range [1, 5] must be chosen
+  var expBE {.noInit.}: array[(ebits + 7) div 8, byte]
+  expBE.exportRawUint(exponent, bigEndian)
+
+  montyPowUnsafeExponent(a, expBE, M, one, negInvModWord, windowSize, canUseNoCarryMontyMul, canUseNoCarryMontySquare)
+
+{.pop.} # inline
+{.pop.} # raises no exceptions
--- a/constantine/arithmetic/finite_fields.nim
+++ b/constantine/arithmetic/finite_fields.nim
@ -29,7 +29,7 @@
 import
  ../primitives,
  ../config/[common, type_ff, curves],
-  ./bigints, ./limbs_montgomery
+  ./bigints, ./bigints_montgomery

 when UseASM_X86_64:
  import ./assembly/limbs_asm_modular_x86
--- a/constantine/arithmetic/limbs.nim
+++ b/constantine/arithmetic/limbs.nim
@ -71,6 +71,22 @@ func setOne*(a: var Limbs) =
  when a.len > 1:
    zeroMem(a[1].addr, (a.len - 1) * sizeof(SecretWord))

+func setUint*(a: var Limbs, n: SomeUnsignedInt) =
+  ## set ``a`` to an unsigned integer ``n``
+  when sizeof(SecretWord) >= sizeof(n):
+    a[0] = SecretWord(n)
+    when a.len > 1:
+      zeroMem(a[1].addr, (a.len - 1) * sizeof(SecretWord))
+  else:
+    static: doAssert a.len >= 2,
+      "Overflow, trying to store a " & $(sizeof(n)*8) & " integer " &
+      "in ", a.len, " limb of size ", sizeof(SecretWord), "."
+
+    a[0] = SecretWord(n) # Truncate the upper part
+    a[1] = SecretWord(n shr log2(sizeof(SecretWord)))
+    when a.len > 2:
+      zeroMem(a[2].addr, (a.len - 2) * sizeof(SecretWord))
+
 func czero*(a: var Limbs, ctl: SecretBool) =
  ## Set ``a`` to 0 if ``ctl`` is true
  # Only used for FF neg in pure Nim fallback
@ -401,4 +417,26 @@ func prod_high_words*[rLen, aLen, bLen](

  r = z

+# Division
+# ------------------------------------------------------------
+
+func div10*(a: var Limbs): SecretWord =
+  ## Divide `a` by 10 in-place and return the remainder
+  ## TODO constant-time
+  result = Zero
+
+  let clz = WordBitWidth - 1 - log2(10)
+  let norm10 = SecretWord(10) shl clz
+
+  for i in countdown(a.len-1, 0):
+    # dividend = 2^64 * remainder + a[i]
+    var hi = result
+    var lo = a[i]
+    # Normalize
+    hi = (hi shl clz) or (lo shr (WordBitWidth - clz))
+    lo = lo shl clz
+    unsafeDiv2n1n(a[i], result, hi, lo, norm10)
+    # Undo normalization
+    result = result shr clz
+
 {.pop.} # raises no exceptions
--- a/constantine/arithmetic/limbs_modular.nim
+++ b/constantine/arithmetic/limbs_modular.nim
@ -340,7 +340,7 @@ func shlAddMod(a: LimbsViewMut, aLen: int,
  ## Fused modular left-shift + add
  ## Shift input `a` by a word and add `c` modulo `M`
  ##
-  ## With a word W = 2^WordBitSize and a modulus M
+  ## With a word W = 2^WordBitWidth and a modulus M
  ## Does a <- a * W + c (mod M)
  ##
  ## The modulus `M` most-significant bit at `mBits` MUST be set.
--- a/constantine/arithmetic/limbs_montgomery.nim
+++ b/constantine/arithmetic/limbs_montgomery.nim
@ -344,7 +344,7 @@ func redc*(r: var Limbs, a, one, M: Limbs,
  ## to the regular natural representation (mod N)
  ##
  ## with W = M.len
-  ## and R = (2^WordBitSize)^W
+  ## and R = (2^WordBitWidth)^W
  ##
  ## Does "a * R^-1 (mod M)"
  ##
@ -368,13 +368,13 @@ func montyResidue*(r: var Limbs, a, M, r2modM: Limbs,
  ## Montgomery-Multiplication - based
  ##
  ## with W = M.len
-  ## and R = (2^WordBitSize)^W
+  ## and R = (2^WordBitWidth)^W
  ##
  ## Does "a * R (mod M)"
  ##
  ## `a`: The source BigInt in the natural representation. `a` in [0, N) range
  ## `M`: The field modulus. M must be odd.
-  ## `r2modM`: 2^WordBitSize mod `M`. Can be precomputed with `r2mod` function
+  ## `r2modM`: 2^WordBitWidth mod `M`. Can be precomputed with `r2mod` function
  ##
  ## Important: `r` is overwritten
  ## The result `r` buffer size MUST be at least the size of `M` buffer
@ -515,7 +515,7 @@ func montyPow*(
  ##   Use ``exportRawUint`` for conversion
  ## - ``M`` is the modulus
  ## - ``one`` is 1 (mod M) in montgomery representation
-  ## - ``m0ninv`` is the montgomery magic constant "-1/M[0] mod 2^WordBitSize"
+  ## - ``m0ninv`` is the montgomery magic constant "-1/M[0] mod 2^WordBitWidth"
  ## - ``scratchspace`` with k the window bitsize of size up to 5
  ##   This is a buffer that can hold between 2^k + 1 big-ints
  ##   A window of of 1-bit (no window optimization) requires only 2 big-ints
--- a/constantine/config/README.md
+++ b/constantine/config/README.md
@ -3,3 +3,19 @@
 - Low-level logical and physical word definitions
 - Elliptic curve declarations
 - Cipher suites
+
+## Algorithms
+
+
+### Modular inverses mod 2ⁿ
+
+We use "Dumas iterations" to precompute Montgomery magic number `-1/n[0] (mod 2^Wordbitwidth)`
+
+Explanation p11 "Dumas iterations" based on Newton-Raphson:
+- Cetin Kaya Koc (2017), https://eprint.iacr.org/2017/411
+- 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:
+- https://crypto.stackexchange.com/questions/47493/how-to-determine-the-multiplicative-inverse-modulo-64-or-other-power-of-two
+- https://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two
+- http://marc-b-reynolds.github.io/math/2017/09/18/ModInverse.html
--- a/constantine/config/precompute.nim
+++ b/constantine/config/precompute.nim
@ -226,8 +226,11 @@ func shiftRight*(a: var BigInt, k: int) =
 #
 # ############################################################

-func checkOddModulus(M: BigInt) =
-  doAssert bool(BaseType(M.limbs[0]) and 1), "Internal Error: the modulus must be odd to use the Montgomery representation."
+func checkOdd(a: BaseType) =
+  doAssert bool(a and 1), "Internal Error: the modulus must be odd to use the Montgomery representation."
+
+func checkOdd(M: BigInt) =
+  checkOdd(BaseType M.limbs[0])

 func checkValidModulus(M: BigInt) =
  const expectedMsb = M.bits-1 - WordBitWidth * (M.limbs.len - 1)
@ -251,18 +254,7 @@ func useNoCarryMontySquare*(M: BigInt): bool =
  # https://github.com/nim-lang/Nim/issues/9679
  BaseType(M.limbs[^1]) < high(BaseType) shr 2

-func negInvModWord*(M: BigInt): BaseType =
-  ## Returns the Montgomery domain magic constant for the input modulus:
-  ##
-  ##   µ ≡ -1/M[0] (mod SecretWord)
-  ##
-  ## M[0] is the least significant limb of M
-  ## M must be odd and greater than 2.
-  ##
-  ## Assuming 64-bit words:
-  ##
-  ## µ ≡ -1/M[0] (mod 2^64)
-
+func invModBitwidth[T: SomeUnsignedInt](a: T): T =
  # We use BaseType for return value because static distinct type
  # confuses Nim semchecks [UPSTREAM BUG]
  # We don't enforce compile-time evaluation here
@ -278,31 +270,34 @@ func negInvModWord*(M: BigInt): BaseType =
  # - https://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two
  # - http://marc-b-reynolds.github.io/math/2017/09/18/ModInverse.html

-  # For Montgomery magic number, we are in a special case
-  # where a = M and m = 2^WordBitWidth.
+  # We are in a special case
+  # where m = 2^WordBitWidth.
  # For a and m to be coprimes, a must be odd.
-
+  #
  # 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 negate the result x of `ax(2 - ax) ≡ 1 (mod 2^(2k))`
-  # or if k is odd: do ax(2 + ax) ≡ 1 (mod 2^(2k))
-  #
-  # To get the the modular inverse of 2^k' with arbitrary k'
-  # we can do modInv(a, 2^64) mod 2^63 as mentionned in Koc paper.
+  # which grows in O(log(log(a)))
+  checkOdd(a)

-  checkOddModulus(M)
+  let k = log2(T.sizeof() * 8)
+  result = a                 # Start from an inverse of M0 modulo 2, M0 is odd and it's own inverse
+  for _ in 0 ..< k:          # at each iteration we get the inverse mod(2^2k)
+    result *= 2 - a * result # x' = x(2 - ax)
+
+func negInvModWord*(M: BigInt): BaseType =
+  ## Returns the Montgomery domain magic constant for the input modulus:
+  ##
+  ##   µ ≡ -1/M[0] (mod SecretWord)
+  ##
+  ## M[0] is the least significant limb of M
+  ## M must be odd and greater than 2.
+  ##
+  ## Assuming 64-bit words:
+  ##
+  ## µ ≡ -1/M[0] (mod 2^64)
  checkValidModulus(M)

-  let
-    M0 = BaseType(M.limbs[0])
-    k = log2(WordBitWidth.uint32)
-
-  result = M0                 # Start from an inverse of M0 modulo 2, M0 is odd and it's own inverse
-  for _ in 0 ..< k:           # at each iteration we get the inverse mod(2^2k)
-    result *= 2 - M0 * result # x' = x(2 - ax)
-
+  result = invModBitwidth(BaseType M.limbs[0])
  # negate to obtain the negative inverse
  result = not(result) + 1

@ -329,7 +324,7 @@ func r_powmod(n: static int, M: BigInt): BigInt =
  #
  # Thus: C(wn+1) ≡ 2^(wn+1) C0 ≡ 2^(wn + 1) 2^(wn - 1) ≡ 2^(2wn) ≡ (2^wn)^2 ≡ R² (mod M)

-  checkOddModulus(M)
+  checkOdd(M)
  checkValidModulus(M)

  const
@ -380,7 +375,7 @@ func primePlus1div2*(P: BigInt): BigInt =
  ## For use in constant-time modular inversion
  ##
  ## Warning ⚠️: Result is in the canonical domain (not Montgomery)
-  checkOddModulus(P)
+  checkOdd(P)

  # (P+1)/2 = P/2 + 1 if P is odd,
  # this avoids overflowing if the prime uses all bits
@ -442,7 +437,7 @@ func primePlus1Div4_BE*[bits: static int](
  # - (bits + 7 - 1): dividing by 4 means 2 bits are unused
  #                   but we also add 1 to an odd number so using an extra bit
  # => TODO: reduce the output size (to potentially save a byte and corresponding multiplication/squarings)
-  checkOddModulus(P)
+  checkOdd(P)

  # First we do P+1/2 in a way that guarantees no overflow
  var tmp = primePlus1div2(P)
--- a/constantine/elliptic/ec_endomorphism_accel.nim
+++ b/constantine/elliptic/ec_endomorphism_accel.nim
@ -168,7 +168,8 @@ proc `[]=`(recoding: var Recoded,
           digitIdx: int, value: uint8) {.inline.}=
  ## 0 <= digitIdx < LengthInDigits
  ## returns digit ∈ {0, 1}
-  ## This is write-once
+  ## This is write-once, requires zero-init
+  ## This works in BigEndian mode
  const len = Recoded.LengthInDigits
  # assert digitIdx * BitSize < Recoded.LengthInDigits

--- a/constantine/io/README.md
+++ b/constantine/io/README.md
@ -1,5 +1,71 @@
 # I/O, serialization, encoding/decoding

+## Overview
+
+This folder provides serialization, encoding and decoding primitives
+from Constantine internal representation to/from a canonical representation.
+
+**Warning: ⚠️ This folder contains internal APIs**.
+As serialization protocols get added, hardened serialization primitives
+suitable for public use will be provided.
+
+### Internal API
+
+For the time being here is the description of the internal API.
+
+Constant-time APIs only leak the number of bits, of bytes or words of the
+big integer.
+
+The bytes API are constant-time and do not allocate:
+- BigInt or octet-string: fromRawUint, fromUint
+- Machine sized integers: fromUint
+
+If you decide to use the internal hex or decimal API, you SHOULD ensure that the data is well-formatted:
+- Only ['0'..'9'] for decimal
+- Only ['0'..'9'], ['A'..'F'], ['a'..'b'] for hexadecimal
+  An hexadecimal string may start with "0x".
+
+There is no input validation as those are used for configuration, prototyping, research and debugging purposes.
+If the data is too big for the specified BigInt size, the result is undefined.
+
+The internal API is may be constant-time (temporarily) and may allocate.
+
+The hexadecimal API allocates:
+- `toHex` is constant-time
+- `appendHex` is constant-time
+- `fromHex` is **NOT constant-time** (yet), it is intended for debugging or
+  (compile-time) configuration. It does not allocate.
+  In particular it scans spaces and underscores and checks if the string
+  starts with '0x'.
+
+The decimal API allocates:
+- `toDecimal` is **NOT constant-time** (yet) and allocates
+- `fromDecimal` is constant-time and does not allocate.
+
+## Avoiding secret mistakes
+
+Constantine deliberately doesn't define a `$` proc to make directly printing a compiler error or to make the datatype replaced by `...`.
+
+It is recommented that you wrap your own types as distinct type and you can go the extra mile of disabling associated proc:
+
+```
+type SecretKey = distinct BigInt[256]
+
+func toHex*(sk: SecretKey): string {.error: "Someone is about to make a big mistake.".}
+func toDecimal*(sk: SecretKey): string {.error: "Someone is about to make a big mistake.".}
+func `$`*(sk: SecretKey): string {.error: "Someone is about to make a big mistake.".}
+```
+
+Alternatively you can also overload with  dummy procedures:
+
+```
+type SecretKey = distinct BigInt[256]
+
+func toHex*(sk: SecretKey): string = "<secret>"
+func toDecimal*(sk: SecretKey): string = "<smoke screen>"
+func `$`*(sk: SecretKey): string = "<camouflage>"
+```
+
 ## References

 ### Normative references
@ -7,3 +73,15 @@
 - Standards for Efficient Cryptography Group (SECG),\
  "SEC 1: Elliptic Curve Cryptography", May 2009,\
  http://www.secg.org/sec1-v2.pdf
+
+### Algorithms
+
+#### Continued fractions
+
+Continued fractions are used to convert
+
+`size_in_bits <=> size_in_decimal`
+
+for constant-time buffer preallocation when converting to decimal.
+
+- https://en.wikipedia.org/wiki/Continued_fraction
--- a/constantine/io/io_bigints.nim
+++ b/constantine/io/io_bigints.nim
@ -11,7 +11,8 @@
 #   - Burning memory to ensure secrets are not left after dealloc.

 import
-  ../primitives/constant_time,
+  ../primitives,
+  ../arithmetic/bigints,
  ../config/[common, type_bigint],
  ./endians

@ -21,11 +22,18 @@ import
 #
 # ############################################################

+# No exceptions for the byte API
+{.push raises: [].}
+
 # Note: the parsing/serialization routines were initially developed
 #       with an internal representation that used 31 bits out of a uint32
 #       or 63-bits out of an uint64

-# TODO: tag/remove exceptions raised.
+# TODO: the in-place API should return a bool
+#       to indicate success.
+#       the out-of place API are for for configuration,
+#       prototyping, research and debugging purposes,
+#       and can use exceptions.

 func fromRawUintLE(
        dst: var BigInt,
@ -289,15 +297,19 @@ func exportRawUint*(
  else:
    exportRawUintBE(dst, src)

+{.pop.} # {.push raises: [].}
+
 # ############################################################
 #
 #         Conversion helpers
 #
 # ############################################################
+# TODO: constant-time

 func readHexChar(c: char): uint8 {.inline.}=
  ## Converts an hex char to an int
  ## CT: leaks position of invalid input if any.
+  ##     and leaks if 0..9 or a..f or A..F
  case c
  of '0'..'9': result = uint8 ord(c) - ord('0')
  of 'a'..'f': result = uint8 ord(c) - ord('a') + 10
@ -467,3 +479,186 @@ func toHex*(big: BigInt, order: static Endianness = bigEndian): string =
  ## CT:
  ##   - no leaks
  result.appendHex(big, order)
+
+# ############################################################
+#
+#                    Decimal conversion
+#
+# ############################################################
+#
+# We need to convert between the size in binary
+# and the size in decimal. Unlike for the hexadecimal case
+# this is not trivial. We find the following relation.
+#
+# binary_length = log₂(value)
+# decimal_length = log₁₀(value)
+#
+# Hence we have:
+# binary_length = (log₂(value) / log₁₀(value)) * decimal_length
+#
+# the log change of base formula allow us to express
+# log₁₀(value) = log₂(value) / log₂(10)
+#
+# Hence
+# binary_length = log₂(10) * decimal_length
+#
+# Now we need to approximate log₂(10), we can have a best approximation
+# using continued factions:
+# In Sagemath: "continued_fraction(log(10,2)).convergents()[0:10].list()"
+# [3,
+#  10/3,
+#  93/28,
+#  196/59,
+#  485/146,
+#  2136/643,
+#  13301/4004,
+#  28738/8651,
+#  42039/12655,
+#  70777/21306]
+#
+# According to http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/cfCALC.html
+# we have
+# 42039/12655	= [3;3,9,2,2,4,6,2,1]	  = 3.321928091663374	 error -3.223988631617658×10-9 (9.705×10-8%)
+# 70777/21306	= [3;3,9,2,2,4,6,2,1,1]	= 3.3219280953721957 error +4.848330625861763×10-10 (1.459×10-8%
+# as lower and upper bound.
+
+const log2_10_Num = 42039
+const log2_10_Denom = 12655
+
+# No exceptions for the in-place API
+{.push raises: [].}
+
+func hasEnoughBitsForDecimal(bits: uint, decimalLength: uint): bool =
+  ## Check if the decimalLength would fit in a big int of size bits.
+  ## This assumes that bits and decimal length are **public.**
+  ##
+  ## The check uses continued fraction approximation
+  ## In Sagemath: "continued_fraction(log(10,2)).convergents()[0:10].list()"
+  ## which gives 70777/21306
+  ##
+  ## The check might be too lenient by 1 bit.
+  if bits >= high(uint) div log2_10_Num:
+    # The next multiplication would overflow
+    return false
+  # Compute the expected length
+  let maxExpectedBitlength = ((decimalLength * log2_10_Num) div log2_10_Denom)
+
+  # A big int "400....." might take 381 bits and "500....." might take 382
+  let lenientBitlength = maxExpectedBitlength - 1
+
+  result = bits >= lenientBitlength
+
+func fromDecimal*[aBits: static int](a: var BigInt[aBits], s: string): SecretBool =
+  ## Convert a decimal string. The input must be packed
+  ## with no spaces or underscores.
+  ## This assumes that bits and decimal length are **public.**
+  ##
+  ## This function does approximate validation that the BigInt
+  ## can hold the input string.
+  ##
+  ## It is intended for configuration, prototyping, research and debugging purposes.
+  ## You MUST NOT use it for production.
+  ##
+  ## Return true if conversion is successful
+  ##
+  ## Return false if an error occured:
+  ## - There is not enough space in the BigInt
+  ## - An invalid character was found
+
+  if not aBits.hasEnoughBitsForDecimal(s.len.uint):
+    return CtFalse
+
+  a.setZero()
+  result = CtTrue
+
+  for i in 0 ..< s.len:
+    let c = SecretWord(ord(s[i]))
+    result = result and (SecretWord(ord('0')) <= c and c <= SecretWord(ord('9')))
+
+    a += c - SecretWord(ord('0'))
+    if i != s.len - 1:
+      a *= 10
+
+func fromDecimal*(T: type BigInt, s: string): T {.raises: [ValueError].}=
+  ## Convert a decimal string. The input must be packed
+  ## with no spaces or underscores.
+  ## This assumes that bits and decimal length are **public.**
+  ##
+  ## This function does approximate validation that the BigInt
+  ## can hold the input string.
+  ##
+  ## It is intended for configuration, prototyping, research and debugging purposes.
+  ## You MUST NOT use it for production.
+  ##
+  ## This function may raise an exception if input is incorrect
+  ## - There is not enough space in the BigInt
+  ## - An invalid character was found
+  let status = result.fromDecimal(s)
+  if not status.bool:
+    raise newException(ValueError,
+      "BigInt could not be parsed from decimal string." &
+      " <Potentially secret input withheld>")
+
+# Conversion to decimal
+# ----------------------------------------------------------------
+#
+# The first problem to solve is precomputing the final size
+# We also use continued fractions to approximate log₁₀(2)
+#
+# decimal_length = log₁₀(2) * binary_length
+#
+# sage: [(frac, numerical_approx(frac - log(2,10))) for frac in continued_fraction(log(2,10)).convergents()[0:10]]
+# [(0, -0.301029995663981),
+#  (1/3, 0.0323033376693522),
+#  (3/10, -0.00102999566398115),
+#  (28/93, 0.0000452731532231687),
+#  (59/196, -9.58750071583525e-6),
+#  (146/485, 9.32171070389121e-7),
+#  (643/2136, -3.31171646772432e-8),
+#  (4004/13301, 2.08054934391910e-9),
+#  (8651/28738, -5.35579747218407e-10),
+#  (12655/42039, 2.92154800352051e-10)]
+
+const log10_2_Num = 12655
+const log10_2_Denom = 42039
+
+# Then the naive way to serialize is to repeatedly do
+# const intToCharMap = "0123456789"
+# rest = number
+# while rest != 0:
+#   digitToPrint = rest mod 10
+#   result.add intToCharMap[digitToPrint]
+#   rest /= 10
+#
+# For constant-time we:
+# 1. can't compare with 0 as a stopping condition
+# 2. repeatedly add to a buffer
+# 3. can't use naive indexing (cache timing attacks though very unlikely given the small size)
+# 4. need (fast) constant-time division
+#
+# 1 and 2 is solved by precomputing the length and make the number of add be fixed.
+# 3 is easily solved by doing "digitToPrint + ord('0')" instead
+#
+# For 4 for now we use non-constant-time division (TODO)
+
+func decimalLength(bits: static int): int =
+  doAssert bits < (high(uint) div log10_2_Num),
+    "Constantine does not support that many bits to convert to a decimal string: " & $bits
+    # The next multiplication would overflow
+
+  result = 1 + ((bits * log10_2_Num) div log10_2_Denom)
+
+func toDecimal*(a: BigInt): string =
+  ## Convert to a decimal string.
+  ##
+  ## It is intended for configuration, prototyping, research and debugging purposes.
+  ## You MUST NOT use it for production.
+  ##
+  ## This function is NOT constant-time at the moment.
+  const len = decimalLength(BigInt.bits)
+  result = newString(len)
+
+  var a = a
+  for i in countdown(len-1, 0):
+    let c = ord('0') + a.div10().int
+    result[i] = char(c)
--- a/constantine/io/io_fields.nim
+++ b/constantine/io/io_fields.nim
@ -82,9 +82,56 @@ func toHex*(f: FF, order: static Endianness = bigEndian): string =

 func fromHex*(dst: var FF, hexString: string) {.raises: [ValueError].}=
  ## Convert a hex string to a element of Fp or Fr
+  # TODO: review API, should return bool
  let raw {.noinit.} = fromHex(dst.mres.typeof, hexString)
  dst.fromBig(raw)

 func fromHex*(T: type FF, hexString: string): T {.noInit, raises: [ValueError].}=
  ## Convert a hex string to a element of Fp
  result.fromHex(hexString)
+
+func toDecimal*(f: FF): string =
+  ## Convert to a decimal string.
+  ##
+  ## It is intended for configuration, prototyping, research and debugging purposes.
+  ## You MUST NOT use it for production.
+  ##
+  ## This function is NOT constant-time at the moment.
+  # TODO constant-time
+  f.toBig().toDecimal()
+
+func fromDecimal*(dst: var FF, decimalString: string) {.raises: [ValueError].}=
+  ## Convert a decimal string. The input must be packed
+  ## with no spaces or underscores.
+  ## This assumes that bits and decimal length are **public.**
+  ##
+  ## This function does approximate validation that the BigInt
+  ## can hold the input string.
+  ##
+  ## It is intended for configuration, prototyping, research and debugging purposes.
+  ## You MUST NOT use it for production.
+  ##
+  ## Return true if conversion is successful
+  ##
+  ## Return false if an error occured:
+  ## - There is not enough space in the BigInt
+  ## - An invalid character was found
+  # TODO: review API, should return bool
+  let raw {.noinit.} = fromDecimal(dst.mres.typeof, decimalString)
+  dst.fromBig(raw)
+
+func fromDecimal*(T: type FF, hexString: string): T {.noInit, raises: [ValueError].}=
+  ## Convert a decimal string. The input must be packed
+  ## with no spaces or underscores.
+  ## This assumes that bits and decimal length are **public.**
+  ##
+  ## This function does approximate validation that the BigInt
+  ## can hold the input string.
+  ##
+  ## It is intended for configuration, prototyping, research and debugging purposes.
+  ## You MUST NOT use it for production.
+  ##
+  ## This function may raise an exception if input is incorrect
+  ## - There is not enough space in the BigInt
+  ## - An invalid character was found
+  result.fromDecimal(hexString)
--- a/constantine/primitives/bithacks.nim
+++ b/constantine/primitives/bithacks.nim
@ -41,7 +41,7 @@ proc clearBit*[T: SomeInteger](v: T, bit: T): T {.inline.} =
  ## Returns ``v``, with the bit at position ``bit`` set to 0
  v.clearMask(1.T shl bit)

-func log2*(x: uint32): uint32 =
+func log2Impl(x: uint32): uint32 =
  ## Find the log base 2 of a 32-bit or less integer.
  ## using De Bruijn multiplication
  ## Works at compile-time, guaranteed constant-time.
@ -57,7 +57,7 @@ func log2*(x: uint32): uint32 =
  v = v or v shr 16
  lookup[(v * 0x07C4ACDD'u32) shr 27]

-func log2*(x: uint64): uint64 {.inline, noSideEffect.} =
+func log2Impl(x: uint64): uint64 {.inline, noSideEffect.} =
  ## Find the log base 2 of a 32-bit or less integer.
  ## using De Bruijn multiplication
  ## Works at compile-time, guaranteed constant-time.
@ -75,3 +75,11 @@ func log2*(x: uint64): uint64 {.inline, noSideEffect.} =
  v = v or v shr 16
  v = v or v shr 32
  lookup[(v * 0x03F6EAF2CD271461'u64) shr 58]
+
+func log2*[T: SomeUnsignedInt](n: T): T =
+  ## Find the log base 2 of an integer
+  when sizeof(T) == sizeof(uint64):
+    T(log2Impl(uint64(n)))
+  else:
+    static: doAssert sizeof(T) <= sizeof(uint32)
+    T(log2Impl(uint32(n)))
--- a/helpers/prng_unsafe.nim
+++ b/helpers/prng_unsafe.nim
@ -7,7 +7,7 @@
 # at your option. This file may not be copied, modified, or distributed except according to those terms.

 import
-  ../constantine/arithmetic/bigints,
+  ../constantine/arithmetic,
  ../constantine/primitives,
  ../constantine/config/[common, curves, type_ff],
  ../constantine/elliptic/[
--- a/tests/t_bigints_mod_vs_gmp.nim
+++ b/tests/t_bigints_mod_vs_gmp.nim
@ -63,7 +63,7 @@ macro testRandomModSizes(numSizes: static int, aBits, mBits, body: untyped): unt
      if bool(bitSizeRNG.rand(high(int)) and 1):
        bitSizeRNG.sample(CryptoModSizes)
      else:
-        # range 62..1024 to highlight edge effects of the WordBitSize (63)
+        # range 62..1024 to highlight edge effects of the WordBitWidth (63)
        bitSizeRNG.rand(62 .. 1024)

    result.add quote do:
--- a/tests/t_io_bigints.nim
+++ b/tests/t_io_bigints.nim
@ -22,7 +22,7 @@ echo "test_io_bigints xoshiro512** seed: ", seed
 type T = BaseType

 proc main() =
-  suite "IO - BigInt" & " [" & $WordBitwidth & "-bit mode]":
+  suite "IO Hex - BigInt" & " [" & $WordBitwidth & "-bit mode]":
    test "Parsing raw integers":
      block: # Sanity check
        let x = 0'u64
@ -98,4 +98,18 @@ proc main() =

      check: n == h

+  suite "IO Decimal - BigInt" & " [" & $WordBitwidth & "-bit mode]":
+    test "Checks elliptic curve constants":
+      block: # BLS12-381 - https://github.com/ethereum/py_ecc/blob/master/py_ecc/fields/field_properties.py
+        const p = "4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129015664037894272559787"
+        const pHex = "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab"
+        let x = BigInt[381].fromHex(pHex)
+        let dec = x.toDecimal()
+
+        check: p == dec
+
+        let y = BigInt[381].fromDecimal(p)
+        let hex = y.toHex()
+
+        check: pHex == hex
 main()