From 946b85b409d3dfc1bc9dc6125bf270cfc5cc9a87 Mon Sep 17 00:00:00 2001
From: Sun Yimin <emmansun@users.noreply.github.com>
Date: Fri, 28 Apr 2023 15:31:59 +0800
Subject: [PATCH] sm9: reduce mul, improve performance

---
 sm9/bn256/gfp2.go      | 16 +++++++++++++++-
 sm9/bn256/gfp4.go      | 38 +++++++++++++++++++++++---------------
 sm9/bn256/gfp4_test.go | 29 +++++++++++++++++++++++++++++
 3 files changed, 67 insertions(+), 16 deletions(-)

diff --git a/sm9/bn256/gfp2.go b/sm9/bn256/gfp2.go
index 872c0e7..17198cb 100644
--- a/sm9/bn256/gfp2.go
+++ b/sm9/bn256/gfp2.go
@@ -146,7 +146,7 @@ func (e *gfP2) MulU(a, b *gfP2) *gfP2 {
 	gfpSub(ty, ty, v1)
 	gfpAdd(ty, ty, ty)
 	gfpNeg(ty, ty)
-	
+
 	gfpSub(tx, v0, v1)
 	gfpSub(tx, tx, v1)
 
@@ -155,6 +155,20 @@ func (e *gfP2) MulU(a, b *gfP2) *gfP2 {
 	return e
 }
 
+// MulU1: a  * u
+//(a0+a1*u)u=c0+c1*u, where
+//c1 = a0
+//c0 = -2a1
+func (e *gfP2) MulU1(a *gfP2) *gfP2 {
+	t := &gfP{}
+	gfpAdd(t, &a.x, &a.x)
+	gfpNeg(t, t)
+
+	e.x.Set(&a.y)
+	e.y.Set(t)
+	return e
+}
+
 func (e *gfP2) Square(a *gfP2) *gfP2 {
 	// Complex squaring algorithm:
 	// (xu+y)² = y^2-2*x^2 + 2*u*x*y
diff --git a/sm9/bn256/gfp4.go b/sm9/bn256/gfp4.go
index 763159b..96c3c57 100644
--- a/sm9/bn256/gfp4.go
+++ b/sm9/bn256/gfp4.go
@@ -99,15 +99,18 @@ func (e *gfP4) Mul(a, b *gfP4) *gfP4 {
 	//(a0+a1*v)(b0+b1*v)=c0+c1*v, where
 	//c0 = a0*b0 +a1*b1*u
 	//c1 = (a0 + a1)(b0 + b1) - a0*b0 - a1*b1 = a0*b1 + a1*b0
-	tx, t := &gfP2{}, &gfP2{}
-	tx.Mul(&a.x, &b.y)
-	t.Mul(&a.y, &b.x)
-	tx.Add(tx, t)
+	tx, ty, v0, v1 := &gfP2{}, &gfP2{}, &gfP2{}, &gfP2{}
+	v0.Mul(&a.y, &b.y)
+	v1.Mul(&a.x, &b.x)
 
-	ty := &gfP2{}
-	ty.Mul(&a.y, &b.y)
-	t.MulU(&a.x, &b.x)
-	ty.Add(ty, t)
+	tx.Add(&a.x, &a.y)
+	ty.Add(&b.x, &b.y)
+	tx.Mul(tx, ty)
+	tx.Sub(tx, v0)
+	tx.Sub(tx, v1)
+
+	ty.MulU1(v1)
+	ty.Add(ty, v0)
 
 	e.x.Set(tx)
 	e.y.Set(ty)
@@ -121,14 +124,19 @@ func (e *gfP4) Mul(a, b *gfP4) *gfP4 {
 // c0 = a0*b1*u + a1*b0*u
 // c1 = a0*b0 + a1*b1*u
 func (e *gfP4) MulV(a, b *gfP4) *gfP4 {
-	tx, ty, t := &gfP2{}, &gfP2{}, &gfP2{}
-	ty.MulU(&a.y, &b.x)
-	t.MulU(&a.x, &b.y)
-	ty.Add(ty, t)
+	tx, ty, v0, v1 := &gfP2{}, &gfP2{}, &gfP2{}, &gfP2{}
+	v0.Mul(&a.y, &b.y)
+	v1.Mul(&a.x, &b.x)
 
-	tx.Mul(&a.y, &b.y)
-	t.MulU(&a.x, &b.x)
-	tx.Add(tx, t)
+	tx.Add(&a.x, &a.y)
+	ty.Add(&b.x, &b.y)
+	ty.Mul(tx, ty)
+	ty.Sub(ty, v0)
+	ty.Sub(ty, v1)
+	ty.MulU1(ty)
+
+	tx.MulU1(v1)
+	tx.Add(tx, v0)
 
 	e.x.Set(tx)
 	e.y.Set(ty)
diff --git a/sm9/bn256/gfp4_test.go b/sm9/bn256/gfp4_test.go
index 1ab525b..a875f54 100644
--- a/sm9/bn256/gfp4_test.go
+++ b/sm9/bn256/gfp4_test.go
@@ -178,3 +178,32 @@ func Test_gfP4FrobeniusP3(t *testing.T) {
 		t.Errorf("got %v, expected %v", got, expected)
 	}
 }
+
+func BenchmarkGfP4Mul(b *testing.B) {
+	x := &gfP4{
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+	}
+	y := &gfP4{
+		gfP2{
+			*fromBigInt(bigFromHex("17509B092E845C1266BA0D262CBEE6ED0736A96FA347C8BD856DC76B84EBEB96")),
+			*fromBigInt(bigFromHex("A7CF28D519BE3DA65F3170153D278FF247EFBA98A71A08116215BBA5C999A7C7")),
+		},
+		gfP2{
+			*fromBigInt(bigFromHex("85AEF3D078640C98597B6027B441A01FF1DD2C190F5E93C454806C11D8806141")),
+			*fromBigInt(bigFromHex("3722755292130B08D2AAB97FD34EC120EE265948D19C17ABF9B7213BAF82D65B")),
+		},
+	}
+	b.ReportAllocs()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		t := &gfP4{}
+		t.Mul(x, y)
+	}
+}