Fused projective line eval (#96)

* Reorg line functions to allow for Jacobian eval

* 2x faster Miller loop!!! with fused line eval double

* Support Line Double Fusion for D-Twists

* Implement fused line addition

constantine.nimble

@@ -98,6 +98,8 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   # Edge cases highlighted by past bugs
   ("tests/t_ec_shortw_prj_edge_cases.nim", false),
   # Pairing
+  ("tests/t_pairing_bls12_377_line_functions.nim", false),
+  ("tests/t_pairing_bls12_381_line_functions.nim", false),
   ("tests/t_pairing_mul_fp12_by_lines.nim", false),
   ("tests/t_pairing_cyclotomic_fp12.nim", false),
   ("tests/t_pairing_bn254_nogami_optate.nim", false),

constantine/pairing/lines_common.nim

@@ -0,0 +1,58 @@
+# 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
+  std/typetraits,
+  ../primitives,
+  ../config/[common, curves],
+  ../arithmetic,
+  ../towers,
+  ../elliptic/[
+    ec_shortweierstrass_affine,
+    ec_shortweierstrass_projective
+  ],
+  ../io/io_towers
+
+type
+  Line*[F; twist: static SexticTwist] = object
+    ## Packed line representation over a E'(Fp^k/d)
+    ## with k the embedding degree and d the twist degree
+    ## i.e. for a curve with embedding degree 12 and sextic twist
+    ## F is Fp2
+    ##
+    ## Assuming a Sextic Twist
+    ##
+    ## Out of 6 Fp2 coordinates, 3 are 0 and
+    ## the non-zero coordinates depend on the twist kind.
+    ##
+    ## For a D-twist,
+    ##   (x, y, z) corresponds to an sparse element of Fp12
+    ##   with Fp2 coordinates: xy00z0
+    ## For a M-Twist
+    ##   (x, y, z) corresponds to an sparse element of Fp12
+    ##   with Fp2 coordinates: xyz000
+    x*, y*, z*: F
+
+func toHex*(line: Line, order: static Endianness = bigEndian): string =
+  result = static($line.typeof.genericHead() & '(')
+  for fieldName, fieldValue in fieldPairs(line):
+    when fieldName != "x":
+      result.add ", "
+    result.add fieldName & ": "
+    result.appendHex(fieldValue, order)
+  result.add ")"
+
+# Line evaluation
+# --------------------------------------------------
+
+func line_update*(line: var Line, P: ECP_ShortW_Aff) =
+  ## Update the line evaluation with P
+  ## after addition or doubling
+  ## P in G1
+  line.x *= P.y
+  line.z *= P.x

constantine/pairing/lines_projective.nim

@@ -16,7 +16,10 @@ import
     ec_shortweierstrass_affine,
     ec_shortweierstrass_projective
   ],
-  ../io/io_towers
+  ../io/io_towers,
+  ./lines_common
+
+export lines_common
 
 # ############################################################
 #
@@ -34,58 +37,16 @@ import
 #   and Patrick Longa and Jefferson E. Ricardini, 2013
 #   https://eprint.iacr.org/2013/722.pdf
 #   http://sac2013.irmacs.sfu.ca/slides/s1.pdf
-
 #
-# TODO: Implement fused line doubling and addition
+# - Efficient Implementation of Bilinear Pairings on ARM Processors
-#       from Costello2009 or Aranha2010
+#   Gurleen Grewal, Reza Azarderakhsh,
-#       We don't need the complete formulae in the Miller Loop
+#   Patrick Longa, Shi Hu, and David Jao, 2012
+#   https://eprint.iacr.org/2012/408.pdf
 
-type
+# Line evaluation only
-  Line*[F; twist: static SexticTwist] = object
+# -----------------------------------------------------------------------------
-    ## Packed line representation over a E'(Fp^k/d)
-    ## with k the embedding degree and d the twist degree
-    ## i.e. for a curve with embedding degree 12 and sextic twist
-    ## F is Fp2
-    ##
-    ## Assuming a Sextic Twist
-    ##
-    ## Out of 6 Fp2 coordinates, 3 are 0 and
-    ## the non-zero coordinates depend on the twist kind.
-    ##
-    ## For a D-twist,
-    ##   (x, y, z) corresponds to an sparse element of Fp12
-    ##   with Fp2 coordinates: xy00z0
-    ## For a M-Twist
-    ##   (x, y, z) corresponds to an sparse element of Fp12
-    ##   with Fp2 coordinates: xyz000
-    x*, y*, z*: F
 
