diff --git a/benchmarks/bench_pairing_bls12_377.nim b/benchmarks/bench_pairing_bls12_377.nim
index b7ce8b1..ac23eee 100644
--- a/benchmarks/bench_pairing_bls12_377.nim
+++ b/benchmarks/bench_pairing_bls12_377.nim
@@ -36,11 +36,11 @@ proc main() =
     const curve = AvailableCurves[i]
     lineDoubleBench(curve, Iters)
     lineAddBench(curve, Iters)
-    mulFp12byLine_xyz000_Bench(curve, Iters)
-    mulLinebyLine_xyz000_Bench(curve, Iters)
-    mulFp12by_abcdefghij00_Bench(curve, Iters)
-    mulFp12_by_2lines_v1_xyz000_Bench(curve, Iters)
-    mulFp12_by_2lines_v2_xyz000_Bench(curve, Iters)
+    mulFp12byLine_Bench(curve, Iters)
+    mulLinebyLine_Bench(curve, Iters)
+    mulFp12by_prod2lines_Bench(curve, Iters)
+    mulFp12_by_2lines_v1_Bench(curve, Iters)
+    mulFp12_by_2lines_v2_Bench(curve, Iters)
     separator()
     mulBench(curve, Iters)
     sqrBench(curve, Iters)
diff --git a/benchmarks/bench_pairing_bls12_381.nim b/benchmarks/bench_pairing_bls12_381.nim
index 2854b9d..3f6c25c 100644
--- a/benchmarks/bench_pairing_bls12_381.nim
+++ b/benchmarks/bench_pairing_bls12_381.nim
@@ -36,11 +36,11 @@ proc main() =
     const curve = AvailableCurves[i]
     lineDoubleBench(curve, Iters)
     lineAddBench(curve, Iters)
-    mulFp12byLine_xy000z_Bench(curve, Iters)
-    mulLinebyLine_xy000z_Bench(curve, Iters)
-    mulFp12by_abcd00efghij_Bench(curve, Iters)
-    mulFp12_by_2lines_v1_xy000z_Bench(curve, Iters)
-    mulFp12_by_2lines_v2_xy000z_Bench(curve, Iters)
+    mulFp12byLine_Bench(curve, Iters)
+    mulLinebyLine_Bench(curve, Iters)
+    mulFp12by_prod2lines_Bench(curve, Iters)
+    mulFp12_by_2lines_v1_Bench(curve, Iters)
+    mulFp12_by_2lines_v2_Bench(curve, Iters)
     separator()
     mulBench(curve, Iters)
     sqrBench(curve, Iters)
diff --git a/benchmarks/bench_pairing_bn254_nogami.nim b/benchmarks/bench_pairing_bn254_nogami.nim
index ffd0653..a182c82 100644
--- a/benchmarks/bench_pairing_bn254_nogami.nim
+++ b/benchmarks/bench_pairing_bn254_nogami.nim
@@ -36,11 +36,11 @@ proc main() =
     const curve = AvailableCurves[i]
     lineDoubleBench(curve, Iters)
     lineAddBench(curve, Iters)
-    mulFp12byLine_xyz000_Bench(curve, Iters)
-    mulLinebyLine_xyz000_Bench(curve, Iters)
-    mulFp12by_abcdefghij00_Bench(curve, Iters)
-    mulFp12_by_2lines_v1_xyz000_Bench(curve, Iters)
-    mulFp12_by_2lines_v2_xyz000_Bench(curve, Iters)
+    mulFp12byLine_Bench(curve, Iters)
+    mulLinebyLine_Bench(curve, Iters)
+    mulFp12by_prod2lines_Bench(curve, Iters)
+    mulFp12_by_2lines_v1_Bench(curve, Iters)
+    mulFp12_by_2lines_v2_Bench(curve, Iters)
     separator()
     mulBench(curve, Iters)
     sqrBench(curve, Iters)
diff --git a/benchmarks/bench_pairing_bn254_snarks.nim b/benchmarks/bench_pairing_bn254_snarks.nim
index 1625029..186b6b8 100644
--- a/benchmarks/bench_pairing_bn254_snarks.nim
+++ b/benchmarks/bench_pairing_bn254_snarks.nim
@@ -36,11 +36,11 @@ proc main() =
     const curve = AvailableCurves[i]
     lineDoubleBench(curve, Iters)
     lineAddBench(curve, Iters)
-    mulFp12byLine_xyz000_Bench(curve, Iters)
-    mulLinebyLine_xyz000_Bench(curve, Iters)
-    mulFp12by_abcdefghij00_Bench(curve, Iters)
-    mulFp12_by_2lines_v1_xyz000_Bench(curve, Iters)
-    mulFp12_by_2lines_v2_xyz000_Bench(curve, Iters)
+    mulFp12byLine_Bench(curve, Iters)
+    mulLinebyLine_Bench(curve, Iters)
+    mulFp12by_prod2lines_Bench(curve, Iters)
+    mulFp12_by_2lines_v1_Bench(curve, Iters)
+    mulFp12_by_2lines_v2_Bench(curve, Iters)
     separator()
     mulBench(curve, Iters)
     sqrBench(curve, Iters)
diff --git a/benchmarks/bench_pairing_template.nim b/benchmarks/bench_pairing_template.nim
index ee03345..b22812c 100644
--- a/benchmarks/bench_pairing_template.nim
+++ b/benchmarks/bench_pairing_template.nim
@@ -75,63 +75,35 @@ proc lineAddBench*(C: static Curve, iters: int) =
   bench("Line add", C, iters):
     line.line_add(T, Q, P)
 
-proc mulFp12byLine_xyz000_Bench*(C: static Curve, iters: int) =
+proc mulFp12byLine_Bench*(C: static Curve, iters: int) =
   var line: Line[Fp2[C]]
   var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
   let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
   line.line_double(T, P)
   var f = rng.random_unsafe(Fp12[C])
 
-  bench("Mul 𝔽p12 by line xyz000", C, iters):
-    f.mul_sparse_by_line_xyz000(line)
+  bench("Mul 𝔽p12 by line", C, iters):
+    f.mul_by_line(line)
 
-proc mulFp12byLine_xy000z_Bench*(C: static Curve, iters: int) =
-  var line: Line[Fp2[C]]
-  var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
-  let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
-  line.line_double(T, P)
-  var f = rng.random_unsafe(Fp12[C])
-
-  bench("Mul 𝔽p12 by line xy000z", C, iters):
-    f.mul_sparse_by_line_xy000z(line)
-
-proc mulLinebyLine_xyz000_Bench*(C: static Curve, iters: int) =
+proc mulLinebyLine_Bench*(C: static Curve, iters: int) =
   var l0, l1: Line[Fp2[C]]
   var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
   let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
   l0.line_double(T, P)
   l1.line_double(T, P)
-  var f = rng.random_unsafe(Fp12[C])
+  var f {.noInit.}: Fp12[C]
 
-  bench("Mul line xyz000 by line xyz000", C, iters):
-    f.prod_xyz000_xyz000_into_abcdefghij00(l0, l1)
+  bench("Mul line by line", C, iters):
+    f.prod_from_2_lines(l0, l1)
 
-proc mulLinebyLine_xy000z_Bench*(C: static Curve, iters: int) =
-  var l0, l1: Line[Fp2[C]]
-  var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
-  let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
-  l0.line_double(T, P)
-  l1.line_double(T, P)
-  var f = rng.random_unsafe(Fp12[C])
-
-  bench("Mul line xy000z by line xy000z", C, iters):
-    f.prod_xy000z_xy000z_into_abcd00efghij(l0, l1)
-
-proc mulFp12by_abcdefghij00_Bench*(C: static Curve, iters: int) =
+proc mulFp12by_prod2lines_Bench*(C: static Curve, iters: int) =
   var f = rng.random_unsafe(Fp12[C])
   let g = rng.random_unsafe(Fp12[C])
 
-  bench("Mul 𝔽p12 by abcdefghij00", C, iters):
-    f.mul_sparse_by_abcdefghij00(g)
+  bench("Mul 𝔽p12 by product of 2 lines", C, iters):
+    f.mul_by_prod_of_2_lines(g)
 
-proc mulFp12by_abcd00efghij_Bench*(C: static Curve, iters: int) =
-  var f = rng.random_unsafe(Fp12[C])
-  let g = rng.random_unsafe(Fp12[C])
-
-  bench("Mul 𝔽p12 by abcd00efghij", C, iters):
-    f.mul_sparse_by_abcd00efghij(g)
-
-proc mulFp12_by_2lines_v1_xyz000_Bench*(C: static Curve, iters: int) =
+proc mulFp12_by_2lines_v1_Bench*(C: static Curve, iters: int) =
   var l0, l1: Line[Fp2[C]]
   var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
   let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
@@ -140,10 +112,10 @@ proc mulFp12_by_2lines_v1_xyz000_Bench*(C: static Curve, iters: int) =
   var f = rng.random_unsafe(Fp12[C])
 
   bench("mulFp12 by 2 lines v1", C, iters):
-    f.mul_sparse_by_line_xyz000(l0)
-    f.mul_sparse_by_line_xyz000(l1)
+    f.mul_by_line(l0)
+    f.mul_by_line(l1)
 
-proc mulFp12_by_2lines_v2_xyz000_Bench*(C: static Curve, iters: int) =
+proc mulFp12_by_2lines_v2_Bench*(C: static Curve, iters: int) =
   var l0, l1: Line[Fp2[C]]
   var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
   let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
@@ -153,33 +125,8 @@ proc mulFp12_by_2lines_v2_xyz000_Bench*(C: static Curve, iters: int) =
 
   bench("mulFp12 by 2 lines v2", C, iters):
     var f2 {.noInit.}: Fp12[C]
