Update to goff v0.2.0

Fix compat with 32 bit arch

.travis.yml

@@ -4,5 +4,12 @@ language: go
 go:
   - "1.12"
 
+jobs:
+  include:
+    - name: "Unit Tests 64 bit arch"
+      env: GOARCH="amd64"
+    - name: "Unit Test 32 bit arch"
+      env: GOARCH="386"
+
 env:
   - GO111MODULE=on

ff/arith.go

@@ -12,14 +12,19 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-// Code generated by goff DO NOT EDIT
+// Code generated by goff (v0.2.0) DO NOT EDIT
 
+// Package ff contains field arithmetic operations
 package ff
 
 import (
 	"math/bits"
+
+	"golang.org/x/sys/cpu"
 )
 
+var supportAdx = cpu.X86.HasADX && cpu.X86.HasBMI2
+
 func madd(a, b, t, u, v uint64) (uint64, uint64, uint64) {
 	var carry uint64
 	hi, lo := bits.Mul64(a, b)

ff/element.go

@@ -12,29 +12,33 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
-// field modulus q =
-//
-// 21888242871839275222246405745257275088548364400416034343698204186575808495617
-// Code generated by goff DO NOT EDIT
-// goff version:  - build:
-// Element are assumed to be in Montgomery form in all methods
+// Code generated by goff (v0.2.0) DO NOT EDIT
 
-// Package ff (generated by goff) contains field arithmetics operations
+// Package ff contains field arithmetic operations
 package ff
 
+// /!\ WARNING /!\
+// this code has not been audited and is provided as-is. In particular,
+// there is no security guarantees such as constant time implementation
+// or side-channel attack resistance
+// /!\ WARNING /!\
+
 import (
 	"crypto/rand"
 	"encoding/binary"
 	"io"
 	"math/big"
 	"math/bits"
+	"strconv"
 	"sync"
-
 	"unsafe"
 )
 
 // Element represents a field element stored on 4 words (uint64)
 // Element are assumed to be in Montgomery form in all methods
+// field modulus q =
+//
+// 21888242871839275222246405745257275088548364400416034343698204186575808495617
 type Element [4]uint64
 
 // ElementLimbs number of 64 bits words needed to represent Element
@@ -311,6 +315,7 @@ func (z *Element) SetRandom() *Element {
 	z[3] %= 3486998266802970665
 
 	// if z > q --> z -= q
+	// note: this is NOT constant time
 	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
 		var b uint64
 		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
@@ -322,6 +327,38 @@ func (z *Element) SetRandom() *Element {
 	return z
 }
 
+// One returns 1 (in montgommery form)
+func One() Element {
+	var one Element
+	one.SetOne()
+	return one
+}
+
+// FromInterface converts i1 from uint64, int, string, or Element, big.Int into Element
+// panic if provided type is not supported
+func FromInterface(i1 interface{}) Element {
+	var val Element
+
+	switch c1 := i1.(type) {
+	case uint64:
+		val.SetUint64(c1)
+	case int:
+		val.SetString(strconv.Itoa(c1))
+	case string:
+		val.SetString(c1)
+	case big.Int:
+		val.SetBigInt(&c1)
+	case Element:
+		val = c1
+	case *Element:
+		val.Set(c1)
+	default:
+		panic("invalid type")
+	}
+
+	return val
+}
+
 // Add z = x + y mod q
 func (z *Element) Add(x, y *Element) *Element {
 	var carry uint64
@@ -332,6 +369,7 @@ func (z *Element) Add(x, y *Element) *Element {
 	z[3], _ = bits.Add64(x[3], y[3], carry)
 
 	// if z > q --> z -= q
+	// note: this is NOT constant time
 	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
 		var b uint64
 		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
@@ -352,6 +390,7 @@ func (z *Element) AddAssign(x *Element) *Element {
 	z[3], _ = bits.Add64(z[3], x[3], carry)
 
 	// if z > q --> z -= q
+	// note: this is NOT constant time
 	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
 		var b uint64
 		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
@@ -372,6 +411,7 @@ func (z *Element) Double(x *Element) *Element {
 	z[3], _ = bits.Add64(x[3], x[3], carry)
 
 	// if z > q --> z -= q
+	// note: this is NOT constant time
 	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
 		var b uint64
 		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
@@ -416,18 +456,31 @@ func (z *Element) SubAssign(x *Element) *Element {
 	return z
 }
 
-// Exp z = x^e mod q
-func (z *Element) Exp(x Element, e uint64) *Element {
-	if e == 0 {
+// Exp z = x^exponent mod q
+// (not optimized)
+// exponent (non-montgomery form) is ordered from least significant word to most significant word
+func (z *Element) Exp(x Element, exponent ...uint64) *Element {
+	r := 0
+	msb := 0
+	for i := len(exponent) - 1; i >= 0; i-- {
+		if exponent[i] == 0 {
+			r++
+		} else {
+			msb = (i * 64) + bits.Len64(exponent[i])
+			break
+		}
+	}
+	exponent = exponent[:len(exponent)-r]
+	if len(exponent) == 0 {
 		return z.SetOne()
 	}
 
 	z.Set(&x)
 
-	l := bits.Len64(e) - 2
+	l := msb - 2
 	for i := l; i >= 0; i-- {
 		z.Square(z)
-		if e&(1<<uint(i)) != 0 {
+		if exponent[i/64]&(1<<uint(i%64)) != 0 {
 			z.MulAssign(&x)
 		}
 	}
@@ -478,6 +531,7 @@ func (z *Element) FromMont() *Element {
 	}
 
 	// if z > q --> z -= q
+	// note: this is NOT constant time
 	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
 		var b uint64
 		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
@@ -513,15 +567,33 @@ func (z *Element) String() string {
 
 // ToBigInt returns z as a big.Int in Montgomery form
 func (z *Element) ToBigInt(res *big.Int) *big.Int {
-	bits := (*[4]big.Word)(unsafe.Pointer(z))
-	return res.SetBits(bits[:])
+	if bits.UintSize == 64 {
+		bits := (*[4]big.Word)(unsafe.Pointer(z))
+		return res.SetBits(bits[:])
+	} else {
+		var bits [8]big.Word
+		for i := 0; i < len(z); i++ {
+			bits[i*2] = big.Word(z[i])
+			bits[i*2+1] = big.Word(z[i] >> 32)
+		}
+		return res.SetBits(bits[:])
+	}
 }
 
 // ToBigIntRegular returns z as a big.Int in regular form
 func (z Element) ToBigIntRegular(res *big.Int) *big.Int {
 	z.FromMont()
-	bits := (*[4]big.Word)(unsafe.Pointer(&z))
-	return res.SetBits(bits[:])
+	if bits.UintSize == 64 {
+		bits := (*[4]big.Word)(unsafe.Pointer(&z))
+		return res.SetBits(bits[:])
+	} else {
+		var bits [8]big.Word
+		for i := 0; i < len(z); i++ {
+			bits[i*2] = big.Word(z[i])
+			bits[i*2+1] = big.Word(z[i] >> 32)
+		}
+		return res.SetBits(bits[:])
+	}
 }
 
 // SetBigInt sets z to v (regular form) and returns z in Montgomery form
@@ -531,6 +603,19 @@ func (z *Element) SetBigInt(v *big.Int) *Element {
 	zero := big.NewInt(0)
 	q := elementModulusBigInt()
 
+	// fast path
+	c := v.Cmp(q)
+	if c == 0 {
+		return z
+	} else if c != 1 && v.Cmp(zero) != -1 {
+		// v should
+		vBits := v.Bits()
+		for i := 0; i < len(vBits); i++ {
+			z[i] = uint64(vBits[i])
+		}
+		return z.ToMont()
+	}
+
 	// copy input
 	vv := new(big.Int).Set(v)
 
@@ -548,8 +633,18 @@ func (z *Element) SetBigInt(v *big.Int) *Element {
 	}
 	// v should
 	vBits := vv.Bits()
-	for i := 0; i < len(vBits); i++ {
-		z[i] = uint64(vBits[i])
+	if bits.UintSize == 64 {
+		for i := 0; i < len(vBits); i++ {
+			z[i] = uint64(vBits[i])
+		}
+	} else {
+		for i := 0; i < len(vBits); i++ {
+			if i%2 == 0 {
+				z[i/2] = uint64(vBits[i])
+			} else {
+				z[i/2] |= uint64(vBits[i]) << 32
+			}
+		}
 	}
 	return z.ToMont()
 }
@@ -563,202 +658,97 @@ func (z *Element) SetString(s string) *Element {
 	return z.SetBigInt(x)
 }
 
-// Mul z = x * y mod q
-func (z *Element) Mul(x, y *Element) *Element {
+// Legendre returns the Legendre symbol of z (either +1, -1, or 0.)
+func (z *Element) Legendre() int {
+	var l Element
+	// z^((q-1)/2)
+	l.Exp(*z,
+		11669102379873075200,
+		10671829228508198984,
+		15863968012492123182,
+		1743499133401485332,
+	)
 
-	var t [4]uint64
-	var c [3]uint64
-	{
-		// round 0
-		v := x[0]
-		c[1], c[0] = bits.Mul64(v, y[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd1(v, y[1], c[1])
-		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd1(v, y[2], c[1])
-		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd1(v, y[3], c[1])
-		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
-	}
-	{
-		// round 1
-		v := x[1]
-		c[1], c[0] = madd1(v, y[0], t[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd2(v, y[1], c[1], t[1])
-		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd2(v, y[2], c[1], t[2])
-		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd2(v, y[3], c[1], t[3])
-		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
-	}
-	{
-		// round 2
-		v := x[2]
-		c[1], c[0] = madd1(v, y[0], t[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd2(v, y[1], c[1], t[1])
-		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd2(v, y[2], c[1], t[2])
-		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd2(v, y[3], c[1], t[3])
-		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
-	}
-	{
-		// round 3
-		v := x[3]
-		c[1], c[0] = madd1(v, y[0], t[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd2(v, y[1], c[1], t[1])
-		c[2], z[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd2(v, y[2], c[1], t[2])
-		c[2], z[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd2(v, y[3], c[1], t[3])
-		z[3], z[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	if l.IsZero() {
+		return 0
 	}
 
-	// if z > q --> z -= q
-	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
-		var b uint64
-		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
-		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
-		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
-		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
+	// if l == 1
+	if (l[3] == 1011752739694698287) && (l[2] == 7381016538464732718) && (l[1] == 3962172157175319849) && (l[0] == 12436184717236109307) {
+		return 1
 	}
-	return z
+	return -1
 }
 
-// MulAssign z = z * x mod q
-func (z *Element) MulAssign(x *Element) *Element {
+// Sqrt z = √x mod q
+// if the square root doesn't exist (x is not a square mod q)
+// Sqrt leaves z unchanged and returns nil
+func (z *Element) Sqrt(x *Element) *Element {
+	// q ≡ 1 (mod 4)
+	// see modSqrtTonelliShanks in math/big/int.go
+	// using https://www.maa.org/sites/default/files/pdf/upload_library/22/Polya/07468342.di020786.02p0470a.pdf
 
-	var t [4]uint64
-	var c [3]uint64
-	{
-		// round 0
-		v := z[0]
-		c[1], c[0] = bits.Mul64(v, x[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd1(v, x[1], c[1])
-		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd1(v, x[2], c[1])
-		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd1(v, x[3], c[1])
-		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
-	}
-	{
-		// round 1
-		v := z[1]
-		c[1], c[0] = madd1(v, x[0], t[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd2(v, x[1], c[1], t[1])
-		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd2(v, x[2], c[1], t[2])
-		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd2(v, x[3], c[1], t[3])
-		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
-	}
-	{
-		// round 2
-		v := z[2]
-		c[1], c[0] = madd1(v, x[0], t[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd2(v, x[1], c[1], t[1])
-		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd2(v, x[2], c[1], t[2])
-		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd2(v, x[3], c[1], t[3])
-		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
-	}
-	{
-		// round 3
-		v := z[3]
-		c[1], c[0] = madd1(v, x[0], t[0])
-		m := c[0] * 14042775128853446655
-		c[2] = madd0(m, 4891460686036598785, c[0])
-		c[1], c[0] = madd2(v, x[1], c[1], t[1])
-		c[2], z[0] = madd2(m, 2896914383306846353, c[2], c[0])
-		c[1], c[0] = madd2(v, x[2], c[1], t[2])
-		c[2], z[1] = madd2(m, 13281191951274694749, c[2], c[0])
-		c[1], c[0] = madd2(v, x[3], c[1], t[3])
-		z[3], z[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
-	}
+	var y, b, t, w Element
+	// w = x^((s-1)/2))
+	w.Exp(*x,
+		14829091926808964255,
+		867720185306366531,
+		688207751544974772,
+		6495040407,
+	)
 
-	// if z > q --> z -= q
-	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
-		var b uint64
-		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
-		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
-		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
-		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
+	// y = x^((s+1)/2)) = w * x
+	y.Mul(x, &w)
+
+	// b = x^s = w * w * x = y * x
+	b.Mul(&w, &y)
+
+	// g = nonResidue ^ s
+	var g = Element{
+		7164790868263648668,
+		11685701338293206998,
+		6216421865291908056,
+		1756667274303109607,
+	}
+	r := uint64(28)
+
+	// compute legendre symbol
+	// t = x^((q-1)/2) = r-1 squaring of x^s
+	t = b
+	for i := uint64(0); i < r-1; i++ {
+		t.Square(&t)
+	}
+	if t.IsZero() {
+		return z.SetZero()
+	}
+	if !((t[3] == 1011752739694698287) && (t[2] == 7381016538464732718) && (t[1] == 3962172157175319849) && (t[0] == 12436184717236109307)) {
+		// t != 1, we don't have a square root
+		return nil
+	}
+	for {
+		var m uint64
+		t = b
+
+		// for t != 1
+		for !((t[3] == 1011752739694698287) && (t[2] == 7381016538464732718) && (t[1] == 3962172157175319849) && (t[0] == 12436184717236109307)) {
+			t.Square(&t)
+			m++
+		}
+
+		if m == 0 {
+			return z.Set(&y)
+		}
+		// t = g^(2^(r-m-1)) mod q
+		ge := int(r - m - 1)
+		t = g
+		for ge > 0 {
+			t.Square(&t)
+			ge--
+		}
+
+		g.Square(&t)
+		y.MulAssign(&t)
+		b.MulAssign(&g)
+		r = m
 	}
-	return z
-}
-
-// Square z = x * x mod q
-func (z *Element) Square(x *Element) *Element {
-
-	var p [4]uint64
-
-	var u, v uint64
-	{
-		// round 0
-		u, p[0] = bits.Mul64(x[0], x[0])
-		m := p[0] * 14042775128853446655
-		C := madd0(m, 4891460686036598785, p[0])
-		var t uint64
-		t, u, v = madd1sb(x[0], x[1], u)
-		C, p[0] = madd2(m, 2896914383306846353, v, C)
-		t, u, v = madd1s(x[0], x[2], t, u)
-		C, p[1] = madd2(m, 13281191951274694749, v, C)
-		_, u, v = madd1s(x[0], x[3], t, u)
-		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
-	}
-	{
-		// round 1
-		m := p[0] * 14042775128853446655
-		C := madd0(m, 4891460686036598785, p[0])
-		u, v = madd1(x[1], x[1], p[1])
-		C, p[0] = madd2(m, 2896914383306846353, v, C)
-		var t uint64
-		t, u, v = madd2sb(x[1], x[2], p[2], u)
-		C, p[1] = madd2(m, 13281191951274694749, v, C)
-		_, u, v = madd2s(x[1], x[3], p[3], t, u)
-		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
-	}
-	{
-		// round 2
-		m := p[0] * 14042775128853446655
-		C := madd0(m, 4891460686036598785, p[0])
-		C, p[0] = madd2(m, 2896914383306846353, p[1], C)
-		u, v = madd1(x[2], x[2], p[2])
-		C, p[1] = madd2(m, 13281191951274694749, v, C)
-		_, u, v = madd2sb(x[2], x[3], p[3], u)
-		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
-	}
-	{
-		// round 3
-		m := p[0] * 14042775128853446655
-		C := madd0(m, 4891460686036598785, p[0])
-		C, z[0] = madd2(m, 2896914383306846353, p[1], C)
-		C, z[1] = madd2(m, 13281191951274694749, p[2], C)
-		u, v = madd1(x[3], x[3], p[3])
-		z[3], z[2] = madd3(m, 3486998266802970665, v, C, u)
-	}
-
-	// if z > q --> z -= q
-	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
-		var b uint64
-		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
-		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
-		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
-		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
-	}
-	return z
 }

ff/element_mul.go

@@ -0,0 +1,170 @@
+// +build !amd64
+
+// Copyright 2020 ConsenSys AG
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by goff (v0.2.0) DO NOT EDIT
+
+// Package ff contains field arithmetic operations
+package ff
+
+// /!\ WARNING /!\
+// this code has not been audited and is provided as-is. In particular,
+// there is no security guarantees such as constant time implementation
+// or side-channel attack resistance
+// /!\ WARNING /!\
+
+import "math/bits"
+
+// Mul z = x * y mod q
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) Mul(x, y *Element) *Element {
+
+	var t [4]uint64
+	var c [3]uint64
+	{
+		// round 0
+		v := x[0]
+		c[1], c[0] = bits.Mul64(v, y[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd1(v, y[1], c[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd1(v, y[2], c[1])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd1(v, y[3], c[1])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 1
+		v := x[1]
+		c[1], c[0] = madd1(v, y[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, y[1], c[1], t[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, y[2], c[1], t[2])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, y[3], c[1], t[3])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 2
+		v := x[2]
+		c[1], c[0] = madd1(v, y[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, y[1], c[1], t[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, y[2], c[1], t[2])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, y[3], c[1], t[3])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 3
+		v := x[3]
+		c[1], c[0] = madd1(v, y[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, y[1], c[1], t[1])
+		c[2], z[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, y[2], c[1], t[2])
+		c[2], z[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, y[3], c[1], t[3])
+		z[3], z[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
+		var b uint64
+		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
+		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
+		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
+		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
+	}
+	return z
+}
+
+// MulAssign z = z * x mod q
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) MulAssign(x *Element) *Element {
+
+	var t [4]uint64
+	var c [3]uint64
+	{
+		// round 0
+		v := z[0]
+		c[1], c[0] = bits.Mul64(v, x[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd1(v, x[1], c[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd1(v, x[2], c[1])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd1(v, x[3], c[1])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 1
+		v := z[1]
+		c[1], c[0] = madd1(v, x[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, x[1], c[1], t[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, x[2], c[1], t[2])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, x[3], c[1], t[3])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 2
+		v := z[2]
+		c[1], c[0] = madd1(v, x[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, x[1], c[1], t[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, x[2], c[1], t[2])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, x[3], c[1], t[3])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 3
+		v := z[3]
+		c[1], c[0] = madd1(v, x[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, x[1], c[1], t[1])
+		c[2], z[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, x[2], c[1], t[2])
+		c[2], z[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, x[3], c[1], t[3])
+		z[3], z[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
+		var b uint64
+		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
+		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
+		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
+		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
+	}
+	return z
+}

ff/element_mul_amd64.go

@@ -0,0 +1,39 @@
+// Copyright 2020 ConsenSys AG
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by goff (v0.2.0) DO NOT EDIT
+
+// Package ff contains field arithmetic operations
+package ff
+
+// MulAssignElement z = z * x mod q (constant time)
+// calling this instead of z.MulAssign(x) is prefered for performance critical path
+//go:noescape
+func MulAssignElement(res, y *Element)
+
+// Mul z = x * y mod q (constant time)
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) Mul(x, y *Element) *Element {
+	res := *x
+	MulAssignElement(&res, y)
+	z.Set(&res)
+	return z
+}
+
+// MulAssign z = z * x mod q (constant time)
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) MulAssign(x *Element) *Element {
+	MulAssignElement(z, x)
+	return z
+}

ff/element_mul_amd64.s

@@ -0,0 +1,695 @@
+// Code generated by goff (v0.2.0) DO NOT EDIT
+
+#include "textflag.h"
+
+// func MulAssignElement(res,y *Element)
+// montgomery multiplication of res by y
+// stores the result in res
+TEXT ·MulAssignElement(SB), NOSPLIT, $0-16
+    
+    // dereference our parameters
+    MOVQ res+0(FP), DI
+    MOVQ y+8(FP), R8
+    
+    // check if we support adx and mulx
+    CMPB ·supportAdx(SB), $1
+    JNE no_adx
+    
+    // the algorithm is described here
+    // https://hackmd.io/@zkteam/modular_multiplication
+    // however, to benefit from the ADCX and ADOX carry chains
+    // we split the inner loops in 2:
+    // for i=0 to N-1
+    // 		for j=0 to N-1
+    // 		    (A,t[j])  := t[j] + a[j]*b[i] + A
+    // 		m := t[0]*q'[0] mod W
+    // 		C,_ := t[0] + m*q[0]
+    // 		for j=1 to N-1
+    // 		    (C,t[j-1]) := t[j] + m*q[j] + C
+    // 		t[N-1] = C + A
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 0
+    
+    // clear up the carry flags
+    XORQ R9 , R9
+    
+    // R12 = y[0]
+    MOVQ 0(R8), R12
+    
+    // for j=0 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    
+    // DX = res[0]
+    MOVQ 0(DI), DX
+    MULXQ R12, CX ,  R9
+    
+    // DX = res[1]
+    MOVQ 8(DI), DX
+    MOVQ R9, BX
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BX
+    
+    // DX = res[2]
+    MOVQ 16(DI), DX
+    MOVQ R9, BP
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BP
+    
+    // DX = res[3]
+    MOVQ 24(DI), DX
+    MOVQ R9, SI
+    MULXQ R12, AX,  R9
+    ADOXQ AX, SI
+    
+    // add the last carries to R9
+    MOVQ $0, DX
+    ADCXQ DX, R9
+    ADOXQ DX, R9
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, DX
+    MULXQ CX,R11, DX
+    
+    // clear the carry flags
+    XORQ DX, DX
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, DX
+    MULXQ R11, AX, R10
+    ADCXQ CX ,AX
+    
+    // for j=1 to N-1
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    
+    MOVQ $0x2833e84879b97091, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BX, R10
+    ADOXQ AX, R10
+    MOVQ R10, CX
+    MOVQ DX, R10
+    
+    MOVQ $0xb85045b68181585d, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BP, R10
+    ADOXQ AX, R10
+    MOVQ R10, BX
+    MOVQ DX, R10
+    
+    MOVQ $0x30644e72e131a029, DX
+    MULXQ R11, AX, DX
+    ADCXQ  SI, R10
+    ADOXQ AX, R10
+    MOVQ R10, BP
+    MOVQ $0, AX
+    ADCXQ AX, DX
+    ADOXQ DX, R9
+    MOVQ R9, SI
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 1
+    
+    // clear up the carry flags
+    XORQ R9 , R9
+    
+    // R12 = y[1]
+    MOVQ 8(R8), R12
+    
+    // for j=0 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    
+    // DX = res[0]
+    MOVQ 0(DI), DX
+    MULXQ R12, AX,  R9
+    ADOXQ AX, CX
+    
+    // DX = res[1]
+    MOVQ 8(DI), DX
+    ADCXQ R9, BX
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BX
+    
+    // DX = res[2]
+    MOVQ 16(DI), DX
+    ADCXQ R9, BP
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BP
+    
+    // DX = res[3]
+    MOVQ 24(DI), DX
+    ADCXQ R9, SI
+    MULXQ R12, AX,  R9
+    ADOXQ AX, SI
+    
+    // add the last carries to R9
+    MOVQ $0, DX
+    ADCXQ DX, R9
+    ADOXQ DX, R9
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, DX
+    MULXQ CX,R11, DX
+    
+    // clear the carry flags
+    XORQ DX, DX
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, DX
+    MULXQ R11, AX, R10
+    ADCXQ CX ,AX
+    
+    // for j=1 to N-1
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    
+    MOVQ $0x2833e84879b97091, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BX, R10
+    ADOXQ AX, R10
+    MOVQ R10, CX
+    MOVQ DX, R10
+    
+    MOVQ $0xb85045b68181585d, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BP, R10
+    ADOXQ AX, R10
+    MOVQ R10, BX
+    MOVQ DX, R10
+    
+    MOVQ $0x30644e72e131a029, DX
+    MULXQ R11, AX, DX
+    ADCXQ  SI, R10
+    ADOXQ AX, R10
+    MOVQ R10, BP
+    MOVQ $0, AX
+    ADCXQ AX, DX
+    ADOXQ DX, R9
+    MOVQ R9, SI
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 2
+    
+    // clear up the carry flags
+    XORQ R9 , R9
+    
+    // R12 = y[2]
+    MOVQ 16(R8), R12
+    
+    // for j=0 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    
+    // DX = res[0]
+    MOVQ 0(DI), DX
+    MULXQ R12, AX,  R9
+    ADOXQ AX, CX
+    
+    // DX = res[1]
+    MOVQ 8(DI), DX
+    ADCXQ R9, BX
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BX
+    
+    // DX = res[2]
+    MOVQ 16(DI), DX
+    ADCXQ R9, BP
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BP
+    
+    // DX = res[3]
+    MOVQ 24(DI), DX
+    ADCXQ R9, SI
+    MULXQ R12, AX,  R9
+    ADOXQ AX, SI
+    
+    // add the last carries to R9
+    MOVQ $0, DX
+    ADCXQ DX, R9
+    ADOXQ DX, R9
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, DX
+    MULXQ CX,R11, DX
+    
+    // clear the carry flags
+    XORQ DX, DX
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, DX
+    MULXQ R11, AX, R10
+    ADCXQ CX ,AX
+    
+    // for j=1 to N-1
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    
+    MOVQ $0x2833e84879b97091, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BX, R10
+    ADOXQ AX, R10
+    MOVQ R10, CX
+    MOVQ DX, R10
+    
+    MOVQ $0xb85045b68181585d, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BP, R10
+    ADOXQ AX, R10
+    MOVQ R10, BX
+    MOVQ DX, R10
+    
+    MOVQ $0x30644e72e131a029, DX
+    MULXQ R11, AX, DX
+    ADCXQ  SI, R10
+    ADOXQ AX, R10
+    MOVQ R10, BP
+    MOVQ $0, AX
+    ADCXQ AX, DX
+    ADOXQ DX, R9
+    MOVQ R9, SI
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 3
+    
+    // clear up the carry flags
+    XORQ R9 , R9
+    
+    // R12 = y[3]
+    MOVQ 24(R8), R12
+    
+    // for j=0 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    
+    // DX = res[0]
+    MOVQ 0(DI), DX
+    MULXQ R12, AX,  R9
+    ADOXQ AX, CX
+    
+    // DX = res[1]
+    MOVQ 8(DI), DX
+    ADCXQ R9, BX
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BX
+    
+    // DX = res[2]
+    MOVQ 16(DI), DX
+    ADCXQ R9, BP
+    MULXQ R12, AX,  R9
+    ADOXQ AX, BP
+    
+    // DX = res[3]
+    MOVQ 24(DI), DX
+    ADCXQ R9, SI
+    MULXQ R12, AX,  R9
+    ADOXQ AX, SI
+    
+    // add the last carries to R9
+    MOVQ $0, DX
+    ADCXQ DX, R9
+    ADOXQ DX, R9
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, DX
+    MULXQ CX,R11, DX
+    
+    // clear the carry flags
+    XORQ DX, DX
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, DX
+    MULXQ R11, AX, R10
+    ADCXQ CX ,AX
+    
+    // for j=1 to N-1
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    
+    MOVQ $0x2833e84879b97091, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BX, R10
+    ADOXQ AX, R10
+    MOVQ R10, CX
+    MOVQ DX, R10
+    
+    MOVQ $0xb85045b68181585d, DX
+    MULXQ R11, AX, DX
+    ADCXQ  BP, R10
+    ADOXQ AX, R10
+    MOVQ R10, BX
+    MOVQ DX, R10
+    
+    MOVQ $0x30644e72e131a029, DX
+    MULXQ R11, AX, DX
+    ADCXQ  SI, R10
+    ADOXQ AX, R10
+    MOVQ R10, BP
+    MOVQ $0, AX
+    ADCXQ AX, DX
+    ADOXQ DX, R9
+    MOVQ R9, SI
+    
+    reduce:
+    // reduce, constant time version
+    // first we copy registers storing t in a separate set of registers
+    // as SUBQ modifies the 2nd operand
+    MOVQ CX, DX
+    MOVQ BX, R8
+    MOVQ BP, R9
+    MOVQ SI, R10
+    MOVQ $0x43e1f593f0000001, R11
+    SUBQ  R11, DX
+    MOVQ $0x2833e84879b97091, R11
+    SBBQ  R11, R8
+    MOVQ $0xb85045b68181585d, R11
+    SBBQ  R11, R9
+    MOVQ $0x30644e72e131a029, R11
+    SBBQ  R11, R10
+    JCS t_is_smaller // no borrow, we return t
+    
+    // borrow is set, we return u
+    MOVQ DX, (DI)
+    MOVQ R8, 8(DI)
+    MOVQ R9, 16(DI)
+    MOVQ R10, 24(DI)
+    RET
+    t_is_smaller:
+    MOVQ CX, 0(DI)
+    MOVQ BX, 8(DI)
+    MOVQ BP, 16(DI)
+    MOVQ SI, 24(DI)
+    RET
+    
+    no_adx:
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 0
+    
+    // (A,t[0]) := t[0] + x[0]*y[0]
+    MOVQ (DI), AX // x[0]
+    MOVQ 0(R8), R12
+    MULQ R12 // x[0] * y[0]
+    MOVQ DX, R9
+    MOVQ AX, CX
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, R11
+    IMULQ CX , R11
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, AX
+    MULQ R11
+    ADDQ CX ,AX
+    ADCQ $0, DX
+    MOVQ  DX, R10
+    
+    // for j=1 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    MOVQ 8(DI), AX
+    MULQ R12 // x[1] * y[0]
+    MOVQ R9, BX
+    ADDQ AX, BX
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x2833e84879b97091, AX
+    MULQ R11
+    ADDQ  BX, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, CX
+    MOVQ DX, R10
+    MOVQ 16(DI), AX
+    MULQ R12 // x[2] * y[0]
+    MOVQ R9, BP
+    ADDQ AX, BP
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0xb85045b68181585d, AX
+    MULQ R11
+    ADDQ  BP, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BX
+    MOVQ DX, R10
+    MOVQ 24(DI), AX
+    MULQ R12 // x[3] * y[0]
+    MOVQ R9, SI
+    ADDQ AX, SI
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x30644e72e131a029, AX
+    MULQ R11
+    ADDQ  SI, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BP
+    MOVQ DX, R10
+    
+    ADDQ R10, R9
+    MOVQ R9, SI
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 1
+    
+    // (A,t[0]) := t[0] + x[0]*y[1]
+    MOVQ (DI), AX // x[0]
+    MOVQ 8(R8), R12
+    MULQ R12 // x[0] * y[1]
+    ADDQ AX, CX
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, R11
+    IMULQ CX , R11
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, AX
+    MULQ R11
+    ADDQ CX ,AX
+    ADCQ $0, DX
+    MOVQ  DX, R10
+    
+    // for j=1 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    MOVQ 8(DI), AX
+    MULQ R12 // x[1] * y[1]
+    ADDQ R9, BX
+    ADCQ $0, DX
+    ADDQ AX, BX
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x2833e84879b97091, AX
+    MULQ R11
+    ADDQ  BX, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, CX
+    MOVQ DX, R10
+    MOVQ 16(DI), AX
+    MULQ R12 // x[2] * y[1]
+    ADDQ R9, BP
+    ADCQ $0, DX
+    ADDQ AX, BP
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0xb85045b68181585d, AX
+    MULQ R11
+    ADDQ  BP, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BX
+    MOVQ DX, R10
+    MOVQ 24(DI), AX
+    MULQ R12 // x[3] * y[1]
+    ADDQ R9, SI
+    ADCQ $0, DX
+    ADDQ AX, SI
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x30644e72e131a029, AX
+    MULQ R11
+    ADDQ  SI, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BP
+    MOVQ DX, R10
+    
+    ADDQ R10, R9
+    MOVQ R9, SI
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 2
+    
+    // (A,t[0]) := t[0] + x[0]*y[2]
+    MOVQ (DI), AX // x[0]
+    MOVQ 16(R8), R12
+    MULQ R12 // x[0] * y[2]
+    ADDQ AX, CX
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, R11
+    IMULQ CX , R11
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, AX
+    MULQ R11
+    ADDQ CX ,AX
+    ADCQ $0, DX
+    MOVQ  DX, R10
+    
+    // for j=1 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    MOVQ 8(DI), AX
+    MULQ R12 // x[1] * y[2]
+    ADDQ R9, BX
+    ADCQ $0, DX
+    ADDQ AX, BX
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x2833e84879b97091, AX
+    MULQ R11
+    ADDQ  BX, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, CX
+    MOVQ DX, R10
+    MOVQ 16(DI), AX
+    MULQ R12 // x[2] * y[2]
+    ADDQ R9, BP
+    ADCQ $0, DX
+    ADDQ AX, BP
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0xb85045b68181585d, AX
+    MULQ R11
+    ADDQ  BP, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BX
+    MOVQ DX, R10
+    MOVQ 24(DI), AX
+    MULQ R12 // x[3] * y[2]
+    ADDQ R9, SI
+    ADCQ $0, DX
+    ADDQ AX, SI
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x30644e72e131a029, AX
+    MULQ R11
+    ADDQ  SI, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BP
+    MOVQ DX, R10
+    
+    ADDQ R10, R9
+    MOVQ R9, SI
+    
+    // ---------------------------------------------------------------------------------------------
+    // outter loop 3
+    
+    // (A,t[0]) := t[0] + x[0]*y[3]
+    MOVQ (DI), AX // x[0]
+    MOVQ 24(R8), R12
+    MULQ R12 // x[0] * y[3]
+    ADDQ AX, CX
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    // m := t[0]*q'[0] mod W
+    MOVQ $0xc2e1f593efffffff, R11
+    IMULQ CX , R11
+    
+    // C,_ := t[0] + m*q[0]
+    MOVQ $0x43e1f593f0000001, AX
+    MULQ R11
+    ADDQ CX ,AX
+    ADCQ $0, DX
+    MOVQ  DX, R10
+    
+    // for j=1 to N-1
+    //    (A,t[j])  := t[j] + x[j]*y[i] + A
+    //    (C,t[j-1]) := t[j] + m*q[j] + C
+    MOVQ 8(DI), AX
+    MULQ R12 // x[1] * y[3]
+    ADDQ R9, BX
+    ADCQ $0, DX
+    ADDQ AX, BX
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x2833e84879b97091, AX
+    MULQ R11
+    ADDQ  BX, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, CX
+    MOVQ DX, R10
+    MOVQ 16(DI), AX
+    MULQ R12 // x[2] * y[3]
+    ADDQ R9, BP
+    ADCQ $0, DX
+    ADDQ AX, BP
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0xb85045b68181585d, AX
+    MULQ R11
+    ADDQ  BP, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BX
+    MOVQ DX, R10
+    MOVQ 24(DI), AX
+    MULQ R12 // x[3] * y[3]
+    ADDQ R9, SI
+    ADCQ $0, DX
+    ADDQ AX, SI
+    ADCQ $0, DX
+    MOVQ DX, R9
+    
+    MOVQ $0x30644e72e131a029, AX
+    MULQ R11
+    ADDQ  SI, R10
+    ADCQ $0, DX
+    ADDQ AX, R10
+    ADCQ $0, DX
+    
+    MOVQ R10, BP
+    MOVQ DX, R10
+    
+    ADDQ R10, R9
+    MOVQ R9, SI
+    
+    JMP reduce

ff/element_square.go

@@ -0,0 +1,93 @@
+// +build !amd64
+
+// Copyright 2020 ConsenSys AG
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by goff (v0.2.0) DO NOT EDIT
+
+// Package ff contains field arithmetic operations
+package ff
+
+// /!\ WARNING /!\
+// this code has not been audited and is provided as-is. In particular,
+// there is no security guarantees such as constant time implementation
+// or side-channel attack resistance
+// /!\ WARNING /!\
+
+import "math/bits"
+
+// Square z = x * x mod q
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) Square(x *Element) *Element {
+
+	var p [4]uint64
+
+	var u, v uint64
+	{
+		// round 0
+		u, p[0] = bits.Mul64(x[0], x[0])
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		var t uint64
+		t, u, v = madd1sb(x[0], x[1], u)
+		C, p[0] = madd2(m, 2896914383306846353, v, C)
+		t, u, v = madd1s(x[0], x[2], t, u)
+		C, p[1] = madd2(m, 13281191951274694749, v, C)
+		_, u, v = madd1s(x[0], x[3], t, u)
+		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+	{
+		// round 1
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		u, v = madd1(x[1], x[1], p[1])
+		C, p[0] = madd2(m, 2896914383306846353, v, C)
+		var t uint64
+		t, u, v = madd2sb(x[1], x[2], p[2], u)
+		C, p[1] = madd2(m, 13281191951274694749, v, C)
+		_, u, v = madd2s(x[1], x[3], p[3], t, u)
+		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+	{
+		// round 2
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		C, p[0] = madd2(m, 2896914383306846353, p[1], C)
+		u, v = madd1(x[2], x[2], p[2])
+		C, p[1] = madd2(m, 13281191951274694749, v, C)
+		_, u, v = madd2sb(x[2], x[3], p[3], u)
+		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+	{
+		// round 3
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		C, z[0] = madd2(m, 2896914383306846353, p[1], C)
+		C, z[1] = madd2(m, 13281191951274694749, p[2], C)
+		u, v = madd1(x[3], x[3], p[3])
+		z[3], z[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
+		var b uint64
+		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
+		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
+		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
+		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
+	}
+	return z
+
+}

ff/element_square_amd64.go

@@ -0,0 +1,34 @@
+// Copyright 2020 ConsenSys AG
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by goff (v0.2.0) DO NOT EDIT
+
+// Package ff contains field arithmetic operations
+package ff
+
+// SquareElement z = x * x mod q
+// calling this instead of z.Square(x) is prefered for performance critical path
+// go - noescape
+// func SquareElement(res,x *Element)
+
+// Square z = x * x mod q
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) Square(x *Element) *Element {
+	if z != x {
+		z.Set(x)
+	}
+	MulAssignElement(z, x)
+	// SquareElement(z, x)
+	return z
+}

ff/element_test.go

@@ -1,9 +1,26 @@
-// Code generated by goff DO NOT EDIT
+// Copyright 2020 ConsenSys AG
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// Code generated by goff (v0.2.0) DO NOT EDIT
+
+// Package ff contains field arithmetic operations
 package ff
 
 import (
 	"crypto/rand"
 	"math/big"
+	"math/bits"
 	mrand "math/rand"
 	"testing"
 )
@@ -21,7 +38,14 @@ func TestELEMENTCorrectnessAgainstBigInt(t *testing.T) {
 
 	modulusMinusOne.Sub(modulus, &one)
 
-	for i := 0; i < 1000; i++ {
+	var n int
+	if testing.Short() {
+		n = 10
+	} else {
+		n = 500
+	}
+
+	for i := 0; i < n; i++ {
 
 		// sample 2 random big int
 		b1, _ := rand.Int(rand.Reader, modulus)
@@ -57,7 +81,7 @@ func TestELEMENTCorrectnessAgainstBigInt(t *testing.T) {
 
 		rbExp := new(big.Int).SetUint64(rExp)
 
-		var bMul, bAdd, bSub, bDiv, bNeg, bLsh, bInv, bExp, bSquare big.Int
+		var bMul, bAdd, bSub, bDiv, bNeg, bLsh, bInv, bExp, bExp2, bSquare big.Int
 
 		// e1 = mont(b1), e2 = mont(b2)
 		var e1, e2, eMul, eAdd, eSub, eDiv, eNeg, eLsh, eInv, eExp, eSquare, eMulAssign, eSubAssign, eAddAssign Element
@@ -106,12 +130,40 @@ func TestELEMENTCorrectnessAgainstBigInt(t *testing.T) {
 		cmpEandB(&eNeg, &bNeg, "Neg")
 		cmpEandB(&eInv, &bInv, "Inv")
 		cmpEandB(&eExp, &bExp, "Exp")
+
 		cmpEandB(&eLsh, &bLsh, "Lsh")
+
+		// legendre symbol
+		if e1.Legendre() != big.Jacobi(b1, modulus) {
+			t.Fatal("legendre symbol computation failed")
+		}
+		if e2.Legendre() != big.Jacobi(b2, modulus) {
+			t.Fatal("legendre symbol computation failed")
+		}
+
+		// these are slow, killing circle ci
+		if n <= 5 {
+			// sqrt
+			var eSqrt, eExp2 Element
+			var bSqrt big.Int
+			bSqrt.ModSqrt(b1, modulus)
+			eSqrt.Sqrt(&e1)
+			cmpEandB(&eSqrt, &bSqrt, "Sqrt")
+
+			bits := b2.Bits()
+			exponent := make([]uint64, len(bits))
+			for k := 0; k < len(bits); k++ {
+				exponent[k] = uint64(bits[k])
+			}
+			eExp2.Exp(e1, exponent...)
+			bExp2.Exp(b1, b2, modulus)
+			cmpEandB(&eExp2, &bExp2, "Exp multi words")
+		}
 	}
 }
 
 func TestELEMENTIsRandom(t *testing.T) {
-	for i := 0; i < 1000; i++ {
+	for i := 0; i < 50; i++ {
 		var x, y Element
 		x.SetRandom()
 		y.SetRandom()
@@ -125,7 +177,6 @@ func TestELEMENTIsRandom(t *testing.T) {
 // benchmarks
 // most benchmarks are rudimentary and should sample a large number of random inputs
 // or be run multiple times to ensure it didn't measure the fastest path of the function
-// TODO: clean up and push benchmarking branch
 
 var benchResElement Element
 
@@ -219,6 +270,15 @@ func BenchmarkSquareELEMENT(b *testing.B) {
 	}
 }
 
+func BenchmarkSqrtELEMENT(b *testing.B) {
+	var a Element
+	a.SetRandom()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		benchResElement.Sqrt(&a)
+	}
+}
+
 func BenchmarkMulAssignELEMENT(b *testing.B) {
 	x := Element{
 		1997599621687373223,
@@ -232,3 +292,183 @@ func BenchmarkMulAssignELEMENT(b *testing.B) {
 		benchResElement.MulAssign(&x)
 	}
 }
+
+func BenchmarkMulAssignASMELEMENT(b *testing.B) {
+	x := Element{
+		1997599621687373223,
+		6052339484930628067,
+		10108755138030829701,
+		150537098327114917,
+	}
+	benchResElement.SetOne()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		MulAssignElement(&benchResElement, &x)
+	}
+}
+
+func TestELEMENTAsm(t *testing.T) {
+	// ensure ASM implementations matches the ones using math/bits
+	modulus, _ := new(big.Int).SetString("21888242871839275222246405745257275088548364400416034343698204186575808495617", 10)
+	for i := 0; i < 500; i++ {
+		// sample 2 random big int
+		b1, _ := rand.Int(rand.Reader, modulus)
+		b2, _ := rand.Int(rand.Reader, modulus)
+
+		// e1 = mont(b1), e2 = mont(b2)
+		var e1, e2, eTestMul, eMulAssign, eSquare, eTestSquare Element
+		e1.SetBigInt(b1)
+		e2.SetBigInt(b2)
+
+		eTestMul = e1
+		eTestMul.testMulAssign(&e2)
+		eMulAssign = e1
+		eMulAssign.MulAssign(&e2)
+
+		if !eTestMul.Equal(&eMulAssign) {
+			t.Fatal("inconsisntencies between MulAssign and testMulAssign --> check if MulAssign is calling ASM implementaiton on amd64")
+		}
+
+		// square
+		eSquare.Square(&e1)
+		eTestSquare.testSquare(&e1)
+
+		if !eTestSquare.Equal(&eSquare) {
+			t.Fatal("inconsisntencies between Square and testSquare --> check if Square is calling ASM implementaiton on amd64")
+		}
+	}
+}
+
+// this is here for consistency purposes, to ensure MulAssign on AMD64 using asm implementation gives consistent results
+func (z *Element) testMulAssign(x *Element) *Element {
+
+	var t [4]uint64
+	var c [3]uint64
+	{
+		// round 0
+		v := z[0]
+		c[1], c[0] = bits.Mul64(v, x[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd1(v, x[1], c[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd1(v, x[2], c[1])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd1(v, x[3], c[1])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 1
+		v := z[1]
+		c[1], c[0] = madd1(v, x[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, x[1], c[1], t[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, x[2], c[1], t[2])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, x[3], c[1], t[3])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 2
+		v := z[2]
+		c[1], c[0] = madd1(v, x[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, x[1], c[1], t[1])
+		c[2], t[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, x[2], c[1], t[2])
+		c[2], t[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, x[3], c[1], t[3])
+		t[3], t[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+	{
+		// round 3
+		v := z[3]
+		c[1], c[0] = madd1(v, x[0], t[0])
+		m := c[0] * 14042775128853446655
+		c[2] = madd0(m, 4891460686036598785, c[0])
+		c[1], c[0] = madd2(v, x[1], c[1], t[1])
+		c[2], z[0] = madd2(m, 2896914383306846353, c[2], c[0])
+		c[1], c[0] = madd2(v, x[2], c[1], t[2])
+		c[2], z[1] = madd2(m, 13281191951274694749, c[2], c[0])
+		c[1], c[0] = madd2(v, x[3], c[1], t[3])
+		z[3], z[2] = madd3(m, 3486998266802970665, c[0], c[2], c[1])
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
+		var b uint64
+		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
+		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
+		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
+		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
+	}
+	return z
+}
+
+// this is here for consistency purposes, to ensure Square on AMD64 using asm implementation gives consistent results
+func (z *Element) testSquare(x *Element) *Element {
+
+	var p [4]uint64
+
+	var u, v uint64
+	{
+		// round 0
+		u, p[0] = bits.Mul64(x[0], x[0])
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		var t uint64
+		t, u, v = madd1sb(x[0], x[1], u)
+		C, p[0] = madd2(m, 2896914383306846353, v, C)
+		t, u, v = madd1s(x[0], x[2], t, u)
+		C, p[1] = madd2(m, 13281191951274694749, v, C)
+		_, u, v = madd1s(x[0], x[3], t, u)
+		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+	{
+		// round 1
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		u, v = madd1(x[1], x[1], p[1])
+		C, p[0] = madd2(m, 2896914383306846353, v, C)
+		var t uint64
+		t, u, v = madd2sb(x[1], x[2], p[2], u)
+		C, p[1] = madd2(m, 13281191951274694749, v, C)
+		_, u, v = madd2s(x[1], x[3], p[3], t, u)
+		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+	{
+		// round 2
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		C, p[0] = madd2(m, 2896914383306846353, p[1], C)
+		u, v = madd1(x[2], x[2], p[2])
+		C, p[1] = madd2(m, 13281191951274694749, v, C)
+		_, u, v = madd2sb(x[2], x[3], p[3], u)
+		p[3], p[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+	{
+		// round 3
+		m := p[0] * 14042775128853446655
+		C := madd0(m, 4891460686036598785, p[0])
+		C, z[0] = madd2(m, 2896914383306846353, p[1], C)
+		C, z[1] = madd2(m, 13281191951274694749, p[2], C)
+		u, v = madd1(x[3], x[3], p[3])
+		z[3], z[2] = madd3(m, 3486998266802970665, v, C, u)
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[3] < 3486998266802970665 || (z[3] == 3486998266802970665 && (z[2] < 13281191951274694749 || (z[2] == 13281191951274694749 && (z[1] < 2896914383306846353 || (z[1] == 2896914383306846353 && (z[0] < 4891460686036598785))))))) {
+		var b uint64
+		z[0], b = bits.Sub64(z[0], 4891460686036598785, 0)
+		z[1], b = bits.Sub64(z[1], 2896914383306846353, b)
+		z[2], b = bits.Sub64(z[2], 13281191951274694749, b)
+		z[3], _ = bits.Sub64(z[3], 3486998266802970665, b)
+	}
+	return z
+
+}

field/field.go

@@ -1,152 +0,0 @@
-// code originally taken from https://github.com/arnaucube/go-snark (https://github.com/arnaucube/go-snark/blob/master/fields/fq.go), pasted here to ensure compatibility among future changes
-
-package field
-
-import (
-	"bytes"
-	"crypto/rand"
-	"math/big"
-)
-
-// Fq is the Z field over modulus Q
-type Fq struct {
-	Q *big.Int // Q
-}
-
-// NewFq generates a new Fq
-func NewFq(q *big.Int) Fq {
-	return Fq{
-		q,
-	}
-}
-
-// Zero returns a Zero value on the Fq
-func (fq Fq) Zero() *big.Int {
-	return big.NewInt(int64(0))
-}
-
-// One returns a One value on the Fq
-func (fq Fq) One() *big.Int {
-	return big.NewInt(int64(1))
-}
-
-// Add performs an addition on the Fq
-func (fq Fq) Add(a, b *big.Int) *big.Int {
-	r := new(big.Int).Add(a, b)
-	return new(big.Int).Mod(r, fq.Q)
-}
-
-// Double performs a doubling on the Fq
-func (fq Fq) Double(a *big.Int) *big.Int {
-	r := new(big.Int).Add(a, a)
-	return new(big.Int).Mod(r, fq.Q)
-}
-
-// Sub performs a subtraction on the Fq
-func (fq Fq) Sub(a, b *big.Int) *big.Int {
-	r := new(big.Int).Sub(a, b)
-	return new(big.Int).Mod(r, fq.Q)
-}
-
-// Neg performs a negation on the Fq
-func (fq Fq) Neg(a *big.Int) *big.Int {
-	m := new(big.Int).Neg(a)
-	return new(big.Int).Mod(m, fq.Q)
-}
-
-// Mul performs a multiplication on the Fq
-func (fq Fq) Mul(a, b *big.Int) *big.Int {
-	m := new(big.Int).Mul(a, b)
-	return new(big.Int).Mod(m, fq.Q)
-}
-
-func (fq Fq) MulScalar(base, e *big.Int) *big.Int {
-	return fq.Mul(base, e)
-}
-
-// Inverse returns the inverse on the Fq
-func (fq Fq) Inverse(a *big.Int) *big.Int {
-	return new(big.Int).ModInverse(a, fq.Q)
-}
-
-// Div performs the division over the finite field
-func (fq Fq) Div(a, b *big.Int) *big.Int {
-	d := fq.Mul(a, fq.Inverse(b))
-	return new(big.Int).Mod(d, fq.Q)
-}
-
-// Square performs a square operation on the Fq
-func (fq Fq) Square(a *big.Int) *big.Int {
-	m := new(big.Int).Mul(a, a)
-	return new(big.Int).Mod(m, fq.Q)
-}
-
-// Exp performs the exponential over Fq
-func (fq Fq) Exp(base *big.Int, e *big.Int) *big.Int {
-	res := fq.One()
-	rem := fq.Copy(e)
-	exp := base
-
-	for !bytes.Equal(rem.Bytes(), big.NewInt(int64(0)).Bytes()) {
-		if BigIsOdd(rem) {
-			res = fq.Mul(res, exp)
-		}
-		exp = fq.Square(exp)
-		rem = new(big.Int).Rsh(rem, 1)
-	}
-	return res
-}
-
-func (fq Fq) Rand() (*big.Int, error) {
-
-	maxbits := fq.Q.BitLen()
-	b := make([]byte, (maxbits/8)-1)
-	_, err := rand.Read(b)
-	if err != nil {
-		return nil, err
-	}
-	r := new(big.Int).SetBytes(b)
-	rq := new(big.Int).Mod(r, fq.Q)
-
-	// r over q, nil
-	return rq, nil
-}
-
-func (fq Fq) IsZero(a *big.Int) bool {
-	return bytes.Equal(a.Bytes(), fq.Zero().Bytes())
-}
-
-func (fq Fq) Copy(a *big.Int) *big.Int {
-	return new(big.Int).SetBytes(a.Bytes())
-}
-
-func (fq Fq) Affine(a *big.Int) *big.Int {
-	nq := fq.Neg(fq.Q)
-
-	aux := a
-	if aux.Cmp(big.NewInt(int64(0))) == -1 { // negative value
-		if aux.Cmp(nq) != 1 { // aux less or equal nq
-			aux = new(big.Int).Mod(aux, fq.Q)
-		}
-		if aux.Cmp(big.NewInt(int64(0))) == -1 { // negative value
-			aux = new(big.Int).Add(aux, fq.Q)
-		}
-	} else {
-		if aux.Cmp(fq.Q) != -1 { // aux greater or equal nq
-			aux = new(big.Int).Mod(aux, fq.Q)
-		}
-	}
-	return aux
-}
-
-func (fq Fq) Equal(a, b *big.Int) bool {
-	aAff := fq.Affine(a)
-	bAff := fq.Affine(b)
-	return bytes.Equal(aAff.Bytes(), bAff.Bytes())
-}
-
-func BigIsOdd(n *big.Int) bool {
-	one := big.NewInt(int64(1))
-	and := new(big.Int).And(n, one)
-	return bytes.Equal(and.Bytes(), big.NewInt(int64(1)).Bytes())
-}

field/field_test.go

@@ -1,39 +0,0 @@
-// code originally taken from https://github.com/arnaucube/go-snark (https://github.com/arnaucube/go-snark/blob/master/fields/fq.go), pasted here to ensure compatibility among future changes
-
-package field
-
-import (
-	"math/big"
-	"testing"
-
-	"github.com/stretchr/testify/assert"
-)
-
-func iToBig(a int) *big.Int {
-	return big.NewInt(int64(a))
-}
-
-func TestFq1(t *testing.T) {
-	fq1 := NewFq(iToBig(7))
-
-	res := fq1.Add(iToBig(4), iToBig(4))
-	assert.Equal(t, iToBig(1), fq1.Affine(res))
-
-	res = fq1.Double(iToBig(5))
-	assert.Equal(t, iToBig(3), fq1.Affine(res))
-
-	res = fq1.Sub(iToBig(5), iToBig(7))
-	assert.Equal(t, iToBig(5), fq1.Affine(res))
-
-	res = fq1.Neg(iToBig(5))
-	assert.Equal(t, iToBig(2), fq1.Affine(res))
-
-	res = fq1.Mul(iToBig(5), iToBig(11))
-	assert.Equal(t, iToBig(6), fq1.Affine(res))
-
-	res = fq1.Inverse(iToBig(4))
-	assert.Equal(t, iToBig(2), res)
-
-	res = fq1.Square(iToBig(5))
-	assert.Equal(t, iToBig(4), res)
-}

go.mod

@@ -7,4 +7,5 @@ require (
 	github.com/ethereum/go-ethereum v1.8.27
 	github.com/stretchr/testify v1.3.0
 	golang.org/x/crypto v0.0.0-20190621222207-cc06ce4a13d4
+	golang.org/x/sys v0.0.0-20190412213103-97732733099d
 )

mimc7/mimc7.go

@@ -6,7 +6,7 @@ import (
 
 	"github.com/ethereum/go-ethereum/crypto"
 	_constants "github.com/iden3/go-iden3-crypto/constants"
-	"github.com/iden3/go-iden3-crypto/field"
+	"github.com/iden3/go-iden3-crypto/ff"
 	"github.com/iden3/go-iden3-crypto/utils"
 )
 
@@ -15,73 +15,72 @@ const SEED = "mimc"
 var constants = generateConstantsData()
 
 type constantsData struct {
-	maxFieldVal *big.Int
-	seedHash    *big.Int
-	iv          *big.Int
-	fqR         field.Fq
-	nRounds     int
-	cts         []*big.Int
+	seedHash *big.Int
+	iv       *big.Int
+	nRounds  int
+	cts      []*ff.Element
 }
 
 func generateConstantsData() constantsData {
 	var constants constantsData
 
-	fqR := field.NewFq(_constants.Q)
-	constants.fqR = fqR
-
-	// maxFieldVal is the R value of the Finite Field
-	constants.maxFieldVal = constants.fqR.Q
-
 	constants.seedHash = new(big.Int).SetBytes(crypto.Keccak256([]byte(SEED)))
 	c := new(big.Int).SetBytes(crypto.Keccak256([]byte(SEED + "_iv")))
-	constants.iv = new(big.Int).Mod(c, constants.maxFieldVal)
+	constants.iv = new(big.Int).Mod(c, _constants.Q)
 
 	constants.nRounds = 91
-	cts := getConstants(constants.fqR, SEED, constants.nRounds)
+	cts := getConstants(SEED, constants.nRounds)
 	constants.cts = cts
 	return constants
 }
 
-func getConstants(fqR field.Fq, seed string, nRounds int) []*big.Int {
-	cts := make([]*big.Int, nRounds)
-	cts[0] = big.NewInt(int64(0))
+func getConstants(seed string, nRounds int) []*ff.Element {
+	cts := make([]*ff.Element, nRounds)
+	cts[0] = ff.NewElement()
 	c := new(big.Int).SetBytes(crypto.Keccak256([]byte(SEED)))
 	for i := 1; i < nRounds; i++ {
 		c = new(big.Int).SetBytes(crypto.Keccak256(c.Bytes()))
 
-		n := fqR.Affine(c)
-		cts[i] = n
+		n := new(big.Int).Mod(c, _constants.Q)
+		cts[i] = ff.NewElement().SetBigInt(n)
 	}
 	return cts
 }
 
 // MIMC7HashGeneric performs the MIMC7 hash over a *big.Int, in a generic way, where it can be specified the Finite Field over R, and the number of rounds
-func MIMC7HashGeneric(fqR field.Fq, xIn, k *big.Int, nRounds int) *big.Int {
-	cts := getConstants(fqR, SEED, nRounds)
-	var r *big.Int
+func MIMC7HashGeneric(xInBI, kBI *big.Int, nRounds int) *big.Int {
+	xIn := ff.NewElement().SetBigInt(xInBI)
+	k := ff.NewElement().SetBigInt(kBI)
+
+	cts := getConstants(SEED, nRounds)
+	var r *ff.Element
 	for i := 0; i < nRounds; i++ {
-		var t *big.Int
+		var t *ff.Element
 		if i == 0 {
-			t = fqR.Add(xIn, k)
+			t = ff.NewElement().Add(xIn, k)
 		} else {
-			t = fqR.Add(fqR.Add(r, k), cts[i])
+			t = ff.NewElement().Add(ff.NewElement().Add(r, k), cts[i])
 		}
-		t2 := fqR.Square(t)
-		t4 := fqR.Square(t2)
-		r = fqR.Mul(fqR.Mul(t4, t2), t)
+		t2 := ff.NewElement().Square(t)
+		t4 := ff.NewElement().Square(t2)
+		r = ff.NewElement().Mul(ff.NewElement().Mul(t4, t2), t)
 	}
-	return fqR.Affine(fqR.Add(r, k))
+	rE := ff.NewElement().Add(r, k)
+
+	res := big.NewInt(0)
+	rE.ToBigIntRegular(res)
+	return res
 }
 
 // HashGeneric performs the MIMC7 hash over a *big.Int array, in a generic way, where it can be specified the Finite Field over R, and the number of rounds
-func HashGeneric(iv *big.Int, arr []*big.Int, fqR field.Fq, nRounds int) (*big.Int, error) {
+func HashGeneric(iv *big.Int, arr []*big.Int, nRounds int) (*big.Int, error) {
 	if !utils.CheckBigIntArrayInField(arr) {
 		return nil, errors.New("inputs values not inside Finite Field")
 	}
 	r := iv
 	var err error
 	for i := 0; i < len(arr); i++ {
-		r = MIMC7HashGeneric(fqR, r, arr[i], nRounds)
+		r = MIMC7HashGeneric(r, arr[i], nRounds)
 		if err != nil {
 			return r, err
 		}
@@ -90,20 +89,27 @@ func HashGeneric(iv *big.Int, arr []*big.Int, fqR field.Fq, nRounds int) (*big.I
 }
 
 // MIMC7Hash performs the MIMC7 hash over a *big.Int, using the Finite Field over R and the number of rounds setted in the `constants` variable
-func MIMC7Hash(xIn, k *big.Int) *big.Int {
-	var r *big.Int
+func MIMC7Hash(xInBI, kBI *big.Int) *big.Int {
+	xIn := ff.NewElement().SetBigInt(xInBI)
+	k := ff.NewElement().SetBigInt(kBI)
+
+	var r *ff.Element
 	for i := 0; i < constants.nRounds; i++ {
-		var t *big.Int
+		var t *ff.Element
 		if i == 0 {
-			t = constants.fqR.Add(xIn, k)
+			t = ff.NewElement().Add(xIn, k)
 		} else {
-			t = constants.fqR.Add(constants.fqR.Add(r, k), constants.cts[i])
+			t = ff.NewElement().Add(ff.NewElement().Add(r, k), constants.cts[i])
 		}
-		t2 := constants.fqR.Square(t)
-		t4 := constants.fqR.Square(t2)
-		r = constants.fqR.Mul(constants.fqR.Mul(t4, t2), t)
+		t2 := ff.NewElement().Square(t)
+		t4 := ff.NewElement().Square(t2)
+		r = ff.NewElement().Mul(ff.NewElement().Mul(t4, t2), t)
 	}
-	return constants.fqR.Affine(constants.fqR.Add(r, k))
+	rE := ff.NewElement().Add(r, k)
+
+	res := big.NewInt(0)
+	rE.ToBigIntRegular(res)
+	return res
 }
 
 // Hash performs the MIMC7 hash over a *big.Int array
@@ -113,17 +119,18 @@ func Hash(arr []*big.Int, key *big.Int) (*big.Int, error) {
 	}
 	var r *big.Int
 	if key == nil {
-		r = constants.fqR.Zero()
+		r = big.NewInt(0)
 	} else {
 		r = key
 	}
 	for i := 0; i < len(arr); i++ {
-		r = constants.fqR.Add(
-			constants.fqR.Add(
+		r = new(big.Int).Add(
+			new(big.Int).Add(
 				r,
 				arr[i],
 			),
 			MIMC7Hash(arr[i], r))
+		r = new(big.Int).Mod(r, _constants.Q)
 	}
 	return r, nil
 }

mimc7/mimc7_test.go

@@ -6,7 +6,6 @@ import (
 	"testing"
 
 	"github.com/ethereum/go-ethereum/crypto"
-	"github.com/iden3/go-iden3-crypto/field"
 	"github.com/stretchr/testify/assert"
 )
 
@@ -22,16 +21,12 @@ func TestMIMC7Generic(t *testing.T) {
 	b2 := big.NewInt(int64(2))
 	b3 := big.NewInt(int64(3))
 
-	r, ok := new(big.Int).SetString("21888242871839275222246405745257275088548364400416034343698204186575808495617", 10)
-	assert.True(t, ok)
-	fqR := field.NewFq(r)
-
 	bigArray := []*big.Int{b1, b2, b3}
 
 	// Generic Hash
-	mhg := MIMC7HashGeneric(fqR, b1, b2, 91)
+	mhg := MIMC7HashGeneric(b1, b2, 91)
 	assert.Equal(t, "10594780656576967754230020536574539122676596303354946869887184401991294982664", mhg.String())
-	hg, err := HashGeneric(fqR.Zero(), bigArray, fqR, 91)
+	hg, err := HashGeneric(big.NewInt(0), bigArray, 91)
 	assert.Nil(t, err)
 	assert.Equal(t, "6464402164086696096195815557694604139393321133243036833927490113253119343397", (*big.Int)(hg).String())
 }

poseidon/poseidon.go

@@ -20,7 +20,7 @@ var constC []*ff.Element
 var constM [T][T]*ff.Element
 
 func Zero() *ff.Element {
-	return ff.NewElement().SetZero()
+	return ff.NewElement()
 }
 
 func modQ(v *big.Int) {
@@ -71,7 +71,7 @@ func getMDS() [T][T]*ff.Element {
 
 func checkAllDifferent(v []*ff.Element) bool {
 	for i := 0; i < len(v); i++ {
-		if v[i].Equal(ff.NewElement().SetZero()) {
+		if v[i].Equal(ff.NewElement()) {
 			return false
 		}
 		for j := i + 1; j < len(v); j++ {
@@ -92,7 +92,6 @@ func ark(state [T]*ff.Element, c *ff.Element) {
 
 // cubic performs x^5 mod p
 // https://eprint.iacr.org/2019/458.pdf page 8
-// var five = big.NewInt(5)
 
 func cubic(a *ff.Element) {
 	a.Exp(*a, 5)