Fast modular inversion

benchmarks/bench_elliptic_template.nim

@@ -68,10 +68,7 @@ proc mixedAddBench*(T: typedesc, iters: int) =
   let P = rng.random_unsafe(T)
   let Q = rng.random_unsafe(T)
   var Qaff: ECP_ShortW_Aff[T.F, T.G]
-  when Q is ECP_ShortW_Prj:
-    Qaff.affineFromProjective(Q)
-  else:
-    Qaff.affineFromJacobian(Q)
+  Qaff.affine(Q)
   bench("EC Mixed Addition " & G1_or_G2, T, iters):
     r.madd(P, Qaff)
 
@@ -87,14 +84,14 @@ proc affFromProjBench*(T: typedesc, iters: int) =
   var r {.noInit.}: ECP_ShortW_Aff[T.F, T.G]
   let P = rng.random_unsafe(T)
   bench("EC Projective to Affine " & G1_or_G2, T, iters):
-    r.affineFromProjective(P)
+    r.affine(P)
 
 proc affFromJacBench*(T: typedesc, iters: int) =
   const G1_or_G2 = when T.F is Fp: "G1" else: "G2"
   var r {.noInit.}: ECP_ShortW_Aff[T.F, T.G]
   let P = rng.random_unsafe(T)
   bench("EC Jacobian to Affine " & G1_or_G2, T, iters):
-    r.affineFromJacobian(P)
+    r.affine(P)
 
 proc scalarMulGenericBench*(T: typedesc, window: static int, iters: int) =
   const bits = T.F.C.getCurveOrderBitwidth()

benchmarks/bench_fields_template.nim

@@ -95,30 +95,9 @@ proc invBench*(T: typedesc, iters: int) =
   var r: T
   let x = rng.random_unsafe(T)
   preventOptimAway(r)
-  bench("Inversion (constant-time default impl)", T, iters):
+  bench("Inversion (constant-time)", T, iters):
     r.inv(x)
 
-proc invEuclidBench*(T: typedesc, iters: int) =
-  var r: T
-  let x = rng.random_unsafe(T)
-  preventOptimAway(r)
-  bench("Inversion (constant-time Euclid)", T, iters):
-    r.inv_euclid(x)
-
-proc invPowFermatBench*(T: typedesc, iters: int) =
-  let x = rng.random_unsafe(T)
-  const exponent = T.getInvModExponent()
-  bench("Inversion (exponentiation p-2, Little Fermat)", T, iters):
-    var r = x
-    r.powUnsafeExponent(exponent)
-
-proc invAddChainBench*(T: typedesc, iters: int) =
-  var r: T
-  let x = rng.random_unsafe(T)
-  preventOptimAway(r)
-  bench("Inversion (addition chain)", T, iters):
-    r.inv_addchain(x)
-
 proc sqrtBench*(T: typedesc, iters: int) =
   let x = rng.random_unsafe(T)
 

benchmarks/bench_fp.nim

@@ -11,7 +11,7 @@ import
   ../constantine/config/[curves, common],
   ../constantine/arithmetic,
   ../constantine/io/io_bigints,
-  ../constantine/curves/[zoo_inversions, zoo_square_roots],
+  ../constantine/curves/zoo_square_roots,
   # Helpers
   ../helpers/static_for,
   ./bench_fields_template
@@ -50,8 +50,7 @@ proc main() =
     mulBench(Fp[curve], Iters)
     sqrBench(Fp[curve], Iters)
     smallSeparator()
-    invEuclidBench(Fp[curve], ExponentIters)
-    invPowFermatBench(Fp[curve], ExponentIters)
+    invBench(Fp[curve], ExponentIters)
     sqrtBench(Fp[curve], ExponentIters)
     sqrtRatioBench(Fp[curve], ExponentIters)
     # Exponentiation by a "secret" of size ~the curve order

benchmarks/bench_hash_to_curve.nim

@@ -64,7 +64,7 @@ proc bench_BLS12_381_proj_aff_conversion(iters: int) =
   )
 
   bench("Proj->Affine conversion (for pairing)", BLS12_381, iters):