-    f2.prod_xyz000_xyz000_into_abcdefghij00(l0, l1)
-    f.mul_sparse_by_abcdefghij00(f2)
-
-proc mulFp12_by_2lines_v1_xy000z_Bench*(C: static Curve, iters: int) =
-  var l0, l1: Line[Fp2[C]]
-  var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
-  let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
-  l0.line_double(T, P)
-  l1.line_double(T, P)
-  var f = rng.random_unsafe(Fp12[C])
-
-  bench("mulFp12 by 2 lines v1", C, iters):
-    f.mul_sparse_by_line_xy000z(l0)
-    f.mul_sparse_by_line_xy000z(l1)
-
-proc mulFp12_by_2lines_v2_xy000z_Bench*(C: static Curve, iters: int) =
-  var l0, l1: Line[Fp2[C]]
-  var T = rng.random_point(ECP_ShortW_Prj[Fp2[C], G2])
-  let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
-  l0.line_double(T, P)
-  l1.line_double(T, P)
-  var f = rng.random_unsafe(Fp12[C])
-
-  bench("mulFp12 by 2 lines v2", C, iters):
-    var f2 {.noInit.}: Fp12[C]
-    f2.prod_xy000z_xy000z_into_abcd00efghij(l0, l1)
-    f.mul_sparse_by_abcd00efghij(f2)
+    f2.prod_from_2_lines(l0, l1)
+    f.mul_by_prod_of_2_lines(f2)
 
 proc mulBench*(C: static Curve, iters: int) =
   var r: Fp12[C]
diff --git a/constantine/math/curves/bls12_377_constants.nim b/constantine/math/curves/bls12_377_constants.nim
index 9f95ea1..93b72de 100644
--- a/constantine/math/curves/bls12_377_constants.nim
+++ b/constantine/math/curves/bls12_377_constants.nim
@@ -8,7 +8,7 @@
 
 import
   ../config/curves,
-  ../io/io_extfields
+  ../io/[io_fields, io_extfields]
 
 # Curve precomputed parameters
 # -----------------------------------------------------------------
@@ -16,3 +16,7 @@ const BLS12_377_coefB_G2* = Fp2[BLS12_377].fromHex(
   "0x0",
   "0x10222f6db0fd6f343bd03737460c589dc7b4f91cd5fd889129207b63c6bf8000dd39e5c1ccccccd1c9ed9999999999a"
 )
+const BLS12_377_coefB_G2_times_3* = Fp2[BLS12_377].fromHex( 
+  "0x0",
+  "0x1582e9e796a73ef04fc0499f08107627b4f14c2672a760c18c2b4f2fb3aa000126f7dd026666666d0d3cccccccccccd"
+)
diff --git a/constantine/math/curves/bls12_381_constants.nim b/constantine/math/curves/bls12_381_constants.nim
index ab8ca64..8da486c 100644
--- a/constantine/math/curves/bls12_381_constants.nim
+++ b/constantine/math/curves/bls12_381_constants.nim
@@ -8,7 +8,7 @@
 
 import
   ../config/curves,
-  ../io/io_extfields
+  ../io/[io_fields, io_extfields]
 
 # Curve precomputed parameters
 # -----------------------------------------------------------------
@@ -16,3 +16,7 @@ const BLS12_381_coefB_G2* = Fp2[BLS12_381].fromHex(
   "0x4",
   "0x4"
 )
+const BLS12_381_coefB_G2_times_3* = Fp2[BLS12_381].fromHex( 
+  "0xc",
+  "0xc"
+)
diff --git a/constantine/math/curves/bn254_nogami_constants.nim b/constantine/math/curves/bn254_nogami_constants.nim
index 8c2ba96..dd087c6 100644
--- a/constantine/math/curves/bn254_nogami_constants.nim
+++ b/constantine/math/curves/bn254_nogami_constants.nim
@@ -8,7 +8,7 @@
 
 import
   ../config/curves,
-  ../io/io_extfields
+  ../io/[io_fields, io_extfields]
 
 # Curve precomputed parameters
 # -----------------------------------------------------------------
@@ -16,3 +16,7 @@ const BN254_Nogami_coefB_G2* = Fp2[BN254_Nogami].fromHex(
   "0x1",
   "0x2523648240000001ba344d80000000086121000000000013a700000000000012"
 )
+const BN254_Nogami_coefB_G2_times_3* = Fp2[BN254_Nogami].fromHex( 
+  "0x3",
+  "0x2523648240000001ba344d80000000086121000000000013a700000000000010"
+)
diff --git a/constantine/math/curves/bn254_snarks_constants.nim b/constantine/math/curves/bn254_snarks_constants.nim
index 0f53b95..12540aa 100644
--- a/constantine/math/curves/bn254_snarks_constants.nim
+++ b/constantine/math/curves/bn254_snarks_constants.nim
@@ -8,7 +8,7 @@
 
 import
   ../config/curves,
-  ../io/io_extfields
+  ../io/[io_fields, io_extfields]
 
 # Curve precomputed parameters
 # -----------------------------------------------------------------
@@ -16,3 +16,7 @@ const BN254_Snarks_coefB_G2* = Fp2[BN254_Snarks].fromHex(
   "0x2b149d40ceb8aaae81be18991be06ac3b5b4c5e559dbefa33267e6dc24a138e5",
   "0x9713b03af0fed4cd2cafadeed8fdf4a74fa084e52d1852e4a2bd0685c315d2"
 )
+const BN254_Snarks_coefB_G2_times_3* = Fp2[BN254_Snarks].fromHex( 
+  "0x20753adca9c6bfb81499be5e509e8f8ff21b7c8d3cb039cf1ef69c66bce9b021",
+  "0x1c53b10b0d2fc7e67860f09cc8af9ddf5eee18eaf8748f8ade8371391494176"
+)
diff --git a/constantine/math/curves/bw6_761_constants.nim b/constantine/math/curves/bw6_761_constants.nim
index 84fb31a..287b882 100644
--- a/constantine/math/curves/bw6_761_constants.nim
+++ b/constantine/math/curves/bw6_761_constants.nim
@@ -8,9 +8,11 @@
 
 import
   ../config/curves,
-  ../io/io_fields
+  ../io/[io_fields, io_extfields]
 
 # Curve precomputed parameters
 # -----------------------------------------------------------------
 const BW6_761_coefB_G2* = Fp[BW6_761].fromHex( 
   "0x4")
+const BW6_761_coefB_G2_times_3* = Fp[BW6_761].fromHex( 
+  "0xc")
diff --git a/constantine/math/curves/zoo_constants.nim b/constantine/math/curves/zoo_constants.nim
index a9c3933..8449070 100644
--- a/constantine/math/curves/zoo_constants.nim
+++ b/constantine/math/curves/zoo_constants.nim
@@ -25,3 +25,12 @@ macro getCoefB_G2*(C: static Curve): untyped =
   ##   y² = x³ + b*µ (M-Twist)
   ## with µ the non-residue (sextic non-residue with a sextic twist)
   return bindSym($C & "_coefB_G2")
+
+macro getCoefB_G2_times_3*(C: static Curve): untyped =
+  ## A pairing curve has the following equation on G1
+  ##   y² = x³ + b
+  ## and on G2
+  ##   y² = x³ + b/µ (D-Twist)
+  ##   y² = x³ + b*µ (M-Twist)
+  ## with µ the non-residue (sextic non-residue with a sextic twist)
+  return bindSym($C & "_coefB_G2_times_3")
\ No newline at end of file
diff --git a/constantine/math/extension_fields/towers.nim b/constantine/math/extension_fields/towers.nim
index 8c4ac22..92a718a 100644
--- a/constantine/math/extension_fields/towers.nim
+++ b/constantine/math/extension_fields/towers.nim
@@ -724,7 +724,7 @@ func has_large_NR_norm(C: static Curve): bool =
 
   return norm > 5
 
-func has_large_field_elem(C: static Curve): bool =
+func has_large_field_elem*(C: static Curve): bool =
   ## Returns true if field element are large
   ## and necessitate custom routine for assembly in particular
   let a = default(Fp[C])
@@ -1388,9 +1388,13 @@ func mul_sparse_by_0y*(a: var QuadraticExt, sparseB: auto) =
   ## Sparse in-place multiplication
   a.mul_sparse_by_0y(a, sparseB)
 
+func mul_sparse_by_x0*(r: var QuadraticExt, a: QuadraticExt, sparseB: auto) =
+  ## Sparse in-place multiplication
+  r.mul_sparse_generic_by_x0(a, sparseB)
+
 func mul_sparse_by_x0*(a: var QuadraticExt, sparseB: QuadraticExt) =
   ## Sparse in-place multiplication
-  a.mul_sparse_generic_by_x0(a, sparseB)
+  a.mul_sparse_by_x0(a, sparseB)
 
 template mulCheckSparse*(a: var QuadraticExt, b: QuadraticExt) =
   when b.isOne().bool:
diff --git a/constantine/math/pairing/lines_eval.nim b/constantine/math/pairing/lines_eval.nim
index f81a4d1..8c11535 100644
--- a/constantine/math/pairing/lines_eval.nim
+++ b/constantine/math/pairing/lines_eval.nim
@@ -16,7 +16,8 @@ import
     ec_shortweierstrass_affine,
     ec_shortweierstrass_projective
   ],
-  ../io/io_extfields
+  ../io/io_extfields,
+  ../curves/zoo_constants
 
 # No exceptions allowed
 {.push raises: [].}
@@ -38,7 +39,7 @@ import
 #
 # with z = SNR¹ᐟ⁶
 #
