From ede59c5caa02fc7ef18e27a76cf64f9d88de71a9 Mon Sep 17 00:00:00 2001
From: armfazh <armfazh@cloudflare.com>
Date: Sat, 15 Jul 2023 17:30:02 -0700
Subject: [PATCH] Implement Granger-Scott faster squaring in the cyclotomic
 subgroup.

---
 ecc/bls12381/ff/cyclo6.go         | 37 ++++++++++++++++++++++++++++++-
 ecc/bls12381/ff/cyclo6_test.go    | 13 +++++++++++
 ecc/bls12381/ff/fp12cubic.go      |  2 +-
 ecc/bls12381/ff/fp12cubic_test.go |  6 ++---
 ecc/bls12381/ff/fp4.go            |  2 ++
 ecc/bls12381/pair.go              |  2 +-
 6 files changed, 56 insertions(+), 6 deletions(-)

diff --git a/ecc/bls12381/ff/cyclo6.go b/ecc/bls12381/ff/cyclo6.go
index 5396cfa98..5eb8a7483 100644
--- a/ecc/bls12381/ff/cyclo6.go
+++ b/ecc/bls12381/ff/cyclo6.go
@@ -8,9 +8,44 @@ func (z Cyclo6) IsEqual(x *Cyclo6) int    { return (Fp12)(z).IsEqual((*Fp12)(x))
 func (z Cyclo6) IsIdentity() int          { i := &Fp12{}; i.SetOne(); return z.IsEqual((*Cyclo6)(i)) }
 func (z *Cyclo6) Frob(x *Cyclo6)          { (*Fp12)(z).Frob((*Fp12)(x)) }
 func (z *Cyclo6) Mul(x, y *Cyclo6)        { (*Fp12)(z).Mul((*Fp12)(x), (*Fp12)(y)) }
-func (z *Cyclo6) Sqr(x *Cyclo6)           { (*Fp12)(z).Sqr((*Fp12)(x)) }
 func (z *Cyclo6) Inv(x *Cyclo6)           { *z = *x; z[1].Neg() }
 func (z *Cyclo6) exp(x *Cyclo6, n []byte) { (*Fp12)(z).Exp((*Fp12)(x), n) }
+func (z *Cyclo6) Sqr(x *Cyclo6) {
+	// Method of Granger-Scott.
+	// Page 7 of "Faster Squaring in the Cyclotomic Subgroup of Sixth Degree Extensions"
+	// https://www.iacr.org/archive/pkc2010/60560212/60560212.pdf
+	var xx, zz Fp12Cubic
+	xx.FromFp12((*Fp12)(x))
+
+	a, b, c := &xx[0], &xx[1], &xx[2]
+	z0, z1, z2 := &zz[0], &zz[1], &zz[2]
+
+	var aa, bb, cc, tt Fp4
+	aa.Sqr(a)
+	bb.Sqr(b)
+	cc.Sqr(c)
+	cc.mulT(&cc)
+
+	z0.Add(&aa, &aa)
+	z0.Add(z0, &aa)
+	tt.Add(a, a)
+	tt.Cjg()
+	z0.Sub(z0, &tt)
+
+	z1.Add(&cc, &cc)
+	z1.Add(z1, &cc)
+	tt.Add(b, b)
+	tt.Cjg()
+	z1.Add(z1, &tt)
+
+	z2.Add(&bb, &bb)
+	z2.Add(z2, &bb)
+	tt.Add(c, c)
+	tt.Cjg()
+	z2.Sub(z2, &tt)
+
+	(*Fp12)(z).FromFp12Cubic(&zz)
+}
 
 // PowToX computes z = x^paramX, where paramX is the parameter of the BLS curve.
 func (z *Cyclo6) PowToX(x *Cyclo6) {
diff --git a/ecc/bls12381/ff/cyclo6_test.go b/ecc/bls12381/ff/cyclo6_test.go
index e9b73b059..221f0c8b0 100644
--- a/ecc/bls12381/ff/cyclo6_test.go
+++ b/ecc/bls12381/ff/cyclo6_test.go
@@ -82,6 +82,19 @@ func TestCyclo6(t *testing.T) {
 			}
 		}
 	})
