diff --git a/go.mod b/go.mod index cd70e6f6c..caa4a8e03 100644 --- a/go.mod +++ b/go.mod @@ -7,9 +7,9 @@ toolchain go1.22.6 require ( github.com/bits-and-blooms/bitset v1.14.2 github.com/blang/semver/v4 v4.0.0 - github.com/consensys/bavard v0.1.24 + github.com/consensys/bavard v0.1.25 github.com/consensys/compress v0.2.5 - github.com/consensys/gnark-crypto v0.14.1-0.20241217131346-b998989abdbe + github.com/consensys/gnark-crypto v0.14.1-0.20250116204316-e7fd38b0a0a6 github.com/fxamacker/cbor/v2 v2.7.0 github.com/google/go-cmp v0.6.0 github.com/google/pprof v0.0.0-20240727154555-813a5fbdbec8 diff --git a/go.sum b/go.sum index 305b0d264..e20d4967a 100644 --- a/go.sum +++ b/go.sum @@ -57,12 +57,12 @@ github.com/client9/misspell v0.3.4/go.mod h1:qj6jICC3Q7zFZvVWo7KLAzC3yx5G7kyvSDk github.com/cncf/udpa/go v0.0.0-20191209042840-269d4d468f6f/go.mod h1:M8M6+tZqaGXZJjfX53e64911xZQV5JYwmTeXPW+k8Sc= github.com/cncf/udpa/go v0.0.0-20200629203442-efcf912fb354/go.mod h1:WmhPx2Nbnhtbo57+VJT5O0JRkEi1Wbu0z5j0R8u5Hbk= github.com/cncf/udpa/go v0.0.0-20201120205902-5459f2c99403/go.mod h1:WmhPx2Nbnhtbo57+VJT5O0JRkEi1Wbu0z5j0R8u5Hbk= -github.com/consensys/bavard v0.1.24 h1:Lfe+bjYbpaoT7K5JTFoMi5wo9V4REGLvQQbHmatoN2I= -github.com/consensys/bavard v0.1.24/go.mod h1:k/zVjHHC4B+PQy1Pg7fgvG3ALicQw540Crag8qx+dZs= +github.com/consensys/bavard v0.1.25 h1:5YcSBnp03/HvfpKaIQLr/ecspTp2k8YNR5rQLOWvUyc= +github.com/consensys/bavard v0.1.25/go.mod h1:k/zVjHHC4B+PQy1Pg7fgvG3ALicQw540Crag8qx+dZs= github.com/consensys/compress v0.2.5 h1:gJr1hKzbOD36JFsF1AN8lfXz1yevnJi1YolffY19Ntk= github.com/consensys/compress v0.2.5/go.mod h1:pyM+ZXiNUh7/0+AUjUf9RKUM6vSH7T/fsn5LLS0j1Tk= -github.com/consensys/gnark-crypto v0.14.1-0.20241217131346-b998989abdbe h1:WNuXPe50FqynKlUOMdsi1eCzYN8gU4sdCsW3eg3coGA= -github.com/consensys/gnark-crypto v0.14.1-0.20241217131346-b998989abdbe/go.mod h1:ePFa23CZLMRMHxQpY5nMaiAZ3yuEIayaB8ElEvlwLEs= +github.com/consensys/gnark-crypto v0.14.1-0.20250116204316-e7fd38b0a0a6 h1:P4DeR8HYfQGl4Vj6KEv0Eszcokroit/U1dRrUsgt+js= +github.com/consensys/gnark-crypto v0.14.1-0.20250116204316-e7fd38b0a0a6/go.mod h1:q9s22Y0WIHd9UCBfD+xGeW8wDJ7WAGZZpMrLFqzBzrQ= github.com/coreos/go-semver v0.3.0/go.mod h1:nnelYz7RCh+5ahJtPPxZlU+153eP4D4r3EedlOD2RNk= github.com/coreos/go-systemd/v22 v22.3.2/go.mod h1:Y58oyj3AT4RCenI/lSvhwexgC+NSVTIJ3seZv2GcEnc= github.com/coreos/go-systemd/v22 v22.5.0/go.mod h1:Y58oyj3AT4RCenI/lSvhwexgC+NSVTIJ3seZv2GcEnc= diff --git a/internal/tinyfield/element.go b/internal/tinyfield/element.go index 021d9370c..dcaa56b23 100644 --- a/internal/tinyfield/element.go +++ b/internal/tinyfield/element.go @@ -658,7 +658,7 @@ func (z *Element) SetBigInt(v *big.Int) *Element { // v == 0 return z } else if c != 1 && v.Cmp(&zero) != -1 { - // 0 < v < q + // 0 <= v < q return z.setBigInt(v) } @@ -785,6 +785,8 @@ type ByteOrder interface { String() string } +var errInvalidEncoding = errors.New("invalid tinyfield.Element encoding") + // BigEndian is the big-endian implementation of ByteOrder and AppendByteOrder. var BigEndian bigEndian @@ -797,7 +799,7 @@ func (bigEndian) Element(b *[Bytes]byte) (Element, error) { z[0] = binary.BigEndian.Uint64((*b)[0:8]) if !z.smallerThanModulus() { - return Element{}, errors.New("invalid tinyfield.Element encoding") + return Element{}, errInvalidEncoding } z.toMont() @@ -821,7 +823,7 @@ func (littleEndian) Element(b *[Bytes]byte) (Element, error) { z[0] = binary.LittleEndian.Uint64((*b)[0:8]) if !z.smallerThanModulus() { - return Element{}, errors.New("invalid tinyfield.Element encoding") + return Element{}, errInvalidEncoding } z.toMont() diff --git a/std/algebra/emulated/fields_bw6761/e6_pairing.go b/std/algebra/emulated/fields_bw6761/e6_pairing.go index a4a40e840..94e00e4f7 100644 --- a/std/algebra/emulated/fields_bw6761/e6_pairing.go +++ b/std/algebra/emulated/fields_bw6761/e6_pairing.go @@ -423,43 +423,6 @@ func (e *Ext6) mulBy02345(z *E6, x [5]*baseEl) *E6 { } } -// AssertFinalExponentiationIsOne 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 is adapted from 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 Ext6) AssertFinalExponentiationIsOne(x *E6) { - res, err := e.fp.NewHint(finalExpHint, 6, &x.A0, &x.A1, &x.A2, &x.A3, &x.A4, &x.A5) - if err != nil { - // err is non-nil only for invalid number of inputs - panic(err) - } - - residueWitness := E6{ - A0: *res[0], - A1: *res[1], - A2: *res[2], - A3: *res[3], - A4: *res[4], - A5: *res[5], - } - - // Check that x == residueWitness^λ - // where λ = u^3-u^2+1 - (u+1)p, with u the BW6-761 seed - // and residueWitness from the hint. - - // exponentiation by U1=u^3-u^2+1 - t0 := e.ExpByU1(&residueWitness) - // exponentiation by U2=u+1 - t1 := e.ExpByU2(&residueWitness) - - t1 = e.Frobenius(t1) - t0 = e.DivUnchecked(t0, t1) - - e.AssertIsEqual(t0, x) -} - // ExpByU2 set z to z^(x₀+1) in E12 and return z // x₀+1 = 9586122913090633730 func (e Ext6) ExpByU2(z *E6) *E6 { diff --git a/std/algebra/emulated/fields_bw6761/hints.go b/std/algebra/emulated/fields_bw6761/hints.go index 255cc9b2f..994daa0d1 100644 --- a/std/algebra/emulated/fields_bw6761/hints.go +++ b/std/algebra/emulated/fields_bw6761/hints.go @@ -19,7 +19,6 @@ func GetHints() []solver.Hint { divE6Hint, inverseE6Hint, divE6By362880Hint, - finalExpHint, } } @@ -107,41 +106,3 @@ func divE6By362880Hint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.In return nil }) } - -func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { - // This adapted from section 4.3.2 of https://eprint.iacr.org/2024/640.pdf - return emulated.UnwrapHint(nativeInputs, nativeOutputs, - func(mod *big.Int, inputs, outputs []*big.Int) error { - var millerLoop, residueWitness bw6761.E6 - var rInv, mInv big.Int - - millerLoop.B0.A0.SetBigInt(inputs[0]) - millerLoop.B0.A1.SetBigInt(inputs[2]) - millerLoop.B0.A2.SetBigInt(inputs[4]) - millerLoop.B1.A0.SetBigInt(inputs[1]) - millerLoop.B1.A1.SetBigInt(inputs[3]) - millerLoop.B1.A2.SetBigInt(inputs[5]) - - // 1. compute r-th root: - // Exponentiate to rInv where - // rInv = 1/r mod (p^6-1)/r - rInv.SetString("279142441805511726233822077180198394933430419224185936052953462287387912118470357993263103168031788043160461358474005435622327506926362567154401645657309519073154383052970657693950208844465818979551693587858245321454505472049236704031061301292776853925224359757586505231126091244204292668007110271845616234279927419974150119801003450133674289144711275201991607282264849765236206295842916353255855388186086438329721887082685697023028663652777877691341551982676874308309620809049793085180324511691754953492619183755890255644855765188965000691813063771086522132765764526955251054211157804606693386854395171192876178005945476647006847460976477055233044799299417913662363985523123796056692751028712679181978298499780752966303529102009307348414562366180130429432094237007700663759126264893082917308542509779442201840676518234962495304673134599305371982876385622279935346701152286347948653741121231188575146952014672242471261647823749129902237689180055673361938161119768341970519416039779128617354778773830515364777252518313057683396662835013368967463878342754251509207391537635831891662211848811733884861792121210263430418966889668537646457064092991696527814120385172941004264289812969796992647021735186941896252860419364971543301451924917610828019341224722038007513", 10) - residueWitness.Exp(millerLoop, &rInv) - - // 2. compute m-th root: - // where m = (x+1 + x(x^2-x^1-1)q) / r - // Exponentiate to mInv where - // mInv = 1/m mod p^6-1/r - mInv.SetString("420096572758781926988571022578549119077996267041217186563532964653013626327499627643558150289556860284699838191238508062761264485377946319676011525555582097381055209304464769241709045835179375847000286979304653199040198646948595850434830718773056593021324330541604029824826938177546414778934883707126835848724258610612114712835130017082970786784508470382396148858570586085402148355642863720286568566937773459407961735112550507047306343380386401338522186960986251395049985320677251315016812720092326581314645206610216409714397970562842517827716362494341171265008409446148022671451843025093584702610246849007545665518399731546205544005105929880663530772806759681913801835273987094997504640832304570158760940364827187477825525048007459079382410480491250884588399683894539404567701993526561088158396861020181640181843560309670937868772703282755078557149854363818903590441797744966016708880143332350534049482338696654635346189790575286999280892407997722996866724226514621504774811766428733682155766330614074143245300182851212177081558245259537898592443393875891588079021560334726750431309338787970594548465289737362624558256642461612913108676326999205533110217714096123782036214164015261929502119392490941988919030563789520985909704716341786823561745842985678563", 10) - residueWitness.Exp(residueWitness, &mInv) - - residueWitness.B0.A0.BigInt(outputs[0]) - residueWitness.B0.A1.BigInt(outputs[2]) - residueWitness.B0.A2.BigInt(outputs[4]) - residueWitness.B1.A0.BigInt(outputs[1]) - residueWitness.B1.A1.BigInt(outputs[3]) - residueWitness.B1.A2.BigInt(outputs[5]) - - return nil - }) -} diff --git a/std/algebra/emulated/sw_bls12381/hints.go b/std/algebra/emulated/sw_bls12381/hints.go index 676ed3c22..b767113da 100644 --- a/std/algebra/emulated/sw_bls12381/hints.go +++ b/std/algebra/emulated/sw_bls12381/hints.go @@ -14,7 +14,10 @@ func init() { // GetHints returns all hint functions used in the package. func GetHints() []solver.Hint { - return []solver.Hint{finalExpHint} + return []solver.Hint{ + finalExpHint, + pairingCheckHint, + } } func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { @@ -37,69 +40,7 @@ func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) er millerLoop.C1.B2.A0.SetBigInt(inputs[10]) millerLoop.C1.B2.A1.SetBigInt(inputs[11]) - var root, rootPthInverse, root27thInverse, residueWitness, scalingFactor bls12381.E12 - var order3rd, order3rdPower, exponent, exponentInv, finalExpFactor, polyFactor big.Int - // polyFactor = (1-x)/3 - polyFactor.SetString("5044125407647214251", 10) - // finalExpFactor = ((q^12 - 1) / r) / (27 * polyFactor) - finalExpFactor.SetString("2366356426548243601069753987687709088104621721678962410379583120840019275952471579477684846670499039076873213559162845121989217658133790336552276567078487633052653005423051750848782286407340332979263075575489766963251914185767058009683318020965829271737924625612375201545022326908440428522712877494557944965298566001441468676802477524234094954960009227631543471415676620753242466901942121887152806837594306028649150255258504417829961387165043999299071444887652375514277477719817175923289019181393803729926249507024121957184340179467502106891835144220611408665090353102353194448552304429530104218473070114105759487413726485729058069746063140422361472585604626055492939586602274983146215294625774144156395553405525711143696689756441298365274341189385646499074862712688473936093315628166094221735056483459332831845007196600723053356837526749543765815988577005929923802636375670820616189737737304893769679803809426304143627363860243558537831172903494450556755190448279875942974830469855835666815454271389438587399739607656399812689280234103023464545891697941661992848552456326290792224091557256350095392859243101357349751064730561345062266850238821755009430903520645523345000326783803935359711318798844368754833295302563158150573540616830138810935344206231367357992991289265295323280", 10) - - // 1. get pth-root inverse - exponent.Mul(&finalExpFactor, big.NewInt(27)) - root.Exp(millerLoop, &exponent) - if root.IsOne() { - rootPthInverse.SetOne() - } else { - exponentInv.ModInverse(&exponent, &polyFactor) - exponent.Neg(&exponentInv).Mod(&exponent, &polyFactor) - rootPthInverse.Exp(root, &exponent) - } - - // 2.1. get order of 3rd primitive root - var three big.Int - three.SetUint64(3) - exponent.Mul(&polyFactor, &finalExpFactor) - root.Exp(millerLoop, &exponent) - if root.IsOne() { - order3rdPower.SetUint64(0) - } - root.Exp(root, &three) - if root.IsOne() { - order3rdPower.SetUint64(1) - } - root.Exp(root, &three) - if root.IsOne() { - order3rdPower.SetUint64(2) - } - root.Exp(root, &three) - if root.IsOne() { - order3rdPower.SetUint64(3) - } - - // 2.2. get 27th root inverse - if order3rdPower.Uint64() == 0 { - root27thInverse.SetOne() - } else { - order3rd.Exp(&three, &order3rdPower, nil) - exponent.Mul(&polyFactor, &finalExpFactor) - root.Exp(millerLoop, &exponent) - exponentInv.ModInverse(&exponent, &order3rd) - exponent.Neg(&exponentInv).Mod(&exponent, &order3rd) - root27thInverse.Exp(root, &exponent) - } - - // 2.3. shift the Miller loop result so that millerLoop * scalingFactor - // is of order finalExpFactor - scalingFactor.Mul(&rootPthInverse, &root27thInverse) - millerLoop.Mul(&millerLoop, &scalingFactor) - - // 3. get the witness residue - // - // lambda = q - u, the optimal exponent - var lambda big.Int - lambda.SetString("4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129030796414117214202539", 10) - exponent.ModInverse(&lambda, &finalExpFactor) - residueWitness.Exp(millerLoop, &exponent) + residueWitness, scalingFactor := finalExpWitness(&millerLoop) // return the witness residue residueWitness.C0.B0.A0.BigInt(outputs[0]) @@ -126,3 +67,136 @@ func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) er return nil }) } + +func pairingCheckHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { + // This is inspired from https://eprint.iacr.org/2024/640.pdf + // and based on a personal communication with the author Andrija Novakovic. + return emulated.UnwrapHint(nativeInputs, nativeOutputs, + func(mod *big.Int, inputs, outputs []*big.Int) error { + var P bls12381.G1Affine + var Q bls12381.G2Affine + n := len(inputs) + p := make([]bls12381.G1Affine, 0, n/6) + q := make([]bls12381.G2Affine, 0, n/6) + for k := 0; k < n/6+1; k += 2 { + P.X.SetBigInt(inputs[k]) + P.Y.SetBigInt(inputs[k+1]) + p = append(p, P) + } + for k := n / 3; k < n/2+3; k += 4 { + Q.X.A0.SetBigInt(inputs[k]) + Q.X.A1.SetBigInt(inputs[k+1]) + Q.Y.A0.SetBigInt(inputs[k+2]) + Q.Y.A1.SetBigInt(inputs[k+3]) + q = append(q, Q) + } + + lines := make([][2][len(bls12381.LoopCounter) - 1]bls12381.LineEvaluationAff, 0, len(q)) + for _, qi := range q { + lines = append(lines, bls12381.PrecomputeLines(qi)) + } + millerLoop, err := bls12381.MillerLoopFixedQ(p, lines) + if err != nil { + return err + } + millerLoop.Conjugate(&millerLoop) + + residueWitnessInv, scalingFactor := finalExpWitness(&millerLoop) + residueWitnessInv.Inverse(&residueWitnessInv) + + // return the witness residue + residueWitnessInv.C0.B0.A0.BigInt(outputs[0]) + residueWitnessInv.C0.B0.A1.BigInt(outputs[1]) + residueWitnessInv.C0.B1.A0.BigInt(outputs[2]) + residueWitnessInv.C0.B1.A1.BigInt(outputs[3]) + residueWitnessInv.C0.B2.A0.BigInt(outputs[4]) + residueWitnessInv.C0.B2.A1.BigInt(outputs[5]) + residueWitnessInv.C1.B0.A0.BigInt(outputs[6]) + residueWitnessInv.C1.B0.A1.BigInt(outputs[7]) + residueWitnessInv.C1.B1.A0.BigInt(outputs[8]) + residueWitnessInv.C1.B1.A1.BigInt(outputs[9]) + residueWitnessInv.C1.B2.A0.BigInt(outputs[10]) + residueWitnessInv.C1.B2.A1.BigInt(outputs[11]) + + // return the scaling factor + scalingFactor.C0.B0.A0.BigInt(outputs[12]) + scalingFactor.C0.B0.A1.BigInt(outputs[13]) + scalingFactor.C0.B1.A0.BigInt(outputs[14]) + scalingFactor.C0.B1.A1.BigInt(outputs[15]) + scalingFactor.C0.B2.A0.BigInt(outputs[16]) + scalingFactor.C0.B2.A1.BigInt(outputs[17]) + + return nil + + }) + +} + +func finalExpWitness(millerLoop *bls12381.E12) (residueWitness, scalingFactor bls12381.E12) { + + var root, rootPthInverse, root27thInverse bls12381.E12 + var order3rd, order3rdPower, exponent, exponentInv, finalExpFactor, polyFactor big.Int + // polyFactor = (1-x)/3 + polyFactor.SetString("5044125407647214251", 10) + // finalExpFactor = ((q^12 - 1) / r) / (27 * polyFactor) + finalExpFactor.SetString("2366356426548243601069753987687709088104621721678962410379583120840019275952471579477684846670499039076873213559162845121989217658133790336552276567078487633052653005423051750848782286407340332979263075575489766963251914185767058009683318020965829271737924625612375201545022326908440428522712877494557944965298566001441468676802477524234094954960009227631543471415676620753242466901942121887152806837594306028649150255258504417829961387165043999299071444887652375514277477719817175923289019181393803729926249507024121957184340179467502106891835144220611408665090353102353194448552304429530104218473070114105759487413726485729058069746063140422361472585604626055492939586602274983146215294625774144156395553405525711143696689756441298365274341189385646499074862712688473936093315628166094221735056483459332831845007196600723053356837526749543765815988577005929923802636375670820616189737737304893769679803809426304143627363860243558537831172903494450556755190448279875942974830469855835666815454271389438587399739607656399812689280234103023464545891697941661992848552456326290792224091557256350095392859243101357349751064730561345062266850238821755009430903520645523345000326783803935359711318798844368754833295302563158150573540616830138810935344206231367357992991289265295323280", 10) + + // 1. get pth-root inverse + exponent.Mul(&finalExpFactor, big.NewInt(27)) + root.Exp(*millerLoop, &exponent) + if root.IsOne() { + rootPthInverse.SetOne() + } else { + exponentInv.ModInverse(&exponent, &polyFactor) + exponent.Neg(&exponentInv).Mod(&exponent, &polyFactor) + rootPthInverse.Exp(root, &exponent) + } + + // 2.1. get order of 3rd primitive root + var three big.Int + three.SetUint64(3) + exponent.Mul(&polyFactor, &finalExpFactor) + root.Exp(*millerLoop, &exponent) + if root.IsOne() { + order3rdPower.SetUint64(0) + } + root.Exp(root, &three) + if root.IsOne() { + order3rdPower.SetUint64(1) + } + root.Exp(root, &three) + if root.IsOne() { + order3rdPower.SetUint64(2) + } + root.Exp(root, &three) + if root.IsOne() { + order3rdPower.SetUint64(3) + } + + // 2.2. get 27th root inverse + if order3rdPower.Uint64() == 0 { + root27thInverse.SetOne() + } else { + order3rd.Exp(&three, &order3rdPower, nil) + exponent.Mul(&polyFactor, &finalExpFactor) + root.Exp(*millerLoop, &exponent) + exponentInv.ModInverse(&exponent, &order3rd) + exponent.Neg(&exponentInv).Mod(&exponent, &order3rd) + root27thInverse.Exp(root, &exponent) + } + + // 2.3. shift the Miller loop result so that millerLoop * scalingFactor + // is of order finalExpFactor + scalingFactor.Mul(&rootPthInverse, &root27thInverse) + millerLoop.Mul(millerLoop, &scalingFactor) + + // 3. get the witness residue + // + // lambda = q - u, the optimal exponent + var lambda big.Int + lambda.SetString("4002409555221667393417789825735904156556882819939007885332058136124031650490837864442687629129030796414117214202539", 10) + exponent.ModInverse(&lambda, &finalExpFactor) + residueWitness.Exp(*millerLoop, &exponent) + + return residueWitness, scalingFactor +} diff --git a/std/algebra/emulated/sw_bls12381/pairing.go b/std/algebra/emulated/sw_bls12381/pairing.go index b798bf978..f19f592ed 100644 --- a/std/algebra/emulated/sw_bls12381/pairing.go +++ b/std/algebra/emulated/sw_bls12381/pairing.go @@ -100,14 +100,84 @@ 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.MillerLoop(P, Q) + // check input size match + nP := len(P) + nQ := len(Q) + if nP == 0 || nP != nQ { + return errors.New("invalid inputs sizes") + } + // hint the non-residue witness + inputs := make([]*baseEl, 0, 2*nP+4*nQ) + for _, p := range P { + inputs = append(inputs, &p.X, &p.Y) + } + for _, q := range Q { + inputs = append(inputs, &q.P.X.A0, &q.P.X.A1, &q.P.Y.A0, &q.P.Y.A1) + } + hint, err := pr.curveF.NewHint(pairingCheckHint, 18, inputs...) if err != nil { - return err + // err is non-nil only for invalid number of inputs + panic(err) + } + + residueWitnessInv := pr.FromTower([12]*baseEl{hint[0], hint[1], hint[2], hint[3], hint[4], hint[5], hint[6], hint[7], hint[8], hint[9], hint[10], hint[11]}) + // constrain cubicNonResiduePower to be in Fp6 + // that is: a100=a101=a110=a111=a120=a121=0 + // or + // A0 = a000 - a001 + // A1 = 0 + // A2 = a010 - a011 + // A3 = 0 + // A4 = a020 - a021 + // A5 = 0 + // A6 = a001 + // A7 = 0 + // A8 = a011 + // A9 = 0 + // A10 = a021 + // A11 = 0 + scalingFactor := GTEl{ + A0: *pr.curveF.Sub(hint[12], hint[13]), + A1: *pr.curveF.Zero(), + A2: *pr.curveF.Sub(hint[14], hint[15]), + A3: *pr.curveF.Zero(), + A4: *pr.curveF.Sub(hint[16], hint[17]), + A5: *pr.curveF.Zero(), + A6: *hint[13], + A7: *pr.curveF.Zero(), + A8: *hint[15], + A9: *pr.curveF.Zero(), + A10: *hint[17], + A11: *pr.curveF.Zero(), + } + lines := make([]lineEvaluations, nQ) + for i := range Q { + if Q[i].Lines == nil { + Qlines := pr.computeLines(&Q[i].P) + Q[i].Lines = &Qlines + } + lines[i] = *Q[i].Lines } - pr.AssertFinalExponentiationIsOne(f) + res, err := pr.millerLoopLines(P, lines, residueWitnessInv, false) + if err != nil { + return fmt.Errorf("miller loop: %w", err) + } + res = pr.Ext12.Conjugate(res) + // Check that: MillerLoop(P,Q) * scalingFactor * residueWitnessInv^(p-x₀) == 1 + // where u=-0xd201000000010000 is the BLS12-381 seed, and residueWitness, + // scalingFactor from the hint. + // Note that res is already MillerLoop(P,Q) * residueWitnessInv^{-x₀} since + // we initialized the Miller loop accumulator with residueWitnessInv. + // So we only need to check that: + // res * scalingFactor * residueWitnessInv^p == 1 + res = pr.Ext12.Mul(res, &scalingFactor) + t0 := pr.Frobenius(residueWitnessInv) + res = pr.Ext12.Mul(res, t0) + + pr.AssertIsEqual(res, pr.Ext12.One()) return nil } @@ -199,12 +269,12 @@ func (pr Pairing) MillerLoop(P []*G1Affine, Q []*G2Affine) (*GTEl, error) { } lines[i] = *Q[i].Lines } - return pr.millerLoopLines(P, lines) + return pr.millerLoopLines(P, lines, nil, true) } // millerLoopLines computes the multi-Miller loop from points in G1 and precomputed lines in G2 -func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl, error) { +func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations, init *GTEl, first bool) (*GTEl, error) { // check input size match n := len(P) @@ -225,34 +295,48 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl xNegOverY[k] = pr.curveF.Neg(xNegOverY[k]) } + // Compute ∏ᵢ { fᵢ_{x₀,Q}(P) } res := pr.Ext12.One() - // Compute ∏ᵢ { fᵢ_{x₀,Q}(P) } - // i = 62, separately to avoid an E12 Square - // (Square(res) = 1² = 1) - for k := 0; k < n; k++ { - res = pr.MulBy02368(res, - pr.MulByElement(&lines[k][0][62].R1, yInv[k]), - pr.MulByElement(&lines[k][0][62].R0, xNegOverY[k]), - ) - res = pr.MulBy02368(res, - pr.MulByElement(&lines[k][1][62].R1, yInv[k]), - pr.MulByElement(&lines[k][1][62].R0, xNegOverY[k]), - ) + if init != nil { + res = init + } + + j := len(loopCounter) - 2 + if first { + // i = j, separately to avoid an E12 Square + // (Square(res) = 1² = 1) + for k := 0; k < n; k++ { + res = pr.MulBy02368(res, + pr.MulByElement(&lines[k][0][j].R1, yInv[k]), + pr.MulByElement(&lines[k][0][j].R0, xNegOverY[k]), + ) + res = pr.MulBy02368(res, + pr.MulByElement(&lines[k][1][j].R1, yInv[k]), + pr.MulByElement(&lines[k][1][j].R0, xNegOverY[k]), + ) + } + j-- } - for i := 61; i >= 0; i-- { + for i := j; i >= 0; i-- { // mutualize the square among n Miller loops // (∏ᵢfᵢ)² res = pr.Ext12.Square(res) - for k := 0; k < n; k++ { - if loopCounter[i] == 0 { + if loopCounter[i] == 0 { + for k := 0; k < n; k++ { res = pr.MulBy02368(res, pr.MulByElement(&lines[k][0][i].R1, yInv[k]), pr.MulByElement(&lines[k][0][i].R0, xNegOverY[k]), ) - } else { + } + } else { + if init != nil { + // multiply by init when bit=1 + res = pr.Ext12.Mul(res, init) + } + for k := 0; k < n; k++ { res = pr.MulBy02368(res, pr.MulByElement(&lines[k][0][i].R1, yInv[k]), pr.MulByElement(&lines[k][0][i].R0, xNegOverY[k]), @@ -320,7 +404,7 @@ func (pr Pairing) FinalExponentiation(e *GTEl) *GTEl { func (pr Pairing) AssertFinalExponentiationIsOne(x *GTEl) { tower := pr.ToTower(x) - res, err := pr.curveF.NewHint(finalExpHint, 24, tower[0], tower[1], tower[2], tower[3], tower[4], tower[5], tower[6], tower[7], tower[8], tower[9], tower[10], tower[11]) + res, err := pr.curveF.NewHint(finalExpHint, 18, tower[0], tower[1], tower[2], tower[3], tower[4], tower[5], tower[6], tower[7], tower[8], tower[9], tower[10], tower[11]) if err != nil { // err is non-nil only for invalid number of inputs panic(err) diff --git a/std/algebra/emulated/sw_bls12381/pairing_test.go b/std/algebra/emulated/sw_bls12381/pairing_test.go index 30cbb14a0..7020c3240 100644 --- a/std/algebra/emulated/sw_bls12381/pairing_test.go +++ b/std/algebra/emulated/sw_bls12381/pairing_test.go @@ -33,6 +33,47 @@ func randomG1G2Affines() (bls12381.G1Affine, bls12381.G2Affine) { return p, q } +type MillerLoopCircuit struct { + In1G1, In2G1 G1Affine + In1G2, In2G2 G2Affine + Res GTEl +} + +func (c *MillerLoopCircuit) Define(api frontend.API) error { + pairing, err := NewPairing(api) + if err != nil { + return fmt.Errorf("new pairing: %w", err) + } + res, err := pairing.MillerLoop([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2}) + if err != nil { + return fmt.Errorf("pair: %w", err) + } + pairing.AssertIsEqual(res, &c.Res) + return nil +} + +func TestMillerLoopTestSolve(t *testing.T) { + assert := test.NewAssert(t) + p1, q1 := randomG1G2Affines() + p2, q2 := randomG1G2Affines() + lines1 := bls12381.PrecomputeLines(q1) + lines2 := bls12381.PrecomputeLines(q2) + res, err := bls12381.MillerLoopFixedQ( + []bls12381.G1Affine{p1, p2}, + [][2][len(bls12381.LoopCounter) - 1]bls12381.LineEvaluationAff{lines1, lines2}, + ) + assert.NoError(err) + witness := MillerLoopCircuit{ + In1G1: NewG1Affine(p1), + In1G2: NewG2Affine(q1), + In2G1: NewG1Affine(p2), + In2G2: NewG2Affine(q2), + Res: NewGTEl(res), + } + err = test.IsSolved(&MillerLoopCircuit{}, &witness, ecc.BN254.ScalarField()) + assert.NoError(err) +} + type FinalExponentiationCircuit struct { InGt GTEl Res GTEl diff --git a/std/algebra/emulated/sw_bn254/hints.go b/std/algebra/emulated/sw_bn254/hints.go index c73329d82..9c0bd4aaf 100644 --- a/std/algebra/emulated/sw_bn254/hints.go +++ b/std/algebra/emulated/sw_bn254/hints.go @@ -16,19 +16,126 @@ func init() { // GetHints returns all hint functions used in the package. func GetHints() []solver.Hint { return []solver.Hint{ - millerLoopAndCheckFinalExpHint, finalExpHint, + pairingCheckHint, + millerLoopAndCheckFinalExpHint, } } +func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { + // This follows section 4.3.2 of https://eprint.iacr.org/2024/640.pdf + return emulated.UnwrapHint(nativeInputs, nativeOutputs, + func(mod *big.Int, inputs, outputs []*big.Int) error { + + var millerLoop bn254.E12 + + 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]) + + residueWitness, cubicNonResiduePower := finalExpWitness(&millerLoop) + + 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]) + + // we also need to return the cubic non-residue power + cubicNonResiduePower.C0.B0.A0.BigInt(outputs[12]) + cubicNonResiduePower.C0.B0.A1.BigInt(outputs[13]) + cubicNonResiduePower.C0.B1.A0.BigInt(outputs[14]) + cubicNonResiduePower.C0.B1.A1.BigInt(outputs[15]) + cubicNonResiduePower.C0.B2.A0.BigInt(outputs[16]) + cubicNonResiduePower.C0.B2.A1.BigInt(outputs[17]) + + return nil + }) +} + +func pairingCheckHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { + // This follows section 4.3.2 of https://eprint.iacr.org/2024/640.pdf + return emulated.UnwrapHint(nativeInputs, nativeOutputs, + func(mod *big.Int, inputs, outputs []*big.Int) error { + var P bn254.G1Affine + var Q bn254.G2Affine + n := len(inputs) + p := make([]bn254.G1Affine, 0, n/6) + q := make([]bn254.G2Affine, 0, n/6) + for k := 0; k < n/6+1; k += 2 { + P.X.SetBigInt(inputs[k]) + P.Y.SetBigInt(inputs[k+1]) + p = append(p, P) + } + for k := n / 3; k < n/2+3; k += 4 { + Q.X.A0.SetBigInt(inputs[k]) + Q.X.A1.SetBigInt(inputs[k+1]) + Q.Y.A0.SetBigInt(inputs[k+2]) + Q.Y.A1.SetBigInt(inputs[k+3]) + q = append(q, Q) + } + + lines := make([][2][len(bn254.LoopCounter)]bn254.LineEvaluationAff, 0, len(q)) + for _, qi := range q { + lines = append(lines, bn254.PrecomputeLines(qi)) + } + millerLoop, err := bn254.MillerLoopFixedQ(p, lines) + if err != nil { + return err + } + + residueWitnessInv, cubicNonResiduePower := finalExpWitness(&millerLoop) + residueWitnessInv.Inverse(&residueWitnessInv) + + residueWitnessInv.C0.B0.A0.BigInt(outputs[0]) + residueWitnessInv.C0.B0.A1.BigInt(outputs[1]) + residueWitnessInv.C0.B1.A0.BigInt(outputs[2]) + residueWitnessInv.C0.B1.A1.BigInt(outputs[3]) + residueWitnessInv.C0.B2.A0.BigInt(outputs[4]) + residueWitnessInv.C0.B2.A1.BigInt(outputs[5]) + residueWitnessInv.C1.B0.A0.BigInt(outputs[6]) + residueWitnessInv.C1.B0.A1.BigInt(outputs[7]) + residueWitnessInv.C1.B1.A0.BigInt(outputs[8]) + residueWitnessInv.C1.B1.A1.BigInt(outputs[9]) + residueWitnessInv.C1.B2.A0.BigInt(outputs[10]) + residueWitnessInv.C1.B2.A1.BigInt(outputs[11]) + + // we also need to return the cubic non-residue power + cubicNonResiduePower.C0.B0.A0.BigInt(outputs[12]) + cubicNonResiduePower.C0.B0.A1.BigInt(outputs[13]) + cubicNonResiduePower.C0.B1.A0.BigInt(outputs[14]) + cubicNonResiduePower.C0.B1.A1.BigInt(outputs[15]) + cubicNonResiduePower.C0.B2.A0.BigInt(outputs[16]) + cubicNonResiduePower.C0.B2.A1.BigInt(outputs[17]) + + return nil + }) + +} + func millerLoopAndCheckFinalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { // This follows section 4.3.2 of https://eprint.iacr.org/2024/640.pdf return emulated.UnwrapHint(nativeInputs, nativeOutputs, func(mod *big.Int, inputs, outputs []*big.Int) error { - var previous, tmp, x3, cubicNonResiduePower, x, millerLoop, residueWitness, residueWitnessInv, one, root27thOf1 bn254.E12 - var exp1, exp2, rInv, mInv big.Int var P bn254.G1Affine var Q bn254.G2Affine + var previous bn254.E12 P.X.SetBigInt(inputs[0]) P.Y.SetBigInt(inputs[1]) @@ -65,78 +172,7 @@ func millerLoopAndCheckFinalExpHint(nativeMod *big.Int, nativeInputs, nativeOutp millerLoop.Mul(&millerLoop, &previous) - // exp1 = (p^12-1)/3 - exp1.SetString("4030969696062745741797811005853058291874379204406359442560681893891674450106959530046539719647151210908190211459382793062006703141168852426020468083171325367934590379984666859998399967609544754664110191464072930598755441160008826659219834762354786403012110463250131961575955268597858015384895449311534622125256548620283853223733396368939858981844663598065852816056384933498610930035891058807598891752166582271931875150099691598048016175399382213304673796601585080509443902692818733420199004555566113537482054218823936116647313678747500267068559627206777530424029211671772692598157901876223857571299238046741502089890557442500582300718504160740314926185458079985126192563953772118929726791041828902047546977272656240744693339962973939047279285351052107950250121751682659529260304162131862468322644288196213423232132152125277136333208005221619443705106431645884840489295409272576227859206166894626854018093044908314720", 10) - // root27thOf1 = (0, c010, c011, 0, 0, 0, 0, 0, 0, 0, 0, 0) - // is a 27-th root of unity which is necessarily a cubic non-residue - // since h/r = (p^12-1)/r = 27·l and 3 does not divide l. - // it was computed as w^((p^12-1)/27) = c2 * w^2 + c8 * w^8 where - // Fp12 = Fp[w]/w^12-18w^6+82 which is isomorphic to our Fp12 tower - // then c010 = (c2 + 9 * c8) % p and c011 = c8 - root27thOf1.C0.B1.A0.SetString("9483667112135124394372960210728142145589475128897916459350428495526310884707") - root27thOf1.C0.B1.A1.SetString("4534159768373982659291990808346042891252278737770656686799127720849666919525") - - if one.Exp(millerLoop, &exp1).IsOne() { - // residueWitness = millerLoop is a cubic residue - cubicNonResiduePower.SetOne() - residueWitness.Set(&millerLoop) - } else if one.Exp(*millerLoop.Mul(&millerLoop, &root27thOf1), &exp1).IsOne() { - // residueWitness = millerLoop * root27thOf1 is a cubic residue - cubicNonResiduePower.Set(&root27thOf1) - residueWitness.Set(&millerLoop) - } else { - // residueWitness = millerLoop * root27thOf1^2 is a cubic residue - cubicNonResiduePower.Square(&root27thOf1) - residueWitness.Mul(&millerLoop, &root27thOf1) - } - - // 1. compute r-th root: - // Exponentiate to rInv where - // rInv = 1/r mod (p^12-1)/r - rInv.SetString("495819184011867778744231927046742333492451180917315223017345540833046880485481720031136878341141903241966521818658471092566752321606779256340158678675679238405722886654128392203338228575623261160538734808887996935946888297414610216445334190959815200956855428635568184508263913274453942864817234480763055154719338281461936129150171789463489422401982681230261920147923652438266934726901346095892093443898852488218812468761027620988447655860644584419583586883569984588067403598284748297179498734419889699245081714359110559679136004228878808158639412436468707589339209058958785568729925402190575720856279605832146553573981587948304340677613460685405477047119496887534881410757668344088436651291444274840864486870663164657544390995506448087189408281061890434467956047582679858345583941396130713046072603335601764495918026585155498301896749919393", 10) - residueWitness.Exp(residueWitness, &rInv) - - // 2. compute m-th root: - // where m = (6x + 2 + q^3 - q^2 + q)/(3r) - // Exponentiate to mInv where - // mInv = 1/m mod p^12-1 - mInv.SetString("17840267520054779749190587238017784600702972825655245554504342129614427201836516118803396948809179149954197175783449826546445899524065131269177708416982407215963288737761615699967145070776364294542559324079147363363059480104341231360692143673915822421222230661528586799190306058519400019024762424366780736540525310403098758015600523609594113357130678138304964034267260758692953579514899054295817541844330584721967571697039986079722203518034173581264955381924826388858518077894154909963532054519350571947910625755075099598588672669612434444513251495355121627496067454526862754597351094345783576387352673894873931328099247263766690688395096280633426669535619271711975898132416216382905928886703963310231865346128293216316379527200971959980873989485521004596686352787540034457467115536116148612884807380187255514888720048664139404687086409399", 10) - residueWitness.Exp(residueWitness, &mInv) - - // 3. compute cube root: - // since gcd(3, (p^12-1)/r) ≠ 1 we use a modified Toneelli-Shanks algorithm - // see Alg.4 of https://eprint.iacr.org/2024/640.pdf - // Typo in the paper: p^k-1 = 3^n * s instead of p-1 = 3^r * s - // where k=12 and n=3 here and exp2 = (s+1)/3 - residueWitnessInv.Inverse(&residueWitness) - exp2.SetString("149295173928249842288807815031594751550902933496531831205951181255247201855813315927649619246190785589192230054051214557852100116339587126889646966043382421034614458517950624444385183985538694617189266350521219651805757080000326913304438324531658755667115202342597480058368713651772519088329461085612393412046538837788290860138273939590365147475728281409846400594680923462911515927255224400281440435265428973034513894448136725853630228718495637529802733207466114092942366766400693830377740909465411612499335341437923559875826432546203713595131838044695464089778859691547136762894737106526809539677749557286722299625576201574095640767352005953344997266128077036486155280146436004404804695964512181557316554713802082990544197776406442186936269827816744738898152657469728130713344598597476387715653492155415311971560450078713968012341037230430349766855793764662401499603533676762082513303932107208402000670112774382027", 10) - x.Exp(residueWitness, &exp2) - - // 3^t is ord(x^3 / residueWitness) - x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv) - t := 0 - for !x3.IsOne() { - t++ - tmp.Square(&x3) - x3.Mul(&tmp, &x3) - } - - for t != 0 { - x.Mul(&x, tmp.Exp(root27thOf1, &exp2)) - - // 3^t is ord(x^3 / residueWitness) - x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv) - t = 0 - for !x3.IsOne() { - t++ - tmp.Square(&x3) - x3.Mul(&tmp, &x3) - } - } - - // x is now the cube root of residueWitness - residueWitness.Set(&x) - residueWitnessInv.Inverse(&residueWitness) + residueWitness, cubicNonResiduePower := finalExpWitness(&millerLoop) residueWitness.C0.B0.A0.BigInt(outputs[0]) residueWitness.C0.B0.A1.BigInt(outputs[1]) @@ -163,125 +199,80 @@ func millerLoopAndCheckFinalExpHint(nativeMod *big.Int, nativeInputs, nativeOutp }) } -func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { - // This follows section 4.3.2 of https://eprint.iacr.org/2024/640.pdf - return emulated.UnwrapHint(nativeInputs, nativeOutputs, - func(mod *big.Int, inputs, outputs []*big.Int) error { - var tmp, x3, cubicNonResiduePower, x, millerLoop, residueWitness, residueWitnessInv, one, root27thOf1 bn254.E12 - var exp1, exp2, rInv, mInv 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]) - - // exp1 = (p^12-1)/3 - exp1.SetString("4030969696062745741797811005853058291874379204406359442560681893891674450106959530046539719647151210908190211459382793062006703141168852426020468083171325367934590379984666859998399967609544754664110191464072930598755441160008826659219834762354786403012110463250131961575955268597858015384895449311534622125256548620283853223733396368939858981844663598065852816056384933498610930035891058807598891752166582271931875150099691598048016175399382213304673796601585080509443902692818733420199004555566113537482054218823936116647313678747500267068559627206777530424029211671772692598157901876223857571299238046741502089890557442500582300718504160740314926185458079985126192563953772118929726791041828902047546977272656240744693339962973939047279285351052107950250121751682659529260304162131862468322644288196213423232132152125277136333208005221619443705106431645884840489295409272576227859206166894626854018093044908314720", 10) - // root27thOf1 = (0, c010, c011, 0, 0, 0, 0, 0, 0, 0, 0, 0) - // is a 27-th root of unity which is necessarily a cubic non-residue - // since h/r = (p^12-1)/r = 27·l and 3 does not divide l. - // it was computed as w^((p^12-1)/27) = c2 * w^2 + c8 * w^8 where - // Fp12 = Fp[w]/w^12-18w^6+82 which is isomorphic to our Fp12 tower - // then c010 = (c2 + 9 * c8) % p and c011 = c8 - root27thOf1.C0.B1.A0.SetString("9483667112135124394372960210728142145589475128897916459350428495526310884707") - root27thOf1.C0.B1.A1.SetString("4534159768373982659291990808346042891252278737770656686799127720849666919525") - - if one.Exp(millerLoop, &exp1).IsOne() { - // residueWitness = millerLoop is a cubic residue - cubicNonResiduePower.SetOne() - residueWitness.Set(&millerLoop) - } else if one.Exp(*millerLoop.Mul(&millerLoop, &root27thOf1), &exp1).IsOne() { - // residueWitness = millerLoop * root27thOf1 is a cubic residue - cubicNonResiduePower.Set(&root27thOf1) - residueWitness.Set(&millerLoop) - } else { - // residueWitness = millerLoop * root27thOf1^2 is a cubic residue - cubicNonResiduePower.Square(&root27thOf1) - residueWitness.Mul(&millerLoop, &root27thOf1) - } - - // 1. compute r-th root: - // Exponentiate to rInv where - // rInv = 1/r mod (p^12-1)/r - rInv.SetString("495819184011867778744231927046742333492451180917315223017345540833046880485481720031136878341141903241966521818658471092566752321606779256340158678675679238405722886654128392203338228575623261160538734808887996935946888297414610216445334190959815200956855428635568184508263913274453942864817234480763055154719338281461936129150171789463489422401982681230261920147923652438266934726901346095892093443898852488218812468761027620988447655860644584419583586883569984588067403598284748297179498734419889699245081714359110559679136004228878808158639412436468707589339209058958785568729925402190575720856279605832146553573981587948304340677613460685405477047119496887534881410757668344088436651291444274840864486870663164657544390995506448087189408281061890434467956047582679858345583941396130713046072603335601764495918026585155498301896749919393", 10) - residueWitness.Exp(residueWitness, &rInv) - - // 2. compute m-th root: - // where m = (6x + 2 + q^3 - q^2 + q)/(3r) - // Exponentiate to mInv where - // mInv = 1/m mod p^12-1 - mInv.SetString("17840267520054779749190587238017784600702972825655245554504342129614427201836516118803396948809179149954197175783449826546445899524065131269177708416982407215963288737761615699967145070776364294542559324079147363363059480104341231360692143673915822421222230661528586799190306058519400019024762424366780736540525310403098758015600523609594113357130678138304964034267260758692953579514899054295817541844330584721967571697039986079722203518034173581264955381924826388858518077894154909963532054519350571947910625755075099598588672669612434444513251495355121627496067454526862754597351094345783576387352673894873931328099247263766690688395096280633426669535619271711975898132416216382905928886703963310231865346128293216316379527200971959980873989485521004596686352787540034457467115536116148612884807380187255514888720048664139404687086409399", 10) - residueWitness.Exp(residueWitness, &mInv) - - // 3. compute cube root: - // since gcd(3, (p^12-1)/r) ≠ 1 we use a modified Toneelli-Shanks algorithm - // see Alg.4 of https://eprint.iacr.org/2024/640.pdf - // Typo in the paper: p^k-1 = 3^n * s instead of p-1 = 3^r * s - // where k=12 and n=3 here and exp2 = (s+1)/3 - residueWitnessInv.Inverse(&residueWitness) - exp2.SetString("149295173928249842288807815031594751550902933496531831205951181255247201855813315927649619246190785589192230054051214557852100116339587126889646966043382421034614458517950624444385183985538694617189266350521219651805757080000326913304438324531658755667115202342597480058368713651772519088329461085612393412046538837788290860138273939590365147475728281409846400594680923462911515927255224400281440435265428973034513894448136725853630228718495637529802733207466114092942366766400693830377740909465411612499335341437923559875826432546203713595131838044695464089778859691547136762894737106526809539677749557286722299625576201574095640767352005953344997266128077036486155280146436004404804695964512181557316554713802082990544197776406442186936269827816744738898152657469728130713344598597476387715653492155415311971560450078713968012341037230430349766855793764662401499603533676762082513303932107208402000670112774382027", 10) - x.Exp(residueWitness, &exp2) - - // 3^t is ord(x^3 / residueWitness) - x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv) - t := 0 - for !x3.IsOne() { - t++ - tmp.Square(&x3) - x3.Mul(&tmp, &x3) - } - - for t != 0 { - x.Mul(&x, tmp.Exp(root27thOf1, &exp2)) - - // 3^t is ord(x^3 / residueWitness) - x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv) - t = 0 - for !x3.IsOne() { - t++ - tmp.Square(&x3) - x3.Mul(&tmp, &x3) - } - } +func finalExpWitness(millerLoop *bn254.E12) (residueWitness, cubicNonResiduePower bn254.E12) { + var tmp, x3, x, residueWitnessInv, one, root27thOf1 bn254.E12 + var exp1, exp2, rInv, mInv big.Int + // exp1 = (p^12-1)/3 + exp1.SetString("4030969696062745741797811005853058291874379204406359442560681893891674450106959530046539719647151210908190211459382793062006703141168852426020468083171325367934590379984666859998399967609544754664110191464072930598755441160008826659219834762354786403012110463250131961575955268597858015384895449311534622125256548620283853223733396368939858981844663598065852816056384933498610930035891058807598891752166582271931875150099691598048016175399382213304673796601585080509443902692818733420199004555566113537482054218823936116647313678747500267068559627206777530424029211671772692598157901876223857571299238046741502089890557442500582300718504160740314926185458079985126192563953772118929726791041828902047546977272656240744693339962973939047279285351052107950250121751682659529260304162131862468322644288196213423232132152125277136333208005221619443705106431645884840489295409272576227859206166894626854018093044908314720", 10) + // root27thOf1 = (0, c010, c011, 0, 0, 0, 0, 0, 0, 0, 0, 0) + // is a 27-th root of unity which is necessarily a cubic non-residue + // since h/r = (p^12-1)/r = 27·l and 3 does not divide l. + // it was computed as w^((p^12-1)/27) = c2 * w^2 + c8 * w^8 where + // Fp12 = Fp[w]/w^12-18w^6+82 which is isomorphic to our Fp12 tower + // then c010 = (c2 + 9 * c8) % p and c011 = c8 + root27thOf1.C0.B1.A0.SetString("9483667112135124394372960210728142145589475128897916459350428495526310884707") + root27thOf1.C0.B1.A1.SetString("4534159768373982659291990808346042891252278737770656686799127720849666919525") + + if one.Exp(*millerLoop, &exp1).IsOne() { + // residueWitness = millerLoop is a cubic residue + cubicNonResiduePower.SetOne() + residueWitness.Set(millerLoop) + } else if one.Exp(*millerLoop.Mul(millerLoop, &root27thOf1), &exp1).IsOne() { + // residueWitness = millerLoop * root27thOf1 is a cubic residue + cubicNonResiduePower.Set(&root27thOf1) + residueWitness.Set(millerLoop) + } else { + // residueWitness = millerLoop * root27thOf1^2 is a cubic residue + cubicNonResiduePower.Square(&root27thOf1) + residueWitness.Mul(millerLoop, &root27thOf1) + } - // x is now the cube root of residueWitness - residueWitness.Set(&x) + // 1. compute r-th root: + // Exponentiate to rInv where + // rInv = 1/r mod (p^12-1)/r + rInv.SetString("495819184011867778744231927046742333492451180917315223017345540833046880485481720031136878341141903241966521818658471092566752321606779256340158678675679238405722886654128392203338228575623261160538734808887996935946888297414610216445334190959815200956855428635568184508263913274453942864817234480763055154719338281461936129150171789463489422401982681230261920147923652438266934726901346095892093443898852488218812468761027620988447655860644584419583586883569984588067403598284748297179498734419889699245081714359110559679136004228878808158639412436468707589339209058958785568729925402190575720856279605832146553573981587948304340677613460685405477047119496887534881410757668344088436651291444274840864486870663164657544390995506448087189408281061890434467956047582679858345583941396130713046072603335601764495918026585155498301896749919393", 10) + residueWitness.Exp(residueWitness, &rInv) + + // 2. compute m-th root: + // where m = (6x + 2 + q^3 - q^2 + q)/(3r) + // Exponentiate to mInv where + // mInv = 1/m mod p^12-1 + mInv.SetString("17840267520054779749190587238017784600702972825655245554504342129614427201836516118803396948809179149954197175783449826546445899524065131269177708416982407215963288737761615699967145070776364294542559324079147363363059480104341231360692143673915822421222230661528586799190306058519400019024762424366780736540525310403098758015600523609594113357130678138304964034267260758692953579514899054295817541844330584721967571697039986079722203518034173581264955381924826388858518077894154909963532054519350571947910625755075099598588672669612434444513251495355121627496067454526862754597351094345783576387352673894873931328099247263766690688395096280633426669535619271711975898132416216382905928886703963310231865346128293216316379527200971959980873989485521004596686352787540034457467115536116148612884807380187255514888720048664139404687086409399", 10) + residueWitness.Exp(residueWitness, &mInv) + + // 3. compute cube root: + // since gcd(3, (p^12-1)/r) ≠ 1 we use a modified Toneelli-Shanks algorithm + // see Alg.4 of https://eprint.iacr.org/2024/640.pdf + // Typo in the paper: p^k-1 = 3^n * s instead of p-1 = 3^r * s + // where k=12 and n=3 here and exp2 = (s+1)/3 + residueWitnessInv.Inverse(&residueWitness) + exp2.SetString("149295173928249842288807815031594751550902933496531831205951181255247201855813315927649619246190785589192230054051214557852100116339587126889646966043382421034614458517950624444385183985538694617189266350521219651805757080000326913304438324531658755667115202342597480058368713651772519088329461085612393412046538837788290860138273939590365147475728281409846400594680923462911515927255224400281440435265428973034513894448136725853630228718495637529802733207466114092942366766400693830377740909465411612499335341437923559875826432546203713595131838044695464089778859691547136762894737106526809539677749557286722299625576201574095640767352005953344997266128077036486155280146436004404804695964512181557316554713802082990544197776406442186936269827816744738898152657469728130713344598597476387715653492155415311971560450078713968012341037230430349766855793764662401499603533676762082513303932107208402000670112774382027", 10) + x.Exp(residueWitness, &exp2) + + // 3^t is ord(x^3 / residueWitness) + x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv) + t := 0 + for !x3.IsOne() { + t++ + tmp.Square(&x3) + x3.Mul(&tmp, &x3) + } - 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]) + for t != 0 { + x.Mul(&x, tmp.Exp(root27thOf1, &exp2)) + + // 3^t is ord(x^3 / residueWitness) + x3.Square(&x).Mul(&x3, &x).Mul(&x3, &residueWitnessInv) + t = 0 + for !x3.IsOne() { + t++ + tmp.Square(&x3) + x3.Mul(&tmp, &x3) + } + } - // we also need to return the cubic non-residue power - cubicNonResiduePower.C0.B0.A0.BigInt(outputs[12]) - cubicNonResiduePower.C0.B0.A1.BigInt(outputs[13]) - cubicNonResiduePower.C0.B1.A0.BigInt(outputs[14]) - cubicNonResiduePower.C0.B1.A1.BigInt(outputs[15]) - cubicNonResiduePower.C0.B2.A0.BigInt(outputs[16]) - cubicNonResiduePower.C0.B2.A1.BigInt(outputs[17]) - cubicNonResiduePower.C1.B0.A0.BigInt(outputs[18]) - cubicNonResiduePower.C1.B0.A1.BigInt(outputs[19]) - cubicNonResiduePower.C1.B1.A0.BigInt(outputs[20]) - cubicNonResiduePower.C1.B1.A1.BigInt(outputs[21]) - cubicNonResiduePower.C1.B2.A0.BigInt(outputs[22]) - cubicNonResiduePower.C1.B2.A1.BigInt(outputs[23]) + // x is now the cube root of residueWitness + residueWitness.Set(&x) - return nil - }) + return residueWitness, cubicNonResiduePower } diff --git a/std/algebra/emulated/sw_bn254/pairing.go b/std/algebra/emulated/sw_bn254/pairing.go index 4172de432..0e38ae553 100644 --- a/std/algebra/emulated/sw_bn254/pairing.go +++ b/std/algebra/emulated/sw_bn254/pairing.go @@ -148,7 +148,7 @@ func (pr Pairing) FinalExponentiation(e *GTEl) *GTEl { func (pr Pairing) AssertFinalExponentiationIsOne(a *GTEl) { tower := pr.Ext12.ToTower(a) - res, err := pr.curveF.NewHint(finalExpHint, 24, tower[0], tower[1], tower[2], tower[3], tower[4], tower[5], tower[6], tower[7], tower[8], tower[9], tower[10], tower[11]) + res, err := pr.curveF.NewHint(finalExpHint, 18, tower[0], tower[1], tower[2], tower[3], tower[4], tower[5], tower[6], tower[7], tower[8], tower[9], tower[10], tower[11]) if err != nil { // err is non-nil only for invalid number of inputs panic(err) @@ -214,13 +214,89 @@ func (pr Pairing) AssertFinalExponentiationIsOne(a *GTEl) { // This function checks that the Qᵢ are in the correct subgroup, but does not // check Pᵢ. See AssertIsOnG1. func (pr Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error { - f, err := pr.MillerLoop(P, Q) + // check input size match + nP := len(P) + nQ := len(Q) + if nP == 0 || nP != nQ { + return errors.New("invalid inputs sizes") + } + + // hint the non-residue witness + inputs := make([]*baseEl, 0, 2*nP+4*nQ) + for _, p := range P { + inputs = append(inputs, &p.X, &p.Y) + } + for _, q := range Q { + inputs = append(inputs, &q.P.X.A0, &q.P.X.A1, &q.P.Y.A0, &q.P.Y.A1) + } + hint, err := pr.curveF.NewHint(pairingCheckHint, 18, inputs...) if err != nil { - return err + // err is non-nil only for invalid number of inputs + panic(err) + } + residueWitnessInv := pr.Ext12.FromTower([12]*baseEl{hint[0], hint[1], hint[2], hint[3], hint[4], hint[5], hint[6], hint[7], hint[8], hint[9], hint[10], hint[11]}) + + // constrain cubicNonResiduePower to be in Fp6 + // that is: a100=a101=a110=a111=a120=a121=0 + // or + // A0 = a000 - 9 * a001 + // A1 = 0 + // A2 = a010 - 9 * a011 + // A3 = 0 + // A4 = a020 - 9 * a021 + // A5 = 0 + // A6 = a001 + // A7 = 0 + // A8 = a011 + // A9 = 0 + // A10 = a021 + // A11 = 0 + nine := big.NewInt(9) + cubicNonResiduePower := GTEl{ + A0: *pr.curveF.Sub(hint[12], pr.curveF.MulConst(hint[13], nine)), + A1: *pr.curveF.Zero(), + A2: *pr.curveF.Sub(hint[14], pr.curveF.MulConst(hint[15], nine)), + A3: *pr.curveF.Zero(), + A4: *pr.curveF.Sub(hint[16], pr.curveF.MulConst(hint[17], nine)), + A5: *pr.curveF.Zero(), + A6: *hint[13], + A7: *pr.curveF.Zero(), + A8: *hint[15], + A9: *pr.curveF.Zero(), + A10: *hint[17], + A11: *pr.curveF.Zero(), + } + + lines := make([]lineEvaluations, nQ) + for i := range Q { + if Q[i].Lines == nil { + Qlines := pr.computeLines(&Q[i].P) + Q[i].Lines = &Qlines + } + lines[i] = *Q[i].Lines + } + res, err := pr.millerLoopLines(P, lines, residueWitnessInv, false) + if err != nil { + return fmt.Errorf("miller loop: %w", err) } - pr.AssertFinalExponentiationIsOne(f) + // Check that res * cubicNonResiduePower * residueWitnessInv^λ' == 1 + // where λ' = q^3 - q^2 + q, with u the BN254 seed + // and residueWitnessInv, cubicNonResiduePower from the hint. + // Note that res is already MillerLoop(P,Q) * residueWitnessInv^{6x₀+2} since + // we initialized the Miller loop accumulator with residueWitnessInv. + t2 := pr.Ext12.Mul(&cubicNonResiduePower, res) + + t1 := pr.FrobeniusCube(residueWitnessInv) + t0 := pr.FrobeniusSquare(residueWitnessInv) + t1 = pr.Ext12.DivUnchecked(t1, t0) + t0 = pr.Frobenius(residueWitnessInv) + t1 = pr.Ext12.Mul(t1, t0) + + t2 = pr.Ext12.Mul(t2, t1) + + pr.AssertIsEqual(t2, pr.Ext12.One()) return nil } @@ -455,12 +531,12 @@ func (pr Pairing) MillerLoop(P []*G1Affine, Q []*G2Affine) (*GTEl, error) { } lines[i] = *Q[i].Lines } - return pr.millerLoopLines(P, lines) + return pr.millerLoopLines(P, lines, nil, true) } // millerLoopLines computes the multi-Miller loop from points in G1 and precomputed lines in G2 -func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl, error) { +func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations, init *GTEl, first bool) (*GTEl, error) { // check input size match n := len(P) @@ -481,82 +557,115 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl xNegOverY[k] = pr.curveF.Neg(xNegOverY[k]) } + // Compute f_{6x₀+2,Q}(P) var prodLines [10]*baseEl + res := pr.Ext12.One() - // Compute f_{6x₀+2,Q}(P) - // i = 64 - // - // k = 0 - c3 := pr.Ext2.MulByElement(&lines[0][0][64].R0, xNegOverY[0]) - c4 := pr.Ext2.MulByElement(&lines[0][0][64].R1, yInv[0]) - nine := big.NewInt(9) - res := >El{ - A0: *pr.curveF.One(), - A1: *pr.curveF.Sub(&c3.A0, pr.curveF.MulConst(&c3.A1, nine)), - A2: *pr.curveF.Zero(), - A3: *pr.curveF.Sub(&c4.A0, pr.curveF.MulConst(&c4.A1, nine)), - A4: *pr.curveF.Zero(), - A5: *pr.curveF.Zero(), - A6: *pr.curveF.Zero(), - A7: c3.A1, - A8: *pr.curveF.Zero(), - A9: c4.A1, - A10: *pr.curveF.Zero(), - A11: *pr.curveF.Zero(), + var initInv GTEl + if init != nil { + res = init + initInv = *pr.Ext12.Inverse(init) } - if n >= 2 { - // k = 1, separately to avoid MulBy01379 (res × ℓ) - // (res is also a line at this point, so we use Mul01379By01379 ℓ × ℓ) - // line evaluation at P[1] - prodLines = pr.Mul01379By01379( - pr.Ext2.MulByElement(&lines[1][0][64].R0, xNegOverY[1]), - pr.Ext2.MulByElement(&lines[1][0][64].R1, yInv[1]), - c3, - c4, - ) + j := len(loopCounter) - 2 + if first { + // i = j + // k = 0 + c3 := pr.Ext2.MulByElement(&lines[0][0][j].R0, xNegOverY[0]) + c4 := pr.Ext2.MulByElement(&lines[0][0][j].R1, yInv[0]) + nine := big.NewInt(9) res = >El{ - A0: *prodLines[0], - A1: *prodLines[1], - A2: *prodLines[2], - A3: *prodLines[3], - A4: *prodLines[4], + A0: *pr.curveF.One(), + A1: *pr.curveF.Sub(&c3.A0, pr.curveF.MulConst(&c3.A1, nine)), + A2: *pr.curveF.Zero(), + A3: *pr.curveF.Sub(&c4.A0, pr.curveF.MulConst(&c4.A1, nine)), + A4: *pr.curveF.Zero(), A5: *pr.curveF.Zero(), - A6: *prodLines[5], - A7: *prodLines[6], - A8: *prodLines[7], - A9: *prodLines[8], - A10: *prodLines[9], + A6: *pr.curveF.Zero(), + A7: c3.A1, + A8: *pr.curveF.Zero(), + A9: c4.A1, + A10: *pr.curveF.Zero(), A11: *pr.curveF.Zero(), } - } - if n >= 3 { - // k >= 2 - for k := 2; k < n; k++ { - // line evaluation at P[k] - // ℓ × res - res = pr.MulBy01379( - res, - pr.Ext2.MulByElement(&lines[k][0][64].R0, xNegOverY[k]), - pr.Ext2.MulByElement(&lines[k][0][64].R1, yInv[k]), + if n >= 2 { + // k = 1, separately to avoid MulBy01379 (res × ℓ) + // (res is also a line at this point, so we use Mul01379By01379 ℓ × ℓ) + // line evaluation at P[1] + prodLines = pr.Mul01379By01379( + pr.Ext2.MulByElement(&lines[1][0][j].R0, xNegOverY[1]), //nolint: gosec // incorrectly flagged by gosec as out of bounds read (G602) + pr.Ext2.MulByElement(&lines[1][0][j].R1, yInv[1]), //nolint: gosec // incorrectly flagged by gosec as out of bounds read (G602) + c3, + c4, ) + res = >El{ + A0: *prodLines[0], + A1: *prodLines[1], + A2: *prodLines[2], + A3: *prodLines[3], + A4: *prodLines[4], + A5: *pr.curveF.Zero(), + A6: *prodLines[5], + A7: *prodLines[6], + A8: *prodLines[7], + A9: *prodLines[8], + A10: *prodLines[9], + A11: *pr.curveF.Zero(), + } + } + + if n >= 3 { + // k >= 2 + for k := 2; k < n; k++ { + // line evaluation at P[k] + // ℓ × res + res = pr.MulBy01379( + res, + pr.Ext2.MulByElement(&lines[k][0][j].R0, xNegOverY[k]), + pr.Ext2.MulByElement(&lines[k][0][j].R1, yInv[k]), + ) + } } + j-- } - for i := 63; i >= 0; i-- { + for i := j; i >= 0; i-- { res = pr.Ext12.Square(res) - for k := 0; k < n; k++ { - if loopCounter[i] == 0 { + switch loopCounter[i] { + case 0: + for k := 0; k < n; k++ { res = pr.MulBy01379( res, pr.Ext2.MulByElement(&lines[k][0][i].R0, xNegOverY[k]), pr.Ext2.MulByElement(&lines[k][0][i].R1, yInv[k]), ) - } else { + } + case 1: + if init != nil { + // multiply by init when bit=1 + res = pr.Ext12.Mul(res, init) + } + // ℓ × ℓ + for k := 0; k < n; k++ { + prodLines := pr.Mul01379By01379( + pr.Ext2.MulByElement(&lines[k][0][i].R0, xNegOverY[k]), + pr.Ext2.MulByElement(&lines[k][0][i].R1, yInv[k]), + pr.Ext2.MulByElement(&lines[k][1][i].R0, xNegOverY[k]), + pr.Ext2.MulByElement(&lines[k][1][i].R1, yInv[k]), + ) + // (ℓ × ℓ) × res + res = pr.Ext12.MulBy012346789(res, prodLines) + } + case -1: + if init != nil { + // multiply by 1/init when bit=-1 + res = pr.Ext12.Mul(res, &initInv) + } + for k := 0; k < n; k++ { // ℓ × ℓ - prodLines = pr.Mul01379By01379( + prodLines := pr.Mul01379By01379( pr.Ext2.MulByElement(&lines[k][0][i].R0, xNegOverY[k]), pr.Ext2.MulByElement(&lines[k][0][i].R1, yInv[k]), pr.Ext2.MulByElement(&lines[k][1][i].R0, xNegOverY[k]), @@ -565,6 +674,8 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl // (ℓ × ℓ) × res res = pr.Ext12.MulBy012346789(res, prodLines) } + default: + panic(fmt.Sprintf("invalid loop counter value %d", loopCounter[i])) } } @@ -758,7 +869,7 @@ func (pr Pairing) millerLoopAndFinalExpResult(P *G1Affine, Q *G2Affine, previous // err is non-nil only for invalid number of inputs panic(err) } - residueWitness := pr.FromTower([12]*baseEl{hint[0], hint[1], hint[2], hint[3], hint[4], hint[5], hint[6], hint[7], hint[8], hint[9], hint[10], hint[11]}) + residueWitness := pr.Ext12.FromTower([12]*baseEl{hint[0], hint[1], hint[2], hint[3], hint[4], hint[5], hint[6], hint[7], hint[8], hint[9], hint[10], hint[11]}) // constrain cubicNonResiduePower to be in Fp6 // that is: a100=a101=a110=a111=a120=a121=0 @@ -800,67 +911,15 @@ func (pr Pairing) millerLoopAndFinalExpResult(P *G1Affine, Q *G2Affine, previous } lines := *Q.Lines - // precomputations - yInv := pr.curveF.Inverse(&P.Y) - xNegOverY := pr.curveF.Mul(&P.X, yInv) - xNegOverY = pr.curveF.Neg(xNegOverY) - - // init Miller loop accumulator to residueWitnessInv to share the squarings - // of residueWitnessInv^{6x₀+2} - res := residueWitnessInv - - // Compute f_{6x₀+2,Q}(P) - for i := 64; i >= 0; i-- { - res = pr.Ext12.Square(res) - - switch loopCounter[i] { - case 0: - // ℓ × res - res = pr.MulBy01379( - res, - pr.Ext2.MulByElement(&lines[0][i].R0, xNegOverY), - pr.Ext2.MulByElement(&lines[0][i].R1, yInv), - ) - case 1: - // multiply by residueWitnessInv when bit=1 - res = pr.Ext12.Mul(res, residueWitnessInv) - // lines evaluations at P - // and ℓ × ℓ - prodLines := pr.Mul01379By01379( - pr.Ext2.MulByElement(&lines[0][i].R0, xNegOverY), - pr.Ext2.MulByElement(&lines[0][i].R1, yInv), - pr.Ext2.MulByElement(&lines[1][i].R0, xNegOverY), - pr.Ext2.MulByElement(&lines[1][i].R1, yInv), - ) - // (ℓ × ℓ) × res - res = pr.Ext12.MulBy012346789(res, prodLines) - case -1: - // multiply by residueWitness when bit=-1 - res = pr.Ext12.Mul(res, residueWitness) - // lines evaluations at P - // and ℓ × ℓ - prodLines := pr.Mul01379By01379( - pr.Ext2.MulByElement(&lines[0][i].R0, xNegOverY), - pr.Ext2.MulByElement(&lines[0][i].R1, yInv), - pr.Ext2.MulByElement(&lines[1][i].R0, xNegOverY), - pr.Ext2.MulByElement(&lines[1][i].R1, yInv), - ) - // (ℓ × ℓ) × res - res = pr.Ext12.MulBy012346789(res, prodLines) - default: - panic(fmt.Sprintf("invalid loop counter value %d", loopCounter[i])) - } - } - - // Compute ℓ_{[6x₀+2]Q,π(Q)}(P) · ℓ_{[6x₀+2]Q+π(Q),-π²(Q)}(P) - // lines evaluations at P - prodLines := pr.Mul01379By01379( - pr.Ext2.MulByElement(&lines[0][65].R0, xNegOverY), - pr.Ext2.MulByElement(&lines[0][65].R1, yInv), - pr.Ext2.MulByElement(&lines[1][65].R0, xNegOverY), - pr.Ext2.MulByElement(&lines[1][65].R1, yInv), + res, err := pr.millerLoopLines( + []*G1Affine{P}, + []lineEvaluations{lines}, + residueWitnessInv, + false, ) - res = pr.Ext12.MulBy012346789(res, prodLines) + if err != nil { + return nil + } // multiply by previous multi-Miller function res = pr.Ext12.Mul(res, previous) diff --git a/std/algebra/emulated/sw_bn254/pairing_test.go b/std/algebra/emulated/sw_bn254/pairing_test.go index ae2fe87e8..70e97c975 100644 --- a/std/algebra/emulated/sw_bn254/pairing_test.go +++ b/std/algebra/emulated/sw_bn254/pairing_test.go @@ -33,6 +33,47 @@ func randomG1G2Affines() (bn254.G1Affine, bn254.G2Affine) { return p, q } +type MillerLoopCircuit struct { + In1G1, In2G1 G1Affine + In1G2, In2G2 G2Affine + Res GTEl +} + +func (c *MillerLoopCircuit) Define(api frontend.API) error { + pairing, err := NewPairing(api) + if err != nil { + return fmt.Errorf("new pairing: %w", err) + } + res, err := pairing.MillerLoop([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2}) + if err != nil { + return fmt.Errorf("pair: %w", err) + } + pairing.AssertIsEqual(res, &c.Res) + return nil +} + +func TestMillerLoopTestSolve(t *testing.T) { + assert := test.NewAssert(t) + p1, q1 := randomG1G2Affines() + p2, q2 := randomG1G2Affines() + lines1 := bn254.PrecomputeLines(q1) + lines2 := bn254.PrecomputeLines(q2) + res, err := bn254.MillerLoopFixedQ( + []bn254.G1Affine{p1, p2}, + [][2][len(bn254.LoopCounter)]bn254.LineEvaluationAff{lines1, lines2}, + ) + assert.NoError(err) + witness := MillerLoopCircuit{ + In1G1: NewG1Affine(p1), + In1G2: NewG2Affine(q1), + In2G1: NewG1Affine(p2), + In2G2: NewG2Affine(q2), + Res: NewGTEl(res), + } + err = test.IsSolved(&MillerLoopCircuit{}, &witness, ecc.BN254.ScalarField()) + assert.NoError(err) +} + type FinalExponentiation struct { InGt GTEl Res GTEl diff --git a/std/algebra/emulated/sw_bw6761/hints.go b/std/algebra/emulated/sw_bw6761/hints.go new file mode 100644 index 000000000..796651934 --- /dev/null +++ b/std/algebra/emulated/sw_bw6761/hints.go @@ -0,0 +1,111 @@ +package sw_bw6761 + +import ( + "math/big" + + bw6761 "github.com/consensys/gnark-crypto/ecc/bw6-761" + "github.com/consensys/gnark/constraint/solver" + "github.com/consensys/gnark/std/math/emulated" +) + +func init() { + solver.RegisterHint(GetHints()...) +} + +// GetHints returns all hint functions used in the package. +func GetHints() []solver.Hint { + return []solver.Hint{ + finalExpHint, + pairingCheckHint, + } +} + +func finalExpHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { + // This adapted from section 4.3.2 of https://eprint.iacr.org/2024/640.pdf + return emulated.UnwrapHint(nativeInputs, nativeOutputs, + func(mod *big.Int, inputs, outputs []*big.Int) error { + + var millerLoop bw6761.E6 + millerLoop.B0.A0.SetBigInt(inputs[0]) + millerLoop.B0.A1.SetBigInt(inputs[2]) + millerLoop.B0.A2.SetBigInt(inputs[4]) + millerLoop.B1.A0.SetBigInt(inputs[1]) + millerLoop.B1.A1.SetBigInt(inputs[3]) + millerLoop.B1.A2.SetBigInt(inputs[5]) + + // m = (x₀³-x₀²+1-p(x₀+1)) / r + // mInv = 1/m mod p^6-1/r + var mInv big.Int + mInv.SetString("420096572758781926988571022578549119077996267041217186563532964653013626327499627643558150289556860284699838191238508062761264485377946319676011525555582097381055209304464769241709045835179375847000286979304653199040198646948595850434830718773056593021324330541604029824826938177546414778934883707126835848724258610612114712835130017082970786784508470382396148858570586085402148355642863720286568566937773459407961735112550507047306343380386401338522186960986251395049985320677251315016812720092326581314645206610216409714397970562842517827716362494341171265008409446148022671451843025093584702610246849007545665518399731546205544005105929880663530772806759681913801835273987094997504640832304570158760940364827187477825525048007459079382410480491250884588399683894539404567701993526561088158396861020181640181843560309670937868772703282755078557149854363818903590441797744966016708880143332350534049482338696654635346189790575286999280892407997722996866724226514621504774811766428733682155766330614074143245300182851212177081558245259537898592443393875891588079021560334726750431309338787970594548465289737362624558256642461612913108676326999205533110217714096123782036214164015261929502119392490941988919030563789520985909704716341786823561745842985678563", 10) + + residueWitness := finalExpWitness(&millerLoop, &mInv) + + residueWitness.B0.A0.BigInt(outputs[0]) + residueWitness.B0.A1.BigInt(outputs[2]) + residueWitness.B0.A2.BigInt(outputs[4]) + residueWitness.B1.A0.BigInt(outputs[1]) + residueWitness.B1.A1.BigInt(outputs[3]) + residueWitness.B1.A2.BigInt(outputs[5]) + + return nil + }) +} + +func pairingCheckHint(nativeMod *big.Int, nativeInputs, nativeOutputs []*big.Int) error { + // This adapted from section 4.3.2 of https://eprint.iacr.org/2024/640.pdf + return emulated.UnwrapHint(nativeInputs, nativeOutputs, + func(mod *big.Int, inputs, outputs []*big.Int) error { + var P bw6761.G1Affine + var Q bw6761.G2Affine + n := len(inputs) + p := make([]bw6761.G1Affine, 0, n/4) + q := make([]bw6761.G2Affine, 0, n/4) + for k := 0; k < n/4+1; k += 2 { + P.X.SetBigInt(inputs[k]) + P.Y.SetBigInt(inputs[k+1]) + p = append(p, P) + } + for k := n / 2; k < n; k += 2 { + Q.X.SetBigInt(inputs[k]) + Q.Y.SetBigInt(inputs[k+1]) + q = append(q, Q) + } + + millerLoop, err := bw6761.MillerLoopDirect(p, q) + if err != nil { + return err + } + + // m = (x₀+1+p(x₀³-x₀²-x₀)) / r + // mInv = 1/m mod p^6-1/r + var mInv big.Int + mInv.SetString("105300887666978464659709343582542432109497460559010677145223399327335567156593762277982229043678237863242655241846768823344862796112034076814141083092751207576412334798103601349742476585775877619451019850167305863473223932142842098178714149254582966792063312581807532675011404956270444910983750120675327025908192761069674135173328190635728173483753211505851991073745950587829640934449952514784880889959559541546684726344944253403018397996950965921029567425987659358091464001225755716260618839676545930683009926269854751616319103606509390667378268460666742713527948268373325914395974070631687649214144656759247037859773349886114399692016935966157297580328600396352321897692663748248168657388300690175586203114387947411720168269584172401784701771662759756974275902513788431327670950496435721956320875507468132703494465092748348925165286946843554008708392819919707156205920861214337368776935547492934209453494196115576830279851512338758088097719490141268227027970070242059962020992385206924254152017997017283665944910844784993588814611604460594039341562723060932582754994971346320340801549001828241339646153773031765187339622798156846331769418880530957782348437016822638577491500694745694281480857816937650066502281171825041093314285283892479458782481150957342407", 10) + + residueWitnessInv := finalExpWitness(&millerLoop, &mInv) + residueWitnessInv.Inverse(&residueWitnessInv) + + residueWitnessInv.B0.A0.BigInt(outputs[0]) + residueWitnessInv.B0.A1.BigInt(outputs[2]) + residueWitnessInv.B0.A2.BigInt(outputs[4]) + residueWitnessInv.B1.A0.BigInt(outputs[1]) + residueWitnessInv.B1.A1.BigInt(outputs[3]) + residueWitnessInv.B1.A2.BigInt(outputs[5]) + + return nil + }) +} + +func finalExpWitness(millerLoop *bw6761.E6, mInv *big.Int) (residueWitness bw6761.E6) { + + var rInv big.Int + // 1. compute r-th root: + // Exponentiate to rInv where + // rInv = 1/r mod (p^6-1)/r + rInv.SetString("279142441805511726233822077180198394933430419224185936052953462287387912118470357993263103168031788043160461358474005435622327506926362567154401645657309519073154383052970657693950208844465818979551693587858245321454505472049236704031061301292776853925224359757586505231126091244204292668007110271845616234279927419974150119801003450133674289144711275201991607282264849765236206295842916353255855388186086438329721887082685697023028663652777877691341551982676874308309620809049793085180324511691754953492619183755890255644855765188965000691813063771086522132765764526955251054211157804606693386854395171192876178005945476647006847460976477055233044799299417913662363985523123796056692751028712679181978298499780752966303529102009307348414562366180130429432094237007700663759126264893082917308542509779442201840676518234962495304673134599305371982876385622279935346701152286347948653741121231188575146952014672242471261647823749129902237689180055673361938161119768341970519416039779128617354778773830515364777252518313057683396662835013368967463878342754251509207391537635831891662211848811733884861792121210263430418966889668537646457064092991696527814120385172941004264289812969796992647021735186941896252860419364971543301451924917610828019341224722038007513", 10) + residueWitness.Exp(*millerLoop, &rInv) + + // 2. compute m-th root: + residueWitness.Exp(residueWitness, mInv) + + return residueWitness +} diff --git a/std/algebra/emulated/sw_bw6761/pairing.go b/std/algebra/emulated/sw_bw6761/pairing.go index c84170117..70dc2f380 100644 --- a/std/algebra/emulated/sw_bw6761/pairing.go +++ b/std/algebra/emulated/sw_bw6761/pairing.go @@ -115,6 +115,41 @@ func (pr Pairing) FinalExponentiation(z *GTEl) *GTEl { return result } +// AssertFinalExponentiationIsOne 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 is adapted from Section 4 of [On Proving Pairings] paper by A. Novakovic and L. Eagen. +// +// [On Proving Pairings]: https://eprint.iacr.org/2024/640.pdf +func (pr Pairing) AssertFinalExponentiationIsOne(x *GTEl) { + res, err := pr.curveF.NewHint(finalExpHint, 6, &x.A0, &x.A1, &x.A2, &x.A3, &x.A4, &x.A5) + if err != nil { + // err is non-nil only for invalid number of inputs + panic(err) + } + + residueWitness := GTEl{ + A0: *res[0], + A1: *res[1], + A2: *res[2], + A3: *res[3], + A4: *res[4], + A5: *res[5], + } + + // Check that x == residueWitness^Λ + // where Λ = x₀+1+p(x₀^3-x₀^2-x₀) and residueWitness from the hint. + + // exponentiation by U1=x₀^3-x₀^2-x₀ + t0 := pr.Ext6.ExpByU1(&residueWitness) + t0 = pr.Ext6.Frobenius(t0) + // exponentiation by U2=x₀+1 + t1 := pr.Ext6.ExpByU2(&residueWitness) + t0 = pr.Ext6.Mul(t0, t1) + + pr.AssertIsEqual(t0, x) +} + // Pair calculates the reduced pairing for a set of points // ∏ᵢ e(Pᵢ, Qᵢ). // @@ -132,22 +167,58 @@ 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.MillerLoop(P, Q) + // check input size match + nP := len(P) + nQ := len(Q) + if nP == 0 || nP != nQ { + return errors.New("invalid inputs sizes") + } + // hint the non-residue witness + inputs := make([]*baseEl, 0, 2*nP+2*nQ) + for _, p := range P { + inputs = append(inputs, &p.X, &p.Y) + } + for _, q := range Q { + inputs = append(inputs, &q.P.X, &q.P.Y) + } + hint, err := pr.curveF.NewHint(pairingCheckHint, 6, inputs...) if err != nil { - return err + // err is non-nil only for invalid number of inputs + panic(err) + } + residueWitnessInv := >El{ + A0: *hint[0], + A1: *hint[1], + A2: *hint[2], + A3: *hint[3], + A4: *hint[4], + A5: *hint[5], + } + + lines := make([]lineEvaluations, len(Q)) + for i := range Q { + if Q[i].Lines == nil { + Qlines := pr.computeLines(&Q[i].P) + Q[i].Lines = &Qlines + } + lines[i] = *Q[i].Lines } - // We perform the easy part of the final exp to push f to the cyclotomic - // subgroup so that AssertFinalExponentiationIsOne 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.Frobenius(buf) - f = pr.Mul(f, buf) - pr.AssertFinalExponentiationIsOne(f) + res, err := pr.millerLoopLines(P, lines, residueWitnessInv, false) + if err != nil { + return fmt.Errorf("miller loop: %w", err) + } + + // Check that: MillerLoop(P,Q) == residueWitness^Λ + // where Λ = x₀+1+p(x₀³-x₀²-x₀) and residueWitness from the hint. + // + // Note that at this point: + // result = MillerLoop(P,Q) * residueWitnessInv^{x₀+1+p(x₀³-x₀²-x₀)} + // since we initialized the Miller loop accumulator with residueWitnessInv^{p}. + // So we only need to check that: + // result == 1. + pr.AssertIsEqual(res, pr.Ext6.One()) return nil } @@ -273,12 +344,12 @@ func (pr Pairing) MillerLoop(P []*G1Affine, Q []*G2Affine) (*GTEl, error) { } lines[i] = *Q[i].Lines } - return pr.millerLoopLines(P, lines) + return pr.millerLoopLines(P, lines, nil, true) } // millerLoopLines computes the multi-Miller loop from points in G1 and precomputed lines in G2 -func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl, error) { +func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations, init *GTEl, first bool) (*GTEl, error) { // check input size match n := len(P) @@ -299,53 +370,88 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl xNegOverY[k] = pr.curveF.Neg(xNegOverY[k]) } - // f_{x₀+1+λ(x₀³-x₀²-x₀),Q}(P), Q is known in advance + // Compute f_{x₀+1+λ(x₀³-x₀²-x₀),Q}(P) var prodLines [5]*baseEl result := pr.Ext6.One() - // i = 188 - // k = 0 - result = &fields_bw6761.E6{ - A0: *pr.curveF.Mul(&lines[0][0][188].R1, yInv[0]), - A1: result.A1, - A2: *pr.curveF.Mul(&lines[0][0][188].R0, xNegOverY[0]), - A3: *pr.curveF.One(), - A4: result.A4, - A5: result.A5, + var initInv, frobInit, frobInitInv GTEl + if init != nil { + initInv = *pr.Ext6.Inverse(init) + frobInit = *pr.Ext6.Frobenius(init) + frobInitInv = *pr.Ext6.Frobenius(&initInv) + result = &frobInit } - if n >= 2 { - // k = 1, separately to avoid MulBy023 (res × ℓ) - // (res is also a line at this point, so we use Mul023By023 ℓ × ℓ) - prodLines = pr.Mul023By023( - pr.curveF.Mul(&lines[1][0][188].R1, yInv[1]), - pr.curveF.Mul(&lines[1][0][188].R0, xNegOverY[1]), - &result.A0, - &result.A2, - ) - result = &fields_bw6761.E6{ - A0: *prodLines[0], + j := len(loopCounter2) - 2 + if first { + // i = j + // k = 0 + result = >El{ + A0: *pr.curveF.Mul(&lines[0][0][j].R1, yInv[0]), A1: result.A1, - A2: *prodLines[1], - A3: *prodLines[2], - A4: *prodLines[3], - A5: *prodLines[4], + A2: *pr.curveF.Mul(&lines[0][0][j].R0, xNegOverY[0]), + A3: *pr.curveF.One(), + A4: result.A4, + A5: result.A5, } - } - for k := 2; k < n; k++ { - result = pr.MulBy023(result, - pr.curveF.Mul(&lines[k][0][188].R1, yInv[k]), - pr.curveF.Mul(&lines[k][0][188].R0, xNegOverY[k]), - ) + if n >= 2 { + // k = 1, separately to avoid MulBy023 (res × ℓ) + // (res is also a line at this point, so we use Mul023By023 ℓ × ℓ) + prodLines = pr.Mul023By023( + pr.curveF.Mul(&lines[1][0][j].R1, yInv[1]), + pr.curveF.Mul(&lines[1][0][j].R0, xNegOverY[1]), + &result.A0, + &result.A2, + ) + result = >El{ + A0: *prodLines[0], + A1: result.A1, + A2: *prodLines[1], + A3: *prodLines[2], + A4: *prodLines[3], + A5: *prodLines[4], + } + } + + for k := 2; k < n; k++ { + result = pr.MulBy023(result, + pr.curveF.Mul(&lines[k][0][j].R1, yInv[k]), + pr.curveF.Mul(&lines[k][0][j].R0, xNegOverY[k]), + ) + } + j-- } - for i := 187; i >= 0; i-- { + for i := j; i > 0; i-- { // mutualize the square among n Miller loops // (∏ᵢfᵢ)² result = pr.Square(result) - if i > 0 && loopCounter2[i]*3+loopCounter1[i] != 0 { + j := loopCounter1[i] + 3*loopCounter2[i] + switch j { + // cases -4, -2, 2, 4 do not occur, given the static LoopCounters + case -3: + if init != nil { + // mul by frobInitInv to capture -1's in x₀³-x₀²-x₀ + result = pr.Ext6.Mul(result, &frobInitInv) + } + // mul by tangent and line + for k := 0; k < n; k++ { + prodLines = pr.Mul023By023( + pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]), + pr.curveF.Mul(&lines[k][0][i].R0, xNegOverY[k]), + pr.curveF.Mul(&lines[k][1][i].R1, yInv[k]), + pr.curveF.Mul(&lines[k][1][i].R0, xNegOverY[k]), + ) + result = pr.MulBy02345(result, prodLines) + } + case -1: + if init != nil { + // mul by initInv to capture -1's in x₀+1 + result = pr.Ext6.Mul(result, &initInv) + } + // mul by tangent and line for k := 0; k < n; k++ { prodLines = pr.Mul023By023( pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]), @@ -355,8 +461,18 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl ) result = pr.MulBy02345(result, prodLines) } - } else { - // if number of lines is odd, mul last line by res + case 0: + // mul tangents 2-by-2 and then by accumulator + for k := 1; k < n; k += 2 { + prodLines = pr.Mul023By023( + pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]), + pr.curveF.Mul(&lines[k][0][i].R0, xNegOverY[k]), + pr.curveF.Mul(&lines[k-1][0][i].R1, yInv[k-1]), + pr.curveF.Mul(&lines[k-1][0][i].R0, xNegOverY[k-1]), + ) + result = pr.MulBy02345(result, prodLines) + } + // if number of tangents is odd, mul last line by res // works for n=1 as well if n%2 != 0 { // ℓ × res @@ -365,19 +481,68 @@ func (pr Pairing) millerLoopLines(P []*G1Affine, lines []lineEvaluations) (*GTEl pr.curveF.Mul(&lines[n-1][0][i].R0, xNegOverY[n-1]), ) } - // mul lines 2-by-2 - for k := 1; k < n; k += 2 { + case 1: + if init != nil { + // mul by init to capture 1's in x₀+1 + result = pr.Ext6.Mul(result, init) + } + // mul by line and tangent + for k := 0; k < n; k++ { prodLines = pr.Mul023By023( pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]), pr.curveF.Mul(&lines[k][0][i].R0, xNegOverY[k]), - pr.curveF.Mul(&lines[k-1][0][i].R1, yInv[k-1]), - pr.curveF.Mul(&lines[k-1][0][i].R0, xNegOverY[k-1]), + pr.curveF.Mul(&lines[k][1][i].R1, yInv[k]), + pr.curveF.Mul(&lines[k][1][i].R0, xNegOverY[k]), ) result = pr.MulBy02345(result, prodLines) } + case 3: + if init != nil { + // mul by frobInit to capture 1's in x₀³-x₀²-x₀ + result = pr.Ext6.Mul(result, &frobInit) + } + for k := 0; k < n; k++ { + prodLines = pr.Mul023By023( + pr.curveF.Mul(&lines[k][0][i].R1, yInv[k]), + pr.curveF.Mul(&lines[k][0][i].R0, xNegOverY[k]), + pr.curveF.Mul(&lines[k][1][i].R1, yInv[k]), + pr.curveF.Mul(&lines[k][1][i].R0, xNegOverY[k]), + ) + result = pr.MulBy02345(result, prodLines) + } + default: + panic("unknown case for loopCounter") } } + // i = 0 (j = -3) + result = pr.Square(result) + if init != nil { + // mul by frobInitInv to capture -1's in x₀³-x₀²-x₀ + result = pr.Ext6.Mul(result, &frobInitInv) + } + // x₀+1+λ(x₀³-x₀²-x₀) = 0 mod r so accQ = ∞ at the last iteration, + // we only mul by tangent. + // mul tangents 2-by-2 and then by accumulator + for k := 1; k < n; k += 2 { + prodLines = pr.Mul023By023( + pr.curveF.Mul(&lines[k][0][0].R1, yInv[k]), + pr.curveF.Mul(&lines[k][0][0].R0, xNegOverY[k]), + pr.curveF.Mul(&lines[k-1][0][0].R1, yInv[k-1]), + pr.curveF.Mul(&lines[k-1][0][0].R0, xNegOverY[k-1]), + ) + result = pr.MulBy02345(result, prodLines) + } + // if number of tangents is odd, mul last line by res + // works for n=1 as well + if n%2 != 0 { + // ℓ × res + result = pr.MulBy023(result, + pr.curveF.Mul(&lines[n-1][0][0].R1, yInv[n-1]), + pr.curveF.Mul(&lines[n-1][0][0].R0, xNegOverY[n-1]), + ) + } + return result, nil } diff --git a/std/algebra/emulated/sw_bw6761/pairing_test.go b/std/algebra/emulated/sw_bw6761/pairing_test.go index 65f087a55..4c10432ce 100644 --- a/std/algebra/emulated/sw_bw6761/pairing_test.go +++ b/std/algebra/emulated/sw_bw6761/pairing_test.go @@ -33,6 +33,45 @@ func randomG1G2Affines() (bw6761.G1Affine, bw6761.G2Affine) { return p, q } +type MillerLoopCircuit struct { + In1G1, In2G1 G1Affine + In1G2, In2G2 G2Affine + Res GTEl +} + +func (c *MillerLoopCircuit) Define(api frontend.API) error { + pairing, err := NewPairing(api) + if err != nil { + return fmt.Errorf("new pairing: %w", err) + } + res, err := pairing.MillerLoop([]*G1Affine{&c.In1G1, &c.In2G1}, []*G2Affine{&c.In1G2, &c.In2G2}) + if err != nil { + return fmt.Errorf("Miller loop: %w", err) + } + pairing.AssertIsEqual(res, &c.Res) + return nil +} + +func TestMillerLoopTestSolve(t *testing.T) { + assert := test.NewAssert(t) + p1, q1 := randomG1G2Affines() + p2, q2 := randomG1G2Affines() + res, err := bw6761.MillerLoopDirect( + []bw6761.G1Affine{p1, p2}, + []bw6761.G2Affine{q1, q2}, + ) + assert.NoError(err) + witness := MillerLoopCircuit{ + In1G1: NewG1Affine(p1), + In1G2: NewG2Affine(q1), + In2G1: NewG1Affine(p2), + In2G2: NewG2Affine(q2), + Res: NewGTEl(res), + } + err = test.IsSolved(&MillerLoopCircuit{}, &witness, ecc.BN254.ScalarField()) + assert.NoError(err) +} + type FinalExponentiationCircuit struct { InGt GTEl Res GTEl diff --git a/std/algebra/emulated/sw_bw6761/precomputations.go b/std/algebra/emulated/sw_bw6761/precomputations.go index 829dd1b2b..f91db8173 100644 --- a/std/algebra/emulated/sw_bw6761/precomputations.go +++ b/std/algebra/emulated/sw_bw6761/precomputations.go @@ -56,6 +56,9 @@ func (p *Pairing) computeLines(Q *g2AffP) lineEvaluations { panic("unknown case for loopCounter") } } + // i = 0 (case -3) + // x₀+1+λ(x₀³-x₀²-x₀) = 0 mod r so accQ = ∞ at the last iteration, + // we only compute the tangent. cLines[0][0] = p.tangentCompute(accQ) return cLines } diff --git a/std/algebra/native/sw_bls12377/hints.go b/std/algebra/native/sw_bls12377/hints.go index 6a993b909..519e577b9 100644 --- a/std/algebra/native/sw_bls12377/hints.go +++ b/std/algebra/native/sw_bls12377/hints.go @@ -18,6 +18,7 @@ func GetHints() []solver.Hint { scalarMulGLVG1Hint, halfGCDEisenstein, halfGCDEisensteinSigns, + pairingCheckHint, } } @@ -25,6 +26,90 @@ func init() { solver.RegisterHint(GetHints()...) } +func pairingCheckHint(scalarField *big.Int, inputs, outputs []*big.Int) error { + var P bls12377.G1Affine + var Q bls12377.G2Affine + n := len(inputs) + p := make([]bls12377.G1Affine, 0, n/6) + q := make([]bls12377.G2Affine, 0, n/6) + for k := 0; k < n/6+1; k += 2 { + P.X.SetBigInt(inputs[k]) + P.Y.SetBigInt(inputs[k+1]) + p = append(p, P) + } + for k := n / 3; k < n/2+3; k += 4 { + Q.X.A0.SetBigInt(inputs[k]) + Q.X.A1.SetBigInt(inputs[k+1]) + Q.Y.A0.SetBigInt(inputs[k+2]) + Q.Y.A1.SetBigInt(inputs[k+3]) + q = append(q, Q) + } + + lines := make([][2][len(bls12377.LoopCounter) - 1]bls12377.LineEvaluationAff, 0, len(q)) + for _, qi := range q { + lines = append(lines, bls12377.PrecomputeLines(qi)) + } + millerLoop, err := bls12377.MillerLoopFixedQ(p, lines) + if err != nil { + return err + } + + var root, rootPthInverse, residueWitness, scalingFactor bls12377.E12 + var exponent, exponentInv, finalExpFactor, polyFactor big.Int + // polyFactor = 12(x-1) + polyFactor.SetString("115033474957087604736", 10) + // finalExpFactor = ((q^12 - 1) / r) / polyFactor + finalExpFactor.SetString("92351561334497520756349650336409370070948672672207914824247073415859727964231807559847070685040742345026775319680739143654748316009031763764029886042408725311062057776702838555815712331129279611544378217895455619058809454575474763035923260395518532422855090028311239234310116353269618927871828693919559964406939845784130633021661399269804065961999062695977580539176029238189119059338698461832966347603096853909366901376879505972606045770762516580639801134008192256366142553202619529638202068488750102055204336502584141399828818871664747496033599618827160583206926869573005874449182200210044444351826855938563862937638034918413235278166699461287943529570559518592586872860190313088429391521694808994276205429071153237122495989095857292965461625387657577981811772819764071512345106346232882471034669258055302790607847924560040527682025558360106509628206144255667203317787586698694011876342903106644003067103035176245790275561392007119121995936066014208972135762663107247939004517852248103325700169848524693333524025685325993207375736519358185783520948988673594976115901587076295116293065682366935313875411927779217584729138600463438806153265891176654957439524358472291492028580820575807385461119025678550977847392818655362610734928283105671242634809807533919011078145", 10) + + // 1. get pth-root inverse + exponent.Set(&finalExpFactor) + root.Exp(millerLoop, &finalExpFactor) + if root.IsOne() { + rootPthInverse.SetOne() + } else { + exponentInv.ModInverse(&exponent, &polyFactor) + exponent.Neg(&exponentInv).Mod(&exponent, &polyFactor) + rootPthInverse.Exp(root, &exponent) + } + + // 3. shift the Miller loop result so that millerLoop * scalingFactor + // is of order finalExpFactor + scalingFactor.Set(&rootPthInverse) + millerLoop.Mul(&millerLoop, &scalingFactor) + + // 4. get the witness residue + // + // lambda = q - u, the optimal exponent + var lambda big.Int + lambda.SetString("258664426012969094010652733694893533536393512754914660539884262666720468348340822774968888139563774001527230824448", 10) + exponent.ModInverse(&lambda, &finalExpFactor) + residueWitness.Exp(millerLoop, &exponent) + + // return the witness residue + 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 the scaling factor + scalingFactor.C0.B0.A0.BigInt(outputs[12]) + scalingFactor.C0.B0.A1.BigInt(outputs[13]) + scalingFactor.C0.B1.A0.BigInt(outputs[14]) + scalingFactor.C0.B1.A1.BigInt(outputs[15]) + scalingFactor.C0.B2.A0.BigInt(outputs[16]) + scalingFactor.C0.B2.A1.BigInt(outputs[17]) + + return nil +} + func decomposeScalarG1Simple(scalarField *big.Int, inputs []*big.Int, outputs []*big.Int) error { if len(inputs) != 1 { return errors.New("expecting one input") diff --git a/std/algebra/native/sw_bls12377/pairing.go b/std/algebra/native/sw_bls12377/pairing.go index e4cf8f38d..1756dc41a 100644 --- a/std/algebra/native/sw_bls12377/pairing.go +++ b/std/algebra/native/sw_bls12377/pairing.go @@ -240,22 +240,118 @@ func Pair(api frontend.API, P []G1Affine, Q []G2Affine) (GT, error) { // // This function doesn't check that the inputs are in the correct subgroups func PairingCheck(api frontend.API, P []G1Affine, Q []G2Affine) error { - f, err := MillerLoop(api, P, Q) + + // check input size match + nP := len(P) + nQ := len(Q) + if nP == 0 || nP != nQ { + return errors.New("invalid inputs sizes") + } + // hint the non-residue witness + inputs := make([]frontend.Variable, 0, 2*nP+4*nQ) + for _, p := range P { + inputs = append(inputs, p.X, p.Y) + } + for _, q := range Q { + inputs = append(inputs, q.P.X.A0, q.P.X.A1, q.P.Y.A0, q.P.Y.A1) + } + hint, err := api.NewHint(pairingCheckHint, 18, inputs...) if err != nil { - return err + // err is non-nil only for invalid number of inputs + panic(err) + } + var residueWitness GT + residueWitness.C0.B0.A0 = hint[0] + residueWitness.C0.B0.A1 = hint[1] + residueWitness.C0.B1.A0 = hint[2] + residueWitness.C0.B1.A1 = hint[3] + residueWitness.C0.B2.A0 = hint[4] + residueWitness.C0.B2.A1 = hint[5] + residueWitness.C1.B0.A0 = hint[6] + residueWitness.C1.B0.A1 = hint[7] + residueWitness.C1.B1.A0 = hint[8] + residueWitness.C1.B1.A1 = hint[9] + residueWitness.C1.B2.A0 = hint[10] + residueWitness.C1.B2.A1 = hint[11] + + var scalingFactor fields_bls12377.E6 + // constrain cubicNonResiduePower to be in Fp6 + scalingFactor.B0.A0 = hint[12] + scalingFactor.B0.A1 = hint[13] + scalingFactor.B1.A0 = hint[14] + scalingFactor.B1.A1 = hint[15] + scalingFactor.B2.A0 = hint[16] + scalingFactor.B2.A1 = hint[17] + + lines := make([]lineEvaluations, nQ) + for i := range Q { + if Q[i].Lines == nil { + Qlines := computeLines(api, Q[i].P) + Q[i].Lines = Qlines + } + lines[i] = *Q[i].Lines } - // We perform the easy part of the final exp to push f to the cyclotomic - // subgroup so that AssertFinalExponentiationIsOne is carried with optimized - // cyclotomic squaring (e.g. Karabina12345). - // - // f = f^(p⁶-1)(p²+1) - var buf GT - buf.Conjugate(api, f) - buf.DivUnchecked(api, buf, f) - f.FrobeniusSquare(api, buf). - Mul(api, f, buf) - - f.AssertFinalExponentiationIsOne(api) + + // precomputations + yInv := make([]frontend.Variable, nP) + xNegOverY := make([]frontend.Variable, nP) + for k := 0; k < nP; k++ { + yInv[k] = api.DivUnchecked(1, P[k].Y) + xNegOverY[k] = api.Mul(P[k].X, yInv[k]) + xNegOverY[k] = api.Neg(xNegOverY[k]) + } + + // init Miller loop accumulator to residueWitness to share the squarings + // of residueWitness^{x₀} + res := residueWitness + + var prodLines [5]fields_bls12377.E2 + var l0, l1 lineEvaluation + + // Compute ∏ᵢ { fᵢ_{x₀,Q}(P) } + for i := 62; i >= 0; i-- { + // mutualize the square among n Miller loops + // (∏ᵢfᵢ)² + res.Square(api, res) + + if loopCounter[i] == 0 { + for k := 0; k < nP; k++ { + // line evaluation at P + // ℓ × res + res.MulBy034(api, + *l0.R0.MulByFp(api, lines[k][0][i].R0, xNegOverY[k]), + *l0.R1.MulByFp(api, lines[k][0][i].R1, yInv[k]), + ) + } + } else { + // multiply by residueWitness when bit=1 + res.Mul(api, res, residueWitness) + for k := 0; k < nP; k++ { + // lines evaluation at P + // ℓ × ℓ + prodLines = *fields_bls12377.Mul034By034(api, + *l0.R0.MulByFp(api, lines[k][0][i].R0, xNegOverY[k]), + *l0.R1.MulByFp(api, lines[k][0][i].R1, yInv[k]), + *l1.R0.MulByFp(api, lines[k][1][i].R0, xNegOverY[k]), + *l1.R1.MulByFp(api, lines[k][1][i].R1, yInv[k]), + ) + // (ℓ × ℓ) × res + res.MulBy01234(api, prodLines) + } + } + } + + // Check that res * scalingFactor == residueWitness^(q) + // where u=0x8508c00000000001 is the BLS12-377 seed, + // and residueWitness, scalingFactor from the hint. + // Note that res is already MillerLoop(P,Q) * residueWitness^{x₀} since + // we initialized the Miller loop accumulator with residueWitness. + var t0, t1 GT + t1.C0.Mul(api, res.C0, scalingFactor) + t1.C1.Mul(api, res.C1, scalingFactor) + t0.Frobenius(api, residueWitness) + + t0.AssertIsEqual(api, t1) return nil } diff --git a/std/algebra/native/sw_bls12377/pairing2.go b/std/algebra/native/sw_bls12377/pairing2.go index 072e2707b..122c61ff6 100644 --- a/std/algebra/native/sw_bls12377/pairing2.go +++ b/std/algebra/native/sw_bls12377/pairing2.go @@ -313,21 +313,11 @@ func (p *Pairing) PairingCheck(P []*G1Affine, Q []*G2Affine) error { for i := range Q { inQ[i] = *Q[i] } - res, err := MillerLoop(p.api, inP, inQ) + err := PairingCheck(p.api, inP, inQ) if err != nil { return err } - // We perform the easy part of the final exp to push res to the cyclotomic - // subgroup so that AssertFinalExponentiationIsOne is carried with optimized - // cyclotomic squaring (e.g. Karabina12345). - // - // res = res^(p⁶-1)(p²+1) - var buf GT - buf.Conjugate(p.api, res) - buf.DivUnchecked(p.api, buf, res) - res.FrobeniusSquare(p.api, buf).Mul(p.api, res, buf) - - res.AssertFinalExponentiationIsOne(p.api) + return nil }