-# The cubic over quadatric towering
+# The cubic over quadratic towering
 # ---------------------------------
 #
 #   (a₀ + a₁ u) + (a₂ + a₃u) v + (a₄ + a₅u) v²
@@ -55,12 +56,13 @@ import
 # Mapping between towering schemes
 # --------------------------------
 #
-# c₀ <=> a₀ <=> b₀
-# c₁ <=> a₂ <=> b₃
-# c₂ <=> a₄ <=> b₁
-# c₃ <=> a₁ <=> b₄
-# c₄ <=> a₃ <=> b₂
-# c₅ <=> a₅ <=> b₅
+# canonical <=> cubic over quadratic <=> quadratic over cubic
+#    c₀     <=>        a₀            <=>            b₀
+#    c₁     <=>        a₂            <=>            b₃
+#    c₂     <=>        a₄            <=>            b₁
+#    c₃     <=>        a₁            <=>            b₄
+#    c₄     <=>        a₃            <=>            b₂
+#    c₅     <=>        a₅            <=>            b₅
 #
 # See also chapter 6.4
 # - Multiplication and Squaring on Pairing-Friendly Fields
@@ -69,7 +71,7 @@ import
 
 type
   Line*[F] = object
-    ## Packed line representation over a E'(Fpᵏ/d)
+    ## Packed line representation over a E'(𝔽pᵏ/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
@@ -80,16 +82,22 @@ type
     ## the non-zero coordinates depend on the twist kind.
     ##
     ## For a D-twist
-    ##   (x, y, z) corresponds to an sparse element of 𝔽p12
-    ##   with 𝔽p2 coordinates:
-    ##     - xyz000 (cubic over quadratic towering)
-    ##     - x00yz0 (quadratic over cubic towering)
+    ##   in canonical coordinates over sextic polynomial [1, w, w², w³, w⁴, w⁵]
+    ##   when evaluating the line at P(xₚ, yₚ)
+    ##     a.yₚ + b.xₚ w + c w³ 
+    ##     
+    ##   This translates in 𝔽pᵏ to
+    ##     - acb000 (cubic over quadratic towering)
+    ##     - a00bc0 (quadratic over cubic towering)
     ## For a M-Twist
-    ##   (x, y, z) corresponds to an sparse element of 𝔽p12
-    ##   with 𝔽p2 coordinates:
-    ##     - xy000z (cubic over quadratic towering)
-    ##     - xy00z0 (quadratic over cubic towering)
-    x*, y*, z*: F
+    ##   in canonical coordinates over sextic polynomial [1, w, w², w³, w⁴, w⁵]
+    ##   when evaluating the line at the twist ψ(P)(xₚw², yₚw³)
+    ##     a.yₚ w³ + b.xₚ w² + c 
+    ## 
+    ##   This translates in 𝔽pᵏ to
+    ##     - ca00b0 (cubic over quadratic towering)
+    ##     - cb00a0 (quadratic over cubic towering)
+    a*, b*, c*: F
 
   SexticNonResidue* = NonResidue
     ## The Sextic non-residue to build
@@ -124,8 +132,13 @@ func line_update[F1, F2](line: var Line[F2], P: ECP_ShortW_Aff[F1, G1]) =
   ## after addition or doubling
   ## P in G1
   static: doAssert F1.C == F2.C
-  line.x *= P.y
-  line.z *= P.x
+  # D-Twist: line at P(xₚ, yₚ):
+  #   a.yₚ + b.xₚ w + c w³
+  #
+  # M-Twist: line at ψ(P)(xₚw², yₚw³)
+  #   a.yₚ w³ + b.xₚ w² + c 
+  line.a *= P.y
+  line.b *= P.x
 
 # ############################################################
 #
@@ -149,189 +162,121 @@ func line_update[F1, F2](line: var Line[F2], P: ECP_ShortW_Aff[F1, G1]) =
 #   Patrick Longa, Shi Hu, and David Jao, 2012
 #   https://eprint.iacr.org/2012/408.pdf
 
-# Line evaluation only
+# Line evaluation
 # -----------------------------------------------------------------------------
 
-func line_eval_double[F](
-       line: var Line[F],
-       T: ECP_ShortW_Prj[F, G2]) =
-  ## Evaluate the line function for doubling
-  ## i.e. the tangent at T
-  ##
-  ## With T in homogenous projective coordinates (X, Y, Z)
-  ## And ξ the sextic non residue to construct the towering
-  ##
-  ## M-Twist:
-  ##   A = -2ξ Y.Z
-  ##   B = 3bξ Z² - Y²
-  ##   C = 3 X²
-  ##
-  ## D-Twist are scaled by ξ to avoid dividing by ξ:
-  ##   A = -2ξ Y.Z
-  ##   B = 3b Z² - ξY²
-  ##   C = 3ξ X²
-  ##
-  ## Instead of
-  ##   - equation 10 from The Real of pairing, Aranha et al, 2013
-  ##   - or chapter 3 from pairing Implementation Revisited, Scott 2019
-  ##   A = -2 Y.Z
-  ##   B = 3b/ξ Z² - Y²
-  ##   C = 3 X²
-  ##
-  ## A constant factor will be wiped by the final exponentiation
-  ## as for all non-zero α ∈ GF(pᵐ)
-  ## with
-  ## - p odd prime
-  ## - and gcd(α,pᵐ) = 1 (i.e. the extension field pᵐ is using irreducible polynomials)
-  ##
-  ## Little Fermat holds and we have
-  ## α^(pᵐ - 1) ≡ 1 (mod pᵐ)
-  ##
-  ## The final exponent is of the form
-  ## (pᵏ-1)/r
-  ##
-  ## A constant factor on twisted coordinates pᵏᐟᵈ
-  ## is a constant factor on pᵏ with d the twisting degree
-  ## and so will be elminated. QED.
-  var v {.noInit.}: F
-  const b3 = 3 * F.C.getCoefB()
-
-  template A: untyped = line.x
-  template B: untyped = line.y
-  template C: untyped = line.z
-
-  A.prod(T.y, T.z)      # A = Y.Z
-  C.square(T.x)         # C = X²
-  v.square(T.y)         # v = Y²
-  B.square(T.z)         # B = Z²
-
-  A.double()            # A =  2 Y.Z
-  A.neg()               # A = -2 Y.Z
-  A *= SexticNonResidue # A = -2 ξ Y.Z
-
-  B *= b3               # B = 3b Z²
-  C *= 3                # C = 3X²
-  when F.C.getSexticTwist() == M_Twist:
-    B *= SexticNonResidue # B = 3b' Z² = 3bξ Z²
-  elif F.C.getSexticTwist() == D_Twist:
-    v *= SexticNonResidue # v =  ξ Y²
-    C *= SexticNonResidue # C = 3ξ X²
-  else:
-    {.error: "unreachable".}
-
-  B -= v                # B = 3bξ Z² - Y²  (M-twist)
-                        # B = 3b Z² - ξ Y² (D-twist)
-
-func line_eval_add[F](
-       line: var Line[F],
-       T: ECP_ShortW_Prj[F, G2],
-       Q: ECP_ShortW_Aff[F, G2]) =
-  ## Evaluate the line function for addition
-  ## i.e. the line between T and Q
-  ##
-  ## With T in homogenous projective coordinates (X, Y, Z)
-  ## And ξ the sextic non residue to construct the towering
-  ##
-  ## M-Twist:
-  ##   A = ξ (X₁ - Z₁X₂)
-  ##   B = (Y₁ - Z₁Y₂) X₂ - (X₁ - Z₁X₂) Y₂
-  ##   C = - (Y₁ - Z₁Y₂)
-  ##
-  ## D-Twist:
-  ##   A = X₁ - Z₁X₂
-  ##   B = (Y₁ - Z₁Y₂) X₂ - (X₁ - Z₁X₂) Y₂
-  ##   C = - (Y₁ - Z₁Y₂)
-  ##
-  ## Note: There is no need for complete formula as
-  ## we have T ∉ [Q, -Q] in the Miller loop doubling-and-add
-  ## i.e. the line cannot be vertical
-  var v {.noInit.}: F
-
-  template A: untyped = line.x
-  template B: untyped = line.y
-  template C: untyped = line.z
-
-  v.prod(T.z, Q.y) # v = Z₁Y₂
-  B.prod(T.z, Q.x) # B = Z₁X₂
-
-  A.diff(T.x, B)   # A = X₁-Z₁X₂
-  C.diff(T.y, v)   # C = Y₁-Z₁Y₂
-
-  v.prod(A, Q.y)   # v = (X₁-Z₁X₂) Y₂
-  B.prod(C, Q.x)   # B = (Y₁-Z₁Y₂) X₂
-  B -= v           # B = (Y₁-Z₁Y₂) X₂ - (X₁-Z₁X₂) Y₂
-
-  C.neg()          # C = -(Y₁-Z₁Y₂)
-
-  when F.C.getSexticTwist() == M_Twist:
-    A *= SexticNonResidue # A = ξ (X₁ - Z₁X₂)
+# Line for a doubling
+# ===================
+#
+# With T in homogenous projective coordinates (X, Y, Z)
+# And ξ the sextic non residue to construct the towering
+#
+# M-Twist:
+#   A = -2ξ Y.Z      [w³]
+#   B = 3 X²         [w²]
+#   C = 3bξ Z² - Y²  [1]
+#
+# D-Twist may be scaled by ξ to avoid dividing by ξ:
+#   A = -2ξ Y.Z      [1]
+#   C = 3ξ X²        [w]
+#   B = 3b Z² - ξY²  [w³]
+#
+# Instead of
+#   - equation 10 from The Real of pairing, Aranha et al, 2013
+#   - or chapter 3 from pairing Implementation Revisited, Scott 2019
+#   A = -2 Y.Z
+#   B = 3 X²
+#   C = 3b/ξ Z² - Y²
+#
+# Note: This tradeoff a division with 3 multiplication by a non-residue.
+#       This is interesting for ξ has a small norm, but
+#       BN254_Snarks for example is 9+𝑖
+#
+# A constant factor will be wiped by the final exponentiation
+# as for all non-zero α ∈ GF(pᵐ)
+# with
+# - p odd prime
+# - and gcd(α,pᵐ) = 1 (i.e. the extension field pᵐ is using irreducible polynomials)
+#
+# Little Fermat holds and we have
+# α^(pᵐ - 1) ≡ 1 (mod pᵐ)
+#
+# The final exponent is of the form
+# (pᵏ-1)/r
+#
+# A constant factor on twisted coordinates pᵏᐟᵈ
+# is a constant factor on pᵏ with d the twisting degree
+# and so will be elminated. QED.
+#
+# Line for an addition
+# ====================
+#
+# With T in homogenous projective coordinates (X, Y, Z)
+# And ξ the sextic non residue to construct the towering
+#
+# M-Twist:
+#   A = X₁ - Z₁X₂
+#   B = - (Y₁ - Z₁Y₂)
+#   C = (Y₁ - Z₁Y₂) X₂ - (X₁ - Z₁X₂) Y₂
+#
+# D-Twist:
+#   A = X₁ - Z₁X₂
+#   B = - (Y₁ - Z₁Y₂)
+#   C = (Y₁ - Z₁Y₂) X₂ - (X₁ - Z₁X₂) Y₂
+#
+# Note: There is no need for complete formula as
+# we have T ∉ [Q, -Q] in the Miller loop doubling-and-add
+# i.e. the line cannot be vertical
 
 func line_eval_fused_double[Field](
        line: var Line[Field],
        T: var ECP_ShortW_Prj[Field, G2]) =
   ## Fused line evaluation and elliptic point doubling
   # Grewal et al, 2012 adapted to Scott 2019 line notation
