Initial support for Twisted Edwards curves (#167)

* Point decoding: optimized sqrt for p ≡ 5 (mod 8) (Curve25519)

* Implement fused sqrt(u/v) for twisted edwards point deserialization

* Introduce twisted edwards affine

* Allow declaration of curve field elements (and fight against recursive dependencies

* Twisted edwards group law + tests

* Add support for jubjub and bandersnatch #162

* test twisted edwards scalar mul

README.md

@@ -89,8 +89,15 @@ Curves:
 - BN254_Snarks (Zero-Knowledge Proofs, Snarks, Starks, Zcash, Ethereum 1)
 - BLS12-377 (Zexe)
 - BLS12-381 (Algorand, Chia Networks, Dfinity, Ethereum 2, Filecoin, Zcash Sapling)
-- BW6-671 (Celo, EY Blockchain) (Pairings are WIP)
+- BW6-671 (Celo, EY Blockchain) (Pairings are WIP)\
+  BLS12-377 is embedded in BW6-761 for one layer proof composition in zk-SNARKS.
 
+### Other curves
+
+- Curve25519, used in ed25519 and X25519 from TLS 1.3 protocol and the Signal protocol.
+  With Ristretto, it can be used in bulletproofs.
+- Jubjub, a curve embedded in BLS12-381 scalar field to be used in zk-SNARKS circuits.
+- Bandersnatch, a more efficient curve embedded in BLS12-381 scalar field to be used in zk-SNARKS circuits.
 ## Security
 
 Hardening an implementation against all existing and upcoming attack vectors is an extremely complex task.

constantine.nimble

@@ -78,6 +78,11 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   ("tests/t_ec_shortw_jac_g1_mul_vs_ref.nim", false),
   ("tests/t_ec_shortw_jac_g1_mixed_add.nim", false),
 
+  ("tests/t_ec_twedwards_prj_add_double", false),
+  ("tests/t_ec_twedwards_prj_mul_sanity", false),
+  ("tests/t_ec_twedwards_prj_mul_distri", false),
+ 
+
   # Elliptic curve arithmetic G2
   # ----------------------------------------------------------
   # ("tests/t_ec_shortw_prj_g2_add_double_bn254_snarks.nim", false),

constantine/arithmetic/README.md

@@ -35,3 +35,21 @@ where code size is not an issue for example for multi-precision addition.
 - Faster big-integer modular multiplication for most moduli\
   Gautam Botrel, Gus Gutoski, and Thomas Piellard, 2020\
   https://hackmd.io/@zkteam/modular_multiplication
+
+### Square roots
+
+- Probabilistic Primality Testing
+  A. Oliver L. Atkin
+  http://algo.inria.fr/seminars/sem91-92/atkin.pdf
+
+- Square root computation over even extension fields
+  Gora Adj, Francisco Rodríguez-Henríquez, 2012
+  https://eprint.iacr.org/2012/685
+
+- A Complete Generalization of Atkin’s Square Root Algorithm
+  Armand Stefan Rotaru, Sorin Iftene, 2013
+  https://profs.info.uaic.ro/~siftene/fi125(1)04.pdf
+
+- Computing Square Roots Faster than the Tonelli-Shanks/Bernstein Algorithm
+  Palash Sarkar, 2020
+  https://eprint.iacr.org/2020/1407

constantine/arithmetic/finite_fields.nim

@@ -28,7 +28,7 @@
 
 import
   ../primitives,
-  ../config/[common, type_ff, curves],
+  ../config/[common, type_ff, curves_prop_field_core, curves_prop_field_derived],
   ./bigints, ./bigints_montgomery
 
 when UseASM_X86_64:
@@ -138,6 +138,13 @@ func setOne*(a: var FF) =
   #       Check if the compiler optimizes it away
   a.mres = FF.getMontyOne()
 
+func setMinusOne*(a: var FF) =
+  ## Set ``a`` to -1 (mod p)
+  # Note: we need -1 in Montgomery residue form
+  # TODO: Nim codegen is not optimal it uses a temporary
+  #       Check if the compiler optimizes it away
+  a.mres = FF.getMontyPrimeMinus1()
+
 func `+=`*(a: var FF, b: FF) {.meter.} =
   ## In-place addition modulo p
   when UseASM_X86_64 and a.mres.limbs.len <= 6: # TODO: handle spilling

constantine/arithmetic/finite_fields_square_root.nim

@@ -10,7 +10,8 @@ import
   ../primitives,
   ../config/[common, type_ff, curves],
   ../curves/zoo_square_roots,
-  ./bigints, ./finite_fields
+  ./bigints, ./finite_fields,
+  ./finite_fields_inversion
 
 # ############################################################
 #
@@ -65,19 +66,95 @@ func invsqrt_p3mod4*(r: var Fp, a: Fp) =
   ## The square root, if it exist is multivalued,
   ## i.e. both x² == (-x)²
   ## This procedure returns a deterministic result
-  # TODO: deterministic sign
-  #
   # Algorithm
   #
   #
-  # From Euler's criterion:   a^((p-1)/2)) ≡ 1 (mod p) if square
+  # From Euler's criterion: 
+  #    𝛘(a) = a^((p-1)/2)) ≡ 1 (mod p) if square
   # a^((p-1)/2)) * a^-1 ≡ 1/a  (mod p)
   # a^((p-3)/2))        ≡ 1/a  (mod p)
   # a^((p-3)/4))        ≡ 1/√a (mod p)      # Requires p ≡ 3 (mod 4)
-  static: doAssert BaseType(Fp.C.Mod.limbs[0]) mod 4 == 3
+  static: doAssert Fp.C.hasP3mod4_primeModulus()
   r = a
   r.powUnsafeExponent(Fp.getPrimeMinus3div4_BE())
 
+# Specialized routine for p ≡ 5 (mod 8)
+# ------------------------------------------------------------
+
+func hasP5mod8_primeModulus(C: static Curve): static bool =
+  ## Returns true iff p ≡ 5 (mod 8)
+  (BaseType(C.Mod.limbs[0]) and 7) == 5
+
+func invsqrt_p5mod8*(r: var Fp, a: Fp) =
+  ## Compute the inverse square root of ``a``
+  ##
+  ## This requires ``a`` to be a square
+  ## and the prime field modulus ``p``: p ≡ 5 (mod 8)
+  ##
+  ## The result is undefined otherwise
+  ##
+  ## The square root, if it exist is multivalued,
+  ## i.e. both x² == (-x)²
+  ## This procedure returns a deterministic result
+  #
+  # Intuition: Branching algorithm, that requires √-1 (mod p) precomputation
+  #
+  # From Euler's criterion:
+  #    𝛘(a) = a^((p-1)/2)) ≡ 1 (mod p) if square
+  # a^((p-1)/4))² ≡ 1 (mod p)
+  # if a is square, a^((p-1)/4)) ≡ ±1 (mod p)
+  #
+  # Case a^((p-1)/4)) ≡ 1 (mod p)
+  #   a^((p-1)/4)) * a⁻¹ ≡  1/a  (mod p)
+  #   a^((p-5)/4))        ≡  1/a  (mod p)
+  #   a^((p-5)/8))        ≡ ±1/√a (mod p)  # Requires p ≡ 5 (mod 8)
+  #
+  # Case a^((p-1)/4)) ≡ -1 (mod p)
+  #   a^((p-1)/4)) * a⁻¹ ≡  -1/a  (mod p)
+  #   a^((p-5)/4))       ≡  -1/a  (mod p)
+  #   a^((p-5)/8))       ≡ ± √-1/√a (mod p)
+  # as p ≡ 5 (mod 8), hence 𝑖 ∈ Fp with 𝑖² ≡ −1 (mod p)
+  #   a^((p-5)/8)) * 𝑖    ≡ ± 1/√a (mod p)
+  #
+  # Atkin Algorithm: branchless, no precomputation
+  #   Atkin, 1992, http://algo.inria.fr/seminars/sem91-92/atkin.pdf
+  #   Gora Adj 2012, https://eprint.iacr.org/2012/685
+  #   Rotaru, 2013, https://profs.info.uaic.ro/~siftene/fi125(1)04.pdf
+  #
+  # We express √a = αa(β − 1) where β² = −1 and 2aα² = β
+  # confirm that        (αa(β − 1))² = α²a²(β²-2β+1) = α²a²β² - 2a²α²β - a²α²
+  # Which simplifies to (αa(β − 1))² = -aβ² = a
+  #
+  # 𝛘(2) = 2^((p-1)/2) ≡ (-1)^((p²-1)/8) (mod p) hence 2 is QR iff p ≡ ±1 (mod 8)
+  # Here p ≡ 5 (mod 8), so 2 is a QNR, hence 2^((p-1)/2) ≡ -1 (mod 8)
+  #
+  # The product of a quadratic non-residue with quadratic residue is a QNR
+  # as 𝛘(QNR*QR) = 𝛘(QNR).𝛘(QR) = -1*1 = -1, hence:
+  #   (2a)^((p-1)/2) ≡ -1 (mod p)
+  #   (2a)^((p-1)/4) ≡ ± √-1 (mod p)
+  #
+  # Hence we set β = (2a)^((p-1)/4)
+  # and α = (β/2a)⁽¹⸍²⁾= (2a)^(((p-1)/4 - 1)/2) = (2a)^((p-5)/8)
+  static: doAssert Fp.C.hasP5mod8_primeModulus()
+  var alpha{.noInit.}, beta{.noInit.}: Fp
+  
+  # α = (2a)^((p-5)/8)
+  alpha.double(a)
+  beta = alpha
+  alpha.powUnsafeExponent(Fp.getPrimeMinus5div8_BE())
+
+  # Note: if r aliases a, for inverse square root we don't use `a` again
+
+  # β = 2aα²
+  r.square(alpha)
+  beta *= r
+  
+  # √a = αa(β − 1), so 1/√a = α(β − 1)
+  r.setOne()
+  beta -= r
+  r.prod(alpha, beta)
+  
+
 # Specialized routines for addchain-based square roots
 # ------------------------------------------------------------
 
@@ -192,6 +269,8 @@ func invsqrt*[C](r: var Fp[C], a: Fp[C]) =
     r.invsqrt_addchain(a)
   elif C.hasP3mod4_primeModulus():
     r.invsqrt_p3mod4(a)
+  elif C.hasP5mod8_primeModulus():
+    r.invsqrt_p5mod8(a)
   else:
     r.invsqrt_tonelli_shanks(a, useAddChain = C.hasTonelliShanksAddchain())
 
@@ -262,4 +341,88 @@ func invsqrt_if_square*[C](r: var Fp[C], a: Fp[C]): SecretBool =
   result = sqrt_invsqrt_if_square(sqrt, r, a)
 
 {.pop.} # inline
+
+# Fused routines
+# ------------------------------------------------------------
+
+func sqrt_ratio_if_square_p5mod8(r: var Fp, u, v: Fp): SecretBool =
+  ## If u/v is a square, compute √(u/v)
+  ## if not, the result is undefined
+  ## 
+  ## Requires p ≡ 5 (mod 8)
+  ## r must not alias u or v
+  ## 
+  ## The square root, if it exist is multivalued,
+  ## i.e. both (u/v)² == (-u/v)²
+  ## This procedure returns a deterministic result
+  ## This procedure is constant-time
+
+  # References:
+  #   - High-Speed High-Security Signature, Bernstein et al, p15 "Fast decompression", https://ed25519.cr.yp.to/ed25519-20110705.pdf
+  #   - IETF Hash-to-Curve: https://github.com/cfrg/draft-irtf-cfrg-hash-to-curve/blob/9939a07/draft-irtf-cfrg-hash-to-curve.md#optimized-sqrt_ratio-for-q--5-mod-8
+  #   - Pasta curves divsqrt: https://github.com/zcash/pasta/blob/f0f7068/squareroottab.sage#L139-L193
+  #
+  # p ≡ 5 (mod 8), hence 𝑖 ∈ Fp with 𝑖² ≡ −1 (mod p)
+  # if α is a square, with β ≡ α^((p+3)/8) (mod p)
+  # - either β² ≡ α (mod p), hence √α ≡ ± β (mod p)
+  # - or β² ≡ -α (mod p), hence √α ≡ ± 𝑖β (mod p)
+  # (see explanation in invsqrt_p5mod8)
+  #
+  # In our fused division and sqrt case we have
+  # β = (u/v)^((p+3)/8)
+  #   = u^((p+3)/8).v^(p−1−(p+3)/8) via Fermat's little theorem
+  #   = u^((p+3)/8).v^((7p−11)/8)
+  #   = u.u^((p-5)/8).v³.v^((7p−35)/8)
+  #   = uv³.u^((p-5)/8).v^(7(p-5)/8)
+  #   = uv³(uv⁷)^((p−5)/8)
+  #
+  # We can check if β² ≡ -α (mod p)
+  # by checking vβ² ≡ -u (mod p), and then multiply by 𝑖
+  # and if it's neither u or -u it wasn't a square.
+  static: doAssert Fp.C.hasP5mod8_primeModulus()
+  var t {.noInit.}: Fp
+  t.square(v)
+  t *= v
+
+  # r = uv³
+  r.prod(u, t)
+
+  # t = (uv⁷)^((p−5)/8)
+  t *= r
+  t *= v
+  t.powUnsafeExponent(Fp.getPrimeMinus5div8_BE())
+
+  # r = β = uv³(uv⁷)^((p−5)/8)
+  r *= t
+
+  # Check candidate square roots
+  t.square(r)
+  t *= v
+  block:
+    result = t == u
+  block:
+    t.neg()
+    let isSol = t == u
+    result = result or isSol
+    t.prod(r, Fp.C.sqrt_minus_one())
+    r.ccopy(t, isSol)
+
+func sqrt_ratio_if_square*(r: var Fp, u, v: Fp): SecretBool {.inline.} =
+  ## If u/v is a square, compute √(u/v)
+  ## if not, the result is undefined
+  ## 
+  ## r must not alias u or v
+  ## 
+  ## The square root, if it exist is multivalued,
+  ## i.e. both (u/v)² == (-u/v)²
+  ## This procedure returns a deterministic result
+  ## This procedure is constant-time
+  when Fp.C.hasP5mod8_primeModulus():
+    sqrt_ratio_if_square_p5mod8(r, u, v)
+  else:
+    # TODO: Fuse inversion and tonelli-shanks and legendre symbol
+    r.inv(v)
+    r *= u
+    result = r.sqrt_if_square()
+
 {.pop.} # raises no exceptions

constantine/config/curves.nim

@@ -6,5 +6,5 @@
 #   * 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 curves_prop_core, curves_prop_derived
+import curves_prop_field_core, curves_prop_field_derived, curves_prop_curve
-export curves_prop_core, curves_prop_derived
+export curves_prop_field_core, curves_prop_field_derived, curves_prop_curve

constantine/config/curves_declaration.nim

@@ -8,9 +8,9 @@
 
 import
   # Internal
-  ./curves_parser
+  ./curves_parser_field
 
-export CurveFamily
+export CurveFamily, SexticTwist
 
 # ############################################################
 #
@@ -108,9 +108,66 @@ declareCurves:
     embedding_degree: 12
     sexticTwist: D_Twist
 
+  curve BabyJubjub: # Curve embedded in BN254_Snarks scalar field
+    # https://iden3-docs.readthedocs.io/en/latest/_downloads/33717d75ab84e11313cc0d8a090b636f/Baby-Jubjub.pdf
+    bitwidth: 254
+    modulus: "0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001"
+
+    # Montgomery form: y² = x³ + 168698x² + x
+    # Edwards form: x² + y² = 1 + dx²y² with d=168696/168700
+    order: "0x60c89ce5c263405370a08b6d0302b0bab3eedb83920ee0a677297dc392126f1"
+    orderBitwidth: 251
+    cofactor: 8
+    # eq_form: Edwards
+    coef_d: "0x1575bd81821016c07a5fd2dee78446612498beee8e01a829736c2b06fb281473"
+
+  curve Jubjub: # Zcash Sapling curve embedded in BLS12-381 scalar field
+    # https://z.cash/technology/jubjub/
+    bitwidth: 255
+    modulus: "0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001"
+    # Montgomery form: y² = x³ + 40962x² + x 
+    # Twisted Edwards: ax² + y² = 1+dx²y² with a = -1 d=-10240/10241
+    order: "0xe7db4ea6533afa906673b0101343b00a6682093ccc81082d0970e5ed6f72cb7"
+    orderBitwidth: 252
+    cofactor: 8
+    eq_form: TwistedEdwards
+    coef_a: -1
+    coef_d: "0x2a9318e74bfa2b48f5fd9207e6bd7fd4292d7f6d37579d2601065fd6d6343eb1"
+
+  curve Bandersnatch: # Anoma curve embedded in BLS12-381 scalar field
+    # https://eprint.iacr.org/2021/1152
+    bitwidth: 255
+    modulus: "0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001"
+
+    # Weierstrass form: y² = x³ − 3763200000x − 7867596800000
+    # Mongomery form: By² = x³ + Ax² + x
+    #   B=0x300c3385d13bedb7c9e229e185c4ce8b1dd3b71366bb97c30855c0aa41d62727
+    #   A=0x4247698f4e32ad45a293959b4ca17afa4a2d2317e4c6ce5023e1f
+    # Twisted Edwards form: −5x² + y² = 1 + dx²y²
+    #   d = 138827208126141220649022263972958607803 / 171449701953573178309673572579671231137
+    order: "0x1cfb69d4ca675f520cce760202687600ff8f87007419047174fd06b52876e7e1"
+    orderBitwidth: 253
+    cofactor: 4
+    eq_form: TwistedEdwards
+    coef_a: -5
+    coef_d: "6389c12633c267cbc66e3bf86be3b6d8cb66677177e54f92b369f2f5188d58e7"
+
   curve Curve25519: # Bernstein curve
     bitwidth: 255
     modulus: "0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed"
+
+    # Montgomery form:            y² = x³ + 486662x² + x 
+    # Edwards form:          x² + y² = 1+dx²y²           with d = 121665/121666
+    # Twisted Edwards form: ax² + y² = 1+dx²y²           with a = 121666 and d = 121665
+    # or for use in Hisil, Wong, Carter, and Dawson extended coordinates 
+    #                       ax² + y² = 1+dx²y²           with a = -1 d=-121665/121666                     
+    order: "0x1000000000000000000000000000000014def9dea2f79cd65812631a5cf5d3ed"
+    orderBItwidth: 253
+    cofactor: 8
+    eq_form: TwistedEdwards
+    coef_a: -1
+    coef_d: "0x52036cee2b6ffe738cc740797779e89800700a4d4141d8ab75eb4dca135978a3"
+
   curve P256: # secp256r1 / NIST P-256
     bitwidth: 256
     modulus: "0xffffffff00000001000000000000000000000000ffffffffffffffffffffffff"

constantine/config/curves_derived.nim

@@ -123,6 +123,13 @@ macro genDerivedConstants*(mode: static DerivedConstantMode): untyped =
         M
       )
     )
