From 6ecbedbd0944d70c765ea5d0206af00df29b22c3 Mon Sep 17 00:00:00 2001 From: Mamy Ratsimbazafy Date: Sat, 26 Sep 2020 09:16:29 +0200 Subject: [PATCH] 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 | 2 + benchmarks/bench_ec_g2.nim | 2 + benchmarks/bench_elliptic_template.nim | 12 +- constantine.nimble | 1 + .../elliptic/ec_weierstrass_projective.nim | 206 ++++++++++++------ constantine/pairing/lines_projective.nim | 5 +- tests/t_ec_template.nim | 46 +++- tests/t_ec_wstrass_prj_g1_mixed_add.nim | 30 +++ ..._ec_wstrass_prj_g2_mixed_add_bls12_381.nim | 24 ++ ..._wstrass_prj_g2_mixed_add_bn254_snarks.nim | 24 ++ 10 files changed, 275 insertions(+), 77 deletions(-) create mode 100644 tests/t_ec_wstrass_prj_g1_mixed_add.nim create mode 100644 tests/t_ec_wstrass_prj_g2_mixed_add_bls12_381.nim create mode 100644 tests/t_ec_wstrass_prj_g2_mixed_add_bn254_snarks.nim diff --git a/benchmarks/bench_ec_g1.nim b/benchmarks/bench_ec_g1.nim index 7684fdb..1b4f284 100644 --- a/benchmarks/bench_ec_g1.nim +++ b/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) diff --git a/benchmarks/bench_ec_g2.nim b/benchmarks/bench_ec_g2.nim index c37f299..01dd424 100644 --- a/benchmarks/bench_ec_g2.nim +++ b/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) diff --git a/benchmarks/bench_elliptic_template.nim b/benchmarks/bench_elliptic_template.nim index 586a10a..f8440b1 100644 --- a/benchmarks/bench_elliptic_template.nim +++ b/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 diff --git a/constantine.nimble b/constantine.nimble index 0cb970c..44b910c 100644 --- a/constantine.nimble +++ b/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), diff --git a/constantine/elliptic/ec_weierstrass_projective.nim b/constantine/elliptic/ec_weierstrass_projective.nim index 697ffe1..3de0bf0 100644 --- a/constantine/elliptic/ec_weierstrass_projective.nim +++ b/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) - # - (Y1 Z2 + Y2 Z1)(a X1 X2 + 3b(X1 Z2 + X2 Z1) - a² Z1 Z2) - # Y3 = (3 X1 X2 + a Z1 Z2)(a X1 X2 + 3b (X1 Z2 + X2 Z1) - a² Z1 Z2) - # + (Y1 Y2 + a (X1 Z2 + X2 Z1) + 3b Z1 Z2)(Y1 Y2 - a(X1 Z2 + X2 Z1) - 3b Z1 Z2) - # 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) + # X₃ = (X₁Y₂ + X₂Y₁)(Y₁Y₂ - a (X₁Z₂ + X₂Z₁) - 3bZ₁Z₂) + # - (Y₁Z₂ + Y₂Z₁)(aX₁X₂ + 3b(X₁Z₂ + X₂Z₁) - a²Z₁Z₂) + # Y₃ = (3X₁X₂ + aZ₁Z₂)(aX₁X₂ + 3b(X₁Z₂ + X₂Z₁) - a²Z₁Z₂) + # + (Y₁Y₂ + a (X₁Z₂ + X₂Z₁) + 3bZ₁Z₂)(Y₁Y₂ - a(X₁Z₂ + X₂Z₁) - 3bZ₁Z₂) + # 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) - # − 3b(Y1 Z2 + Y2 Z1)(X1 Z2 + X2 Z1) - # Y3 = (Y1 Y2 + 3b Z1 Z2)(Y1 Y2 − 3b Z1 Z2) - # + 9b X1 X2 (X1 Z2 + X2 Z1) - # Z3= (Y1 Z2 + Y2 Z1)(Y1 Y2 + 3b Z1 Z2) + 3 X1 X2 (X1 Y2 + X2 Y1) - t0.prod(P.x, Q.x) # 1. t0 <- X1 X2 - t1.prod(P.y, Q.y) # 2. t1 <- Y1 Y2 - t2.prod(P.z, Q.z) # 3. t2 <- Z1 Z2 - t3.sum(P.x, P.y) # 4. t3 <- X1 + Y1 - t4.sum(Q.x, Q.y) # 5. t4 <- X2 + Y2 - t3 *= t4 # 6. t3 <- t3 * t4 - t4.sum(t0, t1) # 7. t4 <- t0 + t1 - t3 -= t4 # 8. t3 <- t3 - t4 t3 = (X1 + Y1)(X2 + Y2) - (X1 X2 + Y1 Y2) = X1.Y2 + X2.Y1 + # X₃ = (X₁Y₂ + X₂Y₁)(Y₁Y₂ − 3bZ₁Z₂) + # − 3b(Y₁Z₂ + Y₂Z₁)(X₁Z₂ + X₂Z₁) + # Y₃ = (Y₁Y₂ + 3bZ₁Z₂)(Y₁Y₂ − 3bZ₁Z₂) + # + 9bX₁X₂ (X₁Z₂ + X₂Z₁) + # Z₃= (Y₁Z₂ + Y₂Z₁)(Y₁Y₂ + 3bZ₁Z₂) + 3X₁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₂ + t2.prod(P.z, Q.z) # 3. t₂ <- Z₁Z₂ + t3.sum(P.x, P.y) # 4. t₃ <- X₁ + Y₁ + t4.sum(Q.x, Q.y) # 5. t₄ <- X₂ + Y₂ + t3 *= t4 # 6. t₃ <- t₃ * t₄ + t4.sum(t0, t1) # 7. t₄ <- t₀ + t₁ + 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 - r.x.sum(Q.y, Q.z) # 10. X3 <- Y2 + Z2 - t4 *= r.x # 11. t4 <- t4 X3 - r.x.sum(t1, t2) # 12. X3 <- t1 + t2 X3 = Y1 Y2 + Z1 Z2 - t4 -= r.x # 13. t4 <- t4 - X3 t4 = (Y1 + Z1)(Y2 + Z2) - (Y1 Y2 + Z1 Z2) = Y1 Z2 + Y2 Z1 + t4.sum(P.y, P.z) # 9. t₄ <- Y₁ + Z₁ + r.x.sum(Q.y, Q.z) # 10. X₃ <- Y₂ + Z₂ + t4 *= r.x # 11. t₄ <- t₄ X₃ + r.x.sum(t1, t2) # 12. X₃ <- t₁ + t₂ X₃ = Y₁Y₂ + Z₁Z₂ + 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.y.sum(Q.x, Q.z) # 15. Y3 <- X2 + Z2 - r.x *= r.y # 16. X3 <- X3 Y3 X3 = (X1 Z1)(X2 Z2) - r.y.sum(t0, t2) # 17. Y3 <- t0 + t2 Y3 = X1 X2 + Z1 Z2 - 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.x.sum(P.x, P.z) # 14. X₃ <- X₁ + Z₁ + r.y.sum(Q.x, Q.z) # 15. Y₃ <- X₂ + Z₂ + r.x *= r.y # 16. X₃ <- X₃ Y₃ X₃ = (X₁Z₁)(X₂Z₂) + r.y.sum(t0, t2) # 17. Y₃ <- t₀ + t₂ Y₃ = X₁ X₂ + Z₁ Z₂ + 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 - t0 += r.x # 20. t0 <- X3 + t0 t0 = 3 X1 X2 - t2 *= b3 # 21. t2 <- b3 t2 t2 = 3b Z1 Z2 + r.x.double(t0) # 19. X₃ <- t₀ + t₀ X₃ = 2 X₁X₂ + t0 += r.x # 20. t₀ <- X₃ + t₀ t₀ = 3 X₁X₂ + 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 - t1 -= t2 # 23. t1 <- t1 - t2 t1 = Y1 Y2 - 3b Z1 Z2 - r.y *= b3 # 24. Y3 <- b3 Y3 Y3 = 3b(X1 Z2 + X2 Z1) + r.z.sum(t1, t2) # 22. Z₃ <- t₁ + t₂ Z₃ = Y₁Y₂ + 3bZ₁Z₂ + t1 -= t2 # 23. t₁ <- t₁ - t₂ t₁ = Y₁Y₂ - 3bZ₁Z₂ + 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) - t2.prod(t3, t1) # 26. t2 <- t3 t1 t2 = (X1 Y2 + X2 Y1) (Y1 Y2 - 3b Z1 Z2) - 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.y *= t0 # 28. Y3 <- Y3 t0 Y3 = 9b X1 X2 (X1 Z2 + X2 Z1) - t1 *= r.z # 29. t1 <- t1 Z3 t1 = (Y1 Y2 - 3b Z1 Z2)(Y1 Y2 + 3b Z1 Z2) - 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) - t0 *= t3 # 31. t0 <- t0 t3 t0 = 3 X1 X2 (X1.Y2 + X2.Y1) - r.z *= t4 # 32. Z3 <- Z3 t4 Z3 = (Y1 Y2 + 3b Z1 Z2)(Y1 Z2 + Y2 Z1) - 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.x.prod(t4, r.y) # 25. X₃ <- t₄ Y₃ X₃ = 3b(Y₁Z₂ + Y₂Z₁)(X₁Z₂ + X₂Z₁) + t2.prod(t3, t1) # 26. t₂ <- t₃ t₁ t₂ = (X₁Y₂ + X₂Y₁) (Y₁Y₂ - 3bZ₁Z₂) + 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. Y₃ <- Y₃ t₀ Y₃ = 9bX₁X₂ (X₁Z₂ + X₂Z₁) + t1 *= r.z # 29. t₁ <- t₁ Z₃ t₁ = (Y₁Y₂ - 3bZ₁Z₂)(Y₁Y₂ + 3bZ₁Z₂) + 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. t₀ <- t₀ t₃ t₀ = 3X₁X₂ (X₁Y₂ + X₂Y₁) + r.z *= t4 # 32. Z₃ <- Z₃ t₄ Z₃ = (Y₁Y₂ + 3bZ₁Z₂)(Y₁Z₂ + Y₂Z₁) + 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²) - # Y3 = (Y² - 9bZ²)(Y² + 3bZ²) + 24bY²Z² - # Z3 = 8Y³Z + # X₃ = 2XY(Y² - 9bZ²) + # Y₃ = (Y² - 9bZ²)(Y² + 3bZ²) + 24bY²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 - r.z.double(t0) # 2. Z3 <- t0 + t0 - r.z.double() # 3. Z3 <- Z3 + Z3 - r.z.double() # 4. Z3 <- Z3 + Z3 Z3 = 8Y² - t1.prod(snrY, P.z) # 5. t1 <- Y Z - t2.square(P.z) # 6. t2 <- Z Z - t2 *= b3 # 7. t2 <- b3 t2 + t0.square(P.y) # 1. t₀ <- Y Y + r.z.double(t0) # 2. Z₃ <- t₀ + t₀ + r.z.double() # 3. Z₃ <- Z₃ + Z₃ + r.z.double() # 4. Z₃ <- Z₃ + Z₃ Z₃ = 8Y² + t1.prod(snrY, P.z) # 5. t₁ <- Y Z + t2.square(P.z) # 6. t₂ <- Z Z + 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.y.sum(t0, t2) # 9. Y3 <- t0 + t2 - r.z *= t1 # 10. Z3 <- t1 Z3 - t1.double(t2) # 11. t1 <- t2 + t2 - t2 += t1 # 12. t2 <- t1 + t2 - t0 -= t2 # 13. t0 <- t0 - t2 - r.y *= t0 # 14. Y3 <- t0 Y3 - r.y += r.x # 15. Y3 <- X3 + Y3 - t1.prod(P.x, snrY) # 16. t1 <- X Y - r.x.prod(t0, t1) # 17. X3 <- t0 t1 - r.x.double() # 18. X3 <- X3 + X3 + r.x.prod(t2, r.z) # 8. X₃ <- t₂ Z₃ + r.y.sum(t0, t2) # 9. Y₃ <- t₀ + t₂ + r.z *= t1 # 10. Z₃ <- t₁ Z₃ + t1.double(t2) # 11. t₁ <- t₂ + t₂ + t2 += t1 # 12. t₂ <- t₁ + t₂ + t0 -= t2 # 13. t₀ <- t₀ - t₂ + r.y *= t0 # 14. Y₃ <- t₀ Y₃ + r.y += r.x # 15. Y₃ <- X₃ + Y₃ + t1.prod(P.x, snrY) # 16. t₁ <- X Y + r.x.prod(t0, t1) # 17. X₃ <- t₀ t₁ + 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) diff --git a/constantine/pairing/lines_projective.nim b/constantine/pairing/lines_projective.nim index 9e115d2..28ab263 100644 --- a/constantine/pairing/lines_projective.nim +++ b/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 - var QProj {.noInit.}: ECP_SWei_Proj[Fp2[C]] - QProj.projectiveFromAffine(Q) - T += QProj + T += Q diff --git a/tests/t_ec_template.nim b/tests/t_ec_template.nim index 0840475..b1e53f1 100644 --- a/tests/t_ec_template.nim +++ b/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_projective, ec_scalar_mul], + ../constantine/elliptic/[ec_weierstrass_affine, 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) diff --git a/tests/t_ec_wstrass_prj_g1_mixed_add.nim b/tests/t_ec_wstrass_prj_g1_mixed_add.nim new file mode 100644 index 0000000..df6da2b --- /dev/null +++ b/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 + ) diff --git a/tests/t_ec_wstrass_prj_g2_mixed_add_bls12_381.nim b/tests/t_ec_wstrass_prj_g2_mixed_add_bls12_381.nim new file mode 100644 index 0000000..eb87576 --- /dev/null +++ b/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 + ) diff --git a/tests/t_ec_wstrass_prj_g2_mixed_add_bn254_snarks.nim b/tests/t_ec_wstrass_prj_g2_mixed_add_bn254_snarks.nim new file mode 100644 index 0000000..067c6fe --- /dev/null +++ b/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 + )