+
   var A {.noInit.}, B {.noInit.}, C {.noInit.}: Field
   var E {.noInit.}, F {.noInit.}, G {.noInit.}: Field
-  template H: untyped = line.x
-  const b3 = 3*Field.C.getCoefB()
-
-  var snrY = T.y
-  when Field.C.getSexticTwist() == D_Twist:
-    snrY *= SexticNonResidue
-
-  A.prod(T.x, snrY)
+  
+  template H: untyped = line.a
+  template I: untyped = line.b
+  template J: untyped = line.c
+  
+  A.prod(T.x, T.y)
   A.div2()          # A = XY/2
   B.square(T.y)     # B = Y²
   C.square(T.z)     # C = Z²
 
-  var snrB = B
-  when Field.C.getSexticTwist() == D_Twist:
-    snrB *= SexticNonResidue
+  # E = 3b'Z² = 3bξ Z² (M-Twist) or 3b/ξ Z² (D-Twist)
+  when Field.C.getSexticTwist() == M_Twist and E.fromComplexExtension():
+    const b3 = 3*Field.C.getCoefB()
+    E.prod(C, b3)
+    E *= SexticNonResidue
+  else:
+    E = C
+    E.mulCheckSparse(Field.C.getCoefB_G2_times_3())
 
-  E.prod(C, b3)
-  when Field.C.getSexticTwist() == M_Twist:
-    E *= SexticNonResidue # E = 3b'Z² = 3bξ Z²
-
-  F.prod(E, 3)      # F = 3E = 9bZ²
-  G.sum(snrB, F)
+  F.prod(E, 3)      # F = 3E = 9b'Z²
+  G.sum(B, 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
+  H -= C            # H = (Y+Z)²-(B+C)= 2YZ
 
-  line.z.square(T.x)
-  line.z *= 3       # lz = 3X²
-  when Field.C.getSexticTwist() == D_Twist:
-    line.z *= SexticNonResidue
+  I.square(T.x)
+  I *= 3            # I = 3X²
 
-  line.y.diff(E, snrB) # ly = E-B = 3b'Z² - Y²
+  J.diff(E, B)      # J = E-B = 3b'Z² - Y²
 
   # In-place modification: invalidates `T.` calls
-  T.x.diff(snrB, F)
+  T.x.diff(B, 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 Field.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 Field.C.getSexticTwist() == M_Twist:
-    H *= SexticNonResidue
-    # else: the SNR is already integrated in H
+  T.z.prod(B, H)    # Z₃ = BH = Y²((Y+Z)² - (Y²+Z²)) = 2Y³Z
+  H.neg()
 
 func line_eval_fused_add[Field](
        line: var Line[Field],
@@ -350,13 +295,13 @@ func line_eval_fused_add[Field](
     H {.noInit.}: Field
     I {.noInit.}: Field
 
-  template lambda: untyped = line.x
-  template theta: untyped = line.z
-  template J: untyped = line.y
+  template lambda: untyped = line.a
+  template theta: untyped = line.b
+  template J: untyped = line.c
 
   A.prod(Q.y, T.z)
   B.prod(Q.x, T.z)
-  theta.diff(T.y, A)  # θ = Y₁ - Z₁X₂
+  theta.diff(T.y, A)  # θ = Y₁ - Z₁Y₂
   lambda.diff(T.x, B) # λ = X₁ - Z₁X₂
   C.square(theta)
   D.square(lambda)
@@ -383,8 +328,6 @@ func line_eval_fused_add[Field](
 
   # Line evaluation
   theta.neg()
-  when Field.C.getSexticTwist() == M_Twist:
-    lambda *= SexticNonResidue # A = ξ (X₁ - Z₁X₂)
 
 # Public line evaluation procedures
 # -----------------------------------------------------------------------------
@@ -398,13 +341,8 @@ func line_double*[F1, F2](
   ##
   ## Compute lt,t(P)
   static: doAssert F1.C == F2.C
-  when true:
-    line_eval_fused_double(line, T)
-    line.line_update(P)
-  else:
-    line_eval_double(line, T)
-    line.line_update(P)
-    T.double()
+  line_eval_fused_double(line, T)
+  line.line_update(P)
 
 func line_add*[F1, F2](
        line: var Line[F2],
@@ -416,14 +354,8 @@ func line_add*[F1, F2](
   ##
   ## Compute lt,q(P)
   static: doAssert F1.C == F2.C
-  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
-
+  line_eval_fused_add(line, T, Q)
+  line.line_update(P)
 
 # ############################################################
 #
@@ -461,7 +393,7 @@ func line_add*[F1, F2](
 func mul_by_line_xy0*[Fpkdiv2, Fpkdiv6](
        r: var Fpkdiv2,
        a: Fpkdiv2,
-       b: Line[Fpkdiv6]) =
+       x, y: Fpkdiv6) =
   ## Sparse multiplication of an 𝔽pᵏᐟ² element
   ## with coordinates (a₀, a₁, a₂) by a line (x, y, 0)
   ## The z coordinates in the line will be ignored.
@@ -475,26 +407,27 @@ func mul_by_line_xy0*[Fpkdiv2, Fpkdiv6](
     v0 {.noInit.}: Fpkdiv6
     v1 {.noInit.}: Fpkdiv6
 
-  v0.prod(a.c0, b.x)
-  v1.prod(a.c1, b.y)
-  r.c0.prod(a.c2, b.y)
+  v0.prod(a.c0, x)
+  v1.prod(a.c1, y)
+  r.c0.prod(a.c2, y)
   r.c0 *= SexticNonResidue
   r.c0 += v0
 
   r.c1.sum(a.c0, a.c1) # Error when r and a alias as r.c0 was updated
-  r.c2.sum(b.x, b.y)
+  r.c2.sum(x, y)
   r.c1 *= r.c2
   r.c1 -= v0
   r.c1 -= v1
 
-  r.c2.prod(a.c2, b.x)
+  r.c2.prod(a.c2, x)
   r.c2 += v1
 
-func mul_sparse_by_line_xy00z0*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
+func mul_sparse_by_line_ab00c0*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
   ## Sparse multiplication of an 𝔽pᵏ element
-  ## by a sparse 𝔽pᵏ element coming from an D-Twist line function.
+  ## by a sparse 𝔽pᵏ element coming from an D-Twist line function
+  ## With a quadratic over cubic towering (Fp2 -> Fp6 -> Fp12)
   ## The sparse element is represented by a packed Line type
-  ## with coordinate (x,y,z) matching 𝔽pᵏ coordinates xy00z0
+  ## with coordinate (a,b,c) matching 𝔽pᵏ coordinates ab00c0
 
   static:
     doAssert Fpk.C.getSexticTwist() == D_Twist
@@ -509,13 +442,13 @@ func mul_sparse_by_line_xy00z0*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
     v2 {.noInit.}: Line[Fpkdiv6]
     v3 {.noInit.}: Fpkdiv2
 
-  v0.mul_by_line_xy0(f.c0, l)
-  v1.mul_sparse_by_0y0(f.c1, l.z)
+  v0.mul_by_line_xy0(f.c0, l.a, l.b)
+  v1.mul_sparse_by_0y0(f.c1, l.c)
 
-  v2.x = l.x
-  v2.y.sum(l.y, l.z)
+  v2.x = l.a
+  v2.y.sum(l.b, l.c)
   f.c1 += f.c0
-  v3.mul_by_line_xy0(f.c1, v2)
+  v3.mul_by_line_xy0(f.c1, v2.x, v2.y)
   v3 -= v0
   f.c1.diff(v3, v1)
 
@@ -534,11 +467,12 @@ func mul_sparse_by_line_xy00z0*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
 # D-Twist
 # ------------------------------------------------------------
 
-func mul_sparse_by_line_xyz000*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
+func mul_sparse_by_line_acb000*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
   ## Sparse multiplication of an 𝔽pᵏ element
   ## by a sparse 𝔽pᵏ element coming from an D-Twist line function.
+  ## with a cubic over quadratic towering (Fp2 -> Fp4 -> Fp12)
   ## The sparse element is represented by a packed Line type
-  ## with coordinates (x,y,z) matching 𝔽pᵏ coordinates xyz000
+  ## with coordinate (a,b,c) matching 𝔽pᵏ coordinates acb000
 
   static:
     doAssert Fpk.C.getSexticTwist() == D_Twist
@@ -549,12 +483,12 @@ func mul_sparse_by_line_xyz000*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
 
   # In the following equations (taken from cubic extension implementation)
   # a = f
-  # b0 = (x, y)
-  # b1 = (z, 0)
+  # b0 = (a, c)
+  # b1 = (b, 0)
   # b2 = (0, 0)
   #
-  # v0 = a0 b0 = (f00, f01).(x, y)
-  # v1 = a1 b1 = (f10, f11).(z, 0)
+  # v0 = a0 b0 = (f00, f01).(a, c)
+  # v1 = a1 b1 = (f10, f11).(b, 0)
   # v2 = a2 b2 = (f20, f21).(0, 0)
   #
   # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
@@ -566,25 +500,25 @@ func mul_sparse_by_line_xyz000*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
   #    = a0 b0 + a2 b0 - v0 + v1
   #    = a2 b0 + v1
 
-  when false:
-    var b0 {.noInit.}, v0{.noInit.}, v1{.noInit.}, t{.noInit.}: Fp4[C]
+  when Fpk.C.has_large_field_elem():
+    var b0 {.noInit.}, v0{.noInit.}, v1{.noInit.}, t{.noInit.}: Fpkdiv3
 
-    b0.c0 = l.x
-    b0.c1 = l.y
+    b0.c0 = l.a
+    b0.c1 = l.c
 
     v0.prod(f.c0, b0)
-    v1.mul_sparse_by_x0(f.c1, l.z)
+    v1.mul_sparse_by_x0(f.c1, l.b)
 
     # r1 = (a0 + a1) * (b0 + b1) - v0 - v1
-    f.c1 += f.c0 # r1 = a0 + a1
+    f.c1 += f.c0  # r1 = a0 + a1
     t = b0
-    t.c0 += l.z  # t = b0 + b1
-    f.c1 *= t    # r2 = (a0 + a1)(b0 + b1)
+    t.c0 += l.b   # t = b0 + b1
+    f.c1 *= t     # r2 = (a0 + a1)(b0 + b1)
     f.c1 -= v0
-    f.c1 -= v1   # r2 = (a0 + a1)(b0 + b1) - v0 - v1
+    f.c1 -= v1    # r2 = (a0 + a1)(b0 + b1) - v0 - v1
 
     # r0 = ξ a2 b1 + v0
-    f.c0.mul_sparse_by_x0(f.c2, l.z)
+    f.c0.mul_sparse_by_x0(f.c2, l.b)
     f.c0 *= SexticNonResidue
     f.c0 += v0
 
@@ -596,41 +530,41 @@ func mul_sparse_by_line_xyz000*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
     var V0{.noInit.}, V1{.noInit.}, f2x{.noInit.}: doublePrec(Fpkdiv3)
     var t{.noInit.}: Fpkdiv6
 
-    V0.prod2x_disjoint(f.c0, l.x, l.y)
-    V1.mul2x_sparse_by_x0(f.c1, l.z)
+    V0.prod2x_disjoint(f.c0, l.a, l.c)
+    V1.mul2x_sparse_by_x0(f.c1, l.b)
 
     # r1 = (a0 + a1) * (b0 + b1) - v0 - v1
     f.c1.sum(f.c1, f.c0)
-    t.sum(l.x, l.z)                 # b0 is (x, y)
-    f2x.prod2x_disjoint(f.c1, t, l.y) # b1 is (z, 0)
+    t.sum(l.a, l.b)                  # b0 is (x, y)
+    f2x.prod2x_disjoint(f.c1, t, l.c) # b1 is (z, 0)
     f2x.diff2xMod(f2x, V0)
     f2x.diff2xMod(f2x, V1)
     f.c1.redc2x(f2x)
 
     # r0 = ξ a2 b1 + v0
-    f2x.mul2x_sparse_by_x0(f.c2, l.z)
+    f2x.mul2x_sparse_by_x0(f.c2, l.b)
     f2x.prod2x(f2x, SexticNonResidue)
     f2x.sum2xMod(f2x, V0)
     f.c0.redc2x(f2x)
 
     # r2 = a2 b0 + v1
-    f2x.prod2x_disjoint(f.c2, l.x, l.y)
+    f2x.prod2x_disjoint(f.c2, l.a, l.c)
     f2x.sum2xMod(f2x, V1)
     f.c2.redc2x(f2x)
 
-func prod_xyz000_xyz000_into_abcdefghij00*[Fpk, Fpkdiv6](f: var Fpk, l0, l1: Line[Fpkdiv6]) =
+func prod_xzy000_xzy000_into_abcdefghij00*[Fpk, Fpkdiv6](f: var Fpk, l0, l1: Line[Fpkdiv6]) =
   ## Multiply 2 lines together
   ## The result is sparse in f.c1.c1
   # In the following equations (taken from cubic extension implementation)
-  # a0 = (x0, y0)
-  # a1 = (z0,  0)
+  # a0 = (x0, z0)
+  # a1 = (y0,  0)
   # a2 = ( 0,  0)
-  # b0 = (x1, y1)
-  # b1 = (z1,  0)
+  # b0 = (x1, z1)
+  # b1 = (y1,  0)
   # b2 = ( 0,  0)
   #
-  # v0 = a0 b0 = (x0, y0).(x1, y1)
-  # v1 = a1 b1 = (z0,  0).(z1,  0)
+  # v0 = a0 b0 = (x0, z0).(x1, z1)
+  # v1 = a1 b1 = (y0,  0).(y1,  0)
   # v2 = a2 b2 = ( 0,  0).( 0,  0)
   #
   # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
@@ -652,13 +586,13 @@ func prod_xyz000_xyz000_into_abcdefghij00*[Fpk, Fpkdiv6](f: var Fpk, l0, l1: Lin
   var V0{.noInit.}, f2x{.noInit.}: doublePrec(Fpkdiv3)
   var V1{.noInit.}: doublePrec(Fpkdiv6)
 
-  V0.prod2x_disjoint(l0.x, l0.y, l1.x, l1.y) # a0 b0 = (x0, y0).(x1, y1)
-  V1.prod2x(l0.z, l1.z)                      # a1 b1 = (z0,  0).(z1,  0)
+  V0.prod2x_disjoint(l0.a, l0.c, l1.a, l1.c) # a0 b0 = (x0, z0).(x1, z1)
+  V1.prod2x(l0.b, l1.b)                      # a1 b1 = (y0,  0).(y1,  0)
 
   # r1 = (a0 + a1) * (b0 + b1) - v0 - v1
-  f.c1.c0.sum(l0.x, l0.z)                           # x0 + z0
-  f.c1.c1.sum(l1.x, l1.z)                           # x1 + z1
-  f2x.prod2x_disjoint(f.c1.c0, l0.y, f.c1.c1, l1.y) # (x0 + z0, y0)(x1 + z1, y1) = (a0 + a1) * (b0 + b1)
+  f.c1.c0.sum(l0.a, l0.b)                           # x0 + y0
+  f.c1.c1.sum(l1.a, l1.b)                           # x1 + y1
+  f2x.prod2x_disjoint(f.c1.c0, l0.c, f.c1.c1, l1.c) # (x0 + y0, z0)(x1 + y1, z1) = (a0 + a1) * (b0 + b1)
   f2x.diff2xMod(f2x, V0)
   f2x.c0.diff2xMod(f2x.c0, V1)
   f.c1.redc2x(f2x)
@@ -744,25 +678,29 @@ func mul_sparse_by_abcdefghij00*[Fpk](
 # M-Twist
 # ------------------------------------------------------------
 
-func mul_sparse_by_line_xy000z*[Fpk, Fpkdiv6](
-       f: var Fpk, l: Line[Fpkdiv6]) =
+func mul_sparse_by_line_ca00b0*[Fpk, Fpkdiv6](f: var Fpk, l: Line[Fpkdiv6]) =
+  ## Sparse multiplication of an 𝔽pᵏ element
+  ## by a sparse 𝔽pᵏ element coming from an M-Twist line function
+  ## with a cubic over quadratic towering (Fp2 -> Fp4 -> Fp12)
+  ## The sparse element is represented by a packed Line type
+  ## with coordinate (a,b,c) matching 𝔽pᵏ coordinates ca00b0
 
   static:
     doAssert Fpk.C.getSexticTwist() == M_Twist
     doAssert f is CubicExt, "This assumes 𝔽pᵏ as a cubic extension of 𝔽pᵏᐟ³"
-    doAssert f.c0 is QuadraticExt, "This assumes 𝔽pᵏᐟ³ as a quadratic extension of 𝔽pᵏᐟ⁶"
+    doAssert f.c0 is QuadraticExt, "This assumes 𝔽pᵏᐟ³ as a cubic extension of 𝔽pᵏᐟ⁶"
 
   type Fpkdiv3 = typeof(f.c0)
 
   # In the following equations (taken from cubic extension implementation)
   # a = f
-  # b0 = (x, y)
+  # b0 = (c, a)
   # b1 = (0, 0)
-  # b2 = (0, z)
+  # b2 = (b, 0)
   #
-  # v0 = a0 b0 = (f00, f01).(x, y)
+  # v0 = a0 b0 = (f00, f01).(c, a)
   # v1 = a1 b1 = (f10, f11).(0, 0)
-  # v2 = a2 b2 = (f20, f21).(0, z)
+  # v2 = a2 b2 = (f20, f21).(b, 0)
   #
   # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
   #    = ξ (a1 b2 + a2 b2 - v2) + v0
@@ -773,25 +711,25 @@ func mul_sparse_by_line_xy000z*[Fpk, Fpkdiv6](
   # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1
   #    = (a0 + a2) * (b0 + b2) - v0 - v2
 
-  when false:
+  when Fpk.C.has_large_field_elem():
     var b0 {.noInit.}, v0{.noInit.}, v2{.noInit.}, t{.noInit.}: Fpkdiv3
 
-    b0.c0 = l.x
-    b0.c1 = l.y
+    b0.c0 = l.c
+    b0.c1 = l.a
 
     v0.prod(f.c0, b0)
-    v2.mul_sparse_by_0y(f.c2, l.z)
+    v2.mul_sparse_by_x0(f.c2, l.b)
 
     # r2 = (a0 + a2) * (b0 + b2) - v0 - v2
     f.c2 += f.c0 # r2 = a0 + a2
     t = b0
-    t.c1 += l.z  # t = b0 + b2
+    t.c0 += l.b # t = b0 + b2
     f.c2 *= t    # r2 = (a0 + a2)(b0 + b2)
     f.c2 -= v0
     f.c2 -= v2   # r2 = (a0 + a2)(b0 + b2) - v0 - v2
 
     # r0 = ξ a1 b2 + v0
-    f.c0.mul_sparse_by_0y(f.c1, l.z)
+    f.c0.mul_sparse_by_x0(f.c1, l.b)
     f.c0 *= SexticNonResidue
     f.c0 += v0
 
@@ -804,43 +742,43 @@ func mul_sparse_by_line_xy000z*[Fpk, Fpkdiv6](
     var V0{.noInit.}, V2{.noInit.}, f2x{.noInit.}: doublePrec(Fpkdiv3)
     var t{.noInit.}: Fpkdiv6
 
-    V0.prod2x_disjoint(f.c0, l.x, l.y)
-    V2.mul2x_sparse_by_0y(f.c2, l.z)
+    V0.prod2x_disjoint(f.c0, l.c, l.a)
+    V2.mul2x_sparse_by_x0(f.c2, l.b)
 
     # r2 = (a0 + a2) * (b0 + b2) - v0 - v2
     f.c2.sum(f.c2, f.c0)
-    t.sum(l.y, l.z)                 # b0 is (x, y)
-    f2x.prod2x_disjoint(f.c2, l.x, t) # b2 is (0, z)
+    t.sum(l.c, l.b)                  # b0 + b2 = (c+b, a)
+    f2x.prod2x_disjoint(f.c2, t, l.a)
     f2x.diff2xMod(f2x, V0)
     f2x.diff2xMod(f2x, V2)
     f.c2.redc2x(f2x)
 
     # r0 = ξ a1 b2 + v0
-    f2x.mul2x_sparse_by_0y(f.c1, l.z)
+    f2x.mul2x_sparse_by_x0(f.c1, l.b)
     f2x.prod2x(f2x, SexticNonResidue)
     f2x.sum2xMod(f2x, V0)
     f.c0.redc2x(f2x)
 
     # r1 = a1 b0 + ξ v2
-    f2x.prod2x_disjoint(f.c1, l.x, l.y)
+    f2x.prod2x_disjoint(f.c1, l.c, l.a)
     V2.prod2x(V2, SexticNonResidue)
     f2x.sum2xMod(f2x, V2)
     f.c1.redc2x(f2x)
 
-func prod_xy000z_xy000z_into_abcd00efghij*[Fpk, Fpkdiv6](f: var Fpk, l0, l1: Line[Fpkdiv6]) =
+func prod_zx00y0_zx00y0_into_abcd00efghij*[Fpk, Fpkdiv6](f: var Fpk, l0, l1: Line[Fpkdiv6]) =
   ## Multiply 2 lines together
-  ## The result is sparse in f.c1.c0
+  ## The result is sparse in f.c1.c1
   # In the following equations (taken from cubic extension implementation)
-  # a0 = (x0, y0)
+  # a0 = (z0, x0)
   # a1 = ( 0,  0)
-  # a2 = ( 0, z0)
-  # b0 = (x1, y1)
+  # a2 = (y0,  0)
+  # b0 = (z1, x1)
   # b1 = ( 0,  0)
-  # b2 = ( 0, z1)
+  # b2 = (y0,  0)
   #
-  # v0 = a0 b0 = (x0, y0).(x1, y1)
+  # v0 = a0 b0 = (z0, x0).(z1, x1)
   # v1 = a1 b1 = ( 0,  0).( 0,  0)
-  # v2 = a2 b2 = ( 0, z0).( 0, z1)
+  # v2 = a2 b2 = (y0,  0).(y1,  0)
   #
   # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0
   #    = ξ (a1 b2 + a2 b2 - v2) + v0
@@ -861,14 +799,13 @@ func prod_xy000z_xy000z_into_abcd00efghij*[Fpk, Fpkdiv6](f: var Fpk, l0, l1: Lin
   var V0{.noInit.}, f2x{.noInit.}: doublePrec(Fpkdiv3)
   var V2{.noInit.}: doublePrec(Fpkdiv6)
 
-  V0.prod2x_disjoint(l0.x, l0.y, l1.x, l1.y) # a0 b0 = (x0, y0).(x1, y1)
-  V2.prod2x(l0.z, l1.z)                      # a2 b2 = ( 0, z0).( 0, z1)
-  V2.prod2x(V2, NonResidue)
+  V0.prod2x_disjoint(l0.c, l0.a, l1.c, l1.a) # a0 b0 = (z0, x0).(z1, x1)
+  V2.prod2x(l0.b, l1.b)                      # a2 b2 = (y0,  0).(y1,  0)
 
   # r2 = (a0 + a2) * (b0 + b2) - v0 - v2
-  f.c2.c0.sum(l0.y, l0.z)                           # y0 + z0
-  f.c2.c1.sum(l1.y, l1.z)                           # y1 + z1
-  f2x.prod2x_disjoint(l0.x, f.c2.c0, l1.x, f.c2.c1) # (x0, y0 + z0).(x1, y1 + z1) = (a0 + a2) * (b0 + b2)
+  f.c2.c0.sum(l0.b, l0.c)                           # y0 + z0
+  f.c2.c1.sum(l1.b, l1.c)                           # y1 + z1
+  f2x.prod2x_disjoint(f.c2.c0, l0.a, f.c2.c1, l1.a) # (z0 + y0, x0).(z1 + y1, x1) = (a0 + a2) * (b0 + b2)
   f2x.diff2xMod(f2x, V0)                            # (a0 + a2) * (b0 + b2) - v0
   f2x.c0.diff2xMod(f2x.c0, V2)                      # (a0 + a2) * (b0 + b2) - v0 - v2
   f.c2.redc2x(f2x)
@@ -955,34 +892,37 @@ func mul_sparse_by_abcd00efghij*[Fpk](
 # ------------------------------------------------------------
 
 func mul_by_line*[Fpk, Fpkdiv6](f: var Fpk, line: Line[Fpkdiv6]) {.inline.} =
-  ## Multiply an element of Fp12 by a sparse line function (xyz000 or xy000z)
+  ## Multiply an element of Fp12 by a sparse line function
   when Fpk.C.getSexticTwist() == D_Twist:
-    f.mul_sparse_by_line_xyz000(line)
+    f.mul_sparse_by_line_acb000(line)
   elif Fpk.C.getSexticTwist() == M_Twist:
-    f.mul_sparse_by_line_xy000z(line)
+    f.mul_sparse_by_line_ca00b0(line)
   else:
     {.error: "A line function assumes that the curve has a twist".}
 
 func prod_from_2_lines*[Fpk, Fpkdiv6](f: var Fpk, line0, line1: Line[Fpkdiv6]) {.inline.} =
-  ## Multiply 2 lines function (xyz000 or xy000z)
+  ## Multiply 2 lines function
   ## and store the result in f
   ## f is overwritten
   when Fpk.C.getSexticTwist() == D_Twist:
-    f.prod_xyz000_xyz000_into_abcdefghij00(line0, line1)
+    f.prod_xzy000_xzy000_into_abcdefghij00(line0, line1)
   elif Fpk.C.getSexticTwist() == M_Twist:
-    f.prod_xy000z_xy000z_into_abcd00efghij(line0, line1)
+    f.prod_zx00y0_zx00y0_into_abcd00efghij(line0, line1)
+  else:
+    {.error: "A line function assumes that the curve has a twist".}
+
+func mul_by_prod_of_2_lines*[Fpk](f: var Fpk, g: Fpk) {.inline.} =
+  ## Multiply f by the somewhat sparse product of 2 lines
+  when Fpk.C.getSexticTwist() == D_Twist:
+    f.mul_sparse_by_abcdefghij00(g)
+  elif Fpk.C.getSexticTwist() == M_Twist:
+    f.mul_sparse_by_abcd00efghij(g)
   else:
     {.error: "A line function assumes that the curve has a twist".}
 
 func mul_by_2_lines*[Fpk, Fpkdiv6](f: var Fpk, line0, line1: Line[Fpkdiv6]) {.inline.} =
-  ## Multiply f*line0*line1 with lines (xyz000 or xy000z)
+  ## Multiply f*line0*line1 with lines
   ## f is updated with the result
-  var t{.noInit.}: typeof(f)
-  when Fpk.C.getSexticTwist() == D_Twist:
-    t.prod_xyz000_xyz000_into_abcdefghij00(line0, line1)
-    f.mul_sparse_by_abcdefghij00(t)
-  elif Fpk.C.getSexticTwist() == M_Twist:
-    t.prod_xy000z_xy000z_into_abcd00efghij(line0, line1)
-    f.mul_sparse_by_abcd00efghij(t)
-  else:
-    {.error: "A line function assumes that the curve has a twist".}
+  var t{.noInit.}: Fpk
+  t.prod_from_2_lines(line0, line1)
+  f.mul_by_prod_of_2_lines(t)
diff --git a/sage/g2_params.sage b/sage/g2_params.sage
index 63f4201..dc717c9 100644
--- a/sage/g2_params.sage
+++ b/sage/g2_params.sage
@@ -129,6 +129,10 @@ def gen_coef_b_on_G2(curve_name, curve_config):
   buf += field_to_nim(G2B, G2_field, curve_name)
   buf += '\n'
 
+  buf += f'const {curve_name}_coefB_G2_times_3* = '
+  buf += field_to_nim(3*G2B, G2_field, curve_name)
+  buf += '\n'
+
   return buf
 
 # CLI
@@ -156,17 +160,17 @@ if __name__ == "__main__":
   else:
     G2B = gen_coef_b_on_G2(curve, Curves)
 
-    with open(f'{curve.lower()}_precomputed_params.nim', 'w') as f:
+    with open(f'{curve.lower()}_constants.nim', 'w') as f:
       f.write(copyright())
       f.write('\n\n')
 
       f.write(inspect.cleandoc("""
           import
             ../config/curves,
-            ../io/io_extfields
+            ../io/[io_fields, io_extfields]
       """))
 
       f.write('\n\n')
       f.write(G2B)
 
-    print(f'Successfully created {curve.lower()}_precomputed_params.nim')
+    print(f'Successfully created {curve.lower()}_constants.nim')
diff --git a/tests/math/t_pairing_mul_fp12_by_lines.nim b/tests/math/t_pairing_mul_fp12_by_lines.nim
index 7d2e142..61eec8d 100644
--- a/tests/math/t_pairing_mul_fp12_by_lines.nim
+++ b/tests/math/t_pairing_mul_fp12_by_lines.nim
@@ -94,12 +94,11 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
         let x = rng.random_elem(Fp2[C], gen)
         let y = rng.random_elem(Fp2[C], gen)
         let b = Fp6[C](coords: [x, y, Fp2[C]()])
-        let line = Line[Fp2[C]](x: x, y: y)
 
         var r {.noInit.}, r2 {.noInit.}: Fp6[C]
 
         r.prod(a, b)
-        r2.mul_by_line_xy0(a, line)
+        r2.mul_by_line_xy0(a, x, y)
 
         check: bool(r == r2)
 
@@ -109,8 +108,8 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
       test_fp6_xy0(curve, gen = Long01Sequence)
 
   when Fp12[BN254_Snarks]().c0.typeof is Fp6:
-    test "Sparse 𝔽p12/𝔽p6 resulting from xy00z0 line function":
-      proc test_fp12_xy00z0(C: static Curve, gen: static RandomGen) =
+    test "Sparse 𝔽p12/𝔽p6 resulting from ab00c0 line function":
+      proc test_fp12_ab00c0(C: static Curve, gen: static RandomGen) =
         for _ in 0 ..< Iters:
           var a = rng.random_elem(Fp12[C], gen)
           var a2 = a
@@ -119,14 +118,14 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
           var y = rng.random_elem(Fp2[C], gen)
           var z = rng.random_elem(Fp2[C], gen)
 
-          let line = Line[Fp2[C]](x: x, y: y, z: z)
+          let line = Line[Fp2[C]](a: x, b: y, c: z)
           let b = Fp12[C](
             c0: Fp6[C](coords: [       x, y, Fp2[C]()]),
             c1: Fp6[C](coords: [Fp2[C](), z, Fp2[C]()])
           )
 
           a *= b
-          a2.mul_sparse_by_line_xy00z0(line)
+          a2.mul_sparse_by_line_ab00c0(line)
 
           check: bool(a == a2)
 
@@ -135,8 +134,8 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
         test_fp12_xy00z0(curve, gen = HighHammingWeight)
         test_fp12_xy00z0(curve, gen = Long01Sequence)
 
-    test "Sparse 𝔽p12/𝔽p6 resulting from xyz000 line function":
-      proc test_fp12_xyz000(C: static Curve, gen: static RandomGen) =
+    test "Sparse 𝔽p12/𝔽p6 resulting from acb000 line function":
+      proc test_fp12_acb000(C: static Curve, gen: static RandomGen) =
         for _ in 0 ..< Iters:
           var a = rng.random_elem(Fp12[C], gen)
           var a2 = a
@@ -145,9 +144,9 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
           var y = rng.random_elem(Fp2[C], gen)
           var z = rng.random_elem(Fp2[C], gen)
 
-          let line = Line[Fp2[C]](x: x, y: y, z: z)
+          let line = Line[Fp2[C]](a: x, b: y, c: z)
           let b = Fp12[C](
-            c0: Fp6[C](coords: [x, y, z])
+            c0: Fp6[C](coords: [x, z, y])
           )
 
           a *= b
@@ -156,14 +155,14 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
           check: bool(a == a2)
 
       staticFor(curve, TestCurves):
-        test_fp12_xyz000(curve, gen = Uniform)
-        test_fp12_xyz000(curve, gen = HighHammingWeight)
-        test_fp12_xyz000(curve, gen = Long01Sequence)
+        test_fp12_acb000(curve, gen = Uniform)
+        test_fp12_acb000(curve, gen = HighHammingWeight)
+        test_fp12_acb000(curve, gen = Long01Sequence)
   else:
     static: doAssert Fp12[BN254_Snarks]().c0.typeof is Fp4
 
-    test "Sparse 𝔽p12/𝔽p4 resulting from xy000z line function (M-twist only)":
-      proc test_fp12_xy000z(C: static Curve, gen: static RandomGen) =
+    test "Sparse 𝔽p12/𝔽p4 resulting from ca00b0 line function (M-twist only)":
+      proc test_fp12_ca00b0(C: static Curve, gen: static RandomGen) =
         when C.getSexticTwist() == M_Twist:
           for _ in 0 ..< Iters:
             var a = rng.random_elem(Fp12[C], gen)
@@ -173,27 +172,27 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             var y = rng.random_elem(Fp2[C], gen)
             var z = rng.random_elem(Fp2[C], gen)
 
-            let line = Line[Fp2[C]](x: x, y: y, z: z)
+            let line = Line[Fp2[C]](a: x, b: y, c: z)
             let b = Fp12[C](
               coords: [
-                Fp4[C](coords: [x, y]),
+                Fp4[C](coords: [z, x]),
                 Fp4[C](),
-                Fp4[C](coords: [Fp2[C](), z])
+                Fp4[C](coords: [y, Fp2[C]()])
               ]
             )
 
             a *= b
-            a2.mul_sparse_by_line_xy000z(line)
+            a2.mul_sparse_by_line_ca00b0(line)
 
             check: bool(a == a2)
 
       staticFor(curve, TestCurves):
-        test_fp12_xy000z(curve, gen = Uniform)
-        test_fp12_xy000z(curve, gen = HighHammingWeight)
-        test_fp12_xy000z(curve, gen = Long01Sequence)
+        test_fp12_ca00b0(curve, gen = Uniform)
+        test_fp12_ca00b0(curve, gen = HighHammingWeight)
+        test_fp12_ca00b0(curve, gen = Long01Sequence)
 
     test "Sparse 𝔽p12/𝔽p4 resulting from xyz000 line function (D-twist only)":
-      proc test_fp12_xyz000(C: static Curve, gen: static RandomGen) =
+      proc test_fp12_acb000(C: static Curve, gen: static RandomGen) =
         when C.getSexticTwist() == D_Twist:
           for _ in 0 ..< Iters:
             var a = rng.random_elem(Fp12[C], gen)
@@ -203,39 +202,39 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             var y = rng.random_elem(Fp2[C], gen)
             var z = rng.random_elem(Fp2[C], gen)
 
-            let line = Line[Fp2[C]](x: x, y: y, z: z)
+            let line = Line[Fp2[C]](a: x, b: y, c: z)
             let b = Fp12[C](
               coords: [
-                Fp4[C](coords: [x, y]),
-                Fp4[C](coords: [z, Fp2[C]()]),
+                Fp4[C](coords: [x, z]),
+                Fp4[C](coords: [y, Fp2[C]()]),
                 Fp4[C]()
               ]
             )
 
             a *= b
-            a2.mul_sparse_by_line_xyz000(line)
+            a2.mul_sparse_by_line_acb000(line)
 
             check: bool(a == a2)
 
       staticFor(curve, TestCurves):
-        test_fp12_xyz000(curve, gen = Uniform)
-        test_fp12_xyz000(curve, gen = HighHammingWeight)
-        test_fp12_xyz000(curve, gen = Long01Sequence)
+        test_fp12_acb000(curve, gen = Uniform)
+        test_fp12_acb000(curve, gen = HighHammingWeight)
+        test_fp12_acb000(curve, gen = Long01Sequence)
 
-    test "Somewhat-sparse 𝔽p12/𝔽p4 resulting from xy000z*xy000z line functions (M-twist only)":
-      proc test_fp12_xy000z_xy000z(C: static Curve, gen: static RandomGen) =
+    test "Somewhat-sparse 𝔽p12/𝔽p4 resulting from ca00b0*ca00b0 line functions (M-twist only)":
+      proc test_fp12_ca00b0_ca00b0(C: static Curve, gen: static RandomGen) =
         when C.getSexticTwist() == M_Twist:
           for _ in 0 ..< Iters:
             var x0 = rng.random_elem(Fp2[C], gen)
             var y0 = rng.random_elem(Fp2[C], gen)
             var z0 = rng.random_elem(Fp2[C], gen)
 
-            let line0 = Line[Fp2[C]](x: x0, y: y0, z: z0)
+            let line0 = Line[Fp2[C]](a: x0, b: y0, c: z0)
             let f0 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x0, y0]),
+                Fp4[C](coords: [z0, x0]),
                 Fp4[C](),
-                Fp4[C](coords: [Fp2[C](), z0])
+                Fp4[C](coords: [y0, Fp2[C]()])
               ]
             )
 
@@ -243,12 +242,12 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             var y1 = rng.random_elem(Fp2[C], gen)
             var z1 = rng.random_elem(Fp2[C], gen)
 
-            let line1 = Line[Fp2[C]](x: x1, y: y1, z: z1)
+            let line1 = Line[Fp2[C]](a: x1, b: y1, c: z1)
             let f1 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x1, y1]),
+                Fp4[C](coords: [z1, x1]),
                 Fp4[C](),
-                Fp4[C](coords: [Fp2[C](), z1])
+                Fp4[C](coords: [y1, Fp2[C]()])
               ]
             )
 
@@ -256,23 +255,28 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             r.prod(f0, f1)
 
             var rl: Fp12[C]
-            rl.prod_xy000z_xy000z_into_abcd00efghij(line0, line1)
+            rl.prod_zx00y0_zx00y0_into_abcd00efghij(line0, line1)
 
             check: bool(r == rl)
 
-    test "Somewhat-sparse 𝔽p12/𝔽p4 resulting from xyz000*xyz000 line functions (D-twist only)":
-      proc test_fp12_xyz000_xyz000(C: static Curve, gen: static RandomGen) =
+      staticFor(curve, TestCurves):
+        test_fp12_ca00b0_ca00b0(curve, gen = Uniform)
+        test_fp12_ca00b0_ca00b0(curve, gen = HighHammingWeight)
+        test_fp12_ca00b0_ca00b0(curve, gen = Long01Sequence)
+
+    test "Somewhat-sparse 𝔽p12/𝔽p4 resulting from acb000*acb000 line functions (D-twist only)":
+      proc test_fp12_acb000_acb000(C: static Curve, gen: static RandomGen) =
         when C.getSexticTwist() == D_Twist:
           for _ in 0 ..< Iters:
             var x0 = rng.random_elem(Fp2[C], gen)
             var y0 = rng.random_elem(Fp2[C], gen)
             var z0 = rng.random_elem(Fp2[C], gen)
 
-            let line0 = Line[Fp2[C]](x: x0, y: y0, z: z0)
+            let line0 = Line[Fp2[C]](a: x0, b: y0, c: z0)
             let f0 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x0, y0]),
-                Fp4[C](coords: [z0, Fp2[C]()]),
+                Fp4[C](coords: [x0, z0]),
+                Fp4[C](coords: [y0, Fp2[C]()]),
                 Fp4[C]()
               ]
             )
@@ -281,11 +285,11 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             var y1 = rng.random_elem(Fp2[C], gen)
             var z1 = rng.random_elem(Fp2[C], gen)
 
-            let line1 = Line[Fp2[C]](x: x1, y: y1, z: z1)
+            let line1 = Line[Fp2[C]](a: x1, b: y1, c: z1)
             let f1 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x1, y1]),
-                Fp4[C](coords: [z1, Fp2[C]()]),
+                Fp4[C](coords: [x1, z1]),
+                Fp4[C](coords: [y1, Fp2[C]()]),
                 Fp4[C]()
               ]
             )
