kzg-commitments-study/arithmetic.go
arithmetic.goRaw
1package kzg
2
3import (
4 "bytes"
5 "crypto/rand"
6 "fmt"
7 "math/big"
8 "strconv"
9
10 bn256 "github.com/ethereum/go-ethereum/crypto/bn256/cloudflare"
11)
12
13// Q is the order of the integer field (Zq) that fits inside the snark
14var Q, _ = new(big.Int).SetString(
15 "21888242871839275222246405745257275088696311157297823662689037894645226208583", 10)
16
17// R is the mod of the finite field
18var R, _ = new(big.Int).SetString(
19 "21888242871839275222246405745257275088548364400416034343698204186575808495617", 10)
20
21func randBigInt() (*big.Int, error) {
22 maxbits := R.BitLen()
23 b := make([]byte, (maxbits/8)-1)
24 _, err := rand.Read(b)
25 if err != nil {
26 return nil, err
27 }
28 r := new(big.Int).SetBytes(b)
29 rq := new(big.Int).Mod(r, R)
30
31 return rq, nil
32}
33
34func arrayOfZeroes(n int) []*big.Int {
35 r := make([]*big.Int, n)
36 for i := 0; i < n; i++ {
37 r[i] = new(big.Int).SetInt64(0)
38 }
39 return r[:]
40}
41
42//nolint:deadcode,unused
43func arrayOfZeroesG1(n int) []*bn256.G1 {
44 r := make([]*bn256.G1, n)
45 for i := 0; i < n; i++ {
46 r[i] = new(bn256.G1).ScalarBaseMult(big.NewInt(0))
47 }
48 return r[:]
49}
50
51//nolint:deadcode,unused
52func arrayOfZeroesG2(n int) []*bn256.G2 {
53 r := make([]*bn256.G2, n)
54 for i := 0; i < n; i++ {
55 r[i] = new(bn256.G2).ScalarBaseMult(big.NewInt(0))
56 }
57 return r[:]
58}
59
60func compareBigIntArray(a, b []*big.Int) bool {
61 if len(a) != len(b) {
62 return false
63 }
64 for i := 0; i < len(a); i++ {
65 if a[i] != b[i] {
66 return false
67 }
68 }
69 return true
70}
71
72//nolint:deadcode,unused
73func checkArrayOfZeroes(a []*big.Int) bool {
74 z := arrayOfZeroes(len(a))
75 return compareBigIntArray(a, z)
76}
77
78func fAdd(a, b *big.Int) *big.Int {
79 ab := new(big.Int).Add(a, b)
80 return ab.Mod(ab, R)
81}
82
83func fSub(a, b *big.Int) *big.Int {
84 ab := new(big.Int).Sub(a, b)
85 return new(big.Int).Mod(ab, R)
86}
87
88func fMul(a, b *big.Int) *big.Int {
89 ab := new(big.Int).Mul(a, b)
90 return ab.Mod(ab, R)
91}
92
93func fDiv(a, b *big.Int) *big.Int {
94 ab := new(big.Int).Mul(a, new(big.Int).ModInverse(b, R))
95 return new(big.Int).Mod(ab, R)
96}
97
98func fNeg(a *big.Int) *big.Int {
99 return new(big.Int).Mod(new(big.Int).Neg(a), R)
100}
101
102//nolint:deadcode,unused
103func fInv(a *big.Int) *big.Int {
104 return new(big.Int).ModInverse(a, R)
105}
106
107func fExp(base *big.Int, e *big.Int) *big.Int {
108 res := big.NewInt(1)
109 rem := new(big.Int).Set(e)
110 exp := base
111
112 for !bytes.Equal(rem.Bytes(), big.NewInt(int64(0)).Bytes()) {
113 // if BigIsOdd(rem) {
114 if rem.Bit(0) == 1 { // .Bit(0) returns 1 when is odd
115 res = fMul(res, exp)
116 }
117 exp = fMul(exp, exp)
118 rem.Rsh(rem, 1)
119 }
120 return res
121}
122
123func max(a, b int) int {
124 if a > b {
125 return a
126 }
127 return b
128}
129
130func polynomialAdd(a, b []*big.Int) []*big.Int {
131 r := arrayOfZeroes(max(len(a), len(b)))
132 for i := 0; i < len(a); i++ {
133 r[i] = fAdd(r[i], a[i])
134 }
135 for i := 0; i < len(b); i++ {
136 r[i] = fAdd(r[i], b[i])
137 }
138 return r
139}
140
141func polynomialSub(a, b []*big.Int) []*big.Int {
142 r := arrayOfZeroes(max(len(a), len(b)))
143 for i := 0; i < len(a); i++ {
144 r[i] = fAdd(r[i], a[i])
145 }
146 for i := 0; i < len(b); i++ {
147 r[i] = fSub(r[i], b[i])
148 }
149 return r
150}
151
152func polynomialMul(a, b []*big.Int) []*big.Int {
153 r := arrayOfZeroes(len(a) + len(b) - 1)
154 for i := 0; i < len(a); i++ {
155 for j := 0; j < len(b); j++ {
156 r[i+j] = fAdd(r[i+j], fMul(a[i], b[j]))
157 }
158 }
159 return r
160}
161
162func polynomialDiv(a, b []*big.Int) ([]*big.Int, []*big.Int) {
163 // https://en.wikipedia.org/wiki/Division_algorithm
164 r := arrayOfZeroes(len(a) - len(b) + 1)
165 rem := a
166 for len(rem) >= len(b) {
167 l := fDiv(rem[len(rem)-1], b[len(b)-1])
168 pos := len(rem) - len(b)
169 r[pos] = l
170 aux := arrayOfZeroes(pos)
171 aux1 := append(aux, l)
172 aux2 := polynomialSub(rem, polynomialMul(b, aux1))
173 rem = aux2[:len(aux2)-1]
174 }
175 return r, rem
176}
177
178func polynomialMulByConstant(a []*big.Int, c *big.Int) []*big.Int {
179 for i := 0; i < len(a); i++ {
180 a[i] = fMul(a[i], c)
181 }
182 return a
183}
184func polynomialDivByConstant(a []*big.Int, c *big.Int) []*big.Int {
185 for i := 0; i < len(a); i++ {
186 a[i] = fDiv(a[i], c)
187 }
188 return a
189}
190
191// polynomialEval evaluates the polinomial over the Finite Field at the given value x
192func polynomialEval(p []*big.Int, x *big.Int) *big.Int {
193 r := big.NewInt(int64(0))
194 for i := 0; i < len(p); i++ {
195 xi := fExp(x, big.NewInt(int64(i)))
196 elem := fMul(p[i], xi)
197 r = fAdd(r, elem)
198 }
199 return r
200}
201
202// newPolZeroAt generates a new polynomial that has value zero at the given value
203func newPolZeroAt(pointPos, totalPoints int, height *big.Int) []*big.Int {
204 fac := 1
205 for i := 1; i < totalPoints+1; i++ {
206 if i != pointPos {
207 fac = fac * (pointPos - i)
208 }
209 }
210 facBig := big.NewInt(int64(fac))
211 hf := fDiv(height, facBig)
212 r := []*big.Int{hf}
213 for i := 1; i < totalPoints+1; i++ {
214 if i != pointPos {
215 ineg := big.NewInt(int64(-i))
216 b1 := big.NewInt(int64(1))
217 r = polynomialMul(r, []*big.Int{ineg, b1})
218 }
219 }
220 return r
221}
222
223// zeroPolynomial returns the zero polynomial:
224// z(x) = (x - z_0) (x - z_1) ... (x - z_{k-1})
225func zeroPolynomial(zs []*big.Int) []*big.Int {
226 z := []*big.Int{fNeg(zs[0]), big.NewInt(1)} // (x - z0)
227 for i := 1; i < len(zs); i++ {
228 z = polynomialMul(z, []*big.Int{fNeg(zs[i]), big.NewInt(1)}) // (x - zi)
229 }
230 return z
231}
232
233var sNums = map[string]string{
234 "0": "⁰",
235 "1": "¹",
236 "2": "²",
237 "3": "³",
238 "4": "⁴",
239 "5": "⁵",
240 "6": "⁶",
241 "7": "⁷",
242 "8": "⁸",
243 "9": "⁹",
244}
245
246func intToSNum(n int) string {
247 s := strconv.Itoa(n)
248 sN := ""
249 for i := 0; i < len(s); i++ {
250 sN += sNums[string(s[i])]
251 }
252 return sN
253}
254
255// PolynomialToString converts a polynomial represented by a *big.Int array,
256// into its string human readable representation
257func PolynomialToString(p []*big.Int) string {
258 s := ""
259 for i := len(p) - 1; i >= 1; i-- {
260 if bytes.Equal(p[i].Bytes(), big.NewInt(1).Bytes()) {
261 s += fmt.Sprintf("x%s + ", intToSNum(i))
262 } else if !bytes.Equal(p[i].Bytes(), big.NewInt(0).Bytes()) {
263 s += fmt.Sprintf("%sx%s + ", p[i], intToSNum(i))
264 }
265 }
266 s += p[0].String()
267 return s
268}
269
270// LagrangeInterpolation implements the Lagrange interpolation:
271// https://en.wikipedia.org/wiki/Lagrange_polynomial
272func LagrangeInterpolation(x, y []*big.Int) ([]*big.Int, error) {
273 // p(x) will be the interpoled polynomial
274 // var p []*big.Int
275 if len(x) != len(y) {
276 return nil, fmt.Errorf("len(x)!=len(y): %d, %d", len(x), len(y))
277 }
278 p := arrayOfZeroes(len(x))
279 k := len(x)
280
281 for j := 0; j < k; j++ {
282 // jPol is the Lagrange basis polynomial for each point
283 var jPol []*big.Int
284 for m := 0; m < k; m++ {
285 // if x[m] == x[j] {
286 if m == j {
287 continue
288 }
289 // numerator & denominator of the current iteration
290 num := []*big.Int{fNeg(x[m]), big.NewInt(1)} // (x^1 - x_m)
291 den := fSub(x[j], x[m]) // x_j-x_m
292 mPol := polynomialDivByConstant(num, den)
293 if len(jPol) == 0 {
294 // first j iteration
295 jPol = mPol
296 continue
297 }
298 jPol = polynomialMul(jPol, mPol)
299 }
300 p = polynomialAdd(p, polynomialMulByConstant(jPol, y[j]))
301 }
302
303 return p, nil
304}
305
306// TODO add method to 'clean' the polynomial, to remove right-zeroes