+	t.Run("sqr_sqrfasr", func(t *testing.T) {
+		var want, got Cyclo6
+		for i := 0; i < testTimes; i++ {
+			x := randomCyclo6(t)
+
+			// Specialized square in cyclotomic vs Generic Square in Fp12
+			got.Sqr(x)
+			(*Fp12)(&want).Sqr((*Fp12)(x))
+			if got.IsEqual(&want) == 0 {
+				test.ReportError(t, got, want, x)
+			}
+		}
+	})
 
 	t.Run("invFp12_vs_invCyclo6", func(t *testing.T) {
 		var want, got Fp12
diff --git a/ecc/bls12381/ff/fp12cubic.go b/ecc/bls12381/ff/fp12cubic.go
index 8bcb43a3d..05125837c 100644
--- a/ecc/bls12381/ff/fp12cubic.go
+++ b/ecc/bls12381/ff/fp12cubic.go
@@ -29,7 +29,7 @@ func (z *Fp12Cubic) FromFp12(x *Fp12) {
 	z[2][1] = x[1][2] // w^5
 }
 
-func (z *Fp12) FromFp12Alt(x *Fp12Cubic) {
+func (z *Fp12) FromFp12Cubic(x *Fp12Cubic) {
 	z[0][0] = x[0][0] // w^0
 	z[1][0] = x[1][0] // w^1
 	z[0][1] = x[2][0] // w^2
diff --git a/ecc/bls12381/ff/fp12cubic_test.go b/ecc/bls12381/ff/fp12cubic_test.go
index 9741e867a..899e3b91b 100644
--- a/ecc/bls12381/ff/fp12cubic_test.go
+++ b/ecc/bls12381/ff/fp12cubic_test.go
@@ -17,7 +17,7 @@ func TestFP12CubicAdd(t *testing.T) {
 		yalt.FromFp12(y)
 		zalt.Add(&xalt, &yalt)
 		z.Add(x, y)
-		zcmp.FromFp12Alt(&zalt)
+		zcmp.FromFp12Cubic(&zalt)
 		if z.IsEqual(&zcmp) == 0 {
 			test.ReportError(t, z, zcmp, x, y)
 		}
@@ -35,7 +35,7 @@ func TestFP12CubicMul(t *testing.T) {
 		yalt.FromFp12(y)
 		zalt.Mul(&xalt, &yalt)
 		z.Mul(x, y)
-		zcmp.FromFp12Alt(&zalt)
+		zcmp.FromFp12Cubic(&zalt)
 		if z.IsEqual(&zcmp) == 0 {
 			test.ReportError(t, z, zcmp, x, y)
 		}
@@ -51,7 +51,7 @@ func TestFP12AltSqr(t *testing.T) {
 		xalt.FromFp12(x)
 		zalt.Sqr(&xalt)
 		z.Sqr(x)
-		zcmp.FromFp12Alt(&zalt)
+		zcmp.FromFp12Cubic(&zalt)
 		if z.IsEqual(&zcmp) == 0 {
 			test.ReportError(t, z, zcmp, x)
 		}
diff --git a/ecc/bls12381/ff/fp4.go b/ecc/bls12381/ff/fp4.go
index 5336b4237..e230b4214 100644
--- a/ecc/bls12381/ff/fp4.go
+++ b/ecc/bls12381/ff/fp4.go
@@ -23,6 +23,8 @@ func (z *Fp4) IsEqual(x *Fp4) int {
 	return z[0].IsEqual(&x[0]) & z[1].IsEqual(&x[1])
 }
 
+func (z *Fp4) Cjg() { z[1].Neg() }
+
 func (z *Fp4) Neg() {
 	z[0].Neg()
 	z[1].Neg()
diff --git a/ecc/bls12381/pair.go b/ecc/bls12381/pair.go
index 3b3d71c8d..ead99d03c 100644
--- a/ecc/bls12381/pair.go
+++ b/ecc/bls12381/pair.go
@@ -32,7 +32,7 @@ func miller(f *ff.Fp12, P *G1, Q *G2) {
 			acc.MulLine(acc, g)
 		}
 	}
-	f.FromFp12Alt(acc)
+	f.FromFp12Cubic(acc)
 	f.Cjg() // inverts f as paramX is negative.
 }