@@ -294,16 +298,16 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             r.prod(f0, f1)
 
             var rl: Fp12[C]
-            rl.prod_xyz000_xyz000_into_abcdefghij00(line0, line1)
+            rl.prod_xzy000_xzy000_into_abcdefghij00(line0, line1)
 
             check: bool(r == rl)
 
       staticFor(curve, TestCurves):
-        test_fp12_xyz000_xyz000(curve, gen = Uniform)
-        test_fp12_xyz000_xyz000(curve, gen = HighHammingWeight)
-        test_fp12_xyz000_xyz000(curve, gen = Long01Sequence)
+        test_fp12_acb000_acb000(curve, gen = Uniform)
+        test_fp12_acb000_acb000(curve, gen = HighHammingWeight)
+        test_fp12_acb000_acb000(curve, gen = Long01Sequence)
 
-    test "Somewhat-sparse 𝔽p12/𝔽p4 mul by the product (xyz000*xyz000) of line functions (D-twist only)":
+    test "Somewhat-sparse 𝔽p12/𝔽p4 mul by the product (acb000*acb000) of line functions (D-twist only)":
       proc test_fp12_abcdefghij00(C: static Curve, gen: static RandomGen) =
         when C.getSexticTwist() == D_Twist:
           for _ in 0 ..< Iters:
