1package field
  2
  3import (
  4	"math/bits"
  5
  6	"github.com/consensys/gnark-crypto/field/goldilocks"
  7	"github.com/consensys/gnark/frontend"
  8)
  9
 10const D = 2
 11
 12type QuadraticExtension = [2]F
 13type QEAlgebra = [D]QuadraticExtension
 14
 15type QuadraticExtensionAPI struct {
 16	api      frontend.API
 17	fieldAPI FieldAPI
 18
 19	W        F
 20	DTH_ROOT F
 21
 22	ONE_QE  QuadraticExtension
 23	ZERO_QE QuadraticExtension
 24
 25	ZERO_QE_ALGEBRA QEAlgebra
 26}
 27
 28func NewQuadraticExtensionAPI(api frontend.API, fieldAPI FieldAPI) *QuadraticExtensionAPI {
 29	var ZERO_QE = QuadraticExtension{ZERO_F, ZERO_F}
 30
 31	var ZERO_QE_ALGEBRA QEAlgebra
 32
 33	for i := 0; i < D; i++ {
 34		ZERO_QE_ALGEBRA[i] = ZERO_QE
 35	}
 36
 37	return &QuadraticExtensionAPI{
 38		api:      api,
 39		fieldAPI: fieldAPI,
 40
 41		W:        NewFieldConst(7),
 42		DTH_ROOT: NewFieldConst(18446744069414584320),
 43
 44		ONE_QE:  QuadraticExtension{ONE_F, ZERO_F},
 45		ZERO_QE: ZERO_QE,
 46
 47		ZERO_QE_ALGEBRA: ZERO_QE_ALGEBRA,
 48	}
 49}
 50
 51func (c *QuadraticExtensionAPI) SquareExtension(a QuadraticExtension) QuadraticExtension {
 52	return c.MulExtension(a, a)
 53}
 54
 55func (c *QuadraticExtensionAPI) MulExtension(a QuadraticExtension, b QuadraticExtension) QuadraticExtension {
 56	c_0 := c.fieldAPI.Add(c.fieldAPI.Mul(a[0], b[0]), c.fieldAPI.Mul(c.fieldAPI.Mul(c.W, a[1]), b[1]))
 57	c_1 := c.fieldAPI.Add(c.fieldAPI.Mul(a[0], b[1]), c.fieldAPI.Mul(a[1], b[0]))
 58	return QuadraticExtension{c_0, c_1}
 59}
 60
 61func (c *QuadraticExtensionAPI) AddExtension(a QuadraticExtension, b QuadraticExtension) QuadraticExtension {
 62	c_0 := c.fieldAPI.Add(a[0], b[0])
 63	c_1 := c.fieldAPI.Add(a[1], b[1])
 64	return QuadraticExtension{c_0, c_1}
 65}
 66
 67func (c *QuadraticExtensionAPI) SubExtension(a QuadraticExtension, b QuadraticExtension) QuadraticExtension {
 68	c_0 := c.fieldAPI.Sub(a[0], b[0])
 69	c_1 := c.fieldAPI.Sub(a[1], b[1])
 70	return QuadraticExtension{c_0, c_1}
 71}
 72
 73func (c *QuadraticExtensionAPI) DivExtension(a QuadraticExtension, b QuadraticExtension) QuadraticExtension {
 74	return c.MulExtension(a, c.InverseExtension(b))
 75}
 76
 77func (c *QuadraticExtensionAPI) IsZero(a QuadraticExtension) frontend.Variable {
 78	return c.api.Mul(IsZero(c.api, c.fieldAPI, a[0]), IsZero(c.api, c.fieldAPI, a[1]))
 79}
 80
 81// TODO: Instead of calculating the inverse within the circuit, can witness the
 82// inverse and assert that a_inverse * a = 1.  Should reduce # of constraints.
 83func (c *QuadraticExtensionAPI) InverseExtension(a QuadraticExtension) QuadraticExtension {
 84	// First assert that a doesn't have 0 value coefficients
 85	a0_is_zero := IsZero(c.api, c.fieldAPI, a[0])
 86	a1_is_zero := IsZero(c.api, c.fieldAPI, a[1])
 87
 88	// assert that a0_is_zero OR a1_is_zero == false
 89	c.api.AssertIsEqual(c.api.Mul(a0_is_zero, a1_is_zero), frontend.Variable(0))
 90
 91	a_pow_r_minus_1 := QuadraticExtension{a[0], c.fieldAPI.Mul(a[1], c.DTH_ROOT)}
 92	a_pow_r := c.MulExtension(a_pow_r_minus_1, a)
 93	return c.ScalarMulExtension(a_pow_r_minus_1, c.fieldAPI.Inverse(a_pow_r[0]))
 94}
 95
 96func (c *QuadraticExtensionAPI) ScalarMulExtension(a QuadraticExtension, scalar F) QuadraticExtension {
 97	return QuadraticExtension{c.fieldAPI.Mul(a[0], scalar), c.fieldAPI.Mul(a[1], scalar)}
 98}
 99