-    Paff.affineFromProjective(P)
+    Paff.affine(P)
 
 const Iters = 1000
 

benchmarks/bench_pairing_template.nim

@@ -51,9 +51,9 @@ func clearCofactor[F; G: static Subgroup](
        ec: var ECP_ShortW_Aff[F, G]) =
   # For now we don't have any affine operation defined
   var t {.noInit.}: ECP_ShortW_Prj[F, G]
-  t.projectiveFromAffine(ec)
+  t.fromAffine(ec)
   t.clearCofactor()
-  ec.affineFromProjective(t)
+  ec.affine(t)
 
 func random_point*(rng: var RngState, EC: typedesc): EC {.noInit.} =
   result = rng.random_unsafe(EC)

benchmarks/bench_summary_template.nim

@@ -86,9 +86,9 @@ func clearCofactorReference[F; G: static Subgroup](
        ec: var ECP_ShortW_Aff[F, G]) =
   # For now we don't have any affine operation defined
   var t {.noInit.}: ECP_ShortW_Prj[F, G]
-  t.projectiveFromAffine(ec)
+  t.fromAffine(ec)
   t.clearCofactorReference()
-  ec.affineFromProjective(t)
+  ec.affine(t)
 
 func random_point*(rng: var RngState, EC: typedesc): EC {.noInit.} =
   result = rng.random_unsafe(EC)
@@ -136,10 +136,7 @@ proc mixedAddBench*(T: typedesc, iters: int) =
   let P = rng.random_unsafe(T)
   let Q = rng.random_unsafe(T)
   var Qaff: ECP_ShortW_Aff[T.F, T.G]
-  when Q is ECP_ShortW_Prj:
-    Qaff.affineFromProjective(Q)
-  else:
-    Qaff.affineFromJacobian(Q)
+  Qaff.affine(Q)
   bench("EC Mixed Addition " & G1_or_G2, T, iters):
     r.madd(P, Qaff)
 

constantine.nimble

@@ -28,6 +28,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # Big ints
   # ----------------------------------------------------------
   ("tests/t_io_bigints.nim", false),
+  ("tests/t_io_unsaturated.nim", false),
   ("tests/t_bigints.nim", false),
   ("tests/t_bigints_multimod.nim", false),
   ("tests/t_bigints_mod_vs_gmp.nim", true),
@@ -64,6 +65,11 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   ("tests/t_fp4_frobenius.nim", false),
   ("tests/t_fp6_frobenius.nim", false),
   ("tests/t_fp12_frobenius.nim", false),
+
+  # Elliptic curve arithmetic 
+  # ----------------------------------------------------------
+  ("tests/t_ec_conversion.nim", false),
+
   # Elliptic curve arithmetic G1
   # ----------------------------------------------------------
   # ("tests/t_ec_shortw_prj_g1_add_double.nim", false),

constantine/arithmetic.nim

@@ -10,7 +10,6 @@ import
   arithmetic/[bigints, bigints_montgomery],
   arithmetic/[
     finite_fields,
-    finite_fields_inversion,
     finite_fields_square_root,
     finite_fields_double_precision
   ]
@@ -19,6 +18,5 @@ export
   bigints,
   bigints_montgomery,
   finite_fields,
-  finite_fields_inversion,
   finite_fields_square_root,
   finite_fields_double_precision

constantine/arithmetic/assembly/limbs_asm_montmul_x86_adx_bmi2.nim

@@ -147,7 +147,7 @@ proc partialRedx(
     ctx.comment "  Reduction"
     ctx.comment "  m = t[0] * m0ninv mod 2^w"
     ctx.mov  rdx, t[0]
-    ctx.mulx S, rdx, m0ninv, rdx # (S, RDX) <- m0ninv * RDX
+    ctx.imul rdx, m0ninv
 
     # Clear carry flags - TODO: necessary?
     ctx.`xor` S, S

constantine/arithmetic/bigints.nim

@@ -11,8 +11,8 @@ import
   ../primitives,
   ./limbs,
   ./limbs_extmul,
-  ./limbs_modular,
-  ./limbs_montgomery
+  ./limbs_invmod,
+  ./limbs_division
 
 export BigInt
 
@@ -79,9 +79,9 @@ 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) =
+func csetZero*(a: var BigInt, ctl: SecretBool) =
   ## Set ``a`` to 0 if ``ctl`` is true
