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