@@ -311,11 +315,11 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             var y0 = rng.random_elem(Fp2[C], gen)
             var z0 = rng.random_elem(Fp2[C], gen)
 
-            let line0 = Line[Fp2[C]](x: x0, y: y0, z: z0)
+            let line0 = Line[Fp2[C]](a: x0, b: y0, c: z0)
             let f0 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x0, y0]),
-                Fp4[C](coords: [z0, Fp2[C]()]),
+                Fp4[C](coords: [x0, z0]),
+                Fp4[C](coords: [y0, Fp2[C]()]),
                 Fp4[C]()
               ]
             )
@@ -324,17 +328,20 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             var y1 = rng.random_elem(Fp2[C], gen)
             var z1 = rng.random_elem(Fp2[C], gen)
 
-            let line1 = Line[Fp2[C]](x: x1, y: y1, z: z1)
+            let line1 = Line[Fp2[C]](a: x1, b: y1, c: z1)
             let f1 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x1, y1]),
-                Fp4[C](coords: [z1, Fp2[C]()]),
+                Fp4[C](coords: [x1, z1]),
+                Fp4[C](coords: [y1, Fp2[C]()]),
                 Fp4[C]()
               ]
             )
 
+            var r: Fp12[C]
+            r.prod(f0, f1)
+
             var rl: Fp12[C]
