From 95e23339b21ab8beac5b5c36bfd9bb197f7d80f9 Mon Sep 17 00:00:00 2001 From: Mamy Ratsimbazafy Date: Fri, 29 Jan 2021 20:42:36 +0100 Subject: [PATCH] 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] --- constantine/arithmetic.nim | 3 +- constantine/arithmetic/README.md | 4 +- constantine/arithmetic/bigints.nim | 167 ++------------- constantine/arithmetic/bigints_montgomery.nim | 167 +++++++++++++++ constantine/arithmetic/finite_fields.nim | 2 +- constantine/arithmetic/limbs.nim | 38 ++++ constantine/arithmetic/limbs_modular.nim | 2 +- constantine/arithmetic/limbs_montgomery.nim | 8 +- constantine/config/README.md | 16 ++ constantine/config/precompute.nim | 67 +++--- .../elliptic/ec_endomorphism_accel.nim | 3 +- constantine/io/README.md | 78 +++++++ constantine/io/io_bigints.nim | 199 +++++++++++++++++- constantine/io/io_fields.nim | 47 +++++ constantine/primitives/bithacks.nim | 12 +- helpers/prng_unsafe.nim | 2 +- tests/t_bigints_mod_vs_gmp.nim | 2 +- tests/t_io_bigints.nim | 16 +- 18 files changed, 631 insertions(+), 202 deletions(-) create mode 100644 constantine/arithmetic/bigints_montgomery.nim diff --git a/constantine/arithmetic.nim b/constantine/arithmetic.nim index ed0284b..7e5339a 100644 --- 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, diff --git a/constantine/arithmetic/README.md b/constantine/arithmetic/README.md index b1c20f9..48befc5 100644 --- 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. diff --git a/constantine/arithmetic/bigints.nim b/constantine/arithmetic/bigints.nim index bf3c7f6..cb7dcbd 100644 --- 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 diff --git a/constantine/arithmetic/bigints_montgomery.nim b/constantine/arithmetic/bigints_montgomery.nim new file mode 100644 index 0000000..489d610 --- /dev/null +++ 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 diff --git a/constantine/arithmetic/finite_fields.nim b/constantine/arithmetic/finite_fields.nim index fdd5961..7760908 100644 --- 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 diff --git a/constantine/arithmetic/limbs.nim b/constantine/arithmetic/limbs.nim index 153975c..cb46557 100644 --- 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 diff --git a/constantine/arithmetic/limbs_modular.nim b/constantine/arithmetic/limbs_modular.nim index cffe429..82f6d80 100644 --- 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. diff --git a/constantine/arithmetic/limbs_montgomery.nim b/constantine/arithmetic/limbs_montgomery.nim index 31ba98d..b8afe32 100644 --- 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 diff --git a/constantine/config/README.md b/constantine/config/README.md index 75266ba..b045dd0 100644 --- 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 diff --git a/constantine/config/precompute.nim b/constantine/config/precompute.nim index 5844555..a5cf9f5 100644 --- 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) diff --git a/constantine/elliptic/ec_endomorphism_accel.nim b/constantine/elliptic/ec_endomorphism_accel.nim index 1877827..37e245f 100644 --- 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 diff --git a/constantine/io/README.md b/constantine/io/README.md index d13da22..4d7d059 100644 --- 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 = "" +func toDecimal*(sk: SecretKey): string = "" +func `$`*(sk: SecretKey): string = "" +``` + ## 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 diff --git a/constantine/io/io_bigints.nim b/constantine/io/io_bigints.nim index 11d8cac..9a6cf6a 100644 --- 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." & + " ") + +# 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) diff --git a/constantine/io/io_fields.nim b/constantine/io/io_fields.nim index 3e2cbbd..5144d2a 100644 --- 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) diff --git a/constantine/primitives/bithacks.nim b/constantine/primitives/bithacks.nim index 160ff2c..cf22287 100644 --- 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))) diff --git a/helpers/prng_unsafe.nim b/helpers/prng_unsafe.nim index ebaaeb1..b156d79 100644 --- 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/[ diff --git a/tests/t_bigints_mod_vs_gmp.nim b/tests/t_bigints_mod_vs_gmp.nim index 241df7c..4ea6ea7 100644 --- 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: diff --git a/tests/t_io_bigints.nim b/tests/t_io_bigints.nim index 3f66984..d31bbb0 100644 --- 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()