+    # const MyCurve_PrimeMinus5div8_BE = primeMinus5div8_BE(MyCurve_Modulus)
+    result.add newConstStmt(
+      used(curve & ff & "_PrimeMinus5div8_BE"), newCall(
+        bindSym"primeMinus5div8_BE",
+        M
+      )
+    )
     # const MyCurve_PrimePlus1div4_BE = primePlus1div4_BE(MyCurve_Modulus)
     result.add newConstStmt(
       used(curve & ff & "_PrimePlus1div4_BE"), newCall(

constantine/config/curves_parser_curve.nim

@@ -0,0 +1,126 @@
+# 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
+  # Standard library
+  std/macros,
+  # Internal
+  ./type_bigint, ./type_ff,
+  ../io/[io_bigints, io_fields],
+  ./curves_declaration, ./curves_parser_field
+
+export CurveFamily, Curve, SexticTwist
+
+# ############################################################
+#
+#                   Curve properties generator
+#
+# ############################################################
+
+template getCoef(c: CurveCoef, curveName: untyped): untyped {.dirty.}=
+  case c.kind
+  of NoCoef:
+    error "Unreachable"
+    nnkDiscardStmt.newTree(newLit "Dummy")
+  of Small:
+    newLit c.coef
+  of Large:
+    newCall(
+      bindSym"fromHex",
+      nnkBracketExpr.newTree(bindSym"Fp", curveName),
+      newLit c.coefHex
+    )
+
+proc genCurveConstants(defs: seq[CurveParams]): NimNode =
+  ## Generate curves main constants
+
+  # MapCurveBitWidth & MapCurveOrderBitWidth
+  # are workaround for https://github.com/nim-lang/Nim/issues/16774
+
+  var MapCurveFamily = nnkBracket.newTree()
+  var curveEllipticStmts = newStmtList()
+
+  for curveDef in defs:
+    curveDef.name.expectKind(nnkIdent)
+    curveDef.bitWidth.expectKind(nnkIntLit)
+    curveDef.modulus.expectKind(nnkStrLit)
+
+    let curve = curveDef.name
+    let family = curveDef.family
+
+    MapCurveFamily.add nnkExprColonExpr.newTree(
+        curve, newLit(family)
+    )
+
+    # Curve equation
+    # -----------------------------------------------
+    curveEllipticStmts.add newConstStmt(
+      exported($curve & "_equation_form"),
+      newLit curveDef.eq_form
+    )
+
+    if curveDef.eq_form == ShortWeierstrass and
+         curveDef.coef_A.kind != NoCoef and curveDef.coef_B.kind != NoCoef:
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_coef_A"),
+        curveDef.coef_A.getCoef(curve)
+      )
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_coef_B"),
+        curveDef.coef_B.getCoef(curve)
+      )
+
+      # Towering
+      # -----------------------------------------------
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_nonresidue_fp"),
+        curveDef.nonresidue_fp
+      )
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_nonresidue_fp2"),
+        curveDef.nonresidue_fp2
+      )
+
+      # Pairing
+      # -----------------------------------------------
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_embedding_degree"),
+        newLit curveDef.embedding_degree
+      )
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_sexticTwist"),
+        newLit curveDef.sexticTwist
+      )
+    
+    if curveDef.eq_form == TwistedEdwards and
+         curveDef.coef_A.kind != NoCoef and curveDef.coef_D.kind != NoCoef:
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_coef_A"),
+        curveDef.coef_A.getCoef(curve)
+      )
+      curveEllipticStmts.add newConstStmt(
+        exported($curve & "_coef_D"),
+        curveDef.coef_D.getCoef(curve)
+      )
+
+  # end for ---------------------------------------------------
+
+  result = newStmtList()
+
+  # const CurveFamily: array[Curve, CurveFamily] = ...
+  result.add newConstStmt(
+    exported("CurveFamilies"), MapCurveFamily
+  )
+
+  result.add curveEllipticStmts
+
+  # echo result.toStrLit()
+
+macro setupCurves(): untyped =
+  result = genCurveConstants(curvesDefinitions)
+
+setupCurves()

