snarkjs/src/polfield.js
src/polfield.jsRaw
1/*
2 Copyright 2018 0kims association
3
4 This file is part of zksnark javascript library.
5
6 zksnark javascript library is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 zksnark javascript library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with zksnark javascript library. If not, see <https://www.gnu.org/licenses/>.
18*/
19
20/*
21 This library do operations on polynomials where their coefficients are in field F
22
23 The polynomial P(x) = p0 + p1 * x + p2 * x^2 + p3 * x^3, ...
24 is represented by the array [ p0, p1, p2, p3, ... ]
25 */
26
27const bigInt = require("./bigint.js");
28
29class PolField {
30 constructor (F) {
31 this.F = F;
32
33 const q = this.F.q;
34 let rem = q.sub(bigInt(1));
35 let s = 0;
36 while (!rem.isOdd()) {
37 s ++;
38 rem = rem.shr(1);
39 }
40
41 const five = this.F.add(this.F.add(this.F.two, this.F.two), this.F.one);
42
43 this.w = new Array(s+1);
44 this.wi = new Array(s+1);
45 this.w[s] = this.F.exp(five, rem);
46 this.wi[s] = this.F.inverse(this.w[s]);
47
48 let n=s-1;
49 while (n>=0) {
50 this.w[n] = this.F.square(this.w[n+1]);
51 this.wi[n] = this.F.square(this.wi[n+1]);
52 n--;
53 }
54
55 }
56
57 add(a, b) {
58 const m = Math.max(a.length, b.length);
59 const res = new Array(m);
60 for (let i=0; i<m; i++) {
61 res[i] = this.F.add(a[i] || this.F.zero, b[i] || this.F.zero);
62 }
63 return this.reduce(res);
64 }
65
66 double(a) {
67 return this.add(a,a);
68 }
69
70 sub(a, b) {
71 const m = Math.max(a.length, b.length);
72 const res = new Array(m);
73 for (let i=0; i<m; i++) {
74 res[i] = this.F.sub(a[i] || this.F.zero, b[i] || this.F.zero);
75 }
76 return this.reduce(res);
77 }
78
79 mulScalar(a, b) {
80 if (this.F.isZero(b)) return [];
81 if (this.F.equals(b, this.F.one)) return a;
82 const res = new Array(a.length);
83 for (let i=0; i<a.length; i++) {
84 res[i] = this.F.mul(a[i], b);
85 }
86 return res;
87 }
88
89
90
91 mul(a, b) {
92 if (a.length == 0) return [];
93 if (b.length == 0) return [];
94 if (a.length == 1) return this.mulScalar(b, a[0]);
95 if (b.length == 1) return this.mulScalar(a, b[0]);
96
97 if (b.length > a.length) {
98 [b, a] = [a, b];
99 }
100
101 if ((b.length <= 2) || (b.length < log2(a.length))) {
102 return this.mulNormal(a,b);
103 } else {
104 return this.mulFFT(a,b);
105 }
106 }
107
108 mulNormal(a, b) {
109 let res = [];
110 b = this.affine(b);
111 for (let i=0; i<b.length; i++) {
112 res = this.add(res, this.scaleX(this.mulScalar(a, b[i]), i) );
113 }
114 return res;
115 }
116
117 mulFFT(a,b) {
118 const longestN = Math.max(a.length, b.length);
119 const bitsResult = log2(longestN-1)+2;
120 const m = 1 << bitsResult;
121 const ea = this.extend(a,m);
122 const eb = this.extend(b,m);
123
124 const ta = this._fft(ea, bitsResult, 0, 1, false);
125 const tb = this._fft(eb, bitsResult, 0, 1, false);
126
127 const tres = new Array(m);
128
129 for (let i=0; i<m; i++) {
130 tres[i] = this.F.mul(ta[i], tb[i]);
131 }
132
133 const res = this._fft(tres, bitsResult, 0, 1, true);
134
135 const twoinvm = this.F.inverse( this.F.mulScalar(this.F.one, m) );
136 const resn = new Array(m);
137 for (let i=0; i<m; i++) {
138 resn[i] = this.F.mul(res[(m-i)%m], twoinvm);
139 }
140
141 return this.reduce(this.affine(resn));
142 }
143
144
145
146 square(a) {
147 return this.mul(a,a);
148 }
149
150 scaleX(p, n) {
151 if (n==0) {
152 return p;
153 } else if (n>0) {
154 const z = new Array(n).fill(this.F.zero);
155 return z.concat(p);
156 } else {
157 if (-n >= p.length) return [];
158 return p.slice(-n);
159 }
160 }
161
162 eval(p, x) {
163 let v = this.F.zero;
164 let ix = this.F.one;
165 for (let i=0; i<p.length; i++) {
166 v = this.F.add(v, this.F.mul(p[i], ix));
167 ix = this.F.mul(ix, x);
168 }
169 return v;
170 }
171
172 lagrange(points) {
173 let roots = [this.F.one];
174 for (let i=0; i<points.length; i++) {
175 roots = this.mul(roots, [this.F.neg(points[i][0]), this.F.one]);
176 }
177
178 let sum = [];
179 for (let i=0; i<points.length; i++) {
180 let mpol = this.ruffini(roots, points[i][0]);
181 const factor =
182 this.F.mul(
183 this.F.inverse(this.eval(mpol, points[i][0])),
184 points[i][1]);
185 mpol = this.mulScalar(mpol, factor);
186 sum = this.add(sum, mpol);
187 }
188 return sum;
189 }
190
191 _fft(pall, bits, offset, step) {
192
193 const n = 1 << bits;
194 if (n==1) {
195 return [ pall[offset] ];
196 }
197
198 const ndiv2 = n >> 1;
199 const p1 = this._fft(pall, bits-1, offset, step*2);
200 const p2 = this._fft(pall, bits-1, offset+step, step*2);
201
202 const out = new Array(n);
203
204 let m= this.F.one;
205 for (let i=0; i<ndiv2; i++) {
206 out[i] = this.F.add(p1[i], this.F.mul(m, p2[i]));
207 out[i+ndiv2] = this.F.sub(p1[i], this.F.mul(m, p2[i]));
208 m = this.F.mul(m, this.w[bits]);
209 }
210
211 return out;
212 }
213
214 extend(p, e) {
215 if (e == p.length) return p;
216 const z = new Array(e-p.length).fill(this.F.zero);
217
218 return p.concat(z);
219 }
220
221 reduce(p) {
222 if (p.length == 0) return p;
223 if (! this.F.isZero(p[p.length-1]) ) return p;
224 let i=p.length-1;
225 while( i>0 && this.F.isZero(p[i]) ) i--;
226 return p.slice(0, i+1);
227 }
228
229 affine(p) {
230 for (let i=0; i<p.length; i++) {
231 p[i] = this.F.affine(p[i]);
232 }
233 return p;
234 }
235
236 equals(a, b) {
237 const pa = this.reduce(this.affine(a));
238 const pb = this.reduce(this.affine(b));
239
240 if (pa.length != pb.length) return false;
241 for (let i=0; i<pb.length; i++) {
242 if (!this.F.equals(pa[i], pb[i])) return false;
243 }
244
245 return true;
246 }
247
248 ruffini(p, r) {
249 const res = new Array(p.length-1);
250 res[res.length-1] = p[p.length-1];
251 for (let i = res.length-2; i>=0; i--) {
252 res[i] = this.F.add(this.F.mul(res[i+1], r), p[i+1]);
253 }
254 return res;
255 }
256
257 _next2Power(v) {
258 v--;
259 v |= v >> 1;
260 v |= v >> 2;
261 v |= v >> 4;
262 v |= v >> 8;
263 v |= v >> 16;
264 v++;
265 return v;
266 }
267
268 toString(p) {
269 const ap = this.affine(p);
270 let S = "";
271 for (let i=ap.length-1; i>=0; i--) {
272 if (!this.F.isZero(p[i])) {
273 if (S!="") S += " + ";
274 S = S + p[i].toString(10);
275 if (i>0) {
276 S = S + "x";
277 if (i>1) {
278 S = S + "^" +i;
279 }
280 }
281 }
282 }
283 return S;
284 }
285
286
287 _reciprocal(p, bits) {
288 const k = 1 << bits;
289 if (k==1) {
290 return [ this.F.inverse(p[0]) ];
291 }
292 const np = this.scaleX(p, -k/2);
293 const q = this._reciprocal(np, bits-1);
294 const a = this.scaleX(this.double(q), 3*k/2-2);
295 const b = this.mul( this.square(q), p);
296
297 return this.scaleX(this.sub(a,b), -(k-2));
298 }
299
300 // divides x^m / v
301 _div2(m, v) {
302 const kbits = log2(v.length-1)+1;
303 const k = 1 << kbits;
304
305 const scaleV = k - v.length;
306
307 // rec = x^(k - 2) / v* x^scaleV =>
308 // rec = x^(k-2-scaleV)/ v
309 //
310 // res = x^m/v = x^(m + (2*k-2 - scaleV) - (2*k-2 - scaleV)) /v =>
311 // res = rec * x^(m - (2*k-2 - scaleV)) =>
312 // res = rec * x^(m - 2*k + 2 + scaleV)
313
314 const rec = this._reciprocal(this.scaleX(v, scaleV), kbits);
315 const res = this.scaleX(rec, m - 2*k + 2 + scaleV);
316
317 return res;
318 }
319
320 div(_u, _v) {
321 if (_u.length < _v.length) return [];
322 const kbits = log2(_v.length-1)+1;
323 const k = 1 << kbits;
324
325 const u = this.scaleX(_u, k-_v.length);
326 const v = this.scaleX(_v, k-_v.length);
327
328 const n = v.length-1;
329 let m = u.length-1;
330
331 const s = this._reciprocal(v, kbits);
332 let t;
333 if (m>2*n) {
334 t = this.sub(this.scaleX([this.F.one], 2*n), this.mul(s, v));
335 }
336
337 let q = [];
338 let rem = u;
339 let us, ut;
340 let finish = false;
341
342 while (!finish) {
343 us = this.mul(rem, s);
344 q = this.add(q, this.scaleX(us, -2*n));
345
346 if ( m > 2*n ) {
347 ut = this.mul(rem, t);
348 rem = this.scaleX(ut, -2*n);
349 m = rem.length-1;
350 } else {
351 finish = true;
352 }
353 }
354
355 return q;
356 }
357}
358
359function log2( V )
360{
361 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 ) );
362}
363
364module.exports = PolField;