-func toHex*(line: Line, order: static Endianness = bigEndian): string =
+func line_eval_double(line: var Line, T: ECP_ShortW_Proj) =
-  result = static($line.typeof.genericHead() & '(')
-  for fieldName, fieldValue in fieldPairs(line):
-    when fieldName != "x":
-      result.add ", "
-    result.add fieldName & ": "
-    result.appendHex(fieldValue, order)
-  result.add ")"
-
-# Line evaluation
-# --------------------------------------------------
-
-func `*=`(a: var Fp2, b: Fp) =
-  ## Multiply an element of Fp2 by an element of Fp
-  # TODO: make generic and move to tower_field_extensions
-  a.c0 *= b
-  a.c1 *= b
-
-func line_update(line: var Line, P: ECP_ShortW_Aff) =
-  ## Update the line evaluation with P
-  ## after addition or doubling
-  ## P in G1
-  line.x *= P.y
-  line.z *= P.x
-
-func line_eval_double*(line: var Line, T: ECP_ShortW_Proj) =
   ## Evaluate the line function for doubling
   ## i.e. the tangent at T
   ##
@@ -158,7 +119,7 @@ func line_eval_double*(line: var Line, T: ECP_ShortW_Proj) =
   B -= v                # B = 3bξ Z² - Y²  (M-twist)
                         # B = 3b Z² - ξ Y² (D-twist)
 
-func line_eval_add*(line: var Line, T: ECP_ShortW_Proj, Q: ECP_ShortW_Aff) =
+func line_eval_add(line: var Line, T: ECP_ShortW_Proj, Q: ECP_ShortW_Aff) =
   ## Evaluate the line function for addition
   ## i.e. the line between T and Q
   ##
@@ -206,13 +167,135 @@ func line_eval_add*(line: var Line, T: ECP_ShortW_Proj, Q: ECP_ShortW_Aff) =
 
   C.neg()     # C = -(Y₁-Z₁Y₂)
 