constantine/config/curves_parser_field.nim

@@ -43,13 +43,14 @@ type
     Large
 
   CurveCoef* = object
-    case kind: CurveCoefKind
+    case kind*: CurveCoefKind
     of NoCoef: discard
-    of Small: coef: int
+    of Small: coef*: int
-    of Large: coefHex: string
+    of Large: coefHex*: string
 
   CurveEquationForm* = enum
     ShortWeierstrass
+    TwistedEdwards
 
   SexticTwist* = enum
     ## The sextic twist type of the current elliptic curve
@@ -84,36 +85,37 @@ type
     D_Twist
     M_Twist
 
-  CurveParams = object
+  CurveParams* = object
     ## All the curve parameters that may be defined
     # Note: we don't use case object here, the transition is annoying
     #       and would force use to scan all "kind" field (eq_form, family, ...)
     #       before instantiating the object.
-    name: NimNode
+    name*: NimNode
 
     # Field parameters
-    bitWidth: NimNode # nnkIntLit
+    bitWidth*: NimNode # nnkIntLit
-    modulus: NimNode  # nnkStrLit (hex)
+    modulus*: NimNode  # nnkStrLit (hex)
 
     # Towering
-    nonresidue_fp: NimNode # nnkIntLit
+    nonresidue_fp*: NimNode # nnkIntLit
-    nonresidue_fp2: NimNode # nnkPar(nnkIntLit, nnkIntLit)
+    nonresidue_fp2*: NimNode # nnkPar(nnkIntLit, nnkIntLit)
 
     # Curve parameters
-    eq_form: CurveEquationForm
+    eq_form*: CurveEquationForm
-    coef_A: CurveCoef
+    coef_A*: CurveCoef
-    coef_B: CurveCoef
+    coef_B*: CurveCoef
-    order: NimNode # nnkStrLit (hex)
+    coef_D*: CurveCoef
-    orderBitwidth: NimNode # nnkIntLit
+    order*: NimNode # nnkStrLit (hex)
+    orderBitwidth*: NimNode # nnkIntLit
 
-    embedding_degree: int
+    embedding_degree*: int
-    sexticTwist: SexticTwist
+    sexticTwist*: SexticTwist
 
-    family: CurveFamily
+    family*: CurveFamily
 
-var curvesDefinitions {.compileTime.}: seq[CurveParams]
+var curvesDefinitions* {.compileTime.}: seq[CurveParams]
 
-proc parseCurveDecls(defs: var seq[CurveParams], curves: NimNode) =
+proc parseCurveDecls*(defs: var seq[CurveParams], curves: NimNode) =
   ## Parse the curve declarations and store them in the curve definitions
   curves.expectKind(nnkStmtList)
 
@@ -180,6 +182,10 @@ proc parseCurveDecls(defs: var seq[CurveParams], curves: NimNode) =
       elif sectionId.eqIdent"coef_a":
         if sectionVal.kind == nnkIntLit:
           params.coef_A = CurveCoef(kind: Small, coef: sectionVal.intVal.int)
+        elif sectionVal.kind == nnkPrefix: # Got -1
+          sectionVal[0].expectIdent"-"
+          sectionVal[1].expectKind(nnkIntLit)
+          params.coef_A = CurveCoef(kind: Small, coef: -sectionVal[1].intVal.int)
         else:
           params.coef_A = CurveCoef(kind: Large, coefHex: sectionVal.strVal)
       elif sectionId.eqIdent"coef_b":
@@ -191,6 +197,15 @@ proc parseCurveDecls(defs: var seq[CurveParams], curves: NimNode) =
           params.coef_B = CurveCoef(kind: Small, coef: -sectionVal[1].intVal.int)
         else:
           params.coef_B = CurveCoef(kind: Large, coefHex: sectionVal.strVal)
+      elif sectionId.eqIdent"coef_d":
+        if sectionVal.kind == nnkIntLit:
+          params.coef_D = CurveCoef(kind: Small, coef: sectionVal.intVal.int)
+        elif sectionVal.kind == nnkPrefix: # Got -1
+          sectionVal[0].expectIdent"-"
+          sectionVal[1].expectKind(nnkIntLit)
+          params.coef_D = CurveCoef(kind: Small, coef: -sectionVal[1].intVal.int)
+        else:
+          params.coef_D = CurveCoef(kind: Large, coefHex: sectionVal.strVal)
       elif sectionId.eqIdent"order":
         params.order = sectionVal
       elif sectionId.eqIdent"orderBitwidth":
@@ -211,28 +226,14 @@ proc parseCurveDecls(defs: var seq[CurveParams], curves: NimNode) =
 
     defs.add params
 
-proc exported(id: string): NimNode =
+proc exported*(id: string): NimNode =
   nnkPostfix.newTree(
     ident"*",
     ident(id)
   )
 
-template getCoef(c: CurveCoef, width: NimNode): untyped {.dirty.}=
+proc genFieldsConstants(defs: seq[CurveParams]): NimNode =
-  case c.kind
+  ## Generate fields main constants
-  of NoCoef:
-    error "Unreachable"
-    nnkDiscardStmt.newTree(newLit "Dummy")
-  of Small:
-    newLit c.coef
-  of Large:
-    newCall(
-      bindSym"fromHex",
-      nnkBracketExpr.newTree(bindSym"BigInt", width),
-      newLit c.coefHex
-    )
-
-proc genMainConstants(defs: var seq[CurveParams]): NimNode =
-  ## Generate curves and fields main constants
 
   # MapCurveBitWidth & MapCurveOrderBitWidth
   # are workaround for https://github.com/nim-lang/Nim/issues/16774
@@ -240,12 +241,10 @@ proc genMainConstants(defs: var seq[CurveParams]): NimNode =
   var Curves: seq[NimNode]
   var MapCurveBitWidth = nnkBracket.newTree()
   var MapCurveOrderBitWidth = nnkBracket.newTree()
-  var MapCurveFamily = nnkBracket.newTree()
   var curveModStmts = newStmtList()
-  var curveEllipticStmts = newStmtList()
-  var curveExtraStmts = newStmtList()
 
   for curveDef in defs:
+
     curveDef.name.expectKind(nnkIdent)
     curveDef.bitWidth.expectKind(nnkIntLit)
     curveDef.modulus.expectKind(nnkStrLit)
@@ -253,9 +252,10 @@ proc genMainConstants(defs: var seq[CurveParams]): NimNode =
     let curve = curveDef.name
     let bitWidth = curveDef.bitWidth
     let modulus = curveDef.modulus
-    let family = curveDef.family
 
     Curves.add curve
+
+    # Field Fp
     # "BN254_Snarks: 254" array construction expression
     MapCurveBitWidth.add nnkExprColonExpr.newTree(
       curve, bitWidth
@@ -271,19 +271,10 @@ proc genMainConstants(defs: var seq[CurveParams]): NimNode =
       )
     )
 
-    MapCurveFamily.add nnkExprColonExpr.newTree(
+    # Field Fr
-        curve, newLit(family)
-    )
-
-    # Curve equation
-    # -----------------------------------------------
-    curveEllipticStmts.add newConstStmt(
-      exported($curve & "_equation_form"),
-      newLit curveDef.eq_form
-    )
     if not curveDef.order.isNil:
       curveDef.orderBitwidth.expectKind(nnkIntLit)
