diff --git a/constantine.nimble b/constantine.nimble index b291c7b..84b5964 100644 --- a/constantine.nimble +++ b/constantine.nimble @@ -21,3 +21,4 @@ task test, "Run all tests": test "", "tests/test_bigints.nim" test "", "tests/test_bigints_multimod.nim" test "", "tests/test_bigints_vs_gmp.nim" + test "", "tests/test_finite_fields.nim" diff --git a/constantine/config/curves.nim b/constantine/config/curves.nim index 585bf1f..5009287 100644 --- a/constantine/config/curves.nim +++ b/constantine/config/curves.nim @@ -18,57 +18,51 @@ import # # ############################################################ -func montyMagic(M: static BigInt): static Word {.inline.} = +# Note: The declareCurve macro builds an exported montyMagic*(C: static Curve) on top of the one below +func montyMagic(M: BigInt): BaseType = ## Returns the Montgomery domain magic constant for the input modulus: ## -1/M[0] mod LimbSize ## M[0] is the least significant limb of M ## M must be odd and greater than 2. + # We use BaseType for return value because static distinct type + # confuses Nim semchecks [UPSTREAM BUG] + # We don't enforce compile-time evaluation here + # because static BigInt[bits] also causes semcheck troubles [UPSTREAM BUG] - # Test vectors: https://www.researchgate.net/publication/4107322_Montgomery_modular_multiplication_architecture_for_public_key_cryptosystems - # on p354 - # Reference C impl: http://www.hackersdelight.org/hdcodetxt/mont64.c.txt - # Explanation: https://eprint.iacr.org/2017/1057.pdf - - # ###################################################################### - # Implementation of modular multiplicative inverse - # Assuming 2 positive integers a and m the modulo - # - # We are looking for z that solves `az ≡ 1 mod m` - # - # References: - # - Knuth, The Art of Computer Programming, Vol2 p342 - # - Menezes, Handbook of Applied Cryptography (HAC), p610 - # http://cacr.uwaterloo.ca/hac/about/chap14.pdf - - # Starting from the extended GCD formula (Bezout identity), - # `ax + by = gcd(x,y)` with input x,y and outputs a, b, gcd - # We assume a and m are coprimes, i.e. gcd is 1, otherwise no inverse - # `ax + my = 1` <=> `ax + my ≡ 1 mod m` <=> `ax ≡ 1 mod m` + # Modular inverse algorithm: + # 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 # For Montgomery magic number, we are in a special case # where a = M and m = 2^LimbSize. # For a and m to be coprimes, a must be odd. - # `m` (2^LimbSize) being a power of 2 greatly simplifies computation: - # - https://crypto.stackexchange.com/questions/47493/how-to-determine-the-multiplicative-inverse-modulo-64-or-other-power-of-two - # - http://groups.google.com/groups?selm=1994Apr6.093116.27805%40mnemosyne.cs.du.edu - # - https://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two - # - https://eprint.iacr.org/2017/411 - # We have the following relation # ax ≡ 1 (mod 2^k) <=> ax(2 - ax) ≡ 1 (mod 2^(2k)) # # To get -1/M0 mod LimbSize # we can either negate the resulting x of `ax(2 - ax) ≡ 1 (mod 2^(2k))` # or do ax(2 + ax) ≡ 1 (mod 2^(2k)) + # + # To get the the modular inverse of 2^k' with arbitrary k' (like k=63 in our case) + # we can do modInv(a, 2^64) mod 2^63 as mentionned in Koc paper. - const - M0 = M.limbs[0] - k = log2(WordBitSize) + let + M0 = BaseType(M.limbs[0]) + k = log2(WordPhysBitSize) - result = M0 # Start from an inverse of M0 modulo 2, M0 is odd and it's own inverse - for _ in 0 ..< k: - result *= 2 + M * result # x' = x(2 + ax) (`+` to avoid negating at the end) + 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) (`+` to avoid negating at the end) + + # Our actual word size is 2^63 not 2^64 + result = result and BaseType(MaxWord) # ############################################################ # diff --git a/constantine/config/curves_parser.nim b/constantine/config/curves_parser.nim index 1087b67..1edf1b5 100644 --- a/constantine/config/curves_parser.nim +++ b/constantine/config/curves_parser.nim @@ -139,7 +139,7 @@ macro declareCurves*(curves: untyped): untyped = let C = ident"C" let matchingBigInt = ident"matchingBigInt" result.add newProc( - name = matchingBigInt, + name = nnkPostFix.newTree(ident"*", matchingBigInt), params = [ident"untyped", newIdentDefs(C, nnkStaticTy.newTree(Curve))], body = nnkBracketExpr.newTree(bindSym"BigInt", nnkBracketExpr.newTree(cbs, C)), procType = nnkTemplateDef @@ -158,16 +158,13 @@ macro declareCurves*(curves: untyped): untyped = nnkGenericParams.newTree(newIdentDefs( C, nnkStaticTy.newTree(Curve), newEmptyNode() )), + # TODO: where should I put the nnkCommentStmt? nnkObjectTy.newTree( newEmptyNode(), newEmptyNode(), nnkRecList.newTree( - nnkCommentStmt.newTree "All operations on a field are modulo P", - nnkCommentStmt.newTree "P being the prime modulus of the Curve C", - nnkCommentStmt.newTree "Internally, data is stored in Montgomery n-residue form", - nnkCommentStmt.newTree "with the magic constant chosen for convenient division (a power of 2 depending on P bitsize)" newIdentDefs( - nnkPostfix.newTree(ident"*", ident"value"), + nnkPostfix.newTree(ident"*", ident"nres"), newCall(matchingBigInt, C) ) ) @@ -203,7 +200,7 @@ macro declareCurves*(curves: untyped): untyped = # proc MontyMagic(curve: static Curve): static Word result.add newProc( - name = nnkPostfix.newTree(ident"*", ident"MontyMagic"), + name = nnkPostfix.newTree(ident"*", ident"montyMagic"), params = [ ident"auto", newIdentDefs( @@ -213,10 +210,14 @@ macro declareCurves*(curves: untyped): untyped = ], body = newCall( ident"montyMagic", - newCall(ident"Mod", ident"curve") + # curve.Mod().nres + nnkDotExpr.newTree( + newCall(ident"Mod", ident"curve"), + ident"nres" + ) + ), - procType = nnkFuncDef, - pragmas = nnkPragma.newTree(ident"compileTime") + procType = nnkFuncDef ) # echo result.toStrLit() diff --git a/constantine/io/io_bigints.nim b/constantine/io/io_bigints.nim index da2bf44..bf0be99 100644 --- a/constantine/io/io_bigints.nim +++ b/constantine/io/io_bigints.nim @@ -181,13 +181,14 @@ func serializeRawUint*( dstEndianness: static Endianness) = ## Serialize a bigint into its canonical big-endian or little endian ## representation. - ## A destination buffer of size "BigInt.bits div 8" at minimum is needed. + ## A destination buffer of size "(BigInt.bits + 7) div 8" at minimum is needed, + ## i.e. bits -> byte conversion rounded up ## ## If the buffer is bigger, output will be zero-padded left for big-endian ## or zero-padded right for little-endian. ## I.e least significant bit is aligned to buffer boundary - assert dst.len >= static(BigInt.bits div 8), "BigInt -> Raw int conversion: destination buffer is too small" + assert dst.len >= (BigInt.bits + 7) div 8, "BigInt -> Raw int conversion: destination buffer is too small" when BigInt.bits == 0: zeroMem(dst, dst.len) @@ -339,16 +340,7 @@ func toHex*(big: BigInt, order: static Endianness = bigEndian): string = ## Note. Leading zeros are not removed. ## Result is prefixed with 0x ## - ## This is a raw memory dump. Output will be padded with 0 - ## if the big int does not use the full memory allocated for it. - ## - ## Regardless of the machine endianness the output will be big-endian hex. - ## - ## For example a BigInt representing 10 will be - ## - 0x0a for BigInt[8] - ## - 0x000a for BigInt[16] - ## - 0x00000000_0000000a for BigInt[64] - ## (underscore added for docuentation readability only) + ## Output will be padded with 0s to maintain constant-time. ## ## CT: ## - no leaks diff --git a/constantine/io/io_fields.nim b/constantine/io/io_fields.nim index c2b0c27..33d7822 100644 --- a/constantine/io/io_fields.nim +++ b/constantine/io/io_fields.nim @@ -8,7 +8,8 @@ import ./io_bigints, - ../math/finite_fields + ../config/curves, + ../math/[bigints_checked, finite_fields] # ############################################################ # @@ -16,8 +17,34 @@ import # # ############################################################ -func fromUint*(dst: var Fp, +func fromUint*(dst: var Fq, src: SomeUnsignedInt) = ## Parse a regular unsigned integer ## and store it into a BigInt of size `bits` - dst.value.fromRawUint(cast[array[sizeof(src), byte]](src), cpuEndian) + dst.nres.fromRawUint(cast[array[sizeof(src), byte]](src), cpuEndian) + dst.nres.unsafeMontgomeryResidue(Fq.C.Mod.nres) + +func serializeRawUint*(dst: var openarray[byte], + src: Fq, + dstEndianness: static Endianness) = + ## Serialize a finite field element to its canonical big-endian or little-endian + ## representation + ## With `bits` the number of bits of the field modulus + ## a buffer of size "(bits + 7) div 8" at minimum is needed + ## i.e. bits -> byte conversion rounded up + ## + ## If the buffer is bigger, output will be zero-padded left for big-endian + ## or zero-padded right for little-endian. + ## I.e least significant bit is aligned to buffer boundary + serializeRawUint(dst, src.toBig(), dstEndianness) + +func toHex*(f: Fq, order: static Endianness = bigEndian): string = + ## Stringify a finite field element to hex. + ## Note. Leading zeros are not removed. + ## Result is prefixed with 0x + ## + ## Output will be padded with 0s to maintain constant-time. + ## + ## CT: + ## - no leaks + result = f.toBig().toHex() diff --git a/constantine/math/bigints_checked.nim b/constantine/math/bigints_checked.nim index ab6243a..4185683 100644 --- a/constantine/math/bigints_checked.nim +++ b/constantine/math/bigints_checked.nim @@ -110,3 +110,24 @@ func reduce*[aBits, mBits](r: var BigInt[mBits], a: BigInt[aBits], M: BigInt[mBi # but we don't want to inline it as it would increase codesize, better have Nim # pass a pointer+length to a fixed session of the BSS. reduce(r.view, a.view, M.view) + +func unsafeMontgomeryResidue*[mBits](nres: var BigInt[mBits], N: BigInt[mBits]) = + ## Convert a BigInt from its natural representation + ## to the Montgomery n-residue form + ## + ## `nres` is modified in-place + ## + ## Caller must take care of properly switching between + ## the natural and montgomery domain. + ## Nesting Montgomery form is possible by applying this function twice. + montgomeryResidue(nres.view, N.view) + +func unsafeRedc*[mBits](nres: var BigInt[mBits], N: BigInt[mBits], montyMagic: static BaseType) = + ## Convert a BigInt from its Montgomery n-residue form + ## to the natural representation + ## + ## `nres` is modified in-place + ## + ## Caller must take care of properly switching between + ## the natural and montgomery domain. + redc(nres.view, N.view, Word(montyMagic)) diff --git a/constantine/math/bigints_raw.nim b/constantine/math/bigints_raw.nim index a1d78f1..31f6802 100644 --- a/constantine/math/bigints_raw.nim +++ b/constantine/math/bigints_raw.nim @@ -181,6 +181,12 @@ template checkValidModulus(m: BigIntViewConst) = debug: assert not m[^1].isZero.bool, "Internal Error: the modulus must use all declared bits" +template checkOddModulus(m: BigIntViewConst) = + ## CHeck that the modulus is odd + ## and valid for use in the Montgomery n-residue representation + debug: + assert bool(BaseType(m[0]) and 1), "Internal Error: the modulus must be odd to use the Montgomery representation." + debug: func `$`*(a: BigIntViewAny): string = let len = a.numLimbs() @@ -375,3 +381,56 @@ func reduce*(r: BigIntViewMut, a: BigIntViewAny, M: BigIntViewConst) = # Now shift-left the copied words while adding the new word modulo M for i in countdown(aOffset, 0): r.shlAddMod(a[i], M) + +func montgomeryResidue*(a: BigIntViewMut, N: BigIntViewConst) = + ## Transform a bigint ``a`` from it's natural representation (mod N) + ## to a the Montgomery n-residue representation + ## i.e. Does "a * (2^LimbSize)^W (mod N), where W is the number + ## of words needed to represent n in base 2^LimbSize + ## + ## `a`: The source BigInt in the natural representation. `a` in [0, N) range + ## `N`: The field modulus. N must be odd. + ## + ## Important: `a` is overwritten + # Reference: https://eprint.iacr.org/2017/1057.pdf + checkValidModulus(N) + checkOddModulus(N) + checkMatchingBitlengths(a, N) + + let nLen = N.numLimbs() + for i in countdown(nLen, 1): + a.shlAddMod(Zero, N) + +func redc*(a: BigIntViewMut, N: BigIntViewConst, montyMagic: Word) = + ## Transform a bigint ``a`` from it's Montgomery N-residue representation (mod N) + ## to the regular natural representation (mod N) + ## + ## i.e. Does "a * ((2^LimbSize)^W)^-1 (mod N), where W is the number + ## of words needed to represent n in base 2^LimbSize + ## + ## This is called a Montgomery Reduction + ## The Montgomery Magic Constant is µ = -1/N mod N + ## is used internally and can be precomputed with montyMagic(Curve) + # References: + # - https://eprint.iacr.org/2017/1057.pdf (Montgomery) + # page: Radix-r interleaved multiplication algorithm + # - https://en.wikipedia.org/wiki/Montgomery_modular_multiplication#Montgomery_arithmetic_on_multiprecision_(variable-radix)_integers + # - http://langevin.univ-tln.fr/cours/MLC/extra/montgomery.pdf + # Montgomery original paper + # + checkValidModulus(N) + checkOddModulus(N) + checkMatchingBitlengths(a, N) + + let nLen = N.numLimbs() + for i in 0 ..< nLen: + let z0 = Word(BaseType(a[0]) * BaseType(montyMagic)) and MaxWord + + var carry = DoubleWord(0) + for j in 0 ..< nLen: + let z = DoubleWord(a[i]) + unsafeExtPrecMul(z0, N[i]) + carry + carry = z shr WordBitSize + if j != 0: + a[j] = Word(z) and MaxWord + + a[^1] = Word(carry) diff --git a/constantine/math/finite_fields.nim b/constantine/math/finite_fields.nim index 572577f..db2a7a5 100644 --- a/constantine/math/finite_fields.nim +++ b/constantine/math/finite_fields.nim @@ -42,16 +42,21 @@ debug: # No exceptions allowed {.push raises: [].} -func toMonty*[C: static Curve](a: Fq[C]): Montgomery[C] = - ## Convert a big integer over Fq to its montgomery representation - ## over Fq. - ## i.e. Does "a * (2^LimbSize)^W (mod p), where W is the number - ## of words needed to represent p in base 2^LimbSize +# ############################################################ +# +# Conversion +# +# ############################################################ - result = a - for i in static(countdown(C.Mod.limbs.high, 1)): - shiftAdd(result, 0) +func fromBig*(T: type Fq, src: BigInt): T = + ## Convert a BigInt to its Montgomery form + result.nres = src + result.nres.unsafeMontgomeryResidue(Fq.C.Mod) +func toBig*[C: static Curve](src: Fq[C]): auto = + ## Convert a finite-field element to a BigInt in natral representation + result = src.nres + result.unsafeRedC(C.Mod.nres, montyMagic(C)) # ############################################################ # @@ -66,7 +71,7 @@ template add(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] = ## ## a and b MAY be the same buffer ## a and b MUST have the same announced bitlength (i.e. `bits` static parameters) - add(a.value, b.value, ctl) + add(a.nres, b.nres, ctl) template sub(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] = ## Constant-time big integer in-place optional substraction @@ -75,7 +80,7 @@ template sub(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] = ## ## a and b MAY be the same buffer ## a and b MUST have the same announced bitlength (i.e. `bits` static parameters) - sub(a.value, b.value, ctl) + sub(a.nres, b.nres, ctl) # ############################################################ # @@ -83,9 +88,6 @@ template sub(a: var Fq, b: Fq, ctl: CTBool[Word]): CTBool[Word] = # # ############################################################ -# No exceptions allowed -{.push raises: [].} - func `+=`*(a: var Fq, b: Fq) = ## Addition over Fq var ctl = add(a, b, CtTrue) diff --git a/tests/test_finite_fields.nim b/tests/test_finite_fields.nim index 3b1bfd3..c54d342 100644 --- a/tests/test_finite_fields.nim +++ b/tests/test_finite_fields.nim @@ -11,6 +11,8 @@ import unittest, random, ../constantine/io/io_fields, ../constantine/config/curves +import ../constantine/io/io_bigints + static: doAssert defined(testingCurves), "This modules requires the -d:testingCurves compile option" proc main() = @@ -24,7 +26,15 @@ proc main() = z.fromUint(90'u32) x += y - check: bool(z == x) + + var x_bytes: array[8, byte] + x_bytes.serializeRawUint(x, cpuEndian) + + check: + # Check equality in the Montgomery domain + bool(z == x) + # Check equality when converting back to natural domain + 90'u64 == cast[uint64](x_bytes) block: var x, y, z: Fq[Fake101] @@ -44,7 +54,15 @@ proc main() = z.fromUint(1'u32) x += y - check: bool(z == x) + + var x_bytes: array[8, byte] + x_bytes.serializeRawUint(x, cpuEndian) + + check: + # Check equality in the Montgomery domain + bool(z == x) + # Check equality when converting back to natural domain + 1'u64 == cast[uint64](x_bytes) test "Substraction mod 101": block: @@ -55,7 +73,15 @@ proc main() = z.fromUint(70'u32) x -= y - check: bool(z == x) + + var x_bytes: array[8, byte] + x_bytes.serializeRawUint(x, cpuEndian) + + check: + # Check equality in the Montgomery domain + bool(z == x) + # Check equality when converting back to natural domain + 70'u64 == cast[uint64](x_bytes) block: var x, y, z: Fq[Fake101] @@ -65,7 +91,15 @@ proc main() = z.fromUint(0'u32) x -= y - check: bool(z == x) + + var x_bytes: array[8, byte] + x_bytes.serializeRawUint(x, cpuEndian) + + check: + # Check equality in the Montgomery domain + bool(z == x) + # Check equality when converting back to natural domain + 0'u64 == cast[uint64](x_bytes) block: var x, y, z: Fq[Fake101] @@ -75,6 +109,14 @@ proc main() = z.fromUint(100'u32) x -= y - check: bool(z == x) + + var x_bytes: array[8, byte] + x_bytes.serializeRawUint(x, cpuEndian) + + check: + # Check equality in the Montgomery domain + bool(z == x) + # Check equality when converting back to natural domain + 100'u64 == cast[uint64](x_bytes) main()