+func line_eval_fused_double(line: var Line, T: var ECP_ShortW_Proj) =
+  ## Fused line evaluation and elliptic point doubling
+  # Grewal et al, 2012 adapted to Scott 2019 line notation
+  var A {.noInit.}, B {.noInit.}, C {.noInit.}: Line.F
+  var E {.noInit.}, F {.noInit.}, G {.noInit.}: Line.F
+  template H: untyped = line.x
+  const b3 = 3*Line.F.C.getCoefB()
+
+  var snrY = T.y
+  when Line.F.C.getSexticTwist() == D_Twist:
+    snrY *= SexticNonResidue
+
+  A.prod(T.x, snrY)
+  A.div2()          # A = XY/2
+  B.square(T.y)     # B = Y²
+  C.square(T.z)     # C = Z²
+
+  var snrB = B
+  when Line.F.C.getSexticTwist() == D_Twist:
+    snrB *= SexticNonResidue
+
+  E = C
+  E *= b3
+  when Line.F.C.getSexticTwist() == M_Twist:
+    E *= SexticNonResidue # E = 3b'Z² = 3bξ Z²
+
+  F = E
+  F *= 3            # F = 3E = 9bZ²
+  G.sum(snrB, F)
+  G.div2()          # G = (B+F)/2
+  H.sum(T.y, T.z)
+  H.square()
+  H -= B
+  H -= C            # lx = H = (Y+Z)²-(B+C)= 2YZ
+
+  line.z.square(T.x)
+  line.z *= 3       # lz = 3X²
+  when Line.F.C.getSexticTwist() == D_Twist:
+    line.z *= SexticNonResidue
+
+  line.y.diff(E, snrB) # ly = E-B = 3b'Z² - Y²
+
+  # In-place modification: invalidates `T.` calls
+  T.x.diff(snrB, F)
+  T.x *= A          # X₃ = A(B-F) = XY/2.(Y²-9b'Z²)
+                    # M-twist: XY/2.(Y²-9bξZ²)
+                    # D-Twist: ξXY/2.(Y²ξ-9bZ²)
+
+  T.y.square(G)
+  E.square()
+  E *= 3
+  T.y -= E          # Y₃ = G² - 3E² = (Y²+9b'Z²)²/4 - 3*(3b'Z²)²
+                    # M-twist: (Y²+9bξZ²)²/4 - 3*(3bξZ²)²
+                    # D-Twist: (ξY²+9bZ²)²/4 - 3*(3bZ²)²
+
+  when Line.F.C.getSexticTwist() == D_Twist:
+    H *= SexticNonResidue
+  T.z.prod(snrB, H) # Z₃ = BH = Y²((Y+Z)² - (Y²+Z²)) = 2Y³Z
+                    # M-twist: 2Y³Z
+                    # D-twist: 2ξ²Y³Z
+
+  # Correction for Fp4 towering
+  H.neg()           # lx = -H
+  when Line.F.C.getSexticTwist() == M_Twist:
+    H *= SexticNonResidue
+    # else: the SNR is already integrated in H
+
+func line_eval_fused_add(line: var Line, T: var ECP_ShortW_Proj, Q: ECP_ShortW_Aff) =
+  ## Fused line evaluation and elliptic point addition
+  # Grewal et al, 2012 adapted to Scott 2019 line notation
+  var
+    A {.noInit.}: Line.F
+    B {.noInit.}: Line.F
+    C {.noInit.}: Line.F
+    D {.noInit.}: Line.F
+    E {.noInit.}: Line.F
+    F {.noInit.}: Line.F
+    G {.noInit.}: Line.F
+    H {.noInit.}: Line.F
+    I {.noInit.}: Line.F
+
+  template lambda: untyped = line.x
+  template theta: untyped = line.z
+  template J: untyped = line.y
+
+  A.prod(Q.y, T.z)
+  B.prod(Q.x, T.z)
+  theta.diff(T.y, A)  # θ = Y₁ - Z₁X₂
+  lambda.diff(T.x, B) # λ = X₁ - Z₁X₂
+  C.square(theta)
+  D.square(lambda)
+  E.prod(D, lambda)
+  F.prod(T.z, C)
+  G.prod(T.x, D)
+  H.double(G)
+  H.diffAlias(F, H)
+  H += E
+  I.prod(T.y, E)
+
+  T.x.prod(theta, Q.x)
+  T.y.prod(lambda, Q.y)
+  J.diff(T.x, T.y)
+
+  # EC addition
+  T.x.prod(lambda, H)
+
+  T.y.diff(G, H)
+  T.y *= theta
+  T.y -= I
+
+  T.z *= E
+
+  # Line evaluation
+  theta.neg()
+  when ECP_ShortW_Proj.F.C.getSexticTwist() == M_Twist:
+    lambda *= SexticNonResidue # A = ξ (X₁ - Z₁X₂)
+
+# Public proc
+# -----------------------------------------------------------------------------
+
 func line_double*(line: var Line, T: var ECP_ShortW_Proj, P: ECP_ShortW_Aff) =
   ## Doubling step of the Miller loop
   ## T in G2, P in G1
   ##
   ## Compute lt,t(P)
-  ##
+  when true:
-  # TODO fused line doubling from Costello 2009, Grewal 2012, Aranha 2013
+    line_eval_fused_double(line, T)
+    line.line_update(P)
+  else:
     line_eval_double(line, T)
     line.line_update(P)
     T.double()
