Mixed addition (#90)

* ptrettier comments

* Implement mixed addition on G1

* Test for mixed addition in G2 and use it for Miller Loop

benchmarks/bench_ec_g1.nim

@@ -49,6 +49,8 @@ proc main() =
     const curve = AvailableCurves[i]
     addBench(ECP_SWei_Proj[Fp[curve]], Iters)
     separator()
+    mixedAddBench(ECP_SWei_Proj[Fp[curve]], Iters)
+    separator()
     doublingBench(ECP_SWei_Proj[Fp[curve]], Iters)
     separator()
     scalarMulUnsafeDoubleAddBench(ECP_SWei_Proj[Fp[curve]], MulIters)

benchmarks/bench_ec_g2.nim

@@ -50,6 +50,8 @@ proc main() =
     const curve = AvailableCurves[i]
     addBench(ECP_SWei_Proj[Fp2[curve]], Iters)
     separator()
+    mixedAddBench(ECP_SWei_Proj[Fp2[curve]], Iters)
+    separator()
     doublingBench(ECP_SWei_Proj[Fp2[curve]], Iters)
     separator()
     scalarMulUnsafeDoubleAddBench(ECP_SWei_Proj[Fp2[curve]], MulIters)

benchmarks/bench_elliptic_template.nim

@@ -17,7 +17,7 @@ import
   ../constantine/config/[curves, common],
   ../constantine/arithmetic,
   ../constantine/io/io_bigints,
-  ../constantine/elliptic/[ec_weierstrass_projective, ec_scalar_mul, ec_endomorphism_accel],
+  ../constantine/elliptic/[ec_weierstrass_affine, ec_weierstrass_projective, ec_scalar_mul, ec_endomorphism_accel],
   # Helpers
   ../helpers/[prng_unsafe, static_for],
   ./platforms,
@@ -135,6 +135,16 @@ proc addBench*(T: typedesc, iters: int) =
   bench("EC Add " & G1_or_G2, T, iters):
     r.sum(P, Q)
 
+proc mixedAddBench*(T: typedesc, iters: int) =
+  const G1_or_G2 = when T.F is Fp: "G1" else: "G2"
+  var r {.noInit.}: T
+  let P = rng.random_unsafe(T)
+  let Q = rng.random_unsafe(T)
+  var Qaff: ECP_SWei_Aff[T.F]
+  Qaff.affineFromProjective(Q)
+  bench("EC Mixed Addition " & G1_or_G2, T, iters):
+    r.madd(P, Qaff)
+
 proc doublingBench*(T: typedesc, iters: int) =
   const G1_or_G2 = when T.F is Fp: "G1" else: "G2"
   var r {.noInit.}: T

constantine.nimble

@@ -54,6 +54,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   ("tests/t_ec_wstrass_prj_g1_mul_sanity.nim", false),
   ("tests/t_ec_wstrass_prj_g1_mul_distri.nim", false),
   ("tests/t_ec_wstrass_prj_g1_mul_vs_ref.nim", false),
+  ("tests/t_ec_wstrass_prj_g1_mixed_add.nim", false),
   # Elliptic curve arithmetic G2
   ("tests/t_ec_wstrass_prj_g2_add_double_bn254_snarks.nim", false),
   ("tests/t_ec_wstrass_prj_g2_mul_sanity_bn254_snarks.nim", false),

constantine/elliptic/ec_weierstrass_projective.nim

@@ -153,11 +153,11 @@ func sum*[F](
   #
   # with the indices 1 corresponding to ``P``, 2 to ``Q`` and 3 to the result ``r``
   #
-  # X3 = (X1 Y2 + X2 Y1)(Y1 Y2 - a(X1 Z2 + X2 Z1) - 3b Z1 Z2)
+  # X₃ = (X₁Y₂ + X₂Y₁)(Y₁Y₂ - a (X₁Z₂ + X₂Z₁) - 3bZ₁Z₂)
-  #      - (Y1 Z2 + Y2 Z1)(a X1 X2 + 3b(X1 Z2 + X2 Z1) - a² Z1 Z2)
+  #      - (Y₁Z₂ + Y₂Z₁)(aX₁X₂ + 3b(X₁Z₂ + X₂Z₁) - a²Z₁Z₂)
-  # Y3 = (3 X1 X2 + a Z1 Z2)(a X1 X2 + 3b (X1 Z2 + X2 Z1) - a² Z1 Z2)
+  # Y₃ = (3X₁X₂ + aZ₁Z₂)(aX₁X₂ + 3b(X₁Z₂ + X₂Z₁) - a²Z₁Z₂)
-  #      + (Y1 Y2 + a (X1 Z2 + X2 Z1) + 3b Z1 Z2)(Y1 Y2 - a(X1 Z2 + X2 Z1) - 3b Z1 Z2)
+  #      + (Y₁Y₂ + a (X₁Z₂ + X₂Z₁) + 3bZ₁Z₂)(Y₁Y₂ - a(X₁Z₂ + X₂Z₁) - 3bZ₁Z₂)
-  # Z3 = (Y1 Z2 + Y2 Z1)(Y1 Y2 + a(X1 Z2 + X2 Z1) + 3b Z1 Z2) + (X1 Y2 + X2 Y1)(3 X1 X2 + a Z1 Z2)
+  # Z₃ = (Y₁Z₂ + Y₂Z₁)(Y₁Y₂ + a(X₁Z₂ + X₂Z₁) + 3bZ₁Z₂) + (X₁Y₂ + X₂Y₁)(3X₁X₂ + aZ₁Z₂)
   #
   # Cost: 12M + 3 mul(a) + 2 mul(3b) + 23 a
 
@@ -170,55 +170,117 @@ func sum*[F](
     # Algorithm 7 for curves: y² = x³ + b
     # 12M + 2 mul(3b) + 19A
     #
-    # X3 = (X1 Y2 + X2 Y1)(Y1 Y2 − 3b Z1 Z2)
+    # X₃ = (X₁Y₂ + X₂Y₁)(Y₁Y₂ − 3bZ₁Z₂)
-    #     − 3b(Y1 Z2 + Y2 Z1)(X1 Z2 + X2 Z1)
+    #     − 3b(Y₁Z₂ + Y₂Z₁)(X₁Z₂ + X₂Z₁)
-    # Y3 = (Y1 Y2 + 3b Z1 Z2)(Y1 Y2 − 3b Z1 Z2)
+    # Y₃ = (Y₁Y₂ + 3bZ₁Z₂)(Y₁Y₂ − 3bZ₁Z₂)
-    #     + 9b X1 X2 (X1 Z2 + X2 Z1)
+    #     + 9bX₁X₂ (X₁Z₂ + X₂Z₁)
-    # Z3= (Y1 Z2 + Y2 Z1)(Y1 Y2 + 3b Z1 Z2) + 3 X1 X2 (X1 Y2 + X2 Y1)
+    # Z₃= (Y₁Z₂ + Y₂Z₁)(Y₁Y₂ + 3bZ₁Z₂) + 3X₁X₂ (X₁Y₂ + X₂Y₁)
-    t0.prod(P.x, Q.x)         # 1.  t0 <- X1 X2
+    t0.prod(P.x, Q.x)         # 1.  t₀ <- X₁X₂
-    t1.prod(P.y, Q.y)         # 2.  t1 <- Y1 Y2
+    t1.prod(P.y, Q.y)         # 2.  t₁ <- Y₁Y₂
-    t2.prod(P.z, Q.z)         # 3.  t2 <- Z1 Z2
+    t2.prod(P.z, Q.z)         # 3.  t₂ <- Z₁Z₂
-    t3.sum(P.x, P.y)          # 4.  t3 <- X1 + Y1
+    t3.sum(P.x, P.y)          # 4.  t₃ <- X₁ + Y₁
-    t4.sum(Q.x, Q.y)          # 5.  t4 <- X2 + Y2
+    t4.sum(Q.x, Q.y)          # 5.  t₄ <- X₂ + Y₂
-    t3 *= t4                  # 6.  t3 <- t3 * t4
+    t3 *= t4                  # 6.  t₃ <- t₃ * t₄
-    t4.sum(t0, t1)            # 7.  t4 <- t0 + t1
+    t4.sum(t0, t1)            # 7.  t₄ <- t₀ + t₁
-    t3 -= t4                  # 8.  t3 <- t3 - t4   t3 = (X1 + Y1)(X2 + Y2) - (X1 X2 + Y1 Y2) = X1.Y2 + X2.Y1
+    t3 -= t4                  # 8.  t₃ <- t₃ - t₄   t₃ = (X₁ + Y₁)(X₂ + Y₂) - (X₁X₂ + Y₁Y₂) = X₁Y₂ + X₂Y₁
     when F is Fp2 and F.C.getSexticTwist() == D_Twist:
       t3 *= SexticNonResidue
-    t4.sum(P.y, P.z)          # 9.  t4 <- Y1 + Z1
+    t4.sum(P.y, P.z)          # 9.  t₄ <- Y₁ + Z₁
-    r.x.sum(Q.y, Q.z)         # 10. X3 <- Y2 + Z2
+    r.x.sum(Q.y, Q.z)         # 10. X₃ <- Y₂ + Z₂
-    t4 *= r.x                 # 11. t4 <- t4 X3
+    t4 *= r.x                 # 11. t₄ <- t₄ X₃
-    r.x.sum(t1, t2)           # 12. X3 <- t1 + t2   X3 = Y1 Y2 + Z1 Z2
+    r.x.sum(t1, t2)           # 12. X₃ <- t₁ + t₂   X₃ = Y₁Y₂ + Z₁Z₂
-    t4 -= r.x                 # 13. t4 <- t4 - X3   t4 = (Y1 + Z1)(Y2 + Z2) - (Y1 Y2 + Z1 Z2) = Y1 Z2 + Y2 Z1
+    t4 -= r.x                 # 13. t₄ <- t₄ - X₃   t₄ = (Y₁ + Z₁)(Y₂ + Z₂) - (Y₁Y₂ + Z₁Z₂) = Y₁Z₂ + Y₂Z₁
     when F is Fp2 and F.C.getSexticTwist() == D_Twist:
       t4 *= SexticNonResidue
-    r.x.sum(P.x, P.z)         # 14. X3 <- X1 + Z1
+    r.x.sum(P.x, P.z)         # 14. X₃ <- X₁ + Z₁
-    r.y.sum(Q.x, Q.z)         # 15. Y3 <- X2 + Z2
+    r.y.sum(Q.x, Q.z)         # 15. Y₃ <- X₂ + Z₂
-    r.x *= r.y                # 16. X3 <- X3 Y3     X3 = (X1 Z1)(X2 Z2)
+    r.x *= r.y                # 16. X₃ <- X₃ Y₃     X₃ = (X₁Z₁)(X₂Z₂)
-    r.y.sum(t0, t2)           # 17. Y3 <- t0 + t2   Y3 = X1 X2 + Z1 Z2
+    r.y.sum(t0, t2)           # 17. Y₃ <- t₀ + t₂   Y₃ = X₁ X₂ + Z₁ Z₂
-    r.y.diffAlias(r.x, r.y)   # 18. Y3 <- X3 - Y3   Y3 = (X1 + Z1)(X2 + Z2) - (X1 X2 + Z1 Z2) = X1 Z2 + X2 Z1
+    r.y.diffAlias(r.x, r.y)   # 18. Y₃ <- X₃ - Y₃   Y₃ = (X₁ + Z₁)(X₂ + Z₂) - (X₁ X₂ + Z₁ Z₂) = X₁Z₂ + X₂Z₁
     when F is Fp2 and F.C.getSexticTwist() == D_Twist:
       t0 *= SexticNonResidue
       t1 *= SexticNonResidue
-    r.x.double(t0)            # 19. X3 <- t0 + t0   X3 = 2 X1 X2
+    r.x.double(t0)            # 19. X₃ <- t₀ + t₀   X₃ = 2 X₁X₂
-    t0 += r.x                 # 20. t0 <- X3 + t0   t0 = 3 X1 X2
+    t0 += r.x                 # 20. t₀ <- X₃ + t₀   t₀ = 3 X₁X₂
-    t2 *= b3                  # 21. t2 <- b3 t2     t2 = 3b Z1 Z2
+    t2 *= b3                  # 21. t₂ <- 3b t₂     t₂ = 3bZ₁Z₂
     when F is Fp2 and F.C.getSexticTwist() == M_Twist:
       t2 *= SexticNonResidue
-    r.z.sum(t1, t2)           # 22. Z3 <- t1 + t2   Z3 = Y1 Y2 + 3b Z1 Z2
+    r.z.sum(t1, t2)           # 22. Z₃ <- t₁ + t₂   Z₃ = Y₁Y₂ + 3bZ₁Z₂
-    t1 -= t2                  # 23. t1 <- t1 - t2   t1 = Y1 Y2 - 3b Z1 Z2
+    t1 -= t2                  # 23. t₁ <- t₁ - t₂   t₁ = Y₁Y₂ - 3bZ₁Z₂
-    r.y *= b3                 # 24. Y3 <- b3 Y3     Y3 = 3b(X1 Z2 + X2 Z1)
+    r.y *= b3                 # 24. Y₃ <- 3b Y₃     Y₃ = 3b(X₁Z₂ + X₂Z₁)
     when F is Fp2 and F.C.getSexticTwist() == M_Twist:
       r.y *= SexticNonResidue
-    r.x.prod(t4, r.y)         # 25. X3 <- t4 Y3     X3 = 3b(Y1 Z2 + Y2 Z1)(X1 Z2 + X2 Z1)
+    r.x.prod(t4, r.y)         # 25. X₃ <- t₄ Y₃     X₃ = 3b(Y₁Z₂ + Y₂Z₁)(X₁Z₂ + X₂Z₁)
-    t2.prod(t3, t1)           # 26. t2 <- t3 t1     t2 = (X1 Y2 + X2 Y1) (Y1 Y2 - 3b Z1 Z2)
+    t2.prod(t3, t1)           # 26. t₂ <- t₃ t₁     t₂ = (X₁Y₂ + X₂Y₁) (Y₁Y₂ - 3bZ₁Z₂)
-    r.x.diffAlias(t2, r.x)    # 27. X3 <- t2 - X3   X3 = (X1 Y2 + X2 Y1) (Y1 Y2 - 3b Z1 Z2) - 3b(Y1 Z2 + Y2 Z1)(X1 Z2 + X2 Z1)
+    r.x.diffAlias(t2, r.x)    # 27. X₃ <- t₂ - X₃   X₃ = (X₁Y₂ + X₂Y₁) (Y₁Y₂ - 3bZ₁Z₂) - 3b(Y₁Z₂ + Y₂Z₁)(X₁Z₂ + X₂Z₁)
-    r.y *= t0                 # 28. Y3 <- Y3 t0     Y3 = 9b X1 X2 (X1 Z2 + X2 Z1)
+    r.y *= t0                 # 28. Y₃ <- Y₃ t₀     Y₃ = 9bX₁X₂ (X₁Z₂ + X₂Z₁)
-    t1 *= r.z                 # 29. t1 <- t1 Z3     t1 = (Y1 Y2 - 3b Z1 Z2)(Y1 Y2 + 3b Z1 Z2)
+    t1 *= r.z                 # 29. t₁ <- t₁ Z₃     t₁ = (Y₁Y₂ - 3bZ₁Z₂)(Y₁Y₂ + 3bZ₁Z₂)
-    r.y += t1                 # 30. Y3 <- t1 + Y3   Y3 = (Y1 Y2 + 3b Z1 Z2)(Y1 Y2 - 3b Z1 Z2) + 9b X1 X2 (X1 Z2 + X2 Z1)
+    r.y += t1                 # 30. Y₃ <- t₁ + Y₃   Y₃ = (Y₁Y₂ + 3bZ₁Z₂)(Y₁Y₂ - 3bZ₁Z₂) + 9bX₁X₂ (X₁Z₂ + X₂Z₁)
-    t0 *= t3                  # 31. t0 <- t0 t3     t0 = 3 X1 X2 (X1.Y2 + X2.Y1)
+    t0 *= t3                  # 31. t₀ <- t₀ t₃     t₀ = 3X₁X₂ (X₁Y₂ + X₂Y₁)
-    r.z *= t4                 # 32. Z3 <- Z3 t4     Z3 = (Y1 Y2 + 3b Z1 Z2)(Y1 Z2 + Y2 Z1)
+    r.z *= t4                 # 32. Z₃ <- Z₃ t₄     Z₃ = (Y₁Y₂ + 3bZ₁Z₂)(Y₁Z₂ + Y₂Z₁)
-    r.z += t0                 # 33. Z3 <- Z3 + t0   Z3 = (Y1 Z2 + Y2 Z1)(Y1 Y2 + 3b Z1 Z2) + 3 X1 X2 (X1.Y2 + X2.Y1)
+    r.z += t0                 # 33. Z₃ <- Z₃ + t₀   Z₃ = (Y₁Z₂ + Y₂Z₁)(Y₁Y₂ + 3bZ₁Z₂) + 3X₁X₂ (X₁Y₂ + X₂Y₁)
+  else:
+    {.error: "Not implemented.".}
+
+func madd*[F](
+       r: var ECP_SWei_Proj[F],
+       P: ECP_SWei_Proj[F], Q: ECP_SWei_Aff[F]
+     ) =
+  ## Elliptic curve mixed addition for Short Weierstrass curves
+  ## with p in Projective coordinates and Q in affine coordinates
+  ##
+  ##   R = P + Q
+
+  # TODO: static doAssert odd order
+
+  when F.C.getCoefA() == 0:
+    var t0 {.noInit.}, t1 {.noInit.}, t2 {.noInit.}, t3 {.noInit.}, t4 {.noInit.}: F
+    const b3 = 3 * F.C.getCoefB()
+
+    # Algorithm 8 for curves: y² = x³ + b
+    # X₃ = (X₁Y₂ + X₂Y₁)(Y₁Y₂ − 3bZ₁)
+    #     − 3b(Y₁ + Y₂Z₁)(X₁ + X₂Z₁)
+    # Y₃ = (Y₁Y₂ + 3bZ₁)(Y₁Y₂ − 3bZ₁)
+    #     + 9bX₁X₂ (X₁ + X₂Z₁)
+    # Z₃= (Y₁ + Y₂Z₁)(Y₁Y₂ + 3bZ₁) + 3 X₁X₂ (X₁Y₂ + X₂Y₁)
+    t0.prod(P.x, Q.x)         # 1.  t₀ <- X₁ X₂
+    t1.prod(P.y, Q.y)         # 2.  t₁ <- Y₁ Y₂
+    t3.sum(P.x, P.y)          # 3.  t₃ <- X₁ + Y₁ ! error in paper
+    t4.sum(Q.x, Q.y)          # 4.  t₄ <- X₂ + Y₂ ! error in paper
+    t3 *= t4                  # 5.  t₃ <- t₃ * t₄
+    t4.sum(t0, t1)            # 6.  t₄ <- t₀ + t₁
+    t3 -= t4                  # 7.  t₃ <- t₃ - t₄, t₃ = (X₁ + Y₁)(X₂ + Y₂) - (X₁ X₂ + Y₁ Y₂) = X₁Y₂ + X₂Y₁
+    when F is Fp2 and F.C.getSexticTwist() == D_Twist:
+      t3 *= SexticNonResidue
+    t4.prod(Q.y, P.z)         # 8.  t₄ <- Y₂ Z₁
+    t4 += P.y                 # 9.  t₄ <- t₄ + Y₁, t₄ = Y₁+Y₂Z₁
+    when F is Fp2 and F.C.getSexticTwist() == D_Twist:
+      t4 *= SexticNonResidue
+    r.y.prod(Q.x, P.z)        # 10. Y₃ <- X₂ Z₁
+    r.y += P.x                # 11. Y₃ <- Y₃ + X₁, Y₃ = X₁ + X₂Z₁
+    when F is Fp2 and F.C.getSexticTwist() == D_Twist:
+      t0 *= SexticNonResidue
+      t1 *= SexticNonResidue
+    r.x.double(t0)            # 12. X₃ <- t₀ + t₀
+    t0 += r.x                 # 13. t₀ <- X₃ + t₀, t₀ = 3X₁X₂
+    t2 = P.z
+    t2 *= b3                  # 14. t₂ <- 3bZ₁
+    when F is Fp2 and F.C.getSexticTwist() == M_Twist:
+      t2 *= SexticNonResidue
+    r.z.sum(t1, t2)           # 15. Z₃ <- t₁ + t₂, Z₃ = Y₁Y₂ + 3bZ₁
+    t1 -= t2                  # 16. t₁ <- t₁ - t₂, t₁ = Y₁Y₂ - 3bZ₁
+    r.y *= b3                 # 17. Y₃ <- 3bY₃,    Y₃ = 3b(X₁ + X₂Z₁)
+    when F is Fp2 and F.C.getSexticTwist() == M_Twist:
+      r.y *= SexticNonResidue
+    r.x.prod(t4, r.y)         # 18. X₃ <- t₄ Y₃,   X₃ = (Y₁ + Y₂Z₁) 3b(X₁ + X₂Z₁)
+    t2.prod(t3, t1)           # 19. t₂ <- t₃ t₁,   t₂ = (X₁Y₂ + X₂Y₁)(Y₁Y₂ - 3bZ₁)
+    r.x.diffAlias(t2, r.x)    # 20. X₃ <- t₂ - X₃, X₃ = (X₁Y₂ + X₂Y₁)(Y₁Y₂ - 3bZ₁) - 3b(Y₁ + Y₂Z₁)(X₁ + X₂Z₁)
+    r.y *= t0                 # 21. Y₃ <- Y₃ t₀,   Y₃ = 9bX₁X₂ (X₁ + X₂Z₁)
+    t1 *= r.z                 # 22. t₁ <- t₁ Z₃,   t₁ = (Y₁Y₂ - 3bZ₁)(Y₁Y₂ + 3bZ₁)
+    r.y += t1                 # 23. Y₃ <- t₁ + Y₃, Y₃ = (Y₁Y₂ + 3bZ₁)(Y₁Y₂ - 3bZ₁) + 9bX₁X₂ (X₁ + X₂Z₁)
+    t0 *= t3                  # 31. t₀ <- t₀ t₃,   t₀ = 3X₁X₂ (X₁Y₂ + X₂Y₁)
+    r.z *= t4                 # 32. Z₃ <- Z₃ t₄,   Z₃ = (Y₁Y₂ + 3bZ₁)(Y₁ + Y₂Z₁)
+    r.z += t0                 # 33. Z₃ <- Z₃ + t₀, Z₃ = (Y₁ + Y₂Z₁)(Y₁Y₂ + 3bZ₁) + 3X₁X₂ (X₁Y₂ + X₂Y₁)
   else:
     {.error: "Not implemented.".}
 
@@ -249,11 +311,11 @@ func double*[F](
   #   Joost Renes and Craig Costello and Lejla Batina, 2015
   #   https://eprint.iacr.org/2015/1060
   #
-  # X3 = 2XY (Y² - 2aXZ - 3bZ²)
+  # X₃ = 2XY (Y² - 2aXZ - 3bZ²)
   #      - 2YZ (aX² + 6bXZ - a²Z²)
-  # Y3 = (Y² + 2aXZ + 3bZ²)(Y² - 2aXZ - 3bZ²)
+  # Y₃ = (Y² + 2aXZ + 3bZ²)(Y² - 2aXZ - 3bZ²)
   #      + (3X² + aZ²)(aX² + 6bXZ - a²Z²)
-  # Z3 = 8Y³Z
+  # Z₃ = 8Y³Z
   #
   # Cost: 8M + 3S + 3 mul(a) + 2 mul(3b) + 15a
 
@@ -264,43 +326,50 @@ func double*[F](
     # Algorithm 9 for curves:
     # 6M + 2S + 1 mul(3b) + 9a
     #
-    # X3 = 2XY(Y² - 9bZ²)
+    # X₃ = 2XY(Y² - 9bZ²)
-    # Y3 = (Y² - 9bZ²)(Y² + 3bZ²) + 24bY²Z²
+    # Y₃ = (Y² - 9bZ²)(Y² + 3bZ²) + 24bY²Z²
-    # Z3 = 8Y³Z
+    # Z₃ = 8Y³Z
     snrY = P.y
     when F is Fp2 and F.C.getSexticTwist() == D_Twist:
       snrY *= SexticNonResidue
       t0.square(P.y)
       t0 *= SexticNonResidue
     else:
-      t0.square(P.y)          # 1.  t0 <- Y Y
+      t0.square(P.y)          # 1.  t₀ <- Y Y
-    r.z.double(t0)            # 2.  Z3 <- t0 + t0
+    r.z.double(t0)            # 2.  Z₃ <- t₀ + t₀
-    r.z.double()              # 3.  Z3 <- Z3 + Z3
+    r.z.double()              # 3.  Z₃ <- Z₃ + Z₃
-    r.z.double()              # 4.  Z3 <- Z3 + Z3   Z3 = 8Y²
+    r.z.double()              # 4.  Z₃ <- Z₃ + Z₃   Z₃ = 8Y²
-    t1.prod(snrY, P.z)        # 5.  t1 <- Y Z
+    t1.prod(snrY, P.z)        # 5.  t₁ <- Y Z
-    t2.square(P.z)            # 6.  t2 <- Z Z
+    t2.square(P.z)            # 6.  t₂ <- Z Z
-    t2 *= b3                  # 7.  t2 <- b3 t2
+    t2 *= b3                  # 7.  t₂ <- 3b t₂
     when F is Fp2 and F.C.getSexticTwist() == M_Twist:
       t2 *= SexticNonResidue
-    r.x.prod(t2, r.z)         # 8.  X3 <- t2 Z3
+    r.x.prod(t2, r.z)         # 8.  X₃ <- t₂ Z₃
-    r.y.sum(t0, t2)           # 9.  Y3 <- t0 + t2
+    r.y.sum(t0, t2)           # 9.  Y₃ <- t₀ + t₂
-    r.z *= t1                 # 10. Z3 <- t1 Z3
+    r.z *= t1                 # 10. Z₃ <- t₁ Z₃
-    t1.double(t2)             # 11. t1 <- t2 + t2
+    t1.double(t2)             # 11. t₁ <- t₂ + t₂
-    t2 += t1                  # 12. t2 <- t1 + t2
+    t2 += t1                  # 12. t₂ <- t₁ + t₂
-    t0 -= t2                  # 13. t0 <- t0 - t2
+    t0 -= t2                  # 13. t₀ <- t₀ - t₂
-    r.y *= t0                 # 14. Y3 <- t0 Y3
+    r.y *= t0                 # 14. Y₃ <- t₀ Y₃
-    r.y += r.x                # 15. Y3 <- X3 + Y3
+    r.y += r.x                # 15. Y₃ <- X₃ + Y₃
-    t1.prod(P.x, snrY)        # 16. t1 <- X Y
+    t1.prod(P.x, snrY)        # 16. t₁ <- X Y
-    r.x.prod(t0, t1)          # 17. X3 <- t0 t1
+    r.x.prod(t0, t1)          # 17. X₃ <- t₀ t₁
-    r.x.double()              # 18. X3 <- X3 + X3
+    r.x.double()              # 18. X₃ <- X₃ + X₃
   else:
     {.error: "Not implemented.".}
 
 func `+=`*[F](P: var ECP_SWei_Proj[F], Q: ECP_SWei_Proj[F]) =
+  ## In-place point addition
+  # TODO test for aliasing support
   var tmp {.noInit.}: ECP_SWei_Proj[F]
   tmp.sum(P, Q)
   P = tmp
 
+func `+=`*[F](P: var ECP_SWei_Proj[F], Q: ECP_SWei_Aff[F]) =
+  ## In-place mixed point addition
+  # used in line_addition
+  P.madd(P, Q)
+
 func double*[F](P: var ECP_SWei_Proj[F]) =
   var tmp {.noInit.}: ECP_SWei_Proj[F]
   tmp.double(P)
@@ -316,9 +385,6 @@ func diff*[F](r: var ECP_SWei_Proj[F],
   r.sum(P, nQ)
 
 func affineFromProjective*[F](aff: var ECP_SWei_Aff[F], proj: ECP_SWei_Proj) =
-  # TODO: for now we reuse projective coordinate backend
-  # TODO: simultaneous inversion
-
   var invZ {.noInit.}: F
   invZ.inv(proj.z)
 

constantine/pairing/lines_projective.nim

@@ -228,7 +228,4 @@ func line_add*[C](
   # TODO fused line addition from Costello 2009, Grewal 2012, Aranha 2013
   line_eval_add(line, T, Q)
   line.line_update(P)
-  # TODO: mixed addition
+  T += Q
-  var QProj {.noInit.}: ECP_SWei_Proj[Fp2[C]]
-  QProj.projectiveFromAffine(Q)
-  T += QProj

tests/t_ec_template.nim

@@ -19,8 +19,7 @@ import
   ../constantine/config/[common, curves],
   ../constantine/arithmetic,
   ../constantine/towers,
-  ../constantine/io/io_bigints,
+  ../constantine/elliptic/[ec_weierstrass_affine, ec_weierstrass_projective, ec_scalar_mul],
-  ../constantine/elliptic/[ec_weierstrass_projective, ec_scalar_mul],
   # Test utilities
   ../helpers/prng_unsafe,
   ./support/ec_reference_scalar_mult
@@ -409,3 +408,46 @@ proc run_EC_mul_vs_ref_impl*(
       test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = true, gen = HighHammingWeight)
       test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = false, gen = Long01Sequence)
       test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = true, gen = Long01Sequence)
+
+proc run_EC_mixed_add_impl*(
+       ec: typedesc,
+       Iters: static int,
+       moduleName: string
+     ) =
+
+  # Random seed for reproducibility
+  var rng: RngState
+  let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+  rng.seed(seed)
+  echo "\n------------------------------------------------------\n"
+  echo moduleName, " xoshiro512** seed: ", seed
+
+  when ec.F is Fp:
+    const G1_or_G2 = "G1"
+  else:
+    const G1_or_G2 = "G2"
+
+  const testSuiteDesc = "Elliptic curve mixed addition for Short Weierstrass form"
+
+  suite testSuiteDesc & " - " & $ec & " - [" & $WordBitwidth & "-bit mode]":
+    test "EC " & G1_or_G2 & " mixed addition is consistent with general addition":
+      proc test(EC: typedesc, bits: static int, randZ: bool, gen: RandomGen) =
+        for _ in 0 ..< Iters:
+          let a = rng.random_point(EC, randZ, gen)
+          let b = rng.random_point(EC, randZ, gen)
+          var bAff: ECP_SWei_Aff[EC.F]
+          bAff.affineFromProjective(b)
+
+          var r_generic, r_mixed: EC
+
+          r_generic.sum(a, b)
+          r_mixed.madd(a, bAff)
+
+          check: bool(r_generic == r_mixed)
+
+      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = false, gen = Uniform)
+      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = true, gen = Uniform)
+      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = false, gen = HighHammingWeight)
+      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = true, gen = HighHammingWeight)
+      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = false, gen = Long01Sequence)
+      test(ec, bits = ec.F.C.getCurveOrderBitwidth(), randZ = true, gen = Long01Sequence)

