perf(bw6): use more efficient addchains

std/algebra/emulated/fields_bw6761/e3.go

@@ -173,17 +173,17 @@ func (e Ext3) MulBy01(z *E3, c0, c1 *baseEl) *E3 {
 }
 
 // MulBy1 multiplication of E6 by sparse element (0, c1, 0)
-func (e Ext3) MulBy1(z *E3, c1 baseEl) *E3 {
+func (e Ext3) MulBy1(z *E3, c1 *baseEl) *E3 {
 
-	b := e.fp.Mul(&z.A1, &c1)
+	b := e.fp.Mul(&z.A1, c1)
 
 	tmp := e.fp.Add(&z.A1, &z.A2)
-	t0 := e.fp.Mul(&c1, tmp)
+	t0 := e.fp.Mul(c1, tmp)
 	t0 = e.fp.Sub(t0, b)
 	t0 = mulFpByNonResidue(e.fp, t0)
 
 	tmp = e.fp.Add(&z.A0, &z.A1)
-	t1 := e.fp.Mul(&c1, tmp)
+	t1 := e.fp.Mul(c1, tmp)
 	t1 = e.fp.Sub(t1, b)
 
 	return &E3{

std/algebra/emulated/fields_bw6761/e6_pairing.go

@@ -7,102 +7,122 @@ func (e Ext6) nSquareCompressed(z *E6, n int) *E6 {
 	return z
 }
 
-// Expt set x to x^t in E6 and return x
-func (e Ext6) Expt(x *E6) *E6 {
-	x = e.Reduce(x)
-
-	// const tAbsVal uint64 = 9586122913090633729
-	// tAbsVal in binary: 1000010100001000110000000000000000000000000000000000000000000001
-	// drop the low 46 bits (all 0 except the least significant bit): 100001010000100011 = 136227
-	// Shortest addition chains can be found at https://wwwhomes.uni-bielefeld.de/achim/addition_chain.html
-
-	// a shortest addition chain for 136227
-	result := e.Copy(x)
+// ExpX0Minus1 set z to z^{x₀-1} in E6 and return z
+// x₀-1 = 91893752504881257682351033800651177983
+func (e Ext6) ExpX0Minus1(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.Copy(z)
 	result = e.nSquareCompressed(result, 5)
 	result = e.DecompressKarabina(result)
-	result = e.Mul(result, x)
-	x33 := e.Copy(result)
+	result = e.Mul(result, z)
+	z33 := e.Copy(result)
 	result = e.nSquareCompressed(result, 7)
 	result = e.DecompressKarabina(result)
-	result = e.Mul(result, x33)
+	result = e.Mul(result, z33)
 	result = e.nSquareCompressed(result, 4)
 	result = e.DecompressKarabina(result)
-	result = e.Mul(result, x)
+	result = e.Mul(result, z)
 	result = e.CyclotomicSquare(result)
-	result = e.Mul(result, x)
-
-	// the remaining 46 bits
+	result = e.Mul(result, z)
 	result = e.nSquareCompressed(result, 46)
 	result = e.DecompressKarabina(result)
-	result = e.Mul(result, x)
 
 	return result
 }
 
-// ExptMinus1 set x to x^(t-1) in E6 and return x
-func (e Ext6) ExptMinus1(x *E6) *E6 {
-	result := e.Expt(x)
-	result = e.Mul(result, e.Conjugate(x))
+// ExpX0Minus1Square set z to z^{(x₀-1)²} in E6 and return z
+// (x₀-1)² = 91893752504881257682351033800651177984
+func (e Ext6) ExpX0Minus1Square(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.Copy(z)
+	result = e.CyclotomicSquare(result)
+	t0 := e.Mul(z, result)
+	t1 := e.CyclotomicSquare(t0)
+	t0 = e.Mul(t0, t1)
+	result = e.Mul(result, t0)
+	t1 = e.Mul(t1, result)
+	t0 = e.Mul(t0, t1)
+	t2 := e.CyclotomicSquare(t0)
+	t2 = e.Mul(t1, t2)
+	t0 = e.Mul(t0, t2)
+	t2 = e.nSquareCompressed(t2, 7)
+	t2 = e.DecompressKarabina(t2)
+	t1 = e.Mul(t1, t2)
+	t1 = e.nSquareCompressed(t1, 11)
+	t1 = e.DecompressKarabina(t1)
+	t1 = e.Mul(t0, t1)
+	t1 = e.nSquareCompressed(t1, 9)
+	t1 = e.DecompressKarabina(t1)
+	t0 = e.Mul(t0, t1)
+	t0 = e.CyclotomicSquare(t0)
+	result = e.Mul(result, t0)
+	result = e.nSquareCompressed(result, 92)
+	result = e.DecompressKarabina(result)
+
 	return result
+
 }
 
-// ExptPlus1 set x to x^(t+1) in E6 and return x
-func (e Ext6) ExptPlus1(x *E6) *E6 {
-	result := e.Expt(x)
-	result = e.Mul(result, x)
+// ExpX0Plus1 set z to z^(x₀+1) in E6 and return z
+// x₀+1 = 91893752504881257682351033800651177985
+func (e Ext6) ExpX0Plus1(z *E6) *E6 {
+	result := e.ExpX0Minus1(z)
+	t := e.CyclotomicSquare(z)
+	result = e.Mul(result, t)
 	return result
 }
 
-// ExptMinus1Div3 set x to x^(t-1)/3 in E6 and return x
-func (e Ext6) ExptMinus1Div3(x *E6) *E6 {
-	x = e.Reduce(x)
-	result := e.Copy(x)
+// ExpX0Minus1Div3 set z to z^(x₀-1)/3 in E6 and return z
+// (x₀-1)/3 = 3195374304363544576
+func (e Ext6) ExptMinus1Div3(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.Copy(z)
 	result = e.CyclotomicSquare(result)
-	result = e.Mul(result, x)
-	t0 := e.Mul(result, x)
+	result = e.Mul(result, z)
+	t0 := e.Mul(result, z)
 	t0 = e.CyclotomicSquare(t0)
 	result = e.Mul(result, t0)
-	t0 = e.CyclotomicSquare(result)
-	t0 = e.nSquareCompressed(t0, 6)
+	t0 = result
+	t0 = e.nSquareCompressed(t0, 7)
 	t0 = e.DecompressKarabina(t0)
 	result = e.Mul(result, t0)
 	result = e.nSquareCompressed(result, 5)
 	result = e.DecompressKarabina(result)
-	result = e.Mul(result, x)
+	result = e.Mul(result, z)
 	result = e.nSquareCompressed(result, 46)
 	result = e.DecompressKarabina(result)
 
 	return result
 }
 
-// Expc1 set x to x^c2 in E6 and return x
+// ExpC1 set z to z^C1 in E6 and return z
 // ht, hy = 13, 9
-// c1 = (ht+hy)/2 = 11
-func (e Ext6) Expc1(x *E6) *E6 {
-	x = e.Reduce(x)
-	result := e.CyclotomicSquare(x)
-	result = e.Mul(result, x)
-	t0 := e.Mul(x, result)
+// C1 = (ht+hy)/2 = 11
+func (e Ext6) ExpC1(z *E6) *E6 {
+	z = e.Reduce(z)
+	result := e.CyclotomicSquare(z)
+	result = e.Mul(result, z)
+	t0 := e.Mul(z, result)
 	t0 = e.CyclotomicSquare(t0)
 	result = e.Mul(result, t0)
 
 	return result
 }
 
-// Expc2 set x to x^c1 in E6 and return x
+// ExpC2 set z to z^C2 in E6 and return z
 // ht, hy = 13, 9
-// c1 = (ht**2+3*hy**2)/4 = 103
-func (e Ext6) Expc2(x *E6) *E6 {
-	x = e.Reduce(x)
+// C2 = (ht**2+3*hy**2)/4 = 103
+func (e Ext6) ExpC2(z *E6) *E6 {
+	z = e.Reduce(z)
 
-	result := e.CyclotomicSquare(x)
-	result = e.Mul(result, x)
-	t0 := e.CyclotomicSquare(result)
-	t0 = e.nSquareCompressed(t0, 3)
+	result := e.CyclotomicSquare(z)
+	result = e.Mul(result, z)
+	t0 := result
+	t0 = e.nSquareCompressed(t0, 4)
 	t0 = e.DecompressKarabina(t0)
 	result = e.Mul(result, t0)
 	result = e.CyclotomicSquare(result)
-	result = e.Mul(result, x)
+	result = e.Mul(result, z)
 
 	return result
 }

std/algebra/emulated/sw_bw6761/pairing.go

@@ -66,32 +66,30 @@ func (pr Pairing) FinalExponentiation(z *GTEl) *GTEl {
 	// 2. Hard part (up to permutation)
 	// (x₀+1)(p²-p+1)/r
 	// Algorithm 4.4 from https://yelhousni.github.io/phd.pdf
-	a := pr.ExptMinus1(result)
-	a = pr.ExptMinus1(a)
+	a := pr.ExpX0Minus1Square(result)
 	a = pr.Mul(a, pr.Frobenius(result))
-	b := pr.ExptPlus1(a)
+	b := pr.ExpX0Plus1(a)
 	b = pr.Mul(b, pr.Conjugate(result))
 	t := pr.CyclotomicSquare(a)
 	a = pr.Mul(a, t)
 	c := pr.ExptMinus1Div3(b)
-	d := pr.ExptMinus1(c)
-	e := pr.ExptMinus1(d)
-	e = pr.ExptMinus1(e)
+	d := pr.ExpX0Minus1(c)
+	e := pr.ExpX0Minus1Square(d)
 	e = pr.Mul(e, d)
 	d = pr.Conjugate(d)
 	f := pr.Mul(d, b)
-	g := pr.ExptPlus1(e)
+	g := pr.ExpX0Plus1(e)
 	g = pr.Mul(g, f)
 	h := pr.Mul(g, c)
 	i := pr.Mul(g, d)
-	i = pr.ExptPlus1(i)
+	i = pr.ExpX0Plus1(i)
 	i = pr.Mul(i, pr.Conjugate(f))
-	j := pr.Expc1(h)
+	j := pr.ExpC1(h)
 	j = pr.Mul(j, e)
 	k := pr.CyclotomicSquare(j)
 	k = pr.Mul(k, j)
 	k = pr.Mul(k, b)
-	t = pr.Expc2(i)
+	t = pr.ExpC2(i)
 	k = pr.Mul(k, t)
 	result = pr.Mul(a, k)
 

std/algebra/emulated/sw_bw6761/pairing_test.go

@@ -45,7 +45,7 @@ func (c *FinalExponentiationCircuit) Define(api frontend.API) error {
 	}
 	res := pairing.FinalExponentiation(&c.InGt)
 	// gnark-crypto exponent is 12(x₀+1) and gnark is (x₀-1)
-	tmp := pairing.Expc1(res)
+	tmp := pairing.ExpC1(res)
 	res = pairing.Mul(res, tmp)
 	pairing.AssertIsEqual(res, &c.Res)
 	return nil
@@ -77,7 +77,7 @@ func (c *PairCircuit) Define(api frontend.API) error {
 	}
 	res, err := pairing.Pair([]*G1Affine{&c.InG1}, []*G2Affine{&c.InG2})
 	// gnark-crypto exponent is 12(x₀+1) and gnark is (x₀-1)
-	tmp := pairing.Expc1(res)
+	tmp := pairing.ExpC1(res)
 	res = pairing.Mul(res, tmp)
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)
@@ -118,7 +118,7 @@ func (c *MultiPairCircuit) Define(api frontend.API) error {
 		Q = append(Q, &c.InG2)
 	}
 	res, err := pairing.Pair(P, Q)
-	tmp := pairing.Expc1(res)
+	tmp := pairing.ExpC1(res)
 	res = pairing.Mul(res, tmp)
 	if err != nil {
 		return fmt.Errorf("pair: %w", err)