@@ -225,7 +308,10 @@ func line_add*[C](
   ## T and Q in G2, P in G1
   ##
   ## Compute lt,q(P)
-  # TODO fused line addition from Costello 2009, Grewal 2012, Aranha 2013
+  when true:
+    line_eval_fused_add(line, T, Q)
+    line.line_update(P)
+  else:
     line_eval_add(line, T, Q)
     line.line_update(P)
     T += Q

constantine/pairing/mul_fp12_by_lines.nim

@@ -41,43 +41,6 @@ import
 # 𝔽p12 by line - Sparse functions
 # ----------------------------------------------------------------
 
-func mul_sparse_by_y0*[C: static Curve](r: var Fp4[C], a: Fp4[C], b: Fp2[C]) =
-  ## Sparse multiplication of an Fp4 element
-  ## with coordinates (a₀, a₁) by (b₀, 0)
-  r.c0.prod(a.c0, b)
-  r.c1.prod(a.c1, b)
-
-func mul_sparse_by_0y*[C: static Curve](r: var Fp4[C], a: Fp4[C], b: Fp2[C]) =
-  ## Sparse multiplication of an Fp4 element
-  ## with coordinates (a₀, a₁) by (0, b₁)
-  r.c0.prod(a.c1, b)
-  r.c0 *= NonResidue
-  r.c1.prod(a.c0, b)
-
-func mul_sparse_by_0y0*[C: static Curve](r: var Fp6[C], a: Fp6[C], b: Fp2[C]) =
-  ## Sparse multiplication of an Fp6 element
-  ## with coordinates (a₀, a₁, a₂) by (0, b₁, 0)
-  # TODO: make generic and move to tower_field_extensions
-
-  # v0 = a0 b0 = 0
-  # v1 = a1 b1
-  # v2 = a2 b2 = 0
-  #
-  # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
-  #    = ξ (a1 b1 + a2 b1 - v1)
-  #    = ξ a2 b1
-  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2
-  #    = a0 b1 + a1 b1 - v1
-  #    = a0 b1
-  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
-  #    = v1
-  #    = a1 b1
-
-  r.c0.prod(a.c2, b)
-  r.c0 *= ξ
-  r.c1.prod(a.c0, b)
-  r.c2.prod(a.c1, b)
-
 func mul_by_line_xy0*[C: static Curve, twist: static SexticTwist](
        r: var Fp6[C],
        a: Fp6[C],

constantine/towers.nim

@@ -197,3 +197,49 @@ func `*`*(_: typedesc[γ], a: Fp4): Fp4 {.noInit, inline.} =
 
 func `*=`*(a: var Fp4, _: typedesc[γ]) {.inline.} =
   a = γ * a
+
+# Sparse functions
+# ----------------------------------------------------------------
+
+func `*=`*(a: var Fp2, b: Fp) =
+  ## Multiply an element of Fp2 by an element of Fp
+  # TODO: make generic and move to tower_field_extensions
+  a.c0 *= b
+  a.c1 *= b
+
+func mul_sparse_by_y0*[C: static Curve](r: var Fp4[C], a: Fp4[C], b: Fp2[C]) =
+  ## Sparse multiplication of an Fp4 element
+  ## with coordinates (a₀, a₁) by (b₀, 0)
+  r.c0.prod(a.c0, b)
+  r.c1.prod(a.c1, b)
+
+func mul_sparse_by_0y*[C: static Curve](r: var Fp4[C], a: Fp4[C], b: Fp2[C]) =
+  ## Sparse multiplication of an Fp4 element
+  ## with coordinates (a₀, a₁) by (0, b₁)
+  r.c0.prod(a.c1, b)
+  r.c0 *= NonResidue
+  r.c1.prod(a.c0, b)
+
+func mul_sparse_by_0y0*[C: static Curve](r: var Fp6[C], a: Fp6[C], b: Fp2[C]) =
+  ## Sparse multiplication of an Fp6 element
+  ## with coordinates (a₀, a₁, a₂) by (0, b₁, 0)
+  # TODO: make generic and move to tower_field_extensions
+
+  # v0 = a0 b0 = 0
+  # v1 = a1 b1
+  # v2 = a2 b2 = 0
+  #
+  # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
+  #    = ξ (a1 b1 + a2 b1 - v1)
+  #    = ξ a2 b1
+  # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2
+  #    = a0 b1 + a1 b1 - v1
+  #    = a0 b1
+  # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
+  #    = v1
+  #    = a1 b1
+
+  r.c0.prod(a.c2, b)
+  r.c0 *= ξ
+  r.c1.prod(a.c0, b)
+  r.c2.prod(a.c1, b)

tests/t_pairing_bls12_377_line_functions.nim

@@ -0,0 +1,114 @@
+# 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
+  # Standard library
+  std/[tables, unittest, times],
+  # Internals
+  ../constantine/config/common,
+  ../constantine/[arithmetic, primitives],
+  ../constantine/towers,
+  ../constantine/config/curves,
+  ../constantine/io/io_towers,
+  ../constantine/elliptic/[
+    ec_shortweierstrass_affine,
+    ec_shortweierstrass_projective,
+    ec_scalar_mul],
+  ../constantine/pairing/lines_projective,
+  # Test utilities
+  ../helpers/[prng_unsafe, static_for]
+
+const
+  Iters = 4
+  TestCurves = [
+    BLS12_377
+  ]
+
+type
+  RandomGen = enum
+    Uniform
+    HighHammingWeight
+    Long01Sequence
+
+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 "test_pairing_bls12_377_line_functions xoshiro512** seed: ", seed
+
+func random_point*(rng: var RngState, EC: typedesc, gen: RandomGen): EC {.noInit.} =
+  if gen == Uniform:
+    result = rng.random_unsafe(EC)
+  elif gen == HighHammingWeight:
+    result = rng.random_highHammingWeight(EC)
+  else:
+    result = rng.random_long01Seq(EC)
+
+func random_point*(rng: var RngState, EC: typedesc, randZ: bool, gen: RandomGen): EC {.noInit.} =
+  if not randZ:
+    if gen == Uniform:
+      result = rng.random_unsafe(EC)
+    elif gen == HighHammingWeight:
+      result = rng.random_highHammingWeight(EC)
+    else:
+      result = rng.random_long01Seq(EC)
+  else:
+    if gen == Uniform:
+      result = rng.random_unsafe_with_randZ(EC)
+    elif gen == HighHammingWeight:
+      result = rng.random_highHammingWeight_with_randZ(EC)
+    else:
+      result = rng.random_long01Seq_with_randZ(EC)
+
+suite "Pairing - Line Functions on BLS12-377" & " [" & $WordBitwidth & "-bit mode]":
+  test "Line double - lt,t(P)":
+    proc test_line_double(C: static Curve, randZ: bool, gen: RandomGen) =
+      for _ in 0 ..< Iters:
+        let P = rng.random_point(ECP_ShortW_Aff[Fp[C]], gen)
+        var T = rng.random_point(ECP_ShortW_Proj[Fp2[C]], randZ, gen)
+        let Q = rng.random_point(ECP_ShortW_Proj[Fp2[C]], randZ, gen)
+        var l: Line[Fp2[C], C.getSexticTwist()]
+
+        var T2: typeof(Q)
+        T2.double(T)
+        l.line_double(T, P)
+
+        doAssert: bool(T == T2)
+
+    staticFor(curve, TestCurves):
+      test_line_double(curve, randZ = false, gen = Uniform)
+      test_line_double(curve, randZ = true, gen = Uniform)
+      test_line_double(curve, randZ = false, gen = HighHammingWeight)
+      test_line_double(curve, randZ = true, gen = HighHammingWeight)
+      test_line_double(curve, randZ = false, gen = Long01Sequence)
+      test_line_double(curve, randZ = true, gen = Long01Sequence)
+
+  test "Line add - lt,q(P)":
+    proc test_line_add(C: static Curve, randZ: bool, gen: RandomGen) =
+      for _ in 0 ..< Iters:
+        let P = rng.random_point(ECP_ShortW_Aff[Fp[C]], gen)
+        let Q = rng.random_point(ECP_ShortW_Proj[Fp2[C]], randZ, gen)
+        var T = rng.random_point(ECP_ShortW_Proj[Fp2[C]], randZ, gen)
+        var l: Line[Fp2[C], C.getSexticTwist()]
+
+        var TQ{.noInit.}: typeof(T)
+        TQ.sum(T, Q)
+
+        var Qaff{.noInit.}: ECP_ShortW_Aff[Fp2[C]]
+        Qaff.affineFromProjective(Q)
+        l.line_add(T, Qaff, P)
+
+        doAssert: bool(T == TQ)
+
+    staticFor(curve, TestCurves):
+      test_line_add(curve, randZ = false, gen = Uniform)
+      test_line_add(curve, randZ = true, gen = Uniform)
+      test_line_add(curve, randZ = false, gen = HighHammingWeight)
+      test_line_add(curve, randZ = true, gen = HighHammingWeight)
+      test_line_add(curve, randZ = false, gen = Long01Sequence)
+      test_line_add(curve, randZ = true, gen = Long01Sequence)