Use optimized goldilocks in codebase (#26)

* gl * stage 1 optimizations * working optimized poseidon * Fix posedion tests * in progress gate type refactor * working gates * working e2e * hm' * hm2 * debug saga continues * more debugging cry * more debug * it finally works * optimizations * more optimizations * new changes * more optimizations * more cleanup * some refactoring * new files * flattening of packages * working commit * more refactor * more flattening * more flattening * more more refactor * more optimizations * more optimizations * more optimizations * plonk benchmark * plonk * fix r1cs * resolve kevin's comments * Update goldilocks/base.go Co-authored-by: Kevin Jue <kjue235@gmail.com> * Update goldilocks/base.go Co-authored-by: Kevin Jue <kjue235@gmail.com> * Update goldilocks/base.go Co-authored-by: Kevin Jue <kjue235@gmail.com> * Update goldilocks/quadratic_extension.go Co-authored-by: Kevin Jue <kjue235@gmail.com> * fix: resolve kevin's confusion --------- Co-authored-by: Kevin Jue <kjue235@gmail.com>
2026-01-12 09:01:32 +01:00 · 2023-07-24 16:08:17 -07:00
parent 103c7ca47d
commit b670530e7f
61 changed files with 3506 additions and 3211 deletions
--- a/goldilocks/base.go
+++ b/goldilocks/base.go
@@ -0,0 +1,362 @@
+// This package implements efficient Golidlocks arithmetic operations within Gnark. We do not use
+// the emulated field arithmetic API, because it is too slow for our purposes. Instead, we use
+// an efficient reduction method that leverages the fact that the modulus is a simple
+// linear combination of powers of two.
+package goldilocks
+
+// In general, methods whose name do not contain `NoReduce` can be used without any extra mental
+// overhead. These methods act exactly as you would expect a normal field would operate.
+//
+// However, if you want to aggressively optimize the number of constraints in your circuit, it can
+// be very beneficial to use the no reduction methods and keep track of the maximum number of bits
+// your computation uses.
+
+import (
+	"fmt"
+	"math/big"
+
+	"github.com/consensys/gnark-crypto/field/goldilocks"
+	"github.com/consensys/gnark/constraint/solver"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/std/math/bits"
+	"github.com/consensys/gnark/std/math/emulated"
+)
+
+// The multiplicative group generator of the field.
+var MULTIPLICATIVE_GROUP_GENERATOR goldilocks.Element = goldilocks.NewElement(7)
+
+// The two adicity of the field.
+var TWO_ADICITY uint64 = 32
+
+// The power of two generator of the field.
+var POWER_OF_TWO_GENERATOR goldilocks.Element = goldilocks.NewElement(1753635133440165772)
+
+// The modulus of the field.
+var MODULUS *big.Int = emulated.Goldilocks{}.Modulus()
+
+// The threshold maximum number of bits at which we must reduce the element.
+var REDUCE_NB_BITS_THRESHOLD uint8 = 254 - 64
+
+// The number of bits to use for range checks on inner products of field elements.
+var RANGE_CHECK_NB_BITS int = 140
+
+// Registers the hint functions with the solver.
+func init() {
+	solver.RegisterHint(MulAddHint)
+	solver.RegisterHint(ReduceHint)
+	solver.RegisterHint(InverseHint)
+}
+
+// A type alias used to represent Goldilocks field elements.
+type Variable struct {
+	Limb frontend.Variable
+}
+
+// Creates a new Goldilocks field element from an existing variable. Assumes that the element is
+// already reduced.
+func NewVariable(x frontend.Variable) Variable {
+	return Variable{Limb: x}
+}
+
+// The zero element in the Golidlocks field.
+func Zero() Variable {
+	return NewVariable(0)
+}
+
+// The one element in the Goldilocks field.
+func One() Variable {
+	return NewVariable(1)
+}
+
+// The negative one element in the Goldilocks field.
+func NegOne() Variable {
+	return NewVariable(MODULUS.Uint64() - 1)
+}
+
+// The chip used for Goldilocks field operations.
+type Chip struct {
+	api frontend.API
+}
+
+// Creates a new Goldilocks chip.
+func NewChip(api frontend.API) *Chip {
+	return &Chip{api: api}
+}
+
+// Adds two field elements such that x + y = z within the Golidlocks field.
+func (p *Chip) Add(a Variable, b Variable) Variable {
+	return p.MulAdd(a, NewVariable(1), b)
+}
+
+// Adds two field elements such that x + y = z within the Golidlocks field without reducing.
+func (p *Chip) AddNoReduce(a Variable, b Variable) Variable {
+	return NewVariable(p.api.Add(a.Limb, b.Limb))
+}
+
+// Subtracts two field elements such that x + y = z within the Golidlocks field.
+func (p *Chip) Sub(a Variable, b Variable) Variable {
+	return p.MulAdd(b, NewVariable(MODULUS.Uint64()-1), a)
+}
+
+// Subtracts two field elements such that x + y = z within the Golidlocks field without reducing.
+func (p *Chip) SubNoReduce(a Variable, b Variable) Variable {
+	return NewVariable(p.api.Add(a.Limb, p.api.Mul(b.Limb, MODULUS.Uint64()-1)))
+}
+
+// Multiplies two field elements such that x * y = z within the Golidlocks field.
+func (p *Chip) Mul(a Variable, b Variable) Variable {
+	return p.MulAdd(a, b, Zero())
+}
+
+// Multiplies two field elements such that x * y = z within the Golidlocks field without reducing.
+func (p *Chip) MulNoReduce(a Variable, b Variable) Variable {
+	return NewVariable(p.api.Mul(a.Limb, b.Limb))
+}
+
+// Multiplies two field elements and adds a field element such that x * y + z = c within the
+// Golidlocks field.
+func (p *Chip) MulAdd(a Variable, b Variable, c Variable) Variable {
+	result, err := p.api.Compiler().NewHint(MulAddHint, 2, a.Limb, b.Limb, c.Limb)
+	if err != nil {
+		panic(err)
+	}
+
+	quotient := NewVariable(result[0])
+	remainder := NewVariable(result[1])
+
+	lhs := p.api.Mul(a.Limb, b.Limb)
+	lhs = p.api.Add(lhs, c.Limb)
+	rhs := p.api.Add(p.api.Mul(quotient.Limb, MODULUS), remainder.Limb)
+	p.api.AssertIsEqual(lhs, rhs)
+
+	p.RangeCheck(quotient)
+	p.RangeCheck(remainder)
+	return remainder
+}
+
+// Multiplies two field elements and adds a field element such that x * y + z = c within the
+// Golidlocks field without reducing.
+func (p *Chip) MulAddNoReduce(a Variable, b Variable, c Variable) Variable {
+	return p.AddNoReduce(p.MulNoReduce(a, b), c)
+}
+
+// The hint used to compute MulAdd.
+func MulAddHint(_ *big.Int, inputs []*big.Int, results []*big.Int) error {
+	if len(inputs) != 3 {
+		panic("MulAddHint expects 3 input operands")
+	}
+
+	for _, operand := range inputs {
+		if operand.Cmp(MODULUS) >= 0 {
+			panic(fmt.Sprintf("%s is not in the field", operand.String()))
+		}
+	}
+
+	product := new(big.Int).Mul(inputs[0], inputs[1])
+	sum := new(big.Int).Add(product, inputs[2])
+	quotient := new(big.Int).Div(sum, MODULUS)
+	remainder := new(big.Int).Rem(sum, MODULUS)
+
+	results[0] = quotient
+	results[1] = remainder
+
+	return nil
+}
+
+// Reduces a field element x such that x % MODULUS = y.
+func (p *Chip) Reduce(x Variable) Variable {
+	// Witness a `quotient` and `remainder` such that:
+	//
+	// 		MODULUS * quotient + remainder = x
+	//
+	// Must check that offset \in [0, MODULUS) and carry \in [0, 2^RANGE_CHECK_NB_BITS) to ensure
+	// that this computation does not overflow. We use 2^RANGE_CHECK_NB_BITS to reduce the cost of the range check
+	//
+	// In other words, we assume that we at most compute a a dot product with dimension at most RANGE_CHECK_NB_BITS - 128.
+
+	result, err := p.api.Compiler().NewHint(ReduceHint, 2, x.Limb)
+	if err != nil {
+		panic(err)
+	}
+
+	quotient := result[0]
+	rangeCheckNbBits := RANGE_CHECK_NB_BITS
+	p.api.ToBinary(quotient, rangeCheckNbBits)
+
+	remainder := NewVariable(result[1])
+	p.RangeCheck(remainder)
+	return remainder
+}
+
+// Reduces a field element x such that x % MODULUS = y.
+func (p *Chip) ReduceWithMaxBits(x Variable, maxNbBits uint64) Variable {
+	// Witness a `quotient` and `remainder` such that:
+	//
+	// 		MODULUS * quotient + remainder = x
+	//
+	// Must check that remainder \in [0, MODULUS) and quotient \in [0, 2^maxNbBits) to ensure that this
+	// computation does not overflow.
+
+	result, err := p.api.Compiler().NewHint(ReduceHint, 2, x.Limb)
+	if err != nil {
+		panic(err)
+	}
+
+	quotient := result[0]
+	p.api.ToBinary(quotient, int(maxNbBits))
+
+	remainder := NewVariable(result[1])
+	p.RangeCheck(remainder)
+	return remainder
+}
+
+// The hint used to compute Reduce.
+func ReduceHint(_ *big.Int, inputs []*big.Int, results []*big.Int) error {
+	if len(inputs) != 1 {
+		panic("ReduceHint expects 1 input operand")
+	}
+	input := inputs[0]
+	quotient := new(big.Int).Div(input, MODULUS)
+	remainder := new(big.Int).Rem(input, MODULUS)
+	results[0] = quotient
+	results[1] = remainder
+	return nil
+}
+
+// Computes the inverse of a field element x such that x * x^-1 = 1.
+func (p *Chip) Inverse(x Variable) Variable {
+	result, err := p.api.Compiler().NewHint(InverseHint, 1, x.Limb)
+	if err != nil {
+		panic(err)
+	}
+
+	inverse := NewVariable(result[0])
+	product := p.Mul(inverse, x)
+	p.api.AssertIsEqual(product.Limb, frontend.Variable(1))
+	return inverse
+}
+
+// The hint used to compute Inverse.
+func InverseHint(_ *big.Int, inputs []*big.Int, results []*big.Int) error {
+	if len(inputs) != 1 {
+		panic("InverseHint expects 1 input operand")
+	}
+
+	input := inputs[0]
+	if input.Cmp(MODULUS) == 0 || input.Cmp(MODULUS) == 1 {
+		panic("Input is not in the field")
+	}
+
+	inputGl := goldilocks.NewElement(input.Uint64())
+	resultGl := goldilocks.NewElement(0)
+	resultGl.Inverse(&inputGl)
+
+	result := big.NewInt(0)
+	results[0] = resultGl.BigInt(result)
+
+	return nil
+}
+
+// Computes a field element raised to some power.
+func (p *Chip) Exp(x Variable, k *big.Int) Variable {
+	if k.IsUint64() && k.Uint64() == 0 {
+		return One()
+	}
+
+	e := k
+	if k.Sign() == -1 {
+		panic("Unsupported negative exponent. Need to implement inversion.")
+	}
+
+	z := x
+	for i := e.BitLen() - 2; i >= 0; i-- {
+		z = p.Mul(z, z)
+		if e.Bit(i) == 1 {
+			z = p.Mul(z, x)
+		}
+	}
+
+	return z
+}
+
+// Range checks a field element x to be less than the Golidlocks modulus 2 ^ 64 - 2 ^ 32 + 1.
+func (p *Chip) RangeCheck(x Variable) {
+	// The Goldilocks' modulus is 2^64 - 2^32 + 1, which is:
+	//
+	// 		1111111111111111111111111111111100000000000000000000000000000001
+	//
+	// in big endian binary. This function will first verify that x is at most 64 bits wide. Then it
+	// checks that if the bits[0:31] (in big-endian) are all 1, then bits[32:64] are all zero.
+
+	// First decompose x into 64 bits.  The bits will be in little-endian order.
+	bits := bits.ToBinary(p.api, x.Limb, bits.WithNbDigits(64))
+
+	// Those bits should compose back to x.
+	reconstructedX := frontend.Variable(0)
+	c := uint64(1)
+	for i := 0; i < 64; i++ {
+		reconstructedX = p.api.Add(reconstructedX, p.api.Mul(bits[i], c))
+		c = c << 1
+		p.api.AssertIsBoolean(bits[i])
+	}
+	p.api.AssertIsEqual(x.Limb, reconstructedX)
+
+	mostSigBits32Sum := frontend.Variable(0)
+	for i := 32; i < 64; i++ {
+		mostSigBits32Sum = p.api.Add(mostSigBits32Sum, bits[i])
+	}
+
+	leastSigBits32Sum := frontend.Variable(0)
+	for i := 0; i < 32; i++ {
+		leastSigBits32Sum = p.api.Add(leastSigBits32Sum, bits[i])
+	}
+
+	// If mostSigBits32Sum < 32, then we know that:
+	//
+	// 		x < (2^63 + ... + 2^32 + 0 * 2^31 + ... + 0 * 2^0)
+	//
+	// which equals to 2^64 - 2^32. So in that case, we don't need to do any more checks. If
+	// mostSigBits32Sum == 32, then we need to check that x == 2^64 - 2^32 (max GL value).
+	shouldCheck := p.api.IsZero(p.api.Sub(mostSigBits32Sum, 32))
+	p.api.AssertIsEqual(
+		p.api.Select(
+			shouldCheck,
+			leastSigBits32Sum,
+			frontend.Variable(0),
+		),
+		frontend.Variable(0),
+	)
+}
+
+func (p *Chip) AssertIsEqual(x, y Variable) {
+	p.api.AssertIsEqual(x.Limb, y.Limb)
+}
+
+// Computes the n'th primitive root of unity for the Goldilocks field.
+func PrimitiveRootOfUnity(nLog uint64) goldilocks.Element {
+	if nLog > TWO_ADICITY {
+		panic("nLog is greater than TWO_ADICITY")
+	}
+	res := goldilocks.NewElement(POWER_OF_TWO_GENERATOR.Uint64())
+	for i := 0; i < int(TWO_ADICITY-nLog); i++ {
+		res.Square(&res)
+	}
+	return res
+}
+
+func TwoAdicSubgroup(nLog uint64) []goldilocks.Element {
+	if nLog > TWO_ADICITY {
+		panic("nLog is greater than GOLDILOCKS_TWO_ADICITY")
+	}
+
+	var res []goldilocks.Element
+	rootOfUnity := PrimitiveRootOfUnity(nLog)
+	res = append(res, goldilocks.NewElement(1))
+
+	for i := 0; i < (1 << nLog); i++ {
+		lastElement := res[len(res)-1]
+		res = append(res, *lastElement.Mul(&lastElement, &rootOfUnity))
+	}
+
+	return res
+}
--- a/goldilocks/base_test.go
+++ b/goldilocks/base_test.go
@@ -0,0 +1,73 @@
+package goldilocks
+
+import (
+	"math/big"
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	"github.com/consensys/gnark/backend"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/test"
+)
+
+type TestGoldilocksRangeCheckCircuit struct {
+	X frontend.Variable
+}
+
+func (c *TestGoldilocksRangeCheckCircuit) Define(api frontend.API) error {
+	chip := NewChip(api)
+	chip.RangeCheck(NewVariable(c.X))
+	return nil
+}
+func TestGoldilocksRangeCheck(t *testing.T) {
+	assert := test.NewAssert(t)
+
+	var circuit, witness TestGoldilocksRangeCheckCircuit
+
+	witness.X = 1
+	assert.ProverSucceeded(&circuit, &witness, test.WithCurves(ecc.BN254), test.WithBackends(backend.GROTH16), test.NoSerialization())
+
+	witness.X = 0
+	assert.ProverSucceeded(&circuit, &witness, test.WithCurves(ecc.BN254), test.WithBackends(backend.GROTH16), test.NoSerialization())
+
+	witness.X = MODULUS
+	assert.ProverFailed(&circuit, &witness, test.WithCurves(ecc.BN254), test.WithBackends(backend.GROTH16), test.NoSerialization())
+
+	one := big.NewInt(1)
+	maxValidVal := new(big.Int).Sub(MODULUS, one)
+	witness.X = maxValidVal
+	assert.ProverSucceeded(&circuit, &witness, test.WithCurves(ecc.BN254), test.WithBackends(backend.GROTH16))
+}
+
+type TestGoldilocksMulAddCircuit struct {
+	X, Y, Z        frontend.Variable
+	ExpectedResult frontend.Variable
+}
+
+func (c *TestGoldilocksMulAddCircuit) Define(api frontend.API) error {
+	chip := NewChip(api)
+	calculateValue := chip.MulAdd(NewVariable(c.X), NewVariable(c.Y), NewVariable(c.Z))
+	api.AssertIsEqual(calculateValue.Limb, c.ExpectedResult)
+	return nil
+}
+
+func TestGoldilocksMulAdd(t *testing.T) {
+	assert := test.NewAssert(t)
+
+	var circuit, witness TestGoldilocksMulAddCircuit
+
+	witness.X = 1
+	witness.Y = 2
+	witness.Z = 3
+	witness.ExpectedResult = 5
+	assert.ProverSucceeded(&circuit, &witness, test.WithCurves(ecc.BN254), test.WithBackends(backend.GROTH16), test.NoFuzzing())
+
+	bigOperand := new(big.Int).SetUint64(9223372036854775808)
+	expectedValue, _ := new(big.Int).SetString("18446744068340842500", 10)
+
+	witness.X = bigOperand
+	witness.Y = bigOperand
+	witness.Z = 3
+	witness.ExpectedResult = expectedValue
+	assert.ProverSucceeded(&circuit, &witness, test.WithCurves(ecc.BN254), test.WithBackends(backend.GROTH16), test.NoFuzzing())
+}
--- a/goldilocks/quadratic_extension.go
+++ b/goldilocks/quadratic_extension.go
@@ -0,0 +1,234 @@
+package goldilocks
+
+import (
+	"math/bits"
+
+	"github.com/consensys/gnark/frontend"
+)
+
+const W uint64 = 7
+const DTH_ROOT uint64 = 18446744069414584320
+
+type QuadraticExtensionVariable [2]Variable
+
+func NewQuadraticExtensionVariable(x Variable, y Variable) QuadraticExtensionVariable {
+	return QuadraticExtensionVariable{x, y}
+}
+
+func (p Variable) ToQuadraticExtension() QuadraticExtensionVariable {
+	return NewQuadraticExtensionVariable(p, Zero())
+}
+
+func ZeroExtension() QuadraticExtensionVariable {
+	return Zero().ToQuadraticExtension()
+}
+
+func OneExtension() QuadraticExtensionVariable {
+	return One().ToQuadraticExtension()
+}
+
+// Adds two quadratic extension variables in the Goldilocks field.
+func (p *Chip) AddExtension(a, b QuadraticExtensionVariable) QuadraticExtensionVariable {
+	c0 := p.Add(a[0], b[0])
+	c1 := p.Add(a[1], b[1])
+	return NewQuadraticExtensionVariable(c0, c1)
+}
+
+// Adds two quadratic extension variables in the Goldilocks field without reducing.
+func (p *Chip) AddExtensionNoReduce(a, b QuadraticExtensionVariable) QuadraticExtensionVariable {
+	c0 := p.AddNoReduce(a[0], b[0])
+	c1 := p.AddNoReduce(a[1], b[1])
+	return NewQuadraticExtensionVariable(c0, c1)
+}
+
+// Subtracts two quadratic extension variables in the Goldilocks field.
+func (p *Chip) SubExtension(a, b QuadraticExtensionVariable) QuadraticExtensionVariable {
+	c0 := p.Sub(a[0], b[0])
+	c1 := p.Sub(a[1], b[1])
+	return NewQuadraticExtensionVariable(c0, c1)
+}
+
+// Subtracts two quadratic extension variables in the Goldilocks field without reducing.
+func (p *Chip) SubExtensionNoReduce(a, b QuadraticExtensionVariable) QuadraticExtensionVariable {
+	c0 := p.SubNoReduce(a[0], b[0])
+	c1 := p.SubNoReduce(a[1], b[1])
+	return NewQuadraticExtensionVariable(c0, c1)
+}
+
+// Multiplies quadratic extension variable in the Goldilocks field.
+func (p *Chip) MulExtension(a, b QuadraticExtensionVariable) QuadraticExtensionVariable {
+	product := p.MulExtensionNoReduce(a, b)
+	product[0] = p.Reduce(product[0])
+	product[1] = p.Reduce(product[1])
+	return product
+}
+
+// Multiplies quadratic extension variable in the Goldilocks field without reducing.
+func (p *Chip) MulExtensionNoReduce(a, b QuadraticExtensionVariable) QuadraticExtensionVariable {
+	c0o0 := p.MulNoReduce(a[0], b[0])
+	c0o1 := p.MulNoReduce(p.MulNoReduce(NewVariable(7), a[1]), b[1])
+	c0 := p.AddNoReduce(c0o0, c0o1)
+	c1 := p.AddNoReduce(p.MulNoReduce(a[0], b[1]), p.MulNoReduce(a[1], b[0]))
+	return NewQuadraticExtensionVariable(c0, c1)
+}
+
+// Multiplies two operands a and b and adds to c in the Goldilocks extension field. a * b + c must
+// be less than RANGE_CHECK_NB_BITS bits.
+func (p *Chip) MulAddExtension(a, b, c QuadraticExtensionVariable) QuadraticExtensionVariable {
+	product := p.MulExtensionNoReduce(a, b)
+	sum := p.AddExtensionNoReduce(product, c)
+	sum[0] = p.Reduce(sum[0])
+	sum[1] = p.Reduce(sum[1])
+	return sum
+}
+
+func (p *Chip) MulAddExtensionNoReduce(a, b, c QuadraticExtensionVariable) QuadraticExtensionVariable {
+	product := p.MulExtensionNoReduce(a, b)
+	sum := p.AddExtensionNoReduce(product, c)
+	return sum
+}
+
+// Multiplies two operands a and b and subtracts to c in the Goldilocks extension field. a * b - c must
+// be less than RANGE_CHECK_NB_BITS bits.
+func (p *Chip) SubMulExtension(a, b, c QuadraticExtensionVariable) QuadraticExtensionVariable {
+	difference := p.SubExtensionNoReduce(a, b)
+	product := p.MulExtensionNoReduce(difference, c)
+	product[0] = p.Reduce(product[0])
+	product[1] = p.Reduce(product[1])
+	return product
+}
+
+// Multiplies quadratic extension variable in the Goldilocks field by a scalar.
+func (p *Chip) ScalarMulExtension(
+	a QuadraticExtensionVariable,
+	b Variable,
+) QuadraticExtensionVariable {
+	return NewQuadraticExtensionVariable(
+		p.Mul(a[0], b),
+		p.Mul(a[1], b),
+	)
+}
+
+// Computes an inner product over quadratic extension variable vectors in the Goldilocks field.
+func (p *Chip) InnerProductExtension(
+	constant Variable,
+	startingAcc QuadraticExtensionVariable,
+	pairs [][2]QuadraticExtensionVariable,
+) QuadraticExtensionVariable {
+	acc := startingAcc
+	for i := 0; i < len(pairs); i++ {
+		a := pairs[i][0]
+		b := pairs[i][1]
+		mul := p.ScalarMulExtension(a, constant)
+		acc = p.MulAddExtensionNoReduce(mul, b, acc)
+	}
+	return p.ReduceExtension(acc)
+}
+
+// Computes the inverse of a quadratic extension variable in the Goldilocks field.
+func (p *Chip) InverseExtension(a QuadraticExtensionVariable) QuadraticExtensionVariable {
+	a0IsZero := p.api.IsZero(a[0].Limb)
+	a1IsZero := p.api.IsZero(a[1].Limb)
+	p.api.AssertIsEqual(p.api.Mul(a0IsZero, a1IsZero), frontend.Variable(0))
+	aPowRMinus1 := QuadraticExtensionVariable{
+		a[0],
+		p.Mul(a[1], NewVariable(DTH_ROOT)),
+	}
+	aPowR := p.MulExtension(aPowRMinus1, a)
+	return p.ScalarMulExtension(aPowRMinus1, p.Inverse(aPowR[0]))
+}
+
+// Divides two quadratic extension variables in the Goldilocks field.
+func (p *Chip) DivExtension(a, b QuadraticExtensionVariable) QuadraticExtensionVariable {
+	return p.MulExtension(a, p.InverseExtension(b))
+}
+
+// Exponentiates a quadratic extension variable to some exponent in the Golidlocks field.
+func (p *Chip) ExpExtension(
+	a QuadraticExtensionVariable,
+	exponent uint64,
+) QuadraticExtensionVariable {
+	switch exponent {
+	case 0:
+		return OneExtension()
+	case 1:
+		return a
+	case 2:
+		return p.MulExtension(a, a)
+	default:
+	}
+
+	current := a
+	product := OneExtension()
+
+	for i := 0; i < bits.Len64(exponent); i++ {
+		if i != 0 {
+			current = p.MulExtension(current, current)
+		}
+		if (exponent >> i & 1) != 0 {
+			product = p.MulExtension(product, current)
+		}
+	}
+
+	return product
+}
+
+func (p *Chip) ReduceExtension(x QuadraticExtensionVariable) QuadraticExtensionVariable {
+	return NewQuadraticExtensionVariable(p.Reduce(x[0]), p.Reduce(x[1]))
+}
+
+// Reduces a list of extension field terms with a scalar power in the Goldilocks field.
+func (p *Chip) ReduceWithPowers(
+	terms []QuadraticExtensionVariable,
+	scalar QuadraticExtensionVariable,
+) QuadraticExtensionVariable {
+	sum := ZeroExtension()
+	for i := len(terms) - 1; i >= 0; i-- {
+		sum = p.AddExtensionNoReduce(
+			p.MulExtensionNoReduce(
+				sum,
+				scalar,
+			),
+			terms[i],
+		)
+		sum = p.ReduceExtension(sum)
+	}
+	return sum
+}
+
+// Outputs whether the quadratic extension variable is zero.
+func (p *Chip) IsZero(x QuadraticExtensionVariable) frontend.Variable {
+	x0IsZero := p.api.IsZero(x[0].Limb)
+	x1IsZero := p.api.IsZero(x[1].Limb)
+	return p.api.Mul(x0IsZero, x1IsZero)
+}
+
+// Lookup is similar to select, but returns the first variable if the bit is zero and vice-versa.
+func (p *Chip) Lookup(
+	b frontend.Variable,
+	x, y QuadraticExtensionVariable,
+) QuadraticExtensionVariable {
+	c0 := p.api.Select(b, y[0].Limb, x[0].Limb)
+	c1 := p.api.Select(b, y[1].Limb, x[1].Limb)
+	return NewQuadraticExtensionVariable(NewVariable(c0), NewVariable(c1))
+}
+
+// Lookup2 is similar to select2, but returns the first variable if the bit is zero and vice-versa.
+func (p *Chip) Lookup2(
+	b0 frontend.Variable,
+	b1 frontend.Variable,
+	qe0, qe1, qe2, qe3 QuadraticExtensionVariable,
+) QuadraticExtensionVariable {
+	c0 := p.Lookup(b0, qe0, qe1)
+	c1 := p.Lookup(b0, qe2, qe3)
+	return p.Lookup(b1, c0, c1)
+}
+
+// Asserts that two quadratic extension variables are equal.
+func (p *Chip) AssertIsEqualExtension(
+	a QuadraticExtensionVariable,
+	b QuadraticExtensionVariable,
+) {
+	p.AssertIsEqual(a[0], b[0])
+	p.AssertIsEqual(a[1], b[1])
+}
--- a/goldilocks/quadratic_extension_algebra.go
+++ b/goldilocks/quadratic_extension_algebra.go
@@ -0,0 +1,125 @@
+package goldilocks
+
+import "github.com/consensys/gnark-crypto/field/goldilocks"
+
+const D = 2
+
+type QuadraticExtensionAlgebraVariable = [D]QuadraticExtensionVariable
+
+func NewQuadraticExtensionAlgebraVariable(
+	a QuadraticExtensionVariable,
+	b QuadraticExtensionVariable,
+) QuadraticExtensionAlgebraVariable {
+	return QuadraticExtensionAlgebraVariable{a, b}
+}
+
+func (p QuadraticExtensionVariable) ToQuadraticExtensionAlgebra() QuadraticExtensionAlgebraVariable {
+	return [2]QuadraticExtensionVariable{p, ZeroExtension()}
+}
+
+func ZeroExtensionAlgebra() QuadraticExtensionAlgebraVariable {
+	return ZeroExtension().ToQuadraticExtensionAlgebra()
+}
+
+func OneExtensionAlgebra() QuadraticExtensionAlgebraVariable {
+	return OneExtension().ToQuadraticExtensionAlgebra()
+}
+
+func (p *Chip) AddExtensionAlgebra(
+	a QuadraticExtensionAlgebraVariable,
+	b QuadraticExtensionAlgebraVariable,
+) QuadraticExtensionAlgebraVariable {
+	var sum QuadraticExtensionAlgebraVariable
+	for i := 0; i < D; i++ {
+		sum[i] = p.AddExtension(a[i], b[i])
+	}
+	return sum
+}
+
+func (p *Chip) SubExtensionAlgebra(
+	a QuadraticExtensionAlgebraVariable,
+	b QuadraticExtensionAlgebraVariable,
+) QuadraticExtensionAlgebraVariable {
+	var diff QuadraticExtensionAlgebraVariable
+	for i := 0; i < D; i++ {
+		diff[i] = p.SubExtension(a[i], b[i])
+	}
+	return diff
+}
+
+func (p Chip) MulExtensionAlgebra(
+	a QuadraticExtensionAlgebraVariable,
+	b QuadraticExtensionAlgebraVariable,
+) QuadraticExtensionAlgebraVariable {
+	var inner [D][]QuadraticExtensionAlgebraVariable
+	var innerW [D][]QuadraticExtensionAlgebraVariable
+
+	for i := 0; i < D; i++ {
+		for j := 0; j < D-i; j++ {
+			idx := (i + j) % D
+			inner[idx] = append(inner[idx], QuadraticExtensionAlgebraVariable{a[i], b[j]})
+		}
+		for j := D - i; j < D; j++ {
+			idx := (i + j) % D
+			innerW[idx] = append(innerW[idx], QuadraticExtensionAlgebraVariable{a[i], b[j]})
+		}
+	}
+
+	var product QuadraticExtensionAlgebraVariable
+	for i := 0; i < D; i++ {
+		acc := p.InnerProductExtension(NewVariable(W), ZeroExtension(), innerW[i])
+		product[i] = p.InnerProductExtension(One(), acc, inner[i])
+	}
+
+	return product
+}
+
+func (p *Chip) ScalarMulExtensionAlgebra(
+	a QuadraticExtensionVariable,
+	b QuadraticExtensionAlgebraVariable,
+) QuadraticExtensionAlgebraVariable {
+	var product QuadraticExtensionAlgebraVariable
+	for i := 0; i < D; i++ {
+		product[i] = p.MulExtension(a, b[i])
+	}
+	return product
+}
+
+func (p *Chip) PartialInterpolateExtAlgebra(
+	domain []goldilocks.Element,
+	values []QuadraticExtensionAlgebraVariable,
+	barycentricWeights []goldilocks.Element,
+	point QuadraticExtensionAlgebraVariable,
+	initialEval QuadraticExtensionAlgebraVariable,
+	initialPartialProd QuadraticExtensionAlgebraVariable,
+) (QuadraticExtensionAlgebraVariable, QuadraticExtensionAlgebraVariable) {
+	n := len(values)
+	if n == 0 {
+		panic("Cannot interpolate with no values")
+	}
+	if n != len(domain) {
+		panic("Domain and values must have the same length")
+	}
+	if n != len(barycentricWeights) {
+		panic("Domain and barycentric weights must have the same length")
+	}
+
+	newEval := initialEval
+	newPartialProd := initialPartialProd
+	for i := 0; i < n; i++ {
+		val := values[i]
+		x := domain[i]
+		xField := NewVariable(x)
+		xQE := xField.ToQuadraticExtension()
+		xQEAlgebra := xQE.ToQuadraticExtensionAlgebra()
+		weight := NewVariable(barycentricWeights[i].Uint64()).ToQuadraticExtension()
+		term := p.SubExtensionAlgebra(point, xQEAlgebra)
+		weightedVal := p.ScalarMulExtensionAlgebra(weight, val)
+		newEval = p.MulExtensionAlgebra(newEval, term)
+		tmp := p.MulExtensionAlgebra(weightedVal, newPartialProd)
+		newEval = p.AddExtensionAlgebra(newEval, tmp)
+		newPartialProd = p.MulExtensionAlgebra(newPartialProd, term)
+	}
+
+	return newEval, newPartialProd
+}
--- a/goldilocks/quadratic_extension_algebra_test.go
+++ b/goldilocks/quadratic_extension_algebra_test.go
@@ -0,0 +1 @@
+package goldilocks
--- a/goldilocks/quadratic_extension_test.go
+++ b/goldilocks/quadratic_extension_test.go
@@ -0,0 +1,94 @@
+package goldilocks
+
+import (
+	"testing"
+
+	"github.com/consensys/gnark-crypto/ecc"
+	"github.com/consensys/gnark/frontend"
+	"github.com/consensys/gnark/test"
+)
+
+type TestQuadraticExtensionMulCircuit struct {
+	Operand1       QuadraticExtensionVariable
+	Operand2       QuadraticExtensionVariable
+	ExpectedResult QuadraticExtensionVariable
+}
+
+func (c *TestQuadraticExtensionMulCircuit) Define(api frontend.API) error {
+	glApi := NewChip(api)
+	actualRes := glApi.MulExtension(c.Operand1, c.Operand2)
+	glApi.AssertIsEqual(actualRes[0], c.ExpectedResult[0])
+	glApi.AssertIsEqual(actualRes[1], c.ExpectedResult[1])
+	return nil
+}
+
+func TestQuadraticExtensionMul4(t *testing.T) {
+	assert := test.NewAssert(t)
+	operand1 := QuadraticExtensionVariable{
+		NewVariable("4994088319481652598"),
+		NewVariable("16489566008211790727"),
+	}
+	operand2 := QuadraticExtensionVariable{
+		NewVariable("3797605683985595697"),
+		NewVariable("13424401189265534004"),
+	}
+	expectedResult := QuadraticExtensionVariable{
+		NewVariable("15052319864161058789"),
+		NewVariable("16841416332519902625"),
+	}
+	circuit := TestQuadraticExtensionMulCircuit{
+		Operand1:       operand1,
+		Operand2:       operand2,
+		ExpectedResult: expectedResult,
+	}
+	witness := TestQuadraticExtensionMulCircuit{
+		Operand1:       operand1,
+		Operand2:       operand2,
+		ExpectedResult: expectedResult,
+	}
+	err := test.IsSolved(&circuit, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
+
+// Test for quadratic extension division
+type TestQuadraticExtensionDivCircuit struct {
+	Operand1       QuadraticExtensionVariable
+	Operand2       QuadraticExtensionVariable
+	ExpectedResult QuadraticExtensionVariable
+}
+
+func (c *TestQuadraticExtensionDivCircuit) Define(api frontend.API) error {
+	glAPI := NewChip(api)
+	actualRes := glAPI.DivExtension(c.Operand1, c.Operand2)
+	glAPI.AssertIsEqual(actualRes[0], c.ExpectedResult[0])
+	glAPI.AssertIsEqual(actualRes[1], c.ExpectedResult[1])
+	return nil
+}
+
+func TestQuadraticExtensionDiv(t *testing.T) {
+	assert := test.NewAssert(t)
+	operand1 := QuadraticExtensionVariable{
+		NewVariable("4994088319481652598"),
+		NewVariable("16489566008211790727"),
+	}
+	operand2 := QuadraticExtensionVariable{
+		NewVariable("7166004739148609569"),
+		NewVariable("14655965871663555016"),
+	}
+	expectedResult := QuadraticExtensionVariable{
+		NewVariable("15052319864161058789"),
+		NewVariable("16841416332519902625"),
+	}
+	circuit := TestQuadraticExtensionDivCircuit{
+		Operand1:       operand1,
+		Operand2:       operand2,
+		ExpectedResult: expectedResult,
+	}
+	witness := TestQuadraticExtensionDivCircuit{
+		Operand1:       operand1,
+		Operand2:       operand2,
+		ExpectedResult: expectedResult,
+	}
+	err := test.IsSolved(&circuit, &witness, ecc.BN254.ScalarField())
+	assert.NoError(err)
+}
--- a/goldilocks/utils.go
+++ b/goldilocks/utils.go
@@ -0,0 +1,45 @@
+package goldilocks
+
+import (
+	"math/big"
+
+	"github.com/consensys/gnark/frontend"
+)
+
+func StrArrayToBigIntArray(input []string) []big.Int {
+	var output []big.Int
+	for i := 0; i < len(input); i++ {
+		a := new(big.Int)
+		a, _ = a.SetString(input[i], 10)
+		output = append(output, *a)
+	}
+	return output
+}
+
+func StrArrayToFrontendVariableArray(input []string) []frontend.Variable {
+	var output []frontend.Variable
+	for i := 0; i < len(input); i++ {
+		output = append(output, frontend.Variable(input[i]))
+	}
+	return output
+}
+
+func Uint64ArrayToVariableArray(input []uint64) []Variable {
+	var output []Variable
+	for i := 0; i < len(input); i++ {
+		output = append(output, NewVariable(input[i]))
+	}
+	return output
+}
+
+func Uint64ArrayToQuadraticExtension(input []uint64) QuadraticExtensionVariable {
+	return NewQuadraticExtensionVariable(NewVariable(input[0]), NewVariable(input[1]))
+}
+
+func Uint64ArrayToQuadraticExtensionArray(input [][]uint64) []QuadraticExtensionVariable {
+	var output []QuadraticExtensionVariable
+	for i := 0; i < len(input); i++ {
+		output = append(output, NewQuadraticExtensionVariable(NewVariable(input[i][0]), NewVariable(input[i][1])))
+	}
+	return output
+}