From b380bca1d2bb82fe0ba68195980e00fcce6d0224 Mon Sep 17 00:00:00 2001
From: Evgenii <evgeniy@nspcc.ru>
Date: Tue, 29 Jan 2019 14:11:08 +0300
Subject: [PATCH] add inversion function

---
 gf127/gf127.go      | 70 ++++++++++++++++++++++++++++++++++++++++++---
 gf127/gf127_test.go | 20 +++++++++++--
 2 files changed, 83 insertions(+), 7 deletions(-)

diff --git a/gf127/gf127.go b/gf127/gf127.go
index 5d2cf0a..1dadb3b 100644
--- a/gf127/gf127.go
+++ b/gf127/gf127.go
@@ -11,16 +11,24 @@ import (
 	"encoding/binary"
 	"encoding/hex"
 	"errors"
+	"math/bits"
 )
 
 // GF127 represents element of GF(2^127)
 type GF127 [2]uint64
 
-const msb64 = 0x8000000000000000
-const byteSize = 16
+const (
+	msb64    = 0x8000000000000000
+	byteSize = 16
+)
 
-// x127x64 represents x^127 + x^63. Used in assembly file.
-var x127x63 = GF127{msb64, msb64}
+var (
+	// x127x63 represents x^127 + x^63. Used in assembly file.
+	x127x63 = GF127{msb64, msb64}
+
+	// x126x631 is reduction polynomial x^127+x^63+1
+	x127x631 = GF127{msb64 + 1, msb64}
+)
 
 // New constructs new element of GF(2^127) as hi*x^64 + lo.
 // It is assumed that hi has zero MSB.
@@ -66,6 +74,60 @@ func (c *GF127) UnmarshalBinary(data []byte) error {
 	return nil
 }
 
+// Inv sets b to a^-1
+// Algorithm is based on Extended Euclidean Algorithm
+// and is described by Hankerson, Hernandez, Menezes in
+// https://link.springer.com/content/pdf/10.1007/3-540-44499-8_1.pdf
+func Inv(a, b *GF127) {
+	var (
+		v    = x127x631
+		u    = *a
+		c, d = &GF127{1, 0}, &GF127{0, 0}
+		t    = new(GF127)
+		x    *GF127
+	)
+
+	// degree of polynomial is a position of most significant bit
+	for du, dv := msb(&u), msb(&v); du != 0; du, dv = msb(&u), msb(&v) {
+		if du < dv {
+			v, u = u, v
+			dv, du = du, dv
+			d, c = c, d
+		}
+
+		x = xN(du - dv)
+
+		Mul(x, &v, t)
+		Add(&u, t, &u)
+
+		// becasuse mul performs reduction on t, we need
+		// manually reduce u at first step
+		if msb(&u) == 127 {
+			Add(&u, &x127x63, &u)
+			Add(&u, &GF127{1, 0}, &u)
+		}
+
+		Mul(x, d, t)
+		Add(c, t, c)
+	}
+	*b = *c
+}
+
+func xN(n int) *GF127 {
+	if n < 64 {
+		return &GF127{1 << uint(n), 0}
+	}
+	return &GF127{0, 1 << (uint(n) >> 8)}
+}
+
+func msb(a *GF127) (x int) {
+	x = bits.LeadingZeros64(a[1])
+	if x == 64 {
+		x = bits.LeadingZeros64(a[0]) + 64
+	}
+	return 127 - x
+}
+
 // Mul sets c to the product a*b and returns c.
 func (c *GF127) Mul(a, b *GF127) *GF127 {
 	Mul(a, b, c)
diff --git a/gf127/gf127_test.go b/gf127/gf127_test.go
index 1304e6e..6d0c7ae 100644
--- a/gf127/gf127_test.go
+++ b/gf127/gf127_test.go
@@ -30,7 +30,7 @@ var testCasesMul = [][3]*GF127{
 }
 
 func TestMul(t *testing.T) {
-	c := &GF127{0, 0}
+	c := new(GF127)
 	for _, tc := range testCasesMul {
 		if Mul(tc[0], tc[1], c); !c.Equals(tc[2]) {
 			t.Errorf("expected (%s), got (%s)", c.String(), tc[2].String())
@@ -45,7 +45,7 @@ var testCasesMul10 = [][2]*GF127{
 }
 
 func TestMul10(t *testing.T) {
-	c := &GF127{0, 0}
+	c := new(GF127)
 	for _, tc := range testCasesMul10 {
 		if Mul10(tc[0], c); !c.Equals(tc[1]) {
 			t.Errorf("expected (%s), got (%s)", tc[1].String(), c.String())
@@ -60,10 +60,24 @@ var testCasesMul11 = [][2]*GF127{
 }
 
 func TestMul11(t *testing.T) {
-	c := &GF127{0, 0}
+	c := new(GF127)
 	for _, tc := range testCasesMul11 {
 		if Mul11(tc[0], c); !c.Equals(tc[1]) {
 			t.Errorf("expected (%s), got (%s)", tc[1].String(), c.String())
 		}
 	}
 }
+
+var testCasesInv = [][2]*GF127{
+	{&GF127{1, 0}, &GF127{1, 0}},
+	{&GF127{54321, 12345}, &GF127{8230555108620784737, 3929873967650665114}},
+}
+
+func TestInv(t *testing.T) {
+	c := new(GF127)
+	for _, tc := range testCasesInv {
+		if Inv(tc[0], c); !c.Equals(tc[1]) {
+			t.Errorf("expected (%s), got (%s)", tc[1].String(), c.String())
+		}
+	}
+}