Implement binary EEA inversion for faster BNADD precompile

crypto/bn256/cloudflare/bn256.go

@@ -236,7 +236,7 @@ func (e *G2) Marshal() []byte {
 		e.p = &twistPoint{}
 	}
 
-	e.p.MakeAffine()
+	e.p.MakeAffineConstantTime()
 	ret := make([]byte, numBytes*4)
 	if e.p.IsInfinity() {
 		return ret

crypto/bn256/cloudflare/bn256_test.go

@@ -92,6 +92,73 @@ func TestTripartiteDiffieHellman(t *testing.T) {
 	}
 }
 
+func TestBinaryEAA(t *testing.T) {
+	for i := 0; i < 10000; i++ {
+		_, Ga, err := RandomG1(rand.Reader)
+		if err != nil {
+			t.Fatal(err)
+		}
+		tmpLittleFermat := &gfP{}
+		tmpLittleFermat.InvertConstantTime(&Ga.p.x)
+
+		tmpBinaryEAA := &gfP{}
+		tmpBinaryEAA.InvertVariableTime(&Ga.p.x)
+
+		tmpBinaryEAASelfSet := &gfP{}
+		tmpBinaryEAASelfSet.Set(&Ga.p.x)
+		tmpBinaryEAASelfSet.InvertVariableTime(tmpBinaryEAASelfSet)
+
+		eq := equals(tmpLittleFermat, tmpBinaryEAA)
+		if eq == false {
+			t.Fatalf("results of different inversion do not agree")
+		}
+
+		eq = equals(tmpLittleFermat, tmpBinaryEAASelfSet)
+		if eq == false {
+			t.Fatalf("self-assigned inversion is invalid")
+		}
+	}
+}
+
+func BenchmarkLittleFermatInversion(b *testing.B) {
+	el := gfP{0x0, 0x97816a916871ca8d, 0xb85045b68181585d, 0x30644e72e131a029}
+
+	b.ResetTimer()
+
+	tmp := &gfP{}
+	for i := 0; i < b.N; i++ {
+		tmp.InvertConstantTime(&el)
+	}
+}
+
+func BenchmarkBinaryEEAInversion(b *testing.B) {
+	el := gfP{0x0, 0x97816a916871ca8d, 0xb85045b68181585d, 0x30644e72e131a029}
+
+	b.ResetTimer()
+
+	tmp := &gfP{}
+	for i := 0; i < b.N; i++ {
+		tmp.InvertVariableTime(&el)
+	}
+}
+
+func BenchmarkG1AddAndMakeAffine(b *testing.B) {
+	_, Ga, err := RandomG1(rand.Reader)
+	if err != nil {
+		b.Fatal(err)
+	}
+	_, Gb, err := RandomG1(rand.Reader)
+	if err != nil {
+		b.Fatal(err)
+	}
+	b.ResetTimer()
+
+	for i := 0; i < b.N; i++ {
+		e := new(G1).Add(Ga, Gb)
+		e.p.MakeAffine()
+	}
+}
+
 func BenchmarkG1(b *testing.B) {
 	x, _ := rand.Int(rand.Reader, Order)
 	b.ResetTimer()

crypto/bn256/cloudflare/curve.go

@@ -217,7 +217,31 @@ func (c *curvePoint) MakeAffine() {
 	}
 
 	zInv := &gfP{}
-	zInv.Invert(&c.z)
+	zInv.InvertVariableTime(&c.z)
+
+	t, zInv2 := &gfP{}, &gfP{}
+	gfpMul(t, &c.y, zInv)
+	gfpMul(zInv2, zInv, zInv)
+
+	gfpMul(&c.x, &c.x, zInv2)
+	gfpMul(&c.y, t, zInv2)
+
+	c.z = *newGFp(1)
+	c.t = *newGFp(1)
+}
+
+func (c *curvePoint) MakeAffineConstantTime() {
+	if c.z == *newGFp(1) {
+		return
+	} else if c.z == *newGFp(0) {
+		c.x = gfP{0}
+		c.y = *newGFp(1)
+		c.t = gfP{0}
+		return
+	}
+
+	zInv := &gfP{}
+	zInv.InvertConstantTime(&c.z)
 
 	t, zInv2 := &gfP{}, &gfP{}
 	gfpMul(t, &c.y, zInv)

crypto/bn256/cloudflare/gfp.go

@@ -3,6 +3,7 @@ package bn256
 import (
 	"errors"
 	"fmt"
+	"math/bits"
 )
 
 type gfP [4]uint64
@@ -30,7 +31,7 @@ func (e *gfP) Set(f *gfP) {
 	e[3] = f[3]
 }
 
-func (e *gfP) Invert(f *gfP) {
+func (e *gfP) InvertConstantTime(f *gfP) {
 	bits := [4]uint64{0x3c208c16d87cfd45, 0x97816a916871ca8d, 0xb85045b68181585d, 0x30644e72e131a029}
 
 	sum, power := &gfP{}, &gfP{}
@@ -79,3 +80,130 @@ func (e *gfP) Unmarshal(in []byte) error {
 
 func montEncode(c, a *gfP) { gfpMul(c, a, r2) }
 func montDecode(c, a *gfP) { gfpMul(c, a, &gfP{1}) }
+
+func isZero(a *gfP) bool {
+	return a[0] == 0 && a[1] == 0 && a[2] == 0 && a[3] == 0
+}
+
+func isEven(a *gfP) bool {
+	return bits.TrailingZeros64((a[0])) > 0
+}
+
+func div2(a *gfP) {
+	a[0] = a[0]>>1 | a[1]<<63
+	a[1] = a[1]>>1 | a[2]<<63
+	a[2] = a[2]>>1 | a[3]<<63
+	a[3] = a[3] >> 1
+}
+
+func (e *gfP) addNocarry(f *gfP) {
+	carry := uint64(0)
+	e[0], carry = bits.Add64(e[0], f[0], carry)
+	e[1], carry = bits.Add64(e[1], f[1], carry)
+	e[2], carry = bits.Add64(e[2], f[2], carry)
+	e[3], _ = bits.Add64(e[3], f[3], carry)
+}
+
+func (e *gfP) subNoborrow(f *gfP) {
+	borrow := uint64(0)
+	e[0], borrow = bits.Sub64(e[0], f[0], borrow)
+	e[1], borrow = bits.Sub64(e[1], f[1], borrow)
+	e[2], borrow = bits.Sub64(e[2], f[2], borrow)
+	e[3], _ = bits.Sub64(e[3], f[3], borrow)
+}
+
+func gte(a, b *gfP) bool {
+	// subtract b from a. If no borrow occures then a >= b
+	borrow := uint64(0)
+	_, borrow = bits.Sub64(a[0], b[0], borrow)
+	_, borrow = bits.Sub64(a[1], b[1], borrow)
+	_, borrow = bits.Sub64(a[2], b[2], borrow)
+	_, borrow = bits.Sub64(a[3], b[3], borrow)
+
+	return borrow == 0
+}
+
+func equals(a, b *gfP) bool {
+	return a[0] == b[0] && a[1] == b[1] && a[2] == b[2] && a[3] == b[3]
+}
+
+// Performs inversion of the field element using binary EEA.
+// If element is zero (no inverse exists) then set `e` to zero
+func (e *gfP) InvertVariableTime(f *gfP) {
+	if isZero(f) {
+		e.Set(&gfP{0, 0, 0, 0})
+		return
+	}
+
+	// Guajardo Kumar Paar Pelzl
+	// Efficient Software-Implementation of Finite Fields with Applications to Cryptography
+	// Algorithm 16 (BEA for Inversion in Fp)
+
+	one := gfP{1, 0, 0, 0}
+
+	u, b := gfP{}, gfP{}
+	u.Set(f)
+	b.Set(r2)
+
+	v := gfP{p2[0], p2[1], p2[2], p2[3]}
+	c := gfP{0, 0, 0, 0}
+	modulus := gfP{p2[0], p2[1], p2[2], p2[3]}
+
+	for {
+		if equals(&u, &one) || equals(&v, &one) {
+			break
+		}
+
+		// while u is even
+		for {
+			if !isEven(&u) {
+				break
+			}
+
+			div2(&u)
+			if isEven(&b) {
+				div2(&b)
+			} else {
+				// we will not overflow a modulus here,
+				// so we can use specialized function
+				// do perform addition without reduction
+				b.addNocarry(&modulus)
+				div2(&b)
+			}
+		}
+
+		// while v is even
+		for {
+			if !isEven(&v) {
+				break
+			}
+
+			div2(&v)
+			if isEven(&c) {
+				div2(&c)
+			} else {
+				// we will not overflow a modulus here,
+				// so we can use specialized function
+				// do perform addition without reduction
+				c.addNocarry(&modulus)
+				div2(&c)
+			}
+		}
+
+		if gte(&v, &u) {
+			// v >= u
+			v.subNoborrow(&u)
+			gfpSub(&c, &c, &b)
+		} else {
+			// if v < u
+			u.subNoborrow(&v)
+			gfpSub(&b, &b, &c)
+		}
+	}
+
+	if equals(&u, &one) {
+		e.Set(&b)
+	} else {
+		e.Set(&c)
+	}
+}

crypto/bn256/cloudflare/gfp12.go

@@ -143,15 +143,31 @@ func (e *gfP12) Square(a *gfP12) *gfP12 {
 	return e
 }
 
-func (e *gfP12) Invert(a *gfP12) *gfP12 {
+func (e *gfP12) InvertVariableTime(a *gfP12) *gfP12 {
 	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
 	// ftp://136.206.11.249/pub/crypto/pairings.pdf
 	t1, t2 := &gfP6{}, &gfP6{}
 
 	t1.Square(&a.x)
 	t2.Square(&a.y)
 	t1.MulTau(t1).Sub(t2, t1)
-	t2.Invert(t1)
+	t2.InvertVariableTime(t1)
+
+	e.x.Neg(&a.x)
+	e.y.Set(&a.y)
+	e.MulScalar(e, t2)
+	return e
+}
+
+func (e *gfP12) InvertConstantTime(a *gfP12) *gfP12 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+	t1, t2 := &gfP6{}, &gfP6{}
+
+	t1.Square(&a.x)
+	t2.Square(&a.y)
+	t1.MulTau(t1).Sub(t2, t1)
+	t2.InvertConstantTime(t1)
 
 	e.x.Neg(&a.x)
 	e.y.Set(&a.y)

crypto/bn256/cloudflare/gfp2.go

@@ -137,7 +137,7 @@ func (e *gfP2) Square(a *gfP2) *gfP2 {
 	return e
 }
 
-func (e *gfP2) Invert(a *gfP2) *gfP2 {
+func (e *gfP2) InvertVariableTime(a *gfP2) *gfP2 {
 	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
 	// ftp://136.206.11.249/pub/crypto/pairings.pdf
 	t1, t2 := &gfP{}, &gfP{}
@@ -146,7 +146,25 @@ func (e *gfP2) Invert(a *gfP2) *gfP2 {
 	gfpAdd(t1, t1, t2)
 
 	inv := &gfP{}
-	inv.Invert(t1)
+	inv.InvertVariableTime(t1)
+
+	gfpNeg(t1, &a.x)
+
+	gfpMul(&e.x, t1, inv)
+	gfpMul(&e.y, &a.y, inv)
+	return e
+}
+
+func (e *gfP2) InvertConstantTime(a *gfP2) *gfP2 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+	t1, t2 := &gfP{}, &gfP{}
+	gfpMul(t1, &a.x, &a.x)
+	gfpMul(t2, &a.y, &a.y)
+	gfpAdd(t1, t1, t2)
+
+	inv := &gfP{}
+	inv.InvertConstantTime(t1)
 
 	gfpNeg(t1, &a.x)
 

crypto/bn256/cloudflare/gfp6.go

@@ -168,7 +168,7 @@ func (e *gfP6) Square(a *gfP6) *gfP6 {
 	return e
 }
 
-func (e *gfP6) Invert(a *gfP6) *gfP6 {
+func (e *gfP6) InvertConstantTime(a *gfP6) *gfP6 {
 	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
 	// ftp://136.206.11.249/pub/crypto/pairings.pdf
 
@@ -204,7 +204,51 @@ func (e *gfP6) Invert(a *gfP6) *gfP6 {
 	t1.Mul(B, &a.x).MulXi(t1)
 	F.Add(F, t1)
 
-	F.Invert(F)
+	F.InvertVariableTime(F)
+
+	e.x.Mul(C, F)
+	e.y.Mul(B, F)
+	e.z.Mul(A, F)
+	return e
+}
+
+func (e *gfP6) InvertVariableTime(a *gfP6) *gfP6 {
+	// See "Implementing cryptographic pairings", M. Scott, section 3.2.
+	// ftp://136.206.11.249/pub/crypto/pairings.pdf
+
+	// Here we can give a short explanation of how it works: let j be a cubic root of
+	// unity in GF(p²) so that 1+j+j²=0.
+	// Then (xτ² + yτ + z)(xj²τ² + yjτ + z)(xjτ² + yj²τ + z)
+	// = (xτ² + yτ + z)(Cτ²+Bτ+A)
+	// = (x³ξ²+y³ξ+z³-3ξxyz) = F is an element of the base field (the norm).
+	//
+	// On the other hand (xj²τ² + yjτ + z)(xjτ² + yj²τ + z)
+	// = τ²(y²-ξxz) + τ(ξx²-yz) + (z²-ξxy)
+	//
+	// So that's why A = (z²-ξxy), B = (ξx²-yz), C = (y²-ξxz)
+	t1 := (&gfP2{}).Mul(&a.x, &a.y)
+	t1.MulXi(t1)
+
+	A := (&gfP2{}).Square(&a.z)
+	A.Sub(A, t1)
+
+	B := (&gfP2{}).Square(&a.x)
+	B.MulXi(B)
+	t1.Mul(&a.y, &a.z)
+	B.Sub(B, t1)
+
+	C := (&gfP2{}).Square(&a.y)
+	t1.Mul(&a.x, &a.z)
+	C.Sub(C, t1)
+
+	F := (&gfP2{}).Mul(C, &a.y)
+	F.MulXi(F)
+	t1.Mul(A, &a.z)
+	F.Add(F, t1)
+	t1.Mul(B, &a.x).MulXi(t1)
+	F.Add(F, t1)
+
+	F.InvertConstantTime(F)
 
 	e.x.Mul(C, F)
 	e.y.Mul(B, F)