-      curveEllipticStmts.add newConstStmt(
+      curveModStmts.add newConstStmt(
         exported($curve & "_Order"),
         newCall(
           bindSym"fromHex",
@@ -295,7 +286,7 @@ proc genMainConstants(defs: var seq[CurveParams]): NimNode =
         curve, curveDef.orderBitwidth
       )
     else: # Dummy
-      curveEllipticStmts.add newConstStmt(
+      curveModStmts.add newConstStmt(
         exported($curve & "_Order"),
         newCall(
           bindSym"fromHex",
@@ -307,38 +298,6 @@ proc genMainConstants(defs: var seq[CurveParams]): NimNode =
         curve, newLit 1
       )
 
-    if curveDef.coef_A.kind != NoCoef and curveDef.coef_B.kind != NoCoef:
-      curveEllipticStmts.add newConstStmt(
-        exported($curve & "_coef_A"),
-        curveDef.coef_A.getCoef(bitWidth)
-      )
-      curveEllipticStmts.add newConstStmt(
-        exported($curve & "_coef_B"),
-        curveDef.coef_B.getCoef(bitWidth)
-      )
-
-      # Towering
-      # -----------------------------------------------
-      curveEllipticStmts.add newConstStmt(
-        exported($curve & "_nonresidue_fp"),
-        curveDef.nonresidue_fp
-      )
-      curveEllipticStmts.add newConstStmt(
-        exported($curve & "_nonresidue_fp2"),
-        curveDef.nonresidue_fp2
-      )
-
-      # Pairing
-      # -----------------------------------------------
-      curveEllipticStmts.add newConstStmt(
-        exported($curve & "_embedding_degree"),
-        newLit curveDef.embedding_degree
-      )
-      curveEllipticStmts.add newConstStmt(
-        exported($curve & "_sexticTwist"),
-        newLit curveDef.sexticTwist
-      )
-
   # end for ---------------------------------------------------
 
   result = newStmtList()
@@ -356,19 +315,12 @@ proc genMainConstants(defs: var seq[CurveParams]): NimNode =
   result.add newConstStmt(
     exported("CurveBitWidth"), MapCurveBitWidth
   )
-  # const CurveFamily: array[Curve, CurveFamily] = ...
+  result.add curveModStmts
-  result.add newConstStmt(
-    exported("CurveFamilies"), MapCurveFamily
-  )
   # const CurveOrderBitSize: array[Curve, int] = ...
   result.add newConstStmt(
     exported("CurveOrderBitWidth"), MapCurveOrderBitWidth
   )
 
-  result.add curveModStmts
-  result.add curveEllipticStmts
-  result.add curveExtraStmts
-
   # echo result.toStrLit()
 
 macro declareCurves*(curves: untyped): untyped =
@@ -389,4 +341,4 @@ macro declareCurves*(curves: untyped): untyped =
 
   curves.expectKind(nnkStmtList)
   curvesDefinitions.parseCurveDecls(curves)
-  result = curvesDefinitions.genMainConstants()
+  result = curvesDefinitions.genFieldsConstants()

constantine/config/curves_prop_curve.nim

@@ -1,4 +1,3 @@
-# Constantine
 # Copyright (c) 2018-2019    Status Research & Development GmbH
 # Copyright (c) 2020-Present Mamy André-Ratsimbazafy
 # Licensed and distributed under either of
@@ -10,44 +9,19 @@ import
   # Standard library
   std/macros,
   # Internal
-  ./type_bigint, ./common,
+  ./type_bigint,
-  ./curves_declaration, ./curves_parser
+  ./curves_declaration, ./curves_parser_curve
 
 export CurveFamily, Curve, SexticTwist
 
-# ############################################################
-#
-#                 Field properties
-#
-# ############################################################
-
-{.experimental: "dynamicBindSym".}
-
-macro Mod*(C: static Curve): untyped =
-  ## Get the Modulus associated to a curve
-  result = bindSym($C & "_Modulus")
-
-template getCurveBitwidth*(C: Curve): int =
-  ## Returns the number of bits taken by the curve modulus
-  CurveBitWidth[C]
-
-template matchingBigInt*(C: static Curve): untyped =
-  # Workaround: https://github.com/nim-lang/Nim/issues/16774
-  BigInt[CurveBitWidth[C]]
-
-template family*(C: Curve): CurveFamily =
-  CurveFamilies[C]
-
-template matchingLimbs2x*(C: Curve): untyped =
-  const N2 = wordsRequired(getCurveBitwidth(C)) * 2 # TODO upstream, not precomputing N2 breaks semcheck
-  array[N2, SecretWord] # TODO upstream, using Limbs[N2] breaks semcheck
-
 # ############################################################
 #
 #                   Curve properties
 #
 # ############################################################
 
+{.experimental: "dynamicBindSym".}
+
 macro getCurveOrder*(C: static Curve): untyped =
   ## Get the curve order `r`
   ## i.e. the number of points on the elliptic curve
@@ -65,6 +39,9 @@ template matchingOrderBigInt*(C: static Curve): untyped =
   # Workaround: https://github.com/nim-lang/Nim/issues/16774
   BigInt[CurveOrderBitWidth[C]]
 
+template family*(C: Curve): CurveFamily =
+  CurveFamilies[C]
+
 macro getEquationForm*(C: static Curve): untyped =
   ## Returns the equation form
   ## (ShortWeierstrass, Montgomery, Twisted Edwards, Weierstrass, ...)
@@ -82,6 +59,12 @@ macro getCoefB*(C: static Curve): untyped =
   ## or a bigInt depending on the curve
   result = bindSym($C & "_coef_B")
 
+macro getCoefD*(C: static Curve): untyped =
+  ## Returns the D coefficient of the curve
+  ## The return type is polymorphic, it can be an int
+  ## or a bigInt depending on the curve
+  result = bindSym($C & "_coef_D")
+
 macro getNonResidueFp*(C: static Curve): untyped =
   ## Returns the tower extension (and twist) non-residue for 𝔽p
   ## Depending on the curve it might be:

constantine/config/curves_prop_field_core.nim

@@ -0,0 +1,39 @@
+# 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
+  # Standard library
+  std/macros,
+  # Internal
+  ./type_bigint, ./common,
+  ./curves_declaration
+
+export Curve
+
+# ############################################################
+#
+#                 Field properties
+#
+# ############################################################
+
+{.experimental: "dynamicBindSym".}
+
+macro Mod*(C: static Curve): untyped =
+  ## Get the Modulus associated to a curve
+  result = bindSym($C & "_Modulus")
+
+template getCurveBitwidth*(C: Curve): int =
+  ## Returns the number of bits taken by the curve modulus
+  CurveBitWidth[C]
+
+template matchingBigInt*(C: static Curve): untyped =
+  # Workaround: https://github.com/nim-lang/Nim/issues/16774
+  BigInt[CurveBitWidth[C]]
+
+template matchingLimbs2x*(C: Curve): untyped =
+  const N2 = wordsRequired(getCurveBitwidth(C)) * 2 # TODO upstream, not precomputing N2 breaks semcheck
+  array[N2, SecretWord] # TODO upstream, using Limbs[N2] breaks semcheck

constantine/config/curves_prop_field_derived.nim

@@ -11,7 +11,7 @@ import
   std/macros,
   # Internal
   ./type_bigint, ./type_ff, ./common,
-  ./curves_declaration, ./curves_prop_core, ./curves_derived
+  ./curves_declaration, ./curves_prop_field_core, ./curves_derived
 
 # ############################################################
 #
@@ -108,9 +108,13 @@ macro getPrimeMinus1div2_BE*(ff: type FF): untyped =
   result = bindConstant(ff, "PrimeMinus1div2_BE")
 
 macro getPrimeMinus3div4_BE*(ff: type FF): untyped =
-  ## Get (P-3) / 2 in big-endian serialized format
+  ## Get (P-3) / 4 in big-endian serialized format
   result = bindConstant(ff, "PrimeMinus3div4_BE")
 
+macro getPrimeMinus5div8_BE*(ff: type FF): untyped =
+  ## Get (P-5) / 8 in big-endian serialized format
+  result = bindConstant(ff, "PrimeMinus5div8_BE")
+
 macro getPrimePlus1div4_BE*(ff: type FF): untyped =
   ## Get (P+1) / 4 for an odd prime in big-endian serialized format
   result = bindConstant(ff, "PrimePlus1div4_BE")

constantine/config/precompute.nim

@@ -236,7 +236,9 @@ func checkValidModulus(M: BigInt) =
   const expectedMsb = M.bits-1 - WordBitWidth * (M.limbs.len - 1)
   let msb = log2(BaseType(M.limbs[^1]))
 
-  doAssert msb == expectedMsb, "Internal Error: the modulus must use all declared bits and only those"
+  doAssert msb == expectedMsb, "Internal Error: the modulus must use all declared bits and only those:\n" &
+    "    Modulus '" & M.toHex() & "' is declared with " & $M.bits &
+    " bits but uses " & $(msb + WordBitWidth * (M.limbs.len - 1)) & " bits."
 
 func countSpareBits*(M: BigInt): int =
   ## Count the number of extra bits
@@ -434,6 +436,25 @@ func primeMinus3div4_BE*[bits: static int](
 
   result.exportRawUint(tmp, bigEndian)
 
+func primeMinus5div8_BE*[bits: static int](
+       P: BigInt[bits]
+     ): array[(bits+7) div 8, byte] {.noInit.} =
+  ## For an input prime `p`, compute (p-5)/8
+  ## and return the result as a canonical byte array / octet string
+  ## For use to check if a number is a square (quadratic residue)
+  ## and if so compute the square root in a fused manner
+  ##
+  # Output size:
+  # - (bits + 7) div 8: bits => byte conversion rounded up
+  # - (bits + 7 - 3): dividing by 8 means 3 bits is unused
+  # => TODO: reduce the output size (to potentially save a byte and corresponding multiplication/squarings)
+
+  var tmp = P
+  discard tmp.sub(5)
+  tmp.shiftRight(3)
+
+  result.exportRawUint(tmp, bigEndian)
+
 func primePlus1Div4_BE*[bits: static int](
        P: BigInt[bits]
      ): array[(bits+7) div 8, byte] {.noInit.} =

constantine/config/type_ff.nim

@@ -9,7 +9,7 @@
 import
   ./common,
   ./curves_declaration,
-  ./curves_prop_core
+  ./curves_prop_field_core
 
 export matchingBigInt
 

constantine/curves/bandersnatch_sqrt.nim

@@ -0,0 +1,19 @@
+# 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/[curves, type_bigint, type_ff],
+  ../io/[io_bigints, io_fields],
+  ../arithmetic/finite_fields
+
+const
+  # with e = 2adicity
+  # p == s * 2^e + 1
+  # root_of_unity = smallest_quadratic_nonresidue^s
+  # exponent = (p-1-2^e)/2^e / 2
+  Bandersnatch_TonelliShanks_exponent* = BigInt[222].fromHex"0x39f6d3a994cebea4199cec0404d0ec02a9ded2017fff2dff7fffffff"
+  Bandersnatch_TonelliShanks_twoAdicity* = 32
+  Bandersnatch_TonelliShanks_root_of_unity* = Fp[Bandersnatch].fromHex"0x212d79e5b416b6f0fd56dc8d168d6c0c4024ff270b3e0941b788f500b912f1f"

constantine/curves/curve25519_sqrt.nim

@@ -0,0 +1,27 @@
+# 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/[curves, type_bigint, type_ff],
+  ../io/[io_bigints, io_fields],
+  ../arithmetic/finite_fields
+
+# p ≡ 5 (mod 8), hence 𝑖 ∈ Fp with 𝑖² ≡ −1 (mod p)
+# Hence if α is a square
+# with β ≡ α^((p+3)/8) (mod p)
+# - either β² ≡ α (mod p), hence √α ≡ ±β (mod p)
+# - or β² ≡ -α (mod p), hence √α ≡ ±𝑖β (mod p)
+
+# Sage:
+#   p = Integer('0x7fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffed')
+#   Fp = GF(p)
+#   sqrt_minus1 = Fp(-1).sqrt()
+#   print(Integer(sqrt_minus1).hex())
+const Curve25519_sqrt_minus_one* = Fp[Curve25519].fromHex(
+    "0x2b8324804fc1df0b2b4d00993dfbd7a72f431806ad2fe478c4ee1b274a0ea0b0"
+)

constantine/curves/jubjub_sqrt.nim

@@ -0,0 +1,19 @@
+# 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/[curves, type_bigint, type_ff],
+  ../io/[io_bigints, io_fields],
+  ../arithmetic/finite_fields
+
+const
+  # with e = 2adicity
+  # p == s * 2^e + 1
+  # root_of_unity = smallest_quadratic_nonresidue^s
+  # exponent = (p-1-2^e)/2^e / 2
+  Jubjub_TonelliShanks_exponent* = BigInt[222].fromHex"0x39f6d3a994cebea4199cec0404d0ec02a9ded2017fff2dff7fffffff"
+  Jubjub_TonelliShanks_twoAdicity* = 32
+  Jubjub_TonelliShanks_root_of_unity* = Fp[Jubjub].fromHex"0x212d79e5b416b6f0fd56dc8d168d6c0c4024ff270b3e0941b788f500b912f1f"

constantine/curves/zoo_square_roots.nim

@@ -13,14 +13,18 @@ import
   ./bls12_381_sqrt,
   ./bn254_nogami_sqrt,
   ./bn254_snarks_sqrt,
-  ./bw6_761_sqrt
+  ./bw6_761_sqrt,
+  ./curve25519_sqrt,
+  ./jubjub_sqrt,
+  ./bandersnatch_sqrt
 
 export
   bls12_377_sqrt,
   bls12_381_sqrt,
   bn254_nogami_sqrt,
   bn254_snarks_sqrt,
-  bw6_761_sqrt
+  bw6_761_sqrt,
+  curve25519_sqrt
 
 func hasSqrtAddchain*(C: static Curve): static bool =
   when C in {BLS12_381, BN254_Nogami, BN254_Snarks, BW6_761}:
@@ -39,3 +43,7 @@ func hasTonelliShanksAddchain*(C: static Curve): static bool =
     true
   else:
     false
+
+macro sqrt_minus_one*(C: static Curve): untyped = 
+  ## Return 𝑖 ∈ Fp with 𝑖² ≡ −1 (mod p)
+  return bindSym($C & "_sqrt_minus_one")

constantine/elliptic/ec_shortweierstrass_affine.nim

@@ -17,7 +17,7 @@ import
 # ############################################################
 #
 #             Elliptic Curve in Short Weierstrass form
-#                 with Projective Coordinates
+#                 with Affine Coordinates
 #
 # ############################################################
 
@@ -57,10 +57,10 @@ func curve_eq_rhs*[F](y2: var F, x: F, Tw: static Twisted) =
 
   when Tw == NotOnTwist:
     when F.C.getCoefB() >= 0:
-      y2.fromInt F.C.getCoefB()
+      y2.fromUint uint F.C.getCoefB()
       y2 += t
     else:
-      y2.fromInt -F.C.getCoefB()
+      y2.fromUint uint -F.C.getCoefB()
       y2.diff(t, y2)
   else:
     y2.sum(F.C.getCoefB_G2, t)
@@ -86,7 +86,6 @@ func trySetFromCoordX*[F, Tw](
   ## Try to create a point the elliptic curve
   ## y² = x³ + a x + b     (affine coordinate)
   ##
-  ## The `Z` coordinates is set to 1
   ##
   ## return true and update `P` if `x` leads to a valid point
   ## return false otherwise, in that case `P` is undefined.

constantine/elliptic/ec_shortweierstrass_projective.nim

@@ -157,6 +157,7 @@ func sum*[F; Tw: static Twisted](
   ##   y² = x³ + a x + b
   ##
   ## ``r`` is initialized/overwritten with the sum
+  ## ``r`` may alias P
   ##
   ## Implementation is constant-time, in particular it will not expose
   ## that P == Q or P == -Q or P or Q are the infinity points
@@ -252,6 +253,8 @@ func madd*[F; Tw: static Twisted](
   ## with p in Projective coordinates and Q in affine coordinates
   ##
   ##   R = P + Q
+  ## 
+  ## ``r`` may alias P
 
   # TODO: static doAssert odd order
 
@@ -321,6 +324,7 @@ func double*[F; Tw: static Twisted](
   ##   y² = x³ + a x + b
   ##
   ## ``r`` is initialized/overwritten with the sum
+  ## ``r`` may alias P
   ##
   ## Implementation is constant-time, in particular it will not expose
   ## that `P` is an infinity point.

constantine/elliptic/ec_twistededwards_affine.nim

@@ -0,0 +1,149 @@
+# 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
+  ../primitives,
+  ../config/[common, curves],
+  ../arithmetic,
+  ../towers,
+  ../io/[io_fields, io_towers]
+
+# ############################################################
+#
+#             Elliptic Curve in Twisted Edwards form
+#                 with Affine Coordinates
+#
+# ############################################################
+
+type ECP_TwEdwards_Aff*[F] = object
+  ## Elliptic curve point for a curve in Twisted Edwards form
+  ##   ax²+y²=1+dx²y²
+  ## with a, d ≠ 0 and a ≠ d
+  ##
+  ## over a field F
+  x*, y*: F
+
+func `==`*(P, Q: ECP_TwEdwards_Aff): SecretBool =
+  ## Constant-time equality check
+  result = P.x == Q.x
+  result = result and (P.y == Q.y)
+
+
+func isInf*(P: ECP_TwEdwards_Aff): SecretBool =
+  ## Returns true if P is an infinity point
+  ## and false otherwise
+  result = P.x.isZero() and P.y.isOne()
+
+
+func isOnCurve*[F](x, y: F): SecretBool =
+  ## Returns true if the (x, y) coordinates
+  ## represents a point of the Twisted Edwards elliptic curve
+  ## with equation ax²+y²=1+dx²y²
+  var t0{.noInit.}, t1{.noInit.}, t2{.noInit.}: F
+  t0.square(x)
+  t1.square(y)
+  
+  # ax²+y²
+  when F.C.getCoefA() is int:
+    when F.C.getCoefA() == -1:
+      t2.diff(t1, t0)
+    else:
+      t2.prod(t0, F.C.getCoefA())
+      t2 += t1
+  else:
+    t2.prod(F.C.getCoefA(), t0)
+    t2 += t1
+
+  # dx²y²
+  t0 *= t1
+  when F.C.getCoefD() is int:
+    when F.C.getCoefD >= 0:
+      t1.fromUint uint F.C.getCoefD()
+      t0 *= t1
+    else:
+      t1.fromUint uint F.C.getCoefD()
+      t0 *= t1
+      t0.neg()
+  else:
+    t0 *= F.C.getCoefD()
+
+  # ax²+y² - dx²y² =? 1
+  t2 -= t0
+  return t2.isOne()
+
+func trySetFromCoordY*[F](P: var ECP_TwEdwards_Aff[F], y: F): SecretBool =
+  ## Try to create a point the elliptic curve
+  ##   ax²+y²=1+dx²y²    (affine coordinate)
+  ##
+  ##
+  ## return true and update `P` if `y` leads to a valid point
+  ## return false otherwise, in that case `P` is undefined.
+  ##
+  ## Note: Dedicated robust procedures for hashing-to-curve
+  ##       will be provided, this is intended for testing purposes.
+  ## 
+  ##       For **test case generation only**,
+  ##       this is preferred to generating random point
+  ##       via random scalar multiplication of the curve generator
+  ##       as the latter assumes:
+  ##       - point addition, doubling work
+  ##       - scalar multiplication works
+  ##       - a generator point is defined
+  ##       i.e. you can't test unless everything is already working
+
+  # https://eprint.iacr.org/2015/677.pdf
+  # p2: Encoding and parsing curve points.
+  # x² = (y² − 1)/(dy² − a)
+  var t {.noInit.}: F
+  t.square(y)
+
+  # (dy² − a)
+  when F.C.getCoefD() is int:
+    when F.C.getCoefD() >= 0:
+      P.y.fromUint uint F.C.getCoefD()
+    else:
+      P.y.fromUint uint -F.C.getCoefD()
+      P.y.neg()
+  else:
+    P.y = F.C.getCoefD()
+  P.y *= t
+  when F.C.getCoefA() is int:
+    when F.C.getCoefA == -1:
+      P.x.setOne()
+      P.y += P.x
+    elif F.C.getCoefA >= 0:
+      P.x.fromUint uint F.C.getCoefA()
+      P.y -= P.x
+    else:
+      P.x.fromUint uint -F.C.getCoefA()
+      P.y += P.x
+  else:
+    P.y -= F.C.getCoefA()
+
+  # y² − 1
+  P.x.setMinusOne()
+  P.x += t
+
+  # √((y² − 1)/(dy² − a))
+  result = sqrt_ratio_if_square(t, P.x, P.y)
+  P.x = t
+  P.y = y
+
+func neg*(P: var ECP_TwEdwards_Aff, Q: ECP_TwEdwards_Aff) =
+  ## Negate ``P``
+  P.x.neg(Q.x)
+  P.y = Q.y
+
+func neg*(P: var ECP_TwEdwards_Aff) =
+  ## Negate ``P``
+  P.x.neg()
+
+func cneg*(P: var ECP_TwEdwards_Aff, ctl: CTBool) =
+  ## Conditional negation.
+  ## Negate if ``ctl`` is true
+  P.x.cneg(ctl)

constantine/elliptic/ec_twistededwards_projective.nim

@@ -0,0 +1,306 @@
+# 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
+  ../primitives,
+  ../config/[common, curves],
+  ../arithmetic,
+  ../towers,
+  ./ec_twistededwards_affine
+
+# ############################################################
+#
+#             Elliptic Curve in Twisted Edwards form
+#                 with Projective Coordinates
+#
+# ############################################################
+
+type ECP_TwEdwards_Prj*[F] = object
+  ## Elliptic curve point for a curve in Twisted Edwards form
+  ##   ax²+y²=1+dx²y²
+  ## with a, d ≠ 0 and a ≠ d
+  ##
+  ## over a field F
+  ## 
+  ## in projective coordinate (X, Y, Z)
+  ## with x = X/Z and y = Y/Z
+  ## hence (aX² + Y²)Z² = Z⁴ + dX²Y²
+  x*, y*, z*: F
+
+func `==`*(P, Q: ECP_TwEdwards_Prj): SecretBool =
+  ## Constant-time equality check
+  ## This is a costly operation
+  # Reminder: the representation is not unique
+  var a{.noInit.}, b{.noInit.}: ECP_TwEdwards_Prj.F
+
+  a.prod(P.x, Q.z)
+  b.prod(Q.x, P.z)
+  result = a == b
+
+  a.prod(P.y, Q.z)
+  b.prod(Q.y, P.z)
+  result = result and a == b
+
+func isInf*(P: ECP_TwEdwards_Prj): SecretBool {.inline.} =
+  ## Returns true if P is an infinity point
+  ## and false otherwise
+  result = P.x.isZero() and (P.y == P.z)
+
+func setInf*(P: var ECP_TwEdwards_Prj) {.inline.} =
+  ## Set ``P`` to infinity
+  P.x.setZero()
+  P.y.setOne()
+  P.z.setOne()
+
+func ccopy*(P: var ECP_TwEdwards_Prj, Q: ECP_TwEdwards_Prj, ctl: SecretBool) {.inline.} =
+  ## Constant-time conditional copy
+  ## If ctl is true: Q is copied into P
+  ## if ctl is false: Q is not copied and P is unmodified
+  ## Time and memory accesses are the same whether a copy occurs or not
+  for fP, fQ in fields(P, Q):
+    ccopy(fP, fQ, ctl)
+
+func trySetFromCoordY*[F](
+       P: var ECP_TwEdwards_Prj[F],
+       y: F): SecretBool =
+  ## Try to create a point the elliptic curve
+  ##   ax²+y²=1+dx²y²     (affine coordinate)
+  ##
+  ## The `Z` coordinates is set to 1
+  ##
+  ## return true and update `P` if `y` leads to a valid point
+  ## return false otherwise, in that case `P` is undefined.
+  ##
+  ## Note: Dedicated robust procedures for hashing-to-curve
+  ##       will be provided, this is intended for testing purposes.
+  ## 
+  ##       For **test case generation only**,
+  ##       this is preferred to generating random point
+  ##       via random scalar multiplication of the curve generator
+  ##       as the latter assumes:
+  ##       - point addition, doubling work
+  ##       - scalar multiplication works
+  ##       - a generator point is defined
+  ##       i.e. you can't test unless everything is already working
+  
+  var Q{.noInit.}: ECP_TwEdwards_Aff[F]
+  result = Q.trySetFromCoordY(y)
+
+  P.x = Q.x
+  P.y = Q.y
+  P.z.setOne()
+
+func trySetFromCoordsYandZ*[F](
+       P: var ECP_TwEdwards_Prj[F],
+       y, z: F): SecretBool =
+  ## Try to create a point the elliptic curve
+  ##   ax²+y²=1+dx²y²     (affine coordinate)
+  ##
+  ## return true and update `P` if `y` leads to a valid point
+  ## return false otherwise, in that case `P` is undefined.
+  ##
+  ## Note: Dedicated robust procedures for hashing-to-curve
+  ##       will be provided, this is intended for testing purposes.
+  ## 
+  ##       For **test case generation only**,
+  ##       this is preferred to generating random point
+  ##       via random scalar multiplication of the curve generator
+  ##       as the latter assumes:
+  ##       - point addition, doubling work
+  ##       - scalar multiplication works
+  ##       - a generator point is defined
+  ##       i.e. you can't test unless everything is already working
+  
+  var Q{.noInit.}: ECP_TwEdwards_Aff[F]
+  result = Q.trySetFromCoordY(y)
+
+  P.x.prod(Q.x, z)
+  P.y.prod(Q.y, z)
+  P.z = z
+
+func neg*(P: var ECP_TwEdwards_Prj, Q: ECP_TwEdwards_Prj) {.inline.} =
+  ## Negate ``P``
+  P.x.neg(Q.x)
+  P.y = Q.y
+  P.z = Q.z
+
+func neg*(P: var ECP_TwEdwards_Prj) {.inline.} =
+  ## Negate ``P``
+  P.x.neg()
+
+func cneg*(P: var ECP_TwEdwards_Prj, ctl: CTBool) {.inline.} =
+  ## Conditional negation.
+  ## Negate if ``ctl`` is true
+  P.x.cneg(ctl)
+
+func sum*[Field](
+       r: var ECP_TwEdwards_Prj[Field],
+       P, Q: ECP_TwEdwards_Prj[Field]
+     ) =
+  ## Elliptic curve point addition for Twisted Edwards curves in projective coordinates
+  ##
+  ##   R = P + Q
+  ##
+  ## Twisted Edwards curves have the following equation in projective coordinates
+  ##   (aX² + Y²)Z² = Z⁴ + dX²Y²
+  ## from the affine equation
+  ##   ax²+y²=1+dx²y²
+  ##
+  ## ``r`` is initialized/overwritten with the sum
+  ## ``r`` may alias P
+  ##
+  ## Implementation is constant-time, in particular it will not expose
+  ## that P == Q or P == -Q or P or Q are the infinity points
+  ## to simple side-channel attacks (SCA)
+  ## This is done by using a "complete" or "exception-free" addition law.
+  #
+  # https://www.hyperelliptic.org/EFD/g1p/auto-twisted-projective.html#addition-add-2008-bbjlp
+  # Cost: 10M + 1S + 1*a + 1*d + 7add.
+  #   A = Z1*Z2
+  #   B = A²
+  #   C = X1*X2
+  #   D = Y1*Y2
+  #   E = d*C*D
+  #   F = B-E
+  #   G = B+E
+  #   X3 = A*F*((X1+Y1)*(X2+Y2)-C-D)
+  #   Y3 = A*G*(D-a*C)
+  #   Z3 = F*G
+  var
+    A{.noInit.}, B{.noInit.}, C{.noInit.}: Field
+    D{.noInit.}, E{.noInit.}, F{.noInit.}: Field
+    G{.noInit.}: Field
+
+  A.prod(P.z, Q.z)
+  B.square(A)
+  C.prod(P.x, Q.x)
+  D.prod(P.y, Q.y)
+  E.prod(C, D)
+  when Field.C.getCoefD() is int:
+    # conversion at compile-time
+    const coefD = block:
+      var d: Field
+      d.fromInt F.C.getCoefD() 
+      d
+    E *= coefD
+  else:
+    E *= Field.C.getCoefD()
+  F.diff(B, E)
+  G.sum(B, E)
+
+  # Aliasing: B and E are unused
+  # We store (P.x+P.y)*(Q.x+Q.y)
+  # so that using r.x or r.y is safe even in case of aliasing
+
+  B.sum(P.x, P.y)
+  E.sum(Q.x, Q.y)
+  B *= E          # B = (X1+Y1)*(X2+Y2)
+  E.sum(C, D)     # E = C+D
+  
+  # Y3 = A*G*(D-a*C)
+  when Field.C.getCoefA() == -1:
+    r.y = E       # (D-a*C) = D+C
+  else:
+    r.y.prod(C, Field.C.getCoefA())
+    r.y.diff(D, r.y)
+  r.y *= A
+  r.y *= G
+
+  # X3 = A*F*((X1+Y1)*(X2+Y2)-C-D)
+  B -= E
+  r.x.prod(A, F)
+  r.x *= B
+
+  # Z3 = F*G
+  r.z.prod(F, G)
+
+func double*[Field](
+       r: var ECP_TwEdwards_Prj[Field],
+       P: ECP_TwEdwards_Prj[Field]
+     ) =
+  ## Elliptic curve point doubling for Twisted Edwards curves in projective coordinates
+  ##
+  ##   R = [2] P
+  ##
+  ## Twisted Edwards curves have the following equation in projective coordinates
+  ##   (aX² + Y²)Z² = Z⁴ + dX²Y²
+  ## from the affine equation
+  ##   ax²+y²=1+dx²y²
+  ##
+  ## ``r`` is initialized/overwritten with the sum
+  ## ``r`` may alias P
+  ##
+  ## Implementation is constant-time, in particular it will not expose
+  ## that P == Q or P == -Q or P or Q are the infinity points
+  ## to simple side-channel attacks (SCA)
+  ## This is done by using a "complete" or "exception-free" addition law.
+  #
+  # https://www.hyperelliptic.org/EFD/g1p/auto-twisted-projective.html#addition-add-2008-bbjlp
+  # Cost: 3M + 4S + 1*a + 6add + 1*2. 
+  #  B = (X1+Y1)²
+  #  C = X1²
+  #  D = Y1²
+  #  E = a*C
+  #  F = E+D
+  #  H = Z1²
+  #  J = F-2*H
+  #  X3 = (B-C-D)*J
+  #  Y3 = F*(E-D)
+  #  Z3 = F*J
+
+  var
+    D{.noInit.}, E{.noInit.}: Field
+    H{.noInit.}, J{.noInit.}: FIeld
+  
+  # (B-C-D) => 2X1Y1, but With squaring and 2 substractions instead of mul + addition
+  # In practice, squaring is not cheap enough to compasate the extra substraction cost.
+  r.x.prod(P.x, P.y)
+  r.x.double()
+
+  D.square(P.y)
+  E.square(P.x)
+  E *= Field.C.getCoefA()
+
+  r.y.sum(E, D)    # Ry stores F = E+D
+  H.square(P.z)
+  H.double()
+  J.diff(r.y, H)   # J = F-2H
+
+  r.x *= J         # X3 = (B-C-D)*J
+  r.z.prod(r.y, J) # Z3 = F*J
+  E -= D           # C stores E-D
+  r.y *= E
+
+func double*(P: var ECP_TwEdwards_Prj) {.inline.} =
+  ## In-place EC doubling
+  P.double(P)
+
+func diff*(r: var ECP_TwEdwards_Prj,
+              P, Q: ECP_TwEdwards_Prj
+     ) {.inline.} =
+  ## r = P - Q
+  ## Can handle r and Q aliasing
+  var nQ {.noInit.}: typeof(Q)
+  nQ.neg(Q)
+  r.sum(P, nQ)
+
+func affineFromProjective*[F](
+       aff: var ECP_TwEdwards_Prj[F],
+       proj: ECP_TwEdwards_Prj[F]) =
+  var invZ {.noInit.}: F
+  invZ.inv(proj.z)
+
+  aff.x.prod(proj.x, invZ)
+  aff.y.prod(proj.y, invZ)
+
+func projectiveFromAffine*[F](
+       proj: var ECP_TwEdwards_Prj[F],
+       aff: ECP_TwEdwards_Prj[F]) {.inline.} =
+  proj.x = aff.x
+  proj.y = aff.y
+  proj.z.setOne()

constantine/io/io_fields.nim

@@ -8,7 +8,7 @@
 
 import
   ./io_bigints,
-  ../config/[common, curves],
+  ../config/[common, type_ff],
   ../arithmetic/finite_fields,
   ../primitives
 
@@ -45,6 +45,7 @@ func fromInt*(dst: var FF,
     let src = isNeg.mux(SecretWord -src, SecretWord src)
     let raw {.noinit.} = (type dst.mres).fromRawUint(cast[array[sizeof(src), byte]](src), cpuEndian)
     dst.fromBig(raw)
+    dst.cneg(isNeg)
 
 func exportRawUint*(dst: var openarray[byte],
                        src: FF,

helpers/prng_unsafe.nim

@@ -13,7 +13,9 @@ import
   ../constantine/elliptic/[
     ec_shortweierstrass_affine,
     ec_shortweierstrass_projective,
-    ec_shortweierstrass_jacobian],
+    ec_shortweierstrass_jacobian,
+    ec_twistededwards_affine,
+    ec_twistededwards_projective],
   ../constantine/io/io_bigints,
   ../constantine/tower_field_extensions/extension_fields
 
@@ -233,7 +235,24 @@ func random_long01Seq(rng: var RngState, a: var ExtensionField) =
 # Elliptic curves
 # ------------------------------------------------------------
 
-func random_unsafe(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Aff or ECP_ShortW_Jac)) =
+type ECP = ECP_ShortW_Aff or ECP_ShortW_Prj or ECP_ShortW_Jac or
+           ECP_TwEdwards_Aff or ECP_TwEdwards_Prj
+type ECP_ext = ECP_ShortW_Prj or ECP_ShortW_Jac or
+               ECP_TwEdwards_Prj
+
+template trySetFromCoord[F](a: ECP, fieldElem: F): SecretBool =
+  when a is (ECP_ShortW_Aff or ECP_ShortW_Prj or ECP_ShortW_Jac):
+    trySetFromCoordX(a, fieldElem)
+  else:
+    trySetFromCoordY(a, fieldElem)
+
+template trySetFromCoords[F](a: ECP, fieldElem, scale: F): SecretBool =
+  when a is (ECP_ShortW_Prj or ECP_ShortW_Jac):
+    trySetFromCoordsXandZ(a, fieldElem, scale)
+  else:
+    trySetFromCoordsYandZ(a, fieldElem, scale)
+
+func random_unsafe(rng: var RngState, a: var ECP) =
   ## Initialize a random curve point with Z coordinate == 1
   ## Unsafe: for testing and benchmarking purposes only
   var fieldElem {.noInit.}: a.F
@@ -243,9 +262,9 @@ func random_unsafe(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Aff o
     # Euler's criterion: there are (p-1)/2 squares in a field with modulus `p`
     #                    so we have a probability of ~0.5 to get a good point
     rng.random_unsafe(fieldElem)
-    success = trySetFromCoordX(a, fieldElem)
+    success = trySetFromCoord(a, fieldElem)
 
-func random_unsafe_with_randZ(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Jac)) =
+func random_unsafe_with_randZ(rng: var RngState, a: var ECP_ext) =
   ## Initialize a random curve point with Z coordinate being random
   ## Unsafe: for testing and benchmarking purposes only
   var Z{.noInit.}: a.F
@@ -256,9 +275,9 @@ func random_unsafe_with_randZ(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_S
 
   while not bool(success):
     rng.random_unsafe(fieldElem)
-    success = trySetFromCoordsXandZ(a, fieldElem, Z)
+    success = trySetFromCoords(a, fieldElem, Z)
 
-func random_highHammingWeight(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Aff or ECP_ShortW_Jac)) =
+func random_highHammingWeight(rng: var RngState, a: var ECP) =
   ## Initialize a random curve point with Z coordinate == 1
   ## This will be generated with a biaised RNG with high Hamming Weight
   ## to trigger carry bugs
@@ -269,9 +288,9 @@ func random_highHammingWeight(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_S
     # Euler's criterion: there are (p-1)/2 squares in a field with modulus `p`
     #                    so we have a probability of ~0.5 to get a good point
     rng.random_highHammingWeight(fieldElem)
-    success = trySetFromCoordX(a, fieldElem)
+    success = trySetFromCoord(a, fieldElem)
 
-func random_highHammingWeight_with_randZ(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Jac)) =
+func random_highHammingWeight_with_randZ(rng: var RngState, a: var ECP_ext) =
   ## Initialize a random curve point with Z coordinate == 1
   ## This will be generated with a biaised RNG with high Hamming Weight
   ## to trigger carry bugs
@@ -283,9 +302,9 @@ func random_highHammingWeight_with_randZ(rng: var RngState, a: var (ECP_ShortW_P
 
   while not bool(success):
     rng.random_highHammingWeight(fieldElem)
-    success = trySetFromCoordsXandZ(a, fieldElem, Z)
+    success = trySetFromCoords(a, fieldElem, Z)
 
-func random_long01Seq(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Aff or ECP_ShortW_Jac)) =
+func random_long01Seq(rng: var RngState, a: var ECP) =
   ## Initialize a random curve point with Z coordinate == 1
   ## This will be generated with a biaised RNG
   ## that produces long bitstrings of 0 and 1
@@ -297,9 +316,9 @@ func random_long01Seq(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Af
     # Euler's criterion: there are (p-1)/2 squares in a field with modulus `p`
     #                    so we have a probability of ~0.5 to get a good point
     rng.random_long01Seq(fieldElem)
-    success = trySetFromCoordX(a, fieldElem)
+    success = trySetFromCoord(a, fieldElem)
 
-func random_long01Seq_with_randZ(rng: var RngState, a: var (ECP_ShortW_Prj or ECP_ShortW_Jac)) =
+func random_long01Seq_with_randZ(rng: var RngState, a: var ECP_ext) =
   ## Initialize a random curve point with Z coordinate == 1
   ## This will be generated with a biaised RNG
   ## that produces long bitstrings of 0 and 1
@@ -312,7 +331,7 @@ func random_long01Seq_with_randZ(rng: var RngState, a: var (ECP_ShortW_Prj or EC
 
   while not bool(success):
     rng.random_long01Seq(fieldElem)
-    success = trySetFromCoordsXandZ(a, fieldElem, Z)
+    success = trySetFromCoords(a, fieldElem, Z)
 
 # Generic over any Constantine type
 # ------------------------------------------------------------
@@ -320,7 +339,7 @@ func random_long01Seq_with_randZ(rng: var RngState, a: var (ECP_ShortW_Prj or EC
 func random_unsafe*(rng: var RngState, T: typedesc): T =
   ## Create a random Field or Extension Field or Curve Element
   ## Unsafe: for testing and benchmarking purposes only
-  when T is (ECP_ShortW_Prj or ECP_ShortW_Aff or ECP_ShortW_Jac):
+  when T is ECP:
     rng.random_unsafe(result)
   elif T is SomeNumber:
     cast[T](rng.next()) # TODO: Rely on casting integer actually converting in C (i.e. uint64->uint32 is valid)
@@ -329,7 +348,7 @@ func random_unsafe*(rng: var RngState, T: typedesc): T =
   else: # Fields
     rng.random_unsafe(result)
 
-func random_unsafe_with_randZ*(rng: var RngState, T: typedesc[ECP_ShortW_Prj or ECP_ShortW_Jac]): T =
+func random_unsafe_with_randZ*(rng: var RngState, T: typedesc[ECP_ext]): T =
   ## Create a random curve element with a random Z coordinate
   ## Unsafe: for testing and benchmarking purposes only
   rng.random_unsafe_with_randZ(result)
@@ -337,7 +356,7 @@ func random_unsafe_with_randZ*(rng: var RngState, T: typedesc[ECP_ShortW_Prj or
 func random_highHammingWeight*(rng: var RngState, T: typedesc): T =
   ## Create a random Field or Extension Field or Curve Element
   ## Skewed towards high Hamming Weight
-  when T is (ECP_ShortW_Prj or ECP_ShortW_Aff or ECP_ShortW_Jac):
+  when T is ECP:
     rng.random_highHammingWeight(result)
   elif T is SomeNumber:
     cast[T](rng.next()) # TODO: Rely on casting integer actually converting in C (i.e. uint64->uint32 is valid)
@@ -346,7 +365,7 @@ func random_highHammingWeight*(rng: var RngState, T: typedesc): T =
   else: # Fields
     rng.random_highHammingWeight(result)
 
-func random_highHammingWeight_with_randZ*(rng: var RngState, T: typedesc[ECP_ShortW_Prj or ECP_ShortW_Jac]): T =
+func random_highHammingWeight_with_randZ*(rng: var RngState, T: typedesc[ECP_ext]): T =
   ## Create a random curve element with a random Z coordinate
   ## Skewed towards high Hamming Weight
   rng.random_highHammingWeight_with_randZ(result)
@@ -354,7 +373,7 @@ func random_highHammingWeight_with_randZ*(rng: var RngState, T: typedesc[ECP_Sho
 func random_long01Seq*(rng: var RngState, T: typedesc): T =
   ## Create a random Field or Extension Field or Curve Element
   ## Skewed towards long bitstrings of 0 or 1
-  when T is (ECP_ShortW_Prj or ECP_ShortW_Aff or ECP_ShortW_Jac):
+  when T is ECP:
     rng.random_long01Seq(result)
   elif T is SomeNumber:
     cast[T](rng.next()) # TODO: Rely on casting integer actually converting in C (i.e. uint64->uint32 is valid)
@@ -363,7 +382,7 @@ func random_long01Seq*(rng: var RngState, T: typedesc): T =
   else: # Fields
     rng.random_long01Seq(result)
 
-func random_long01Seq_with_randZ*(rng: var RngState, T: typedesc[ECP_ShortW_Prj or ECP_ShortW_Jac]): T =
+func random_long01Seq_with_randZ*(rng: var RngState, T: typedesc[ECP_ext]): T =
   ## Create a random curve element with a random Z coordinate
   ## Skewed towards long bitstrings of 0 or 1
   rng.random_long01Seq_with_randZ(result)

tests/t_ec_shortw_prj_g1_add_double.nim

@@ -14,7 +14,7 @@ import
   ./t_ec_template
 
 const
-  Iters = 1
+  Iters = 8
 
 run_EC_addition_tests(
     ec = ECP_ShortW_Prj[Fp[BN254_Snarks], NotOnTwist],

tests/t_ec_template.nim

@@ -23,6 +23,8 @@ import
     ec_shortweierstrass_affine,
     ec_shortweierstrass_jacobian,
     ec_shortweierstrass_projective,
+    ec_twistededwards_affine,
+    ec_twistededwards_projective,
     ec_scalar_mul],
     ../constantine/io/[io_bigints, io_fields],
   # Test utilities
@@ -51,6 +53,15 @@ func random_point*(rng: var RngState, EC: typedesc, randZ: bool, gen: RandomGen)
     else:
       result = rng.random_long01Seq_with_randZ(EC)
 
+template pairingGroup(EC: typedesc): string =
+  when EC is (ECP_ShortW_Aff or ECP_ShortW_Prj or ECP_ShortW_Jac):
+    when EC.Tw == NotOnTwist:
+      "G1"
+    else:
+      "G2"
+  else:
+    ""
+
 proc run_EC_addition_tests*(
        ec: typedesc,
        Iters: static int,
@@ -64,12 +75,9 @@ proc run_EC_addition_tests*(
   echo "\n------------------------------------------------------\n"
   echo moduleName, " xoshiro512** seed: ", seed
 
-  when ec.Tw == NotOnTwist:
+  const G1_or_G2 = pairingGroup(ec)
-    const G1_or_G2 = "G1"
-  else:
-    const G1_or_G2 = "G2"
 
-  const testSuiteDesc = "Elliptic curve in Short Weierstrass form with projective coordinates"
+  const testSuiteDesc = "Elliptic curve in " & $ec.F.C.getEquationForm() & " form with projective coordinates"
 
   suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
     test "The infinity point is the neutral element w.r.t. to EC " & G1_or_G2 & " addition":
@@ -215,12 +223,9 @@ proc run_EC_mul_sanity_tests*(
   echo "\n------------------------------------------------------\n"
   echo moduleName, " xoshiro512** seed: ", seed
 
-  when ec.Tw == NotOnTwist:
+  const G1_or_G2 = pairingGroup(ec)
-    const G1_or_G2 = "G1"
-  else:
-    const G1_or_G2 = "G2"
 
-  const testSuiteDesc = "Elliptic curve in Short Weierstrass form"
+  const testSuiteDesc = "Elliptic curve in " & $ec.F.C.getEquationForm() & " form"
 
   suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
     test "EC " & G1_or_G2 & " mul [0]P == Inf":
@@ -313,12 +318,9 @@ proc run_EC_mul_distributive_tests*(
   echo "\n------------------------------------------------------\n"
   echo moduleName, " xoshiro512** seed: ", seed
 
-  when ec.Tw == NotOnTwist:
+  const G1_or_G2 = pairingGroup(ec)
-    const G1_or_G2 = "G1"
-  else:
-    const G1_or_G2 = "G2"
 
-  const testSuiteDesc = "Elliptic curve in Short Weierstrass form"
+  const testSuiteDesc = "Elliptic curve in " & $ec.F.C.getEquationForm() & " form"
 
   suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
 
@@ -383,12 +385,9 @@ proc run_EC_mul_vs_ref_impl*(
   echo "\n------------------------------------------------------\n"
   echo moduleName, " xoshiro512** seed: ", seed
 
-  when ec.Tw == NotOnTwist:
+  const G1_or_G2 = pairingGroup(ec)
-    const G1_or_G2 = "G1"
-  else:
-    const G1_or_G2 = "G2"
 
-  const testSuiteDesc = "Elliptic curve in Short Weierstrass form"
+  const testSuiteDesc = "Elliptic curve in " & $ec.F.C.getEquationForm() & " form"
 
   suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
     test "EC " & G1_or_G2 & " mul constant-time is equivalent to a simple double-and-add algorithm":

tests/t_ec_twedwards_prj_add_double.nim

@@ -0,0 +1,35 @@
+# 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
+  # Internals
+  ../constantine/config/[type_ff, curves],
+  ../constantine/elliptic/ec_twistededwards_projective,
+  # Test utilities
+  ./t_ec_template
+
+const
+  Iters = 8
+
+run_EC_addition_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Curve25519]],
+    Iters = Iters,
+    moduleName = "test_ec_twistededwards_projective_add_double_" & $Curve25519
+  )
+
+run_EC_addition_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Jubjub]],
+    Iters = Iters,
+    moduleName = "test_ec_twistededwards_projective_add_double_" & $Jubjub
+  )
+
+run_EC_addition_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Bandersnatch]],
+    Iters = Iters,
+    moduleName = "test_ec_twistededwards_projective_add_double_" & $Bandersnatch
+  )

tests/t_ec_twedwards_prj_mul_distri.nim

@@ -0,0 +1,36 @@
+# 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
+  # Internals
+  ../constantine/config/[type_ff, curves],
+  ../constantine/elliptic/ec_twistededwards_projective,
+  # Test utilities
+  ./t_ec_template
+
+const
+  Iters = 12
+  ItersMul = Iters div 4
+
+run_EC_mul_distributive_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Curve25519]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_distributive_" & $Curve25519
+  )
+
+run_EC_mul_distributive_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Jubjub]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_distributive_" & $Jubjub
+  )
+
+run_EC_mul_distributive_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Bandersnatch]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_distributive_" & $Bandersnatch
+  )

tests/t_ec_twedwards_prj_mul_sanity.nim

@@ -0,0 +1,36 @@
+# 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
+  # Internals
+  ../constantine/config/[type_ff, curves],
+  ../constantine/elliptic/ec_twistededwards_projective,
+  # Test utilities
+  ./t_ec_template
+  
+const
+  Iters = 12
+  ItersMul = Iters div 4
+
+run_EC_mul_sanity_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Curve25519]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_sanity_" & $Curve25519
+  )
+
+run_EC_mul_sanity_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Jubjub]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_sanity_" & $Jubjub
+  )
+
+run_EC_mul_sanity_tests(
+    ec = ECP_TwEdwards_Prj[Fp[Bandersnatch]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_sanity_" & $Bandersnatch
+  )

tests/t_ec_twedwards_prj_mul_vs_ref.nim

@@ -0,0 +1,36 @@
+# 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
+  # Internals
+  ../constantine/config/[type_ff, curves],
+  ../constantine/elliptic/ec_twistededwards_projective,
+  # Test utilities
+  ./t_ec_template
+
+const
+  Iters = 12
+  ItersMul = Iters div 4
+
+run_EC_mul_vs_ref_impl(
+    ec = ECP_TwEdwards_Prj[Fp[Curve25519]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_vs_ref_" & $Curve25519
+  )
+
+run_EC_mul_vs_ref_impl(
+    ec = ECP_TwEdwards_Prj[Fp[Jubjub]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_vs_ref_" & $Jubjub
+  )
+
+run_EC_mul_vs_ref_impl(
+    ec = ECP_TwEdwards_Prj[Fp[Bandersnatch]],
+    ItersMul = ItersMul,
+    moduleName = "test_ec_twistededwards_projective_mul_vs_ref_" & $Bandersnatch
+  )

tests/t_finite_fields_sqrt.nim

@@ -81,7 +81,7 @@ proc exhaustiveCheck(C: static Curve, modulus: static int) =
           bool not a.isSquare()
           bool not a.sqrt_if_square()
 
-template testImpl(a: untyped): untyped {.dirty.} =
+template testSqrtImpl(a: untyped): untyped {.dirty.} =
   var na{.noInit.}: typeof(a)
   na.neg(a)
 
@@ -103,17 +103,46 @@ template testImpl(a: untyped): untyped {.dirty.} =
 
 proc randomSqrtCheck(C: static Curve) =
   test "Random square root check for " & $Curve(C):
-    for _ in 0 ..< 1: # Iters:
+    for _ in 0 ..< Iters:
       let a = rng.random_unsafe(Fp[C])
-      testImpl(a)
+      testSqrtImpl(a)
 
     for _ in 0 ..< Iters:
       let a = rng.randomHighHammingWeight(Fp[C])
-      testImpl(a)
+      testSqrtImpl(a)
 
     for _ in 0 ..< Iters:
       let a = rng.random_long01Seq(Fp[C])
-      testImpl(a)
+      testSqrtImpl(a)
+
+template testSqrtRatioImpl(u, v: untyped): untyped {.dirty.} =
+  var u_over_v, r{.noInit.}: typeof(v)
+  u_over_v.inv(v)
+  u_over_v *= u
+
+  let qr = r.sqrt_ratio_if_square(u, v)
+  check: bool(qr) == bool(u_over_v.isSquare())
+
+  if bool(qr):
+    r.square()
+    check: bool(r == u_over_v)
+
+proc randomSqrtRatioCheck(C: static Curve) =
+  test "Random square root check for " & $Curve(C):
+    for _ in 0 ..< Iters:
+      let u = rng.random_unsafe(Fp[C])
+      let v = rng.random_unsafe(Fp[C])
+      testSqrtRatioImpl(u, v)
+
+    for _ in 0 ..< Iters:
+      let u = rng.randomHighHammingWeight(Fp[C])
+      let v = rng.randomHighHammingWeight(Fp[C])
+      testSqrtRatioImpl(u, v)
+
+    for _ in 0 ..< Iters:
+      let u = rng.random_long01Seq(Fp[C])
+      let v = rng.random_long01Seq(Fp[C])
+      testSqrtRatioImpl(u, v)
 
 proc main() =
   suite "Modular square root" & " [" & $WordBitwidth & "-bit mode]":
@@ -125,6 +154,14 @@ proc main() =
     randomSqrtCheck BLS12_377 # p ≢ 3 (mod 4)
     randomSqrtCheck BLS12_381
     randomSqrtCheck BW6_761
+    randomSqrtCheck Curve25519
+    randomSqrtCheck Jubjub
+    randomSqrtCheck Bandersnatch
+  
+  suite "Modular sqrt(u/v)" & " [" & $WordBitwidth & "-bit mode]":
+    randomSqrtRatioCheck Curve25519
+    randomSqrtRatioCheck Jubjub
+    randomSqrtRatioCheck Bandersnatch
 
   suite "Modular square root - 32-bit bugs highlighted by property-based testing " & " [" & $WordBitwidth & "-bit mode]":
     # test "FKM12_447 - #30": - Deactivated, we don't support the curve as no one uses it.
@@ -150,11 +187,11 @@ proc main() =
     test "Fp[2^127 - 1] - #61":
       var a: Fp[Mersenne127]
       a.fromHex"0x75bfffefbfffffff7fd9dfd800000000"
-      testImpl(a)
+      testSqrtImpl(a)
 
     test "Fp[2^127 - 1] - #62":
       var a: Fp[Mersenne127]
       a.fromHex"0x7ff7ffffffffffff1dfb7fafc0000000"
-      testImpl(a)
+      testSqrtImpl(a)
 
 main()