status-im/constantine
2
0
Fork 0
mirror of https://github.com/status-im/constantine.git synced 2025-02-20 16:08:34 +00:00
Code Issues Packages Projects Releases Wiki Activity

Pairings optimizations (#178)

Browse Source 
* bench for cyclotomic square, exp and rename cyclotomic exp + multipairings for BLS12-377

* refactor/unify lines and cyclotomic functions

* Add Karabina's compressed squaring

* Use compressed squarings in final exponentiation

* Weighted addchain for bn254_snarks

* Add new towering options and cost functions

* Rearrange bench summaries

* fix BW6-761
This commit is contained in:
Mamy Ratsimbazafy 2022-02-20 20:15:20 +01:00 committed by GitHub
parent f2d51a3b6e
commit dc73c71801
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
47 changed files with 2208 additions and 1852 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
benchmarks/bench_fields_template.nim
View File

@ -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)

3
benchmarks/bench_fp.nim
View File

@ -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()

5
benchmarks/bench_fp2.nim
View File

@ -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()

5
benchmarks/bench_fp4.nim
View File

@ -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()

5
benchmarks/bench_fp6.nim
View File

@ -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()

15
benchmarks/bench_pairing_bls12_377.nim
View File

@ -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()

14
benchmarks/bench_pairing_bls12_381.nim
View File

@ -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()

8
benchmarks/bench_pairing_bn254_nogami.nim
View File

@ -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()

8
benchmarks/bench_pairing_bn254_snarks.nim
View File

@ -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()

57
benchmarks/bench_pairing_template.nim
View File

@ -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) =

15
benchmarks/bench_summary_bls12_377.nim
View File

@ -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)

16
benchmarks/bench_summary_bls12_381.nim
View File

@ -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)

15
benchmarks/bench_summary_bn254_nogami.nim
View File

@ -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)

15
benchmarks/bench_summary_bn254_snarks.nim
View File

@ -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)

2
benchmarks/bench_summary_template.nim
View File

@ -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
],

2
constantine.nimble
View File

@ -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),

11
constantine/arithmetic/assembly/limbs_asm_mul_mont_x86_adx_bmi2.nim
View File

@ -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
# ------------------------------------------------------------

46
constantine/curves/bls12_377_pairing.nim
View File

@ -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

51
constantine/curves/bls12_381_pairing.nim
View File

@ -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

17
constantine/curves/bn254_nogami_pairing.nim
View File

@ -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²)

131
constantine/curves/bn254_snarks_pairing.nim
View File

@ -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²)

16
constantine/curves/bw6_761_pairing.nim
View File

@ -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⁴

2
constantine/curves/zoo_pairings.nim
View File

@ -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

9
constantine/io/io_ec.nim
View File

@ -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].}=

18
constantine/io/io_towers.nim
View File

@ -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

233
constantine/pairing/cyclotomic_fp12.nim
View File

@ -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

233
constantine/pairing/cyclotomic_fp6.nim
View File

@ -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

680
constantine/pairing/cyclotomic_subgroup.nim Normal file
View File

@ -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

74
constantine/pairing/lines_common.nim
View File

@ -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

988
constantine/pairing/lines_eval.nim Normal file
View File

@ -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".}

320
constantine/pairing/lines_projective.nim
View File

@ -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

47
constantine/pairing/miller_loops.nim
View File

@ -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)

564
constantine/pairing/mul_fp12_by_lines.nim
View File

@ -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".}

143
constantine/pairing/mul_fp6_by_lines.nim
View File

@ -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".}

16
constantine/pairing/pairing_bls12.nim
View File

@ -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))

11
constantine/pairing/pairing_bn.nim
View File

@ -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):

15
constantine/pairing/pairing_bw6_761.nim
View File

@ -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

2
constantine/protocols/ethereum_evm_precompiles.nim
View File

@ -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]

8
constantine/tower_field_extensions/assembly/fp2_asm_x86_adx_bmi2.nim
View File

@ -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],

164
constantine/tower_field_extensions/extension_fields.nim
View File

@ -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``

2
docs/optimizations.md
View File

@ -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

2
research/kzg_poly_commit/polynomials.nim
View File

@ -10,7 +10,7 @@ import
../../constantine/pairing/[
pairing_bls12,
miller_loops,
cyclotomic_fp12
cyclotomic_subgroup
]
type

2
tests/t_pairing_bls12_377_line_functions.nim
View File

@ -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]

2
tests/t_pairing_bls12_381_line_functions.nim
View File

@ -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]

44
tests/t_pairing_cyclotomic_fp12.nim → tests/t_pairing_cyclotomic_subgroup.nim
View File

@ -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)

5
tests/t_pairing_mul_fp12_by_lines.nim
View File

@ -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]

4
tests/t_pairing_template.nim
View File

@ -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