Add multipairing for BN curves (#194)

benchmarks/bench_pairing_bn254_nogami.nim

@@ -58,6 +58,12 @@ proc main() =
     separator()
     pairingBNBench(curve, Iters)
     separator()
+    staticFor j, 2, 4:
+      pairing_multisingle_BNBench(curve, j, Iters div j)
+      pairing_multipairing_BNBench(curve, j, Iters div j)
+    separator()
+    staticFor j, 4, 9:
+      pairing_multipairing_BNBench(curve, j, Iters div j)
 
 main()
 notes()

benchmarks/bench_pairing_bn254_snarks.nim

@@ -58,6 +58,12 @@ proc main() =
     separator()
     pairingBNBench(curve, Iters)
     separator()
+    staticFor j, 2, 4:
+      pairing_multisingle_BNBench(curve, j, Iters div j)
+      pairing_multipairing_BNBench(curve, j, Iters div j)
+    separator()
+    staticFor j, 4, 9:
+      pairing_multipairing_BNBench(curve, j, Iters div j)
 
 main()
 notes()

benchmarks/bench_pairing_template.nim

@@ -233,10 +233,6 @@ proc pairingBLS12Bench*(C: static Curve, iters: int) =
     f.pairing_bls12(P, Q)
 
 proc pairing_multisingle_BLS12Bench*(C: static Curve, N: static int, iters: int) =
-  let
-    P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1])
-    Q = rng.random_point(ECP_ShortW_Aff[Fp2[C], G2])
-
   var
     Ps {.noInit.}: array[N, ECP_ShortW_Aff[Fp[C], G1]]
     Qs {.noInit.}: array[N, ECP_ShortW_Aff[Fp2[C], G2]]
@@ -277,3 +273,36 @@ proc pairingBNBench*(C: static Curve, iters: int) =
   var f: Fp12[C]
   bench("Pairing BN", C, iters):
     f.pairing_bn(P, Q)
+
+proc pairing_multisingle_BNBench*(C: static Curve, N: static int, iters: int) =
+  var
+    Ps {.noInit.}: array[N, ECP_ShortW_Aff[Fp[C], G1]]
+    Qs {.noInit.}: array[N, ECP_ShortW_Aff[Fp2[C], G2]]
+
+    GTs {.noInit.}: array[N, Fp12[C]]
+
+  for i in 0 ..< N:
+    Ps[i] = rng.random_unsafe(typeof(Ps[0]))
+    Qs[i] = rng.random_unsafe(typeof(Qs[0]))
+
+  var f: Fp12[C]
+  bench("Pairing BN non-batched: " & $N, C, iters):
+    for i in 0 ..< N:
+      GTs[i].pairing_bn(Ps[i], Qs[i])
+
+    f = GTs[0]
+    for i in 1 ..< N:
+      f *= GTs[i]
+
+proc pairing_multipairing_BNBench*(C: static Curve, N: static int, iters: int) =
+  var
+    Ps {.noInit.}: array[N, ECP_ShortW_Aff[Fp[C], G1]]
+    Qs {.noInit.}: array[N, ECP_ShortW_Aff[Fp2[C], G2]]
+
+  for i in 0 ..< N:
+    Ps[i] = rng.random_unsafe(typeof(Ps[0]))
+    Qs[i] = rng.random_unsafe(typeof(Qs[0]))
+
+  var f: Fp12[C]
+  bench("Pairing BN batched:     " & $N, C, iters):
+    f.pairing_bn(Ps, Qs)

constantine.nimble