100func (c *QuadraticExtensionAPI) VarToQE(a frontend.Variable) QuadraticExtension {
101	return c.FieldToQE(c.fieldAPI.NewElement(a))
102}
103
104func (c *QuadraticExtensionAPI) FieldToQE(a F) QuadraticExtension {
105	return QuadraticExtension{a, ZERO_F}
106}
107
108// / Exponentiate `base` to the power of a known `exponent`.
109func (c *QuadraticExtensionAPI) ExpU64Extension(a QuadraticExtension, exponent uint64) QuadraticExtension {
110	switch exponent {
111	case 0:
112		return c.ONE_QE
113	case 1:
114		return a
115	case 2:
116		return c.SquareExtension(a)
117	default:
118	}
119
120	current := a
121	product := c.ONE_QE
122
123	for i := 0; i < bits.Len64(exponent); i++ {
124		if i != 0 {
125			current = c.SquareExtension(current)
126		}
127
128		if (exponent >> i & 1) != 0 {
129			product = c.MulExtension(product, current)
130		}
131	}
132
133	return product
134}
135
136func (c *QuadraticExtensionAPI) ReduceWithPowers(terms []QuadraticExtension, scalar QuadraticExtension) QuadraticExtension {
137	sum := c.ZERO_QE
138
139	for i := len(terms) - 1; i >= 0; i-- {
140		sum = c.AddExtension(
141			c.MulExtension(
142				sum,
143				scalar,
144			),
145			terms[i],
146		)
147	}
148
149	return sum
150}
151
152func (c *QuadraticExtensionAPI) Select(b frontend.Variable, qe0, qe1 QuadraticExtension) QuadraticExtension {
153	var retQE QuadraticExtension
154
155	for i := 0; i < 2; i++ {
156		retQE[i] = c.fieldAPI.Select(b, qe0[i], qe1[i])
157	}
158
159	return retQE
160}
161
162func (c *QuadraticExtensionAPI) Lookup2(b0 frontend.Variable, b1 frontend.Variable, qe0, qe1, qe2, qe3 QuadraticExtension) QuadraticExtension {
163	var retQE QuadraticExtension
164
165	for i := 0; i < 2; i++ {
166		retQE[i] = c.fieldAPI.Lookup2(b0, b1, qe0[i], qe1[i], qe2[i], qe3[i])
167	}
168
169	return retQE
170}
171
172func (c *QuadraticExtensionAPI) AssertIsEqual(a, b QuadraticExtension) {
173	for i := 0; i < 2; i++ {
174		c.fieldAPI.AssertIsEqual(a[i], b[i])
175	}
176}
177
178func (c *QuadraticExtensionAPI) InnerProductExtension(constant F, startingAcc QuadraticExtension, pairs [][2]QuadraticExtension) QuadraticExtension {
179	acc := startingAcc
180
181	for i := 0; i < len(pairs); i++ {
182		a := pairs[i][0]
183		b := pairs[i][1]
184		mul := c.ScalarMulExtension(a, constant)
185		mul = c.MulExtension(mul, b)
186		acc = c.AddExtension(acc, mul)
187	}
188
189	return acc
190}
191
192/*
193func (c *QuadraticExtensionAPI) Println(a QuadraticExtension) {
194	fmt.Print("Degree 0 coefficient")
195	c.fieldAPI.Println(a[0])
196
197	fmt.Print("Degree 1 coefficient")
198	c.fieldAPI.Println(a[1])
199}
200*/
201
202func (c *QuadraticExtensionAPI) MulExtensionAlgebra(a, b QEAlgebra) QEAlgebra {
203	var inner [D][][2]QuadraticExtension
204	var inner_w [D][][2]QuadraticExtension
205	for i := 0; i < D; i++ {
206		for j := 0; j < D-i; j++ {
207			idx := (i + j) % D
208			inner[idx] = append(inner[idx], [2]QuadraticExtension{a[i], b[j]})
209		}
210		for j := D - i; j < D; j++ {
211			idx := (i + j) % D
212			inner_w[idx] = append(inner_w[idx], [2]QuadraticExtension{a[i], b[j]})
213		}
214	}
215
216	var product QEAlgebra
217	for i := 0; i < D; i++ {
218		acc := c.InnerProductExtension(c.W, c.ZERO_QE, inner_w[i])
219		product[i] = c.InnerProductExtension(ONE_F, acc, inner[i])
220	}
221
222	return product
223}
224
225func (c *QuadraticExtensionAPI) ScalarMulExtensionAlgebra(a QuadraticExtension, b QEAlgebra) QEAlgebra {
226	var product QEAlgebra
227	for i := 0; i < D; i++ {
228		product[i] = c.MulExtension(a, b[i])
229	}
230
231	return product
232}
233
234func (c *QuadraticExtensionAPI) AddExtensionAlgebra(a, b QEAlgebra) QEAlgebra {
235	var sum QEAlgebra
236	for i := 0; i < D; i++ {
237		sum[i] = c.AddExtension(a[i], b[i])
238	}
239
240	return sum
241}
242
243func (c *QuadraticExtensionAPI) SubExtensionAlgebra(a, b QEAlgebra) QEAlgebra {
244	var diff QEAlgebra
245	for i := 0; i < D; i++ {
246		diff[i] = c.SubExtension(a[i], b[i])
247	}
248
249	return diff
250}
251
252func (c *QuadraticExtensionAPI) PartialInterpolateExtAlgebra(
253	domain []goldilocks.Element,
254	values []QEAlgebra,
255	barycentricWeights []goldilocks.Element,
256	point QEAlgebra,
257	initialEval QEAlgebra,
258	initialPartialProd QEAlgebra,
259) (QEAlgebra, QEAlgebra) {
260	n := len(values)
261	if n == 0 {
262		panic("Cannot interpolate with no values")
263	}
264	if n != len(domain) {
265		panic("Domain and values must have the same length")
266	}
267	if n != len(barycentricWeights) {
268		panic("Domain and barycentric weights must have the same length")
269	}
270
271	newEval := initialEval
272	newPartialProd := initialPartialProd
273	for i := 0; i < n; i++ {
274		val := values[i]
275		x := domain[i]
276		xField := NewFieldConst(x.Uint64())
277		xQE := QuadraticExtension{xField, ZERO_F}
278		xQEAlgebra := QEAlgebra{xQE, c.ZERO_QE}
279		weight := QuadraticExtension{NewFieldConst(barycentricWeights[i].Uint64()), ZERO_F}
280		term := c.SubExtensionAlgebra(point, xQEAlgebra)
281		weightedVal := c.ScalarMulExtensionAlgebra(weight, val)
282		newEval = c.MulExtensionAlgebra(newEval, term)
283		tmp := c.MulExtensionAlgebra(weightedVal, newPartialProd)
284		newEval = c.AddExtensionAlgebra(newEval, tmp)
285		newPartialProd = c.MulExtensionAlgebra(newPartialProd, term)
286	}
287
288	return newEval, newPartialProd
289}