From 842b934d7549c89aac5da5ca13633753f2e000f2 Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 14 Jun 2024 13:15:40 +0100 Subject: [PATCH 1/2] perf(bls12-381): eliminate finalexp ~naively --- .../emulated/fields_bls12381/e12_pairing.go | 184 ++++++++++++++++++ std/algebra/emulated/fields_bls12381/hints.go | 44 +++++ std/algebra/emulated/sw_bls12381/pairing.go | 5 +- 3 files changed, 230 insertions(+), 3 deletions(-) diff --git a/std/algebra/emulated/fields_bls12381/e12_pairing.go b/std/algebra/emulated/fields_bls12381/e12_pairing.go index 4cd3dc5613..2c2fe1d73b 100644 --- a/std/algebra/emulated/fields_bls12381/e12_pairing.go +++ b/std/algebra/emulated/fields_bls12381/e12_pairing.go @@ -2,6 +2,122 @@ package fields_bls12381 import "github.com/consensys/gnark/std/math/emulated" +func (e Ext12) nSquare(z *E12, n int) *E12 { + for i := 0; i < n; i++ { + z = e.Square(z) + } + return z +} + +// Expt sets z to x^t in E12 and return z +// where t = -u = 15132376222941642752 +func (e Ext12) Expt(x *E12) *E12 { + // Expt computation is derived from the addition chain: + // + // _10 = 2*1 + // _11 = 1 + _10 + // _1100 = _11 << 2 + // _1101 = 1 + _1100 + // _1101000 = _1101 << 3 + // _1101001 = 1 + _1101000 + // return ((_1101001 << 9 + 1) << 32 + 1) << 15 + // + // Operations: 62 squares 5 multiplies + // + // Generated by github.com/mmcloughlin/addchain v0.4.0. + + z := e.Square(x) + z = e.Mul(x, z) + z = e.Square(z) + z = e.Square(z) + z = e.Mul(x, z) + z = e.nSquare(z, 3) + z = e.Mul(x, z) + z = e.nSquare(z, 9) + z = e.Mul(x, z) + z = e.nSquare(z, 32) + z = e.Mul(x, z) + z = e.nSquare(z, 15) + z = e.Square(z) + + return z +} + +// ExpByU sets z to x^U in E12 and return z +// where U = (u-1)^2/3 = 76329603384216526031706109802092473003 +func (e Ext12) ExpByU(x *E12) *E12 { + // ExpByU computation is derived from the addition chain: + // + // _10 = 2*1 + // _11 = 1 + _10 + // _110 = 2*_11 + // _111 = 1 + _110 + // _1000 = 1 + _111 + // _100000 = _1000 << 2 + // _100011 = _11 + _100000 + // _101010 = _111 + _100011 + // _1010100 = 2*_101010 + // _1010101 = 1 + _1010100 + // _1010110 = 1 + _1010101 + // _1011001 = _11 + _1010110 + // _1100001 = _1000 + _1011001 + // _10101011 = _1010101 + _1010110 + // _11000010 = 2*_1100001 + // _11100101 = _100011 + _11000010 + // i49 = ((_11100101 << 7 + _1011001) << 5 + _11) << 18 + // i68 = ((_1010101 + i49) << 9 + _10101011) << 7 + _1100001 + // i85 = (i68 << 9 + _10 + _10101011) << 5 + _11 + // i136 = ((i85 << 17 + _1010101) << 9 + _10101011) << 23 + // return (_1010101 + i136) << 9 + _10101011 + // + // Operations: 124 squares 23 multiplies + // + // Generated by github.com/mmcloughlin/addchain v0.4.0. + + t2 := e.Square(x) + t1 := e.Mul(x, t2) + z := e.Square(t1) + z = e.Mul(x, z) + t3 := e.Mul(x, z) + t0 := e.Square(t3) + t0 = e.Square(t0) + t5 := e.Mul(t1, t0) + z = e.Mul(z, t5) + z = e.Square(z) + t0 = e.Mul(x, z) + z = e.Mul(x, t0) + t4 := e.Mul(t1, z) + t3 = e.Mul(t3, t4) + z = e.Mul(t0, z) + t6 := e.Square(t3) + t5 = e.Mul(t5, t6) + t5 = e.nSquare(t5, 7) + t4 = e.Mul(t4, t5) + t4 = e.nSquare(t4, 5) + t4 = e.Mul(t1, t4) + t4 = e.nSquare(t4, 18) + t4 = e.Mul(t0, t4) + t4 = e.nSquare(t4, 9) + t4 = e.Mul(z, t4) + t4 = e.nSquare(t4, 7) + t3 = e.Mul(t3, t4) + t3 = e.nSquare(t3, 9) + t2 = e.Mul(t2, t3) + t2 = e.Mul(z, t2) + t2 = e.nSquare(t2, 5) + t1 = e.Mul(t1, t2) + t1 = e.nSquare(t1, 17) + t1 = e.Mul(t0, t1) + t1 = e.nSquare(t1, 9) + t1 = e.Mul(z, t1) + t1 = e.nSquare(t1, 23) + t0 = e.Mul(t0, t1) + t0 = e.nSquare(t0, 9) + z = e.Mul(z, t0) + + return z +} + func (e Ext12) nSquareTorus(z *E6, n int) *E6 { for i := 0; i < n; i++ { z = e.SquareTorus(z) @@ -268,3 +384,71 @@ func (e Ext12) FrobeniusSquareTorus(y *E6) *E6 { return &E6{B0: *t0, B1: *t1, B2: *t2} } + +// FinalExponentiationCheck checks that a Miller function output x lies in the +// same equivalence class as the reduced pairing. This replaces the final +// exponentiation step in-circuit. +// The method follows Section 4 of [On Proving Pairings] paper by A. Novakovic and L. Eagen. +// +// [On Proving Pairings]: https://eprint.iacr.org/2024/640.pdf +func (e Ext12) FinalExponentiationCheck(x *E12) *E12 { + res, err := e.fp.NewHint(finalExpHint, 12, &x.C0.B0.A0, &x.C0.B0.A1, &x.C0.B1.A0, &x.C0.B1.A1, &x.C0.B2.A0, &x.C0.B2.A1, &x.C1.B0.A0, &x.C1.B0.A1, &x.C1.B1.A0, &x.C1.B1.A1, &x.C1.B2.A0, &x.C1.B2.A1) + if err != nil { + // err is non-nil only for invalid number of inputs + panic(err) + } + + residueWitness := E12{ + C0: E6{ + B0: E2{A0: *res[0], A1: *res[1]}, + B1: E2{A0: *res[2], A1: *res[3]}, + B2: E2{A0: *res[4], A1: *res[5]}, + }, + C1: E6{ + B0: E2{A0: *res[6], A1: *res[7]}, + B1: E2{A0: *res[8], A1: *res[9]}, + B2: E2{A0: *res[10], A1: *res[11]}, + }, + } + + // Check that x == residueWitness^r by checking that: + // x^k == residueWitness^(q-u) + // where k = (u-1)^2/3, u=-0xd201000000010000 the BLS12-381 seed + // and residueWitness from the hint. + t0 := e.Frobenius(&residueWitness) + // exponentiation by -u + t1 := e.Expt(&residueWitness) + t0 = e.Mul(t0, t1) + // exponentiation by U=(u-1)^2/3 + t1 = e.ExpByU(x) + + e.AssertIsEqual(t0, t1) + + return nil +} + +func (e Ext12) Frobenius(x *E12) *E12 { + t0 := e.Ext2.Conjugate(&x.C0.B0) + t1 := e.Ext2.Conjugate(&x.C0.B1) + t2 := e.Ext2.Conjugate(&x.C0.B2) + t3 := e.Ext2.Conjugate(&x.C1.B0) + t4 := e.Ext2.Conjugate(&x.C1.B1) + t5 := e.Ext2.Conjugate(&x.C1.B2) + t1 = e.Ext2.MulByNonResidue1Power2(t1) + t2 = e.Ext2.MulByNonResidue1Power4(t2) + t3 = e.Ext2.MulByNonResidue1Power1(t3) + t4 = e.Ext2.MulByNonResidue1Power3(t4) + t5 = e.Ext2.MulByNonResidue1Power5(t5) + return &E12{ + C0: E6{ + B0: *t0, + B1: *t1, + B2: *t2, + }, + C1: E6{ + B0: *t3, + B1: *t4, + B2: *t5, + }, + } +} diff --git a/std/algebra/emulated/fields_bls12381/hints.go b/std/algebra/emulated/fields_bls12381/hints.go index fdc9700504..320aaacb9d 100644 --- a/std/algebra/emulated/fields_bls12381/hints.go +++ b/std/algebra/emulated/fields_bls12381/hints.go @@ -27,6 +27,7 @@ func GetHints() []solver.Hint { // E12 divE12Hint, inverseE12Hint, + finalExpHint, } } @@ -268,3 +269,46 @@ func divE12Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) erro return nil }) } + +func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { + // This follows section 4.1 of https://eprint.iacr.org/2024/640.pdf (Th. 1) + return emulated.UnwrapHint(nativeInputs, nativeOutputs, + func(mod *big.Int, inputs, outputs []*big.Int) error { + var millerLoop, residueWitness bls12381.E12 + var rInv big.Int + + millerLoop.C0.B0.A0.SetBigInt(inputs[0]) + millerLoop.C0.B0.A1.SetBigInt(inputs[1]) + millerLoop.C0.B1.A0.SetBigInt(inputs[2]) + millerLoop.C0.B1.A1.SetBigInt(inputs[3]) + millerLoop.C0.B2.A0.SetBigInt(inputs[4]) + millerLoop.C0.B2.A1.SetBigInt(inputs[5]) + millerLoop.C1.B0.A0.SetBigInt(inputs[6]) + millerLoop.C1.B0.A1.SetBigInt(inputs[7]) + millerLoop.C1.B1.A0.SetBigInt(inputs[8]) + millerLoop.C1.B1.A1.SetBigInt(inputs[9]) + millerLoop.C1.B2.A0.SetBigInt(inputs[10]) + millerLoop.C1.B2.A1.SetBigInt(inputs[11]) + + // compute r-th root: + // Exponentiate to rInv where + // rInv = 1/r mod (p^12-1)/r + rInv.SetString("169662389312441398885310937191698694666993326870281216192803558492181163400934408837135364582394949149589560242411491538960982200559697133935443307582773537814554128992403254243871087441488619811839498788505657962013599019994544063402394719913759780901881538869078447034832302535303591303383830742161317593225991746471557492001710830538428792119562309446698444646787667517629943447802199824630112988907247336627481159245442124709621313522294197747687500252452962523217400829932174349352696726049683687654879009114460723993703760367089269403767790334911644010940272722630305066645230222732316445557889124653426141642271480304669447694344127599708992364443461893123938202386892312748211835322692697497854107961493711137028209148238339237355911496376520814450515612396561384525661635220451168152178239892009375229296874955612623691164738926395993739297557487207643426168321070539996994036837992284584225139752716615623194417718962478029165908544042568334172107008712033983002554672734519081879196926275059798317879322062358113986901925780890205936071364647548199159506709147492864081514759663116291487638998943660232689862634717010538047493292265992334130695994203833154950619462266484292385471162124464248375625748097868775829652908052615424796255913420292818674303286242639225711610323988077268116737", 10) + residueWitness.Exp(millerLoop, &rInv) + + residueWitness.C0.B0.A0.BigInt(outputs[0]) + residueWitness.C0.B0.A1.BigInt(outputs[1]) + residueWitness.C0.B1.A0.BigInt(outputs[2]) + residueWitness.C0.B1.A1.BigInt(outputs[3]) + residueWitness.C0.B2.A0.BigInt(outputs[4]) + residueWitness.C0.B2.A1.BigInt(outputs[5]) + residueWitness.C1.B0.A0.BigInt(outputs[6]) + residueWitness.C1.B0.A1.BigInt(outputs[7]) + residueWitness.C1.B1.A0.BigInt(outputs[8]) + residueWitness.C1.B1.A1.BigInt(outputs[9]) + residueWitness.C1.B2.A0.BigInt(outputs[10]) + residueWitness.C1.B2.A1.BigInt(outputs[11]) + + return nil + }) +} diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go index bfe908cbd7..47c0a185aa 100644 --- a/std/algebra/emulated/sw_bls12381/pairing.go +++ b/std/algebra/emulated/sw_bls12381/pairing.go @@ -245,13 +245,12 @@ func (pr Pairing) Pair(P []*G1Affine, Q []*G2Affine) (*GTEl, error) { // // This function doesn't check that the inputs are in the correct subgroups. func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error { - f, err := pr.Pair(P, Q) + f, err := pr.MillerLoop(P, Q) if err != nil { return err } - one := pr.One() - pr.AssertIsEqual(f, one) + pr.FinalExponentiationCheck(f) return nil } From ce537580ab5dc6885c28bdcf370d1c894679205e Mon Sep 17 00:00:00 2001 From: Youssef El Housni Date: Fri, 14 Jun 2024 17:00:26 +0100 Subject: [PATCH 2/2] perf(bls12-381): use cyclotomic group in finel exp check --- std/algebra/emulated/fields_bls12381/e12.go | 59 +++++++++++++++++ .../emulated/fields_bls12381/e12_pairing.go | 66 +++++++++++-------- std/algebra/emulated/sw_bls12381/pairing.go | 10 +++ 3 files changed, 108 insertions(+), 27 deletions(-) diff --git a/std/algebra/emulated/fields_bls12381/e12.go b/std/algebra/emulated/fields_bls12381/e12.go index 84728728c2..dbc394db1a 100644 --- a/std/algebra/emulated/fields_bls12381/e12.go +++ b/std/algebra/emulated/fields_bls12381/e12.go @@ -114,6 +114,65 @@ func (e Ext12) Square(x *E12) *E12 { } } +// Granger--Scott cyclotomic square +func (e Ext12) CyclotomicSquare(x *E12) *E12 { + t0 := e.Ext2.Square(&x.C1.B1) + t1 := e.Ext2.Square(&x.C0.B0) + t6 := e.Ext2.Add(&x.C1.B1, &x.C0.B0) + t6 = e.Ext2.Square(t6) + t6 = e.Ext2.Sub(t6, t0) + t6 = e.Ext2.Sub(t6, t1) + t2 := e.Ext2.Square(&x.C0.B2) + t3 := e.Ext2.Square(&x.C1.B0) + t7 := e.Ext2.Add(&x.C0.B2, &x.C1.B0) + t7 = e.Ext2.Square(t7) + t7 = e.Ext2.Sub(t7, t2) + t7 = e.Ext2.Sub(t7, t3) + t4 := e.Ext2.Square(&x.C1.B2) + t5 := e.Ext2.Square(&x.C0.B1) + t8 := e.Ext2.Add(&x.C1.B2, &x.C0.B1) + t8 = e.Ext2.Square(t8) + t8 = e.Ext2.Sub(t8, t4) + t8 = e.Ext2.Sub(t8, t5) + t8 = e.Ext2.MulByNonResidue(t8) + t0 = e.Ext2.MulByNonResidue(t0) + t0 = e.Ext2.Add(t0, t1) + t2 = e.Ext2.MulByNonResidue(t2) + t2 = e.Ext2.Add(t2, t3) + t4 = e.Ext2.MulByNonResidue(t4) + t4 = e.Ext2.Add(t4, t5) + z00 := e.Ext2.Sub(t0, &x.C0.B0) + z00 = e.Ext2.Double(z00) + z00 = e.Ext2.Add(z00, t0) + z01 := e.Ext2.Sub(t2, &x.C0.B1) + z01 = e.Ext2.Double(z01) + z01 = e.Ext2.Add(z01, t2) + z02 := e.Ext2.Sub(t4, &x.C0.B2) + z02 = e.Ext2.Double(z02) + z02 = e.Ext2.Add(z02, t4) + z10 := e.Ext2.Add(t8, &x.C1.B0) + z10 = e.Ext2.Double(z10) + z10 = e.Ext2.Add(z10, t8) + z11 := e.Ext2.Add(t6, &x.C1.B1) + z11 = e.Ext2.Double(z11) + z11 = e.Ext2.Add(z11, t6) + z12 := e.Ext2.Add(t7, &x.C1.B2) + z12 = e.Ext2.Double(z12) + z12 = e.Ext2.Add(z12, t7) + return &E12{ + C0: E6{ + B0: *z00, + B1: *z01, + B2: *z02, + }, + C1: E6{ + B0: *z10, + B1: *z11, + B2: *z12, + }, + } +} + func (e Ext12) AssertIsEqual(x, y *E12) { e.Ext6.AssertIsEqual(&x.C0, &y.C0) e.Ext6.AssertIsEqual(&x.C1, &y.C1) diff --git a/std/algebra/emulated/fields_bls12381/e12_pairing.go b/std/algebra/emulated/fields_bls12381/e12_pairing.go index 2c2fe1d73b..7bbb60b7e5 100644 --- a/std/algebra/emulated/fields_bls12381/e12_pairing.go +++ b/std/algebra/emulated/fields_bls12381/e12_pairing.go @@ -2,9 +2,9 @@ package fields_bls12381 import "github.com/consensys/gnark/std/math/emulated" -func (e Ext12) nSquare(z *E12, n int) *E12 { +func (e Ext12) nSquareGS(z *E12, n int) *E12 { for i := 0; i < n; i++ { - z = e.Square(z) + z = e.CyclotomicSquare(z) } return z } @@ -26,19 +26,18 @@ func (e Ext12) Expt(x *E12) *E12 { // // Generated by github.com/mmcloughlin/addchain v0.4.0. - z := e.Square(x) + z := e.CyclotomicSquare(x) z = e.Mul(x, z) - z = e.Square(z) - z = e.Square(z) + z = e.nSquareGS(z, 2) z = e.Mul(x, z) - z = e.nSquare(z, 3) + z = e.nSquareGS(z, 3) z = e.Mul(x, z) - z = e.nSquare(z, 9) + z = e.nSquareGS(z, 9) z = e.Mul(x, z) - z = e.nSquare(z, 32) + z = e.nSquareGS(z, 32) z = e.Mul(x, z) - z = e.nSquare(z, 15) - z = e.Square(z) + z = e.nSquareGS(z, 15) + z = e.CyclotomicSquare(z) return z } @@ -74,45 +73,45 @@ func (e Ext12) ExpByU(x *E12) *E12 { // // Generated by github.com/mmcloughlin/addchain v0.4.0. - t2 := e.Square(x) + t2 := e.CyclotomicSquare(x) t1 := e.Mul(x, t2) - z := e.Square(t1) + z := e.CyclotomicSquare(t1) z = e.Mul(x, z) t3 := e.Mul(x, z) - t0 := e.Square(t3) - t0 = e.Square(t0) + t0 := e.CyclotomicSquare(t3) + t0 = e.CyclotomicSquare(t0) t5 := e.Mul(t1, t0) z = e.Mul(z, t5) - z = e.Square(z) + z = e.CyclotomicSquare(z) t0 = e.Mul(x, z) z = e.Mul(x, t0) t4 := e.Mul(t1, z) t3 = e.Mul(t3, t4) z = e.Mul(t0, z) - t6 := e.Square(t3) + t6 := e.CyclotomicSquare(t3) t5 = e.Mul(t5, t6) - t5 = e.nSquare(t5, 7) + t5 = e.nSquareGS(t5, 7) t4 = e.Mul(t4, t5) - t4 = e.nSquare(t4, 5) + t4 = e.nSquareGS(t4, 5) t4 = e.Mul(t1, t4) - t4 = e.nSquare(t4, 18) + t4 = e.nSquareGS(t4, 18) t4 = e.Mul(t0, t4) - t4 = e.nSquare(t4, 9) + t4 = e.nSquareGS(t4, 9) t4 = e.Mul(z, t4) - t4 = e.nSquare(t4, 7) + t4 = e.nSquareGS(t4, 7) t3 = e.Mul(t3, t4) - t3 = e.nSquare(t3, 9) + t3 = e.nSquareGS(t3, 9) t2 = e.Mul(t2, t3) t2 = e.Mul(z, t2) - t2 = e.nSquare(t2, 5) + t2 = e.nSquareGS(t2, 5) t1 = e.Mul(t1, t2) - t1 = e.nSquare(t1, 17) + t1 = e.nSquareGS(t1, 17) t1 = e.Mul(t0, t1) - t1 = e.nSquare(t1, 9) + t1 = e.nSquareGS(t1, 9) t1 = e.Mul(z, t1) - t1 = e.nSquare(t1, 23) + t1 = e.nSquareGS(t1, 23) t0 = e.Mul(t0, t1) - t0 = e.nSquare(t0, 9) + t0 = e.nSquareGS(t0, 9) z = e.Mul(z, t0) return z @@ -452,3 +451,16 @@ func (e Ext12) Frobenius(x *E12) *E12 { }, } } + +func (e Ext12) FrobeniusSquare(x *E12) *E12 { + z00 := &x.C0.B0 + z01 := e.Ext2.MulByNonResidue2Power2(&x.C0.B1) + z02 := e.Ext2.MulByNonResidue2Power4(&x.C0.B2) + z10 := e.Ext2.MulByNonResidue2Power1(&x.C1.B0) + z11 := e.Ext2.MulByNonResidue2Power3(&x.C1.B1) + z12 := e.Ext2.MulByNonResidue2Power5(&x.C1.B2) + return &E12{ + C0: E6{B0: *z00, B1: *z01, B2: *z02}, + C1: E6{B0: *z10, B1: *z11, B2: *z12}, + } +} diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go index 47c0a185aa..c5f96e62d8 100644 --- a/std/algebra/emulated/sw_bls12381/pairing.go +++ b/std/algebra/emulated/sw_bls12381/pairing.go @@ -250,6 +250,16 @@ func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error { return err } + // We perform the easy part of the final exp to push f to the cyclotomic + // subgroup so that FinalExponentiationCheck is carried with optimized + // cyclotomic squaring (e.g. Karabina12345). + // + // f = f^(p⁶-1)(p²+1) + buf := pr.Conjugate(f) + buf = pr.DivUnchecked(buf, f) + f = pr.FrobeniusSquare(buf) + f = pr.Mul(f, buf) + pr.FinalExponentiationCheck(f) return nil