From 68727e5c8db780fa1df0dede05b543ab0acf58e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mamy=20Andr=C3=A9-Ratsimbazafy?= Date: Sat, 22 Feb 2020 19:50:24 +0100 Subject: [PATCH] Add modular inversion + test vs GMP --- constantine/config/curves.nim | 14 +++++ constantine/io/io_bigints.nim | 22 +++++--- constantine/math/bigints_checked.nim | 26 +++++++++ constantine/math/finite_fields.nim | 24 ++++++++- constantine/math/precomputed.nim | 15 +++++- tests/test_finite_fields_powinv.nim | 14 +++++ tests/test_finite_fields_vs_gmp.nim | 81 ++++++++++++++++++++++++++-- 7 files changed, 184 insertions(+), 12 deletions(-) diff --git a/constantine/config/curves.nim b/constantine/config/curves.nim index bc8427d..db45200 100644 --- a/constantine/config/curves.nim +++ b/constantine/config/curves.nim @@ -136,6 +136,16 @@ macro genMontyMagics(T: typed): untyped = ) ) ) + # const MyCurve_InvModExponent = primeMinus2_BE(MyCurve_Modulus) + result.add newConstStmt( + ident($curve & "_InvModExponent"), newCall( + bindSym"primeMinus2_BE", + nnkDotExpr.newTree( + bindSym($curve & "_Modulus"), + ident"mres" + ) + ) + ) # echo result.toStrLit @@ -153,6 +163,10 @@ macro getMontyOne*(C: static Curve): untyped = ## Get one in Montgomery representation (i.e. R mod P) result = bindSym($C & "_MontyOne") +macro getInvModExponent*(C: static Curve): untyped = + ## Get modular inversion exponent (Modulus-2 in canonical representation) + result = bindSym($C & "_InvModExponent") + # ############################################################ # # Debug info printed at compile-time diff --git a/constantine/io/io_bigints.nim b/constantine/io/io_bigints.nim index 10f03b8..83bad0e 100644 --- a/constantine/io/io_bigints.nim +++ b/constantine/io/io_bigints.nim @@ -151,6 +151,15 @@ func fromUint*( # ############################################################ import strutils +template toByte(x: SomeUnsignedInt): byte = + ## At compile-time, conversion to bytes checks the range + ## we want to ensure this is done at the register level + ## at runtime in a single "mov byte" instruction + when nimvm: + byte(x and 0xFF) + else: + byte(x) + template blobFrom(dst: var openArray[byte], src: SomeUnsignedInt, startIdx: int, endian: static Endianness) = ## Write an integer into a raw binary blob ## Swapping endianness if needed @@ -159,10 +168,10 @@ template blobFrom(dst: var openArray[byte], src: SomeUnsignedInt, startIdx: int, when endian == cpuEndian: for i in 0 ..< sizeof(src): - dst[startIdx+i] = byte((src shr (i * 8))) + dst[startIdx+i] = toByte((src shr (i * 8))) else: for i in 0 ..< sizeof(src): - dst[startIdx+sizeof(src)-1-i] = byte((src shr (i * 8))) + dst[startIdx+sizeof(src)-1-i] = toByte((src shr (i * 8))) func exportRawUintLE( dst: var openarray[byte], @@ -203,11 +212,11 @@ func exportRawUintLE( # we can just copy each byte # tail is inclusive for i in 0 ..< tail: - dst[dst_idx+i] = byte(lo shr (i*8)) + dst[dst_idx+i] = toByte(lo shr (i*8)) else: # TODO check this # We need to copy from the end for i in 0 ..< tail: - dst[dst_idx+i] = byte(lo shr ((tail-i)*8)) + dst[dst_idx+i] = toByte(lo shr ((tail-i)*8)) return func exportRawUintBE( @@ -251,12 +260,13 @@ func exportRawUintBE( # When requesting little-endian on little-endian platform # we can just copy each byte # tail is inclusive + debugEcho "tail: " for i in 0 ..< tail: - dst[tail-i] = byte(lo shr (i*8)) + dst[tail-1-i] = toByte(lo shr (i*8)) else: # TODO check this # We need to copy from the end for i in 0 ..< tail: - dst[tail-i] = byte(lo shr ((tail-i)*8)) + dst[tail-1-i] = toByte(lo shr ((tail-i)*8)) return func exportRawUint*( diff --git a/constantine/math/bigints_checked.nim b/constantine/math/bigints_checked.nim index 0696708..a7f92a7 100644 --- a/constantine/math/bigints_checked.nim +++ b/constantine/math/bigints_checked.nim @@ -215,3 +215,29 @@ func montyPowUnsafeExponent*[mBits, eBits: static int]( scratchPtrs[i] = scratchSpace[i].view() montyPowUnsafeExponent(a.view, expBE, M.view, one.view, Word(negInvModWord), scratchPtrs) + +func montyPowUnsafeExponent*[mBits: static int]( + a: var BigInt[mBits], exponent: openarray[byte], + M, one: BigInt[mBits], negInvModWord: static BaseType, windowSize: static int) = + ## Compute a <- a^exponent (mod M) + ## ``a`` in the Montgomery domain + ## ``exponent`` is a BigInt in canonical 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, BigInt[mBits]] + var scratchPtrs {.noInit.}: array[scratchLen, BigIntViewMut] + for i in 0 ..< scratchLen: + scratchPtrs[i] = scratchSpace[i].view() + + montyPowUnsafeExponent(a.view, exponent, M.view, one.view, Word(negInvModWord), scratchPtrs) diff --git a/constantine/math/finite_fields.nim b/constantine/math/finite_fields.nim index abbcbe9..9fcf273 100644 --- a/constantine/math/finite_fields.nim +++ b/constantine/math/finite_fields.nim @@ -134,7 +134,11 @@ func pow*(a: var Fq, exponent: BigInt) = ## ``a``: a field element to be exponentiated ## ``exponent``: a big integer const windowSize = 5 # TODO: find best window size for each curves - a.mres.montyPow(exponent, Fq.C.Mod.mres, Fq.C.getMontyOne(), Fq.C.getNegInvModWord(), windowSize) + a.mres.montyPow( + exponent, + Fq.C.Mod.mres, Fq.C.getMontyOne(), + Fq.C.getNegInvModWord(), windowSize + ) func powUnsafeExponent*(a: var Fq, exponent: BigInt) = ## Exponentiation over Fq @@ -148,4 +152,20 @@ func powUnsafeExponent*(a: var Fq, exponent: BigInt) = ## - power analysis ## - timing analysis const windowSize = 5 # TODO: find best window size for each curves - a.mres.montyPowUnsafeExponent(exponent, Fq.C.Mod.mres, Fq.C.getMontyOne(), Fq.C.getNegInvModWord(), windowSize) + a.mres.montyPowUnsafeExponent( + exponent, + Fq.C.Mod.mres, Fq.C.getMontyOne(), + Fq.C.getNegInvModWord(), windowSize + ) + +func inv*(a: var Fq) = + ## Modular inversion + ## Warning ⚠️ : + ## - This assumes that `Fq` is a prime field + + const windowSize = 5 # TODO: find best window size for each curves + a.mres.montyPowUnsafeExponent( + Fq.C.getInvModExponent(), + Fq.C.Mod.mres, Fq.C.getMontyOne(), + Fq.C.getNegInvModWord(), windowSize + ) diff --git a/constantine/math/precomputed.nim b/constantine/math/precomputed.nim index 011aeb2..c760903 100644 --- a/constantine/math/precomputed.nim +++ b/constantine/math/precomputed.nim @@ -9,7 +9,8 @@ import ./bigints_checked, ../primitives/constant_time, - ../config/common + ../config/common, + ../io/io_bigints # Precomputed constants # ############################################################ @@ -187,3 +188,15 @@ func montyOne*(M: BigInt): BigInt = ## Returns "1 (mod M)" in the Montgomery domain. ## This is equivalent to R (mod M) in the natural domain r_powmod(1, M) + +func primeMinus2_BE*[bits: static int]( + P: BigInt[bits] + ): array[(bits+7) div 8, byte] {.noInit.} = + ## Compute an input prime-2 + ## and return the result as a canonical byte array / octet string + ## For use to precompute modular inverse exponent. + + var tmp = P + discard tmp.sub(BigInt[bits].fromRawUint([byte 2], bigEndian), true) + + result.exportRawUint(tmp, bigEndian) diff --git a/tests/test_finite_fields_powinv.nim b/tests/test_finite_fields_powinv.nim index 1bb2ef3..3c710b8 100644 --- a/tests/test_finite_fields_powinv.nim +++ b/tests/test_finite_fields_powinv.nim @@ -125,4 +125,18 @@ proc main() = check: computed == expected + suite "Modular inversion over prime fields": + test "x^(-1) mod p": + var x: Fq[BLS12_381] + + # BN254 field modulus + x.fromHex("0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47") + + let expected = "0x0636759a0f3034fa47174b2c0334902f11e9915b7bd89c6a2b3082b109abbc9837da17201f6d8286fe6203caa1b9d4c8" + x.inv() + let computed = x.toHex() + + check: + computed == expected + main() diff --git a/tests/test_finite_fields_vs_gmp.nim b/tests/test_finite_fields_vs_gmp.nim index a327b0f..8ff77a0 100644 --- a/tests/test_finite_fields_vs_gmp.nim +++ b/tests/test_finite_fields_vs_gmp.nim @@ -47,7 +47,7 @@ const # https://gmplib.org/manual/Integer-Import-and-Export.html GMP_MostSignificantWordFirst = 1'i32 GMP_LeastSignificantWordFirst = -1'i32 -proc main() = +proc mainMul() = var gmpRng: gmp_randstate_t gmp_randinit_mt(gmpRng) # The GMP seed varies between run so that @@ -66,7 +66,7 @@ proc main() = randomTests(128, curve): # echo "--------------------------------------------------------------------------------" - echo "Testing: random input on ", $curve + echo "Testing: random modular multiplication on ", $curve const bits = CurveParams[curve][0] @@ -127,4 +127,79 @@ proc main() = " GMP: " & rGMP.toHex() & "\n" & " Constantine: " & rConstantine.toHex() -main() +proc mainInv() = + var gmpRng: gmp_randstate_t + gmp_randinit_mt(gmpRng) + # The GMP seed varies between run so that + # test coverage increases as the library gets tested. + # This requires to dump the seed in the console or the function inputs + # to be able to reproduce a bug + let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32 + echo "GMP seed: ", seed + gmp_randseed_ui(gmpRng, seed) + + var a, p, r: mpz_t + mpz_init(a) + mpz_init(p) + mpz_init(r) + + randomTests(128, curve): + # echo "--------------------------------------------------------------------------------" + echo "Testing: random modular inversion on ", $curve + + const bits = CurveParams[curve][0] + + # Generate random value in the range 0 ..< 2^(bits-1) + mpz_urandomb(a, gmpRng, uint bits) + # Set modulus to curve modulus + let err = mpz_set_str(p, CurveParams[curve][1], 0) + doAssert err == 0 + + ######################################################### + # Conversion buffers + const len = csize (bits + 7) div 8 + + # Note: GMP does not pad right the bigendian numbers if there is extra space + var aBuf: array[len, byte] + + var aW: csize # Word written by GMP + discard mpz_export(aBuf[0].addr, aW.addr, GMP_MostSignificantWordFirst, 1, GMP_WordNativeEndian, 0, a) + + # Since the modulus is using all bits, it's we can test for exact amount copy + doAssert len >= aW, "Expected at most " & $len & " bytes but wrote " & $aW & " for " & toHex(aBuf) & " (little-endian)" + + # Build the bigint - TODO more fields codecs + let aTest = Fq[curve].fromBig BigInt[bits].fromRawUint(aBuf[0 ..< aW], bigEndian) + + ######################################################### + # Modular inversion + let exist = mpz_invert(r, a, p) + doAssert exist != 0 + + var rTest = aTest + rTest.inv() + + ######################################################### + # Check + var rGMP: array[len, byte] + var rW: csize # Word written by GMP + discard mpz_export(rGMP[0].addr, rW.addr, GMP_MostSignificantWordFirst, 1, GMP_WordNativeEndian, 0, r) + + var rConstantine: array[len, byte] + exportRawUint(rConstantine, rTest, bigEndian) + + # echo "rGMP: ", rGMP.toHex() + # echo "rConstantine: ", rConstantine.toHex() + + doAssert rGMP[0 ..< rW] == rConstantine[^rW..^1], block: + # Reexport as bigEndian for debugging + discard mpz_export(aBuf[0].addr, aW.addr, GMP_MostSignificantWordFirst, 1, GMP_WordNativeEndian, 0, a) + "\nModular Inversion on curve " & $curve & " with operand\n" & + " a: 0x" & aBuf.toHex & "\n" & + " p: " & CurveParams[curve][1] & "\n" & + "failed:" & "\n" & + " GMP: " & rGMP.toHex() & "\n" & + " Constantine: " & rConstantine.toHex() + +mainMul() +mainInv()