tests/t_ec_wstrass_prj_g1_mixed_add.nim

@@ -0,0 +1,30 @@
+# 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/curves,
+  ../constantine/elliptic/ec_weierstrass_projective,
+  ../constantine/arithmetic,
+  # Test utilities
+  ./t_ec_template
+
+const
+  Iters = 12
+
+run_EC_mixed_add_impl(
+    ec = ECP_SWei_Proj[Fp[BN254_Snarks]],
+    Iters = Iters,
+    moduleName = "test_ec_weierstrass_projective_mixed_add_" & $BN254_Snarks
+  )
+
+run_EC_mixed_add_impl(
+    ec = ECP_SWei_Proj[Fp[BLS12_381]],
+    Iters = Iters,
+    moduleName = "test_ec_weierstrass_projective_mixed_add_" & $BLS12_381
+  )

tests/t_ec_wstrass_prj_g2_mixed_add_bls12_381.nim

@@ -0,0 +1,24 @@
+# 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/curves,
+  ../constantine/elliptic/ec_weierstrass_projective,
+  ../constantine/towers,
+  # Test utilities
+  ./t_ec_template
+
+const
+  Iters = 12
+
+run_EC_mixed_add_impl(
+    ec = ECP_SWei_Proj[Fp2[BLS12_381]],
+    Iters = Iters,
+    moduleName = "test_ec_weierstrass_projective_mixed_add_" & $BLS12_381
+  )

tests/t_ec_wstrass_prj_g2_mixed_add_bn254_snarks.nim

@@ -0,0 +1,24 @@
+# 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/curves,
+  ../constantine/elliptic/ec_weierstrass_projective,
+  ../constantine/towers,
+  # Test utilities
+  ./t_ec_template
+
+const
+  Iters = 12
+
+run_EC_mixed_add_impl(
+    ec = ECP_SWei_Proj[Fp2[BN254_Snarks]],
+    Iters = Iters,
+    moduleName = "test_ec_weierstrass_projective_mixed_add_" & $BN254_Snarks
+  )