1  (function() {
  2      var me = {};
  3  
  4      function ZeroDivisionError() {
  5          if (!this) return new ZeroDivisionError();
  6          this.message = "division by zero";
  7          this.name = "ZeroDivisionError";
  8      }
  9      me.ZeroDivisionError = ZeroDivisionError;
 10  
 11      // per-byte 2^8 Galois field
 12      // Note that this imposes a hard limit that the number of extended chunks can
 13      // be at most 256 along each dimension
 14      function galoistpl(a) {
 15          // 2 is not a primitive root, so we have to use 3 as our logarithm base
 16          var r = a ^ (a<<1); // a * (x+1)
 17          if (r > 0xff) { // overflow?
 18              r = r ^ 0x11b;
 19          }
 20          return r;
 21      }
 22  
 23      // Precomputing a multiplication and XOR table for increased speed
 24      var glogtable = new Array(256);
 25      var gexptable = [];
 26      (function() {
 27          var v = 1;
 28          for (var i = 0; i < 255; i++) {
 29              glogtable[v] = i;
 30              gexptable.push(v);
 31              v = galoistpl(v);
 32          }
 33      })();
 34      me.glogtable = glogtable;
 35      me.gexptable = gexptable;
 36  
 37      function Galois(val) {
 38          if (!(this instanceof Galois)) return new Galois(val);
 39          if (val instanceof Galois) {
 40              this.val = val.val;
 41          } else {
 42              this.val = val;
 43          }
 44          if (typeof Object.freeze == 'function') {
 45              Object.freeze(this);
 46          }
 47      }
 48      me.Galois = Galois;
 49      Galois.prototype.add = Galois.prototype.sub = function(other) {
 50          return new Galois(this.val ^ other.val);
 51      };
 52      Galois.prototype.mul = function(other) {
 53          if (this.val == 0 || other.val == 0) {
 54              return new Galois(0);
 55          }
 56          return new Galois(gexptable[(glogtable[this.val] +
 57                      glogtable[other.val]) % 255]);
 58      };
 59      Galois.prototype.div = function(other) {
 60          if (other.val == 0) {
 61              throw new ZeroDivisionError();
 62          }
 63          if (this.val == 0) {
 64              return new Galois(0);
 65          }
 66          return new Galois(gexptable[(glogtable[this.val] + 255 -
 67                      glogtable[other.val]) % 255]);
 68      };
 69      Galois.prototype.inspect = function() {
 70          return ""+this.val;
 71      };
 72  
 73      function powmod(b, e, m) {
 74          var r = 1;
 75          while (e > 0) {
 76              if (e & 1) r = (r * b) % m;
 77              b = (b * b) % m;
 78              e = e >> 1;
 79          }
 80          return r;
 81      }
 82  
 83  
 84      // Modular division class
 85      function mkModuloClass(n) {
 86          if (n <= 2) throw new Error("n must be prime!");
 87          for (var divisor = 2; divisor * divisor <= n; divisor++) {
 88              if (n % divisor == 0) {
 89                  throw new Error("n must be prime!");
 90              }
 91          }
 92  
 93          function Mod(val) {
 94              if (!(this instanceof Mod)) return new Mod(val);
 95              if (val instanceof Mod) {
 96                  this.val = val.val;
 97              } else {
 98                  this.val = val;
 99              }
100              if (typeof Object.freeze == 'function') {
101                  Object.freeze(this);
102              }
103          }
104          Mod.modulo = n;
105          Mod.prototype.add = function(other) {
106              return new Mod((this.val + other.val) % n);
107          };
108          Mod.prototype.sub = function(other) {
109              return new Mod((this.val + n - other.val) % n);
110          };
111          Mod.prototype.mul = function(other) {
112              return new Mod((this.val * other.val) % n);
113          };
114          Mod.prototype.div = function(other) {
115              return new Mod((this.val * powmod(other.val, n-2, n)) % n);
116          };
117          Mod.prototype.inspect = function() {
118              return ""+this.val;
119          };
120  
121          return Mod;
122      }
123      me.mkModuloClass = mkModuloClass;
124  
125      // Evaluates a polynomial in little-endian form, eg. x^2 + 3x + 2 = [2, 3, 1]
126      // (normally I hate little-endian, but in this case dealing with polynomials
127      // it's justified, since you get the nice property that p[n] is the nth degree
128      // term of p) at coordinate x, eg. eval_poly_at([2, 3, 1], 5) = 42 if you are
129      // using float as your arithmetic
130      function eval_poly_at(p, x) {
131          var arithmetic = p[0].constructor;
132          var y = new arithmetic(0);
133          var x_to_the_i = new arithmetic(1);
134          for (var i = 0; i < p.length; i++) {
135              y = y.add(x_to_the_i.mul(p[i]))
136                  x_to_the_i = x_to_the_i.mul(x);
137          }
138          return y;
139      }
140      me.eval_poly_at = eval_poly_at;
141  
142      // Given p+1 y values and x values with no errors, recovers the original
143      // p+1 degree polynomial. For example,
144      // lagrange_interp([51.0, 59.0, 66.0], [1, 3, 4]) = [50.0, 0, 1.0]
145      // if you are using float as your arithmetic
146      function lagrange_interp(pieces, xs) {
147          var arithmetic = pieces[0].constructor;
148          var zero = new arithmetic(0);
149          var one = new arithmetic(1);
150          // Generate master numerator polynomial
151          var root = [one];
152          var i, j;
153          for (i = 0; i < xs.length; i++) {
154              root.unshift(zero);
155              for (j = 0; j < root.length - 1; j++) {
156                  root[j] = root[j].sub(root[j+1].mul(xs[i]));
157              }
158          }
159          // Generate per-value numerator polynomials by dividing the master
160          // polynomial back by each x coordinate
161          var nums = [];
162          for (i = 0; i < xs.length; i++) {
163              var output = [];
164              var last = one;
165              for (j = 2; j < root.length+1; j++) {
166                  output.unshift(last);
167                  if (j != root.length) {
168                      last = root[root.length-j].add(last.mul(xs[i]));
169                  }
170              }
171              nums.push(output);
172          }
173          // Generate denominators by evaluating numerator polys at their x
174          var denoms = [];
175          for (i = 0; i < xs.length; i++) {
176              var denom = zero;
177              var x_to_the_j = one;
178              for (j = 0; j < nums[i].length; j++) {
179                  denom = denom.add(x_to_the_j.mul(nums[i][j]));
180                  x_to_the_j = x_to_the_j.mul(xs[i]);
181              }
182              denoms.push(denom);
183          }
184          // Generate output polynomial
185          var b = [];
186          for (i = 0; i < pieces.length; i++) {
187              b[i] = zero;
188          }
189          for (i = 0; i < xs.length; i++) {
190              var yslice = pieces[i].div(denoms[i]);
191              for (j = 0; j < pieces.length; j++) {
192                  b[j] = b[j].add(nums[i][j].mul(yslice));
193              }
194          }
195          return b;
196      }
197      me.lagrange_interp = lagrange_interp;
198  
199      // Compresses two linear equations of length n into one
200      // equation of length n-1
201      // Format:
202      // 3x + 4y = 80 (ie. 3x + 4y - 80 = 0) -> a = [3,4,-80]
203      // 5x + 2y = 70 (ie. 5x + 2y - 70 = 0) -> b = [5,2,-70]
204      function elim(a, b) {
205          var c = [];
206          for (var i = 1; i < a.length; i++) {
207              c[i-1] = a[i].mul(b[0]).sub(b[i].mul(a[0]));
208          }
209          return c;
210      }
211  
212      // Linear equation solver
213      // Format:
214      // 3x + 4y = 80, y = 5 (ie. 3x + 4y - 80z = 0, y = 5, z = 1)
215      //      -> coeffs = [3,4,-80], vals = [5,1]
216      function evaluate(coeffs, vals) {
217          var arithmetic = coeffs[0].constructor;
218          var tot = new arithmetic(0);
219          for (var i = 0; i < vals.length; i++) {
220              tot = tot.sub(coeffs[i+1].mul(vals[i]));
221          }
222          if (coeffs[0].val == 0) {
223              throw new ZeroDivisionError();
224          }
225          return tot.div(coeffs[0]);
226      }
227  
228      // Linear equation system solver
229      // Format:
230      // ax + by + c = 0, dx + ey + f = 0
231      // -> [[a, b, c], [d, e, f]]
232      // eg.
233      // [[3.0, 5.0, -13.0], [9.0, 1.0, -11.0]] -> [1.0, 2.0]
234      function sys_solve(eqs) {
235          var arithmetic = eqs[0][0].constructor;
236          var one = new arithmetic(1);
237          var back_eqs = [eqs[0]];
238          var i;
239          while (eqs.length > 1) {
240              var neweqs = [];
241              for (i = 0; i < eqs.length - 1; i++) {
242                  neweqs.push(elim(eqs[i], eqs[i+1]));
243              }
244              eqs = neweqs;
245              i = 0;
246              while (i < eqs.length - 1 && eqs[i][0].val == 0) {
247                  i++;
248              }
249              back_eqs.unshift(eqs[i]);
250          }
251          var kvals = [one];
252          for (i = 0; i < back_eqs.length; i++) {
253              kvals.unshift(evaluate(back_eqs[i], kvals));
254          }
255          return kvals.slice(0, -1);
256      }
257      me.sys_solve = sys_solve;
258  
259      function polydiv(Q, E) {
260          var qpoly = Q.slice();
261          var epoly = E.slice();
262          var div = [];
263          while (qpoly.length >= epoly.length) {
264              div.unshift(qpoly[qpoly.length-1].div(epoly[epoly.length-1]));
265              for (var i = 2; i < epoly.length + 1; i++) {
266                  qpoly[qpoly.length-i] =
267                      qpoly[qpoly.length-i].sub(div[0].mul(epoly[epoly.length-i]));
268              }
269              qpoly.pop();
270          }
271          return div;
272      }
273      me.polydiv = polydiv;
274  
275      // Given a set of y coordinates and x coordinates, and the degree of the
276      // original polynomial, determines the original polynomial even if some of
277      // the y coordinates are wrong. If m is the minimal number of pieces (ie.
278      // degree + 1), t is the total number of pieces provided, then the algo can
279      // handle up to (t-m)/2 errors. See:
280      // http://en.wikipedia.org/wiki/Berlekamp%E2%80%93Welch_algorithm#Example
281      // (just skip to my example, the rest of the article sucks imo)
282      function berlekamp_welch_attempt(pieces, xs, master_degree) {
283          var error_locator_degree = Math.floor((pieces.length - master_degree - 1) / 2);
284          var arithmetic = pieces[0].constructor;
285          var zero = new arithmetic(0);
286          var one = new arithmetic(1);
287          // Set up the equations for y[i]E(x[i]) = Q(x[i])
288          // degree(E) = error_locator_degree
289          // degree(Q) = master_degree + error_locator_degree - 1
290          var eqs = [];
291          var i,j;
292          for (i = 0; i < 2 * error_locator_degree + master_degree + 1; i++) {
293              eqs.push([]);
294          }
295          for (i = 0; i < 2 * error_locator_degree + master_degree + 1; i++) {
296              var neg_x_to_the_j = zero.sub(one);
297              for (j = 0; j < error_locator_degree + master_degree + 1; j++) {
298                  eqs[i].push(neg_x_to_the_j);
299                  neg_x_to_the_j = neg_x_to_the_j.mul(xs[i]);
300              }
301              var x_to_the_j = one;
302              for (j = 0; j < error_locator_degree + 1; j++) {
303                  eqs[i].push(x_to_the_j.mul(pieces[i]));
304                  x_to_the_j = x_to_the_j.mul(xs[i]);
305              }
306          }
307          // Solve 'em
308          // Assume the top error polynomial term to be one
309          var errors = error_locator_degree;
310          var ones = 1;
311          var polys;
312          while (errors >= 0) {
313              try {
314                  polys = sys_solve(eqs);
315                  for (i = 0; i < ones; i++) polys.push(one);
316                  break;
317              } catch (e) {
318                  if (e instanceof ZeroDivisionError) {
319                      eqs.pop();
320                      for (i = 0; i < eqs.length; i++) {
321                          var eq = eqs[i];
322                          eq[eq.length-2] = eq[eq.length-2].add(eq[eq.length-1]);
323                          eq.pop();
324                      }
325                      errors--;
326                      ones++;
327                  } else {
328                      throw e;
329                  }
330              }
331          }
332          if (errors < 0) {
333              throw new Error("Not enough data!");
334          }
335          // Divide the polynomials
336          var qpoly = polys.slice(0, error_locator_degree + master_degree + 1);
337          var epoly = polys.slice(error_locator_degree + master_degree + 1);
338          var div = polydiv(qpoly, epoly);
339          // Check
340          var corrects = 0;
341          for (i = 0; i < xs.length; i++) {
342              if (eval_poly_at(div, xs[i]).val == pieces[i].val) {
343                  corrects++;
344              }
345          }
346          if (corrects < master_degree + errors) {
347              throw new Error("Answer doesn't match (too many errors)!");
348          }
349          return div;
350      }
351      me.berlekamp_welch_attempt = berlekamp_welch_attempt;
352  
353      // Extends a list of integers in [0 ... 255] (if using Galois arithmetic) by
354      // adding n redundant error-correction values
355      function extend(data, n, arithmetic) {
356          arithmetic = arithmetic || Galois;
357          function mk(x) { return new arithmetic(x); }
358          var data2 = data.map(mk);
359          var data3 = data.slice();
360          var xs = [];
361          var i;
362          for (i = 0; i < data.length; i++) {
363              xs.push(new arithmetic(i));
364          }
365          var poly = berlekamp_welch_attempt(data2, xs, data.length - 1);
366          for (i = 0; i < n; i++) {
367              data3.push(eval_poly_at(poly, new arithmetic(data.length + i)).val);
368          }
369          return data3;
370      }
371      me.extend = extend;
372  
373      // Repairs a list of integers in [0 ... 255]. Some integers can be
374      // erroneous, and you can put null (or undefined) in place of an integer if
375      // you know that a certain value is defective or missing. Uses the
376      // Berlekamp-Welch algorithm to do error-correction
377      function repair(data, datasize, arithmetic) {
378          arithmetic = arithmetic || Galois;
379          var vs = [];
380          var xs = [];
381          var i;
382          for (var i = 0; i < data.length; i++) {
383              if (data[i] != null) {
384                  vs.push(new arithmetic(data[i]));
385                  xs.push(new arithmetic(i));
386              }
387          }
388          var poly = berlekamp_welch_attempt(vs, xs, datasize - 1);
389          var result = [];
390          for (i = 0; i < data.length; i++) {
391              result.push(eval_poly_at(poly, new arithmetic(i)).val);
392          }
393          return result;
394      }
395      me.repair = repair;
396  
397      function transpose(xs) {
398          var ys = [];
399          for (var i = 0; i < xs[0].length; i++) {
400              var y = [];
401              for (var j = 0; j < xs.length; j++) {
402                  y.push(xs[j][i]);
403              }
404              ys.push(y);
405          }
406          return ys;
407      }
408  
409      // Extends a list of bytearrays
410      // eg. extend_chunks([map(ord, 'hello'), map(ord, 'world')], 2)
411      // n is the number of redundant error-correction chunks to add
412      function extend_chunks(data, n, arithmetic) {
413          arithmetic = arithmetic || Galois;
414          var o = [];
415          for (var i = 0; i < data[0].length; i++) {
416              o.push(extend(data.map(function(x) { return x[i]; }), n, arithmetic));
417          }
418          return transpose(o);
419      }
420      me.extend_chunks = extend_chunks;
421  
422      // Repairs a list of bytearrays. Use null in place of a missing array.
423      // Individual arrays can contain some missing or erroneous data.
424      function repair_chunks(data, datasize, arithmetic) {
425          arithmetic = arithmetic || Galois;
426          var first_nonzero = 0;
427          while (data[first_nonzero] == null) {
428              first_nonzero++;
429          }
430          var i;
431          for (i = 0; i < data.length; i++) {
432              if (data[i] == null) {
433                  data[i] = new Array(data[first_nonzero].length);
434              }
435          }
436          var o = [];
437          for (i = 0; i < data[0].length; i++) {
438              o.push(repair(data.map(function(x) { return x[i]; }), datasize, arithmetic));
439          }
440          return transpose(o);
441      }
442      me.repair_chunks = repair_chunks;
443  
444      // Extends either a bytearray or a list of bytearrays or a list of lists...
445      // Used in the cubify method to expand a cube in all dimensions
446      function deep_extend_chunks(data, n, arithmetic) {
447          arithmetic = arithmetic || Galois;
448          if (!(data[0] instanceof Array)) {
449              return extend(data, n, arithmetic)
450          } else {
451              var o = [];
452              for (var i = 0; i < data[0].length; i++) {
453                  o.push(deep_extend_chunks(
454                              data.map(function(x) { return x[i]; }), n, arithmetic));
455              }
456              return transpose(o);
457          }
458      }
459      me.deep_extend_chunks = deep_extend_chunks;
460  
461      function isObject(o) {
462          return typeof o == 'object' || typeof o == 'function';
463      }
464      if (typeof define == 'function' && typeof define.amd == 'object' && define.amd) {
465          define(function() {
466              return me;
467          });
468      } else {
469          done = 0
470          try {
471              if (isObject(module)) { module.exports = me; }
472              else (isObject(window) ? window : this).Erasure = me;
473          }
474          catch(e) {
475              (isObject(window) ? window : this).Erasure = me;
476          }
477      }
478  }.call(this));