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