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