-            rl.prod_xyz000_xyz000_into_abcdefghij00(line0, line1)
+            rl.prod_xzy000_xzy000_into_abcdefghij00(line0, line1)
 
             var f = rng.random_elem(Fp12[C], gen)
             var f2 = f
@@ -349,20 +356,20 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
         test_fp12_abcdefghij00(curve, gen = HighHammingWeight)
         test_fp12_abcdefghij00(curve, gen = Long01Sequence)
 
-    test "Somewhat-sparse 𝔽p12/𝔽p4 mul by the product (xy000z*xy000z) of line functions (M-twist only)":
-      proc test_fp12_abcd00efghij(C: static Curve, gen: static RandomGen) =
+    test "Somewhat-sparse 𝔽p12/𝔽p4 mul by the product (ca00b0*ca00b0) of line functions (M-twist only)":
+      proc test_fp12_abcdef00ghij(C: static Curve, gen: static RandomGen) =
         when C.getSexticTwist() == M_Twist:
           for _ in 0 ..< Iters:
             var x0 = rng.random_elem(Fp2[C], gen)
             var y0 = rng.random_elem(Fp2[C], gen)
             var z0 = rng.random_elem(Fp2[C], gen)
 
-            let line0 = Line[Fp2[C]](x: x0, y: y0, z: z0)
+            let line0 = Line[Fp2[C]](a: x0, b: y0, c: z0)
             let f0 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x0, y0]),
