snarkjs/src/polfield.js
src/polfield.jsRaw
1/*
2 Copyright 2018 0kims association.
3
4 This file is part of snarkjs.
5
6 snarkjs is a free software: you can redistribute it and/or
7 modify it under the terms of the GNU General Public License as published by the
8 Free Software Foundation, either version 3 of the License, or (at your option)
9 any later version.
10
11 snarkjs is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
14 more details.
15
16 You should have received a copy of the GNU General Public License along with
17 snarkjs. If not, see <https://www.gnu.org/licenses/>.
18*/
19
20/*
21 This library does operations on polynomials with coefficients in a field F.
22
23 A polynomial P(x) = p0 + p1 * x + p2 * x^2 + ... + pn * x^n is represented
24 by the array [ p0, p1, p2, ... , pn ].
25 */
26
27const bigInt = require("./bigint.js");
28const assert = require("assert");
29
30class PolField {
31 constructor (F) {
32 this.F = F;
33
34 const q = this.F.q;
35 let rem = q.sub(bigInt(1));
36 let s = 0;
37 while (!rem.isOdd()) {
38 s ++;
39 rem = rem.shr(1);
40 }
41
42 const five = this.F.add(this.F.add(this.F.two, this.F.two), this.F.one);
43
44 this.w = new Array(s+1);
45 this.wi = new Array(s+1);
46 this.w[s] = this.F.exp(five, rem);
47 this.wi[s] = this.F.inverse(this.w[s]);
48
49 let n=s-1;
50 while (n>=0) {
51 this.w[n] = this.F.square(this.w[n+1]);
52 this.wi[n] = this.F.square(this.wi[n+1]);
53 n--;
54 }
55
56
57 this.roots = [];
58/* for (let i=0; i<16; i++) {
59 let r = this.F.one;
60 n = 1 << i;
61 const rootsi = new Array(n);
62 for (let j=0; j<n; j++) {
63 rootsi[j] = r;
64 r = this.F.mul(r, this.w[i]);
65 }
66
67 this.roots.push(rootsi);
68 }
69 */
70 this._setRoots(15);
71 }
72
73 _setRoots(n) {
74 for (let i=n; (i>=0) && (!this.roots[i]); i--) {
75 let r = this.F.one;
76 const nroots = 1 << i;
77 const rootsi = new Array(nroots);
78 for (let j=0; j<nroots; j++) {
79 rootsi[j] = r;
80 r = this.F.mul(r, this.w[i]);
81 }
82
83 this.roots[i] = rootsi;
84 }
85 }
86
87 add(a, b) {
88 const m = Math.max(a.length, b.length);
89 const res = new Array(m);
90 for (let i=0; i<m; i++) {
91 res[i] = this.F.add(a[i] || this.F.zero, b[i] || this.F.zero);
92 }
93 return this.reduce(res);
94 }
95
96 double(a) {
97 return this.add(a,a);
98 }
99
100 sub(a, b) {
101 const m = Math.max(a.length, b.length);
102 const res = new Array(m);
103 for (let i=0; i<m; i++) {
104 res[i] = this.F.sub(a[i] || this.F.zero, b[i] || this.F.zero);
105 }
106 return this.reduce(res);
107 }
108
109 mulScalar(p, b) {
110 if (this.F.isZero(b)) return [];
111 if (this.F.equals(b, this.F.one)) return p;
112 const res = new Array(p.length);
113 for (let i=0; i<p.length; i++) {
114 res[i] = this.F.mul(p[i], b);
115 }
116 return res;
117 }
118
119
120
121 mul(a, b) {
122 if (a.length == 0) return [];
123 if (b.length == 0) return [];
124 if (a.length == 1) return this.mulScalar(b, a[0]);
125 if (b.length == 1) return this.mulScalar(a, b[0]);
126
127 if (b.length > a.length) {
128 [b, a] = [a, b];
129 }
130
131 if ((b.length <= 2) || (b.length < log2(a.length))) {
132 return this.mulNormal(a,b);
133 } else {
134 return this.mulFFT(a,b);
135 }
136 }
137
138 mulNormal(a, b) {
139 let res = [];
140 b = this.affine(b);
141 for (let i=0; i<b.length; i++) {
142 res = this.add(res, this.scaleX(this.mulScalar(a, b[i]), i) );
143 }
144 return res;
145 }
146
147 mulFFT(a,b) {
148 const longestN = Math.max(a.length, b.length);
149 const bitsResult = log2(longestN-1)+2;
150 this._setRoots(bitsResult);
151
152 const m = 1 << bitsResult;
153 const ea = this.extend(a,m);
154 const eb = this.extend(b,m);
155
156 const ta = __fft(this, ea, bitsResult, 0, 1, false);
157 const tb = __fft(this, eb, bitsResult, 0, 1, false);
158
159 const tres = new Array(m);
160
161 for (let i=0; i<m; i++) {
162 tres[i] = this.F.mul(ta[i], tb[i]);
163 }
164
165 const res = __fft(this, tres, bitsResult, 0, 1, true);
166
167 const twoinvm = this.F.inverse( this.F.mulScalar(this.F.one, m) );
168 const resn = new Array(m);
169 for (let i=0; i<m; i++) {
170 resn[i] = this.F.mul(res[(m-i)%m], twoinvm);
171 }
172
173 return this.reduce(this.affine(resn));
174 }
175
176
177
178 square(a) {
179 return this.mul(a,a);
180 }
181
182 scaleX(p, n) {
183 if (n==0) {
184 return p;
185 } else if (n>0) {
186 const z = new Array(n).fill(this.F.zero);
187 return z.concat(p);
188 } else {
189 if (-n >= p.length) return [];
190 return p.slice(-n);
191 }
192 }
193
194 eval2(p, x) {
195 let v = this.F.zero;
196 let ix = this.F.one;
197 for (let i=0; i<p.length; i++) {
198 v = this.F.add(v, this.F.mul(p[i], ix));
199 ix = this.F.mul(ix, x);
200 }
201 return v;
202 }
203
204 eval(p,x) {
205 const F = this.F;
206 if (p.length == 0) return F.zero;
207 const m = this._next2Power(p.length);
208 const ep = this.extend(p, m);
209
210 return _eval(ep, x, 0, 1, m);
211
212 function _eval(p, x, offset, step, n) {
213 if (n==1) return p[offset];
214 const newX = F.square(x);
215 const res= F.add(
216 _eval(p, newX, offset, step << 1, n >> 1),
217 F.mul(
218 x,
219 _eval(p, newX, offset+step , step << 1, n >> 1)));
220 return res;
221 }
222 }
223
224 lagrange(points) {
225 let roots = [this.F.one];
226 for (let i=0; i<points.length; i++) {
227 roots = this.mul(roots, [this.F.neg(points[i][0]), this.F.one]);
228 }
229
230 let sum = [];
231 for (let i=0; i<points.length; i++) {
232 let mpol = this.ruffini(roots, points[i][0]);
233 const factor =
234 this.F.mul(
235 this.F.inverse(this.eval(mpol, points[i][0])),
236 points[i][1]);
237 mpol = this.mulScalar(mpol, factor);
238 sum = this.add(sum, mpol);
239 }
240 return sum;
241 }
242
243
244 fft(p) {
245 if (p.length <= 1) return p;
246 const bits = log2(p.length-1)+1;
247 this._setRoots(bits);
248
249 const m = 1 << bits;
250 const ep = this.extend(p, m);
251 const res = __fft(this, ep, bits, 0, 1);
252 return res;
253 }
254
255 ifft(p) {
256
257 if (p.length <= 1) return p;
258 const bits = log2(p.length-1)+1;
259 this._setRoots(bits);
260 const m = 1 << bits;
261 const ep = this.extend(p, m);
262 const res = __fft(this, ep, bits, 0, 1);
263
264 const twoinvm = this.F.inverse( this.F.mulScalar(this.F.one, m) );
265 const resn = new Array(m);
266 for (let i=0; i<m; i++) {
267 resn[i] = this.F.mul(res[(m-i)%m], twoinvm);
268 }
269
270 return resn;
271
272 }
273
274
275 _fft(pall, bits, offset, step) {
276
277 const n = 1 << bits;
278 if (n==1) {
279 return [ pall[offset] ];
280 }
281
282 const ndiv2 = n >> 1;
283 const p1 = this._fft(pall, bits-1, offset, step*2);
284 const p2 = this._fft(pall, bits-1, offset+step, step*2);
285
286 const out = new Array(n);
287
288 let m= this.F.one;
289 for (let i=0; i<ndiv2; i++) {
290 out[i] = this.F.add(p1[i], this.F.mul(m, p2[i]));
291 out[i+ndiv2] = this.F.sub(p1[i], this.F.mul(m, p2[i]));
292 m = this.F.mul(m, this.w[bits]);
293 }
294
295 return out;
296 }
297
298 extend(p, e) {
299 if (e == p.length) return p;
300 const z = new Array(e-p.length).fill(this.F.zero);
301
302 return p.concat(z);
303 }
304
305 reduce(p) {
306 if (p.length == 0) return p;
307 if (! this.F.isZero(p[p.length-1]) ) return p;
308 let i=p.length-1;
309 while( i>0 && this.F.isZero(p[i]) ) i--;
310 return p.slice(0, i+1);
311 }
312
313 affine(p) {
314 for (let i=0; i<p.length; i++) {
315 p[i] = this.F.affine(p[i]);
316 }
317 return p;
318 }
319
320 equals(a, b) {
321 const pa = this.reduce(this.affine(a));
322 const pb = this.reduce(this.affine(b));
323
324 if (pa.length != pb.length) return false;
325 for (let i=0; i<pb.length; i++) {
326 if (!this.F.equals(pa[i], pb[i])) return false;
327 }
328
329 return true;
330 }
331
332 ruffini(p, r) {
333 const res = new Array(p.length-1);
334 res[res.length-1] = p[p.length-1];
335 for (let i = res.length-2; i>=0; i--) {
336 res[i] = this.F.add(this.F.mul(res[i+1], r), p[i+1]);
337 }
338 return res;
339 }
340
341 _next2Power(v) {
342 v--;
343 v |= v >> 1;
344 v |= v >> 2;
345 v |= v >> 4;
346 v |= v >> 8;
347 v |= v >> 16;
348 v++;
349 return v;
350 }
351
352 toString(p) {
353 const ap = this.affine(p);
354 let S = "";
355 for (let i=ap.length-1; i>=0; i--) {
356 if (!this.F.isZero(p[i])) {
357 if (S!="") S += " + ";
358 S = S + p[i].toString(10);
359 if (i>0) {
360 S = S + "x";
361 if (i>1) {
362 S = S + "^" +i;
363 }
364 }
365 }
366 }
367 return S;
368 }
369
370
371 _reciprocal(p, bits) {
372 const k = 1 << bits;
373 if (k==1) {
374 return [ this.F.inverse(p[0]) ];
375 }
376 const np = this.scaleX(p, -k/2);
377 const q = this._reciprocal(np, bits-1);
378 const a = this.scaleX(this.double(q), 3*k/2-2);
379 const b = this.mul( this.square(q), p);
380
381 return this.scaleX(this.sub(a,b), -(k-2));
382 }
383
384 // divides x^m / v
385 _div2(m, v) {
386 const kbits = log2(v.length-1)+1;
387 const k = 1 << kbits;
388
389 const scaleV = k - v.length;
390
391 // rec = x^(k - 2) / v* x^scaleV =>
392 // rec = x^(k-2-scaleV)/ v
393 //
394 // res = x^m/v = x^(m + (2*k-2 - scaleV) - (2*k-2 - scaleV)) /v =>
395 // res = rec * x^(m - (2*k-2 - scaleV)) =>
396 // res = rec * x^(m - 2*k + 2 + scaleV)
397
398 const rec = this._reciprocal(this.scaleX(v, scaleV), kbits);
399 const res = this.scaleX(rec, m - 2*k + 2 + scaleV);
400
401 return res;
402 }
403
404 div(_u, _v) {
405 if (_u.length < _v.length) return [];
406 const kbits = log2(_v.length-1)+1;
407 const k = 1 << kbits;
408
409 const u = this.scaleX(_u, k-_v.length);
410 const v = this.scaleX(_v, k-_v.length);
411
412 const n = v.length-1;
413 let m = u.length-1;
414
415 const s = this._reciprocal(v, kbits);
416 let t;
417 if (m>2*n) {
418 t = this.sub(this.scaleX([this.F.one], 2*n), this.mul(s, v));
419 }
420
421 let q = [];
422 let rem = u;
423 let us, ut;
424 let finish = false;
425
426 while (!finish) {
427 us = this.mul(rem, s);
428 q = this.add(q, this.scaleX(us, -2*n));
429
430 if ( m > 2*n ) {
431 ut = this.mul(rem, t);
432 rem = this.scaleX(ut, -2*n);
433 m = rem.length-1;
434 } else {
435 finish = true;
436 }
437 }
438
439 return q;
440 }
441
442
443 // returns the ith nth-root of one
444 oneRoot(n, i) {
445 let nbits = log2(n-1)+1;
446 let res = this.F.one;
447 let r = i;
448
449 assert(i<n);
450 assert(1<<nbits === n);
451
452 while (r>0) {
453 if (r & 1 == 1) {
454 res = this.F.mul(res, this.w[nbits]);
455 }
456 r = r >> 1;
457 nbits --;
458 }
459 return res;
460 }
461
462 computeVanishingPolinomial(bits, t) {
463 const m = 1 << bits;
464 return this.F.sub(this.F.exp(t, bigInt(m)), this.F.one);
465 }
466
467 evaluateLagrangePolynomials(bits, t) {
468 const m= 1 << bits;
469 const tm = this.F.exp(t, bigInt(m));
470 const u= new Array(m).fill(this.F.zero);
471 this._setRoots(bits);
472 const omega = this.w[bits];
473
474 if (this.F.equals(tm, this.F.one)) {
475 for (let i = 0; i < m; i++) {
476 if (this.F.equals(this.roots[bits][0],t)) { // i.e., t equals omega^i
477 u[i] = this.F.one;
478 return u;
479 }
480 }
481 }
482
483 const z = this.F.sub(tm, this.F.one);
484// let l = this.F.mul(z, this.F.exp(this.F.twoinv, m));
485 let l = this.F.mul(z, this.F.inverse(bigInt(m)));
486 for (let i = 0; i < m; i++) {
487 u[i] = this.F.mul(l, this.F.inverse(this.F.sub(t,this.roots[bits][i])));
488 l = this.F.mul(l, omega);
489 }
490
491 return u;
492 }
493
494 log2(V) {
495 return log2(V);
496 }
497}
498
499function log2( V )
500{
501 return( ( ( V & 0xFFFF0000 ) !== 0 ? ( V &= 0xFFFF0000, 16 ) : 0 ) | ( ( V & 0xFF00FF00 ) !== 0 ? ( V &= 0xFF00FF00, 8 ) : 0 ) | ( ( V & 0xF0F0F0F0 ) !== 0 ? ( V &= 0xF0F0F0F0, 4 ) : 0 ) | ( ( V & 0xCCCCCCCC ) !== 0 ? ( V &= 0xCCCCCCCC, 2 ) : 0 ) | ( ( V & 0xAAAAAAAA ) !== 0 ) );
502}
503
504
505function __fft(PF, pall, bits, offset, step) {
506
507 const n = 1 << bits;
508 if (n==1) {
509 return [ pall[offset] ];
510 } else if (n==2) {
511 return [
512 PF.F.add(pall[offset], pall[offset + step]),
513 PF.F.sub(pall[offset], pall[offset + step])];
514 }
515
516 const ndiv2 = n >> 1;
517 const p1 = __fft(PF, pall, bits-1, offset, step*2);
518 const p2 = __fft(PF, pall, bits-1, offset+step, step*2);
519
520 const out = new Array(n);
521
522 for (let i=0; i<ndiv2; i++) {
523 out[i] = PF.F.add(p1[i], PF.F.mul(PF.roots[bits][i], p2[i]));
524 out[i+ndiv2] = PF.F.sub(p1[i], PF.F.mul(PF.roots[bits][i], p2[i]));
525 }
526
527 return out;
528}
529
530
531module.exports = PolField;