-  a.limbs.czero(ctl)
+  a.limbs.csetZero(ctl)
 
 # Copy
 # ------------------------------------------------------------
@@ -450,35 +450,36 @@ func reduce*[aBits, mBits](r: var BigInt[mBits], a: BigInt[aBits], M: BigInt[mBi
   # pass a pointer+length to a fixed session of the BSS.
   reduce(r.limbs, a.limbs, aBits, M.limbs, mBits)
 
-func div2_modular*[bits](a: var BigInt[bits], mp1div2: BigInt[bits]) =
-  ## Compute a <- a/2 (mod M)
-  ## `mp1div2` is the modulus (M+1)/2
-  ##
-  ## Normally if `a` is odd we add the modulus before dividing by 2
-  ## but this may overflow and we might lose a bit before shifting.
-  ## Instead we shift first and then add half the modulus rounded up
+func invmod*[bits](
+       r: var BigInt[bits],
+       a, F, M: BigInt[bits]) =
+  ## Compute the modular inverse of ``a`` modulo M
+  ## r ≡ F.a⁻¹ (mod M)
   ##
-  ## Assuming M is odd, `mp1div2` can be precomputed without
-  ## overflowing the "Limbs" by dividing by 2 first
-  ## and add 1
-  ## Otherwise `mp1div2` should be M/2
-  a.limbs.div2_modular(mp1div2.limbs)
-
-func steinsGCD*[bits](r: var BigInt[bits], a, F, M, mp1div2: BigInt[bits]) =
-  ## Compute F multiplied the modular inverse of ``a`` modulo M
-  ## r ≡ F . a^-1 (mod M)
+  ## M MUST be odd, M does not need to be prime.
+  ## ``a`` MUST be less than M.
+  r.limbs.invmod(a.limbs, F.limbs, M.limbs, bits)
+
+func invmod*[bits](
+       r: var BigInt[bits],
+       a: BigInt[bits],
+       F, M: static BigInt[bits]) =
+  ## Compute the modular inverse of ``a`` modulo M
+  ## r ≡ F.a⁻¹ (mod M)
+  ## 
+  ## with F and M known at compile-time
   ##
   ## M MUST be odd, M does not need to be prime.
   ## ``a`` MUST be less than M.
-  r.limbs.steinsGCD(a.limbs, F.limbs, M.limbs, bits, mp1div2.limbs)
+  r.limbs.invmod(a.limbs, F.limbs, M.limbs, bits)
 
-func invmod*[bits](r: var BigInt[bits], a, M, mp1div2: BigInt[bits]) =
+func invmod*[bits](r: var BigInt[bits], a, M: BigInt[bits]) =
   ## Compute the modular inverse of ``a`` modulo M
   ##
   ## The modulus ``M`` MUST be odd
   var one {.noInit.}: BigInt[bits]
   one.setOne()
-  r.steinsGCD(a, one, M, mp1div2)
+  r.invmod(a, one, M)
 
 {.pop.} # inline
 {.pop.} # raises no exceptions

constantine/arithmetic/bigints_montgomery.nim

@@ -11,7 +11,6 @@ import
   ../primitives,
   ../io/io_bigints,
   ./limbs,
-  ./limbs_modular,
   ./limbs_montgomery,
   ./bigints
 

constantine/arithmetic/finite_fields.nim

@@ -231,25 +231,29 @@ func neg*(r: var FF, a: FF) {.meter.} =
     var t {.noInit.}: FF
     let isZero = a.isZero()
     discard t.mres.diff(FF.fieldMod(), a.mres)
-    t.mres.czero(isZero)
+    t.mres.csetZero(isZero)
     r = t
 
 func neg*(a: var FF) {.meter.} =
   ## Negate modulo p
   a.neg(a)
 
-func div2*(a: var FF) {.meter.} =
-  ## Modular division by 2
-  a.mres.div2_modular(FF.getPrimePlus1div2())
-
 # ############################################################
 #
 #                Field arithmetic conditional
 #
 # ############################################################
 
+func csetZero*(a: var FF, ctl: SecretBool) =
+  ## Set ``a`` to 0 if ``ctl`` is true
+  a.mres.csetZero(ctl)
+
+func csetOne*(a: var FF, ctl: SecretBool) =
+  ## Set ``a`` to 1 if ``ctl`` is true
+  a.mres.ccopy(FF.getMontyOne(), ctl)
+
 func cneg*(r: var FF, a: FF, ctl: SecretBool) {.meter.} =
-  ## Constant-time in-place conditional negation
+  ## Constant-time out-of-place conditional negation
   ## The negation is only performed if ctl is "true"
   r.neg(a)
   r.ccopy(a, not ctl)
@@ -261,19 +265,68 @@ func cneg*(a: var FF, ctl: SecretBool) {.meter.} =
   a.cneg(t, ctl)
 
 func cadd*(a: var FF, b: FF, ctl: SecretBool) {.meter.} =
-  ## Constant-time in-place conditional addition
+  ## Constant-time out-place conditional addition
   ## The addition is only performed if ctl is "true"
   var t = a
   t += b
   a.ccopy(t, ctl)
 
 func csub*(a: var FF, b: FF, ctl: SecretBool) {.meter.} =
-  ## Constant-time in-place conditional substraction
+  ## Constant-time out-place conditional substraction
   ## The substraction is only performed if ctl is "true"
   var t = a
   t -= b
   a.ccopy(t, ctl)
 
+# ############################################################
+#
+#           Field arithmetic division and inversion
+#
+# ############################################################
+
+func div2*(a: var FF) {.meter.} =
+  ## Modular division by 2
+  ## `a` will be divided in-place
+  #
+  # Normally if `a` is odd we add the modulus before dividing by 2
+  # but this may overflow and we might lose a bit before shifting.
+  # Instead we shift first and then add half the modulus rounded up
+  #
+  # Assuming M is odd, `mp1div2` can be precomputed without
+  # overflowing the "Limbs" by dividing by 2 first
+  # and add 1
+  # Otherwise `mp1div2` should be M/2
+
+  # if a.isOdd:
+  #   a += M
+  # a = a shr 1
+  let wasOdd = a.mres.isOdd()
+  a.mres.shiftRight(1)
+  let carry {.used.} = a.mres.cadd(FF.getPrimePlus1div2(), wasOdd)
+
+  # a < M -> a/2 <= M/2:
+  #   a/2 + M/2 < M if a is odd
+  #   a/2       < M if a is even
+  debug: doAssert not carry.bool
+
+func inv*(r: var FF, a: FF) =
+  ## Inversion modulo p
+  ##
+  ## The inverse of 0 is 0.
+  ## Incidentally this avoids extra check
+  ## to convert Jacobian and Projective coordinates
+  ## to affine for elliptic curve
+  r.mres.invmod(a.mres, FF.getR2modP(), FF.fieldMod())
+
+func inv*(a: var FF) =
+  ## Inversion modulo p
+  ##
+  ## The inverse of 0 is 0.
+  ## Incidentally this avoids extra check
+  ## to convert Jacobian and Projective coordinates
+  ## to affine for elliptic curve
+  a.inv(a)
+
 # ############################################################
 #
 #               Field arithmetic exponentiation

constantine/arithmetic/finite_fields_double_precision.nim

@@ -162,7 +162,7 @@ func neg2xMod*(r: var FpDbl, a: FpDbl) =
       subB(borrow, t.limbs2x[i], M.limbs[i-N], a.limbs2x[i], borrow)
 
     # Zero the result if input was zero
-    t.limbs2x.czero(isZero)
+    t.limbs2x.csetZero(isZero)
     r = t
 
 func prod2xImpl(