@@ -175,6 +175,11 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[
   ("tests/math/t_pairing_bn254_snarks_optate.nim", false),
   ("tests/math/t_pairing_bls12_377_optate.nim", false),
   ("tests/math/t_pairing_bls12_381_optate.nim", false),
+
+  # Multi-Pairing
+  # ----------------------------------------------------------
+  ("tests/math/t_pairing_bn254_nogami_multi.nim", false),
+  ("tests/math/t_pairing_bn254_snarks_multi.nim", false),
   ("tests/math/t_pairing_bls12_381_multi.nim", false),
 
   # Prime order fields

constantine/math/curves/bls12_381_pairing.nim

@@ -81,8 +81,6 @@ func millerLoopAddchain*[N: static int](
   f.miller_accum_double_then_add(Ts, Qs, Ps, 32)              # 0b110100100000000100000000000000000000000000000001
   f.miller_accum_double_then_add(Ts, Qs, Ps, 16, add = false) # 0b1101001000000001000000000000000000000000000000010000000000000000
 
-  # TODO: what is the threshold for Karabina's compressed squarings?
-
 func cycl_exp_by_curve_param_div2*(
        r: var Fp12[BLS12_381], a: Fp12[BLS12_381],
        invert = BLS12_381_pairing_ate_param_isNeg) =

constantine/math/curves/bn254_nogami_pairing.nim

@@ -59,6 +59,27 @@ func millerLoopAddchain*(
   # Ate pairing for BN curves needs adjustment after basic Miller loop
   f.millerCorrectionBN(T, Q, P, BN254_Nogami_pairing_ate_param_isNeg)
 
+func millerLoopAddchain*[N: static int](
+       f: var Fp12[BN254_Nogami],
+       Qs: array[N, ECP_ShortW_Aff[Fp2[BN254_Nogami], G2]],
+       Ps: array[N, ECP_ShortW_Aff[Fp[BN254_Nogami], G1]]
+     ) =
+  ## Miller Loop for BN254-Nogami curve
+  ## Computes f{6u+2,Q}(P) with u the BLS curve parameter
+  var Ts {.noInit.}: array[N, ECP_ShortW_Prj[Fp2[BN254_Nogami], G2]]
+
+  f.miller_init_double_then_add(Ts, Qs, Ps, 1)                # 0b11
+  f.miller_accum_double_then_add(Ts, Qs, Ps, 6)               # 0b11000001
+  f.miller_accum_double_then_add(Ts, Qs, Ps, 1)               # 0b110000011
+  f.miller_accum_double_then_add(Ts, Qs, Ps, 54)              # 0b110000011000000000000000000000000000000000000000000000000000001
+  f.miller_accum_double_then_add(Ts, Qs, Ps, 2, add = false)  # 0b11000001100000000000000000000000000000000000000000000000000000100
+
+  # Negative AteParam
+  f.conj()
+
+  for i in 0 ..< N:
+    f.millerCorrectionBN(Ts[i], Qs[i], Ps[i], BN254_Nogami_pairing_ate_param_isNeg)
+
 func cycl_exp_by_curve_param*(
        r: var Fp12[BN254_Nogami], a: Fp12[BN254_Nogami],
        invert = BN254_Nogami_pairing_ate_param_isNeg) =

constantine/math/pairing/miller_loops.nim

@@ -77,16 +77,16 @@ func millerCorrectionBN*[FT, F1, F2](
   when ate_param_isNeg:
     T.neg()
   var V {.noInit.}: typeof(Q)
-  var line {.noInit.}: Line[F2]
+  var line1 {.noInit.}, line2 {.noInit.}: Line[F2]
 
   V.frobenius_psi(Q)
-  line.line_add(T, V, P)
+  line1.line_add(T, V, P)
-  f.mul_by_line(line)
 
   V.frobenius_psi(Q, 2)
   V.neg()
-  line.line_add(T, V, P)
+  line2.line_add(T, V, P)
-  f.mul_by_line(line)
+
+  f.mul_by_2_lines(line1, line2)
 
 # ############################################################
 #                                                            #
@@ -265,6 +265,43 @@ func add_jToN[N: static int, FT, F1, F2](
 
   {.pop.}
 
+template basicMillerLoop*[FT, F1, F2; N: static int](
+       f: var FT,
+       Ts: var array[N, ECP_ShortW_Prj[F2, G2]],
+       line0, line1: var Line[F2],
+       Ps: array[N, ECP_ShortW_Aff[F1, G1]],
+       Qs, nQs: array[N, ECP_ShortW_Aff[F2, G2]],
+       ate_param: untyped,
+       ate_param_isNeg: untyped
+    ) =
+  ## Basic Miller loop iterations
+  # TODO: recompute nQ on-the-fly to save stack space
+  mixin pairing # symbol from zoo_pairings
+
+  static:
+    doAssert FT.C == F1.C
+    doAssert FT.C == F2.C
+
+  f.setOne()
+
+  template u: untyped = pairing(C, ate_param)
+  var u3 = pairing(C, ate_param)
+  u3 *= 3
+  for i in countdown(u3.bits - 2, 1):
+    f.square()
+    double_jToN(f, j=0, line0, line1, Ts, Ps)
+
+    let naf = bit(u3, i).int8 - bit(u, i).int8 # This can throw exception
+    if naf == 1:
+      add_jToN(f, j=0, line0, line1, Ts, Qs, Ps)
+    elif naf == -1:
+      add_jToN(f, j=0, line0, line1, Ts, nQs, Ps)
+
+  when pairing(C, ate_param_isNeg):
+    # In GT, x^-1 == conjugate(x)
+    # Remark 7.1, chapter 7.1.1 of Guide to Pairing-Based Cryptography, El Mrabet, 2017
+    conj(f)
+
 func miller_init_double_then_add*[N: static int, FT, F1, F2](
        f: var FT,
        Ts: var array[N, ECP_ShortW_Prj[F2, G2]],

constantine/math/pairing/pairing_bn.nim

@@ -7,6 +7,7 @@
 # at your option. This file may not be copied, modified, or distributed except according to those terms.
 
 import
+  ../../platforms/abstractions,
   ../config/curves,
   ../extension_fields,
   ../elliptic/[
@@ -51,7 +52,7 @@ func millerLoopGenericBN*[C](
        f: var Fp12[C],
        P: ECP_ShortW_Aff[Fp[C], G1],
        Q: ECP_ShortW_Aff[Fp2[C], G2]
-     ) =
+     ) {.meter.} =
   ## Generic Miller Loop for BN curves
   ## Computes f{6u+2,Q}(P) with u the BN curve parameter
 
@@ -75,6 +76,34 @@ func millerLoopGenericBN*[C](
     pairing(C, ate_param_isNeg)
   )
 
+func millerLoopGenericBN*[N: static int, C](
+       f: var Fp12[C],
+       Ps: array[N, ECP_ShortW_Aff[Fp[C], G1]],
+       Qs: array[N, ECP_ShortW_Aff[Fp2[C], G2]]
+     ) {.meter.} =
+  ## Generic Miller Loop for BN curves
+  ## Computes f{6u+2,Q}(P) with u the BN curve parameter
+
+  var
+    Ts {.noInit.}: array[N, ECP_ShortW_Prj[Fp2[C], G2]]
+    line0 {.noInit.}, line1 {.noInit.}: Line[Fp2[C]]
+    nQs{.noInit.}: typeof(Qs)
+
+  for i in 0 ..< N:
+    Ts[i].fromAffine(Qs[i])
+  for i in 0 ..< N:
+    nQs[i].neg(Qs[i])
+
+  basicMillerLoop(
+    f, Ts, line0, line1,
+    Ps, Qs, nQs,
+    ate_param, ate_param_isNeg
+  )
+
+  # Ate pairing for BN curves needs adjustment after basic Miller loop
+  for i in 0 ..< N:
+    f.millerCorrectionBN(Ts[i], Qs[i], Ps[i], pairing(C, ate_param_isNeg))
+
 func finalExpGeneric[C: static Curve](f: var Fp12[C]) =
   ## A generic and slow implementation of final exponentiation
   ## for sanity checks purposes.
@@ -95,7 +124,7 @@ func pairing_bn_reference*[C](
 # Optimized pairing implementation
 # ----------------------------------------------------------------
 
-func finalExpHard_BN*[C: static Curve](f: var Fp12[C]) =
+func finalExpHard_BN*[C: static Curve](f: var Fp12[C]) {.meter.} =
   ## Hard part of the final exponentiation
   ## Specialized for BN curves
   ##
@@ -150,8 +179,8 @@ func finalExpHard_BN*[C: static Curve](f: var Fp12[C]) =
 func pairing_bn*[C](
        gt: var Fp12[C],
        P: ECP_ShortW_Aff[Fp[C], G1],
-       Q: ECP_ShortW_Aff[Fp2[C], G2]) =
+       Q: ECP_ShortW_Aff[Fp2[C], G2]) {.meter.} =
-  ## Compute the optimal Ate Pairing for BLS12 curves
+  ## Compute the optimal Ate Pairing for BN curves
   ## Input: P ∈ G1, Q ∈ G2
   ## Output: e(P, Q) ∈ Gt
   when C == BN254_Nogami:
@@ -160,3 +189,19 @@ func pairing_bn*[C](
     gt.millerLoopGenericBN(P, Q)
   gt.finalExpEasy()
   gt.finalExpHard_BN()
+
+func pairing_bn*[N: static int, C](
+       gt: var Fp12[C],
+       Ps: array[N, ECP_ShortW_Aff[Fp[C], G1]],
+       Qs: array[N, ECP_ShortW_Aff[Fp2[C], G2]]) {.meter.} =
+  ## Compute the optimal Ate Pairing for BLS12 curves
+  ## Input: an array of Ps ∈ G1 and Qs ∈ G2
+  ## Output:
+  ##   The product of pairings
+  ##   e(P₀, Q₀) * e(P₁, Q₁) * e(P₂, Q₂) * ... * e(Pₙ, Qₙ) ∈ Gt
+  when C == BN254_Nogami:
+    gt.millerLoopAddChain(Qs, Ps)
+  else:
+    gt.millerLoopGenericBN(Ps, Qs)
+  gt.finalExpEasy()
+  gt.finalExpHard_BN()

tests/math/t_pairing_bn254_nogami_multi.nim

@@ -0,0 +1,62 @@
+# 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/[os, times, strformat],
+  # Internals
+  ../../constantine/platforms/abstractions,
+  ../../constantine/math/[arithmetic, extension_fields, ec_shortweierstrass],
+  ../../constantine/math/io/io_extfields,
+  ../../constantine/math/config/curves,
+  ../../constantine/math/pairing/pairing_bn,
+  # Test utilities
+  ../../helpers/prng_unsafe
+
+# Testing multipairing
+# ----------------------------------------------
+
+var rng: RngState
+let timeseed = uint32(toUnix(getTime()) and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+seed(rng, timeseed)
+echo "\n------------------------------------------------------\n"
+echo "test_pairing_bn254_nogami_multi xoshiro512** seed: ", timeseed
+
+proc testMultiPairing(rng: var RngState, N: static int) =
+  var
+    Ps {.noInit.}: array[N, ECP_ShortW_Aff[Fp[BN254_Nogami], G1]]
+    Qs {.noInit.}: array[N, ECP_ShortW_Aff[Fp2[BN254_Nogami], G2]]
+
+    GTs {.noInit.}: array[N, Fp12[BN254_Nogami]]
+
+  for i in 0 ..< N:
+    Ps[i] = rng.random_unsafe(typeof(Ps[0]))
+    Qs[i] = rng.random_unsafe(typeof(Qs[0]))
+
+  # Simple pairing
+  let clockSimpleStart = cpuTime()
+  var GTsimple {.noInit.}: Fp12[BN254_Nogami]
+  for i in 0 ..< N:
+    GTs[i].pairing_bn(Ps[i], Qs[i])
+
+  GTsimple = GTs[0]
+  for i in 1 ..< N:
+    GTsimple *= GTs[i]
+  let clockSimpleStop = cpuTime()
+
+  # Multipairing
+  let clockMultiStart = cpuTime()
+  var GTmulti {.noInit.}: Fp12[BN254_Nogami]
+  GTmulti.pairing_bn(Ps, Qs)
+  let clockMultiStop = cpuTime()
+
+  echo &"N={N}, Simple: {clockSimpleStop - clockSimpleStart:>4.4f}s, Multi: {clockMultiStop - clockMultiStart:>4.4f}s"
+  doAssert bool GTsimple == GTmulti
+
+staticFor i, 1, 17:
+  rng.testMultiPairing(N = i)

tests/math/t_pairing_bn254_snarks_multi.nim

@@ -0,0 +1,62 @@
+# 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/[os, times, strformat],
+  # Internals
+  ../../constantine/platforms/abstractions,
+  ../../constantine/math/[arithmetic, extension_fields, ec_shortweierstrass],
+  ../../constantine/math/io/io_extfields,
+  ../../constantine/math/config/curves,
+  ../../constantine/math/pairing/pairing_bn,
+  # Test utilities
+  ../../helpers/prng_unsafe
+
+# Testing multipairing
+# ----------------------------------------------
+
+var rng: RngState
+let timeseed = uint32(toUnix(getTime()) and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+seed(rng, timeseed)
+echo "\n------------------------------------------------------\n"
+echo "test_pairing_bn254_snarks_multi xoshiro512** seed: ", timeseed
+
+proc testMultiPairing(rng: var RngState, N: static int) =
+  var
+    Ps {.noInit.}: array[N, ECP_ShortW_Aff[Fp[BN254_Snarks], G1]]
+    Qs {.noInit.}: array[N, ECP_ShortW_Aff[Fp2[BN254_Snarks], G2]]
+
+    GTs {.noInit.}: array[N, Fp12[BN254_Snarks]]
+
+  for i in 0 ..< N:
+    Ps[i] = rng.random_unsafe(typeof(Ps[0]))
+    Qs[i] = rng.random_unsafe(typeof(Qs[0]))
+
+  # Simple pairing
+  let clockSimpleStart = cpuTime()
+  var GTsimple {.noInit.}: Fp12[BN254_Snarks]
+  for i in 0 ..< N:
+    GTs[i].pairing_bn(Ps[i], Qs[i])
+
+  GTsimple = GTs[0]
+  for i in 1 ..< N:
+    GTsimple *= GTs[i]
+  let clockSimpleStop = cpuTime()
+
+  # Multipairing
+  let clockMultiStart = cpuTime()
+  var GTmulti {.noInit.}: Fp12[BN254_Snarks]
+  GTmulti.pairing_bn(Ps, Qs)
+  let clockMultiStop = cpuTime()
+
+  echo &"N={N}, Simple: {clockSimpleStop - clockSimpleStart:>4.4f}s, Multi: {clockMultiStop - clockMultiStart:>4.4f}s"
+  doAssert bool GTsimple == GTmulti
+
+staticFor i, 1, 17:
+  rng.testMultiPairing(N = i)