diff --git a/benchmarks/bench_fields_template.nim b/benchmarks/bench_fields_template.nim index ef16909..67aa0bb 100644 --- a/benchmarks/bench_fields_template.nim +++ b/benchmarks/bench_fields_template.nim @@ -91,6 +91,21 @@ proc sqrBench*(T: typedesc, iters: int) = bench("Squaring", T, iters): r.square(x) +proc mulUnrBench*(T: typedesc, iters: int) = + var r: doublePrec(T) + let x = rng.random_unsafe(T) + let y = rng.random_unsafe(T) + preventOptimAway(r) + bench("Multiplication unreduced", T, iters): + r.prod2x(x, y) + +proc sqrUnrBench*(T: typedesc, iters: int) = + var r: doublePrec(T) + let x = rng.random_unsafe(T) + preventOptimAway(r) + bench("Squaring unreduced", T, iters): + r.square2x(x) + proc toBigBench*(T: typedesc, iters: int) = var r: matchingBigInt(T.C) let x = rng.random_unsafe(T) diff --git a/benchmarks/bench_fp.nim b/benchmarks/bench_fp.nim index 9ca81ce..4e97804 100644 --- a/benchmarks/bench_fp.nim +++ b/benchmarks/bench_fp.nim @@ -50,6 +50,9 @@ proc main() = mulBench(Fp[curve], Iters) sqrBench(Fp[curve], Iters) smallSeparator() + mulUnrBench(Fp[curve], Iters) + sqrUnrBench(Fp[curve], Iters) + smallSeparator() toBigBench(Fp[curve], Iters) toFieldBench(Fp[curve], Iters) smallSeparator() diff --git a/benchmarks/bench_fp2.nim b/benchmarks/bench_fp2.nim index 267a8f0..78e9da7 100644 --- a/benchmarks/bench_fp2.nim +++ b/benchmarks/bench_fp2.nim @@ -42,8 +42,13 @@ proc main() = negBench(Fp2[curve], Iters) ccopyBench(Fp2[curve], Iters) div2Bench(Fp2[curve], Iters) + smallSeparator() mulBench(Fp2[curve], Iters) sqrBench(Fp2[curve], Iters) + smallSeparator() + mulUnrBench(Fp2[curve], Iters) + sqrUnrBench(Fp2[curve], Iters) + smallSeparator() invBench(Fp2[curve], InvIters) sqrtBench(Fp2[curve], InvIters) separator() diff --git a/benchmarks/bench_fp4.nim b/benchmarks/bench_fp4.nim index a056cf0..c3fa0c2 100644 --- a/benchmarks/bench_fp4.nim +++ b/benchmarks/bench_fp4.nim @@ -42,8 +42,13 @@ proc main() = negBench(Fp4[curve], Iters) ccopyBench(Fp4[curve], Iters) div2Bench(Fp4[curve], Iters) + smallSeparator() mulBench(Fp4[curve], Iters) sqrBench(Fp4[curve], Iters) + smallSeparator() + mulUnrBench(Fp2[curve], Iters) + sqrUnrBench(Fp2[curve], Iters) + smallSeparator() invBench(Fp4[curve], InvIters) separator() diff --git a/benchmarks/bench_fp6.nim b/benchmarks/bench_fp6.nim index 6693860..00151e0 100644 --- a/benchmarks/bench_fp6.nim +++ b/benchmarks/bench_fp6.nim @@ -40,8 +40,13 @@ proc main() = addBench(Fp6[curve], Iters) subBench(Fp6[curve], Iters) negBench(Fp6[curve], Iters) + smallSeparator() mulBench(Fp6[curve], Iters) sqrBench(Fp6[curve], Iters) + smallSeparator() + mulUnrBench(Fp6[curve], Iters) + sqrUnrBench(Fp6[curve], Iters) + smallSeparator() invBench(Fp6[curve], InvIters) separator() diff --git a/benchmarks/bench_pairing_bls12_377.nim b/benchmarks/bench_pairing_bls12_377.nim index e36272c..d5897d1 100644 --- a/benchmarks/bench_pairing_bls12_377.nim +++ b/benchmarks/bench_pairing_bls12_377.nim @@ -42,6 +42,14 @@ proc main() = mulFp12_by_2lines_v1_xyz000_Bench(curve, Iters) mulFp12_by_2lines_v2_xyz000_Bench(curve, Iters) separator() + mulBench(curve, Iters) + sqrBench(curve, Iters) + separator() + cyclotomicSquare_Bench(curve, Iters) + cyclotomicSquareCompressed_Bench(curve, Iters) + cyclotomicDecompression_Bench(curve, Iters) + expCurveParamBench(curve, Iters) + separator() finalExpEasyBench(curve, Iters) finalExpHardBLS12Bench(curve, Iters) separator() @@ -49,7 +57,14 @@ proc main() = finalExpBLS12Bench(curve, Iters) separator() pairingBLS12Bench(curve, Iters) + pairing_multipairing_BLS12Bench(curve, 1, Iters) separator() + staticFor j, 2, 4: + pairing_multisingle_BLS12Bench(curve, j, Iters div j) + pairing_multipairing_BLS12Bench(curve, j, Iters div j) + separator() + staticFor j, 4, 9: + pairing_multipairing_BLS12Bench(curve, j, Iters div j) main() notes() diff --git a/benchmarks/bench_pairing_bls12_381.nim b/benchmarks/bench_pairing_bls12_381.nim index cdf19eb..84d27a5 100644 --- a/benchmarks/bench_pairing_bls12_381.nim +++ b/benchmarks/bench_pairing_bls12_381.nim @@ -42,6 +42,14 @@ proc main() = mulFp12_by_2lines_v1_xy000z_Bench(curve, Iters) mulFp12_by_2lines_v2_xy000z_Bench(curve, Iters) separator() + mulBench(curve, Iters) + sqrBench(curve, Iters) + separator() + cyclotomicSquare_Bench(curve, Iters) + cyclotomicSquareCompressed_Bench(curve, Iters) + cyclotomicDecompression_Bench(curve, Iters) + expCurveParamBench(curve, Iters) + separator() finalExpEasyBench(curve, Iters) finalExpHardBLS12Bench(curve, Iters) separator() @@ -51,9 +59,11 @@ proc main() = pairingBLS12Bench(curve, Iters) pairing_multipairing_BLS12Bench(curve, 1, Iters) separator() - staticFor j, 2, 17: + staticFor j, 2, 4: pairing_multisingle_BLS12Bench(curve, j, Iters div j) pairing_multipairing_BLS12Bench(curve, j, Iters div j) - + separator() + staticFor j, 4, 9: + pairing_multipairing_BLS12Bench(curve, j, Iters div j) main() notes() diff --git a/benchmarks/bench_pairing_bn254_nogami.nim b/benchmarks/bench_pairing_bn254_nogami.nim index 7d134fa..7177581 100644 --- a/benchmarks/bench_pairing_bn254_nogami.nim +++ b/benchmarks/bench_pairing_bn254_nogami.nim @@ -42,6 +42,14 @@ proc main() = mulFp12_by_2lines_v1_xyz000_Bench(curve, Iters) mulFp12_by_2lines_v2_xyz000_Bench(curve, Iters) separator() + mulBench(curve, Iters) + sqrBench(curve, Iters) + separator() + cyclotomicSquare_Bench(curve, Iters) + cyclotomicSquareCompressed_Bench(curve, Iters) + cyclotomicDecompression_Bench(curve, Iters) + expCurveParamBench(curve, Iters) + separator() finalExpEasyBench(curve, Iters) finalExpHardBNBench(curve, Iters) separator() diff --git a/benchmarks/bench_pairing_bn254_snarks.nim b/benchmarks/bench_pairing_bn254_snarks.nim index c8be167..a5cae37 100644 --- a/benchmarks/bench_pairing_bn254_snarks.nim +++ b/benchmarks/bench_pairing_bn254_snarks.nim @@ -42,6 +42,14 @@ proc main() = mulFp12_by_2lines_v1_xyz000_Bench(curve, Iters) mulFp12_by_2lines_v2_xyz000_Bench(curve, Iters) separator() + mulBench(curve, Iters) + sqrBench(curve, Iters) + separator() + cyclotomicSquare_Bench(curve, Iters) + cyclotomicSquareCompressed_Bench(curve, Iters) + cyclotomicDecompression_Bench(curve, Iters) + expCurveParamBench(curve, Iters) + separator() finalExpEasyBench(curve, Iters) finalExpHardBNBench(curve, Iters) separator() diff --git a/benchmarks/bench_pairing_template.nim b/benchmarks/bench_pairing_template.nim index 71352aa..ff16507 100644 --- a/benchmarks/bench_pairing_template.nim +++ b/benchmarks/bench_pairing_template.nim @@ -20,9 +20,8 @@ import ../constantine/ec_shortweierstrass, ../constantine/curves/zoo_subgroups, ../constantine/pairing/[ - cyclotomic_fp12, - lines_projective, - mul_fp12_by_lines, + cyclotomic_subgroup, + lines_eval, pairing_bls12, pairing_bn ], @@ -181,6 +180,54 @@ proc mulFp12_by_2lines_v2_xy000z_Bench*(C: static Curve, iters: int) = f2.prod_xy000z_xy000z_into_abcd00efghij(l0, l1) f.mul_sparse_by_abcd00efghij(f2) +proc mulBench*(C: static Curve, iters: int) = + var r: Fp12[C] + let x = rng.random_unsafe(Fp12[C]) + let y = rng.random_unsafe(Fp12[C]) + preventOptimAway(r) + bench("Multiplication 𝔽p12", C, iters): + r.prod(x, y) + +proc sqrBench*(C: static Curve, iters: int) = + var r: Fp12[C] + let x = rng.random_unsafe(Fp12[C]) + preventOptimAway(r) + bench("Squaring 𝔽p12", C, iters): + r.square(x) + +proc cyclotomicSquare_Bench*(C: static Curve, iters: int) = + var f = rng.random_unsafe(Fp12[C]) + + bench("Squaring 𝔽p12 in cyclotomic subgroup", C, iters): + f.cyclotomic_square() + +proc expCurveParamBench*(C: static Curve, iters: int) = + var f = rng.random_unsafe(Fp12[C]) + + bench("Cyclotomic Exp by curve parameter", C, iters): + f.cycl_exp_by_curve_param(f) + +proc cyclotomicSquareCompressed_Bench*(C: static Curve, iters: int) = + var f = rng.random_unsafe(Fp12[C]) + var g: G2345[Fp2[C]] + g.fromFpk(f) + + bench("Cyclotomic Compressed Squaring 𝔽p12", C, iters): + g.cyclotomic_square_compressed() + +proc cyclotomicDecompression_Bench*(C: static Curve, iters: int) = + var f = rng.random_unsafe(Fp12[C]) + var gs: array[1, G2345[Fp2[C]]] + gs[0].fromFpk(f) + + var g1s_ratio: array[1, tuple[g1_num, g1_den: Fp2[C]]] + var g0s, g1s: array[1, Fp2[C]] + + bench("Cyclotomic Decompression 𝔽p12", C, iters): + recover_g1(g1s_ratio[0].g1_num, g1s_ratio[0].g1_den, gs[0]) + g1s.batch_ratio_g1s(g1s_ratio) + g0s[0].recover_g0(g1s[0], gs[0]) + proc millerLoopBLS12Bench*(C: static Curve, iters: int) = let P = rng.random_point(ECP_ShortW_Aff[Fp[C], G1]) @@ -253,7 +300,7 @@ proc pairing_multisingle_BLS12Bench*(C: static Curve, N: static int, iters: int) Qs[i] = rng.random_unsafe(typeof(Qs[0])) var f: Fp12[C] - bench("Pairing BLS12 multi-single " & $N & " pairings", C, iters): + bench("Pairing BLS12 non-batched: " & $N, C, iters): for i in 0 ..< N: GTs[i].pairing_bls12(Ps[i], Qs[i]) @@ -271,7 +318,7 @@ proc pairing_multipairing_BLS12Bench*(C: static Curve, N: static int, iters: int Qs[i] = rng.random_unsafe(typeof(Qs[0])) var f: Fp12[C] - bench("Pairing BLS12 multipairing " & $N & " pairings", C, iters): + bench("Pairing BLS12 batched: " & $N, C, iters): f.pairing_bls12(Ps, Qs) proc pairingBNBench*(C: static Curve, iters: int) = diff --git a/benchmarks/bench_summary_bls12_377.nim b/benchmarks/bench_summary_bls12_377.nim index a06fb49..b7d465c 100644 --- a/benchmarks/bench_summary_bls12_377.nim +++ b/benchmarks/bench_summary_bls12_377.nim @@ -49,29 +49,30 @@ proc main() = invBench(Fp2[curve], Iters) sqrtBench(Fp2[curve], Iters) separator() + mulBench(Fp12[curve], Iters) + sqrBench(Fp12[curve], Iters) + invBench(Fp12[curve], Iters) + separator() addBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) separator() addBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() - mulBench(Fp12[curve], Iters) - sqrBench(Fp12[curve], Iters) - invBench(Fp12[curve], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() millerLoopBLS12Bench(curve, Iters) finalExpBLS12Bench(curve, Iters) diff --git a/benchmarks/bench_summary_bls12_381.nim b/benchmarks/bench_summary_bls12_381.nim index 09de80e..56d4eb2 100644 --- a/benchmarks/bench_summary_bls12_381.nim +++ b/benchmarks/bench_summary_bls12_381.nim @@ -36,6 +36,7 @@ proc main() = staticFor i, 0, AvailableCurves.len: const curve = AvailableCurves[i] + mulBench(Fr[curve], Iters) sqrBench(Fr[curve], Iters) separator() @@ -49,29 +50,30 @@ proc main() = invBench(Fp2[curve], Iters) sqrtBench(Fp2[curve], Iters) separator() + mulBench(Fp12[curve], Iters) + sqrBench(Fp12[curve], Iters) + invBench(Fp12[curve], Iters) + separator() addBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) separator() addBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() - mulBench(Fp12[curve], Iters) - sqrBench(Fp12[curve], Iters) - invBench(Fp12[curve], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() millerLoopBLS12Bench(curve, Iters) finalExpBLS12Bench(curve, Iters) diff --git a/benchmarks/bench_summary_bn254_nogami.nim b/benchmarks/bench_summary_bn254_nogami.nim index 1f85d02..f968dfc 100644 --- a/benchmarks/bench_summary_bn254_nogami.nim +++ b/benchmarks/bench_summary_bn254_nogami.nim @@ -49,29 +49,30 @@ proc main() = invBench(Fp2[curve], Iters) sqrtBench(Fp2[curve], Iters) separator() + mulBench(Fp12[curve], Iters) + sqrBench(Fp12[curve], Iters) + invBench(Fp12[curve], Iters) + separator() addBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) separator() addBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() - mulBench(Fp12[curve], Iters) - sqrBench(Fp12[curve], Iters) - invBench(Fp12[curve], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() millerLoopBNBench(curve, Iters) finalExpBNBench(curve, Iters) diff --git a/benchmarks/bench_summary_bn254_snarks.nim b/benchmarks/bench_summary_bn254_snarks.nim index ae7741c..be20e68 100644 --- a/benchmarks/bench_summary_bn254_snarks.nim +++ b/benchmarks/bench_summary_bn254_snarks.nim @@ -49,29 +49,30 @@ proc main() = invBench(Fp2[curve], Iters) sqrtBench(Fp2[curve], Iters) separator() + mulBench(Fp12[curve], Iters) + sqrBench(Fp12[curve], Iters) + invBench(Fp12[curve], Iters) + separator() addBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) mixedAddBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) doublingBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) separator() addBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) separator() addBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) mixedAddBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) doublingBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) - scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() - mulBench(Fp12[curve], Iters) - sqrBench(Fp12[curve], Iters) - invBench(Fp12[curve], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp[curve], G1], Iters) + scalarMulBench(ECP_ShortW_Prj[Fp2[curve], G2], Iters) + scalarMulBench(ECP_ShortW_Jac[Fp2[curve], G2], Iters) separator() millerLoopBNBench(curve, Iters) finalExpBNBench(curve, Iters) diff --git a/benchmarks/bench_summary_template.nim b/benchmarks/bench_summary_template.nim index d40ca52..2fdf9b7 100644 --- a/benchmarks/bench_summary_template.nim +++ b/benchmarks/bench_summary_template.nim @@ -24,7 +24,7 @@ import ../constantine/curves/zoo_subgroups, ../constantine/hash_to_curve/hash_to_curve, ../constantine/pairing/[ - cyclotomic_fp12, + cyclotomic_subgroup, pairing_bls12, pairing_bn ], diff --git a/constantine.nimble b/constantine.nimble index 7b149a3..0f58dc1 100644 --- a/constantine.nimble +++ b/constantine.nimble @@ -168,7 +168,7 @@ const testDesc: seq[tuple[path: string, useGMP: bool]] = @[ # ("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_cyclotomic_subgroup.nim", false), ("tests/t_pairing_bn254_nogami_optate.nim", false), ("tests/t_pairing_bn254_snarks_optate.nim", false), ("tests/t_pairing_bls12_377_optate.nim", false), diff --git a/constantine/arithmetic/assembly/limbs_asm_mul_mont_x86_adx_bmi2.nim b/constantine/arithmetic/assembly/limbs_asm_mul_mont_x86_adx_bmi2.nim index 2a1811b..b845500 100644 --- a/constantine/arithmetic/assembly/limbs_asm_mul_mont_x86_adx_bmi2.nim +++ b/constantine/arithmetic/assembly/limbs_asm_mul_mont_x86_adx_bmi2.nim @@ -279,7 +279,7 @@ macro mulMont_CIOS_sparebit_adx_gen[N: static int]( result.add ctx.generate -func mulMont_CIOS_sparebit_asm_adx*(r: var Limbs, a, b, M: Limbs, m0ninv: BaseType, skipFinalSub: static bool = false) = +func mulMont_CIOS_sparebit_asm_adx_inline*(r: var Limbs, a, b, M: Limbs, m0ninv: BaseType, skipFinalSub: static bool = false) {.inline.} = ## Constant-time Montgomery multiplication ## If "skipFinalSub" is set ## the result is in the range [0, 2M) @@ -288,6 +288,15 @@ func mulMont_CIOS_sparebit_asm_adx*(r: var Limbs, a, b, M: Limbs, m0ninv: BaseTy ## This procedure can only be called if the modulus doesn't use the full bitwidth of its underlying representation r.mulMont_CIOS_sparebit_adx_gen(a, b, M, m0ninv, skipFinalSub) +func mulMont_CIOS_sparebit_asm_adx*(r: var Limbs, a, b, M: Limbs, m0ninv: BaseType, skipFinalSub: static bool = false) = + ## Constant-time Montgomery multiplication + ## If "skipFinalSub" is set + ## the result is in the range [0, 2M) + ## otherwise the result is in the range [0, M) + ## + ## This procedure can only be called if the modulus doesn't use the full bitwidth of its underlying representation + r.mulMont_CIOS_sparebit_asm_adx_inline(a, b, M, m0ninv, skipFinalSub) + # Montgomery Squaring # ------------------------------------------------------------ diff --git a/constantine/curves/bls12_377_pairing.nim b/constantine/curves/bls12_377_pairing.nim index 2072851..0430619 100644 --- a/constantine/curves/bls12_377_pairing.nim +++ b/constantine/curves/bls12_377_pairing.nim @@ -11,7 +11,7 @@ import ../io/io_bigints, ../towers, ../elliptic/[ec_shortweierstrass_affine, ec_shortweierstrass_projective], - ../pairing/[cyclotomic_fp12, miller_loops], + ../pairing/[cyclotomic_subgroup, miller_loops], ../isogeny/frobenius # Slow generic implementation @@ -58,29 +58,37 @@ func millerLoopAddchain*( f.miller_accum_double_then_add(T, Q, P, 1) # 0b100001010000100011 f.miller_accum_double_then_add(T, Q, P, 46, add = true) # 0b1000010100001000110000000000000000000000000000000000000000000001 -func pow_x*(r: var Fp12[BLS12_377], a: Fp12[BLS12_377], invert = BLS12_377_pairing_ate_param_isNeg) = +func millerLoopAddchain*[N: static int]( + f: var Fp12[BLS12_377], + Qs: array[N, ECP_ShortW_Aff[Fp2[BLS12_377], G2]], + Ps: array[N, ECP_ShortW_Aff[Fp[BLS12_377], G1]] + ) = + ## Miller Loop for BLS12-377 curve + ## Computes f{u,Q}(P) with u the BLS curve parameter + + var Ts {.noInit.}: array[N, ECP_ShortW_Prj[Fp2[BLS12_377], G2]] + + f.miller_init_double_then_add(Ts, Qs, Ps, 5) # 0b100001 + f.miller_accum_double_then_add(Ts, Qs, Ps, 2) # 0b10000101 + f.miller_accum_double_then_add(Ts, Qs, Ps, 5) # 0b1000010100001 + f.miller_accum_double_then_add(Ts, Qs, Ps, 4) # 0b10000101000010001 + f.miller_accum_double_then_add(Ts, Qs, Ps, 1) # 0b100001010000100011 + f.miller_accum_double_then_add(Ts, Qs, Ps, 46, add = true) # 0b1000010100001000110000000000000000000000000000000000000000000001 + +func cycl_exp_by_curve_param*(r: var Fp12[BLS12_377], a: Fp12[BLS12_377], invert = BLS12_377_pairing_ate_param_isNeg) = ## f^x with x the curve parameter ## For BLS12_377 f^0x8508c00000000001 - r.cyclotomic_square(a) + r.cycl_sqr_repeated(a, 5) + r *= a + let t{.noInit.} = r + r.cycl_sqr_repeated(7) + r *= t + r.cycl_sqr_repeated(4) r *= a r.cyclotomic_square() r *= a - let t111 = r - r.cycl_sqr_repeated(2) - let t111000 = r - - r *= t111 - let t100011 = r - - r.cyclotomic_square() - r *= t100011 - r *= t111000 - - r.cycl_sqr_repeated(10) - r *= t100011 - - r.cycl_sqr_repeated(46) # TODO: Karabina's compressed squarings + r.cyclotomic_exp_compressed(r, [46]) r *= a if invert: @@ -95,6 +103,6 @@ func isInPairingSubgroup*(a: Fp12[BLS12_377]): SecretBool = # a is in the GT subgroup iff a^p == a^u var t0{.noInit.}, t1{.noInit.}: Fp12[BLS12_377] t0.frobenius_map(a) - t1.pow_x(a) + t1.cycl_exp_by_curve_param(a) return t0 == t1 \ No newline at end of file diff --git a/constantine/curves/bls12_381_pairing.nim b/constantine/curves/bls12_381_pairing.nim index 3effbf3..4f22cb5 100644 --- a/constantine/curves/bls12_381_pairing.nim +++ b/constantine/curves/bls12_381_pairing.nim @@ -11,7 +11,7 @@ import ../io/io_bigints, ../towers, ../elliptic/[ec_shortweierstrass_affine, ec_shortweierstrass_projective], - ../pairing/[cyclotomic_fp12, miller_loops], + ../pairing/[cyclotomic_subgroup, miller_loops], ../isogeny/frobenius # Slow generic implementation @@ -82,30 +82,45 @@ func millerLoopAddchain*[N: static int]( # TODO: what is the threshold for Karabina's compressed squarings? -func pow_xdiv2*(r: var Fp12[BLS12_381], a: Fp12[BLS12_381], invert = BLS12_381_pairing_ate_param_isNeg) = +func cycl_exp_by_curve_param_div2*( + r: var Fp12[BLS12_381], a: Fp12[BLS12_381], + invert = BLS12_381_pairing_ate_param_isNeg) = ## f^(x/2) with x the curve parameter - ## For BLS12_381 f^-0xd201000000010000 + ## For BLS12_381 f^-0xd201000000010000 = 0b1101001000000001000000000000000000000000000000010000000000000000 + + # Squarings accumulator + var s{.noInit.}: Fp12[BLS12_381] - r.cyclotomic_square(a) - r *= a - r.cycl_sqr_repeated(2) - r *= a - r.cycl_sqr_repeated(3) - r *= a - r.cycl_sqr_repeated(9) - r *= a - r.cycl_sqr_repeated(32) # TODO: use Karabina? - r *= a - r.cycl_sqr_repeated(16-1) # Don't do the last iteration + r.cyclotomic_exp_compressed(s, a, [16-1, 32, 9]) + s.cycl_sqr_repeated(3) + r *= s + s.cycl_sqr_repeated(2) + r *= s + s.cyclotomic_square() + r *= s if invert: r.cyclotomic_inv() -func pow_x*(r: var Fp12[BLS12_381], a: Fp12[BLS12_381], invert = BLS12_381_pairing_ate_param_isNeg) = +func cycl_exp_by_curve_param*( + r: var Fp12[BLS12_381], a: Fp12[BLS12_381], + invert = BLS12_381_pairing_ate_param_isNeg) = ## f^x with x the curve parameter ## For BLS12_381 f^-0xd201000000010000 - r.pow_xdiv2(a, invert) - r.cyclotomic_square() + + # Squarings accumulator + var s{.noInit.}: Fp12[BLS12_381] + + r.cyclotomic_exp_compressed(s, a, [16, 32, 9]) + s.cycl_sqr_repeated(3) + r *= s + s.cycl_sqr_repeated(2) + r *= s + s.cyclotomic_square() + r *= s + + if invert: + r.cyclotomic_inv() func isInPairingSubgroup*(a: Fp12[BLS12_381]): SecretBool = ## Returns true if a is in GT subgroup, i.e. a is an element of order r @@ -116,6 +131,6 @@ func isInPairingSubgroup*(a: Fp12[BLS12_381]): SecretBool = # P is in the G1 subgroup iff a^p == a^u var t0{.noInit.}, t1{.noInit.}: Fp12[BLS12_381] t0.frobenius_map(a) - t1.pow_x(a) + t1.cycl_exp_by_curve_param(a) return t0 == t1 \ No newline at end of file diff --git a/constantine/curves/bn254_nogami_pairing.nim b/constantine/curves/bn254_nogami_pairing.nim index 57245ef..ae21841 100644 --- a/constantine/curves/bn254_nogami_pairing.nim +++ b/constantine/curves/bn254_nogami_pairing.nim @@ -11,7 +11,7 @@ import ../io/io_bigints, ../towers, ../elliptic/[ec_shortweierstrass_affine, ec_shortweierstrass_projective], - ../pairing/[cyclotomic_fp12, miller_loops], + ../pairing/[cyclotomic_subgroup, miller_loops], ../isogeny/frobenius # Slow generic implementation @@ -58,13 +58,12 @@ 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 pow_u*(r: var Fp12[BN254_Nogami], a: Fp12[BN254_Nogami], invert = 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) = ## f^u with u the curve parameter - ## For BN254_Nogami f^-0x4080000000000001 - r.cyclotomic_square(a) - r.cycl_sqr_repeated(6) - r *= a - r.cycl_sqr_repeated(55) + ## For BN254_Nogami f^-0x4080000000000001 = 0b100000010000000000000000000000000000000000000000000000000000001 + r.cyclotomic_exp_compressed(a, [55, 7]) r *= a if invert: @@ -78,8 +77,8 @@ func isInPairingSubgroup*(a: Fp12[BN254_Nogami]): SecretBool = # on BLS pairing-friendly curves # P is in the G1 subgroup iff a^p == a^(6u²) var t0{.noInit.}, t1{.noInit.}: Fp12[BN254_Nogami] - t0.pow_u(a) # a^p - t1.pow_u(t0) # a^(p²) + t0.cycl_exp_by_curve_param(a) # a^p + t1.cycl_exp_by_curve_param(t0) # a^(p²) t0.square(t1) # a^(2p²) t0 *= t1 # a^(3p²) t0.square() # a^(6p²) diff --git a/constantine/curves/bn254_snarks_pairing.nim b/constantine/curves/bn254_snarks_pairing.nim index dc681b8..3f99510 100644 --- a/constantine/curves/bn254_snarks_pairing.nim +++ b/constantine/curves/bn254_snarks_pairing.nim @@ -11,7 +11,7 @@ import ../io/io_bigints, ../towers, ../elliptic/[ec_shortweierstrass_affine, ec_shortweierstrass_projective], - ../pairing/[cyclotomic_fp12, miller_loops], + ../pairing/[cyclotomic_subgroup, miller_loops], ../isogeny/frobenius # Slow generic implementation @@ -41,88 +41,65 @@ const BN254_Snarks_pairing_finalexponent* = block: # it would requires saving accumulators to actually save # operations compared to NAF, and can we combine the saved EC[Fp2] accumulators? -func pow_u*(r: var Fp12[BN254_Snarks], a: Fp12[BN254_Snarks], invert = BN254_Snarks_pairing_ate_param_isNeg) = +func cycl_exp_by_curve_param*( + r: var Fp12[BN254_Snarks], a: Fp12[BN254_Snarks], + invert = BN254_Snarks_pairing_ate_param_isNeg) = ## f^u with u the curve parameter ## For BN254_Snarks f^0x44e992b44a6909f1 - when false: - cyclotomic_exp( - r, a, - BigInt[63].fromHex("0x44e992b44a6909f1"), - invert - ) - else: - var # Hopefully the compiler optimizes away unused Fp12 - # because those are huge - x10 {.noInit.}: Fp12[BN254_Snarks] - x11 {.noInit.}: Fp12[BN254_Snarks] - x100 {.noInit.}: Fp12[BN254_Snarks] - x110 {.noInit.}: Fp12[BN254_Snarks] - x1100 {.noInit.}: Fp12[BN254_Snarks] - x1111 {.noInit.}: Fp12[BN254_Snarks] - x10010 {.noInit.}: Fp12[BN254_Snarks] - x10110 {.noInit.}: Fp12[BN254_Snarks] - x11100 {.noInit.}: Fp12[BN254_Snarks] - x101110 {.noInit.}: Fp12[BN254_Snarks] - x1001010 {.noInit.}: Fp12[BN254_Snarks] - x1111000 {.noInit.}: Fp12[BN254_Snarks] - x10001110 {.noInit.}: Fp12[BN254_Snarks] + # https://github.com/mmcloughlin/addchain + # Addchain weighted by Fp12 mul and cyclotomic square cycle costs + # addchain search -add 3622 -double 1696 "0x44e992b44a6909f1" + var # Hopefully the compiler optimizes away unused Fp12 + # because those are huge + x10 {.noInit.}: Fp12[BN254_Snarks] + x100 {.noInit.}: Fp12[BN254_Snarks] + x1000 {.noInit.}: Fp12[BN254_Snarks] + x10000 {.noInit.}: Fp12[BN254_Snarks] + x10001 {.noInit.}: Fp12[BN254_Snarks] + x10011 {.noInit.}: Fp12[BN254_Snarks] + x10100 {.noInit.}: Fp12[BN254_Snarks] + x11001 {.noInit.}: Fp12[BN254_Snarks] + x100010 {.noInit.}: Fp12[BN254_Snarks] + x100111 {.noInit.}: Fp12[BN254_Snarks] + x101001 {.noInit.}: Fp12[BN254_Snarks] - x10 .cyclotomic_square(a) - x11 .prod(x10, a) - x100 .prod(x11, a) - x110 .prod(x10, x100) - x1100 .cyclotomic_square(x110) - x1111 .prod(x11, x1100) - x10010 .prod(x11, x1111) - x10110 .prod(x100, x10010) - x11100 .prod(x110, x10110) - x101110 .prod(x10010, x11100) - x1001010 .prod(x11100, x101110) - x1111000 .prod(x101110, x1001010) - x10001110 .prod(x10110, x1111000) + x10 .cyclotomic_square(a) + x100 .cyclotomic_square(x10) + x1000 .cyclotomic_square(x100) + x10000 .cyclotomic_square(x1000) + x10001 .prod(x10000, a) + x10011 .prod(x10001, x10) + x10100 .prod(x10011, a) + x11001 .prod(x1000, x10001) + x100010 .cyclotomic_square(x10001) + x100111 .prod(x10011, x10100) + x101001 .prod(x10, x100111) - var - r15 {.noInit.}: Fp12[BN254_Snarks] - r16 {.noInit.}: Fp12[BN254_Snarks] - r17 {.noInit.}: Fp12[BN254_Snarks] - r18 {.noInit.}: Fp12[BN254_Snarks] - r20 {.noInit.}: Fp12[BN254_Snarks] - r21 {.noInit.}: Fp12[BN254_Snarks] - r22 {.noInit.}: Fp12[BN254_Snarks] - r26 {.noInit.}: Fp12[BN254_Snarks] - r27 {.noInit.}: Fp12[BN254_Snarks] - r61 {.noInit.}: Fp12[BN254_Snarks] + r.cycl_sqr_repeated(x100010, 6) + r *= x100 + r *= x11001 + r.cycl_sqr_repeated(7) + r *= x11001 - r15.cyclotomic_square(x10001110) - r15 *= x1001010 - r16.prod(x10001110, r15) - r17.prod(x1111, r16) - r18.prod(r16, r17) + r.cycl_sqr_repeated(8) + r *= x101001 + r *= x10 + r.cycl_sqr_repeated(6) + r *= x10001 - r20.cyclotomic_square(r18) - r20 *= r17 - r21.prod(x1111000, r20) - r22.prod(r15, r21) + r.cycl_sqr_repeated(8) + r *= x101001 + r.cycl_sqr_repeated(6) + r *= x101001 + r.cycl_sqr_repeated(10) - r26.cyclotomic_square(r22) - r26.cyclotomic_square() - r26 *= r22 - r26 *= r18 + r *= x100111 + r.cycl_sqr_repeated(6) + r *= x101001 + r *= x1000 - r27.prod(r22, r26) - - r61.prod(r26, r27) - r61.cycl_sqr_repeated(17) - r61 *= r27 - r61.cycl_sqr_repeated(14) - r61 *= r21 - - r = r61 - r.cycl_sqr_repeated(16) - r *= r20 - - if invert: - r.cyclotomic_inv() + if invert: + r.cyclotomic_inv() func isInPairingSubgroup*(a: Fp12[BN254_Snarks]): SecretBool = ## Returns true if a is in GT subgroup, i.e. a is an element of order r @@ -132,8 +109,8 @@ func isInPairingSubgroup*(a: Fp12[BN254_Snarks]): SecretBool = # on BLS pairing-friendly curves # P is in the G1 subgroup iff a^p == a^(6u²) var t0{.noInit.}, t1{.noInit.}: Fp12[BN254_Snarks] - t0.pow_u(a) # a^p - t1.pow_u(t0) # a^(p²) + t0.cycl_exp_by_curve_param(a) # a^p + t1.cycl_exp_by_curve_param(t0) # a^(p²) t0.square(t1) # a^(2p²) t0 *= t1 # a^(3p²) t0.square() # a^(6p²) diff --git a/constantine/curves/bw6_761_pairing.nim b/constantine/curves/bw6_761_pairing.nim index e14434b..989ecba 100644 --- a/constantine/curves/bw6_761_pairing.nim +++ b/constantine/curves/bw6_761_pairing.nim @@ -10,7 +10,7 @@ import ../config/[common, curves, type_bigint], ../io/io_bigints, ../towers, - ../pairing/cyclotomic_fp6, + ../pairing/cyclotomic_subgroup, ../isogeny/frobenius # Slow generic implementation @@ -64,7 +64,7 @@ const BW6_761_pairing_finalexponent_hard* = block: # Addition chain # ------------------------------------------------------------ -func pow_u*(r: var Fp6[BW6_761], a: Fp6[BW6_761], invert = false) = +func cycl_exp_by_curve_param*(r: var Fp6[BW6_761], a: Fp6[BW6_761], invert = false) = ## f^u with u the curve parameter ## For BLS12_377/BW6_761 f^0x8508c00000000001 r.cyclotomic_square(a) @@ -111,8 +111,8 @@ func isInPairingSubgroup*(a: Fp6[BW6_761]): SecretBool = u0.cyclotomic_square() # u0 = 18 u0 *= a # u0 = 19 - u1.pow_u(a) # u1 = u - u4.pow_u(u1) # u4 = u² + u1.cycl_exp_by_curve_param(a) # u1 = u + u4.cycl_exp_by_curve_param(u1) # u4 = u² u3.cyclotomic_square(u1) # u3 = 2u u3.cyclotomic_square() # u3 = 4u @@ -121,7 +121,7 @@ func isInPairingSubgroup*(a: Fp6[BW6_761]): SecretBool = u0 *= u1 # u0 = 10u + 19 - u1.pow_u(u4) # u1 = u³ + u1.cycl_exp_by_curve_param(u4) # u1 = u³ u3.cyclotomic_square(u1) # u3 = 2u³ u3.cycl_sqr_repeated(3) # u3 = 16u³ u3 *= u1 # u3 = 17u³ @@ -129,9 +129,9 @@ func isInPairingSubgroup*(a: Fp6[BW6_761]): SecretBool = u3 *= u1 # u3 = 35u³ u0 *= u3 # u0 = 35u³ + 10u + 19 - u4.pow_u(u1) # u4 = u⁴ - u5.pow_u(u4) # u5 = u⁵ - u6.pow_u(u5) # u6 = u⁶ + u4.cycl_exp_by_curve_param(u1) # u4 = u⁴ + u5.cycl_exp_by_curve_param(u4) # u5 = u⁵ + u6.cycl_exp_by_curve_param(u5) # u6 = u⁶ u1.cyclotomic_square(u4) # u1 = 2u⁴ u1.cycl_sqr_repeated(2) # u1 = 8u⁴ diff --git a/constantine/curves/zoo_pairings.nim b/constantine/curves/zoo_pairings.nim index 077f654..cd97162 100644 --- a/constantine/curves/zoo_pairings.nim +++ b/constantine/curves/zoo_pairings.nim @@ -21,4 +21,4 @@ macro pairing*(C: static Curve, value: untyped): untyped = ## Get pairing related constants return bindSym($C & "_pairing_" & $value) -export pow_x, pow_xdiv2, pow_u, millerLoopAddchain, isInPairingSubgroup +export cycl_exp_by_curve_param, cycl_exp_by_curve_param_div2, millerLoopAddchain, isInPairingSubgroup diff --git a/constantine/io/io_ec.nim b/constantine/io/io_ec.nim index 3352b43..8203056 100644 --- a/constantine/io/io_ec.nim +++ b/constantine/io/io_ec.nim @@ -27,11 +27,6 @@ import # # ############################################################ -proc spaces(num: int): string = - result = newString(num) - for i in 0 ..< num: - result[i] = ' ' - func toHex*[EC: ECP_ShortW_Prj or ECP_ShortW_Jac or ECP_ShortW_Aff](P: EC, indent: static int = 0): string = ## Stringify an elliptic curve point to Hex ## Note. Leading zeros are not removed. @@ -53,9 +48,9 @@ func toHex*[EC: ECP_ShortW_Prj or ECP_ShortW_Jac or ECP_ShortW_Aff](P: EC, inden const sp = spaces(indent) result = sp & $EC & "(\n" & sp & " x: " - result.appendHex(aff.x, bigEndian) + result.appendHex(aff.x) result &= ",\n" & sp & " y: " - result.appendHex(aff.y, bigEndian) + result.appendHex(aff.y) result &= "\n" & sp & ")" func fromHex*(dst: var (ECP_ShortW_Prj or ECP_ShortW_Jac), x, y: string): bool {.raises: [ValueError].}= diff --git a/constantine/io/io_towers.nim b/constantine/io/io_towers.nim index 00131f2..280f1d0 100644 --- a/constantine/io/io_towers.nim +++ b/constantine/io/io_towers.nim @@ -28,17 +28,25 @@ export extension_fields # # ############################################################ -func appendHex*(accum: var string, f: Fp2 or Fp4 or Fp6 or Fp12, order: static Endianness = bigEndian) = +proc spaces*(num: int): string = + result = newString(num) + for i in 0 ..< num: + result[i] = ' ' + +func appendHex*(accum: var string, f: Fp2 or Fp4 or Fp6 or Fp12, indent = 0, order: static Endianness = bigEndian) = ## Hex accumulator accum.add static($f.typeof.genericHead() & '(') staticFor i, 0, f.coords.len: when i != 0: accum.add ", " - accum.add "c" & $i & ": " - accum.appendHex(f.coords[i], order) + accum.add "\n" & spaces(indent+2) & "c" & $i & ": " + when f is Fp2: + accum.appendHex(f.coords[i], order) + else: + accum.appendHex(f.coords[i], indent+2, order) accum.add ")" -func toHex*(f: Fp2 or Fp4 or Fp6 or Fp12, order: static Endianness = bigEndian): string = +func toHex*(f: Fp2 or Fp4 or Fp6 or Fp12, indent = 0, order: static Endianness = bigEndian): string = ## Stringify a tower field element to hex. ## Note. Leading zeros are not removed. ## Result is prefixed with 0x @@ -47,7 +55,7 @@ func toHex*(f: Fp2 or Fp4 or Fp6 or Fp12, order: static Endianness = bigEndian): ## ## CT: ## - no leaks - result.appendHex(f, order) + result.appendHex(f, indent, order) func fromHex*(dst: var Fp2, c0, c1: string) {.raises: [ValueError].}= ## Convert 2 coordinates to an element of 𝔽p2 diff --git a/constantine/pairing/cyclotomic_fp12.nim b/constantine/pairing/cyclotomic_fp12.nim deleted file mode 100644 index 8d2a142..0000000 --- a/constantine/pairing/cyclotomic_fp12.nim +++ /dev/null @@ -1,233 +0,0 @@ -# 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 - ../primitives, - ../config/[common, curves], - ../arithmetic, - ../towers, - ../isogeny/frobenius - -# No exceptions allowed -{.push raises: [].} - -# ############################################################ -# -# Gϕ₁₂, Cyclotomic subgroup of Fp12 -# with GΦₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) ≡ 1 (mod pⁿ)} -# -# ############################################################ - -# - Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions -# Granger, Scott, 2009 -# https://eprint.iacr.org/2009/565.pdf -# -# - On the final exponentiation for calculating -# pairings on ordinary elliptic curves -# Scott, Benger, Charlemagne, Perez, Kachisa, 2008 -# https://eprint.iacr.org/2008/490.pdf - -# 𝔽p12 -> Gϕ₁₂ - Mapping to Cyclotomic group -# ---------------------------------------------------------------- -func finalExpEasy*[C: static Curve](f: var Fp12[C]) {.meter.} = - ## Easy part of the final exponentiation - ## - ## This maps the result of the Miller loop into the cyclotomic subgroup Gϕ₁₂ - ## - ## We need to clear the Gₜ cofactor to obtain - ## an unique Gₜ representation - ## (reminder, Gₜ is a multiplicative group hence we exponentiate by the cofactor) - ## - ## i.e. Fp¹² --> (fexp easy) --> Gϕ₁₂ --> (fexp hard) --> Gₜ - ## - ## The final exponentiation is fexp = f^((p¹² - 1) / r) - ## It is separated into: - ## f^((p¹² - 1) / r) = (p¹² - 1) / ϕ₁₂(p) * ϕ₁₂(p) / r - ## - ## with the cyclotomic polynomial ϕ₁₂(p) = (p⁴-p²+1) - ## - ## With an embedding degree of 12, the easy part of final exponentiation is - ## - ## f^(p⁶−1)(p²+1) - ## - ## And properties are - ## 0. f^(p⁶) ≡ conj(f) (mod p¹²) for all f in Fp12 - ## - ## After g = f^(p⁶−1) the result g is on the cyclotomic subgroup - ## 1. g^(-1) ≡ g^(p⁶) (mod p¹²) - ## 2. Inversion can be done with conjugate - ## 3. g is unitary, its norm |g| (the product of conjugates) is 1 - ## 4. Squaring has a fast compressed variant. - # - # Proofs: - # - # Fp12 can be defined as a quadratic extension over Fp⁶ - # with g = g₀ + x g₁ with x a quadratic non-residue - # - # with q = p⁶, q² = p¹² - # The frobenius map f^q ≡ (f₀ + x f₁)^q (mod q²) - # ≡ f₀^q + x^q f₁^q (mod q²) - # ≡ f₀ + x^q f₁ (mod q²) - # ≡ f₀ - x f₁ (mod q²) - # hence - # f^p⁶ ≡ conj(f) (mod p¹²) - # Q.E.D. of (0) - # - # ---------------- - # - # p¹² - 1 = (p⁶−1)(p⁶+1) = (p⁶−1)(p²+1)(p⁴-p²+1) - # by Fermat's little theorem we have - # f^(p¹² - 1) ≡ 1 (mod p¹²) - # - # Hence f^(p⁶−1)(p⁶+1) ≡ 1 (mod p¹²) - # - # We call g = f^(p⁶−1) we have - # g^(p⁶+1) ≡ 1 (mod p¹²) <=> g^(p⁶) * g ≡ 1 (mod p¹²) - # hence g^(-1) ≡ g^(p⁶) (mod p¹²) - # Q.E.D. of (1) - # - # -- - # - # From (1) g^(-1) ≡ g^(p⁶) (mod p¹²) for g = f^(p⁶−1) - # and (0) f^p⁶ ≡ conj(f) (mod p¹²) for all f in fp12 - # - # so g^(-1) ≡ conj(g) (mod p¹²) for g = f^(p⁶−1) - # Q.E.D. of (2) - # - # -- - # - # f^(p¹² - 1) ≡ 1 (mod p¹²) by Fermat's Little Theorem - # f^(p⁶−1)(p⁶+1) ≡ 1 (mod p¹²) - # g^(p⁶+1) ≡ 1 (mod p¹²) - # g * g^p⁶ ≡ 1 (mod p¹²) - # g * conj(g) ≡ 1 (mod p¹²) - # Q.E.D. of (3) - var g {.noinit.}: typeof(f) - g.inv(f) # g = f^-1 - conj(f) # f = f^p⁶ - g *= f # g = f^(p⁶-1) - f.frobenius_map(g, 2) # f = f^((p⁶-1) p²) - f *= g # f = f^((p⁶-1) (p²+1)) - -# Gϕ₁₂ - Cyclotomic functions -# ---------------------------------------------------------------- -# A cyclotomic group is a subgroup of Fp^n defined by -# -# GΦₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) = 1} -# -# The result of any pairing is in a cyclotomic subgroup - -func cyclotomic_inv*(a: var Fp12) {.meter.} = - ## Fast inverse for a - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - a.conj() - -func cyclotomic_inv*(r: var Fp12, a: Fp12) {.meter.} = - ## Fast inverse for a - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - r.conj(a) - -func cyclotomic_square*[C](r: var Fp12[C], a: Fp12[C]) {.meter.} = - ## Square `a` into `r` - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - # - # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions - # Granger, Scott, 2009 - # https://eprint.iacr.org/2009/565.pdf - - when a.c0 is Fp4: - # Cubic over quadratic - # A = 3a² − 2 ̄a - # B = 3 √i c² + 2 ̄b - # C = 3b² − 2 ̄c - var t0{.noinit.}, t1{.noinit.}: Fp4[C] - - t0.square(a.c0) # t0 = a² - t1.double(t0) # t1 = 2a² - t1 += t0 # t1 = 3a² - - t0.conj(a.c0) # t0 = ̄a - t0.double() # t0 = 2 ̄a - r.c0.diff(t1, t0) # r0 = 3a² − 2 ̄a - - # Aliasing: a.c0 unused - - t0.square(a.c2) # t0 = c² - t0 *= NonResidue # t0 = √i c² - t1.double(t0) # t1 = 2 √i c² - t0 += t1 # t0 = 3 √i c² - - t1.square(a.c1) # t1 = b² - - r.c1.conj(a.c1) # r1 = ̄b - r.c1.double() # r1 = 2 ̄b - r.c1 += t0 # r1 = 3 √i c² + 2 ̄b - - # Aliasing: a.c1 unused - - t0.double(t1) # t0 = 2b² - t0 += t1 # t0 = 3b² - - t1.conj(a.c2) # r2 = ̄c - t1.double() # r2 = 2 ̄c - r.c2.diff(t0, t1) # r2 = 3b² - 2 ̄c - - else: - {.error: "Not implemented".} - -func cyclotomic_square*[C](a: var Fp12[C]) {.meter.} = - ## Square `a` into `r` - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - # - # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions - # Granger, Scott, 2009 - # https://eprint.iacr.org/2009/565.pdf - a.cyclotomic_square(a) - -func cycl_sqr_repeated*(f: var Fp12, num: int) {.inline, meter.} = - ## Repeated cyclotomic squarings - for _ in 0 ..< num: - f.cyclotomic_square() - -iterator unpack(scalarByte: byte): bool = - yield bool((scalarByte and 0b10000000) shr 7) - yield bool((scalarByte and 0b01000000) shr 6) - yield bool((scalarByte and 0b00100000) shr 5) - yield bool((scalarByte and 0b00010000) shr 4) - yield bool((scalarByte and 0b00001000) shr 3) - yield bool((scalarByte and 0b00000100) shr 2) - yield bool((scalarByte and 0b00000010) shr 1) - yield bool( scalarByte and 0b00000001) - -func cyclotomic_exp*[C](r: var Fp12[C], a: Fp12[C], exponent: BigInt, invert: bool) {.meter.} = - var eBytes: array[(exponent.bits+7) div 8, byte] - eBytes.exportRawUint(exponent, bigEndian) - - r.setOne() - for b in eBytes: - for bit in unpack(b): - r.cyclotomic_square() - if bit: - r *= a - if invert: - r.cyclotomic_inv() - -func isInCyclotomicSubgroup*[C](a: Fp12[C]): SecretBool = - ## Check if a ∈ Fpⁿ: a^Φₙ(p) = 1 - ## Φ₁₂(p) = p⁴-p²+1 - var t{.noInit.}, p2{.noInit.}: Fp12[C] - - p2.frobenius_map(a, 2) # a^(p²) - t.frobenius_map(p2, 2) # a^(p⁴) - t *= a # a^(p⁴+1) - - return t == p2 \ No newline at end of file diff --git a/constantine/pairing/cyclotomic_fp6.nim b/constantine/pairing/cyclotomic_fp6.nim deleted file mode 100644 index fbe6476..0000000 --- a/constantine/pairing/cyclotomic_fp6.nim +++ /dev/null @@ -1,233 +0,0 @@ -# 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 - ../primitives, - ../config/[common, curves], - ../arithmetic, - ../towers, - ../isogeny/frobenius - -# No exceptions allowed -{.push raises: [].} - -# ############################################################ -# -# Gϕ₆, Cyclotomic subgroup of Fp6 -# with GΦₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) ≡ 1 (mod pⁿ)} -# -# ############################################################ - -# - Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions -# Granger, Scott, 2009 -# https://eprint.iacr.org/2009/565.pdf -# -# - On the final exponentiation for calculating -# pairings on ordinary elliptic curves -# Scott, Benger, Charlemagne, Perez, Kachisa, 2008 -# https://eprint.iacr.org/2008/490.pdf - -# 𝔽p6 -> Gϕ₆ - Mapping to Cyclotomic group -# ---------------------------------------------------------------- -func finalExpEasy*[C: static Curve](f: var Fp6[C]) {.meter.} = - ## Easy part of the final exponentiation - ## - ## This maps the result of the Miller loop into the cyclotomic subgroup Gϕ₆ - ## - ## We need to clear the Gₜ cofactor to obtain - ## an unique Gₜ representation - ## (reminder, Gₜ is a multiplicative group hence we exponentiate by the cofactor) - ## - ## i.e. Fp⁶ --> (fexp easy) --> Gϕ₆ --> (fexp hard) --> Gₜ - ## - ## The final exponentiation is fexp = f^((p⁶ - 1) / r) - ## It is separated into: - ## f^((p⁶ - 1) / r) = (p⁶ - 1) / ϕ₆(p) * ϕ₆(p) / r - ## - ## with the cyclotomic polynomial ϕ₆(p) = (p²-p+1) - ## - ## With an embedding degree of 6, the easy part of final exponentiation is - ## - ## f^(p³−1)(p+1) - ## - ## And properties are - ## 0. f^(p³) ≡ conj(f) (mod p⁶) for all f in Fp6 - ## - ## After g = f^(p³−1) the result g is on the cyclotomic subgroup - ## 1. g^(-1) ≡ g^(p³) (mod p⁶) - ## 2. Inversion can be done with conjugate - ## 3. g is unitary, its norm |g| (the product of conjugates) is 1 - ## 4. Squaring has a fast compressed variant. - # - # Proofs: - # - # Fp6 can be defined as a quadratic extension over Fp³ - # with g = g₀ + x g₁ with x a quadratic non-residue - # - # with q = p³, q² = p⁶ - # The frobenius map f^q ≡ (f₀ + x f₁)^q (mod q²) - # ≡ f₀^q + x^q f₁^q (mod q²) - # ≡ f₀ + x^q f₁ (mod q²) - # ≡ f₀ - x f₁ (mod q²) - # hence - # f^p³ ≡ conj(f) (mod p⁶) - # Q.E.D. of (0) - # - # ---------------- - # - # p⁶ - 1 = (p³−1)(p³+1) = (p³−1)(p+1)(p²-p+1) - # by Fermat's little theorem we have - # f^(p⁶ - 1) ≡ 1 (mod p⁶) - # - # Hence f^(p³−1)(p³+1) ≡ 1 (mod p⁶) - # - # We call g = f^(p³−1) we have - # g^(p³+1) ≡ 1 (mod p⁶) <=> g^(p³) * g ≡ 1 (mod p⁶) - # hence g^(-1) ≡ g^(p³) (mod p⁶) - # Q.E.D. of (1) - # - # -- - # - # From (1) g^(-1) ≡ g^(p³) (mod p⁶) for g = f^(p³−1) - # and (0) f^p³ ≡ conj(f) (mod p⁶) for all f in fp12 - # - # so g^(-1) ≡ conj(g) (mod p⁶) for g = f^(p³−1) - # Q.E.D. of (2) - # - # -- - # - # f^(p⁶ - 1) ≡ 1 (mod p⁶) by Fermat's Little Theorem - # f^(p³−1)(p³+1) ≡ 1 (mod p⁶) - # g^(p³+1) ≡ 1 (mod p⁶) - # g * g^p³ ≡ 1 (mod p⁶) - # g * conj(g) ≡ 1 (mod p⁶) - # Q.E.D. of (3) - var g {.noinit.}: typeof(f) - g.inv(f) # g = f^-1 - conj(f) # f = f^p³ - g *= f # g = f^(p³-1) - f.frobenius_map(g) # f = f^((p³-1) p) - f *= g # f = f^((p³-1) (p+1)) - -# Gϕ₆ - Cyclotomic functions -# ---------------------------------------------------------------- -# A cyclotomic group is a subgroup of Fp^n defined by -# -# GΦₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) = 1} -# -# The result of any pairing is in a cyclotomic subgroup - -func cyclotomic_inv*(a: var Fp6) {.meter.} = - ## Fast inverse for a - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - a.conj() - -func cyclotomic_inv*(r: var Fp6, a: Fp6) {.meter.} = - ## Fast inverse for a - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - r.conj(a) - -func cyclotomic_square*[C](r: var Fp6[C], a: Fp6[C]) {.meter.} = - ## Square `a` into `r` - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - # - # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions - # Granger, Scott, 2009 - # https://eprint.iacr.org/2009/565.pdf - - when a.c0 is Fp2: - # Cubic over quadratic - # A = 3a² − 2 ̄a - # B = 3 √i c² + 2 ̄b - # C = 3b² − 2 ̄c - var t0{.noinit.}, t1{.noinit.}: Fp2[C] - - t0.square(a.c0) # t0 = a² - t1.double(t0) # t1 = 2a² - t1 += t0 # t1 = 3a² - - t0.conj(a.c0) # t0 = ̄a - t0.double() # t0 = 2 ̄a - r.c0.diff(t1, t0) # r0 = 3a² − 2 ̄a - - # Aliasing: a.c0 unused - - t0.square(a.c2) # t0 = c² - t0 *= NonResidue # t0 = √i c² - t1.double(t0) # t1 = 2 √i c² - t0 += t1 # t0 = 3 √i c² - - t1.square(a.c1) # t1 = b² - - r.c1.conj(a.c1) # r1 = ̄b - r.c1.double() # r1 = 2 ̄b - r.c1 += t0 # r1 = 3 √i c² + 2 ̄b - - # Aliasing: a.c1 unused - - t0.double(t1) # t0 = 2b² - t0 += t1 # t0 = 3b² - - t1.conj(a.c2) # r2 = ̄c - t1.double() # r2 = 2 ̄c - r.c2.diff(t0, t1) # r2 = 3b² - 2 ̄c - - else: - {.error: "Not implemented".} - -func cyclotomic_square*[C](a: var Fp6[C]) {.meter.} = - ## Square `a` into `r` - ## `a` MUST be in the cyclotomic subgroup - ## consequently `a` MUST be unitary - # - # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions - # Granger, Scott, 2009 - # https://eprint.iacr.org/2009/565.pdf - a.cyclotomic_square(a) - -func cycl_sqr_repeated*(f: var Fp6, num: int) {.inline, meter.} = - ## Repeated cyclotomic squarings - for _ in 0 ..< num: - f.cyclotomic_square() - -iterator unpack(scalarByte: byte): bool = - yield bool((scalarByte and 0b10000000) shr 7) - yield bool((scalarByte and 0b01000000) shr 6) - yield bool((scalarByte and 0b00100000) shr 5) - yield bool((scalarByte and 0b00010000) shr 4) - yield bool((scalarByte and 0b00001000) shr 3) - yield bool((scalarByte and 0b00000100) shr 2) - yield bool((scalarByte and 0b00000010) shr 1) - yield bool( scalarByte and 0b00000001) - -func cyclotomic_exp*[C](r: var Fp6[C], a: Fp6[C], exponent: BigInt, invert: bool) {.meter.} = - var eBytes: array[(exponent.bits+7) div 8, byte] - eBytes.exportRawUint(exponent, bigEndian) - - r.setOne() - for b in eBytes: - for bit in unpack(b): - r.cyclotomic_square() - if bit: - r *= a - if invert: - r.cyclotomic_inv() - -func isInCyclotomicSubgroup*[C](a: Fp6[C]): SecretBool = - ## Check if a ∈ Fpⁿ: a^Φₙ(p) = 1 - ## Φ₆(p) = p⁴-p²+1 - var t{.noInit.}, p{.noInit.}: Fp6[C] - - t.frobenius_map(a, 2) # a^(p²) - t *= a # a^(p²+1) - p.frobenius_map(a) # a^(p) - - return t == p \ No newline at end of file diff --git a/constantine/pairing/cyclotomic_subgroup.nim b/constantine/pairing/cyclotomic_subgroup.nim new file mode 100644 index 0000000..fed9b05 --- /dev/null +++ b/constantine/pairing/cyclotomic_subgroup.nim @@ -0,0 +1,680 @@ +# 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 + ../primitives, + ../config/[common, curves], + ../arithmetic, + ../towers, + ../isogeny/frobenius + +# No exceptions allowed +{.push raises: [].} + +# ############################################################ +# +# Gϕₙ, Cyclotomic subgroup of Fpⁿ +# with Gϕₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) ≡ 1 (mod pⁿ)} +# +# ############################################################ + +# - Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions +# Granger, Scott, 2009 +# https://eprint.iacr.org/2009/565.pdf +# +# - On the final exponentiation for calculating +# pairings on ordinary elliptic curves +# Scott, Benger, Charlemagne, Perez, Kachisa, 2008 +# https://eprint.iacr.org/2008/490.pdf + +# 𝔽pⁿ -> Gϕₙ - Mapping to Cyclotomic group +# ---------------------------------------------------------------- + +func finalExpEasy*[C: static Curve](f: var Fp6[C]) {.meter.} = + ## Easy part of the final exponentiation + ## + ## This maps the result of the Miller loop into the cyclotomic subgroup Gϕ₆ + ## + ## We need to clear the Gₜ cofactor to obtain + ## an unique Gₜ representation + ## (reminder, Gₜ is a multiplicative group hence we exponentiate by the cofactor) + ## + ## i.e. Fp⁶ --> (fexp easy) --> Gϕ₆ --> (fexp hard) --> Gₜ + ## + ## The final exponentiation is fexp = f^((p⁶ - 1) / r) + ## It is separated into: + ## f^((p⁶ - 1) / r) = (p⁶ - 1) / ϕ₆(p) * ϕ₆(p) / r + ## + ## with the cyclotomic polynomial ϕ₆(p) = (p²-p+1) + ## + ## With an embedding degree of 6, the easy part of final exponentiation is + ## + ## f^(p³−1)(p+1) + ## + ## And properties are + ## 0. f^(p³) ≡ conj(f) (mod p⁶) for all f in Fp6 + ## + ## After g = f^(p³−1) the result g is on the cyclotomic subgroup + ## 1. g^(-1) ≡ g^(p³) (mod p⁶) + ## 2. Inversion can be done with conjugate + ## 3. g is unitary, its norm |g| (the product of conjugates) is 1 + ## 4. Squaring has a fast compressed variant. + # + # Proofs: + # + # Fp6 can be defined as a quadratic extension over Fp³ + # with g = g₀ + x g₁ with x a quadratic non-residue + # + # with q = p³, q² = p⁶ + # The frobenius map f^q ≡ (f₀ + x f₁)^q (mod q²) + # ≡ f₀^q + x^q f₁^q (mod q²) + # ≡ f₀ + x^q f₁ (mod q²) + # ≡ f₀ - x f₁ (mod q²) + # hence + # f^p³ ≡ conj(f) (mod p⁶) + # Q.E.D. of (0) + # + # ---------------- + # + # p⁶ - 1 = (p³−1)(p³+1) = (p³−1)(p+1)(p²-p+1) + # by Fermat's little theorem we have + # f^(p⁶ - 1) ≡ 1 (mod p⁶) + # + # Hence f^(p³−1)(p³+1) ≡ 1 (mod p⁶) + # + # We call g = f^(p³−1) we have + # g^(p³+1) ≡ 1 (mod p⁶) <=> g^(p³) * g ≡ 1 (mod p⁶) + # hence g^(-1) ≡ g^(p³) (mod p⁶) + # Q.E.D. of (1) + # + # -- + # + # From (1) g^(-1) ≡ g^(p³) (mod p⁶) for g = f^(p³−1) + # and (0) f^p³ ≡ conj(f) (mod p⁶) for all f in fp12 + # + # so g^(-1) ≡ conj(g) (mod p⁶) for g = f^(p³−1) + # Q.E.D. of (2) + # + # -- + # + # f^(p⁶ - 1) ≡ 1 (mod p⁶) by Fermat's Little Theorem + # f^(p³−1)(p³+1) ≡ 1 (mod p⁶) + # g^(p³+1) ≡ 1 (mod p⁶) + # g * g^p³ ≡ 1 (mod p⁶) + # g * conj(g) ≡ 1 (mod p⁶) + # Q.E.D. of (3) + var g {.noinit.}: typeof(f) + g.inv(f) # g = f^-1 + conj(f) # f = f^p³ + g *= f # g = f^(p³-1) + f.frobenius_map(g) # f = f^((p³-1) p) + f *= g # f = f^((p³-1) (p+1)) + +func finalExpEasy*[C: static Curve](f: var Fp12[C]) {.meter.} = + ## Easy part of the final exponentiation + ## + ## This maps the result of the Miller loop into the cyclotomic subgroup Gϕ₁₂ + ## + ## We need to clear the Gₜ cofactor to obtain + ## an unique Gₜ representation + ## (reminder, Gₜ is a multiplicative group hence we exponentiate by the cofactor) + ## + ## i.e. Fp¹² --> (fexp easy) --> Gϕ₁₂ --> (fexp hard) --> Gₜ + ## + ## The final exponentiation is fexp = f^((p¹² - 1) / r) + ## It is separated into: + ## f^((p¹² - 1) / r) = (p¹² - 1) / ϕ₁₂(p) * ϕ₁₂(p) / r + ## + ## with the cyclotomic polynomial ϕ₁₂(p) = (p⁴-p²+1) + ## + ## With an embedding degree of 12, the easy part of final exponentiation is + ## + ## f^(p⁶−1)(p²+1) + ## + ## And properties are + ## 0. f^(p⁶) ≡ conj(f) (mod p¹²) for all f in Fp12 + ## + ## After g = f^(p⁶−1) the result g is on the cyclotomic subgroup + ## 1. g^(-1) ≡ g^(p⁶) (mod p¹²) + ## 2. Inversion can be done with conjugate + ## 3. g is unitary, its norm |g| (the product of conjugates) is 1 + ## 4. Squaring has a fast compressed variant. + # + # Proofs: + # + # Fp12 can be defined as a quadratic extension over Fp⁶ + # with g = g₀ + x g₁ with x a quadratic non-residue + # + # with q = p⁶, q² = p¹² + # The frobenius map f^q ≡ (f₀ + x f₁)^q (mod q²) + # ≡ f₀^q + x^q f₁^q (mod q²) + # ≡ f₀ + x^q f₁ (mod q²) + # ≡ f₀ - x f₁ (mod q²) + # hence + # f^p⁶ ≡ conj(f) (mod p¹²) + # Q.E.D. of (0) + # + # ---------------- + # + # p¹² - 1 = (p⁶−1)(p⁶+1) = (p⁶−1)(p²+1)(p⁴-p²+1) + # by Fermat's little theorem we have + # f^(p¹² - 1) ≡ 1 (mod p¹²) + # + # Hence f^(p⁶−1)(p⁶+1) ≡ 1 (mod p¹²) + # + # We call g = f^(p⁶−1) we have + # g^(p⁶+1) ≡ 1 (mod p¹²) <=> g^(p⁶) * g ≡ 1 (mod p¹²) + # hence g^(-1) ≡ g^(p⁶) (mod p¹²) + # Q.E.D. of (1) + # + # -- + # + # From (1) g^(-1) ≡ g^(p⁶) (mod p¹²) for g = f^(p⁶−1) + # and (0) f^p⁶ ≡ conj(f) (mod p¹²) for all f in fp12 + # + # so g^(-1) ≡ conj(g) (mod p¹²) for g = f^(p⁶−1) + # Q.E.D. of (2) + # + # -- + # + # f^(p¹² - 1) ≡ 1 (mod p¹²) by Fermat's Little Theorem + # f^(p⁶−1)(p⁶+1) ≡ 1 (mod p¹²) + # g^(p⁶+1) ≡ 1 (mod p¹²) + # g * g^p⁶ ≡ 1 (mod p¹²) + # g * conj(g) ≡ 1 (mod p¹²) + # Q.E.D. of (3) + var g {.noinit.}: typeof(f) + g.inv(f) # g = f^-1 + conj(f) # f = f^p⁶ + g *= f # g = f^(p⁶-1) + f.frobenius_map(g, 2) # f = f^((p⁶-1) p²) + f *= g # f = f^((p⁶-1) (p²+1)) + +# Gϕₙ - Cyclotomic functions +# ---------------------------------------------------------------- +# A cyclotomic group is a subgroup of Fpⁿ defined by +# +# Gϕₙ(p) = {α ∈ Fpⁿ : α^Φₙ(p) = 1} +# +# The result of any pairing is in a cyclotomic subgroup + +func cyclotomic_inv*[FT](a: var FT) {.meter.} = + ## Fast inverse for a + ## `a` MUST be in the cyclotomic subgroup + ## consequently `a` MUST be unitary + a.conj() + +func cyclotomic_inv*[FT](r: var FT, a: FT) {.meter.} = + ## Fast inverse for a + ## `a` MUST be in the cyclotomic subgroup + ## consequently `a` MUST be unitary + r.conj(a) + +func cyclotomic_square*[FT](r: var FT, a: FT) {.meter.} = + ## Square `a` into `r` + ## `a` MUST be in the cyclotomic subgroup + ## consequently `a` MUST be unitary + # + # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions + # Granger, Scott, 2009 + # https://eprint.iacr.org/2009/565.pdf + + when a is CubicExt: + # Cubic extension field + # A = 3a² − 2 ̄a + # B = 3 √i c² + 2 ̄b + # C = 3b² − 2 ̄c + var t0{.noinit.}, t1{.noinit.}: typeof(a.c0) + + t0.square(a.c0) # t0 = a² + t1.double(t0) # t1 = 2a² + t1 += t0 # t1 = 3a² + + t0.conj(a.c0) # t0 = ̄a + t0.double() # t0 = 2 ̄a + r.c0.diff(t1, t0) # r0 = 3a² − 2 ̄a + + # Aliasing: a.c0 unused + + t0.square(a.c2) # t0 = c² + t0 *= NonResidue # t0 = √i c² + t1.double(t0) # t1 = 2 √i c² + t0 += t1 # t0 = 3 √i c² + + t1.square(a.c1) # t1 = b² + + r.c1.conj(a.c1) # r1 = ̄b + r.c1.double() # r1 = 2 ̄b + r.c1 += t0 # r1 = 3 √i c² + 2 ̄b + + # Aliasing: a.c1 unused + + t0.double(t1) # t0 = 2b² + t0 += t1 # t0 = 3b² + + t1.conj(a.c2) # r2 = ̄c + t1.double() # r2 = 2 ̄c + r.c2.diff(t0, t1) # r2 = 3b² - 2 ̄c + + else: + {.error: "Not implemented".} + +func cyclotomic_square*[FT](a: var FT) {.meter.} = + ## Square `a` into `r` + ## `a` MUST be in the cyclotomic subgroup + ## consequently `a` MUST be unitary + # + # Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions + # Granger, Scott, 2009 + # https://eprint.iacr.org/2009/565.pdf + a.cyclotomic_square(a) + +func cycl_sqr_repeated*[FT](f: var FT, num: int) {.inline, meter.} = + ## Repeated cyclotomic squarings + for _ in 0 ..< num: + f.cyclotomic_square() + +func cycl_sqr_repeated*[FT](r: var FT, a: FT, num: int) {.inline, meter.} = + ## Repeated cyclotomic squarings + r.cyclotomic_square(a) + for _ in 1 ..< num: + r.cyclotomic_square() + +iterator unpack(scalarByte: byte): bool = + yield bool((scalarByte and 0b10000000) shr 7) + yield bool((scalarByte and 0b01000000) shr 6) + yield bool((scalarByte and 0b00100000) shr 5) + yield bool((scalarByte and 0b00010000) shr 4) + yield bool((scalarByte and 0b00001000) shr 3) + yield bool((scalarByte and 0b00000100) shr 2) + yield bool((scalarByte and 0b00000010) shr 1) + yield bool( scalarByte and 0b00000001) + +func cyclotomic_exp*[FT](r: var FT, a: FT, exponent: BigInt, invert: bool) {.meter.} = + var eBytes: array[(exponent.bits+7) div 8, byte] + eBytes.exportRawUint(exponent, bigEndian) + + r.setOne() + for b in eBytes: + for bit in unpack(b): + r.cyclotomic_square() + if bit: + r *= a + if invert: + r.cyclotomic_inv() + +func isInCyclotomicSubgroup*[C](a: Fp6[C]): SecretBool = + ## Check if a ∈ Fpⁿ: a^Φₙ(p) = 1 + ## Φ₆(p) = p⁴-p²+1 + var t{.noInit.}, p{.noInit.}: Fp6[C] + + t.frobenius_map(a, 2) # a^(p²) + t *= a # a^(p²+1) + p.frobenius_map(a) # a^(p) + + return t == p + +func isInCyclotomicSubgroup*[C](a: Fp12[C]): SecretBool = + ## Check if a ∈ Fpⁿ: a^Φₙ(p) = 1 + ## Φ₁₂(p) = p⁴-p²+1 + var t{.noInit.}, p2{.noInit.}: Fp12[C] + + p2.frobenius_map(a, 2) # a^(p²) + t.frobenius_map(p2, 2) # a^(p⁴) + t *= a # a^(p⁴+1) + + return t == p2 + +# ############################################################ +# +# Compressed representations +# +# ############################################################ +# +# The special structure of cyclotomic subgroup allows compression that +# can lead to faster exponentiation: +# +# - Compression in Finite Fields and Torus-Based Cryptography +# Rubin and Silverberg, 2003 +# https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.90.8087&rep=rep1&type=pdf +# +# - Squaring in Cyclotomic Subgroup +# Karabina, 2012 +# https://www.ams.org/journals/mcom/2013-82-281/S0025-5718-2012-02625-1/S0025-5718-2012-02625-1.pdf +# +# Karabina's formula G₂₃₄₅ has the best squaring/decompression cost ratio. +# From a sextic tower Fpᵏᐟ⁶ -> Fpᵏᐟ³ -> Fpᵏ with quadratic-non-residue u and cubic non-residue v +# α = (a₀+a₁u) + (b₀+b₁u)v + (c₀+c₁u)v² +# Compressed α: C(α) = (b₀+b₁u)v + (c₀+c₁u)v² = (g₂+g₃u)v + (g₄+g₅u)v² +# C(α)² = (h₂+h₃u)v + (h₄+h₅u)v² +# with h₂ = 2(g₂ + 3ξg₄g₅) +# h₃ = 3((g₄+g₅)(g₄+ξg₅) - (ξ+1)g₄g₅) - 2g₃ +# h₄ = 3((g₂+g₃)(g₂+ξg₃) - (ξ+1)g₂g₃) - 2g₄ +# h₅ = 2(g₅ + 3(g₂+g₃)(g₂+ξg₃)) +# +# For decompression we can recompute the missing coefficients +# if g₂ != 0 +# g₁ = (g₅²ξ + 3g₄² - 2g₃)/4g₂ g₀ = (2g₁² + g₂g₅ - 3g₃g₄)ξ + 1 +# if g₂ == 0 +# g₁ = 2g₄g₅/g₃ g₀ = (2g₁² - 3g₃g₄)ξ + 1 + +type G2345*[F] = object + ## Compressed representation of cyclotomic subgroup element of a sextic extension + ## (0 + 0u) + (g₂+g₃u)v + (g₄+g₅u)v² + g2, g3, g4, g5: F + +func cyclotomic_square_compressed*[F](g: var G2345[F]) = + ## Karabina's compressed squaring + ## for sextic extension fields + # C(α)² = (h₂+h₃u)v + (h₄+h₅u)v² + # with + # h₂ = 2(g₂ + 3ξg₄g₅) + # h₃ = 3((g₄+g₅)(g₄+ξg₅) - (ξ+1)g₄g₅) - 2g₃ + # h₄ = 3((g₂+g₃)(g₂+ξg₃) - (ξ+1)g₂g₃) - 2g₄ + # h₅ = 2(g₅ + 3(g₂+g₃)(g₂+ξg₃)) + # (4 mul, theorem 3.2 p561) + # + # or + # h₂ = 2g₂ + 6ξg₄g₅ + # h₃ = 3g₄² + 3g₅²ξ - 2g₃ + # h₄ = 3g₂² + 3g₃²ξ - 2g₄ + # h₅ = 2g₅ + 6g₂g₃ + # (2 mul, 4 sqr, section 5.3 p567) + # + # or + # h₂ = 2g₂ + 3ξ((g₄+g₅)²-g₄²-g₅²) + # h₃ = 3(g₄² + g₅²ξ) - 2g₃ + # h₄ = 3(g₂² + g₃²ξ) - 2g₄ + # h₅ = 2g₅ + 3 ((g₂+g₃)²-g₂²-g₃²) + # (6 sqr) + # + # or with quadratic arithmetic + # (h₂+h₃u) = 3u(g₄+g₅u)² + 2(g₂-g₃u) + # (h₄+h₅u) = 3 (g₂+g₃u)² - 2(g₄-g₅u) + # (2x2mul or 2x3sqr, section 5.3 p567) + var g2g3 {.noInit.} = QuadraticExt[F](coords:[g.g2, g.g3]) + var g4g5 {.noInit.} = QuadraticExt[F](coords:[g.g4, g.g5]) + var h2h3 {.noInit.}, h4h5 {.noInit.}: QuadraticExt[F] + + h2h3.square(g4g5) + h2h3 *= NonResidue + h2h3 *= 3 + + h4h5.square(g2g3) + h4h5 *= 3 + + g2g3.double() + g4g5.double() + + g.g2.sum(h2h3.c0, g2g3.c0) + g.g3.diff(h2h3.c1, g2g3.c1) + g.g4.diff(h4h5.c0, g4g5.c0) + g.g5.sum(h4h5.c1, g4g5.c1) + +func recover_g1*[F](g1_num, g1_den: var F, g: G2345[F]) = + ## Compute g₁ from coordinates g₂g₃g₄g₅ + ## of a cyclotomic subgroup element of a sextic extension field + # if g₂ != 0 + # g₁ = (g₅²ξ + 3g₄² - 2g₃)/4g₂ + # if g₂ == 0 + # g₁ = 2g₄g₅/g₃ + # + # Theorem 3.1, this is well-defined for all + # g in Gϕₙ \ {1} + # if g₂=g₃=0 then g₄=g₅=0 as well + # and g₀ = 1 + let g2NonZero = not g.g2.isZero() + var t{.noInit.}: F + + g1_num = g.g4 + t = g.g5 + t.ccopy(g.g4, g2NonZero) + t *= g1_num # g₄² or g₄g₅ + g1_num = t + g1_num.csub(g.g3, g2NonZero) # g₄²- g₃ + g1_num.double() # 2g₄²-2g₃ or 2g₄g₅ + g1_num.cadd(t, g2NonZero) # 3g₄²-2g₃ or 2g₄g₅ + + t.square(g.g5) + t *= NonResidue + g1_num.cadd(t, g2NonZero) # g₅²ξ + 3g₄² - 2g₃ or 2g₄g₅ + + t.prod(g.g2, 4) + g1_den = g.g3 + g1_den.ccopy(t, g2NonZero) # 4g₂ or g₃ + +func batch_ratio_g1s*[N: static int, F]( + dst: var array[N, F], + src: array[N, tuple[g1_num, g1_den: F]]) = + ## Batch inversion of g₁ + ## returns g1_numᵢ/g1_denᵢ + ## This requires that all g1_den != 0 or all g1_den == 0 + ## which is the case if this is used to implement + ## exponentiation in cyclotomic subgroup. + + # Algorithm: Montgomery's batch inversion + # - Speeding the Pollard and Elliptic Curve Methods of Factorization + # Section 10.3.1 + # Peter L. Montgomery + # https://www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866113-7/S0025-5718-1987-0866113-7.pdf + # - Modern Computer Arithmetic + # Section 2.5.1 Several inversions at once + # Richard P. Brent and Paul Zimmermann + # https://members.loria.fr/PZimmermann/mca/mca-cup-0.5.9.pdf + + dst[0] = src[0].g1_den + for i in 1 ..< N: + dst[i].prod(dst[i-1], src[i].g1_den) + + var accInv{.noInit.}: F + accInv.inv(dst[N-1]) + + for i in countdown(N-1, 1): + # Compute inverse + dst[i].prod(accInv, dst[i-1]) + # Apply it + dst[i] *= src[i].g1_num + # Next iteration + accInv *= src[i].g1_den + + dst[0].prod(accInv, src[0].g1_num) + +func recover_g0*[F]( + g0: var F, g1: F, + g: G2345[F]) = + ## Compute g₀ from coordinates g₁g₂g₃g₄g₅ + ## of a cyclotomic subgroup element of a sextic extension field + var t{.noInit.}: F + + t.square(g1) + g0.prod(g.g3, g.g4) + t -= g0 + t.double() + t -= g0 + g0.prod(g.g2, g.g5) + t += g0 + g0.prod(t, NonResidue) + t.setOne() + g0 += t + +func fromFpk*[Fpkdiv6, Fpk]( + g: var G2345[Fpkdiv6], + a: Fpk) = + ## Convert from a sextic extension to the Karabina g₂₃₄₅ + ## representation. + + # GT representations isomorphisms + # =============================== + # + # Given a sextic twist, we can express all elements in terms of z = SNR¹ᐟ⁶ + # + # The canonical direct sextic representation uses coefficients + # + # c₀ + c₁ z + c₂ z² + c₃ z³ + c₄ z⁴ + c₅ z⁵ + # + # with z = SNR¹ᐟ⁶ + # + # The cubic over quadatric towering + # --------------------------------- + # + # (a₀ + a₁ u) + (a₂ + a₃u) v + (a₄ + a₅u) v² + # + # with u = (SNR)¹ᐟ² and v = z = u¹ᐟ³ = (SNR)¹ᐟ⁶ + # + # The quadratic over cubic towering + # --------------------------------- + # + # (b₀ + b₁x + b₂x²) + (b₃ + b₄x + b₅x²)y + # + # with x = (SNR)¹ᐟ³ and y = z = x¹ᐟ² = (SNR)¹ᐟ⁶ + # + # Mapping between towering schemes + # -------------------------------- + # + # g₂₃₄₅ uses the cubic over quadratic representation hence: + # + # c₀ <=> a₀ <=> b₀ <=> g₀ + # c₁ <=> a₂ <=> b₃ <=> g₂ + # c₂ <=> a₄ <=> b₁ <=> g₄ + # c₃ <=> a₁ <=> b₄ <=> g₁ + # c₄ <=> a₃ <=> b₂ <=> g₃ + # c₅ <=> a₅ <=> b₅ <=> g₅ + # + # See also chapter 6.4 + # - Multiplication and Squaring on Pairing-Friendly Fields + # Augusto Jun Devegili and Colm Ó hÉigeartaigh and Michael Scott and Ricardo Dahab, 2006 + # https://eprint.iacr.org/2006/471 + + when a is CubicExt: + when a.c0 is QuadraticExt: + g.g2 = a.c1.c0 + g.g3 = a.c1.c1 + g.g4 = a.c2.c0 + g.g5 = a.c2.c1 + else: + {.error: "a must be a sextic extension field".} + elif a is QuadraticExt: + when a.c0 is CubicExt: + {.error: "𝔽pᵏᐟ⁶ -> 𝔽pᵏᐟ³ -> 𝔽pᵏ towering (quadratic over cubic) is not implemented.".} + else: + {.error: "a must be a sextic extension field".} + else: + {.error: "𝔽pᵏᐟ⁶ -> 𝔽pᵏ towering (direct sextic) is not implemented.".} + +func asFpk*[Fpkdiv6, Fpk]( + a: var Fpk, + g0, g1: Fpkdiv6, + g: G2345[Fpkdiv6]) = + ## Convert from a sextic extension to the Karabina g₂₃₄₅ + ## representation. + when a is CubicExt: + when a.c0 is QuadraticExt: + a.c0.c0 = g0 + a.c0.c1 = g1 + a.c1.c0 = g.g2 + a.c1.c1 = g.g3 + a.c2.c0 = g.g4 + a.c2.c1 = g.g5 + else: + {.error: "a must be a sextic extension field".} + elif a is QuadraticExt: + when a.c0 is CubicExt: + {.error: "𝔽pᵏᐟ⁶ -> 𝔽pᵏᐟ³ -> 𝔽pᵏ towering (quadratic over cubic) is not implemented.".} + else: + {.error: "a must be a sextic extension field".} + else: + {.error: "𝔽pᵏᐟ⁶ -> 𝔽pᵏ towering (direct sextic) is not implemented.".} + +func cyclotomic_exp_compressed*[N: static int, Fpk]( + r: var Fpk, a: Fpk, + squarings: static array[N, int]) = + ## Exponentiation on the cyclotomic subgroup + ## via compressed repeated squarings + ## Exponentiation is done least-signigicant bits first + ## `squarings` represents the number of squarings + ## to do before the next multiplication. + ## + ## `accumSquarings` stores the accumulated squarings so far + ## iff N != 1 + + type Fpkdiv6 = typeof(a.c0.c0) + + var gs {.noInit.}: array[N, G2345[Fpkdiv6]] + + var g {.noInit.}: G2345[Fpkdiv6] + g.fromFpk(a) + + # Compressed squarings + for i in 0 ..< N: + for j in 0 ..< squarings[i]: + g.cyclotomic_square_compressed() + gs[i] = g + + # Batch decompress + var g1s_ratio {.noInit.}: array[N, tuple[g1_num, g1_den: Fpkdiv6]] + for i in 0 ..< N: + recover_g1(g1s_ratio[i].g1_num, g1s_ratio[i].g1_den, gs[i]) + + var g1s {.noInit.}: array[N, Fpkdiv6] + g1s.batch_ratio_g1s(g1s_ratio) + + var g0s {.noInit.}: array[N, Fpkdiv6] + for i in 0 ..< N: + g0s[i].recover_g0(g1s[i], gs[i]) + + r.asFpk(g0s[0], g1s[0], gs[0]) + for i in 1 ..< N: + var t {.noInit.}: Fpk + t.asFpk(g0s[i], g1s[i], gs[i]) + r *= t + +func cyclotomic_exp_compressed*[N: static int, Fpk]( + r, accumSquarings: var Fpk, a: Fpk, + squarings: static array[N, int]) = + ## Exponentiation on the cyclotomic subgroup + ## via compressed repeated squarings + ## Exponentiation is done least-signigicant bits first + ## `squarings` represents the number of squarings + ## to do before the next multiplication. + ## + ## `accumSquarings` stores the accumulated squarings so far + ## iff N != 1 + + type Fpkdiv6 = typeof(a.c0.c0) + + var gs {.noInit.}: array[N, G2345[Fpkdiv6]] + + var g {.noInit.}: G2345[Fpkdiv6] + g.fromFpk(a) + + # Compressed squarings + for i in 0 ..< N: + for j in 0 ..< squarings[i]: + g.cyclotomic_square_compressed() + gs[i] = g + + # Batch decompress + var g1s_ratio {.noInit.}: array[N, tuple[g1_num, g1_den: Fpkdiv6]] + for i in 0 ..< N: + recover_g1(g1s_ratio[i].g1_num, g1s_ratio[i].g1_den, gs[i]) + + var g1s {.noInit.}: array[N, Fpkdiv6] + g1s.batch_ratio_g1s(g1s_ratio) + + var g0s {.noInit.}: array[N, Fpkdiv6] + for i in 0 ..< N: + g0s[i].recover_g0(g1s[i], gs[i]) + + r.asFpk(g0s[0], g1s[0], gs[0]) + for i in 1 ..< N: + var t {.noInit.}: Fpk + t.asFpk(g0s[i], g1s[i], gs[i]) + r *= t + + if i+1 == N: + accumSquarings = t \ No newline at end of file diff --git a/constantine/pairing/lines_common.nim b/constantine/pairing/lines_common.nim deleted file mode 100644 index 48c2722..0000000 --- a/constantine/pairing/lines_common.nim +++ /dev/null @@ -1,74 +0,0 @@ -# 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, - ../arithmetic, - ../towers, - ../elliptic/ec_shortweierstrass_affine, - ../io/io_towers - -# No exceptions allowed -{.push raises: [].} - -type - Line*[F] = object - ## Packed line representation over a E'(Fpᵏ/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 - - SexticNonResidue* = NonResidue - ## The Sextic non-residue to build - ## 𝔽p2 -> 𝔽p12 towering and the G2 sextic twist - ## or - ## 𝔽p -> 𝔽p6 towering and the G2 sextic twist - ## - ## Note: - ## while the non-residues for - ## - 𝔽p2 -> 𝔽p4 - ## - 𝔽p2 -> 𝔽p6 - ## are also sextic non-residues by construction. - ## the non-residues for - ## - 𝔽p4 -> 𝔽p12 - ## - 𝔽p6 -> 𝔽p12 - ## are not. - -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*[F1, F2](line: var Line[F2], P: ECP_ShortW_Aff[F1, G1]) = - ## Update the line evaluation with P - ## after addition or doubling - ## P in G1 - static: doAssert F1.C == F2.C - line.x *= P.y - line.z *= P.x diff --git a/constantine/pairing/lines_eval.nim b/constantine/pairing/lines_eval.nim new file mode 100644 index 0000000..243b9e3 --- /dev/null +++ b/constantine/pairing/lines_eval.nim @@ -0,0 +1,988 @@ +# 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, + ../config/curves, + ../primitives, + ../arithmetic, + ../towers, + ../elliptic/[ + ec_shortweierstrass_affine, + ec_shortweierstrass_projective + ], + ../io/io_towers + +# No exceptions allowed +{.push raises: [].} + +# ############################################################ +# +# Lines functions +# +# ############################################################ + +# GT representations isomorphisms +# =============================== +# +# Given a sextic twist, we can express all elements in terms of z = SNR¹ᐟ⁶ +# +# The canonical direct sextic representation uses coefficients +# +# c₀ + c₁ z + c₂ z² + c₃ z³ + c₄ z⁴ + c₅ z⁵ +# +# with z = SNR¹ᐟ⁶ +# +# The cubic over quadatric towering +# --------------------------------- +# +# (a₀ + a₁ u) + (a₂ + a₃u) v + (a₄ + a₅u) v² +# +# with u = (SNR)¹ᐟ² and v = z = u¹ᐟ³ = (SNR)¹ᐟ⁶ +# +# The quadratic over cubic towering +# --------------------------------- +# +# (b₀ + b₁x + b₂x²) + (b₃ + b₄x + b₅x²)y +# +# with x = (SNR)¹ᐟ³ and y = z = x¹ᐟ² = (SNR)¹ᐟ⁶ +# +# Mapping between towering schemes +# -------------------------------- +# +# 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 +# Augusto Jun Devegili and Colm Ó hÉigeartaigh and Michael Scott and Ricardo Dahab, 2006 +# https://eprint.iacr.org/2006/471 + +type + Line*[F] = object + ## Packed line representation over a E'(Fpᵏ/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 with GT in 𝔽p12 + ## + ## Out of 6 𝔽p2 coordinates, 3 are zeroes and + ## 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) + ## 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 + + SexticNonResidue* = NonResidue + ## The Sextic non-residue to build + ## 𝔽p2 -> 𝔽p12 towering and the G2 sextic twist + ## or + ## 𝔽p -> 𝔽p6 towering and the G2 sextic twist + ## + ## Note: + ## while the non-residues for + ## - 𝔽p2 -> 𝔽p4 + ## - 𝔽p2 -> 𝔽p6 + ## are also sextic non-residues by construction. + ## the non-residues for + ## - 𝔽p4 -> 𝔽p12 + ## - 𝔽p6 -> 𝔽p12 + ## are not. + +func toHex*(line: Line): string = + result = static($line.typeof.genericHead() & '(') + for fieldName, fieldValue in fieldPairs(line): + when fieldName != "x": + result.add ", " + result.add fieldName & ": " + result.appendHex(fieldValue) + result.add ")" + +# Line evaluation +# ----------------------------------------------------------------------- + +func line_update[F1, F2](line: var Line[F2], P: ECP_ShortW_Aff[F1, G1]) = + ## Update the line evaluation with P + ## after addition or doubling + ## P in G1 + static: doAssert F1.C == F2.C + line.x *= P.y + line.z *= P.x + +# ############################################################ +# +# Miller Loop's Line Evaluation +# with projective coordinates +# +# ############################################################ +# +# - Pairing Implementation Revisited +# Michael Scott, 2019 +# https://eprint.iacr.org/2019/077 +# +# - The Realm of the Pairings +# Diego F. Aranha and Paulo S. L. M. Barreto +# and Patrick Longa and Jefferson E. Ricardini, 2013 +# https://eprint.iacr.org/2013/722.pdf +# http://sac2013.irmacs.sfu.ca/slides/s1.pdf +# +# - Efficient Implementation of Bilinear Pairings on ARM Processors +# Gurleen Grewal, Reza Azarderakhsh, +# Patrick Longa, Shi Hu, and David Jao, 2012 +# https://eprint.iacr.org/2012/408.pdf + +# Line evaluation only +# ----------------------------------------------------------------------------- + +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₂) + +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) + 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.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) + 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 Field.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 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 + +func line_eval_fused_add[Field]( + line: var Line[Field], + T: var ECP_ShortW_Prj[Field, G2], + Q: ECP_ShortW_Aff[Field, G2]) = + ## Fused line evaluation and elliptic point addition + # Grewal et al, 2012 adapted to Scott 2019 line notation + var + A {.noInit.}: Field + B {.noInit.}: Field + C {.noInit.}: Field + D {.noInit.}: Field + E {.noInit.}: Field + F {.noInit.}: Field + G {.noInit.}: Field + H {.noInit.}: Field + I {.noInit.}: Field + + 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.diff(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 Field.C.getSexticTwist() == M_Twist: + lambda *= SexticNonResidue # A = ξ (X₁ - Z₁X₂) + +# Public line evaluation procedures +# ----------------------------------------------------------------------------- + +func line_double*[F1, F2]( + line: var Line[F2], + T: var ECP_ShortW_Prj[F2, G2], + P: ECP_ShortW_Aff[F1, G1]) = + ## Doubling step of the Miller loop + ## T in G2, P in G1 + ## + ## 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() + +func line_add*[F1, F2]( + line: var Line[F2], + T: var ECP_ShortW_Prj[F2, G2], + Q: ECP_ShortW_Aff[F2, G2], + P: ECP_ShortW_Aff[F1, G1]) = + ## Addition step of the Miller loop + ## T and Q in G2, P in G1 + ## + ## 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 + + +# ############################################################ +# +# Sparse Multiplication +# by lines for curves with a sextic twist +# +# ############################################################ + +# - Pairing Implementation Revisited +# Michael Scott, 2019 +# https://eprint.iacr.org/2019/077 +# +# - Efficient Implementation of Bilinear Pairings on ARM Processors +# Gurleen Grewal, Reza Azarderakhsh, +# Patrick Longa, Shi Hu, and David Jao, 2012 +# https://eprint.iacr.org/2012/408.pdf +# +# - High-Speed Software Implementation of the Optimal Ate Pairing over Barreto-Naehrig Curves\ +# Jean-Luc Beuchat and Jorge Enrique González Díaz and Shigeo Mitsunari and Eiji Okamoto and Francisco Rodríguez-Henríquez and Tadanori Teruya, 2010\ +# https://eprint.iacr.org/2010/354 +# +# - Faster Pairing Computations on Curves with High-Degree Twists +# Craig Costello, Tanja Lange, and Michael Naehrig, 2009 +# https://eprint.iacr.org/2009/615.pdf + +# ############################################################ +# +# 𝔽pᵏ by line - 𝔽pᵏ quadratic over cubic +# +# ############################################################ + +# D-Twist +# ------------------------------------------------------------ + +func mul_by_line_xy0*[Fpkdiv2, Fpkdiv6]( + r: var Fpkdiv2, + a: Fpkdiv2, + b: Line[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. + ## `r` and `a` must not alias + + static: + doAssert a is CubicExt + doAssert a.c0 is Fpkdiv6 + + var + 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) + 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.c1 *= r.c2 + r.c1 -= v0 + r.c1 -= v1 + + r.c2.prod(a.c2, b.x) + r.c2 += v1 + +func mul_sparse_by_line_xy00z0*[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. + ## The sparse element is represented by a packed Line type + ## with coordinate (x,y,z) matching 𝔽pᵏ coordinates xy00z0 + + static: + doAssert Fpk.C.getSexticTwist() == D_Twist + doAssert f is QuadraticExt, "This assumes 𝔽pᵏ as a quadratic extension of 𝔽pᵏᐟ²" + doAssert f.c0 is CubicExt, "This assumes 𝔽pᵏᐟ² as a cubic extension of 𝔽pᵏᐟ⁶" + + type Fpkdiv2 = typeof(f.c0) + + var + v0 {.noInit.}: Fpkdiv2 + v1 {.noInit.}: Fpkdiv2 + v2 {.noInit.}: Line[Fpkdiv6] + v3 {.noInit.}: Fpkdiv2 + + v0.mul_by_line_xy0(f.c0, l) + v1.mul_sparse_by_0y0(f.c1, l.z) + + v2.x = l.x + v2.y.sum(l.y, l.z) + f.c1 += f.c0 + v3.mul_by_line_xy0(f.c1, v2) + v3 -= v0 + f.c1.diff(v3, v1) + + v3.c0.prod(v1.c2, SexticNonResidue) + v3.c0 += v0.c0 + v3.c1.sum(v0.c1, v1.c0) + v3.c2.sum(v0.c2, v1.c1) + f.c0 = v3 + +# ############################################################ +# +# 𝔽pᵏ by line - 𝔽pᵏ cubic over quadratic +# +# ############################################################ + +# D-Twist +# ------------------------------------------------------------ + +func mul_sparse_by_line_xyz000*[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. + ## The sparse element is represented by a packed Line type + ## with coordinates (x,y,z) matching 𝔽pᵏ coordinates xyz000 + + static: + doAssert Fpk.C.getSexticTwist() == D_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ᵏᐟ⁶" + + type Fpkdiv3 = typeof(f.c0) + + # In the following equations (taken from cubic extension implementation) + # a = f + # b0 = (x, y) + # b1 = (z, 0) + # b2 = (0, 0) + # + # v0 = a0 b0 = (f00, f01).(x, y) + # v1 = a1 b1 = (f10, f11).(z, 0) + # v2 = a2 b2 = (f20, f21).(0, 0) + # + # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 + # = ξ (a1 b1 + a2 b1 - v1) + v0 + # = ξ a2 b1 + v0 + # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 + # = (a0 + a1) * (b0 + b1) - v0 - v1 + # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 + # = a0 b0 + a2 b0 - v0 + v1 + # = a2 b0 + v1 + + when false: + var b0 {.noInit.}, v0{.noInit.}, v1{.noInit.}, t{.noInit.}: Fp4[C] + + b0.c0 = l.x + b0.c1 = l.y + + v0.prod(f.c0, b0) + v1.mul_sparse_by_x0(f.c1, l.z) + + # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + f.c1 += f.c0 # r1 = a0 + a1 + t = b0 + t.c0 += l.z # t = b0 + b1 + f.c1 *= t # r2 = (a0 + a1)(b0 + b1) + f.c1 -= v0 + 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 *= SexticNonResidue + f.c0 += v0 + + # r2 = a2 b0 + v1 + f.c2 *= b0 + f.c2 += v1 + + else: # Lazy reduction + 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) + + # 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) + 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.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.sum2xMod(f2x, V1) + f.c2.redc2x(f2x) + +func prod_xyz000_xyz000_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) + # a2 = ( 0, 0) + # b0 = (x1, y1) + # b1 = (z1, 0) + # b2 = ( 0, 0) + # + # v0 = a0 b0 = (x0, y0).(x1, y1) + # v1 = a1 b1 = (z0, 0).(z1, 0) + # v2 = a2 b2 = ( 0, 0).( 0, 0) + # + # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 + # = ξ (a1 b1 + a2 b1 - v1) + v0 + # = v0 + # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 + # = (a0 + a1) * (b0 + b1) - v0 - v1 + # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 + # = a0 b0 - v0 + v1 + # = v1 + + static: + doAssert Fpk.C.getSexticTwist() == D_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ᵏᐟ⁶" + + type Fpkdiv3 = typeof(f.c0) + + 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) + + # 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) + f2x.diff2xMod(f2x, V0) + f2x.c0.diff2xMod(f2x.c0, V1) + f.c1.redc2x(f2x) + + # r0 = v0 + f.c0.redc2x(V0) + + # r2 = v1 + f.c2.c0.redc2x(V1) + f.c2.c1.setZero() + +func mul_sparse_by_abcdefghij00*[Fpk]( + a: var Fpk, b: Fpk) = + ## Sparse multiplication of an 𝔽pᵏ element + ## by a sparse 𝔽pᵏ element abcdefghij00 + ## with each representing 𝔽pᵏᐟ⁶ coordinate + + static: + doAssert Fpk.C.getSexticTwist() == D_Twist + doAssert a is CubicExt, "This assumes 𝔽pᵏ as a cubic extension of 𝔽pᵏᐟ³" + doAssert a.c0 is QuadraticExt, "This assumes 𝔽pᵏᐟ³ as a quadratic extension of 𝔽pᵏᐟ⁶" + + type Fpkdiv3 = typeof(a.c0) + + # In the following equations (taken from cubic extension implementation) + # b0 = (b00, b01) + # b1 = (b10, b11) + # b2 = (b20, 0) + # + # v0 = a0 b0 = (f00, f01).(b00, b01) + # v1 = a1 b1 = (f10, f11).(b10, b11) + # v2 = a2 b2 = (f20, f21).(b20, 0) + # + # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ + # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ + # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ + + var V0 {.noInit.}, V1 {.noInit.}, V2 {.noinit.}: doublePrec(Fpkdiv3) + var t0 {.noInit.}, t1 {.noInit.}: Fpkdiv3 + var f2x{.noInit.}, g2x {.noinit.}: doublePrec(Fpkdiv3) + + V0.prod2x(a.c0, b.c0) + V1.prod2x(a.c1, b.c1) + V2.mul2x_sparse_by_x0(a.c2, b.c2) + + # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ + t0.sum(a.c1, a.c2) + t1.c0.sum(b.c1.c0, b.c2.c0) # b₂ = (b20, 0) + f2x.prod2x_disjoint(t0, t1.c0, b.c1.c1) # (a₁ + a₂).(b₁ + b₂) + f2x.diff2xMod(f2x, V1) + f2x.diff2xMod(f2x, V2) + f2x.prod2x(f2x, NonResidue) + f2x.sum2xMod(f2x, V0) + + # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + t0.sum(a.c0, a.c1) + t1.sum(b.c0, b.c1) + g2x.prod2x(t0, t1) + g2x.diff2xMod(g2x, V0) + g2x.diff2xMod(g2x, V1) + + # r₂ = (a₀ + a₂) and (b₀ + b₂) + t0.sum(a.c0, a.c2) + t1.c0.sum(b.c0.c0, b.c2.c0) # b₂ = (b20, 0) + + # Now we are aliasing free + + # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ + a.c0.redc2x(f2x) + + # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ + f2x.prod2x(V2, NonResidue) + g2x.sum2xMod(g2x, f2x) + a.c1.redc2x(g2x) + + # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ + f2x.prod2x_disjoint(t0, t1.c0, b.c0.c1) + f2x.diff2xMod(f2x, V0) + f2x.diff2xMod(f2x, V2) + f2x.sum2xMod(f2x, V1) + a.c2.redc2x(f2x) + +# M-Twist +# ------------------------------------------------------------ + +func mul_sparse_by_line_xy000z*[Fpk, Fpkdiv6]( + f: var Fpk, l: Line[Fpkdiv6]) = + + 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ᵏᐟ⁶" + + type Fpkdiv3 = typeof(f.c0) + + # In the following equations (taken from cubic extension implementation) + # a = f + # b0 = (x, y) + # b1 = (0, 0) + # b2 = (0, z) + # + # v0 = a0 b0 = (f00, f01).(x, y) + # v1 = a1 b1 = (f10, f11).(0, 0) + # v2 = a2 b2 = (f20, f21).(0, z) + # + # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 + # = ξ (a1 b2 + a2 b2 - v2) + v0 + # = ξ a1 b2 + v0 + # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 + # = a0 b0 + a1 b0 - v0 + ξ v2 + # = a1 b0 + ξ v2 + # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 + # = (a0 + a2) * (b0 + b2) - v0 - v2 + + when false: + var b0 {.noInit.}, v0{.noInit.}, v2{.noInit.}, t{.noInit.}: Fpkdiv3 + + b0.c0 = l.x + b0.c1 = l.y + + v0.prod(f.c0, b0) + v2.mul_sparse_by_0y(f.c2, l.z) + + # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + f.c2 += f.c0 # r2 = a0 + a2 + t = b0 + t.c1 += l.z # 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 *= SexticNonResidue + f.c0 += v0 + + # r1 = a1 b0 + ξ v2 + f.c1 *= b0 + v2 *= SexticNonResidue + f.c1 += v2 + + else: # Lazy reduction + 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) + + # 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) + 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.prod2x(f2x, SexticNonResidue) + f2x.sum2xMod(f2x, V0) + f.c0.redc2x(f2x) + + # r1 = a1 b0 + ξ v2 + f2x.prod2x_disjoint(f.c1, l.x, l.y) + 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]) = + ## Multiply 2 lines together + ## The result is sparse in f.c1.c0 + # In the following equations (taken from cubic extension implementation) + # a0 = (x0, y0) + # a1 = ( 0, 0) + # a2 = ( 0, z0) + # b0 = (x1, y1) + # b1 = ( 0, 0) + # b2 = ( 0, z1) + # + # v0 = a0 b0 = (x0, y0).(x1, y1) + # v1 = a1 b1 = ( 0, 0).( 0, 0) + # v2 = a2 b2 = ( 0, z0).( 0, z1) + # + # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 + # = ξ (a1 b2 + a2 b2 - v2) + v0 + # = v0 + # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 + # = a0 b0 + a1 b0 - v0 + ξ v2 + # = ξ v2 + # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 + # = (a0 + a2) * (b0 + b2) - v0 - v2 + + 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ᵏᐟ⁶" + + type Fpkdiv3 = typeof(f.c0) + + 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) + + # 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) + 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) + + # r1 = ξ v2 + f.c1.c1.redc2x(V2) + f.c1.c0.setZero() + + # r0 = v0 + f.c0.redc2x(V0) + +func mul_sparse_by_abcd00efghij*[Fpk]( + a: var Fpk, b: Fpk) = + ## Sparse multiplication of an 𝔽pᵏ element + ## by a sparse 𝔽pᵏ element abcd00efghij + ## with each representing 𝔽pᵏᐟ⁶ coordinate + + static: + doAssert Fpk.C.getSexticTwist() == M_Twist + doAssert a is CubicExt, "This assumes 𝔽pᵏ as a cubic extension of 𝔽pᵏᐟ³" + doAssert a.c0 is QuadraticExt, "This assumes 𝔽pᵏᐟ³ as a quadratic extension of 𝔽pᵏᐟ⁶" + + type Fpkdiv3 = typeof(a.c0) + + # In the following equations (taken from cubic extension implementation) + # b0 = (b00, b01) + # b1 = ( 0, b11) + # b2 = (b20, b21) + # + # v0 = a0 b0 = (f00, f01).(b00, b01) + # v1 = a1 b1 = (f10, f11).( 0, b11) + # v2 = a2 b2 = (f20, f21).(b20, b21) + # + # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ + # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ + # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ + + var V0 {.noInit.}, V1 {.noInit.}, V2 {.noinit.}: doublePrec(Fpkdiv3) + var t0 {.noInit.}, t1 {.noInit.}: Fpkdiv3 + var f2x{.noInit.}, g2x {.noinit.}: doublePrec(Fpkdiv3) + + V0.prod2x(a.c0, b.c0) + V1.mul2x_sparse_by_0y(a.c1, b.c1) + V2.prod2x(a.c2, b.c2) + + # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ + t0.sum(a.c1, a.c2) + t1.c1.sum(b.c1.c1, b.c2.c1) # b₁ = ( 0, b11) + f2x.prod2x_disjoint(t0, b.c2.c0, t1.c1) # (a₁ + a₂).(b₁ + b₂) + f2x.diff2xMod(f2x, V1) + f2x.diff2xMod(f2x, V2) + f2x.prod2x(f2x, NonResidue) + f2x.sum2xMod(f2x, V0) + + # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + t0.sum(a.c0, a.c1) + t1.c1.sum(b.c0.c1, b.c1.c1) # b₁ = ( 0, b11) + g2x.prod2x_disjoint(t0, b.c0.c0, t1.c1) # (a₀ + a₁).(b₀ + b₁) + g2x.diff2xMod(g2x, V0) + g2x.diff2xMod(g2x, V1) + + # r₂ = (a₀ + a₂) and (b₀ + b₂) + t0.sum(a.c0, a.c2) + t1.sum(b.c0, b.c2) + + # Now we are aliasing free + + # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ + a.c0.redc2x(f2x) + + # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ + f2x.prod2x(V2, NonResidue) + g2x.sum2xMod(g2x, f2x) + a.c1.redc2x(g2x) + + # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ + f2x.prod2x(t0, t1) + f2x.diff2xMod(f2x, V0) + f2x.diff2xMod(f2x, V2) + f2x.sum2xMod(f2x, V1) + a.c2.redc2x(f2x) + +# Dispatch +# ------------------------------------------------------------ + +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) + when Fpk.C.getSexticTwist() == D_Twist: + f.mul_sparse_by_line_xyz000(line) + elif Fpk.C.getSexticTwist() == M_Twist: + f.mul_sparse_by_line_xy000z(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) + ## and store the result in f + ## f is overwritten + when Fpk.C.getSexticTwist() == D_Twist: + f.prod_xyz000_xyz000_into_abcdefghij00(line0, line1) + elif Fpk.C.getSexticTwist() == M_Twist: + f.prod_xy000z_xy000z_into_abcd00efghij(line0, line1) + 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) + ## 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".} diff --git a/constantine/pairing/lines_projective.nim b/constantine/pairing/lines_projective.nim deleted file mode 100644 index 52ed4e9..0000000 --- a/constantine/pairing/lines_projective.nim +++ /dev/null @@ -1,320 +0,0 @@ -# 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 - ../primitives, - ../config/curves, - ../arithmetic, - ../towers, - ../elliptic/[ - ec_shortweierstrass_affine, - ec_shortweierstrass_projective - ], - ./lines_common - -export lines_common - -# No exceptions allowed -{.push raises: [].} - -# ############################################################ -# -# Miller Loop's Line Evaluation -# with projective coordinates -# -# ############################################################ -# -# - Pairing Implementation Revisited -# Michael Scott, 2019 -# https://eprint.iacr.org/2019/077 -# -# - The Realm of the Pairings -# Diego F. Aranha and Paulo S. L. M. Barreto -# and Patrick Longa and Jefferson E. Ricardini, 2013 -# https://eprint.iacr.org/2013/722.pdf -# http://sac2013.irmacs.sfu.ca/slides/s1.pdf -# -# - Efficient Implementation of Bilinear Pairings on ARM Processors -# Gurleen Grewal, Reza Azarderakhsh, -# Patrick Longa, Shi Hu, and David Jao, 2012 -# https://eprint.iacr.org/2012/408.pdf - -# Line evaluation only -# ----------------------------------------------------------------------------- - -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 𝔽p4 / 𝔽p6 / 𝔽p12 - ## - ## 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 𝔽p4 / 𝔽p6 / 𝔽p12 - ## - ## 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₂) - -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) - 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.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) - 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 Field.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 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 - -func line_eval_fused_add[Field]( - line: var Line[Field], - T: var ECP_ShortW_Prj[Field, G2], - Q: ECP_ShortW_Aff[Field, G2]) = - ## Fused line evaluation and elliptic point addition - # Grewal et al, 2012 adapted to Scott 2019 line notation - var - A {.noInit.}: Field - B {.noInit.}: Field - C {.noInit.}: Field - D {.noInit.}: Field - E {.noInit.}: Field - F {.noInit.}: Field - G {.noInit.}: Field - H {.noInit.}: Field - I {.noInit.}: Field - - 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.diff(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 Field.C.getSexticTwist() == M_Twist: - lambda *= SexticNonResidue # A = ξ (X₁ - Z₁X₂) - -# Public proc -# ----------------------------------------------------------------------------- - -func line_double*[F1, F2]( - line: var Line[F2], - T: var ECP_ShortW_Prj[F2, G2], - P: ECP_ShortW_Aff[F1, G1]) = - ## Doubling step of the Miller loop - ## T in G2, P in G1 - ## - ## 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() - -func line_add*[F1, F2]( - line: var Line[F2], - T: var ECP_ShortW_Prj[F2, G2], - Q: ECP_ShortW_Aff[F2, G2], - P: ECP_ShortW_Aff[F1, G1]) = - ## Addition step of the Miller loop - ## T and Q in G2, P in G1 - ## - ## 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 diff --git a/constantine/pairing/miller_loops.nim b/constantine/pairing/miller_loops.nim index 46c24cc..0d2cef0 100644 --- a/constantine/pairing/miller_loops.nim +++ b/constantine/pairing/miller_loops.nim @@ -12,8 +12,7 @@ import ec_shortweierstrass_projective ], ../isogeny/frobenius, - ./lines_projective, - ./mul_fp6_by_lines, ./mul_fp12_by_lines + ./lines_eval # No exceptions allowed {.push raises: [].} @@ -48,15 +47,15 @@ template basicMillerLoop*[FT, F1, F2]( for i in countdown(u3.bits - 2, 1): square(f) line_double(line, T, P) - mul(f, line) + mul_by_line(f, line) let naf = bit(u3, i).int8 - bit(u, i).int8 # This can throw exception if naf == 1: line_add(line, T, Q, P) - mul(f, line) + mul_by_line(f, line) elif naf == -1: line_add(line, T, nQ, P) - mul(f, line) + mul_by_line(f, line) when pairing(C, ate_param_isNeg): # In GT, x^-1 == conjugate(x) @@ -82,12 +81,12 @@ func millerCorrectionBN*[FT, F1, F2]( V.frobenius_psi(Q) line.line_add(T, V, P) - f.mul(line) + f.mul_by_line(line) V.frobenius_psi(Q, 2) V.neg() line.line_add(T, V, P) - f.mul(line) + f.mul_by_line(line) # ############################################################ # # @@ -131,7 +130,6 @@ func miller_init_double_then_add*[FT, F1, F2]( ## f is overwritten ## T is overwritten by Q static: - doAssert f.c0 is Fp4 doAssert FT.C == F1.C doAssert FT.C == F2.C doAssert numDoublings >= 1 @@ -154,13 +152,13 @@ func miller_init_double_then_add*[FT, F1, F2]( # - The first line is squared (sparse * sparse) # - The second is (somewhat-sparse * sparse) when numDoublings >= 2: - f.prod_sparse_sparse(line, line) + f.prod_from_2_lines(line, line) line.line_double(T, P) - f.mul(line) + f.mul_by_line(line) for _ in 2 ..< numDoublings: f.square() line.line_double(T, P) - f.mul(line) + f.mul_by_line(line) # Addition step: 0b10...01 # ------------------------------------------------ @@ -172,10 +170,10 @@ func miller_init_double_then_add*[FT, F1, F2]( # f *= line <=> f = line for the first iteration var line2 {.noInit.}: Line[F2] line2.line_add(T, Q, P) - f.prod_sparse_sparse(line, line2) + f.prod_from_2_lines(line, line2) else: line.line_add(T, Q, P) - f.mul(line) + f.mul_by_line(line) {.pop.} # No OverflowError or IndexError allowed @@ -204,11 +202,11 @@ func miller_accum_double_then_add*[FT, F1, F2]( for _ in 0 ..< numDoublings: f.square() line.line_double(T, P) - f.mul(line) + f.mul_by_line(line) if add: line.line_add(T, Q, P) - f.mul(line) + f.mul_by_line(line) # Miller Loop - multi-pairing # ---------------------------------------------------------------------------- @@ -234,11 +232,11 @@ func double_jToN[N: static int, FT, F1, F2]( line0.line_double(Ts[i], Ps[i]) line1.line_double(Ts[i+1], Ps[i+1]) - f.mul_3way_sparse_sparse(line0, line1) + f.mul_by_2_lines(line0, line1) when (N and 1) == 1: # N >= 2 and N is odd, there is a leftover line0.line_double(Ts[N-1], Ps[N-1]) - f.mul(line0) + f.mul_by_line(line0) {.pop.} @@ -259,11 +257,11 @@ func add_jToN[N: static int, FT, F1, F2]( line0.line_add(Ts[i], Qs[i], Ps[i]) line1.line_add(Ts[i+1], Qs[i+1], Ps[i+1]) - f.mul_3way_sparse_sparse(line0, line1) + f.mul_by_2_lines(line0, line1) when (N and 1) == 1: # N >= 2 and N is odd, there is a leftover line0.line_add(Ts[N-1], Qs[N-1], Ps[N-1]) - f.mul(line0) + f.mul_by_line(line0) {.pop.} @@ -282,7 +280,6 @@ func miller_init_double_then_add*[N: static int, FT, F1, F2]( ## f is overwritten ## Ts are overwritten by Qs static: - doAssert f.c0 is Fp4 doAssert FT.C == F1.C doAssert FT.C == F2.C @@ -297,16 +294,16 @@ func miller_init_double_then_add*[N: static int, FT, F1, F2]( line0.line_double(Ts[0], Ps[0]) when N >= 2: line1.line_double(Ts[1], Ps[1]) - f.prod_sparse_sparse(line0, line1) + f.prod_from_2_lines(line0, line1) f.double_jToN(j=2, line0, line1, Ts, Ps) # Doubling steps: 0b10...00 # ------------------------------------------------ when numDoublings > 1: # Already did the MSB doubling when N == 1: # f = line0 - f.prod_sparse_sparse(line0, line0) # f.square() - line0.line_double(Ts[1], Ps[1]) - f.mul(line0) + f.prod_from_2_lines(line0, line0) # f.square() + line0.line_double(Ts[0], Ps[0]) + f.mul_by_line(line0) for _ in 2 ..< numDoublings: f.square() f.double_jtoN(j=0, line0, line1, Ts, Ps) @@ -320,7 +317,7 @@ func miller_init_double_then_add*[N: static int, FT, F1, F2]( when numDoublings == 1 and N == 1: # f = line0 line1.line_add(Ts[0], Qs[0], Ps[0]) - f.prod_sparse_sparse(line0, line1) + f.prod_from_2_lines(line0, line1) else: f.add_jToN(j=0,line0, line1, Ts, Qs, Ps) diff --git a/constantine/pairing/mul_fp12_by_lines.nim b/constantine/pairing/mul_fp12_by_lines.nim deleted file mode 100644 index 666e9ac..0000000 --- a/constantine/pairing/mul_fp12_by_lines.nim +++ /dev/null @@ -1,564 +0,0 @@ -# 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 - ../primitives, - ../config/curves, - ../arithmetic, - ../towers, - ./lines_projective - -# No exceptions allowed -{.push raises: [].} - -# ############################################################ -# -# Sparse Multiplication -# by lines for embedding degree 12 -# and sextic twist -# -# ############################################################ - -# - Pairing Implementation Revisited -# Michael Scott, 2019 -# https://eprint.iacr.org/2019/077 -# -# - Efficient Implementation of Bilinear Pairings on ARM Processors -# Gurleen Grewal, Reza Azarderakhsh, -# Patrick Longa, Shi Hu, and David Jao, 2012 -# https://eprint.iacr.org/2012/408.pdf -# -# - High-Speed Software Implementation of the Optimal Ate Pairing over Barreto-Naehrig Curves\ -# Jean-Luc Beuchat and Jorge Enrique González Díaz and Shigeo Mitsunari and Eiji Okamoto and Francisco Rodríguez-Henríquez and Tadanori Teruya, 2010\ -# https://eprint.iacr.org/2010/354 -# -# - Faster Pairing Computations on Curves with High-Degree Twists -# Craig Costello, Tanja Lange, and Michael Naehrig, 2009 -# https://eprint.iacr.org/2009/615.pdf - -# ############################################################ -# -# 𝔽p12 by line - 𝔽p12 quadratic over 𝔽p6 -# -# ############################################################ - -# D-Twist -# ------------------------------------------------------------ - -func mul_by_line_xy0*[C: static Curve]( - r: var Fp6[C], - a: Fp6[C], - b: Line[Fp2[C]]) = - ## Sparse multiplication of an 𝔽p6 - ## with coordinates (a₀, a₁, a₂) by a line (x, y, 0) - ## The z coordinates in the line will be ignored. - ## `r` and `a` must not alias - var - v0 {.noInit.}: Fp2[C] - v1 {.noInit.}: Fp2[C] - - v0.prod(a.c0, b.x) - v1.prod(a.c1, b.y) - r.c0.prod(a.c2, b.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.c1 *= r.c2 - r.c1 -= v0 - r.c1 -= v1 - - r.c2.prod(a.c2, b.x) - r.c2 += v1 - -func mul_sparse_by_line_xy00z0*[C: static Curve]( - f: var Fp12[C], l: Line[Fp2[C]]) = - ## Sparse multiplication of an 𝔽p12 element - ## by a sparse 𝔽p12 element coming from an D-Twist line function. - ## The sparse element is represented by a packed Line type - ## with coordinate (x,y,z) matching 𝔽p12 coordinates xy00z0 (TODO: verify this) - - static: - doAssert C.getSexticTwist() == D_Twist - doAssert f.c0 is Fp6, "This assumes 𝔽p12 as a quadratic extension of 𝔽p6" - - var - v0 {.noInit.}: Fp6[C] - v1 {.noInit.}: Fp6[C] - v2 {.noInit.}: Line[Fp2[C]] - v3 {.noInit.}: Fp6[C] - - v0.mul_by_line_xy0(f.c0, l) - v1.mul_sparse_by_0y0(f.c1, l.z) - - v2.x = l.x - v2.y.sum(l.y, l.z) - f.c1 += f.c0 - v3.mul_by_line_xy0(f.c1, v2) - v3 -= v0 - f.c1.diff(v3, v1) - - v3.c0.prod(v1.c2, SexticNonResidue) - v3.c0 += v0.c0 - v3.c1.sum(v0.c1, v1.c0) - v3.c2.sum(v0.c2, v1.c1) - f.c0 = v3 - -# ############################################################ -# -# 𝔽p12 by line - 𝔽p12 cubic over 𝔽p4 -# -# ############################################################ - -# D-Twist -# ------------------------------------------------------------ - -func mul_sparse_by_line_xyz000*[C: static Curve]( - f: var Fp12[C], l: Line[Fp2[C]]) = - ## Sparse multiplication of an 𝔽p12 element - ## by a sparse 𝔽p12 element coming from an D-Twist line function. - ## The sparse element is represented by a packed Line type - ## with coordinates (x,y,z) matching 𝔽p12 coordinates xyz000 - - static: - doAssert C.getSexticTwist() == D_Twist - doAssert f.c0 is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4" - - # In the following equations (taken from cubic extension implementation) - # a = f - # b0 = (x, y) - # b1 = (z, 0) - # b2 = (0, 0) - # - # v0 = a0 b0 = (f00, f01).(x, y) - # v1 = a1 b1 = (f10, f11).(z, 0) - # v2 = a2 b2 = (f20, f21).(0, 0) - # - # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 - # = ξ (a1 b1 + a2 b1 - v1) + v0 - # = ξ a2 b1 + v0 - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 - # = (a0 + a1) * (b0 + b1) - v0 - v1 - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 - # = a0 b0 + a2 b0 - v0 + v1 - # = a2 b0 + v1 - - when false: - var b0 {.noInit.}, v0{.noInit.}, v1{.noInit.}, t{.noInit.}: Fp4[C] - - b0.c0 = l.x - b0.c1 = l.y - - v0.prod(f.c0, b0) - v1.mul_sparse_by_x0(f.c1, l.z) - - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 - f.c1 += f.c0 # r1 = a0 + a1 - t = b0 - t.c0 += l.z # t = b0 + b1 - f.c1 *= t # r2 = (a0 + a1)(b0 + b1) - f.c1 -= v0 - 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 *= SexticNonResidue - f.c0 += v0 - - # r2 = a2 b0 + v1 - f.c2 *= b0 - f.c2 += v1 - - else: # Lazy reduction - var V0{.noInit.}, V1{.noInit.}, f2x{.noInit.}: doublePrec(Fp4[C]) - var t{.noInit.}: Fp2[C] - - V0.prod2x_disjoint(f.c0, l.x, l.y) - V1.mul2x_sparse_by_x0(f.c1, l.z) - - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 - when false: # TODO: what's the condition? - f.c1.sumUnr(f.c1, f.c0) - t.sumUnr(l.x, l.z) # b0 is (x, y) - else: - 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) - 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.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.sum2xMod(f2x, V1) - f.c2.redc2x(f2x) - -func prod_xyz000_xyz000_into_abcdefghij00*[C: static Curve](f: var Fp12[C], l0, l1: Line[Fp2[C]]) = - ## 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) - # a2 = ( 0, 0) - # b0 = (x1, y1) - # b1 = (z1, 0) - # b2 = ( 0, 0) - # - # v0 = a0 b0 = (x0, y0).(x1, y1) - # v1 = a1 b1 = (z0, 0).(z1, 0) - # v2 = a2 b2 = ( 0, 0).( 0, 0) - # - # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 - # = ξ (a1 b1 + a2 b1 - v1) + v0 - # = v0 - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 - # = (a0 + a1) * (b0 + b1) - v0 - v1 - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 - # = a0 b0 - v0 + v1 - # = v1 - - static: - doAssert C.getSexticTwist() == D_Twist - doAssert f.c0 is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4" - - var V0{.noInit.}, f2x{.noInit.}: doublePrec(Fp4[C]) - var V1{.noInit.}: doublePrec(Fp2[C]) - - 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) - - # 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) - f2x.diff2xMod(f2x, V0) - f2x.c0.diff2xMod(f2x.c0, V1) - f.c1.redc2x(f2x) - - # r0 = v0 - f.c0.redc2x(V0) - - # r2 = v1 - f.c2.c0.redc2x(V1) - f.c2.c1.setZero() - -func mul_sparse_by_abcdefghij00*[C: static Curve]( - a: var Fp12[C], b: Fp12[C]) = - ## Sparse multiplication of an 𝔽p12 element - ## by a sparse 𝔽p12 element abcdefghij00 - ## with each representing 𝔽p2 coordinate - - static: - doAssert C.getSexticTwist() == D_Twist - doAssert a.c0 is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4" - - # In the following equations (taken from cubic extension implementation) - # b0 = (b00, b01) - # b1 = (b10, b11) - # b2 = (b20, 0) - # - # v0 = a0 b0 = (f00, f01).(b00, b01) - # v1 = a1 b1 = (f10, f11).(b10, b11) - # v2 = a2 b2 = (f20, f21).(b20, 0) - # - # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ - # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ - # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ - - var V0 {.noInit.}, V1 {.noInit.}, V2 {.noinit.}: doublePrec(Fp4[C]) - var t0 {.noInit.}, t1 {.noInit.}: Fp4[C] - var f2x{.noInit.}, g2x {.noinit.}: doublePrec(Fp4[C]) - - V0.prod2x(a.c0, b.c0) - V1.prod2x(a.c1, b.c1) - V2.mul2x_sparse_by_x0(a.c2, b.c2) - - # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ - t0.sum(a.c1, a.c2) - t1.c0.sum(b.c1.c0, b.c2.c0) # b₂ = (b20, 0) - f2x.prod2x_disjoint(t0, t1.c0, b.c1.c1) # (a₁ + a₂).(b₁ + b₂) - f2x.diff2xMod(f2x, V1) - f2x.diff2xMod(f2x, V2) - f2x.prod2x(f2x, NonResidue) - f2x.sum2xMod(f2x, V0) - - # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ - t0.sum(a.c0, a.c1) - t1.sum(b.c0, b.c1) - g2x.prod2x(t0, t1) - g2x.diff2xMod(g2x, V0) - g2x.diff2xMod(g2x, V1) - - # r₂ = (a₀ + a₂) and (b₀ + b₂) - t0.sum(a.c0, a.c2) - t1.c0.sum(b.c0.c0, b.c2.c0) # b₂ = (b20, 0) - - # Now we are aliasing free - - # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ - a.c0.redc2x(f2x) - - # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ - f2x.prod2x(V2, NonResidue) - g2x.sum2xMod(g2x, f2x) - a.c1.redc2x(g2x) - - # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ - f2x.prod2x_disjoint(t0, t1.c0, b.c0.c1) - f2x.diff2xMod(f2x, V0) - f2x.diff2xMod(f2x, V2) - f2x.sum2xMod(f2x, V1) - a.c2.redc2x(f2x) - -# M-Twist -# ------------------------------------------------------------ - -func mul_sparse_by_line_xy000z*[C: static Curve]( - f: var Fp12[C], l: Line[Fp2[C]]) = - - static: - doAssert C.getSexticTwist() == M_Twist - doAssert f.c0 is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4" - - # In the following equations (taken from cubic extension implementation) - # a = f - # b0 = (x, y) - # b1 = (0, 0) - # b2 = (0, z) - # - # v0 = a0 b0 = (f00, f01).(x, y) - # v1 = a1 b1 = (f10, f11).(0, 0) - # v2 = a2 b2 = (f20, f21).(0, z) - # - # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 - # = ξ (a1 b2 + a2 b2 - v2) + v0 - # = ξ a1 b2 + v0 - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 - # = a0 b0 + a1 b0 - v0 + ξ v2 - # = a1 b0 + ξ v2 - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 - # = (a0 + a2) * (b0 + b2) - v0 - v2 - - when false: - var b0 {.noInit.}, v0{.noInit.}, v2{.noInit.}, t{.noInit.}: Fp4[C] - - b0.c0 = l.x - b0.c1 = l.y - - v0.prod(f.c0, b0) - v2.mul_sparse_by_0y(f.c2, l.z) - - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 - f.c2 += f.c0 # r2 = a0 + a2 - t = b0 - t.c1 += l.z # 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 *= SexticNonResidue - f.c0 += v0 - - # r1 = a1 b0 + ξ v2 - f.c1 *= b0 - v2 *= SexticNonResidue - f.c1 += v2 - - else: # Lazy reduction - var V0{.noInit.}, V2{.noInit.}, f2x{.noInit.}: doublePrec(Fp4[C]) - var t{.noInit.}: Fp2[C] - - V0.prod2x_disjoint(f.c0, l.x, l.y) - V2.mul2x_sparse_by_0y(f.c2, l.z) - - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 - when false: # TODO: what's the condition - f.c2.sumUnr(f.c2, f.c0) - t.sumUnr(l.y, l.z) # b0 is (x, y) - else: - 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) - 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.prod2x(f2x, SexticNonResidue) - f2x.sum2xMod(f2x, V0) - f.c0.redc2x(f2x) - - # r1 = a1 b0 + ξ v2 - f2x.prod2x_disjoint(f.c1, l.x, l.y) - V2.prod2x(V2, SexticNonResidue) - f2x.sum2xMod(f2x, V2) - f.c1.redc2x(f2x) - -func prod_xy000z_xy000z_into_abcd00efghij*[C: static Curve](f: var Fp12[C], l0, l1: Line[Fp2[C]]) = - ## Multiply 2 lines together - ## The result is sparse in f.c1.c0 - # In the following equations (taken from cubic extension implementation) - # a0 = (x0, y0) - # a1 = ( 0, 0) - # a2 = ( 0, z0) - # b0 = (x1, y1) - # b1 = ( 0, 0) - # b2 = ( 0, z1) - # - # v0 = a0 b0 = (x0, y0).(x1, y1) - # v1 = a1 b1 = ( 0, 0).( 0, 0) - # v2 = a2 b2 = ( 0, z0).( 0, z1) - # - # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 - # = ξ (a1 b2 + a2 b2 - v2) + v0 - # = v0 - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 - # = a0 b0 + a1 b0 - v0 + ξ v2 - # = ξ v2 - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 - # = (a0 + a2) * (b0 + b2) - v0 - v2 - - static: - doAssert C.getSexticTwist() == M_Twist - doAssert f.c0 is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4" - - var V0{.noInit.}, f2x{.noInit.}: doublePrec(Fp4[C]) - var V2{.noInit.}: doublePrec(Fp2[C]) - - 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) - - # 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) - 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) - - # r1 = ξ v2 - f.c1.c1.redc2x(V2) - f.c1.c0.setZero() - - # r0 = v0 - f.c0.redc2x(V0) - -func mul_sparse_by_abcd00efghij*[C: static Curve]( - a: var Fp12[C], b: Fp12[C]) = - ## Sparse multiplication of an 𝔽p12 element - ## by a sparse 𝔽p12 element abcd00efghij - ## with each representing 𝔽p2 coordinate - - static: - doAssert C.getSexticTwist() == M_Twist - doAssert a.c0 is Fp4, "This assumes 𝔽p12 as a cubic extension of 𝔽p4" - - # In the following equations (taken from cubic extension implementation) - # b0 = (b00, b01) - # b1 = ( 0, b11) - # b2 = (b20, b21) - # - # v0 = a0 b0 = (f00, f01).(b00, b01) - # v1 = a1 b1 = (f10, f11).( 0, b11) - # v2 = a2 b2 = (f20, f21).(b20, b21) - # - # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ - # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ - # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ - - var V0 {.noInit.}, V1 {.noInit.}, V2 {.noinit.}: doublePrec(Fp4[C]) - var t0 {.noInit.}, t1 {.noInit.}: Fp4[C] - var f2x{.noInit.}, g2x {.noinit.}: doublePrec(Fp4[C]) - - V0.prod2x(a.c0, b.c0) - V1.mul2x_sparse_by_0y(a.c1, b.c1) - V2.prod2x(a.c2, b.c2) - - # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ - t0.sum(a.c1, a.c2) - t1.c1.sum(b.c1.c1, b.c2.c1) # b₁ = ( 0, b11) - f2x.prod2x_disjoint(t0, b.c2.c0, t1.c1) # (a₁ + a₂).(b₁ + b₂) - f2x.diff2xMod(f2x, V1) - f2x.diff2xMod(f2x, V2) - f2x.prod2x(f2x, NonResidue) - f2x.sum2xMod(f2x, V0) - - # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ - t0.sum(a.c0, a.c1) - t1.c1.sum(b.c0.c1, b.c1.c1) # b₁ = ( 0, b11) - g2x.prod2x_disjoint(t0, b.c0.c0, t1.c1) # (a₀ + a₁).(b₀ + b₁) - g2x.diff2xMod(g2x, V0) - g2x.diff2xMod(g2x, V1) - - # r₂ = (a₀ + a₂) and (b₀ + b₂) - t0.sum(a.c0, a.c2) - t1.sum(b.c0, b.c2) - - # Now we are aliasing free - - # r₀ = ξ ((a₁ + a₂)(b₁ + b₂) - v₁ - v₂) + v₀ - a.c0.redc2x(f2x) - - # r₁ = (a₀ + a₁) * (b₀ + b₁) - v₀ - v₁ + β v₂ - f2x.prod2x(V2, NonResidue) - g2x.sum2xMod(g2x, f2x) - a.c1.redc2x(g2x) - - # r₂ = (a₀ + a₂) * (b₀ + b₂) - v₀ - v₂ + v₁ - f2x.prod2x(t0, t1) - f2x.diff2xMod(f2x, V0) - f2x.diff2xMod(f2x, V2) - f2x.sum2xMod(f2x, V1) - a.c2.redc2x(f2x) - -# Dispatch -# ------------------------------------------------------------ - -func mul*[C](f: var Fp12[C], line: Line[Fp2[C]]) {.inline.} = - ## Multiply an element of Fp12 by a sparse line function (xyz000 or xy000z) - when C.getSexticTwist() == D_Twist: - f.mul_sparse_by_line_xyz000(line) - elif C.getSexticTwist() == M_Twist: - f.mul_sparse_by_line_xy000z(line) - else: - {.error: "A line function assumes that the curve has a twist".} - -func prod_sparse_sparse*[C](f: var Fp12[C], line0, line1: Line[Fp2[C]]) {.inline.} = - ## Multiply 2 lines function (xyz000 or xy000z) - ## and store the result in f - ## f is overwritten - when C.getSexticTwist() == D_Twist: - f.prod_xyz000_xyz000_into_abcdefghij00(line0, line1) - elif C.getSexticTwist() == M_Twist: - f.prod_xy000z_xy000z_into_abcd00efghij(line0, line1) - else: - {.error: "A line function assumes that the curve has a twist".} - -func mul_3way_sparse_sparse*[C](f: var Fp12[C], line0, line1: Line[Fp2[C]]) {.inline.} = - ## Multiply f*line0*line1 with lines (xyz000 or xy000z) - ## f is updated with the result - var t{.noInit.}: typeof(f) - when C.getSexticTwist() == D_Twist: - t.prod_xyz000_xyz000_into_abcdefghij00(line0, line1) - f.mul_sparse_by_abcdefghij00(t) - elif 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".} diff --git a/constantine/pairing/mul_fp6_by_lines.nim b/constantine/pairing/mul_fp6_by_lines.nim deleted file mode 100644 index 6474a6b..0000000 --- a/constantine/pairing/mul_fp6_by_lines.nim +++ /dev/null @@ -1,143 +0,0 @@ -# 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 - ../primitives, - ../config/curves, - ../arithmetic, - ../towers, - ./lines_projective - -# No exceptions allowed -{.push raises: [].} - -# ############################################################ -# -# Sparse Multiplication -# by lines for embedding degree 6 -# and sextic twist -# -# ############################################################ - -# 𝔽p6 by line - Sparse functions -# ---------------------------------------------------------------- - -func mul_sparse_by_line_xyz000*[C: static Curve]( - f: var Fp6[C], l: Line[Fp[C]]) = - ## Sparse multiplication of an 𝔽p6 element - ## by a sparse 𝔽p6 element coming from an D-Twist line function. - ## The sparse element is represented by a packed Line type - ## with coordinates (x,y,z) matching 𝔽p6 coordinates xyz000 - - static: - doAssert C.getSexticTwist() == D_Twist - doAssert f.c0.typeof is Fp2, "This assumes 𝔽p6 as a cubic extension of 𝔽p6" - - # In the following equations (taken from cubic extension implementation) - # a = f - # b0 = (x, y) - # b1 = (z, 0) - # b2 = (0, 0) - # - # v0 = a0 b0 = (f00, f01).(x, y) - # v1 = a1 b1 = (f10, f11).(z, 0) - # v2 = a2 b2 = (f20, f21).(0, 0) - # - # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 - # = ξ (a1 b1 + a2 b1 - v1) + v0 - # = ξ a2 b1 + v0 - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 - # = (a0 + a1) * (b0 + b1) - v0 - v1 - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 - # = a0 b0 + a2 b0 - v0 + v1 - # = a2 b0 + v1 - - var b0 {.noInit.}, v0{.noInit.}, v1{.noInit.}, t{.noInit.}: Fp2[C] - - b0.c0 = l.x - b0.c1 = l.y - - v0.prod(f.c0, b0) - v1.mul_sparse_by_y0(f.c1, l.z) - - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 - f.c1 += f.c0 # r1 = a0 + a1 - t = b0 - t.c0 += l.z # t = b0 + b1 - f.c1 *= t # r2 = (a0 + a1)(b0 + b1) - f.c1 -= v0 - f.c1 -= v1 # r2 = (a0 + a1)(b0 + b1) - v0 - v1 - - # r0 = ξ a2 b1 + v0 - f.c0.mul_sparse_by_y0(f.c2, l.z) - f.c0 *= NonResidue - f.c0 += v0 - - # r2 = a2 b0 + v1 - f.c2 *= b0 - f.c2 += v1 - -func mul_sparse_by_line_xy000z*[C: static Curve]( - f: var Fp6[C], l: Line[Fp[C]]) = - - static: - doAssert C.getSexticTwist() == M_Twist - doAssert f.c0.typeof is Fp2, "This assumes 𝔽p6 as a cubic extension of 𝔽p2" - - # In the following equations (taken from cubic extension implementation) - # a = f - # b0 = (x, y) - # b1 = (0, 0) - # b2 = (0, z) - # - # v0 = a0 b0 = (f00, f01).(x, y) - # v1 = a1 b1 = (f10, f11).(0, 0) - # v2 = a2 b2 = (f20, f21).(0, z) - # - # r0 = ξ ((a1 + a2) * (b1 + b2) - v1 - v2) + v0 - # = ξ (a1 b2 + a2 b2 - v2) + v0 - # = ξ a1 b2 + v0 - # r1 = (a0 + a1) * (b0 + b1) - v0 - v1 + ξ v2 - # = a0 b0 + a1 b0 - v0 + ξ v2 - # = a1 b0 + ξ v2 - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 + v1 - # = (a0 + a2) * (b0 + b2) - v0 - v2 - - var b0 {.noInit.}, v0{.noInit.}, v2{.noInit.}, t{.noInit.}: Fp2[C] - - b0.c0 = l.x - b0.c1 = l.y - - v0.prod(f.c0, b0) - v2.mul_sparse_by_0y(f.c2, l.z) - - # r2 = (a0 + a2) * (b0 + b2) - v0 - v2 - f.c2 += f.c0 # r2 = a0 + a2 - t = b0 - t.c1 += l.z # 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 *= NonResidue - f.c0 += v0 - - # r1 = a1 b0 + ξ v2 - f.c1 *= b0 - v2 *= NonResidue - f.c1 += v2 - -func mul*[C](f: var Fp6[C], line: Line[Fp[C]]) {.inline.} = - when C.getSexticTwist() == D_Twist: - f.mul_sparse_by_line_xyz000(line) - elif C.getSexticTwist() == M_Twist: - f.mul_sparse_by_line_xy000z(line) - else: - {.error: "A line function assumes that the curve has a twist".} diff --git a/constantine/pairing/pairing_bls12.nim b/constantine/pairing/pairing_bls12.nim index 07cc051..6759641 100644 --- a/constantine/pairing/pairing_bls12.nim +++ b/constantine/pairing/pairing_bls12.nim @@ -15,8 +15,8 @@ import ], ../isogeny/frobenius, ../curves/zoo_pairings, - ./cyclotomic_fp12, - ./lines_common, + ./cyclotomic_subgroup, + ./lines_eval, ./miller_loops export zoo_pairings # generic sandwich https://github.com/nim-lang/Nim/issues/11225 @@ -119,17 +119,17 @@ func finalExpHard_BLS12*[C](f: var Fp12[C]) {.meter.} = # (x−1)² when C.pairing(ate_param).isEven.bool: - v0.pow_xdiv2(v2) # v0 = (f²)^(x/2) = f^x + v0.cycl_exp_by_curve_param_div2(v2) # v0 = (f²)^(x/2) = f^x else: - v0.pow_x(f) + v0.cycl_exp_by_curve_param(f) v1.cyclotomic_inv(f) # v1 = f^-1 v0 *= v1 # v0 = f^(x-1) - v1.pow_x(v0) # v1 = (f^(x-1))^x + v1.cycl_exp_by_curve_param(v0) # v1 = (f^(x-1))^x v0.cyclotomic_inv() # v0 = (f^(x-1))^-1 v0 *= v1 # v0 = (f^(x-1))^(x-1) = f^((x-1)*(x-1)) = f^((x-1)²) # (x+p) - v1.pow_x(v0) # v1 = f^((x-1)².x) + v1.cycl_exp_by_curve_param(v0) # v1 = f^((x-1)².x) v0.frobenius_map(v0) # v0 = f^((x-1)².p) v0 *= v1 # v0 = f^((x-1)².(x+p)) @@ -137,8 +137,8 @@ func finalExpHard_BLS12*[C](f: var Fp12[C]) {.meter.} = f *= v2 # f = f³ # (x²+p²−1) - v2.pow_x(v0, invert = false) - v1.pow_x(v2, invert = false) # v1 = f^((x-1)².(x+p).x²) + v2.cycl_exp_by_curve_param(v0, invert = false) + v1.cycl_exp_by_curve_param(v2, invert = false) # v1 = f^((x-1)².(x+p).x²) v2.frobenius_map(v0, 2) # v2 = f^((x-1)².(x+p).p²) v0.cyclotomic_inv() # v0 = f^((x-1)².(x+p).-1) v0 *= v1 # v0 = f^((x-1)².(x+p).(x²-1)) diff --git a/constantine/pairing/pairing_bn.nim b/constantine/pairing/pairing_bn.nim index 623fe52..4e65779 100644 --- a/constantine/pairing/pairing_bn.nim +++ b/constantine/pairing/pairing_bn.nim @@ -15,9 +15,8 @@ import ], ../isogeny/frobenius, ../curves/zoo_pairings, - ./lines_projective, - ./mul_fp12_by_lines, - ./cyclotomic_fp12, + ./lines_eval, + ./cyclotomic_subgroup, ./miller_loops export zoo_pairings # generic sandwich https://github.com/nim-lang/Nim/issues/11225 @@ -114,11 +113,11 @@ func finalExpHard_BN*[C: static Curve](f: var Fp12[C]) = # that requires another addition chain var t0 {.noInit.}, t1 {.noinit.}, t2 {.noinit.}, t3 {.noinit.}, t4 {.noinit.}: Fp12[C] - t0.pow_u(f, invert = false) # t0 = f^|u| + t0.cycl_exp_by_curve_param(f, invert = false) # t0 = f^|u| t0.cyclotomic_square() # t0 = f^2|u| t1.cyclotomic_square(t0) # t1 = f^4|u| t1 *= t0 # t1 = f^6|u| - t2.pow_u(t1, invert = false) # t2 = f^6u² + t2.cycl_exp_by_curve_param(t1, invert = false) # t2 = f^6u² if C.pairing(ate_param_is_Neg): t3.cyclotomic_inv(t1) # t3 = f^6u @@ -126,7 +125,7 @@ func finalExpHard_BN*[C: static Curve](f: var Fp12[C]) = t3 = t1 # t3 = f^6u t1.prod(t2, t3) # t1 = f^6u.f^6u² t3.cyclotomic_square(t2) # t3 = f^12u² - t4.pow_u(t3) # t4 = f^12u³ + t4.cycl_exp_by_curve_param(t3) # t4 = f^12u³ t4 *= t1 # t4 = f^(6u + 6u² + 12u³) = f^λ₂ if not C.pairing(ate_param_is_Neg): diff --git a/constantine/pairing/pairing_bw6_761.nim b/constantine/pairing/pairing_bw6_761.nim index 3d128d8..47a6c1d 100644 --- a/constantine/pairing/pairing_bw6_761.nim +++ b/constantine/pairing/pairing_bw6_761.nim @@ -16,8 +16,7 @@ import ], ../isogeny/frobenius, ../curves/zoo_pairings, - ./lines_projective, - ./mul_fp6_by_lines, + ./lines_eval, ./miller_loops export zoo_pairings # generic sandwich https://github.com/nim-lang/Nim/issues/11225 @@ -100,12 +99,12 @@ func millerLoopBW6_761_opt_to_debug[C]( for i in countdown(u.bits - 2, 1): square(f) line_double(line, T, P) - mul(f, line) + mul_by_line(f, line) let bit = u.bit(i).int8 if bit == 1: line_add(line, T, Q, P) - mul(f, line) + mul_by_line(f, line) # Fixup # ------------------------------ @@ -120,7 +119,7 @@ func millerLoopBW6_761_opt_to_debug[C]( # Drop the vertical line line.line_add(T, Q, P) # TODO: eval without updating T muplusone = mu - muplusone.mul(line) + muplusone.mul_by_line(line) # 2nd part: f_{u²-u-1,Q}(P) # ------------------------------ @@ -133,16 +132,16 @@ func millerLoopBW6_761_opt_to_debug[C]( for i in countdown(u3.bits - 2, 1): square(f) line_double(line, T, P) - mul(f, line) + mul_by_line(f, line) let naf = bit(u3, i).int8 - bit(u, i).int8 # This can throw exception if naf == 1: line_add(line, T, Qu, P) - mul(f, line) + mul_by_line(f, line) f *= mu elif naf == -1: line_add(line, T, nQu, P) - mul(f, line) + mul_by_line(f, line) f *= minvu # Final diff --git a/constantine/protocols/ethereum_evm_precompiles.nim b/constantine/protocols/ethereum_evm_precompiles.nim index 27402b7..62d3816 100644 --- a/constantine/protocols/ethereum_evm_precompiles.nim +++ b/constantine/protocols/ethereum_evm_precompiles.nim @@ -11,7 +11,7 @@ import ../arithmetic, ../towers, ../arithmetic/limbs_montgomery, ../ec_shortweierstrass, - ../pairing/[pairing_bn, miller_loops, cyclotomic_fp12], + ../pairing/[pairing_bn, miller_loops, cyclotomic_subgroup], ../curves/zoo_subgroups, ../io/[io_bigints, io_fields] diff --git a/constantine/tower_field_extensions/assembly/fp2_asm_x86_adx_bmi2.nim b/constantine/tower_field_extensions/assembly/fp2_asm_x86_adx_bmi2.nim index d8fea33..9965b71 100644 --- a/constantine/tower_field_extensions/assembly/fp2_asm_x86_adx_bmi2.nim +++ b/constantine/tower_field_extensions/assembly/fp2_asm_x86_adx_bmi2.nim @@ -43,10 +43,6 @@ func has1extraBit(F: type Fp): bool = ## We construct extensions only on Fp (and not Fr) getSpareBits(F) >= 1 -func has2extraBits(F: type Fp): bool = - ## We construct extensions only on Fp (and not Fr) - getSpareBits(F) >= 2 - # 𝔽p2 squaring # ------------------------------------------------------------ @@ -67,9 +63,9 @@ func sqrx2x_complex_asm_adx*( t0.double(a.c1) t1.sum(a.c0, a.c1) - r.c1.mul_asm_adx_inline(t0, a.c0) + r.c1.limbs2x.mul_asm_adx_inline(t0.mres.limbs, a.c0.mres.limbs) t0.diff(a.c0, a.c1) - r.c0.mul_asm_adx_inline(t0, t1) + r.c0.limbs2x.mul_asm_adx_inline(t0.mres.limbs, t1.mres.limbs) func sqrx_complex_sparebit_asm_adx*( r: var array[2, Fp], diff --git a/constantine/tower_field_extensions/extension_fields.nim b/constantine/tower_field_extensions/extension_fields.nim index 6635f84..a68502f 100644 --- a/constantine/tower_field_extensions/extension_fields.nim +++ b/constantine/tower_field_extensions/extension_fields.nim @@ -293,11 +293,14 @@ template doublePrec*(T: type ExtensionField): type = when T is QuadraticExt: when T.F is QuadraticExt: # Fp4Dbl QuadraticExt2x[QuadraticExt2x[doublePrec(T.F.F)]] + elif T.F is CubicExt: + when T.F.F is QuadraticExt: # Fp12 over Fp6 over Fp2 + QuadraticExt2x[CubicExt2x[QuadraticExt2x[doublePrec(T.F.F.F)]]] elif T.F is Fp: # Fp2Dbl QuadraticExt2x[doublePrec(T.F)] elif T is CubicExt: when T.F is QuadraticExt: # - when T.F.F is QuadraticExt: # Fp12 + when T.F.F is QuadraticExt: # Fp12 over Fp4 over Fp2 CubicExt2x[QuadraticExt2x[QuadraticExt2x[doublePrec(T.F.F.F)]]] elif T.F.F is Fp: # Fp6 CubicExt2x[QuadraticExt2x[doublePrec(T.F.F)]] @@ -684,6 +687,47 @@ func prod2x*( {.pop.} # inline +# ############################################################ +# # +# Cost functions # +# # +# ############################################################ + +func prefer_3sqr_over_2mul(F: type ExtensionField): bool {.compileTime.} = + ## Returns true + ## if time(3sqr) < time(2mul) in the extension fields + + let a = default(F) + # No shortcut in the VM + when a.c0 is ExtensionField: + when a.c0.c0 is ExtensionField: + return true + else: return false + else: return false + +func has_large_NR_norm(C: static Curve): bool = + ## Returns true if the non-residue of the extension fields + ## has a large norm + + const j = C.getNonResidueFp() + const u = C.getNonResidueFp2()[0] + const v = C.getNonResidueFp2()[1] + + const norm2 = u*u + (j*v)*(j*v) + + # Compute integer square root + var norm = 0 + while (norm+1) * (norm+1) <= norm2: + norm += 1 + + return norm > 5 + +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]) + return a.mres.limbs.len > 6 + # ############################################################ # # # Quadratic extensions # @@ -863,17 +907,10 @@ func square_generic(r: var QuadraticExt, a: QuadraticExt) = # 2 c0 c1 <=> (a0 + a1)² - a0² - a1² # # This gives us 3 Sqr and 1 Mul-non-residue - const costlyMul = block: - # No shortcutting in the VM :/ - when a.c0 is ExtensionField: - when a.c0.c0 is ExtensionField: - true - else: - false - else: - false - when QuadraticExt.C == BN254_Snarks or costlyMul: + when QuadraticExt.prefer_3sqr_over_2mul() or + # Other path multiplies twice by non-residue + QuadraticExt.C.has_large_NR_norm(): var v0 {.noInit.}, v1 {.noInit.}: typeof(r.c0) v0.square(a.c0) v1.square(a.c1) @@ -923,9 +960,8 @@ func square2x_disjoint*[Fdbl, F]( var V0 {.noInit.}, V1 {.noInit.}: typeof(r.c0) # Double-precision var t {.noInit.}: F # Single-width - # TODO: which is the best formulation? 3 squarings or 2 Mul? - # It seems like the higher the tower the better squarings are - # So for Fp12 = 2xFp6, prefer squarings. + # square2x is faster than prod2x even on Fp + # hence we always use 3 squarings over 2 products V0.square2x(a0) V1.square2x(a1) t.sum(a0, a1) @@ -1234,14 +1270,20 @@ func inv2xImpl(r: var QuadraticExt, a: QuadraticExt) = func square2x*(r: var QuadraticExt2x, a: QuadraticExt) = when a.fromComplexExtension(): - r.square2x_complex(a) + when UseASM_X86_64 and not QuadraticExt.C.has_large_field_elem(): + if ({.noSideEffect.}: hasAdx()): + r.coords.sqrx2x_complex_asm_adx(a.coords) + else: + r.square2x_complex(a) + else: + r.square2x_complex(a) else: r.square2x_disjoint(a.c0, a.c1) func square*(r: var QuadraticExt, a: QuadraticExt) = when r.fromComplexExtension(): when true: - when UseASM_X86_64 and a.c0.mres.limbs.len <= 6 and r.typeof.has1extraBit(): + when UseASM_X86_64 and not QuadraticExt.C.has_large_field_elem() and r.typeof.has1extraBit(): if ({.noSideEffect.}: hasAdx()): r.coords.sqrx_complex_sparebit_asm_adx(a.coords) else: @@ -1250,7 +1292,7 @@ func square*(r: var QuadraticExt, a: QuadraticExt) = r.square_complex(a) else: # slower var d {.noInit.}: doublePrec(typeof(r)) - d.square2x_complex(a) + d.square2x(a) r.c0.redc2x(d.c0) r.c1.redc2x(d.c1) else: @@ -1279,7 +1321,7 @@ func prod*(r: var QuadraticExt, a, b: QuadraticExt) = when false: r.prod_complex(a, b) else: # faster - when UseASM_X86_64 and a.c0.mres.limbs.len <= 6: + when UseASM_X86_64 and not QuadraticExt.C.has_large_field_elem(): if ({.noSideEffect.}: hasAdx()): r.coords.mul_fp2_complex_asm_adx(a.coords, b.coords) else: @@ -1293,7 +1335,7 @@ func prod*(r: var QuadraticExt, a, b: QuadraticExt) = r.c0.redc2x(d.c0) r.c1.redc2x(d.c1) else: - when r.typeof.F.C == BW6_761 or typeof(r.c0) is Fp: + when r.typeof.F.C.has_large_field_elem(): # BW6-761 requires too many registers for Dbl width path r.prod_generic(a, b) else: @@ -1316,7 +1358,7 @@ func prod2x_disjoint*[Fdbl, F]( func prod2x*(r: var QuadraticExt2x, a, b: QuadraticExt) = ## Double-precision multiplication r <- a*b when a.fromComplexExtension(): - when UseASM_X86_64 and a.c0.mres.limbs.len <= 6: + when UseASM_X86_64 and not QuadraticExt.C.has_large_field_elem(): if ({.noSideEffect.}: hasAdx()): r.coords.mul2x_fp2_complex_asm_adx(a.coords, b.coords) else: @@ -1414,7 +1456,7 @@ func inv*(a: var QuadraticExt) = # Cubic extensions can use specific squaring procedures # beyond Schoolbook and Karatsuba: # - Chung-Hasan (3 different algorithms) -# - Toom-Cook-3x +# - Toom-Cook-3x (only in Miller Loops) # # Chung-Hasan papers # http://cacr.uwaterloo.ca/techreports/2006/cacr2006-24.pdf @@ -1434,11 +1476,11 @@ func inv*(a: var QuadraticExt) = # |-------------|----------|-------------------| # | Schoolbook | 3M + 3S | 6A + 2B | # | Karatsuba | 6S | 13A + 2B | -# | Tom-Cook-3x | 5S | 33A + 2B | +# | Tom-Cook-3x | 5S | 33A + 2B | (Miller loop only, factor 6) # | CH-SQR1 | 3M + 2S | 11A + 2B | # | CH-SQR2 | 2M + 3S | 10A + 2B | # | CH-SQR3 | 1M + 4S | 11A + 2B + 1 Div2 | -# | CH-SQR3x | 1M + 4S | 14A + 2B | +# | CH-SQR3x | 1M + 4S | 14A + 2B | (Miller loop only, factor 2) func square_Chung_Hasan_SQR2(r: var CubicExt, a: CubicExt) {.used.}= ## Returns r = a² @@ -1545,7 +1587,8 @@ func square_Chung_Hasan_SQR3(r: var CubicExt, a: CubicExt) = # ------------------------------------------------------------------- func prodImpl(r: var CubicExt, a, b: CubicExt) = - ## Returns r = a * b # Algorithm is Karatsuba + ## Returns r = a * b + ## Algorithm is Karatsuba var v0{.noInit.}, v1{.noInit.}, v2{.noInit.}: typeof(r.c0) var t0{.noInit.}, t1{.noInit.}, t2{.noInit.}: typeof(r.c0) @@ -1721,6 +1764,58 @@ func invImpl(r: var CubicExt, a: CubicExt) = r.c1.prod(B, t) r.c2.prod(C, t) +func inv2xImpl(r: var CubicExt, a: CubicExt) = + ## Compute the multiplicative inverse of ``a`` + ## via lazy reduction + ## + ## The inverse of 0 is 0. + ## Incidentally this avoids extra check + ## to convert Jacobian and Projective coordinates + ## to affine for elliptic curve + var aa{.noInit.}, bb{.noInit.}, cc{.noInit.}: doublePrec(typeof(r.c0)) + var ab{.noInit.}, bc{.noInit.}, ac{.noInit.}: doublePrec(typeof(r.c0)) + var A {.noInit.}, B {.noInit.}, C {.noInit.}: typeof(r.c0) + var t {.noInit.}, t2{.noInit.}: doublePrec(typeof(r.c0)) + var f {.noInit.}: typeof(r.c0) + + aa.square2x(a.c0) + bb.square2x(a.c1) + cc.square2x(a.c2) + ab.prod2x(a.c0, a.c1) + bc.prod2x(a.c1, a.c2) + ac.prod2x(a.c0, a.c2) + + # A <- a₀² - β a₁ a₂ + t.prod2x(bc, NonResidue) + t.diff2xMod(aa, t) + A.redc2x(t) + + # B <- β a₂² - a₀ a₁ + t.prod2x(cc, NonResidue) + t.diff2xMod(t, ab) + B.redc2x(t) + + # C <- a₁² - a₀ a₂ + t.diff2xMod(bb, ac) + C.redc2x(t) + + # F in t + # F <- β a₁ C + a₀ A + β a₂ B + t.prod2x(B, a.c2) + t2.prod2x(C, a.c1) + t.sum2xMod(t, t2) + t.prod2x(t, NonResidue) + t2.prod2x(A, a.c0) + t.sum2xUnr(t, t2) + f.redc2x(t) + + f.inv() + + # (a0 + a1 v + a2 v²)^-1 = (A + B v + C v²) / F + r.c0.prod(A, f) + r.c1.prod(B, f) + r.c2.prod(C, f) + # Dispatch # ---------------------------------------------------------------------- @@ -1728,8 +1823,7 @@ func invImpl(r: var CubicExt, a: CubicExt) = func square*(r: var CubicExt, a: CubicExt) = ## Returns r = a² - when CubicExt.F.C == BW6_761 or # Too large - CubicExt.F.C == BN254_Snarks: # 50 cycles slower on Fp2->Fp4->Fp12 towering + when CubicExt.C.has_large_NR_norm() or CubicExt.C.has_large_field_elem(): square_Chung_Hasan_SQR3(r, a) else: var d {.noInit.}: doublePrec(typeof(a)) @@ -1747,8 +1841,8 @@ func square2x*(r: var CubicExt2x, a: CubicExt) = square2x_Chung_Hasan_SQR2(r, a) func prod*(r: var CubicExt, a, b: CubicExt) = - ## In-place multiplication - when CubicExt.F.C == BW6_761: # Too large + ## Out-of-place multiplication + when CubicExt.C.has_large_field_elem(): r.prodImpl(a, b) else: var d {.noInit.}: doublePrec(typeof(r)) @@ -1763,14 +1857,7 @@ func prod2x*(r: var CubicExt2x, a, b: CubicExt) = func `*=`*(a: var CubicExt, b: CubicExt) = ## In-place multiplication - when CubicExt.F.C == BW6_761: # Too large - a.prodImpl(a, b) - else: - var d {.noInit.}: doublePrec(typeof(a)) - d.prod2x(a, b) - a.c0.redc2x(d.c0) - a.c1.redc2x(d.c1) - a.c2.redc2x(d.c2) + a.prod(a, b) func inv*(r: var CubicExt, a: CubicExt) = ## Compute the multiplicative inverse of ``a`` @@ -1779,7 +1866,10 @@ func inv*(r: var CubicExt, a: CubicExt) = ## Incidentally this avoids extra check ## to convert Jacobian and Projective coordinates ## to affine for elliptic curve - r.invImpl(a) + when true: + r.invImpl(a) + else: + r.inv2xImpl(a) func inv*(a: var CubicExt) = ## Compute the multiplicative inverse of ``a`` diff --git a/docs/optimizations.md b/docs/optimizations.md index 5fb2ae2..1ce4fcc 100644 --- a/docs/optimizations.md +++ b/docs/optimizations.md @@ -191,7 +191,7 @@ The optimizations can be of algebraic, algorithmic or "implementation details" n - Final exponentiation - [x] Cyclotomic squaring - - [ ] Karabina's compressed cyclotomic squarings + - [x] Karabina's compressed cyclotomic squarings - [x] Addition-chain for exponentiation by curve parameter - [x] BN curves: Fuentes-Castañeda - [ ] BN curves: Duquesne, Ghammam diff --git a/research/kzg_poly_commit/polynomials.nim b/research/kzg_poly_commit/polynomials.nim index dd12ea2..e88a84c 100644 --- a/research/kzg_poly_commit/polynomials.nim +++ b/research/kzg_poly_commit/polynomials.nim @@ -10,7 +10,7 @@ import ../../constantine/pairing/[ pairing_bls12, miller_loops, - cyclotomic_fp12 + cyclotomic_subgroup ] type diff --git a/tests/t_pairing_bls12_377_line_functions.nim b/tests/t_pairing_bls12_377_line_functions.nim index 45771b6..1e50b0b 100644 --- a/tests/t_pairing_bls12_377_line_functions.nim +++ b/tests/t_pairing_bls12_377_line_functions.nim @@ -19,7 +19,7 @@ import ec_shortweierstrass_affine, ec_shortweierstrass_projective, ec_scalar_mul], - ../constantine/pairing/lines_projective, + ../constantine/pairing/lines_eval, # Test utilities ../helpers/[prng_unsafe, static_for] diff --git a/tests/t_pairing_bls12_381_line_functions.nim b/tests/t_pairing_bls12_381_line_functions.nim index 1351404..840b2a7 100644 --- a/tests/t_pairing_bls12_381_line_functions.nim +++ b/tests/t_pairing_bls12_381_line_functions.nim @@ -19,7 +19,7 @@ import ec_shortweierstrass_affine, ec_shortweierstrass_projective, ec_scalar_mul], - ../constantine/pairing/lines_projective, + ../constantine/pairing/lines_eval, # Test utilities ../helpers/[prng_unsafe, static_for] diff --git a/tests/t_pairing_cyclotomic_fp12.nim b/tests/t_pairing_cyclotomic_subgroup.nim similarity index 73% rename from tests/t_pairing_cyclotomic_fp12.nim rename to tests/t_pairing_cyclotomic_subgroup.nim index 845a599..e3e4c77 100644 --- a/tests/t_pairing_cyclotomic_fp12.nim +++ b/tests/t_pairing_cyclotomic_subgroup.nim @@ -14,8 +14,8 @@ import ../constantine/[arithmetic, primitives], ../constantine/towers, ../constantine/config/curves, - ../constantine/io/io_towers, - ../constantine/pairing/cyclotomic_fp12, + ../constantine/io/[io_bigints, io_towers], + ../constantine/pairing/cyclotomic_subgroup, ../constantine/isogeny/frobenius, # Test utilities ../helpers/[prng_unsafe, static_for] @@ -125,3 +125,43 @@ suite "Pairing - Cyclotomic subgroup - GΦ₁₂(p) = {α ∈ Fp¹² : α^Φ₁ test_cycl_squaring_out_place(curve, gen = Uniform) test_cycl_squaring_out_place(curve, gen = HighHammingWeight) test_cycl_squaring_out_place(curve, gen = Long01Sequence) + + test "Compressed cyclotomic squarings": + proc test_compressed_cycl_squarings(C: static Curve, gen: static RandomGen) = + for _ in 0 ..< Iters: + var f = rng.random_elem(Fp12[C], gen) + + f.finalExpEasy() + var g = f + + f.cycl_sqr_repeated(55) + g.cyclotomic_exp_compressed(g, [55]) + + check: bool(f == g) + + staticFor(curve, TestCurves): + test_compressed_cycl_squarings(curve, gen = Uniform) + test_compressed_cycl_squarings(curve, gen = HighHammingWeight) + test_compressed_cycl_squarings(curve, gen = Long01Sequence) + + test "Compressed cyclotomic exponentiation": + proc test_compressed_cycl_exp(C: static Curve, gen: static RandomGen) = + for _ in 0 ..< Iters: + var f = rng.random_elem(Fp12[C], gen) + + f.finalExpEasy() + var g = f + let f2 = f + + # 0b1000000000001000000000000000000000000000000010000000000000000 + const e = BigInt[61].fromHex"0x1001000000010000" + + f.cyclotomic_exp(f2, e, invert = false) + g.cyclotomic_exp_compressed(g, [16, 32, 12]) + + check: bool(f == g) + + staticFor(curve, TestCurves): + test_compressed_cycl_exp(curve, gen = Uniform) + test_compressed_cycl_exp(curve, gen = HighHammingWeight) + test_compressed_cycl_exp(curve, gen = Long01Sequence) \ No newline at end of file diff --git a/tests/t_pairing_mul_fp12_by_lines.nim b/tests/t_pairing_mul_fp12_by_lines.nim index d057142..6563df0 100644 --- a/tests/t_pairing_mul_fp12_by_lines.nim +++ b/tests/t_pairing_mul_fp12_by_lines.nim @@ -15,10 +15,7 @@ import ../constantine/towers, ../constantine/config/curves, ../constantine/io/io_towers, - ../constantine/pairing/[ - lines_projective, - mul_fp12_by_lines - ], + ../constantine/pairing/lines_eval, # Test utilities ../helpers/[prng_unsafe, static_for] diff --git a/tests/t_pairing_template.nim b/tests/t_pairing_template.nim index 48823c8..3804027 100644 --- a/tests/t_pairing_template.nim +++ b/tests/t_pairing_template.nim @@ -16,7 +16,7 @@ import ../constantine/config/curves, ../constantine/elliptic/[ec_shortweierstrass_affine, ec_shortweierstrass_projective], ../constantine/curves/[zoo_subgroups, zoo_pairings], - ../constantine/pairing/[cyclotomic_fp12, cyclotomic_fp6], + ../constantine/pairing/cyclotomic_subgroup, ../constantine/io/io_towers, # Test utilities @@ -27,7 +27,7 @@ export ec_shortweierstrass_affine, ec_shortweierstrass_projective, arithmetic, towers, primitives, io_towers, - cyclotomic_fp12, cyclotomic_fp6 + cyclotomic_subgroup type RandomGen* = enum