+                Fp4[C](coords: [z0, x0]),
                 Fp4[C](),
-                Fp4[C](coords: [Fp2[C](), z0])
+                Fp4[C](coords: [y0, Fp2[C]()])
               ]
             )
 
@@ -370,17 +377,20 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             var y1 = rng.random_elem(Fp2[C], gen)
             var z1 = rng.random_elem(Fp2[C], gen)
 
-            let line1 = Line[Fp2[C]](x: x1, y: y1, z: z1)
+            let line1 = Line[Fp2[C]](a: x1, b: y1, c: z1)
             let f1 = Fp12[C](
               coords: [
-                Fp4[C](coords: [x1, y1]),
+                Fp4[C](coords: [z1, x1]),
                 Fp4[C](),
-                Fp4[C](coords: [Fp2[C](), z1])
+                Fp4[C](coords: [y1, Fp2[C]()])
               ]
             )
 
+            var r: Fp12[C]
+            r.prod(f0, f1)
+
             var rl: Fp12[C]
-            rl.prod_xy000z_xy000z_into_abcd00efghij(line0, line1)
+            rl.prod_zx00y0_zx00y0_into_abcd00efghij(line0, line1)
 
             var f = rng.random_elem(Fp12[C], gen)
             var f2 = f
@@ -391,6 +401,6 @@ suite "Pairing - Sparse 𝔽p12 multiplication by line function is consistent wi
             check: bool(f == f2)
 
       staticFor(curve, TestCurves):
-        test_fp12_abcd00efghij(curve, gen = Uniform)
-        test_fp12_abcd00efghij(curve, gen = HighHammingWeight)
-        test_fp12_abcd00efghij(curve, gen = Long01Sequence)
+        test_fp12_abcdef00ghij(curve, gen = Uniform)
+        test_fp12_abcdef00ghij(curve, gen = HighHammingWeight)
+        test_fp12_abcdef00ghij(curve, gen = Long01Sequence)