Merge pull request #46 from iden3/feature/golden

Golden poseidon hash implementation
2026-02-07 11:36:41 +01:00 · 2022-03-22 05:06:23 -04:00
parent f4972de131 5848bf2918
commit d3e4218fe3
12 changed files with 3221 additions and 0 deletions
--- a/ffg/arith.go
+++ b/ffg/arith.go
@@ -0,0 +1,30 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+package ffg
+
+import (
+	"math/bits"
+)
+
+// madd0 hi = a*b + c (discards lo bits)
+func madd0(a, b, c uint64) (hi uint64) {
+	var carry, lo uint64
+	hi, lo = bits.Mul64(a, b)
+	_, carry = bits.Add64(lo, c, 0)
+	hi, _ = bits.Add64(hi, 0, carry)
+	return
+}
--- a/ffg/asm.go
+++ b/ffg/asm.go
@@ -0,0 +1,24 @@
+//go:build !noadx
+// +build !noadx
+
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+package ffg
+
+import "golang.org/x/sys/cpu"
+
+var supportAdx = cpu.X86.HasADX && cpu.X86.HasBMI2
--- a/ffg/asm_noadx.go
+++ b/ffg/asm_noadx.go
@@ -0,0 +1,25 @@
+//go:build noadx
+// +build noadx
+
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+package ffg
+
+// note: this is needed for test purposes, as dynamically changing supportAdx doesn't flag
+// certain errors (like fatal error: missing stackmap)
+// this ensures we test all asm path.
+var supportAdx = false
--- a/ffg/doc.go
+++ b/ffg/doc.go
@@ -0,0 +1,43 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+// Package ffg contains field arithmetic operations for modulus = 0xffffff...000001.
+//
+// The API is similar to math/big (big.Int), but the operations are significantly faster (up to 20x for the modular multiplication on amd64, see also https://hackmd.io/@zkteam/modular_multiplication)
+//
+// The modulus is hardcoded in all the operations.
+//
+// Field elements are represented as an array, and assumed to be in Montgomery form in all methods:
+// 	type Element [1]uint64
+//
+// Example API signature
+// 	// Mul z = x * y mod q
+// 	func (z *Element) Mul(x, y *Element) *Element
+//
+// and can be used like so:
+// 	var a, b Element
+// 	a.SetUint64(2)
+// 	b.SetString("984896738")
+// 	a.Mul(a, b)
+// 	a.Sub(a, a)
+// 	 .Add(a, b)
+// 	 .Inv(a)
+// 	b.Exp(b, new(big.Int).SetUint64(42))
+//
+// Modulus
+// 	0xffffffff00000001 // base 16
+// 	18446744069414584321 // base 10
+package ffg
--- a/ffg/element.go
+++ b/ffg/element.go
@@ -0,0 +1,759 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+package ffg
+
+// /!\ 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"
+	"errors"
+	"io"
+	"math/big"
+	"math/bits"
+	"reflect"
+	"strconv"
+	"sync"
+)
+
+// Element represents a field element stored on 1 words (uint64)
+// Element are assumed to be in Montgomery form in all methods
+// field modulus q =
+//
+// 18446744069414584321
+type Element [1]uint64
+
+// Limbs number of 64 bits words needed to represent Element
+const Limbs = 1
+
+// Bits number bits needed to represent Element
+const Bits = 64
+
+// Bytes number bytes needed to represent Element
+const Bytes = Limbs * 8
+
+// field modulus stored as big.Int
+var _modulus big.Int
+
+// Modulus returns q as a big.Int
+// q =
+//
+// 18446744069414584321
+func Modulus() *big.Int {
+	return new(big.Int).Set(&_modulus)
+}
+
+// q (modulus)
+var qElement = Element{
+	18446744069414584321,
+}
+
+// rSquare
+var rSquare = Element{
+	18446744065119617025,
+}
+
+var bigIntPool = sync.Pool{
+	New: func() interface{} {
+		return new(big.Int)
+	},
+}
+
+func init() {
+	_modulus.SetString("18446744069414584321", 10)
+}
+
+// NewElement returns a new Element
+func NewElement() *Element {
+	return &Element{}
+}
+
+// NewElementFromUint64 returns a new Element from a uint64 value
+//
+// it is equivalent to
+// 		var v NewElement
+// 		v.SetUint64(...)
+func NewElementFromUint64(v uint64) *Element {
+	z := Element{v}
+	z.Mul(&z, &rSquare)
+	return &z
+}
+
+// SetUint64 z = v, sets z LSB to v (non-Montgomery form) and convert z to Montgomery form
+func (z *Element) SetUint64(v uint64) *Element {
+	*z = Element{v}
+	return z.Mul(z, &rSquare) // z.ToMont()
+}
+
+// Set z = x
+func (z *Element) Set(x *Element) *Element {
+	z[0] = x[0]
+	return z
+}
+
+// SetInterface converts provided interface into Element
+// returns an error if provided type is not supported
+// supported types: Element, *Element, uint64, int, string (interpreted as base10 integer),
+// *big.Int, big.Int, []byte
+func (z *Element) SetInterface(i1 interface{}) (*Element, error) {
+	switch c1 := i1.(type) {
+	case Element:
+		return z.Set(&c1), nil
+	case *Element:
+		return z.Set(c1), nil
+	case uint64:
+		return z.SetUint64(c1), nil
+	case int:
+		return z.SetString(strconv.Itoa(c1)), nil
+	case string:
+		return z.SetString(c1), nil
+	case *big.Int:
+		return z.SetBigInt(c1), nil
+	case big.Int:
+		return z.SetBigInt(&c1), nil
+	case []byte:
+		return z.SetBytes(c1), nil
+	default:
+		return nil, errors.New("can't set ffg.Element from type " + reflect.TypeOf(i1).String())
+	}
+}
+
+// SetZero z = 0
+func (z *Element) SetZero() *Element {
+	z[0] = 0
+	return z
+}
+
+// SetOne z = 1 (in Montgomery form)
+func (z *Element) SetOne() *Element {
+	z[0] = 4294967295
+	return z
+}
+
+// Div z = x*y^-1 mod q
+func (z *Element) Div(x, y *Element) *Element {
+	var yInv Element
+	yInv.Inverse(y)
+	z.Mul(x, &yInv)
+	return z
+}
+
+// Bit returns the i'th bit, with lsb == bit 0.
+// It is the responsability of the caller to convert from Montgomery to Regular form if needed
+func (z *Element) Bit(i uint64) uint64 {
+	j := i / 64
+	if j >= 1 {
+		return 0
+	}
+	return uint64(z[j] >> (i % 64) & 1)
+}
+
+// Equal returns z == x
+func (z *Element) Equal(x *Element) bool {
+	return (z[0] == x[0])
+}
+
+// IsZero returns z == 0
+func (z *Element) IsZero() bool {
+	return (z[0]) == 0
+}
+
+// IsUint64 returns true if z[0] >= 0 and all other words are 0
+func (z *Element) IsUint64() bool {
+	return true
+}
+
+// Cmp compares (lexicographic order) z and x and returns:
+//
+//   -1 if z <  x
+//    0 if z == x
+//   +1 if z >  x
+//
+func (z *Element) Cmp(x *Element) int {
+	_z := *z
+	_x := *x
+	_z.FromMont()
+	_x.FromMont()
+	if _z[0] > _x[0] {
+		return 1
+	} else if _z[0] < _x[0] {
+		return -1
+	}
+	return 0
+}
+
+// LexicographicallyLargest returns true if this element is strictly lexicographically
+// larger than its negation, false otherwise
+func (z *Element) LexicographicallyLargest() bool {
+	// adapted from github.com/zkcrypto/bls12_381
+	// we check if the element is larger than (q-1) / 2
+	// if z - (((q -1) / 2) + 1) have no underflow, then z > (q-1) / 2
+
+	_z := *z
+	_z.FromMont()
+
+	var b uint64
+	_, b = bits.Sub64(_z[0], 9223372034707292161, 0)
+
+	return b == 0
+}
+
+// SetRandom sets z to a random element < q
+func (z *Element) SetRandom() (*Element, error) {
+	var bytes [8]byte
+	if _, err := io.ReadFull(rand.Reader, bytes[:]); err != nil {
+		return nil, err
+	}
+	z[0] = binary.BigEndian.Uint64(bytes[0:8])
+	z[0] %= 18446744069414584321
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[0] < 18446744069414584321) {
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+	}
+
+	return z, nil
+}
+
+// One returns 1 (in montgommery form)
+func One() Element {
+	var one Element
+	one.SetOne()
+	return one
+}
+
+// Halve sets z to z / 2 (mod p)
+func (z *Element) Halve() {
+	var twoInv Element
+	twoInv.SetOne().Double(&twoInv).Inverse(&twoInv)
+	z.Mul(z, &twoInv)
+
+}
+
+// API with assembly impl
+
+// Mul z = x * y mod q
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) Mul(x, y *Element) *Element {
+	mul(z, x, y)
+	return z
+}
+
+// Square z = x * x mod q
+// see https://hackmd.io/@zkteam/modular_multiplication
+func (z *Element) Square(x *Element) *Element {
+	mul(z, x, x)
+	return z
+}
+
+// FromMont converts z in place (i.e. mutates) from Montgomery to regular representation
+// sets and returns z = z * 1
+func (z *Element) FromMont() *Element {
+	fromMont(z)
+	return z
+}
+
+// Add z = x + y mod q
+func (z *Element) Add(x, y *Element) *Element {
+	add(z, x, y)
+	return z
+}
+
+// Double z = x + x mod q, aka Lsh 1
+func (z *Element) Double(x *Element) *Element {
+	double(z, x)
+	return z
+}
+
+// Sub  z = x - y mod q
+func (z *Element) Sub(x, y *Element) *Element {
+	sub(z, x, y)
+	return z
+}
+
+// Neg z = q - x
+func (z *Element) Neg(x *Element) *Element {
+	neg(z, x)
+	return z
+}
+
+// Generic (no ADX instructions, no AMD64) versions of multiplication and squaring algorithms
+
+func _mulGeneric(z, x, y *Element) {
+
+	var t [2]uint64
+	var D uint64
+	var m, C uint64
+	// -----------------------------------
+	// First loop
+
+	C, t[0] = bits.Mul64(y[0], x[0])
+
+	t[1], D = bits.Add64(t[1], C, 0)
+
+	// m = t[0]n'[0] mod W
+	m = t[0] * 18446744069414584319
+
+	// -----------------------------------
+	// Second loop
+	C = madd0(m, 18446744069414584321, t[0])
+
+	t[0], C = bits.Add64(t[1], C, 0)
+	t[1], _ = bits.Add64(0, D, C)
+
+	if t[1] != 0 {
+		// we need to reduce, we have a result on 2 words
+		z[0], _ = bits.Sub64(t[0], 18446744069414584321, 0)
+
+		return
+
+	}
+
+	// copy t into z
+	z[0] = t[0]
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[0] < 18446744069414584321) {
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+	}
+}
+
+func _fromMontGeneric(z *Element) {
+	// the following lines implement z = z * 1
+	// with a modified CIOS montgomery multiplication
+	{
+		// m = z[0]n'[0] mod W
+		m := z[0] * 18446744069414584319
+		C := madd0(m, 18446744069414584321, z[0])
+		z[0] = C
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[0] < 18446744069414584321) {
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+	}
+}
+
+func _addGeneric(z, x, y *Element) {
+	var carry uint64
+
+	z[0], carry = bits.Add64(x[0], y[0], 0)
+	// if we overflowed the last addition, z >= q
+	// if z >= q, z = z - q
+	if carry != 0 {
+		// we overflowed, so z >= q
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+		return
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[0] < 18446744069414584321) {
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+	}
+}
+
+func _doubleGeneric(z, x *Element) {
+	var carry uint64
+
+	z[0], carry = bits.Add64(x[0], x[0], 0)
+	// if we overflowed the last addition, z >= q
+	// if z >= q, z = z - q
+	if carry != 0 {
+		// we overflowed, so z >= q
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+		return
+	}
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[0] < 18446744069414584321) {
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+	}
+}
+
+func _subGeneric(z, x, y *Element) {
+	var b uint64
+	z[0], b = bits.Sub64(x[0], y[0], 0)
+	if b != 0 {
+		z[0], _ = bits.Add64(z[0], 18446744069414584321, 0)
+	}
+}
+
+func _negGeneric(z, x *Element) {
+	if x.IsZero() {
+		z.SetZero()
+		return
+	}
+	z[0], _ = bits.Sub64(18446744069414584321, x[0], 0)
+}
+
+func _reduceGeneric(z *Element) {
+
+	// if z > q --> z -= q
+	// note: this is NOT constant time
+	if !(z[0] < 18446744069414584321) {
+		z[0], _ = bits.Sub64(z[0], 18446744069414584321, 0)
+	}
+}
+
+func mulByConstant(z *Element, c uint8) {
+	switch c {
+	case 0:
+		z.SetZero()
+		return
+	case 1:
+		return
+	case 2:
+		z.Double(z)
+		return
+	case 3:
+		_z := *z
+		z.Double(z).Add(z, &_z)
+	case 5:
+		_z := *z
+		z.Double(z).Double(z).Add(z, &_z)
+	default:
+		var y Element
+		y.SetUint64(uint64(c))
+		z.Mul(z, &y)
+	}
+}
+
+// BatchInvert returns a new slice with every element inverted.
+// Uses Montgomery batch inversion trick
+func BatchInvert(a []Element) []Element {
+	res := make([]Element, len(a))
+	if len(a) == 0 {
+		return res
+	}
+
+	zeroes := make([]bool, len(a))
+	accumulator := One()
+
+	for i := 0; i < len(a); i++ {
+		if a[i].IsZero() {
+			zeroes[i] = true
+			continue
+		}
+		res[i] = accumulator
+		accumulator.Mul(&accumulator, &a[i])
+	}
+
+	accumulator.Inverse(&accumulator)
+
+	for i := len(a) - 1; i >= 0; i-- {
+		if zeroes[i] {
+			continue
+		}
+		res[i].Mul(&res[i], &accumulator)
+		accumulator.Mul(&accumulator, &a[i])
+	}
+
+	return res
+}
+
+func _butterflyGeneric(a, b *Element) {
+	t := *a
+	a.Add(a, b)
+	b.Sub(&t, b)
+}
+
+// BitLen returns the minimum number of bits needed to represent z
+// returns 0 if z == 0
+func (z *Element) BitLen() int {
+	return bits.Len64(z[0])
+}
+
+// Exp z = x^exponent mod q
+func (z *Element) Exp(x Element, exponent *big.Int) *Element {
+	var bZero big.Int
+	if exponent.Cmp(&bZero) == 0 {
+		return z.SetOne()
+	}
+
+	z.Set(&x)
+
+	for i := exponent.BitLen() - 2; i >= 0; i-- {
+		z.Square(z)
+		if exponent.Bit(i) == 1 {
+			z.Mul(z, &x)
+		}
+	}
+
+	return z
+}
+
+// ToMont converts z to Montgomery form
+// sets and returns z = z * r^2
+func (z *Element) ToMont() *Element {
+	return z.Mul(z, &rSquare)
+}
+
+// ToRegular returns z in regular form (doesn't mutate z)
+func (z Element) ToRegular() Element {
+	return *z.FromMont()
+}
+
+// String returns the string form of an Element in Montgomery form
+func (z *Element) String() string {
+	zz := *z
+	zz.FromMont()
+	if zz.IsUint64() {
+		return strconv.FormatUint(zz[0], 10)
+	} else {
+		var zzNeg Element
+		zzNeg.Neg(z)
+		zzNeg.FromMont()
+		if zzNeg.IsUint64() {
+			return "-" + strconv.FormatUint(zzNeg[0], 10)
+		}
+	}
+	vv := bigIntPool.Get().(*big.Int)
+	defer bigIntPool.Put(vv)
+	return zz.ToBigInt(vv).String()
+}
+
+// ToBigInt returns z as a big.Int in Montgomery form
+func (z *Element) ToBigInt(res *big.Int) *big.Int {
+	var b [Limbs * 8]byte
+	binary.BigEndian.PutUint64(b[0:8], z[0])
+
+	return res.SetBytes(b[:])
+}
+
+// ToBigIntRegular returns z as a big.Int in regular form
+func (z Element) ToBigIntRegular(res *big.Int) *big.Int {
+	z.FromMont()
+	return z.ToBigInt(res)
+}
+
+// ToUint64Regular returns z as a uint64 in regular form
+func (z Element) ToUint64Regular() uint64 {
+	z.FromMont()
+	return z[0]
+}
+
+// Bytes returns the regular (non montgomery) value
+// of z as a big-endian byte array.
+func (z *Element) Bytes() (res [Limbs * 8]byte) {
+	_z := z.ToRegular()
+	binary.BigEndian.PutUint64(res[0:8], _z[0])
+
+	return
+}
+
+// Marshal returns the regular (non montgomery) value
+// of z as a big-endian byte slice.
+func (z *Element) Marshal() []byte {
+	b := z.Bytes()
+	return b[:]
+}
+
+// SetBytes interprets e as the bytes of a big-endian unsigned integer,
+// sets z to that value (in Montgomery form), and returns z.
+func (z *Element) SetBytes(e []byte) *Element {
+	// get a big int from our pool
+	vv := bigIntPool.Get().(*big.Int)
+	vv.SetBytes(e)
+
+	// set big int
+	z.SetBigInt(vv)
+
+	// put temporary object back in pool
+	bigIntPool.Put(vv)
+
+	return z
+}
+
+// SetBigInt sets z to v (regular form) and returns z in Montgomery form
+func (z *Element) SetBigInt(v *big.Int) *Element {
+	z.SetZero()
+
+	var zero big.Int
+
+	// fast path
+	c := v.Cmp(&_modulus)
+	if c == 0 {
+		// v == 0
+		return z
+	} else if c != 1 && v.Cmp(&zero) != -1 {
+		// 0 < v < q
+		return z.setBigInt(v)
+	}
+
+	// get temporary big int from the pool
+	vv := bigIntPool.Get().(*big.Int)
+
+	// copy input + modular reduction
+	vv.Set(v)
+	vv.Mod(v, &_modulus)
+
+	// set big int byte value
+	z.setBigInt(vv)
+
+	// release object into pool
+	bigIntPool.Put(vv)
+	return z
+}
+
+// setBigInt assumes 0 <= v < q
+func (z *Element) setBigInt(v *big.Int) *Element {
+	vBits := v.Bits()
+
+	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()
+}
+
+// SetString creates a big.Int with s (in base 10) and calls SetBigInt on z
+func (z *Element) SetString(s string) *Element {
+	// get temporary big int from the pool
+	vv := bigIntPool.Get().(*big.Int)
+
+	if _, ok := vv.SetString(s, 10); !ok {
+		panic("Element.SetString failed -> can't parse number in base10 into a big.Int")
+	}
+	z.SetBigInt(vv)
+
+	// release object into pool
+	bigIntPool.Put(vv)
+
+	return z
+}
+
+var (
+	_bLegendreExponentElement *big.Int
+	_bSqrtExponentElement     *big.Int
+)
+
+func init() {
+	_bLegendreExponentElement, _ = new(big.Int).SetString("7fffffff80000000", 16)
+	const sqrtExponentElement = "7fffffff"
+	_bSqrtExponentElement, _ = new(big.Int).SetString(sqrtExponentElement, 16)
+}
+
+// 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, _bLegendreExponentElement)
+
+	if l.IsZero() {
+		return 0
+	}
+
+	// if l == 1
+	if l[0] == 4294967295 {
+		return 1
+	}
+	return -1
+}
+
+// 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 y, b, t, w Element
+	// w = x^((s-1)/2))
+	w.Exp(*x, _bSqrtExponentElement)
+
+	// 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{
+		15733474329512464024,
+	}
+	r := uint64(32)
+
+	// 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[0] == 4294967295) {
+		// t != 1, we don't have a square root
+		return nil
+	}
+	for {
+		var m uint64
+		t = b
+
+		// for t != 1
+		for !(t[0] == 4294967295) {
+			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.Mul(&y, &t)
+		b.Mul(&b, &g)
+		r = m
+	}
+}
+
+// Inverse z = x^-1 mod q
+// note: allocates a big.Int (math/big)
+func (z *Element) Inverse(x *Element) *Element {
+	var _xNonMont big.Int
+	x.ToBigIntRegular(&_xNonMont)
+	_xNonMont.ModInverse(&_xNonMont, Modulus())
+	z.SetBigInt(&_xNonMont)
+	return z
+}
--- a/ffg/element_fuzz.go
+++ b/ffg/element_fuzz.go
@@ -0,0 +1,113 @@
+//go:build gofuzz
+// +build gofuzz
+
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+package ffg
+
+import (
+	"bytes"
+	"encoding/binary"
+	"io"
+	"math/big"
+	"math/bits"
+)
+
+const (
+	fuzzInteresting = 1
+	fuzzNormal      = 0
+	fuzzDiscard     = -1
+)
+
+// Fuzz arithmetic operations fuzzer
+func Fuzz(data []byte) int {
+	r := bytes.NewReader(data)
+
+	var e1, e2 Element
+	e1.SetRawBytes(r)
+	e2.SetRawBytes(r)
+
+	{
+		// mul assembly
+
+		var c, _c Element
+		a, _a, b, _b := e1, e1, e2, e2
+		c.Mul(&a, &b)
+		_mulGeneric(&_c, &_a, &_b)
+
+		if !c.Equal(&_c) {
+			panic("mul asm != mul generic on Element")
+		}
+	}
+
+	{
+		// inverse
+		inv := e1
+		inv.Inverse(&inv)
+
+		var bInv, b1, b2 big.Int
+		e1.ToBigIntRegular(&b1)
+		bInv.ModInverse(&b1, Modulus())
+		inv.ToBigIntRegular(&b2)
+
+		if b2.Cmp(&bInv) != 0 {
+			panic("inverse operation doesn't match big int result")
+		}
+	}
+
+	{
+		// a + -a == 0
+		a, b := e1, e1
+		b.Neg(&b)
+		a.Add(&a, &b)
+		if !a.IsZero() {
+			panic("a + -a != 0")
+		}
+	}
+
+	return fuzzNormal
+
+}
+
+// SetRawBytes reads up to Bytes (bytes needed to represent Element) from reader
+// and interpret it as big endian uint64
+// used for fuzzing purposes only
+func (z *Element) SetRawBytes(r io.Reader) {
+
+	buf := make([]byte, 8)
+
+	for i := 0; i < len(z); i++ {
+		if _, err := io.ReadFull(r, buf); err != nil {
+			goto eof
+		}
+		z[i] = binary.BigEndian.Uint64(buf[:])
+	}
+eof:
+	z[0] %= qElement[0]
+
+	if z.BiggerModulus() {
+		var b uint64
+		z[0], b = bits.Sub64(z[0], qElement[0], 0)
+	}
+
+	return
+}
+
+func (z *Element) BiggerModulus() bool {
+
+	return z[0] >= qElement[0]
+}
--- a/ffg/element_ops_amd64.go
+++ b/ffg/element_ops_amd64.go
@@ -0,0 +1,17 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+package ffg
--- a/ffg/element_ops_noasm.go
+++ b/ffg/element_ops_noasm.go
@@ -0,0 +1,75 @@
+// Copyright 2020 ConsenSys Software Inc.
+//
+// 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 consensys/gnark-crypto DO NOT EDIT
+
+package ffg
+
+// /!\ 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 /!\
+
+// MulBy3 x *= 3
+func MulBy3(x *Element) {
+	mulByConstant(x, 3)
+}
+
+// MulBy5 x *= 5
+func MulBy5(x *Element) {
+	mulByConstant(x, 5)
+}
+
+// MulBy13 x *= 13
+func MulBy13(x *Element) {
+	mulByConstant(x, 13)
+}
+
+// Butterfly sets
+// a = a + b
+// b = a - b
+func Butterfly(a, b *Element) {
+	_butterflyGeneric(a, b)
+}
+
+func mul(z, x, y *Element) {
+	_mulGeneric(z, x, y)
+}
+
+// FromMont converts z in place (i.e. mutates) from Montgomery to regular representation
+// sets and returns z = z * 1
+func fromMont(z *Element) {
+	_fromMontGeneric(z)
+}
+
+func add(z, x, y *Element) {
+	_addGeneric(z, x, y)
+}
+
+func double(z, x *Element) {
+	_doubleGeneric(z, x)
+}
+
+func sub(z, x, y *Element) {
+	_subGeneric(z, x, y)
+}
+
+func neg(z, x *Element) {
+	_negGeneric(z, x)
+}
+
+func reduce(z *Element) {
+	_reduceGeneric(z)
+}
--- a/ffg/element_test.go
+++ b/ffg/element_test.go
--- a/goldenposeidon/constants.go
+++ b/goldenposeidon/constants.go
@@ -0,0 +1,131 @@
+package poseidon
+
+import "github.com/iden3/go-iden3-crypto/ffg"
+
+const (
+	NROUNDSF = 8  //nolint:golint
+	NROUNDSP = 22 //nolint:golint
+	CAPLEN   = 4  //nolint:golint
+	mLen     = 12
+)
+
+var (
+	mcirc = []uint64{17, 15, 41, 16, 2, 28, 13, 13, 39, 18, 34, 20}
+	mdiag = []uint64{8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
+
+	c = []uint64{
+		0xb585f766f2144405, 0x7746a55f43921ad7, 0xb2fb0d31cee799b4, 0x0f6760a4803427d7,
+		0xe10d666650f4e012, 0x8cae14cb07d09bf1, 0xd438539c95f63e9f, 0xef781c7ce35b4c3d,
+		0xcdc4a239b0c44426, 0x277fa208bf337bff, 0xe17653a29da578a1, 0xc54302f225db2c76,
+		0x86287821f722c881, 0x59cd1a8a41c18e55, 0xc3b919ad495dc574, 0xa484c4c5ef6a0781,
+		0x308bbd23dc5416cc, 0x6e4a40c18f30c09c, 0x9a2eedb70d8f8cfa, 0xe360c6e0ae486f38,
+		0xd5c7718fbfc647fb, 0xc35eae071903ff0b, 0x849c2656969c4be7, 0xc0572c8c08cbbbad,
+		0xe9fa634a21de0082, 0xf56f6d48959a600d, 0xf7d713e806391165, 0x8297132b32825daf,
+		0xad6805e0e30b2c8a, 0xac51d9f5fcf8535e, 0x502ad7dc18c2ad87, 0x57a1550c110b3041,
+		0x66bbd30e6ce0e583, 0x0da2abef589d644e, 0xf061274fdb150d61, 0x28b8ec3ae9c29633,
+		0x92a756e67e2b9413, 0x70e741ebfee96586, 0x019d5ee2af82ec1c, 0x6f6f2ed772466352,
+		0x7cf416cfe7e14ca1, 0x61df517b86a46439, 0x85dc499b11d77b75, 0x4b959b48b9c10733,
+		0xe8be3e5da8043e57, 0xf5c0bc1de6da8699, 0x40b12cbf09ef74bf, 0xa637093ecb2ad631,
+		0x3cc3f892184df408, 0x2e479dc157bf31bb, 0x6f49de07a6234346, 0x213ce7bede378d7b,
+		0x5b0431345d4dea83, 0xa2de45780344d6a1, 0x7103aaf94a7bf308, 0x5326fc0d97279301,
+		0xa9ceb74fec024747, 0x27f8ec88bb21b1a3, 0xfceb4fda1ded0893, 0xfac6ff1346a41675,
+		0x7131aa45268d7d8c, 0x9351036095630f9f, 0xad535b24afc26bfb, 0x4627f5c6993e44be,
+		0x645cf794b8f1cc58, 0x241c70ed0af61617, 0xacb8e076647905f1, 0x3737e9db4c4f474d,
+		0xe7ea5e33e75fffb6, 0x90dee49fc9bfc23a, 0xd1b1edf76bc09c92, 0x0b65481ba645c602,
+		0x99ad1aab0814283b, 0x438a7c91d416ca4d, 0xb60de3bcc5ea751c, 0xc99cab6aef6f58bc,
+		0x69a5ed92a72ee4ff, 0x5e7b329c1ed4ad71, 0x5fc0ac0800144885, 0x32db829239774eca,
+		0x0ade699c5830f310, 0x7cc5583b10415f21, 0x85df9ed2e166d64f, 0x6604df4fee32bcb1,
+		0xeb84f608da56ef48, 0xda608834c40e603d, 0x8f97fe408061f183, 0xa93f485c96f37b89,
+		0x6704e8ee8f18d563, 0xcee3e9ac1e072119, 0x510d0e65e2b470c1, 0xf6323f486b9038f0,
+		0x0b508cdeffa5ceef, 0xf2417089e4fb3cbd, 0x60e75c2890d15730, 0xa6217d8bf660f29c,
+		0x7159cd30c3ac118e, 0x839b4e8fafead540, 0x0d3f3e5e82920adc, 0x8f7d83bddee7bba8,
+		0x780f2243ea071d06, 0xeb915845f3de1634, 0xd19e120d26b6f386, 0x016ee53a7e5fecc6,
+		0xcb5fd54e7933e477, 0xacb8417879fd449f, 0x9c22190be7f74732, 0x5d693c1ba3ba3621,
+		0xdcef0797c2b69ec7, 0x3d639263da827b13, 0xe273fd971bc8d0e7, 0x418f02702d227ed5,
+		0x8c25fda3b503038c, 0x2cbaed4daec8c07c, 0x5f58e6afcdd6ddc2, 0x284650ac5e1b0eba,
+		0x635b337ee819dab5, 0x9f9a036ed4f2d49f, 0xb93e260cae5c170e, 0xb0a7eae879ddb76d,
+		0xd0762cbc8ca6570c, 0x34c6efb812b04bf5, 0x40bf0ab5fa14c112, 0xb6b570fc7c5740d3,
+		0x5a27b9002de33454, 0xb1a5b165b6d2b2d2, 0x8722e0ace9d1be22, 0x788ee3b37e5680fb,
+		0x14a726661551e284, 0x98b7672f9ef3b419, 0xbb93ae776bb30e3a, 0x28fd3b046380f850,
+		0x30a4680593258387, 0x337dc00c61bd9ce1, 0xd5eca244c7a4ff1d, 0x7762638264d279bd,
+		0xc1e434bedeefd767, 0x0299351a53b8ec22, 0xb2d456e4ad251b80, 0x3e9ed1fda49cea0b,
+		0x2972a92ba450bed8, 0x20216dd77be493de, 0xadffe8cf28449ec6, 0x1c4dbb1c4c27d243,
+		0x15a16a8a8322d458, 0x388a128b7fd9a609, 0x2300e5d6baedf0fb, 0x2f63aa8647e15104,
+		0xf1c36ce86ecec269, 0x27181125183970c9, 0xe584029370dca96d, 0x4d9bbc3e02f1cfb2,
+		0xea35bc29692af6f8, 0x18e21b4beabb4137, 0x1e3b9fc625b554f4, 0x25d64362697828fd,
+		0x5a3f1bb1c53a9645, 0xdb7f023869fb8d38, 0xb462065911d4e1fc, 0x49c24ae4437d8030,
+		0xd793862c112b0566, 0xaadd1106730d8feb, 0xc43b6e0e97b0d568, 0xe29024c18ee6fca2,
+		0x5e50c27535b88c66, 0x10383f20a4ff9a87, 0x38e8ee9d71a45af8, 0xdd5118375bf1a9b9,
+		0x775005982d74d7f7, 0x86ab99b4dde6c8b0, 0xb1204f603f51c080, 0xef61ac8470250ecf,
+		0x1bbcd90f132c603f, 0x0cd1dabd964db557, 0x11a3ae5beb9d1ec9, 0xf755bfeea585d11d,
+		0xa3b83250268ea4d7, 0x516306f4927c93af, 0xddb4ac49c9efa1da, 0x64bb6dec369d4418,
+		0xf9cc95c22b4c1fcc, 0x08d37f755f4ae9f6, 0xeec49b613478675b, 0xf143933aed25e0b0,
+		0xe4c5dd8255dfc622, 0xe7ad7756f193198e, 0x92c2318b87fff9cb, 0x739c25f8fd73596d,
+		0x5636cac9f16dfed0, 0xdd8f909a938e0172, 0xc6401fe115063f5b, 0x8ad97b33f1ac1455,
+		0x0c49366bb25e8513, 0x0784d3d2f1698309, 0x530fb67ea1809a81, 0x410492299bb01f49,
+		0x139542347424b9ac, 0x9cb0bd5ea1a1115e, 0x02e3f615c38f49a1, 0x985d4f4a9c5291ef,
+		0x775b9feafdcd26e7, 0x304265a6384f0f2d, 0x593664c39773012c, 0x4f0a2e5fb028f2ce,
+		0xdd611f1000c17442, 0xd8185f9adfea4fd0, 0xef87139ca9a3ab1e, 0x3ba71336c34ee133,
+		0x7d3a455d56b70238, 0x660d32e130182684, 0x297a863f48cd1f43, 0x90e0a736a751ebb7,
+		0x549f80ce550c4fd3, 0x0f73b2922f38bd64, 0x16bf1f73fb7a9c3f, 0x6d1f5a59005bec17,
+		0x02ff876fa5ef97c4, 0xc5cb72a2a51159b0, 0x8470f39d2d5c900e, 0x25abb3f1d39fcb76,
+		0x23eb8cc9b372442f, 0xd687ba55c64f6364, 0xda8d9e90fd8ff158, 0xe3cbdc7d2fe45ea7,
+		0xb9a8c9b3aee52297, 0xc0d28a5c10960bd3, 0x45d7ac9b68f71a34, 0xeeb76e397069e804,
+		0x3d06c8bd1514e2d9, 0x9c9c98207cb10767, 0x65700b51aedfb5ef, 0x911f451539869408,
+		0x7ae6849fbc3a0ec6, 0x3bb340eba06afe7e, 0xb46e9d8b682ea65e, 0x8dcf22f9a3b34356,
+		0x77bdaeda586257a7, 0xf19e400a5104d20d, 0xc368a348e46d950f, 0x9ef1cd60e679f284,
+		0xe89cd854d5d01d33, 0x5cd377dc8bb882a2, 0xa7b0fb7883eee860, 0x7684403ec392950d,
+		0x5fa3f06f4fed3b52, 0x8df57ac11bc04831, 0x2db01efa1e1e1897, 0x54846de4aadb9ca2,
+		0xba6745385893c784, 0x541d496344d2c75b, 0xe909678474e687fe, 0xdfe89923f6c9c2ff,
+		0xece5a71e0cfedc75, 0x5ff98fd5d51fe610, 0x83e8941918964615, 0x5922040b47f150c1,
+		0xf97d750e3dd94521, 0x5080d4c2b86f56d7, 0xa7de115b56c78d70, 0x6a9242ac87538194,
+		0xf7856ef7f9173e44, 0x2265fc92feb0dc09, 0x17dfc8e4f7ba8a57, 0x9001a64209f21db8,
+		0x90004c1371b893c5, 0xb932b7cf752e5545, 0xa0b1df81b6fe59fc, 0x8ef1dd26770af2c2,
+		0x0541a4f9cfbeed35, 0x9e61106178bfc530, 0xb3767e80935d8af2, 0x0098d5782065af06,
+		0x31d191cd5c1466c7, 0x410fefafa319ac9d, 0xbdf8f242e316c4ab, 0x9e8cd55b57637ed0,
+		0xde122bebe9a39368, 0x4d001fd58f002526, 0xca6637000eb4a9f8, 0x2f2339d624f91f78,
+		0x6d1a7918c80df518, 0xdf9a4939342308e9, 0xebc2151ee6c8398c, 0x03cc2ba8a1116515,
+		0xd341d037e840cf83, 0x387cb5d25af4afcc, 0xbba2515f22909e87, 0x7248fe7705f38e47,
+		0x4d61e56a525d225a, 0x262e963c8da05d3d, 0x59e89b094d220ec2, 0x055d5b52b78b9c5e,
+		0x82b27eb33514ef99, 0xd30094ca96b7ce7b, 0xcf5cb381cd0a1535, 0xfeed4db6919e5a7c,
+		0x41703f53753be59f, 0x5eeea940fcde8b6f, 0x4cd1f1b175100206, 0x4a20358574454ec0,
+		0x1478d361dbbf9fac, 0x6f02dc07d141875c, 0x296a202ed8e556a2, 0x2afd67999bf32ee5,
+		0x7acfd96efa95491d, 0x6798ba0c0abb2c6d, 0x34c6f57b26c92122, 0x5736e1bad206b5de,
+		0x20057d2a0056521b, 0x3dea5bd5d0578bd7, 0x16e50d897d4634ac, 0x29bff3ecb9b7a6e3,
+		0x475cd3205a3bdcde, 0x18a42105c31b7e88, 0x023e7414af663068, 0x15147108121967d7,
+		0xe4a3dff1d7d6fef9, 0x01a8d1a588085737, 0x11b4c74eda62beef, 0xe587cc0d69a73346,
+		0x1ff7327017aa2a6e, 0x594e29c42473d06b, 0xf6f31db1899b12d5, 0xc02ac5e47312d3ca,
+		0xe70201e960cb78b8, 0x6f90ff3b6a65f108, 0x42747a7245e7fa84, 0xd1f507e43ab749b2,
+		0x1c86d265f15750cd, 0x3996ce73dd832c1c, 0x8e7fba02983224bd, 0xba0dec7103255dd4,
+		0x9e9cbd781628fc5b, 0xdae8645996edd6a5, 0xdebe0853b1a1d378, 0xa49229d24d014343,
+		0x7be5b9ffda905e1c, 0xa3c95eaec244aa30, 0x0230bca8f4df0544, 0x4135c2bebfe148c6,
+		0x166fc0cc438a3c72, 0x3762b59a8ae83efa, 0xe8928a4c89114750, 0x2a440b51a4945ee5,
+		0x80cefd2b7d99ff83, 0xbb9879c6e61fd62a, 0x6e7c8f1a84265034, 0x164bb2de1bbeddc8,
+		0xf3c12fe54d5c653b, 0x40b9e922ed9771e2, 0x551f5b0fbe7b1840, 0x25032aa7c4cb1811,
+		0xaaed34074b164346, 0x8ffd96bbf9c9c81d, 0x70fc91eb5937085c, 0x7f795e2a5f915440,
+		0x4543d9df5476d3cb, 0xf172d73e004fc90d, 0xdfd1c4febcc81238, 0xbc8dfb627fe558fc,
+	}
+
+	// C is a constant array of element
+	C []*ffg.Element
+	// M is a matrix
+	M [][]*ffg.Element
+)
+
+func init() {
+	for i := 0; i < len(c); i++ {
+		C = append(C, ffg.NewElementFromUint64(c[i]))
+	}
+
+	for i := 0; i < mLen; i++ {
+		var row []*ffg.Element
+		for j := 0; j < mLen; j++ {
+			ele := ffg.NewElementFromUint64(mcirc[(-i+j+mLen)%mLen])
+			if i == j {
+				ele = ffg.NewElementFromUint64(mcirc[0] + mdiag[i])
+			}
+			row = append(row, ele)
+		}
+		M = append(M, row)
+	}
+}
--- a/goldenposeidon/poseidon.go
+++ b/goldenposeidon/poseidon.go
@@ -0,0 +1,77 @@
+package poseidon
+
+import (
+	"math/big"
+
+	"github.com/iden3/go-iden3-crypto/ffg"
+)
+
+func zero() *ffg.Element {
+	return ffg.NewElement()
+}
+
+// exp7 performs x^7 mod p
+func exp7(a *ffg.Element) {
+	a.Exp(*a, big.NewInt(7)) //nolint:gomnd
+}
+
+// exp7state perform exp7 for whole state
+func exp7state(state []*ffg.Element) {
+	for i := 0; i < len(state); i++ {
+		exp7(state[i])
+	}
+}
+
+// ark computes Add-Round Key, from the paper https://eprint.iacr.org/2019/458.pdf
+func ark(state []*ffg.Element, it int) {
+	for i := 0; i < len(state); i++ {
+		state[i].Add(state[i], C[it+i])
+	}
+}
+
+// mix returns [[matrix]] * [vector]
+func mix(state []*ffg.Element) []*ffg.Element {
+	mul := zero()
+	newState := make([]*ffg.Element, mLen)
+	for i := 0; i < mLen; i++ {
+		newState[i] = zero()
+	}
+	for i := 0; i < mLen; i++ {
+		newState[i].SetUint64(0)
+		for j := 0; j < mLen; j++ {
+			mul.Mul(M[i][j], state[j])
+			newState[i].Add(newState[i], mul)
+		}
+	}
+	return newState
+}
+
+// Hash computes the Poseidon hash for the given inputs
+func Hash(inpBI [NROUNDSF]uint64, capBI [CAPLEN]uint64) ([CAPLEN]uint64, error) {
+	state := make([]*ffg.Element, mLen)
+	for i := 0; i < NROUNDSF; i++ {
+		state[i] = ffg.NewElement().SetUint64(inpBI[i])
+	}
+	for i := 0; i < CAPLEN; i++ {
+		state[i+NROUNDSF] = ffg.NewElement().SetUint64(capBI[i])
+	}
+
+	for r := 0; r < NROUNDSF+NROUNDSP; r++ {
+		ark(state, r*mLen)
+
+		if r < NROUNDSF/2 || r >= NROUNDSF/2+NROUNDSP {
+			exp7state(state)
+		} else {
+			exp7(state[0])
+		}
+
+		state = mix(state)
+	}
+
+	return [CAPLEN]uint64{
+		state[0].ToUint64Regular(),
+		state[1].ToUint64Regular(),
+		state[2].ToUint64Regular(),
+		state[3].ToUint64Regular(),
+	}, nil
+}
--- a/goldenposeidon/poseidon_test.go
+++ b/goldenposeidon/poseidon_test.go
@@ -0,0 +1,94 @@
+package poseidon
+
+import (
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+const prime uint64 = 18446744069414584321
+
+func TestPoseidonHash(t *testing.T) {
+	b0 := uint64(0)
+	b1 := uint64(1)
+	bm1 := prime - 1
+	bM := prime
+
+	h, err := Hash([NROUNDSF]uint64{b0, b0, b0, b0, b0, b0, b0, b0}, [CAPLEN]uint64{b0, b0, b0, b0})
+	assert.Nil(t, err)
+	assert.Equal(t,
+		[CAPLEN]uint64{
+			4330397376401421145,
+			14124799381142128323,
+			8742572140681234676,
+			14345658006221440202,
+		}, h,
+	)
+
+	h, err = Hash([NROUNDSF]uint64{b1, b1, b1, b1, b1, b1, b1, b1}, [CAPLEN]uint64{b1, b1, b1, b1})
+	assert.Nil(t, err)
+	assert.Equal(t,
+		[CAPLEN]uint64{
+			16428316519797902711,
+			13351830238340666928,
+			682362844289978626,
+			12150588177266359240,
+		}, h,
+	)
+
+	h, err = Hash([NROUNDSF]uint64{b1, b1, b1, b1, b1, b1, b1, b1}, [CAPLEN]uint64{b1, b1, b1, b1})
+	assert.Nil(t, err)
+	assert.Equal(t,
+		[CAPLEN]uint64{
+			16428316519797902711,
+			13351830238340666928,
+			682362844289978626,
+			12150588177266359240,
+		}, h,
+	)
+
+	h, err = Hash(
+		[NROUNDSF]uint64{bm1, bm1, bm1, bm1, bm1, bm1, bm1, bm1},
+		[CAPLEN]uint64{bm1, bm1, bm1, bm1},
+	)
+	assert.Nil(t, err)
+	assert.Equal(t,
+		[CAPLEN]uint64{
+			13691089994624172887,
+			15662102337790434313,
+			14940024623104903507,
+			10772674582659927682,
+		}, h,
+	)
+
+	h, err = Hash([NROUNDSF]uint64{bM, bM, bM, bM, bM, bM, bM, bM}, [CAPLEN]uint64{b0, b0, b0, b0})
+	assert.Nil(t, err)
+	assert.Equal(t,
+		[CAPLEN]uint64{
+			4330397376401421145,
+			14124799381142128323,
+			8742572140681234676,
+			14345658006221440202,
+		}, h,
+	)
+
+	h, err = Hash([NROUNDSF]uint64{
+		uint64(923978),
+		uint64(235763497586),
+		uint64(9827635653498),
+		uint64(112870),
+		uint64(289273673480943876),
+		uint64(230295874986745876),
+		uint64(6254867324987),
+		uint64(2087),
+	}, [CAPLEN]uint64{b0, b0, b0, b0})
+	assert.Nil(t, err)
+	assert.Equal(t,
+		[CAPLEN]uint64{
+			1892171027578617759,
+			984732815927439256,
+			7866041765487844082,
+			8161503938059336191,
+		}, h,
+	)
+}