numeric.js

"use strict";

var numeric = typeof exports === "undefined" ? function numeric() {} : exports;
if (typeof global !== "undefined") {
  global.numeric = numeric;
}

numeric.version = "1.2.6";

// 1. Utility functions
numeric.bench = function bench(f, interval) {
  var t1, t2, n, i;
  if (typeof interval === "undefined") {
    interval = 15;
  }
  n = 0.5;
  t1 = new Date();
  while (1) {
    n *= 2;
    for (i = n; i > 3; i -= 4) {
      f();
      f();
      f();
      f();
    }
    while (i > 0) {
      f();
      i--;
    }
    t2 = new Date();
    if (t2 - t1 > interval) break;
  }
  for (i = n; i > 3; i -= 4) {
    f();
    f();
    f();
    f();
  }
  while (i > 0) {
    f();
    i--;
  }
  t2 = new Date();
  return (1000 * (3 * n - 1)) / (t2 - t1);
};

numeric._myIndexOf = function _myIndexOf(w) {
  var n = this.length,
    k;
  for (k = 0; k < n; ++k) if (this[k] === w) return k;
  return -1;
};
numeric.myIndexOf = Array.prototype.indexOf
  ? Array.prototype.indexOf
  : numeric._myIndexOf;

numeric.Function = Function;
numeric.precision = 4;
numeric.largeArray = 50;

numeric.prettyPrint = function prettyPrint(x) {
  function fmtnum(x) {
    if (x === 0) {
      return "0";
    }
    if (isNaN(x)) {
      return "NaN";
    }
    if (x < 0) {
      return "-" + fmtnum(-x);
    }
    if (isFinite(x)) {
      var scale = Math.floor(Math.log(x) / Math.log(10));
      var normalized = x / Math.pow(10, scale);
      var basic = normalized.toPrecision(numeric.precision);
      if (parseFloat(basic) === 10) {
        scale++;
        normalized = 1;
        basic = normalized.toPrecision(numeric.precision);
      }
      return parseFloat(basic).toString() + "e" + scale.toString();
    }
    return "Infinity";
  }
  var ret = [];
  function foo(x) {
    var k;
    if (typeof x === "undefined") {
      ret.push(Array(numeric.precision + 8).join(" "));
      return false;
    }
    if (typeof x === "string") {
      ret.push('"' + x + '"');
      return false;
    }
    if (typeof x === "boolean") {
      ret.push(x.toString());
      return false;
    }
    if (typeof x === "number") {
      var a = fmtnum(x);
      var b = x.toPrecision(numeric.precision);
      var c = parseFloat(x.toString()).toString();
      var d = [a, b, c, parseFloat(b).toString(), parseFloat(c).toString()];
      for (k = 1; k < d.length; k++) {
        if (d[k].length < a.length) a = d[k];
      }
      ret.push(Array(numeric.precision + 8 - a.length).join(" ") + a);
      return false;
    }
    if (x === null) {
      ret.push("null");
      return false;
    }
    if (typeof x === "function") {
      ret.push(x.toString());
      var flag = false;
      for (k in x) {
        if (x.hasOwnProperty(k)) {
          if (flag) ret.push(",\n");
          else ret.push("\n{");
          flag = true;
          ret.push(k);
          ret.push(": \n");
          foo(x[k]);
        }
      }
      if (flag) ret.push("}\n");
      return true;
    }
    if (x instanceof Array) {
      if (x.length > numeric.largeArray) {
        ret.push("...Large Array...");
        return true;
      }
      var flag = false;
      ret.push("[");
      for (k = 0; k < x.length; k++) {
        if (k > 0) {
          ret.push(",");
          if (flag) ret.push("\n ");
        }
        flag = foo(x[k]);
      }
      ret.push("]");
      return true;
    }
    ret.push("{");
    var flag = false;
    for (k in x) {
      if (x.hasOwnProperty(k)) {
        if (flag) ret.push(",\n");
        flag = true;
        ret.push(k);
        ret.push(": \n");
        foo(x[k]);
      }
    }
    ret.push("}");
    return true;
  }
  foo(x);
  return ret.join("");
};

numeric.parseDate = function parseDate(d) {
  function foo(d) {
    if (typeof d === "string") {
      return Date.parse(d.replace(/-/g, "/"));
    }
    if (!(d instanceof Array)) {
      throw new Error("parseDate: parameter must be arrays of strings");
    }
    var ret = [],
      k;
    for (k = 0; k < d.length; k++) {
      ret[k] = foo(d[k]);
    }
    return ret;
  }
  return foo(d);
};

numeric.parseFloat = function parseFloat_(d) {
  function foo(d) {
    if (typeof d === "string") {
      return parseFloat(d);
    }
    if (!(d instanceof Array)) {
      throw new Error("parseFloat: parameter must be arrays of strings");
    }
    var ret = [],
      k;
    for (k = 0; k < d.length; k++) {
      ret[k] = foo(d[k]);
    }
    return ret;
  }
  return foo(d);
};

numeric.parseCSV = function parseCSV(t) {
  var foo = t.split("\n");
  var j, k;
  var ret = [];
  var pat = /(([^'",]*)|('[^']*')|("[^"]*")),/g;
  var patnum = /^\s*(([+-]?[0-9]+(\.[0-9]*)?(e[+-]?[0-9]+)?)|([+-]?[0-9]*(\.[0-9]+)?(e[+-]?[0-9]+)?))\s*$/;
  var stripper = function(n) {
    return n.substr(0, n.length - 1);
  };
  var count = 0;
  for (k = 0; k < foo.length; k++) {
    var bar = (foo[k] + ",").match(pat),
      baz;
    if (bar.length > 0) {
      ret[count] = [];
      for (j = 0; j < bar.length; j++) {
        baz = stripper(bar[j]);
        if (patnum.test(baz)) {
          ret[count][j] = parseFloat(baz);
        } else ret[count][j] = baz;
      }
      count++;
    }
  }
  return ret;
};

numeric.toCSV = function toCSV(A) {
  var s = numeric.dim(A);
  var i, j, m, n, row, ret;
  m = s[0];
  n = s[1];
  ret = [];
  for (i = 0; i < m; i++) {
    row = [];
    for (j = 0; j < m; j++) {
      row[j] = A[i][j].toString();
    }
    ret[i] = row.join(", ");
  }
  return ret.join("\n") + "\n";
};

numeric.getURL = function getURL(url) {
  var client = new XMLHttpRequest();
  client.open("GET", url, false);
  client.send();
  return client;
};

numeric.imageURL = function imageURL(img) {
  function base64(A) {
    var n = A.length,
      i,
      x,
      y,
      z,
      p,
      q,
      r,
      s;
    var key =
      "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=";
    var ret = "";
    for (i = 0; i < n; i += 3) {
      x = A[i];
      y = A[i + 1];
      z = A[i + 2];
      p = x >> 2;
      q = ((x & 3) << 4) + (y >> 4);
      r = ((y & 15) << 2) + (z >> 6);
      s = z & 63;
      if (i + 1 >= n) {
        r = s = 64;
      } else if (i + 2 >= n) {
        s = 64;
      }
      ret += key.charAt(p) + key.charAt(q) + key.charAt(r) + key.charAt(s);
    }
    return ret;
  }
  function crc32Array(a, from, to) {
    if (typeof from === "undefined") {
      from = 0;
    }
    if (typeof to === "undefined") {
      to = a.length;
    }
    var table = [
      0x00000000,
      0x77073096,
      0xee0e612c,
      0x990951ba,
      0x076dc419,
      0x706af48f,
      0xe963a535,
      0x9e6495a3,
      0x0edb8832,
      0x79dcb8a4,
      0xe0d5e91e,
      0x97d2d988,
      0x09b64c2b,
      0x7eb17cbd,
      0xe7b82d07,
      0x90bf1d91,
      0x1db71064,
      0x6ab020f2,
      0xf3b97148,
      0x84be41de,
      0x1adad47d,
      0x6ddde4eb,
      0xf4d4b551,
      0x83d385c7,
      0x136c9856,
      0x646ba8c0,
      0xfd62f97a,
      0x8a65c9ec,
      0x14015c4f,
      0x63066cd9,
      0xfa0f3d63,
      0x8d080df5,
      0x3b6e20c8,
      0x4c69105e,
      0xd56041e4,
      0xa2677172,
      0x3c03e4d1,
      0x4b04d447,
      0xd20d85fd,
      0xa50ab56b,
      0x35b5a8fa,
      0x42b2986c,
      0xdbbbc9d6,
      0xacbcf940,
      0x32d86ce3,
      0x45df5c75,
      0xdcd60dcf,
      0xabd13d59,
      0x26d930ac,
      0x51de003a,
      0xc8d75180,
      0xbfd06116,
      0x21b4f4b5,
      0x56b3c423,
      0xcfba9599,
      0xb8bda50f,
      0x2802b89e,
      0x5f058808,
      0xc60cd9b2,
      0xb10be924,
      0x2f6f7c87,
      0x58684c11,
      0xc1611dab,
      0xb6662d3d,
      0x76dc4190,
      0x01db7106,
      0x98d220bc,
      0xefd5102a,
      0x71b18589,
      0x06b6b51f,
      0x9fbfe4a5,
      0xe8b8d433,
      0x7807c9a2,
      0x0f00f934,
      0x9609a88e,
      0xe10e9818,
      0x7f6a0dbb,
      0x086d3d2d,
      0x91646c97,
      0xe6635c01,
      0x6b6b51f4,
      0x1c6c6162,
      0x856530d8,
      0xf262004e,
      0x6c0695ed,
      0x1b01a57b,
      0x8208f4c1,
      0xf50fc457,
      0x65b0d9c6,
      0x12b7e950,
      0x8bbeb8ea,
      0xfcb9887c,
      0x62dd1ddf,
      0x15da2d49,
      0x8cd37cf3,
      0xfbd44c65,
      0x4db26158,
      0x3ab551ce,
      0xa3bc0074,
      0xd4bb30e2,
      0x4adfa541,
      0x3dd895d7,
      0xa4d1c46d,
      0xd3d6f4fb,
      0x4369e96a,
      0x346ed9fc,
      0xad678846,
      0xda60b8d0,
      0x44042d73,
      0x33031de5,
      0xaa0a4c5f,
      0xdd0d7cc9,
      0x5005713c,
      0x270241aa,
      0xbe0b1010,
      0xc90c2086,
      0x5768b525,
      0x206f85b3,
      0xb966d409,
      0xce61e49f,
      0x5edef90e,
      0x29d9c998,
      0xb0d09822,
      0xc7d7a8b4,
      0x59b33d17,
      0x2eb40d81,
      0xb7bd5c3b,
      0xc0ba6cad,
      0xedb88320,
      0x9abfb3b6,
      0x03b6e20c,
      0x74b1d29a,
      0xead54739,
      0x9dd277af,
      0x04db2615,
      0x73dc1683,
      0xe3630b12,
      0x94643b84,
      0x0d6d6a3e,
      0x7a6a5aa8,
      0xe40ecf0b,
      0x9309ff9d,
      0x0a00ae27,
      0x7d079eb1,
      0xf00f9344,
      0x8708a3d2,
      0x1e01f268,
      0x6906c2fe,
      0xf762575d,
      0x806567cb,
      0x196c3671,
      0x6e6b06e7,
      0xfed41b76,
      0x89d32be0,
      0x10da7a5a,
      0x67dd4acc,
      0xf9b9df6f,
      0x8ebeeff9,
      0x17b7be43,
      0x60b08ed5,
      0xd6d6a3e8,
      0xa1d1937e,
      0x38d8c2c4,
      0x4fdff252,
      0xd1bb67f1,
      0xa6bc5767,
      0x3fb506dd,
      0x48b2364b,
      0xd80d2bda,
      0xaf0a1b4c,
      0x36034af6,
      0x41047a60,
      0xdf60efc3,
      0xa867df55,
      0x316e8eef,
      0x4669be79,
      0xcb61b38c,
      0xbc66831a,
      0x256fd2a0,
      0x5268e236,
      0xcc0c7795,
      0xbb0b4703,
      0x220216b9,
      0x5505262f,
      0xc5ba3bbe,
      0xb2bd0b28,
      0x2bb45a92,
      0x5cb36a04,
      0xc2d7ffa7,
      0xb5d0cf31,
      0x2cd99e8b,
      0x5bdeae1d,
      0x9b64c2b0,
      0xec63f226,
      0x756aa39c,
      0x026d930a,
      0x9c0906a9,
      0xeb0e363f,
      0x72076785,
      0x05005713,
      0x95bf4a82,
      0xe2b87a14,
      0x7bb12bae,
      0x0cb61b38,
      0x92d28e9b,
      0xe5d5be0d,
      0x7cdcefb7,
      0x0bdbdf21,
      0x86d3d2d4,
      0xf1d4e242,
      0x68ddb3f8,
      0x1fda836e,
      0x81be16cd,
      0xf6b9265b,
      0x6fb077e1,
      0x18b74777,
      0x88085ae6,
      0xff0f6a70,
      0x66063bca,
      0x11010b5c,
      0x8f659eff,
      0xf862ae69,
      0x616bffd3,
      0x166ccf45,
      0xa00ae278,
      0xd70dd2ee,
      0x4e048354,
      0x3903b3c2,
      0xa7672661,
      0xd06016f7,
      0x4969474d,
      0x3e6e77db,
      0xaed16a4a,
      0xd9d65adc,
      0x40df0b66,
      0x37d83bf0,
      0xa9bcae53,
      0xdebb9ec5,
      0x47b2cf7f,
      0x30b5ffe9,
      0xbdbdf21c,
      0xcabac28a,
      0x53b39330,
      0x24b4a3a6,
      0xbad03605,
      0xcdd70693,
      0x54de5729,
      0x23d967bf,
      0xb3667a2e,
      0xc4614ab8,
      0x5d681b02,
      0x2a6f2b94,
      0xb40bbe37,
      0xc30c8ea1,
      0x5a05df1b,
      0x2d02ef8d
    ];

    var crc = -1,
      y = 0,
      n = a.length,
      i;

    for (i = from; i < to; i++) {
      y = (crc ^ a[i]) & 0xff;
      crc = (crc >>> 8) ^ table[y];
    }

    return crc ^ -1;
  }

  var h = img[0].length,
    w = img[0][0].length,
    s1,
    s2,
    next,
    k,
    length,
    a,
    b,
    i,
    j,
    adler32,
    crc32;
  var stream = [
    137,
    80,
    78,
    71,
    13,
    10,
    26,
    10, //  0: PNG signature
    0,
    0,
    0,
    13, //  8: IHDR Chunk length
    73,
    72,
    68,
    82, // 12: "IHDR"
    (w >> 24) & 255,
    (w >> 16) & 255,
    (w >> 8) & 255,
    w & 255, // 16: Width
    (h >> 24) & 255,
    (h >> 16) & 255,
    (h >> 8) & 255,
    h & 255, // 20: Height
    8, // 24: bit depth
    2, // 25: RGB
    0, // 26: deflate
    0, // 27: no filter
    0, // 28: no interlace
    -1,
    -2,
    -3,
    -4, // 29: CRC
    -5,
    -6,
    -7,
    -8, // 33: IDAT Chunk length
    73,
    68,
    65,
    84, // 37: "IDAT"
    // RFC 1950 header starts here
    8, // 41: RFC1950 CMF
    29 // 42: RFC1950 FLG
  ];
  crc32 = crc32Array(stream, 12, 29);
  stream[29] = (crc32 >> 24) & 255;
  stream[30] = (crc32 >> 16) & 255;
  stream[31] = (crc32 >> 8) & 255;
  stream[32] = crc32 & 255;
  s1 = 1;
  s2 = 0;
  for (i = 0; i < h; i++) {
    if (i < h - 1) {
      stream.push(0);
    } else {
      stream.push(1);
    }
    a = (3 * w + 1 + (i === 0)) & 255;
    b = ((3 * w + 1 + (i === 0)) >> 8) & 255;
    stream.push(a);
    stream.push(b);
    stream.push(~a & 255);
    stream.push(~b & 255);
    if (i === 0) stream.push(0);
    for (j = 0; j < w; j++) {
      for (k = 0; k < 3; k++) {
        a = img[k][i][j];
        if (a > 255) a = 255;
        else if (a < 0) a = 0;
        else a = Math.round(a);
        s1 = (s1 + a) % 65521;
        s2 = (s2 + s1) % 65521;
        stream.push(a);
      }
    }
    stream.push(0);
  }
  adler32 = (s2 << 16) + s1;
  stream.push((adler32 >> 24) & 255);
  stream.push((adler32 >> 16) & 255);
  stream.push((adler32 >> 8) & 255);
  stream.push(adler32 & 255);
  length = stream.length - 41;
  stream[33] = (length >> 24) & 255;
  stream[34] = (length >> 16) & 255;
  stream[35] = (length >> 8) & 255;
  stream[36] = length & 255;
  crc32 = crc32Array(stream, 37);
  stream.push((crc32 >> 24) & 255);
  stream.push((crc32 >> 16) & 255);
  stream.push((crc32 >> 8) & 255);
  stream.push(crc32 & 255);
  stream.push(0);
  stream.push(0);
  stream.push(0);
  stream.push(0);
  //    a = stream.length;
  stream.push(73); // I
  stream.push(69); // E
  stream.push(78); // N
  stream.push(68); // D
  stream.push(174); // CRC1
  stream.push(66); // CRC2
  stream.push(96); // CRC3
  stream.push(130); // CRC4
  return "data:image/png;base64," + base64(stream);
};

// 2. Linear algebra with Arrays.
numeric._dim = function _dim(x) {
  var ret = [];
  while (typeof x === "object") {
    ret.push(x.length);
    x = x[0];
  }
  return ret;
};

numeric.dim = function dim(x) {
  var y, z;
  if (typeof x === "object") {
    y = x[0];
    if (typeof y === "object") {
      z = y[0];
      if (typeof z === "object") {
        return numeric._dim(x);
      }
      return [x.length, y.length];
    }
    return [x.length];
  }
  return [];
};

numeric.mapreduce = function mapreduce(body, init) {
  return Function(
    "x",
    "accum",
    "_s",
    "_k",
    'if(typeof accum === "undefined") accum = ' +
      init +
      ";\n" +
      'if(typeof x === "number") { var xi = x; ' +
      body +
      "; return accum; }\n" +
      'if(typeof _s === "undefined") _s = numeric.dim(x);\n' +
      'if(typeof _k === "undefined") _k = 0;\n' +
      "var _n = _s[_k];\n" +
      "var i,xi;\n" +
      "if(_k < _s.length-1) {\n" +
      "    for(i=_n-1;i>=0;i--) {\n" +
      "        accum = arguments.callee(x[i],accum,_s,_k+1);\n" +
      "    }" +
      "    return accum;\n" +
      "}\n" +
      "for(i=_n-1;i>=1;i-=2) { \n" +
      "    xi = x[i];\n" +
      "    " +
      body +
      ";\n" +
      "    xi = x[i-1];\n" +
      "    " +
      body +
      ";\n" +
      "}\n" +
      "if(i === 0) {\n" +
      "    xi = x[i];\n" +
      "    " +
      body +
      "\n" +
      "}\n" +
      "return accum;"
  );
};
numeric.mapreduce2 = function mapreduce2(body, setup) {
  return Function(
    "x",
    "var n = x.length;\n" +
      "var i,xi;\n" +
      setup +
      ";\n" +
      "for(i=n-1;i!==-1;--i) { \n" +
      "    xi = x[i];\n" +
      "    " +
      body +
      ";\n" +
      "}\n" +
      "return accum;"
  );
};

numeric.same = function same(x, y) {
  var i, n;
  if (!(x instanceof Array) || !(y instanceof Array)) {
    return false;
  }
  n = x.length;
  if (n !== y.length) {
    return false;
  }
  for (i = 0; i < n; i++) {
    if (x[i] === y[i]) {
      continue;
    }
    if (typeof x[i] === "object") {
      if (!same(x[i], y[i])) return false;
    } else {
      return false;
    }
  }
  return true;
};

numeric.rep = function rep(s, v, k) {
  if (typeof k === "undefined") {
    k = 0;
  }
  var n = s[k],
    ret = Array(n),
    i;
  if (k === s.length - 1) {
    for (i = n - 2; i >= 0; i -= 2) {
      ret[i + 1] = v;
      ret[i] = v;
    }
    if (i === -1) {
      ret[0] = v;
    }
    return ret;
  }
  for (i = n - 1; i >= 0; i--) {
    ret[i] = numeric.rep(s, v, k + 1);
  }
  return ret;
};

numeric.dotMMsmall = function dotMMsmall(x, y) {
  var i, j, k, p, q, r, ret, foo, bar, woo, i0, k0, p0, r0;
  p = x.length;
  q = y.length;
  r = y[0].length;
  ret = Array(p);
  for (i = p - 1; i >= 0; i--) {
    foo = Array(r);
    bar = x[i];
    for (k = r - 1; k >= 0; k--) {
      woo = bar[q - 1] * y[q - 1][k];
      for (j = q - 2; j >= 1; j -= 2) {
        i0 = j - 1;
        woo += bar[j] * y[j][k] + bar[i0] * y[i0][k];
      }
      if (j === 0) {
        woo += bar[0] * y[0][k];
      }
      foo[k] = woo;
    }
    ret[i] = foo;
  }
  return ret;
};
numeric._getCol = function _getCol(A, j, x) {
  var n = A.length,
    i;
  for (i = n - 1; i > 0; --i) {
    x[i] = A[i][j];
    --i;
    x[i] = A[i][j];
  }
  if (i === 0) x[0] = A[0][j];
};
numeric.dotMMbig = function dotMMbig(x, y) {
  var gc = numeric._getCol,
    p = y.length,
    v = Array(p);
  var m = x.length,
    n = y[0].length,
    A = new Array(m),
    xj;
  var VV = numeric.dotVV;
  var i, j, k, z;
  --p;
  --m;
  for (i = m; i !== -1; --i) A[i] = Array(n);
  --n;
  for (i = n; i !== -1; --i) {
    gc(y, i, v);
    for (j = m; j !== -1; --j) {
      z = 0;
      xj = x[j];
      A[j][i] = VV(xj, v);
    }
  }
  return A;
};

numeric.dotMV = function dotMV(x, y) {
  var p = x.length,
    q = y.length,
    i;
  var ret = Array(p),
    dotVV = numeric.dotVV;
  for (i = p - 1; i >= 0; i--) {
    ret[i] = dotVV(x[i], y);
  }
  return ret;
};

numeric.dotVM = function dotVM(x, y) {
  var i,
    j,
    k,
    p,
    q,
    r,
    ret,
    foo,
    bar,
    woo,
    i0,
    k0,
    p0,
    r0,
    s1,
    s2,
    s3,
    baz,
    accum;
  p = x.length;
  q = y[0].length;
  ret = Array(q);
  for (k = q - 1; k >= 0; k--) {
    woo = x[p - 1] * y[p - 1][k];
    for (j = p - 2; j >= 1; j -= 2) {
      i0 = j - 1;
      woo += x[j] * y[j][k] + x[i0] * y[i0][k];
    }
    if (j === 0) {
      woo += x[0] * y[0][k];
    }
    ret[k] = woo;
  }
  return ret;
};

numeric.dotVV = function dotVV(x, y) {
  var i,
    n = x.length,
    i1,
    ret = x[n - 1] * y[n - 1];
  for (i = n - 2; i >= 1; i -= 2) {
    i1 = i - 1;
    ret += x[i] * y[i] + x[i1] * y[i1];
  }
  if (i === 0) {
    ret += x[0] * y[0];
  }
  return ret;
};

numeric.dot = function dot(x, y) {
  var d = numeric.dim;
  switch (d(x).length * 1000 + d(y).length) {
    case 2002:
      if (y.length < 10) return numeric.dotMMsmall(x, y);
      else return numeric.dotMMbig(x, y);
    case 2001:
      return numeric.dotMV(x, y);
    case 1002:
      return numeric.dotVM(x, y);
    case 1001:
      return numeric.dotVV(x, y);
    case 1000:
      return numeric.mulVS(x, y);
    case 1:
      return numeric.mulSV(x, y);
    case 0:
      return x * y;
    default:
      throw new Error("numeric.dot only works on vectors and matrices");
  }
};

numeric.diag = function diag(d) {
  var i,
    i1,
    j,
    n = d.length,
    A = Array(n),
    Ai;
  for (i = n - 1; i >= 0; i--) {
    Ai = Array(n);
    i1 = i + 2;
    for (j = n - 1; j >= i1; j -= 2) {
      Ai[j] = 0;
      Ai[j - 1] = 0;
    }
    if (j > i) {
      Ai[j] = 0;
    }
    Ai[i] = d[i];
    for (j = i - 1; j >= 1; j -= 2) {
      Ai[j] = 0;
      Ai[j - 1] = 0;
    }
    if (j === 0) {
      Ai[0] = 0;
    }
    A[i] = Ai;
  }
  return A;
};
numeric.getDiag = function(A) {
  var n = Math.min(A.length, A[0].length),
    i,
    ret = Array(n);
  for (i = n - 1; i >= 1; --i) {
    ret[i] = A[i][i];
    --i;
    ret[i] = A[i][i];
  }
  if (i === 0) {
    ret[0] = A[0][0];
  }
  return ret;
};

numeric.identity = function identity(n) {
  return numeric.diag(numeric.rep([n], 1));
};
numeric.pointwise = function pointwise(params, body, setup) {
  if (typeof setup === "undefined") {
    setup = "";
  }
  var fun = [];
  var k;
  var avec = /\[i\]$/,
    p,
    thevec = "";
  var haveret = false;
  for (k = 0; k < params.length; k++) {
    if (avec.test(params[k])) {
      p = params[k].substring(0, params[k].length - 3);
      thevec = p;
    } else {
      p = params[k];
    }
    if (p === "ret") haveret = true;
    fun.push(p);
  }
  fun[params.length] = "_s";
  fun[params.length + 1] = "_k";
  fun[params.length + 2] =
    'if(typeof _s === "undefined") _s = numeric.dim(' +
    thevec +
    ");\n" +
    'if(typeof _k === "undefined") _k = 0;\n' +
    "var _n = _s[_k];\n" +
    "var i" +
    (haveret ? "" : ", ret = Array(_n)") +
    ";\n" +
    "if(_k < _s.length-1) {\n" +
    "    for(i=_n-1;i>=0;i--) ret[i] = arguments.callee(" +
    params.join(",") +
    ",_s,_k+1);\n" +
    "    return ret;\n" +
    "}\n" +
    setup +
    "\n" +
    "for(i=_n-1;i!==-1;--i) {\n" +
    "    " +
    body +
    "\n" +
    "}\n" +
    "return ret;";
  return Function.apply(null, fun);
};
numeric.pointwise2 = function pointwise2(params, body, setup) {
  if (typeof setup === "undefined") {
    setup = "";
  }
  var fun = [];
  var k;
  var avec = /\[i\]$/,
    p,
    thevec = "";
  var haveret = false;
  for (k = 0; k < params.length; k++) {
    if (avec.test(params[k])) {
      p = params[k].substring(0, params[k].length - 3);
      thevec = p;
    } else {
      p = params[k];
    }
    if (p === "ret") haveret = true;
    fun.push(p);
  }
  fun[params.length] =
    "var _n = " +
    thevec +
    ".length;\n" +
    "var i" +
    (haveret ? "" : ", ret = Array(_n)") +
    ";\n" +
    setup +
    "\n" +
    "for(i=_n-1;i!==-1;--i) {\n" +
    body +
    "\n" +
    "}\n" +
    "return ret;";
  return Function.apply(null, fun);
};
numeric._biforeach = function _biforeach(x, y, s, k, f) {
  if (k === s.length - 1) {
    f(x, y);
    return;
  }
  var i,
    n = s[k];
  for (i = n - 1; i >= 0; i--) {
    _biforeach(
      typeof x === "object" ? x[i] : x,
      typeof y === "object" ? y[i] : y,
      s,
      k + 1,
      f
    );
  }
};
numeric._biforeach2 = function _biforeach2(x, y, s, k, f) {
  if (k === s.length - 1) {
    return f(x, y);
  }
  var i,
    n = s[k],
    ret = Array(n);
  for (i = n - 1; i >= 0; --i) {
    ret[i] = _biforeach2(
      typeof x === "object" ? x[i] : x,
      typeof y === "object" ? y[i] : y,
      s,
      k + 1,
      f
    );
  }
  return ret;
};
numeric._foreach = function _foreach(x, s, k, f) {
  if (k === s.length - 1) {
    f(x);
    return;
  }
  var i,
    n = s[k];
  for (i = n - 1; i >= 0; i--) {
    _foreach(x[i], s, k + 1, f);
  }
};
numeric._foreach2 = function _foreach2(x, s, k, f) {
  if (k === s.length - 1) {
    return f(x);
  }
  var i,
    n = s[k],
    ret = Array(n);
  for (i = n - 1; i >= 0; i--) {
    ret[i] = _foreach2(x[i], s, k + 1, f);
  }
  return ret;
};

/*numeric.anyV = numeric.mapreduce('if(xi) return true;','false');
numeric.allV = numeric.mapreduce('if(!xi) return false;','true');
numeric.any = function(x) { if(typeof x.length === "undefined") return x; return numeric.anyV(x); }
numeric.all = function(x) { if(typeof x.length === "undefined") return x; return numeric.allV(x); }*/

numeric.ops2 = {
  add: "+",
  sub: "-",
  mul: "*",
  div: "/",
  mod: "%",
  and: "&&",
  or: "||",
  eq: "===",
  neq: "!==",
  lt: "<",
  gt: ">",
  leq: "<=",
  geq: ">=",
  band: "&",
  bor: "|",
  bxor: "^",
  lshift: "<<",
  rshift: ">>",
  rrshift: ">>>"
};
numeric.opseq = {
  addeq: "+=",
  subeq: "-=",
  muleq: "*=",
  diveq: "/=",
  modeq: "%=",
  lshifteq: "<<=",
  rshifteq: ">>=",
  rrshifteq: ">>>=",
  bandeq: "&=",
  boreq: "|=",
  bxoreq: "^="
};
numeric.mathfuns = [
  "abs",
  "acos",
  "asin",
  "atan",
  "ceil",
  "cos",
  "exp",
  "floor",
  "log",
  "round",
  "sin",
  "sqrt",
  "tan",
  "isNaN",
  "isFinite"
];
numeric.mathfuns2 = ["atan2", "pow", "max", "min"];
numeric.ops1 = {
  neg: "-",
  not: "!",
  bnot: "~",
  clone: ""
};
numeric.mapreducers = {
  any: ["if(xi) return true;", "var accum = false;"],
  all: ["if(!xi) return false;", "var accum = true;"],
  sum: ["accum += xi;", "var accum = 0;"],
  prod: ["accum *= xi;", "var accum = 1;"],
  norm2Squared: ["accum += xi*xi;", "var accum = 0;"],
  norminf: [
    "accum = max(accum,abs(xi));",
    "var accum = 0, max = Math.max, abs = Math.abs;"
  ],
  norm1: ["accum += abs(xi)", "var accum = 0, abs = Math.abs;"],
  sup: ["accum = max(accum,xi);", "var accum = -Infinity, max = Math.max;"],
  inf: ["accum = min(accum,xi);", "var accum = Infinity, min = Math.min;"]
};

(function() {
  var i, o;
  for (i = 0; i < numeric.mathfuns2.length; ++i) {
    o = numeric.mathfuns2[i];
    numeric.ops2[o] = o;
  }
  for (i in numeric.ops2) {
    if (numeric.ops2.hasOwnProperty(i)) {
      o = numeric.ops2[i];
      var code,
        codeeq,
        setup = "";
      if (numeric.myIndexOf.call(numeric.mathfuns2, i) !== -1) {
        setup = "var " + o + " = Math." + o + ";\n";
        code = function(r, x, y) {
          return r + " = " + o + "(" + x + "," + y + ")";
        };
        codeeq = function(x, y) {
          return x + " = " + o + "(" + x + "," + y + ")";
        };
      } else {
        code = function(r, x, y) {
          return r + " = " + x + " " + o + " " + y;
        };
        if (numeric.opseq.hasOwnProperty(i + "eq")) {
          codeeq = function(x, y) {
            return x + " " + o + "= " + y;
          };
        } else {
          codeeq = function(x, y) {
            return x + " = " + x + " " + o + " " + y;
          };
        }
      }
      numeric[i + "VV"] = numeric.pointwise2(
        ["x[i]", "y[i]"],
        code("ret[i]", "x[i]", "y[i]"),
        setup
      );
      numeric[i + "SV"] = numeric.pointwise2(
        ["x", "y[i]"],
        code("ret[i]", "x", "y[i]"),
        setup
      );
      numeric[i + "VS"] = numeric.pointwise2(
        ["x[i]", "y"],
        code("ret[i]", "x[i]", "y"),
        setup
      );
      numeric[i] = Function(
        "var n = arguments.length, i, x = arguments[0], y;\n" +
          "var VV = numeric." +
          i +
          "VV, VS = numeric." +
          i +
          "VS, SV = numeric." +
          i +
          "SV;\n" +
          "var dim = numeric.dim;\n" +
          "for(i=1;i!==n;++i) { \n" +
          "  y = arguments[i];\n" +
          '  if(typeof x === "object") {\n' +
          '      if(typeof y === "object") x = numeric._biforeach2(x,y,dim(x),0,VV);\n' +
          "      else x = numeric._biforeach2(x,y,dim(x),0,VS);\n" +
          '  } else if(typeof y === "object") x = numeric._biforeach2(x,y,dim(y),0,SV);\n' +
          "  else " +
          codeeq("x", "y") +
          "\n" +
          "}\nreturn x;\n"
      );
      numeric[o] = numeric[i];
      numeric[i + "eqV"] = numeric.pointwise2(
        ["ret[i]", "x[i]"],
        codeeq("ret[i]", "x[i]"),
        setup
      );
      numeric[i + "eqS"] = numeric.pointwise2(
        ["ret[i]", "x"],
        codeeq("ret[i]", "x"),
        setup
      );
      numeric[i + "eq"] = Function(
        "var n = arguments.length, i, x = arguments[0], y;\n" +
          "var V = numeric." +
          i +
          "eqV, S = numeric." +
          i +
          "eqS\n" +
          "var s = numeric.dim(x);\n" +
          "for(i=1;i!==n;++i) { \n" +
          "  y = arguments[i];\n" +
          '  if(typeof y === "object") numeric._biforeach(x,y,s,0,V);\n' +
          "  else numeric._biforeach(x,y,s,0,S);\n" +
          "}\nreturn x;\n"
      );
    }
  }
  for (i = 0; i < numeric.mathfuns2.length; ++i) {
    o = numeric.mathfuns2[i];
    delete numeric.ops2[o];
  }
  for (i = 0; i < numeric.mathfuns.length; ++i) {
    o = numeric.mathfuns[i];
    numeric.ops1[o] = o;
  }
  for (i in numeric.ops1) {
    if (numeric.ops1.hasOwnProperty(i)) {
      setup = "";
      o = numeric.ops1[i];
      if (numeric.myIndexOf.call(numeric.mathfuns, i) !== -1) {
        if (Math.hasOwnProperty(o)) setup = "var " + o + " = Math." + o + ";\n";
      }
      numeric[i + "eqV"] = numeric.pointwise2(
        ["ret[i]"],
        "ret[i] = " + o + "(ret[i]);",
        setup
      );
      numeric[i + "eq"] = Function(
        "x",
        'if(typeof x !== "object") return ' +
          o +
          "x\n" +
          "var i;\n" +
          "var V = numeric." +
          i +
          "eqV;\n" +
          "var s = numeric.dim(x);\n" +
          "numeric._foreach(x,s,0,V);\n" +
          "return x;\n"
      );
      numeric[i + "V"] = numeric.pointwise2(
        ["x[i]"],
        "ret[i] = " + o + "(x[i]);",
        setup
      );
      numeric[i] = Function(
        "x",
        'if(typeof x !== "object") return ' +
          o +
          "(x)\n" +
          "var i;\n" +
          "var V = numeric." +
          i +
          "V;\n" +
          "var s = numeric.dim(x);\n" +
          "return numeric._foreach2(x,s,0,V);\n"
      );
    }
  }
  for (i = 0; i < numeric.mathfuns.length; ++i) {
    o = numeric.mathfuns[i];
    delete numeric.ops1[o];
  }
  for (i in numeric.mapreducers) {
    if (numeric.mapreducers.hasOwnProperty(i)) {
      o = numeric.mapreducers[i];
      numeric[i + "V"] = numeric.mapreduce2(o[0], o[1]);
      numeric[i] = Function(
        "x",
        "s",
        "k",
        o[1] +
          'if(typeof x !== "object") {' +
          "    xi = x;\n" +
          o[0] +
          ";\n" +
          "    return accum;\n" +
          "}" +
          'if(typeof s === "undefined") s = numeric.dim(x);\n' +
          'if(typeof k === "undefined") k = 0;\n' +
          "if(k === s.length-1) return numeric." +
          i +
          "V(x);\n" +
          "var xi;\n" +
          "var n = x.length, i;\n" +
          "for(i=n-1;i!==-1;--i) {\n" +
          "   xi = arguments.callee(x[i]);\n" +
          o[0] +
          ";\n" +
          "}\n" +
          "return accum;\n"
      );
    }
  }
})();

numeric.truncVV = numeric.pointwise(
  ["x[i]", "y[i]"],
  "ret[i] = round(x[i]/y[i])*y[i];",
  "var round = Math.round;"
);
numeric.truncVS = numeric.pointwise(
  ["x[i]", "y"],
  "ret[i] = round(x[i]/y)*y;",
  "var round = Math.round;"
);
numeric.truncSV = numeric.pointwise(
  ["x", "y[i]"],
  "ret[i] = round(x/y[i])*y[i];",
  "var round = Math.round;"
);
numeric.trunc = function trunc(x, y) {
  if (typeof x === "object") {
    if (typeof y === "object") return numeric.truncVV(x, y);
    return numeric.truncVS(x, y);
  }
  if (typeof y === "object") return numeric.truncSV(x, y);
  return Math.round(x / y) * y;
};

numeric.inv = function inv(x) {
  var s = numeric.dim(x),
    abs = Math.abs,
    m = s[0],
    n = s[1];
  var A = numeric.clone(x),
    Ai,
    Aj;
  var I = numeric.identity(m),
    Ii,
    Ij;
  var i, j, k, x;
  for (j = 0; j < n; ++j) {
    var i0 = -1;
    var v0 = -1;
    for (i = j; i !== m; ++i) {
      k = abs(A[i][j]);
      if (k > v0) {
        i0 = i;
        v0 = k;
      }
    }
    Aj = A[i0];
    A[i0] = A[j];
    A[j] = Aj;
    Ij = I[i0];
    I[i0] = I[j];
    I[j] = Ij;
    x = Aj[j];
    for (k = j; k !== n; ++k) Aj[k] /= x;
    for (k = n - 1; k !== -1; --k) Ij[k] /= x;
    for (i = m - 1; i !== -1; --i) {
      if (i !== j) {
        Ai = A[i];
        Ii = I[i];
        x = Ai[j];
        for (k = j + 1; k !== n; ++k) Ai[k] -= Aj[k] * x;
        for (k = n - 1; k > 0; --k) {
          Ii[k] -= Ij[k] * x;
          --k;
          Ii[k] -= Ij[k] * x;
        }
        if (k === 0) Ii[0] -= Ij[0] * x;
      }
    }
  }
  return I;
};

numeric.det = function det(x) {
  var s = numeric.dim(x);
  if (s.length !== 2 || s[0] !== s[1]) {
    throw new Error("numeric: det() only works on square matrices");
  }
  var n = s[0],
    ret = 1,
    i,
    j,
    k,
    A = numeric.clone(x),
    Aj,
    Ai,
    alpha,
    temp,
    k1,
    k2,
    k3;
  for (j = 0; j < n - 1; j++) {
    k = j;
    for (i = j + 1; i < n; i++) {
      if (Math.abs(A[i][j]) > Math.abs(A[k][j])) {
        k = i;
      }
    }
    if (k !== j) {
      temp = A[k];
      A[k] = A[j];
      A[j] = temp;
      ret *= -1;
    }
    Aj = A[j];
    for (i = j + 1; i < n; i++) {
      Ai = A[i];
      alpha = Ai[j] / Aj[j];
      for (k = j + 1; k < n - 1; k += 2) {
        k1 = k + 1;
        Ai[k] -= Aj[k] * alpha;
        Ai[k1] -= Aj[k1] * alpha;
      }
      if (k !== n) {
        Ai[k] -= Aj[k] * alpha;
      }
    }
    if (Aj[j] === 0) {
      return 0;
    }
    ret *= Aj[j];
  }
  return ret * A[j][j];
};

numeric.transpose = function transpose(x) {
  var i,
    j,
    m = x.length,
    n = x[0].length,
    ret = Array(n),
    A0,
    A1,
    Bj;
  for (j = 0; j < n; j++) ret[j] = Array(m);
  for (i = m - 1; i >= 1; i -= 2) {
    A1 = x[i];
    A0 = x[i - 1];
    for (j = n - 1; j >= 1; --j) {
      Bj = ret[j];
      Bj[i] = A1[j];
      Bj[i - 1] = A0[j];
      --j;
      Bj = ret[j];
      Bj[i] = A1[j];
      Bj[i - 1] = A0[j];
    }
    if (j === 0) {
      Bj = ret[0];
      Bj[i] = A1[0];
      Bj[i - 1] = A0[0];
    }
  }
  if (i === 0) {
    A0 = x[0];
    for (j = n - 1; j >= 1; --j) {
      ret[j][0] = A0[j];
      --j;
      ret[j][0] = A0[j];
    }
    if (j === 0) {
      ret[0][0] = A0[0];
    }
  }
  return ret;
};
numeric.negtranspose = function negtranspose(x) {
  var i,
    j,
    m = x.length,
    n = x[0].length,
    ret = Array(n),
    A0,
    A1,
    Bj;
  for (j = 0; j < n; j++) ret[j] = Array(m);
  for (i = m - 1; i >= 1; i -= 2) {
    A1 = x[i];
    A0 = x[i - 1];
    for (j = n - 1; j >= 1; --j) {
      Bj = ret[j];
      Bj[i] = -A1[j];
      Bj[i - 1] = -A0[j];
      --j;
      Bj = ret[j];
      Bj[i] = -A1[j];
      Bj[i - 1] = -A0[j];
    }
    if (j === 0) {
      Bj = ret[0];
      Bj[i] = -A1[0];
      Bj[i - 1] = -A0[0];
    }
  }
  if (i === 0) {
    A0 = x[0];
    for (j = n - 1; j >= 1; --j) {
      ret[j][0] = -A0[j];
      --j;
      ret[j][0] = -A0[j];
    }
    if (j === 0) {
      ret[0][0] = -A0[0];
    }
  }
  return ret;
};

numeric._random = function _random(s, k) {
  var i,
    n = s[k],
    ret = Array(n),
    rnd;
  if (k === s.length - 1) {
    rnd = Math.random;
    for (i = n - 1; i >= 1; i -= 2) {
      ret[i] = rnd();
      ret[i - 1] = rnd();
    }
    if (i === 0) {
      ret[0] = rnd();
    }
    return ret;
  }
  for (i = n - 1; i >= 0; i--) ret[i] = _random(s, k + 1);
  return ret;
};
numeric.random = function random(s) {
  return numeric._random(s, 0);
};

numeric.norm2 = function norm2(x) {
  return Math.sqrt(numeric.norm2Squared(x));
};

numeric.linspace = function linspace(a, b, n) {
  if (typeof n === "undefined") n = Math.max(Math.round(b - a) + 1, 1);
  if (n < 2) {
    return n === 1 ? [a] : [];
  }
  var i,
    ret = Array(n);
  n--;
  for (i = n; i >= 0; i--) {
    ret[i] = (i * b + (n - i) * a) / n;
  }
  return ret;
};

numeric.getBlock = function getBlock(x, from, to) {
  var s = numeric.dim(x);
  function foo(x, k) {
    var i,
      a = from[k],
      n = to[k] - a,
      ret = Array(n);
    if (k === s.length - 1) {
      for (i = n; i >= 0; i--) {
        ret[i] = x[i + a];
      }
      return ret;
    }
    for (i = n; i >= 0; i--) {
      ret[i] = foo(x[i + a], k + 1);
    }
    return ret;
  }
  return foo(x, 0);
};

numeric.setBlock = function setBlock(x, from, to, B) {
  var s = numeric.dim(x);
  function foo(x, y, k) {
    var i,
      a = from[k],
      n = to[k] - a;
    if (k === s.length - 1) {
      for (i = n; i >= 0; i--) {
        x[i + a] = y[i];
      }
    }
    for (i = n; i >= 0; i--) {
      foo(x[i + a], y[i], k + 1);
    }
  }
  foo(x, B, 0);
  return x;
};

numeric.getRange = function getRange(A, I, J) {
  var m = I.length,
    n = J.length;
  var i, j;
  var B = Array(m),
    Bi,
    AI;
  for (i = m - 1; i !== -1; --i) {
    B[i] = Array(n);
    Bi = B[i];
    AI = A[I[i]];
    for (j = n - 1; j !== -1; --j) Bi[j] = AI[J[j]];
  }
  return B;
};

numeric.blockMatrix = function blockMatrix(X) {
  var s = numeric.dim(X);
  if (s.length < 4) return numeric.blockMatrix([X]);
  var m = s[0],
    n = s[1],
    M,
    N,
    i,
    j,
    Xij;
  M = 0;
  N = 0;
  for (i = 0; i < m; ++i) M += X[i][0].length;
  for (j = 0; j < n; ++j) N += X[0][j][0].length;
  var Z = Array(M);
  for (i = 0; i < M; ++i) Z[i] = Array(N);
  var I = 0,
    J,
    ZI,
    k,
    l,
    Xijk;
  for (i = 0; i < m; ++i) {
    J = N;
    for (j = n - 1; j !== -1; --j) {
      Xij = X[i][j];
      J -= Xij[0].length;
      for (k = Xij.length - 1; k !== -1; --k) {
        Xijk = Xij[k];
        ZI = Z[I + k];
        for (l = Xijk.length - 1; l !== -1; --l) ZI[J + l] = Xijk[l];
      }
    }
    I += X[i][0].length;
  }
  return Z;
};

numeric.tensor = function tensor(x, y) {
  if (typeof x === "number" || typeof y === "number") return numeric.mul(x, y);
  var s1 = numeric.dim(x),
    s2 = numeric.dim(y);
  if (s1.length !== 1 || s2.length !== 1) {
    throw new Error("numeric: tensor product is only defined for vectors");
  }
  var m = s1[0],
    n = s2[0],
    A = Array(m),
    Ai,
    i,
    j,
    xi;
  for (i = m - 1; i >= 0; i--) {
    Ai = Array(n);
    xi = x[i];
    for (j = n - 1; j >= 3; --j) {
      Ai[j] = xi * y[j];
      --j;
      Ai[j] = xi * y[j];
      --j;
      Ai[j] = xi * y[j];
      --j;
      Ai[j] = xi * y[j];
    }
    while (j >= 0) {
      Ai[j] = xi * y[j];
      --j;
    }
    A[i] = Ai;
  }
  return A;
};

// 3. The Tensor type T
numeric.T = function T(x, y) {
  this.x = x;
  this.y = y;
};
numeric.t = function t(x, y) {
  return new numeric.T(x, y);
};

numeric.Tbinop = function Tbinop(rr, rc, cr, cc, setup) {
  var io = numeric.indexOf;
  if (typeof setup !== "string") {
    var k;
    setup = "";
    for (k in numeric) {
      if (
        numeric.hasOwnProperty(k) &&
        (rr.indexOf(k) >= 0 ||
          rc.indexOf(k) >= 0 ||
          cr.indexOf(k) >= 0 ||
          cc.indexOf(k) >= 0) &&
        k.length > 1
      ) {
        setup += "var " + k + " = numeric." + k + ";\n";
      }
    }
  }
  return Function(
    ["y"],
    "var x = this;\n" +
      "if(!(y instanceof numeric.T)) { y = new numeric.T(y); }\n" +
      setup +
      "\n" +
      "if(x.y) {" +
      "  if(y.y) {" +
      "    return new numeric.T(" +
      cc +
      ");\n" +
      "  }\n" +
      "  return new numeric.T(" +
      cr +
      ");\n" +
      "}\n" +
      "if(y.y) {\n" +
      "  return new numeric.T(" +
      rc +
      ");\n" +
      "}\n" +
      "return new numeric.T(" +
      rr +
      ");\n"
  );
};

numeric.T.prototype.add = numeric.Tbinop(
  "add(x.x,y.x)",
  "add(x.x,y.x),y.y",
  "add(x.x,y.x),x.y",
  "add(x.x,y.x),add(x.y,y.y)"
);
numeric.T.prototype.sub = numeric.Tbinop(
  "sub(x.x,y.x)",
  "sub(x.x,y.x),neg(y.y)",
  "sub(x.x,y.x),x.y",
  "sub(x.x,y.x),sub(x.y,y.y)"
);
numeric.T.prototype.mul = numeric.Tbinop(
  "mul(x.x,y.x)",
  "mul(x.x,y.x),mul(x.x,y.y)",
  "mul(x.x,y.x),mul(x.y,y.x)",
  "sub(mul(x.x,y.x),mul(x.y,y.y)),add(mul(x.x,y.y),mul(x.y,y.x))"
);

numeric.T.prototype.reciprocal = function reciprocal() {
  var mul = numeric.mul,
    div = numeric.div;
  if (this.y) {
    var d = numeric.add(mul(this.x, this.x), mul(this.y, this.y));
    return new numeric.T(div(this.x, d), div(numeric.neg(this.y), d));
  }
  return new T(div(1, this.x));
};
numeric.T.prototype.div = function div(y) {
  if (!(y instanceof numeric.T)) y = new numeric.T(y);
  if (y.y) {
    return this.mul(y.reciprocal());
  }
  var div = numeric.div;
  if (this.y) {
    return new numeric.T(div(this.x, y.x), div(this.y, y.x));
  }
  return new numeric.T(div(this.x, y.x));
};
numeric.T.prototype.dot = numeric.Tbinop(
  "dot(x.x,y.x)",
  "dot(x.x,y.x),dot(x.x,y.y)",
  "dot(x.x,y.x),dot(x.y,y.x)",
  "sub(dot(x.x,y.x),dot(x.y,y.y)),add(dot(x.x,y.y),dot(x.y,y.x))"
);
numeric.T.prototype.transpose = function transpose() {
  var t = numeric.transpose,
    x = this.x,
    y = this.y;
  if (y) {
    return new numeric.T(t(x), t(y));
  }
  return new numeric.T(t(x));
};
numeric.T.prototype.transjugate = function transjugate() {
  var t = numeric.transpose,
    x = this.x,
    y = this.y;
  if (y) {
    return new numeric.T(t(x), numeric.negtranspose(y));
  }
  return new numeric.T(t(x));
};
numeric.Tunop = function Tunop(r, c, s) {
  if (typeof s !== "string") {
    s = "";
  }
  return Function(
    "var x = this;\n" +
      s +
      "\n" +
      "if(x.y) {" +
      "  " +
      c +
      ";\n" +
      "}\n" +
      r +
      ";\n"
  );
};

numeric.T.prototype.exp = numeric.Tunop(
  "return new numeric.T(ex)",
  "return new numeric.T(mul(cos(x.y),ex),mul(sin(x.y),ex))",
  "var ex = numeric.exp(x.x), cos = numeric.cos, sin = numeric.sin, mul = numeric.mul;"
);
numeric.T.prototype.conj = numeric.Tunop(
  "return new numeric.T(x.x);",
  "return new numeric.T(x.x,numeric.neg(x.y));"
);
numeric.T.prototype.neg = numeric.Tunop(
  "return new numeric.T(neg(x.x));",
  "return new numeric.T(neg(x.x),neg(x.y));",
  "var neg = numeric.neg;"
);
numeric.T.prototype.sin = numeric.Tunop(
  "return new numeric.T(numeric.sin(x.x))",
  "return x.exp().sub(x.neg().exp()).div(new numeric.T(0,2));"
);
numeric.T.prototype.cos = numeric.Tunop(
  "return new numeric.T(numeric.cos(x.x))",
  "return x.exp().add(x.neg().exp()).div(2);"
);
numeric.T.prototype.abs = numeric.Tunop(
  "return new numeric.T(numeric.abs(x.x));",
  "return new numeric.T(numeric.sqrt(numeric.add(mul(x.x,x.x),mul(x.y,x.y))));",
  "var mul = numeric.mul;"
);
numeric.T.prototype.log = numeric.Tunop(
  "return new numeric.T(numeric.log(x.x));",
  "var theta = new numeric.T(numeric.atan2(x.y,x.x)), r = x.abs();\n" +
    "return new numeric.T(numeric.log(r.x),theta.x);"
);
numeric.T.prototype.norm2 = numeric.Tunop(
  "return numeric.norm2(x.x);",
  "var f = numeric.norm2Squared;\n" + "return Math.sqrt(f(x.x)+f(x.y));"
);
numeric.T.prototype.inv = function inv() {
  var A = this;
  if (typeof A.y === "undefined") {
    return new numeric.T(numeric.inv(A.x));
  }
  var n = A.x.length,
    i,
    j,
    k;
  var Rx = numeric.identity(n),
    Ry = numeric.rep([n, n], 0);
  var Ax = numeric.clone(A.x),
    Ay = numeric.clone(A.y);
  var Aix, Aiy, Ajx, Ajy, Rix, Riy, Rjx, Rjy;
  var i, j, k, d, d1, ax, ay, bx, by, temp;
  for (i = 0; i < n; i++) {
    ax = Ax[i][i];
    ay = Ay[i][i];
    d = ax * ax + ay * ay;
    k = i;
    for (j = i + 1; j < n; j++) {
      ax = Ax[j][i];
      ay = Ay[j][i];
      d1 = ax * ax + ay * ay;
      if (d1 > d) {
        k = j;
        d = d1;
      }
    }
    if (k !== i) {
      temp = Ax[i];
      Ax[i] = Ax[k];
      Ax[k] = temp;
      temp = Ay[i];
      Ay[i] = Ay[k];
      Ay[k] = temp;
      temp = Rx[i];
      Rx[i] = Rx[k];
      Rx[k] = temp;
      temp = Ry[i];
      Ry[i] = Ry[k];
      Ry[k] = temp;
    }
    Aix = Ax[i];
    Aiy = Ay[i];
    Rix = Rx[i];
    Riy = Ry[i];
    ax = Aix[i];
    ay = Aiy[i];
    for (j = i + 1; j < n; j++) {
      bx = Aix[j];
      by = Aiy[j];
      Aix[j] = (bx * ax + by * ay) / d;
      Aiy[j] = (by * ax - bx * ay) / d;
    }
    for (j = 0; j < n; j++) {
      bx = Rix[j];
      by = Riy[j];
      Rix[j] = (bx * ax + by * ay) / d;
      Riy[j] = (by * ax - bx * ay) / d;
    }
    for (j = i + 1; j < n; j++) {
      Ajx = Ax[j];
      Ajy = Ay[j];
      Rjx = Rx[j];
      Rjy = Ry[j];
      ax = Ajx[i];
      ay = Ajy[i];
      for (k = i + 1; k < n; k++) {
        bx = Aix[k];
        by = Aiy[k];
        Ajx[k] -= bx * ax - by * ay;
        Ajy[k] -= by * ax + bx * ay;
      }
      for (k = 0; k < n; k++) {
        bx = Rix[k];
        by = Riy[k];
        Rjx[k] -= bx * ax - by * ay;
        Rjy[k] -= by * ax + bx * ay;
      }
    }
  }
  for (i = n - 1; i > 0; i--) {
    Rix = Rx[i];
    Riy = Ry[i];
    for (j = i - 1; j >= 0; j--) {
      Rjx = Rx[j];
      Rjy = Ry[j];
      ax = Ax[j][i];
      ay = Ay[j][i];
      for (k = n - 1; k >= 0; k--) {
        bx = Rix[k];
        by = Riy[k];
        Rjx[k] -= ax * bx - ay * by;
        Rjy[k] -= ax * by + ay * bx;
      }
    }
  }
  return new numeric.T(Rx, Ry);
};
numeric.T.prototype.get = function get(i) {
  var x = this.x,
    y = this.y,
    k = 0,
    ik,
    n = i.length;
  if (y) {
    while (k < n) {
      ik = i[k];
      x = x[ik];
      y = y[ik];
      k++;
    }
    return new numeric.T(x, y);
  }
  while (k < n) {
    ik = i[k];
    x = x[ik];
    k++;
  }
  return new numeric.T(x);
};
numeric.T.prototype.set = function set(i, v) {
  var x = this.x,
    y = this.y,
    k = 0,
    ik,
    n = i.length,
    vx = v.x,
    vy = v.y;
  if (n === 0) {
    if (vy) {
      this.y = vy;
    } else if (y) {
      this.y = undefined;
    }
    this.x = x;
    return this;
  }
  if (vy) {
    if (y) {
      /* ok */
    } else {
      y = numeric.rep(numeric.dim(x), 0);
      this.y = y;
    }
    while (k < n - 1) {
      ik = i[k];
      x = x[ik];
      y = y[ik];
      k++;
    }
    ik = i[k];
    x[ik] = vx;
    y[ik] = vy;
    return this;
  }
  if (y) {
    while (k < n - 1) {
      ik = i[k];
      x = x[ik];
      y = y[ik];
      k++;
    }
    ik = i[k];
    x[ik] = vx;
    if (vx instanceof Array) y[ik] = numeric.rep(numeric.dim(vx), 0);
    else y[ik] = 0;
    return this;
  }
  while (k < n - 1) {
    ik = i[k];
    x = x[ik];
    k++;
  }
  ik = i[k];
  x[ik] = vx;
  return this;
};
numeric.T.prototype.getRows = function getRows(i0, i1) {
  var n = i1 - i0 + 1,
    j;
  var rx = Array(n),
    ry,
    x = this.x,
    y = this.y;
  for (j = i0; j <= i1; j++) {
    rx[j - i0] = x[j];
  }
  if (y) {
    ry = Array(n);
    for (j = i0; j <= i1; j++) {
      ry[j - i0] = y[j];
    }
    return new numeric.T(rx, ry);
  }
  return new numeric.T(rx);
};
numeric.T.prototype.setRows = function setRows(i0, i1, A) {
  var j;
  var rx = this.x,
    ry = this.y,
    x = A.x,
    y = A.y;
  for (j = i0; j <= i1; j++) {
    rx[j] = x[j - i0];
  }
  if (y) {
    if (!ry) {
      ry = numeric.rep(numeric.dim(rx), 0);
      this.y = ry;
    }
    for (j = i0; j <= i1; j++) {
      ry[j] = y[j - i0];
    }
  } else if (ry) {
    for (j = i0; j <= i1; j++) {
      ry[j] = numeric.rep([x[j - i0].length], 0);
    }
  }
  return this;
};
numeric.T.prototype.getRow = function getRow(k) {
  var x = this.x,
    y = this.y;
  if (y) {
    return new numeric.T(x[k], y[k]);
  }
  return new numeric.T(x[k]);
};
numeric.T.prototype.setRow = function setRow(i, v) {
  var rx = this.x,
    ry = this.y,
    x = v.x,
    y = v.y;
  rx[i] = x;
  if (y) {
    if (!ry) {
      ry = numeric.rep(numeric.dim(rx), 0);
      this.y = ry;
    }
    ry[i] = y;
  } else if (ry) {
    ry = numeric.rep([x.length], 0);
  }
  return this;
};

numeric.T.prototype.getBlock = function getBlock(from, to) {
  var x = this.x,
    y = this.y,
    b = numeric.getBlock;
  if (y) {
    return new numeric.T(b(x, from, to), b(y, from, to));
  }
  return new numeric.T(b(x, from, to));
};
numeric.T.prototype.setBlock = function setBlock(from, to, A) {
  if (!(A instanceof numeric.T)) A = new numeric.T(A);
  var x = this.x,
    y = this.y,
    b = numeric.setBlock,
    Ax = A.x,
    Ay = A.y;
  if (Ay) {
    if (!y) {
      this.y = numeric.rep(numeric.dim(this), 0);
      y = this.y;
    }
    b(x, from, to, Ax);
    b(y, from, to, Ay);
    return this;
  }
  b(x, from, to, Ax);
  if (y) b(y, from, to, numeric.rep(numeric.dim(Ax), 0));
};
numeric.T.rep = function rep(s, v) {
  var T = numeric.T;
  if (!(v instanceof T)) v = new T(v);
  var x = v.x,
    y = v.y,
    r = numeric.rep;
  if (y) return new T(r(s, x), r(s, y));
  return new T(r(s, x));
};
numeric.T.diag = function diag(d) {
  if (!(d instanceof numeric.T)) d = new numeric.T(d);
  var x = d.x,
    y = d.y,
    diag = numeric.diag;
  if (y) return new numeric.T(diag(x), diag(y));
  return new numeric.T(diag(x));
};
numeric.T.eig = function eig() {
  if (this.y) {
    throw new Error("eig: not implemented for complex matrices.");
  }
  return numeric.eig(this.x);
};
numeric.T.identity = function identity(n) {
  return new numeric.T(numeric.identity(n));
};
numeric.T.prototype.getDiag = function getDiag() {
  var n = numeric;
  var x = this.x,
    y = this.y;
  if (y) {
    return new n.T(n.getDiag(x), n.getDiag(y));
  }
  return new n.T(n.getDiag(x));
};

// 4. Eigenvalues of real matrices

numeric.house = function house(x) {
  var v = numeric.clone(x);
  var s = x[0] >= 0 ? 1 : -1;
  var alpha = s * numeric.norm2(x);
  v[0] += alpha;
  var foo = numeric.norm2(v);
  if (foo === 0) {
    /* this should not happen */ throw new Error("eig: internal error");
  }
  return numeric.div(v, foo);
};

numeric.toUpperHessenberg = function toUpperHessenberg(me) {
  var s = numeric.dim(me);
  if (s.length !== 2 || s[0] !== s[1]) {
    throw new Error(
      "numeric: toUpperHessenberg() only works on square matrices"
    );
  }
  var m = s[0],
    i,
    j,
    k,
    x,
    v,
    A = numeric.clone(me),
    B,
    C,
    Ai,
    Ci,
    Q = numeric.identity(m),
    Qi;
  for (j = 0; j < m - 2; j++) {
    x = Array(m - j - 1);
    for (i = j + 1; i < m; i++) {
      x[i - j - 1] = A[i][j];
    }
    if (numeric.norm2(x) > 0) {
      v = numeric.house(x);
      B = numeric.getBlock(A, [j + 1, j], [m - 1, m - 1]);
      C = numeric.tensor(v, numeric.dot(v, B));
      for (i = j + 1; i < m; i++) {
        Ai = A[i];
        Ci = C[i - j - 1];
        for (k = j; k < m; k++) Ai[k] -= 2 * Ci[k - j];
      }
      B = numeric.getBlock(A, [0, j + 1], [m - 1, m - 1]);
      C = numeric.tensor(numeric.dot(B, v), v);
      for (i = 0; i < m; i++) {
        Ai = A[i];
        Ci = C[i];
        for (k = j + 1; k < m; k++) Ai[k] -= 2 * Ci[k - j - 1];
      }
      B = Array(m - j - 1);
      for (i = j + 1; i < m; i++) B[i - j - 1] = Q[i];
      C = numeric.tensor(v, numeric.dot(v, B));
      for (i = j + 1; i < m; i++) {
        Qi = Q[i];
        Ci = C[i - j - 1];
        for (k = 0; k < m; k++) Qi[k] -= 2 * Ci[k];
      }
    }
  }
  return { H: A, Q: Q };
};

numeric.epsilon = 2.220446049250313e-16;

numeric.QRFrancis = function(H, maxiter) {
  if (typeof maxiter === "undefined") {
    maxiter = 10000;
  }
  H = numeric.clone(H);
  var H0 = numeric.clone(H);
  var s = numeric.dim(H),
    m = s[0],
    x,
    v,
    a,
    b,
    c,
    d,
    det,
    tr,
    Hloc,
    Q = numeric.identity(m),
    Qi,
    Hi,
    B,
    C,
    Ci,
    i,
    j,
    k,
    iter;
  if (m < 3) {
    return { Q: Q, B: [[0, m - 1]] };
  }
  var epsilon = numeric.epsilon;
  for (iter = 0; iter < maxiter; iter++) {
    for (j = 0; j < m - 1; j++) {
      if (
        Math.abs(H[j + 1][j]) <
        epsilon * (Math.abs(H[j][j]) + Math.abs(H[j + 1][j + 1]))
      ) {
        var QH1 = numeric.QRFrancis(
          numeric.getBlock(H, [0, 0], [j, j]),
          maxiter
        );
        var QH2 = numeric.QRFrancis(
          numeric.getBlock(H, [j + 1, j + 1], [m - 1, m - 1]),
          maxiter
        );
        B = Array(j + 1);
        for (i = 0; i <= j; i++) {
          B[i] = Q[i];
        }
        C = numeric.dot(QH1.Q, B);
        for (i = 0; i <= j; i++) {
          Q[i] = C[i];
        }
        B = Array(m - j - 1);
        for (i = j + 1; i < m; i++) {
          B[i - j - 1] = Q[i];
        }
        C = numeric.dot(QH2.Q, B);
        for (i = j + 1; i < m; i++) {
          Q[i] = C[i - j - 1];
        }
        return { Q: Q, B: QH1.B.concat(numeric.add(QH2.B, j + 1)) };
      }
    }
    a = H[m - 2][m - 2];
    b = H[m - 2][m - 1];
    c = H[m - 1][m - 2];
    d = H[m - 1][m - 1];
    tr = a + d;
    det = a * d - b * c;
    Hloc = numeric.getBlock(H, [0, 0], [2, 2]);
    if (tr * tr >= 4 * det) {
      var s1, s2;
      s1 = 0.5 * (tr + Math.sqrt(tr * tr - 4 * det));
      s2 = 0.5 * (tr - Math.sqrt(tr * tr - 4 * det));
      Hloc = numeric.add(
        numeric.sub(numeric.dot(Hloc, Hloc), numeric.mul(Hloc, s1 + s2)),
        numeric.diag(numeric.rep([3], s1 * s2))
      );
    } else {
      Hloc = numeric.add(
        numeric.sub(numeric.dot(Hloc, Hloc), numeric.mul(Hloc, tr)),
        numeric.diag(numeric.rep([3], det))
      );
    }
    x = [Hloc[0][0], Hloc[1][0], Hloc[2][0]];
    v = numeric.house(x);
    B = [H[0], H[1], H[2]];
    C = numeric.tensor(v, numeric.dot(v, B));
    for (i = 0; i < 3; i++) {
      Hi = H[i];
      Ci = C[i];
      for (k = 0; k < m; k++) Hi[k] -= 2 * Ci[k];
    }
    B = numeric.getBlock(H, [0, 0], [m - 1, 2]);
    C = numeric.tensor(numeric.dot(B, v), v);
    for (i = 0; i < m; i++) {
      Hi = H[i];
      Ci = C[i];
      for (k = 0; k < 3; k++) Hi[k] -= 2 * Ci[k];
    }
    B = [Q[0], Q[1], Q[2]];
    C = numeric.tensor(v, numeric.dot(v, B));
    for (i = 0; i < 3; i++) {
      Qi = Q[i];
      Ci = C[i];
      for (k = 0; k < m; k++) Qi[k] -= 2 * Ci[k];
    }
    var J;
    for (j = 0; j < m - 2; j++) {
      for (k = j; k <= j + 1; k++) {
        if (
          Math.abs(H[k + 1][k]) <
          epsilon * (Math.abs(H[k][k]) + Math.abs(H[k + 1][k + 1]))
        ) {
          var QH1 = numeric.QRFrancis(
            numeric.getBlock(H, [0, 0], [k, k]),
            maxiter
          );
          var QH2 = numeric.QRFrancis(
            numeric.getBlock(H, [k + 1, k + 1], [m - 1, m - 1]),
            maxiter
          );
          B = Array(k + 1);
          for (i = 0; i <= k; i++) {
            B[i] = Q[i];
          }
          C = numeric.dot(QH1.Q, B);
          for (i = 0; i <= k; i++) {
            Q[i] = C[i];
          }
          B = Array(m - k - 1);
          for (i = k + 1; i < m; i++) {
            B[i - k - 1] = Q[i];
          }
          C = numeric.dot(QH2.Q, B);
          for (i = k + 1; i < m; i++) {
            Q[i] = C[i - k - 1];
          }
          return { Q: Q, B: QH1.B.concat(numeric.add(QH2.B, k + 1)) };
        }
      }
      J = Math.min(m - 1, j + 3);
      x = Array(J - j);
      for (i = j + 1; i <= J; i++) {
        x[i - j - 1] = H[i][j];
      }
      v = numeric.house(x);
      B = numeric.getBlock(H, [j + 1, j], [J, m - 1]);
      C = numeric.tensor(v, numeric.dot(v, B));
      for (i = j + 1; i <= J; i++) {
        Hi = H[i];
        Ci = C[i - j - 1];
        for (k = j; k < m; k++) Hi[k] -= 2 * Ci[k - j];
      }
      B = numeric.getBlock(H, [0, j + 1], [m - 1, J]);
      C = numeric.tensor(numeric.dot(B, v), v);
      for (i = 0; i < m; i++) {
        Hi = H[i];
        Ci = C[i];
        for (k = j + 1; k <= J; k++) Hi[k] -= 2 * Ci[k - j - 1];
      }
      B = Array(J - j);
      for (i = j + 1; i <= J; i++) B[i - j - 1] = Q[i];
      C = numeric.tensor(v, numeric.dot(v, B));
      for (i = j + 1; i <= J; i++) {
        Qi = Q[i];
        Ci = C[i - j - 1];
        for (k = 0; k < m; k++) Qi[k] -= 2 * Ci[k];
      }
    }
  }
  throw new Error(
    "numeric: eigenvalue iteration does not converge -- increase maxiter?"
  );
};

numeric.eig = function eig(A, maxiter) {
  var QH = numeric.toUpperHessenberg(A);
  var QB = numeric.QRFrancis(QH.H, maxiter);
  var T = numeric.T;
  var n = A.length,
    i,
    k,
    flag = false,
    B = QB.B,
    H = numeric.dot(QB.Q, numeric.dot(QH.H, numeric.transpose(QB.Q)));
  var Q = new T(numeric.dot(QB.Q, QH.Q)),
    Q0;
  var m = B.length,
    j;
  var a, b, c, d, p1, p2, disc, x, y, p, q, n1, n2;
  var sqrt = Math.sqrt;
  for (k = 0; k < m; k++) {
    i = B[k][0];
    if (i === B[k][1]) {
      // nothing
    } else {
      j = i + 1;
      a = H[i][i];
      b = H[i][j];
      c = H[j][i];
      d = H[j][j];
      if (b === 0 && c === 0) continue;
      p1 = -a - d;
      p2 = a * d - b * c;
      disc = p1 * p1 - 4 * p2;
      if (disc >= 0) {
        if (p1 < 0) x = -0.5 * (p1 - sqrt(disc));
        else x = -0.5 * (p1 + sqrt(disc));
        n1 = (a - x) * (a - x) + b * b;
        n2 = c * c + (d - x) * (d - x);
        if (n1 > n2) {
          n1 = sqrt(n1);
          p = (a - x) / n1;
          q = b / n1;
        } else {
          n2 = sqrt(n2);
          p = c / n2;
          q = (d - x) / n2;
        }
        Q0 = new T([[q, -p], [p, q]]);
        Q.setRows(i, j, Q0.dot(Q.getRows(i, j)));
      } else {
        x = -0.5 * p1;
        y = 0.5 * sqrt(-disc);
        n1 = (a - x) * (a - x) + b * b;
        n2 = c * c + (d - x) * (d - x);
        if (n1 > n2) {
          n1 = sqrt(n1 + y * y);
          p = (a - x) / n1;
          q = b / n1;
          x = 0;
          y /= n1;
        } else {
          n2 = sqrt(n2 + y * y);
          p = c / n2;
          q = (d - x) / n2;
          x = y / n2;
          y = 0;
        }
        Q0 = new T([[q, -p], [p, q]], [[x, y], [y, -x]]);
        Q.setRows(i, j, Q0.dot(Q.getRows(i, j)));
      }
    }
  }
  var R = Q.dot(A).dot(Q.transjugate()),
    n = A.length,
    E = numeric.T.identity(n);
  for (j = 0; j < n; j++) {
    if (j > 0) {
      for (k = j - 1; k >= 0; k--) {
        var Rk = R.get([k, k]),
          Rj = R.get([j, j]);
        if (numeric.neq(Rk.x, Rj.x) || numeric.neq(Rk.y, Rj.y)) {
          x = R.getRow(k).getBlock([k], [j - 1]);
          y = E.getRow(j).getBlock([k], [j - 1]);
          E.set(
            [j, k],
            R.get([k, j])
              .neg()
              .sub(x.dot(y))
              .div(Rk.sub(Rj))
          );
        } else {
          E.setRow(j, E.getRow(k));
          continue;
        }
      }
    }
  }
  for (j = 0; j < n; j++) {
    x = E.getRow(j);
    E.setRow(j, x.div(x.norm2()));
  }
  E = E.transpose();
  E = Q.transjugate().dot(E);
  return { lambda: R.getDiag(), E: E };
};

// 5. Compressed Column Storage matrices
numeric.ccsSparse = function ccsSparse(A) {
  var m = A.length,
    n,
    foo,
    i,
    j,
    counts = [];
  for (i = m - 1; i !== -1; --i) {
    foo = A[i];
    for (j in foo) {
      j = parseInt(j);
      while (j >= counts.length) counts[counts.length] = 0;
      if (foo[j] !== 0) counts[j]++;
    }
  }
  var n = counts.length;
  var Ai = Array(n + 1);
  Ai[0] = 0;
  for (i = 0; i < n; ++i) Ai[i + 1] = Ai[i] + counts[i];
  var Aj = Array(Ai[n]),
    Av = Array(Ai[n]);
  for (i = m - 1; i !== -1; --i) {
    foo = A[i];
    for (j in foo) {
      if (foo[j] !== 0) {
        counts[j]--;
        Aj[Ai[j] + counts[j]] = i;
        Av[Ai[j] + counts[j]] = foo[j];
      }
    }
  }
  return [Ai, Aj, Av];
};
numeric.ccsFull = function ccsFull(A) {
  var Ai = A[0],
    Aj = A[1],
    Av = A[2],
    s = numeric.ccsDim(A),
    m = s[0],
    n = s[1],
    i,
    j,
    j0,
    j1,
    k;
  var B = numeric.rep([m, n], 0);
  for (i = 0; i < n; i++) {
    j0 = Ai[i];
    j1 = Ai[i + 1];
    for (j = j0; j < j1; ++j) {
      B[Aj[j]][i] = Av[j];
    }
  }
  return B;
};
numeric.ccsTSolve = function ccsTSolve(A, b, x, bj, xj) {
  var Ai = A[0],
    Aj = A[1],
    Av = A[2],
    m = Ai.length - 1,
    max = Math.max,
    n = 0;
  if (typeof bj === "undefined") x = numeric.rep([m], 0);
  if (typeof bj === "undefined") bj = numeric.linspace(0, x.length - 1);
  if (typeof xj === "undefined") xj = [];
  function dfs(j) {
    var k;
    if (x[j] !== 0) return;
    x[j] = 1;
    for (k = Ai[j]; k < Ai[j + 1]; ++k) dfs(Aj[k]);
    xj[n] = j;
    ++n;
  }
  var i, j, j0, j1, k, l, l0, l1, a;
  for (i = bj.length - 1; i !== -1; --i) {
    dfs(bj[i]);
  }
  xj.length = n;
  for (i = xj.length - 1; i !== -1; --i) {
    x[xj[i]] = 0;
  }
  for (i = bj.length - 1; i !== -1; --i) {
    j = bj[i];
    x[j] = b[j];
  }
  for (i = xj.length - 1; i !== -1; --i) {
    j = xj[i];
    j0 = Ai[j];
    j1 = max(Ai[j + 1], j0);
    for (k = j0; k !== j1; ++k) {
      if (Aj[k] === j) {
        x[j] /= Av[k];
        break;
      }
    }
    a = x[j];
    for (k = j0; k !== j1; ++k) {
      l = Aj[k];
      if (l !== j) x[l] -= a * Av[k];
    }
  }
  return x;
};
numeric.ccsDFS = function ccsDFS(n) {
  this.k = Array(n);
  this.k1 = Array(n);
  this.j = Array(n);
};
numeric.ccsDFS.prototype.dfs = function dfs(J, Ai, Aj, x, xj, Pinv) {
  var m = 0,
    foo,
    n = xj.length;
  var k = this.k,
    k1 = this.k1,
    j = this.j,
    km,
    k11;
  if (x[J] !== 0) return;
  x[J] = 1;
  j[0] = J;
  k[0] = km = Ai[J];
  k1[0] = k11 = Ai[J + 1];
  while (1) {
    if (km >= k11) {
      xj[n] = j[m];
      if (m === 0) return;
      ++n;
      --m;
      km = k[m];
      k11 = k1[m];
    } else {
      foo = Pinv[Aj[km]];
      if (x[foo] === 0) {
        x[foo] = 1;
        k[m] = km;
        ++m;
        j[m] = foo;
        km = Ai[foo];
        k1[m] = k11 = Ai[foo + 1];
      } else ++km;
    }
  }
};
numeric.ccsLPSolve = function ccsLPSolve(A, B, x, xj, I, Pinv, dfs) {
  var Ai = A[0],
    Aj = A[1],
    Av = A[2],
    m = Ai.length - 1,
    n = 0;
  var Bi = B[0],
    Bj = B[1],
    Bv = B[2];

  var i, i0, i1, j, J, j0, j1, k, l, l0, l1, a;
  i0 = Bi[I];
  i1 = Bi[I + 1];
  xj.length = 0;
  for (i = i0; i < i1; ++i) {
    dfs.dfs(Pinv[Bj[i]], Ai, Aj, x, xj, Pinv);
  }
  for (i = xj.length - 1; i !== -1; --i) {
    x[xj[i]] = 0;
  }
  for (i = i0; i !== i1; ++i) {
    j = Pinv[Bj[i]];
    x[j] = Bv[i];
  }
  for (i = xj.length - 1; i !== -1; --i) {
    j = xj[i];
    j0 = Ai[j];
    j1 = Ai[j + 1];
    for (k = j0; k < j1; ++k) {
      if (Pinv[Aj[k]] === j) {
        x[j] /= Av[k];
        break;
      }
    }
    a = x[j];
    for (k = j0; k < j1; ++k) {
      l = Pinv[Aj[k]];
      if (l !== j) x[l] -= a * Av[k];
    }
  }
  return x;
};
numeric.ccsLUP1 = function ccsLUP1(A, threshold) {
  var m = A[0].length - 1;
  var L = [numeric.rep([m + 1], 0), [], []],
    U = [numeric.rep([m + 1], 0), [], []];
  var Li = L[0],
    Lj = L[1],
    Lv = L[2],
    Ui = U[0],
    Uj = U[1],
    Uv = U[2];
  var x = numeric.rep([m], 0),
    xj = numeric.rep([m], 0);
  var i, j, k, j0, j1, a, e, c, d, K;
  var sol = numeric.ccsLPSolve,
    max = Math.max,
    abs = Math.abs;
  var P = numeric.linspace(0, m - 1),
    Pinv = numeric.linspace(0, m - 1);
  var dfs = new numeric.ccsDFS(m);
  if (typeof threshold === "undefined") {
    threshold = 1;
  }
  for (i = 0; i < m; ++i) {
    sol(L, A, x, xj, i, Pinv, dfs);
    a = -1;
    e = -1;
    for (j = xj.length - 1; j !== -1; --j) {
      k = xj[j];
      if (k <= i) continue;
      c = abs(x[k]);
      if (c > a) {
        e = k;
        a = c;
      }
    }
    if (abs(x[i]) < threshold * a) {
      j = P[i];
      a = P[e];
      P[i] = a;
      Pinv[a] = i;
      P[e] = j;
      Pinv[j] = e;
      a = x[i];
      x[i] = x[e];
      x[e] = a;
    }
    a = Li[i];
    e = Ui[i];
    d = x[i];
    Lj[a] = P[i];
    Lv[a] = 1;
    ++a;
    for (j = xj.length - 1; j !== -1; --j) {
      k = xj[j];
      c = x[k];
      xj[j] = 0;
      x[k] = 0;
      if (k <= i) {
        Uj[e] = k;
        Uv[e] = c;
        ++e;
      } else {
        Lj[a] = P[k];
        Lv[a] = c / d;
        ++a;
      }
    }
    Li[i + 1] = a;
    Ui[i + 1] = e;
  }
  for (j = Lj.length - 1; j !== -1; --j) {
    Lj[j] = Pinv[Lj[j]];
  }
  return { L: L, U: U, P: P, Pinv: Pinv };
};
numeric.ccsDFS0 = function ccsDFS0(n) {
  this.k = Array(n);
  this.k1 = Array(n);
  this.j = Array(n);
};
numeric.ccsDFS0.prototype.dfs = function dfs(J, Ai, Aj, x, xj, Pinv, P) {
  var m = 0,
    foo,
    n = xj.length;
  var k = this.k,
    k1 = this.k1,
    j = this.j,
    km,
    k11;
  if (x[J] !== 0) return;
  x[J] = 1;
  j[0] = J;
  k[0] = km = Ai[Pinv[J]];
  k1[0] = k11 = Ai[Pinv[J] + 1];
  while (1) {
    if (isNaN(km)) throw new Error("Ow!");
    if (km >= k11) {
      xj[n] = Pinv[j[m]];
      if (m === 0) return;
      ++n;
      --m;
      km = k[m];
      k11 = k1[m];
    } else {
      foo = Aj[km];
      if (x[foo] === 0) {
        x[foo] = 1;
        k[m] = km;
        ++m;
        j[m] = foo;
        foo = Pinv[foo];
        km = Ai[foo];
        k1[m] = k11 = Ai[foo + 1];
      } else ++km;
    }
  }
};
numeric.ccsLPSolve0 = function ccsLPSolve0(A, B, y, xj, I, Pinv, P, dfs) {
  var Ai = A[0],
    Aj = A[1],
    Av = A[2],
    m = Ai.length - 1,
    n = 0;
  var Bi = B[0],
    Bj = B[1],
    Bv = B[2];

  var i, i0, i1, j, J, j0, j1, k, l, l0, l1, a;
  i0 = Bi[I];
  i1 = Bi[I + 1];
  xj.length = 0;
  for (i = i0; i < i1; ++i) {
    dfs.dfs(Bj[i], Ai, Aj, y, xj, Pinv, P);
  }
  for (i = xj.length - 1; i !== -1; --i) {
    j = xj[i];
    y[P[j]] = 0;
  }
  for (i = i0; i !== i1; ++i) {
    j = Bj[i];
    y[j] = Bv[i];
  }
  for (i = xj.length - 1; i !== -1; --i) {
    j = xj[i];
    l = P[j];
    j0 = Ai[j];
    j1 = Ai[j + 1];
    for (k = j0; k < j1; ++k) {
      if (Aj[k] === l) {
        y[l] /= Av[k];
        break;
      }
    }
    a = y[l];
    for (k = j0; k < j1; ++k) y[Aj[k]] -= a * Av[k];
    y[l] = a;
  }
};
numeric.ccsLUP0 = function ccsLUP0(A, threshold) {
  var m = A[0].length - 1;
  var L = [numeric.rep([m + 1], 0), [], []],
    U = [numeric.rep([m + 1], 0), [], []];
  var Li = L[0],
    Lj = L[1],
    Lv = L[2],
    Ui = U[0],
    Uj = U[1],
    Uv = U[2];
  var y = numeric.rep([m], 0),
    xj = numeric.rep([m], 0);
  var i, j, k, j0, j1, a, e, c, d, K;
  var sol = numeric.ccsLPSolve0,
    max = Math.max,
    abs = Math.abs;
  var P = numeric.linspace(0, m - 1),
    Pinv = numeric.linspace(0, m - 1);
  var dfs = new numeric.ccsDFS0(m);
  if (typeof threshold === "undefined") {
    threshold = 1;
  }
  for (i = 0; i < m; ++i) {
    sol(L, A, y, xj, i, Pinv, P, dfs);
    a = -1;
    e = -1;
    for (j = xj.length - 1; j !== -1; --j) {
      k = xj[j];
      if (k <= i) continue;
      c = abs(y[P[k]]);
      if (c > a) {
        e = k;
        a = c;
      }
    }
    if (abs(y[P[i]]) < threshold * a) {
      j = P[i];
      a = P[e];
      P[i] = a;
      Pinv[a] = i;
      P[e] = j;
      Pinv[j] = e;
    }
    a = Li[i];
    e = Ui[i];
    d = y[P[i]];
    Lj[a] = P[i];
    Lv[a] = 1;
    ++a;
    for (j = xj.length - 1; j !== -1; --j) {
      k = xj[j];
      c = y[P[k]];
      xj[j] = 0;
      y[P[k]] = 0;
      if (k <= i) {
        Uj[e] = k;
        Uv[e] = c;
        ++e;
      } else {
        Lj[a] = P[k];
        Lv[a] = c / d;
        ++a;
      }
    }
    Li[i + 1] = a;
    Ui[i + 1] = e;
  }
  for (j = Lj.length - 1; j !== -1; --j) {
    Lj[j] = Pinv[Lj[j]];
  }
  return { L: L, U: U, P: P, Pinv: Pinv };
};
numeric.ccsLUP = numeric.ccsLUP0;

numeric.ccsDim = function ccsDim(A) {
  return [numeric.sup(A[1]) + 1, A[0].length - 1];
};
numeric.ccsGetBlock = function ccsGetBlock(A, i, j) {
  var s = numeric.ccsDim(A),
    m = s[0],
    n = s[1];
  if (typeof i === "undefined") {
    i = numeric.linspace(0, m - 1);
  } else if (typeof i === "number") {
    i = [i];
  }
  if (typeof j === "undefined") {
    j = numeric.linspace(0, n - 1);
  } else if (typeof j === "number") {
    j = [j];
  }
  var p,
    p0,
    p1,
    P = i.length,
    q,
    Q = j.length,
    r,
    jq,
    ip;
  var Bi = numeric.rep([n], 0),
    Bj = [],
    Bv = [],
    B = [Bi, Bj, Bv];
  var Ai = A[0],
    Aj = A[1],
    Av = A[2];
  var x = numeric.rep([m], 0),
    count = 0,
    flags = numeric.rep([m], 0);
  for (q = 0; q < Q; ++q) {
    jq = j[q];
    var q0 = Ai[jq];
    var q1 = Ai[jq + 1];
    for (p = q0; p < q1; ++p) {
      r = Aj[p];
      flags[r] = 1;
      x[r] = Av[p];
    }
    for (p = 0; p < P; ++p) {
      ip = i[p];
      if (flags[ip]) {
        Bj[count] = p;
        Bv[count] = x[i[p]];
        ++count;
      }
    }
    for (p = q0; p < q1; ++p) {
      r = Aj[p];
      flags[r] = 0;
    }
    Bi[q + 1] = count;
  }
  return B;
};

numeric.ccsDot = function ccsDot(A, B) {
  var Ai = A[0],
    Aj = A[1],
    Av = A[2];
  var Bi = B[0],
    Bj = B[1],
    Bv = B[2];
  var sA = numeric.ccsDim(A),
    sB = numeric.ccsDim(B);
  var m = sA[0],
    n = sA[1],
    o = sB[1];
  var x = numeric.rep([m], 0),
    flags = numeric.rep([m], 0),
    xj = Array(m);
  var Ci = numeric.rep([o], 0),
    Cj = [],
    Cv = [],
    C = [Ci, Cj, Cv];
  var i, j, k, j0, j1, i0, i1, l, p, a, b;
  for (k = 0; k !== o; ++k) {
    j0 = Bi[k];
    j1 = Bi[k + 1];
    p = 0;
    for (j = j0; j < j1; ++j) {
      a = Bj[j];
      b = Bv[j];
      i0 = Ai[a];
      i1 = Ai[a + 1];
      for (i = i0; i < i1; ++i) {
        l = Aj[i];
        if (flags[l] === 0) {
          xj[p] = l;
          flags[l] = 1;
          p = p + 1;
        }
        x[l] = x[l] + Av[i] * b;
      }
    }
    j0 = Ci[k];
    j1 = j0 + p;
    Ci[k + 1] = j1;
    for (j = p - 1; j !== -1; --j) {
      b = j0 + j;
      i = xj[j];
      Cj[b] = i;
      Cv[b] = x[i];
      flags[i] = 0;
      x[i] = 0;
    }
    Ci[k + 1] = Ci[k] + p;
  }
  return C;
};

numeric.ccsLUPSolve = function ccsLUPSolve(LUP, B) {
  var L = LUP.L,
    U = LUP.U,
    P = LUP.P;
  var Bi = B[0];
  var flag = false;
  if (typeof Bi !== "object") {
    B = [[0, B.length], numeric.linspace(0, B.length - 1), B];
    Bi = B[0];
    flag = true;
  }
  var Bj = B[1],
    Bv = B[2];
  var n = L[0].length - 1,
    m = Bi.length - 1;
  var x = numeric.rep([n], 0),
    xj = Array(n);
  var b = numeric.rep([n], 0),
    bj = Array(n);
  var Xi = numeric.rep([m + 1], 0),
    Xj = [],
    Xv = [];
  var sol = numeric.ccsTSolve;
  var i,
    j,
    j0,
    j1,
    k,
    J,
    N = 0;
  for (i = 0; i < m; ++i) {
    k = 0;
    j0 = Bi[i];
    j1 = Bi[i + 1];
    for (j = j0; j < j1; ++j) {
      J = LUP.Pinv[Bj[j]];
      bj[k] = J;
      b[J] = Bv[j];
      ++k;
    }
    bj.length = k;
    sol(L, b, x, bj, xj);
    for (j = bj.length - 1; j !== -1; --j) b[bj[j]] = 0;
    sol(U, x, b, xj, bj);
    if (flag) return b;
    for (j = xj.length - 1; j !== -1; --j) x[xj[j]] = 0;
    for (j = bj.length - 1; j !== -1; --j) {
      J = bj[j];
      Xj[N] = J;
      Xv[N] = b[J];
      b[J] = 0;
      ++N;
    }
    Xi[i + 1] = N;
  }
  return [Xi, Xj, Xv];
};

numeric.ccsbinop = function ccsbinop(body, setup) {
  if (typeof setup === "undefined") setup = "";
  return Function(
    "X",
    "Y",
    "var Xi = X[0], Xj = X[1], Xv = X[2];\n" +
      "var Yi = Y[0], Yj = Y[1], Yv = Y[2];\n" +
      "var n = Xi.length-1,m = Math.max(numeric.sup(Xj),numeric.sup(Yj))+1;\n" +
      "var Zi = numeric.rep([n+1],0), Zj = [], Zv = [];\n" +
      "var x = numeric.rep([m],0),y = numeric.rep([m],0);\n" +
      "var xk,yk,zk;\n" +
      "var i,j,j0,j1,k,p=0;\n" +
      setup +
      "for(i=0;i<n;++i) {\n" +
      "  j0 = Xi[i]; j1 = Xi[i+1];\n" +
      "  for(j=j0;j!==j1;++j) {\n" +
      "    k = Xj[j];\n" +
      "    x[k] = 1;\n" +
      "    Zj[p] = k;\n" +
      "    ++p;\n" +
      "  }\n" +
      "  j0 = Yi[i]; j1 = Yi[i+1];\n" +
      "  for(j=j0;j!==j1;++j) {\n" +
      "    k = Yj[j];\n" +
      "    y[k] = Yv[j];\n" +
      "    if(x[k] === 0) {\n" +
      "      Zj[p] = k;\n" +
      "      ++p;\n" +
      "    }\n" +
      "  }\n" +
      "  Zi[i+1] = p;\n" +
      "  j0 = Xi[i]; j1 = Xi[i+1];\n" +
      "  for(j=j0;j!==j1;++j) x[Xj[j]] = Xv[j];\n" +
      "  j0 = Zi[i]; j1 = Zi[i+1];\n" +
      "  for(j=j0;j!==j1;++j) {\n" +
      "    k = Zj[j];\n" +
      "    xk = x[k];\n" +
      "    yk = y[k];\n" +
      body +
      "\n" +
      "    Zv[j] = zk;\n" +
      "  }\n" +
      "  j0 = Xi[i]; j1 = Xi[i+1];\n" +
      "  for(j=j0;j!==j1;++j) x[Xj[j]] = 0;\n" +
      "  j0 = Yi[i]; j1 = Yi[i+1];\n" +
      "  for(j=j0;j!==j1;++j) y[Yj[j]] = 0;\n" +
      "}\n" +
      "return [Zi,Zj,Zv];"
  );
};

(function() {
  var k, A, B, C;
  for (k in numeric.ops2) {
    if (isFinite(eval("1" + numeric.ops2[k] + "0")))
      A = "[Y[0],Y[1],numeric." + k + "(X,Y[2])]";
    else A = "NaN";
    if (isFinite(eval("0" + numeric.ops2[k] + "1")))
      B = "[X[0],X[1],numeric." + k + "(X[2],Y)]";
    else B = "NaN";
    if (
      isFinite(eval("1" + numeric.ops2[k] + "0")) &&
      isFinite(eval("0" + numeric.ops2[k] + "1"))
    )
      C = "numeric.ccs" + k + "MM(X,Y)";
    else C = "NaN";
    numeric["ccs" + k + "MM"] = numeric.ccsbinop(
      "zk = xk " + numeric.ops2[k] + "yk;"
    );
    numeric["ccs" + k] = Function(
      "X",
      "Y",
      'if(typeof X === "number") return ' +
        A +
        ";\n" +
        'if(typeof Y === "number") return ' +
        B +
        ";\n" +
        "return " +
        C +
        ";\n"
    );
  }
})();

numeric.ccsScatter = function ccsScatter(A) {
  var Ai = A[0],
    Aj = A[1],
    Av = A[2];
  var n = numeric.sup(Aj) + 1,
    m = Ai.length;
  var Ri = numeric.rep([n], 0),
    Rj = Array(m),
    Rv = Array(m);
  var counts = numeric.rep([n], 0),
    i;
  for (i = 0; i < m; ++i) counts[Aj[i]]++;
  for (i = 0; i < n; ++i) Ri[i + 1] = Ri[i] + counts[i];
  var ptr = Ri.slice(0),
    k,
    Aii;
  for (i = 0; i < m; ++i) {
    Aii = Aj[i];
    k = ptr[Aii];
    Rj[k] = Ai[i];
    Rv[k] = Av[i];
    ptr[Aii] = ptr[Aii] + 1;
  }
  return [Ri, Rj, Rv];
};

numeric.ccsGather = function ccsGather(A) {
  var Ai = A[0],
    Aj = A[1],
    Av = A[2];
  var n = Ai.length - 1,
    m = Aj.length;
  var Ri = Array(m),
    Rj = Array(m),
    Rv = Array(m);
  var i, j, j0, j1, p;
  p = 0;
  for (i = 0; i < n; ++i) {
    j0 = Ai[i];
    j1 = Ai[i + 1];
    for (j = j0; j !== j1; ++j) {
      Rj[p] = i;
      Ri[p] = Aj[j];
      Rv[p] = Av[j];
      ++p;
    }
  }
  return [Ri, Rj, Rv];
};

// The following sparse linear algebra routines are deprecated.

numeric.sdim = function dim(A, ret, k) {
  if (typeof ret === "undefined") {
    ret = [];
  }
  if (typeof A !== "object") return ret;
  if (typeof k === "undefined") {
    k = 0;
  }
  if (!(k in ret)) {
    ret[k] = 0;
  }
  if (A.length > ret[k]) ret[k] = A.length;
  var i;
  for (i in A) {
    if (A.hasOwnProperty(i)) dim(A[i], ret, k + 1);
  }
  return ret;
};

numeric.sclone = function clone(A, k, n) {
  if (typeof k === "undefined") {
    k = 0;
  }
  if (typeof n === "undefined") {
    n = numeric.sdim(A).length;
  }
  var i,
    ret = Array(A.length);
  if (k === n - 1) {
    for (i in A) {
      if (A.hasOwnProperty(i)) ret[i] = A[i];
    }
    return ret;
  }
  for (i in A) {
    if (A.hasOwnProperty(i)) ret[i] = clone(A[i], k + 1, n);
  }
  return ret;
};

numeric.sdiag = function diag(d) {
  var n = d.length,
    i,
    ret = Array(n),
    i1,
    i2,
    i3;
  for (i = n - 1; i >= 1; i -= 2) {
    i1 = i - 1;
    ret[i] = [];
    ret[i][i] = d[i];
    ret[i1] = [];
    ret[i1][i1] = d[i1];
  }
  if (i === 0) {
    ret[0] = [];
    ret[0][0] = d[i];
  }
  return ret;
};

numeric.sidentity = function identity(n) {
  return numeric.sdiag(numeric.rep([n], 1));
};

numeric.stranspose = function transpose(A) {
  var ret = [],
    n = A.length,
    i,
    j,
    Ai;
  for (i in A) {
    if (!A.hasOwnProperty(i)) continue;
    Ai = A[i];
    for (j in Ai) {
      if (!Ai.hasOwnProperty(j)) continue;
      if (typeof ret[j] !== "object") {
        ret[j] = [];
      }
      ret[j][i] = Ai[j];
    }
  }
  return ret;
};

numeric.sLUP = function LUP(A, tol) {
  throw new Error(
    "The function numeric.sLUP had a bug in it and has been removed. Please use the new numeric.ccsLUP function instead."
  );
};

numeric.sdotMM = function dotMM(A, B) {
  var p = A.length,
    q = B.length,
    BT = numeric.stranspose(B),
    r = BT.length,
    Ai,
    BTk;
  var i, j, k, accum;
  var ret = Array(p),
    reti;
  for (i = p - 1; i >= 0; i--) {
    reti = [];
    Ai = A[i];
    for (k = r - 1; k >= 0; k--) {
      accum = 0;
      BTk = BT[k];
      for (j in Ai) {
        if (!Ai.hasOwnProperty(j)) continue;
        if (j in BTk) {
          accum += Ai[j] * BTk[j];
        }
      }
      if (accum) reti[k] = accum;
    }
    ret[i] = reti;
  }
  return ret;
};

numeric.sdotMV = function dotMV(A, x) {
  var p = A.length,
    Ai,
    i,
    j;
  var ret = Array(p),
    accum;
  for (i = p - 1; i >= 0; i--) {
    Ai = A[i];
    accum = 0;
    for (j in Ai) {
      if (!Ai.hasOwnProperty(j)) continue;
      if (x[j]) accum += Ai[j] * x[j];
    }
    if (accum) ret[i] = accum;
  }
  return ret;
};

numeric.sdotVM = function dotMV(x, A) {
  var i, j, Ai, alpha;
  var ret = [],
    accum;
  for (i in x) {
    if (!x.hasOwnProperty(i)) continue;
    Ai = A[i];
    alpha = x[i];
    for (j in Ai) {
      if (!Ai.hasOwnProperty(j)) continue;
      if (!ret[j]) {
        ret[j] = 0;
      }
      ret[j] += alpha * Ai[j];
    }
  }
  return ret;
};

numeric.sdotVV = function dotVV(x, y) {
  var i,
    ret = 0;
  for (i in x) {
    if (x[i] && y[i]) ret += x[i] * y[i];
  }
  return ret;
};

numeric.sdot = function dot(A, B) {
  var m = numeric.sdim(A).length,
    n = numeric.sdim(B).length;
  var k = m * 1000 + n;
  switch (k) {
    case 0:
      return A * B;
    case 1001:
      return numeric.sdotVV(A, B);
    case 2001:
      return numeric.sdotMV(A, B);
    case 1002:
      return numeric.sdotVM(A, B);
    case 2002:
      return numeric.sdotMM(A, B);
    default:
      throw new Error(
        "numeric.sdot not implemented for tensors of order " + m + " and " + n
      );
  }
};

numeric.sscatter = function scatter(V) {
  var n = V[0].length,
    Vij,
    i,
    j,
    m = V.length,
    A = [],
    Aj;
  for (i = n - 1; i >= 0; --i) {
    if (!V[m - 1][i]) continue;
    Aj = A;
    for (j = 0; j < m - 2; j++) {
      Vij = V[j][i];
      if (!Aj[Vij]) Aj[Vij] = [];
      Aj = Aj[Vij];
    }
    Aj[V[j][i]] = V[j + 1][i];
  }
  return A;
};

numeric.sgather = function gather(A, ret, k) {
  if (typeof ret === "undefined") ret = [];
  if (typeof k === "undefined") k = [];
  var n, i, Ai;
  n = k.length;
  for (i in A) {
    if (A.hasOwnProperty(i)) {
      k[n] = parseInt(i);
      Ai = A[i];
      if (typeof Ai === "number") {
        if (Ai) {
          if (ret.length === 0) {
            for (i = n + 1; i >= 0; --i) ret[i] = [];
          }
          for (i = n; i >= 0; --i) ret[i].push(k[i]);
          ret[n + 1].push(Ai);
        }
      } else gather(Ai, ret, k);
    }
  }
  if (k.length > n) k.pop();
  return ret;
};

// 6. Coordinate matrices
numeric.cLU = function LU(A) {
  var I = A[0],
    J = A[1],
    V = A[2];
  var p = I.length,
    m = 0,
    i,
    j,
    k,
    a,
    b,
    c;
  for (i = 0; i < p; i++) if (I[i] > m) m = I[i];
  m++;
  var L = Array(m),
    U = Array(m),
    left = numeric.rep([m], Infinity),
    right = numeric.rep([m], -Infinity);
  var Ui, Uj, alpha;
  for (k = 0; k < p; k++) {
    i = I[k];
    j = J[k];
    if (j < left[i]) left[i] = j;
    if (j > right[i]) right[i] = j;
  }
  for (i = 0; i < m - 1; i++) {
    if (right[i] > right[i + 1]) right[i + 1] = right[i];
  }
  for (i = m - 1; i >= 1; i--) {
    if (left[i] < left[i - 1]) left[i - 1] = left[i];
  }
  var countL = 0,
    countU = 0;
  for (i = 0; i < m; i++) {
    U[i] = numeric.rep([right[i] - left[i] + 1], 0);
    L[i] = numeric.rep([i - left[i]], 0);
    countL += i - left[i] + 1;
    countU += right[i] - i + 1;
  }
  for (k = 0; k < p; k++) {
    i = I[k];
    U[i][J[k] - left[i]] = V[k];
  }
  for (i = 0; i < m - 1; i++) {
    a = i - left[i];
    Ui = U[i];
    for (j = i + 1; left[j] <= i && j < m; j++) {
      b = i - left[j];
      c = right[i] - i;
      Uj = U[j];
      alpha = Uj[b] / Ui[a];
      if (alpha) {
        for (k = 1; k <= c; k++) {
          Uj[k + b] -= alpha * Ui[k + a];
        }
        L[j][i - left[j]] = alpha;
      }
    }
  }
  var Ui = [],
    Uj = [],
    Uv = [],
    Li = [],
    Lj = [],
    Lv = [];
  var p, q, foo;
  p = 0;
  q = 0;
  for (i = 0; i < m; i++) {
    a = left[i];
    b = right[i];
    foo = U[i];
    for (j = i; j <= b; j++) {
      if (foo[j - a]) {
        Ui[p] = i;
        Uj[p] = j;
        Uv[p] = foo[j - a];
        p++;
      }
    }
    foo = L[i];
    for (j = a; j < i; j++) {
      if (foo[j - a]) {
        Li[q] = i;
        Lj[q] = j;
        Lv[q] = foo[j - a];
        q++;
      }
    }
    Li[q] = i;
    Lj[q] = i;
    Lv[q] = 1;
    q++;
  }
  return { U: [Ui, Uj, Uv], L: [Li, Lj, Lv] };
};

numeric.cLUsolve = function LUsolve(lu, b) {
  var L = lu.L,
    U = lu.U,
    ret = numeric.clone(b);
  var Li = L[0],
    Lj = L[1],
    Lv = L[2];
  var Ui = U[0],
    Uj = U[1],
    Uv = U[2];
  var p = Ui.length,
    q = Li.length;
  var m = ret.length,
    i,
    j,
    k;
  k = 0;
  for (i = 0; i < m; i++) {
    while (Lj[k] < i) {
      ret[i] -= Lv[k] * ret[Lj[k]];
      k++;
    }
    k++;
  }
  k = p - 1;
  for (i = m - 1; i >= 0; i--) {
    while (Uj[k] > i) {
      ret[i] -= Uv[k] * ret[Uj[k]];
      k--;
    }
    ret[i] /= Uv[k];
    k--;
  }
  return ret;
};

numeric.cgrid = function grid(n, shape) {
  if (typeof n === "number") n = [n, n];
  var ret = numeric.rep(n, -1);
  var i, j, count;
  if (typeof shape !== "function") {
    switch (shape) {
      case "L":
        shape = function(i, j) {
          return i >= n[0] / 2 || j < n[1] / 2;
        };
        break;
      default:
        shape = function(i, j) {
          return true;
        };
        break;
    }
  }
  count = 0;
  for (i = 1; i < n[0] - 1; i++)
    for (j = 1; j < n[1] - 1; j++)
      if (shape(i, j)) {
        ret[i][j] = count;
        count++;
      }
  return ret;
};

numeric.cdelsq = function delsq(g) {
  var dir = [[-1, 0], [0, -1], [0, 1], [1, 0]];
  var s = numeric.dim(g),
    m = s[0],
    n = s[1],
    i,
    j,
    k,
    p,
    q;
  var Li = [],
    Lj = [],
    Lv = [];
  for (i = 1; i < m - 1; i++)
    for (j = 1; j < n - 1; j++) {
      if (g[i][j] < 0) continue;
      for (k = 0; k < 4; k++) {
        p = i + dir[k][0];
        q = j + dir[k][1];
        if (g[p][q] < 0) continue;
        Li.push(g[i][j]);
        Lj.push(g[p][q]);
        Lv.push(-1);
      }
      Li.push(g[i][j]);
      Lj.push(g[i][j]);
      Lv.push(4);
    }
  return [Li, Lj, Lv];
};

numeric.cdotMV = function dotMV(A, x) {
  var ret,
    Ai = A[0],
    Aj = A[1],
    Av = A[2],
    k,
    p = Ai.length,
    N;
  N = 0;
  for (k = 0; k < p; k++) {
    if (Ai[k] > N) N = Ai[k];
  }
  N++;
  ret = numeric.rep([N], 0);
  for (k = 0; k < p; k++) {
    ret[Ai[k]] += Av[k] * x[Aj[k]];
  }
  return ret;
};

// 7. Splines

numeric.Spline = function Spline(x, yl, yr, kl, kr) {
  this.x = x;
  this.yl = yl;
  this.yr = yr;
  this.kl = kl;
  this.kr = kr;
};
numeric.Spline.prototype._at = function _at(x1, p) {
  var x = this.x;
  var yl = this.yl;
  var yr = this.yr;
  var kl = this.kl;
  var kr = this.kr;
  var x1, a, b, t;
  var add = numeric.add,
    sub = numeric.sub,
    mul = numeric.mul;
  a = sub(mul(kl[p], x[p + 1] - x[p]), sub(yr[p + 1], yl[p]));
  b = add(mul(kr[p + 1], x[p] - x[p + 1]), sub(yr[p + 1], yl[p]));
  t = (x1 - x[p]) / (x[p + 1] - x[p]);
  var s = t * (1 - t);
  return add(
    add(add(mul(1 - t, yl[p]), mul(t, yr[p + 1])), mul(a, s * (1 - t))),
    mul(b, s * t)
  );
};
numeric.Spline.prototype.at = function at(x0) {
  if (typeof x0 === "number") {
    var x = this.x;
    var n = x.length;
    var p,
      q,
      mid,
      floor = Math.floor,
      a,
      b,
      t;
    p = 0;
    q = n - 1;
    while (q - p > 1) {
      mid = floor((p + q) / 2);
      if (x[mid] <= x0) p = mid;
      else q = mid;
    }
    return this._at(x0, p);
  }
  var n = x0.length,
    i,
    ret = Array(n);
  for (i = n - 1; i !== -1; --i) ret[i] = this.at(x0[i]);
  return ret;
};
numeric.Spline.prototype.diff = function diff() {
  var x = this.x;
  var yl = this.yl;
  var yr = this.yr;
  var kl = this.kl;
  var kr = this.kr;
  var n = yl.length;
  var i, dx, dy;
  var zl = kl,
    zr = kr,
    pl = Array(n),
    pr = Array(n);
  var add = numeric.add,
    mul = numeric.mul,
    div = numeric.div,
    sub = numeric.sub;
  for (i = n - 1; i !== -1; --i) {
    dx = x[i + 1] - x[i];
    dy = sub(yr[i + 1], yl[i]);
    pl[i] = div(
      add(mul(dy, 6), mul(kl[i], -4 * dx), mul(kr[i + 1], -2 * dx)),
      dx * dx
    );
    pr[i + 1] = div(
      add(mul(dy, -6), mul(kl[i], 2 * dx), mul(kr[i + 1], 4 * dx)),
      dx * dx
    );
  }
  return new numeric.Spline(x, zl, zr, pl, pr);
};
numeric.Spline.prototype.roots = function roots() {
  function sqr(x) {
    return x * x;
  }
  function heval(y0, y1, k0, k1, x) {
    var A = k0 * 2 - (y1 - y0);
    var B = -k1 * 2 + (y1 - y0);
    var t = (x + 1) * 0.5;
    var s = t * (1 - t);
    return (1 - t) * y0 + t * y1 + A * s * (1 - t) + B * s * t;
  }
  var ret = [];
  var x = this.x,
    yl = this.yl,
    yr = this.yr,
    kl = this.kl,
    kr = this.kr;
  if (typeof yl[0] === "number") {
    yl = [yl];
    yr = [yr];
    kl = [kl];
    kr = [kr];
  }
  var m = yl.length,
    n = x.length - 1,
    i,
    j,
    k,
    y,
    s,
    t;
  var ai,
    bi,
    ci,
    di,
    ret = Array(m),
    ri,
    k0,
    k1,
    y0,
    y1,
    A,
    B,
    D,
    dx,
    cx,
    stops,
    z0,
    z1,
    zm,
    t0,
    t1,
    tm;
  var sqrt = Math.sqrt;
  for (i = 0; i !== m; ++i) {
    ai = yl[i];
    bi = yr[i];
    ci = kl[i];
    di = kr[i];
    ri = [];
    for (j = 0; j !== n; j++) {
      if (j > 0 && bi[j] * ai[j] < 0) ri.push(x[j]);
      dx = x[j + 1] - x[j];
      cx = x[j];
      y0 = ai[j];
      y1 = bi[j + 1];
      k0 = ci[j] / dx;
      k1 = di[j + 1] / dx;
      D = sqr(k0 - k1 + 3 * (y0 - y1)) + 12 * k1 * y0;
      A = k1 + 3 * y0 + 2 * k0 - 3 * y1;
      B = 3 * (k1 + k0 + 2 * (y0 - y1));
      if (D <= 0) {
        z0 = A / B;
        if (z0 > x[j] && z0 < x[j + 1]) stops = [x[j], z0, x[j + 1]];
        else stops = [x[j], x[j + 1]];
      } else {
        z0 = (A - sqrt(D)) / B;
        z1 = (A + sqrt(D)) / B;
        stops = [x[j]];
        if (z0 > x[j] && z0 < x[j + 1]) stops.push(z0);
        if (z1 > x[j] && z1 < x[j + 1]) stops.push(z1);
        stops.push(x[j + 1]);
      }
      t0 = stops[0];
      z0 = this._at(t0, j);
      for (k = 0; k < stops.length - 1; k++) {
        t1 = stops[k + 1];
        z1 = this._at(t1, j);
        if (z0 === 0) {
          ri.push(t0);
          t0 = t1;
          z0 = z1;
          continue;
        }
        if (z1 === 0 || z0 * z1 > 0) {
          t0 = t1;
          z0 = z1;
          continue;
        }
        var side = 0;
        while (1) {
          tm = (z0 * t1 - z1 * t0) / (z0 - z1);
          if (tm <= t0 || tm >= t1) {
            break;
          }
          zm = this._at(tm, j);
          if (zm * z1 > 0) {
            t1 = tm;
            z1 = zm;
            if (side === -1) z0 *= 0.5;
            side = -1;
          } else if (zm * z0 > 0) {
            t0 = tm;
            z0 = zm;
            if (side === 1) z1 *= 0.5;
            side = 1;
          } else break;
        }
        ri.push(tm);
        t0 = stops[k + 1];
        z0 = this._at(t0, j);
      }
      if (z1 === 0) ri.push(t1);
    }
    ret[i] = ri;
  }
  if (typeof this.yl[0] === "number") return ret[0];
  return ret;
};
numeric.spline = function spline(x, y, k1, kn) {
  var n = x.length,
    b = [],
    dx = [],
    dy = [];
  var i;
  var sub = numeric.sub,
    mul = numeric.mul,
    add = numeric.add;
  for (i = n - 2; i >= 0; i--) {
    dx[i] = x[i + 1] - x[i];
    dy[i] = sub(y[i + 1], y[i]);
  }
  if (typeof k1 === "string" || typeof kn === "string") {
    k1 = kn = "periodic";
  }
  // Build sparse tridiagonal system
  var T = [[], [], []];
  switch (typeof k1) {
    case "undefined":
      b[0] = mul(3 / (dx[0] * dx[0]), dy[0]);
      T[0].push(0, 0);
      T[1].push(0, 1);
      T[2].push(2 / dx[0], 1 / dx[0]);
      break;
    case "string":
      b[0] = add(
        mul(3 / (dx[n - 2] * dx[n - 2]), dy[n - 2]),
        mul(3 / (dx[0] * dx[0]), dy[0])
      );
      T[0].push(0, 0, 0);
      T[1].push(n - 2, 0, 1);
      T[2].push(1 / dx[n - 2], 2 / dx[n - 2] + 2 / dx[0], 1 / dx[0]);
      break;
    default:
      b[0] = k1;
      T[0].push(0);
      T[1].push(0);
      T[2].push(1);
      break;
  }
  for (i = 1; i < n - 1; i++) {
    b[i] = add(
      mul(3 / (dx[i - 1] * dx[i - 1]), dy[i - 1]),
      mul(3 / (dx[i] * dx[i]), dy[i])
    );
    T[0].push(i, i, i);
    T[1].push(i - 1, i, i + 1);
    T[2].push(1 / dx[i - 1], 2 / dx[i - 1] + 2 / dx[i], 1 / dx[i]);
  }
  switch (typeof kn) {
    case "undefined":
      b[n - 1] = mul(3 / (dx[n - 2] * dx[n - 2]), dy[n - 2]);
      T[0].push(n - 1, n - 1);
      T[1].push(n - 2, n - 1);
      T[2].push(1 / dx[n - 2], 2 / dx[n - 2]);
      break;
    case "string":
      T[1][T[1].length - 1] = 0;
      break;
    default:
      b[n - 1] = kn;
      T[0].push(n - 1);
      T[1].push(n - 1);
      T[2].push(1);
      break;
  }
  if (typeof b[0] !== "number") b = numeric.transpose(b);
  else b = [b];
  var k = Array(b.length);
  if (typeof k1 === "string") {
    for (i = k.length - 1; i !== -1; --i) {
      k[i] = numeric.ccsLUPSolve(numeric.ccsLUP(numeric.ccsScatter(T)), b[i]);
      k[i][n - 1] = k[i][0];
    }
  } else {
    for (i = k.length - 1; i !== -1; --i) {
      k[i] = numeric.cLUsolve(numeric.cLU(T), b[i]);
    }
  }
  if (typeof y[0] === "number") k = k[0];
  else k = numeric.transpose(k);
  return new numeric.Spline(x, y, y, k, k);
};

// 8. FFT
numeric.fftpow2 = function fftpow2(x, y) {
  var n = x.length;
  if (n === 1) return;
  var cos = Math.cos,
    sin = Math.sin,
    i,
    j;
  var xe = Array(n / 2),
    ye = Array(n / 2),
    xo = Array(n / 2),
    yo = Array(n / 2);
  j = n / 2;
  for (i = n - 1; i !== -1; --i) {
    --j;
    xo[j] = x[i];
    yo[j] = y[i];
    --i;
    xe[j] = x[i];
    ye[j] = y[i];
  }
  fftpow2(xe, ye);
  fftpow2(xo, yo);
  j = n / 2;
  var t,
    k = -6.2831853071795864769252867665590057683943387987502116419 / n,
    ci,
    si;
  for (i = n - 1; i !== -1; --i) {
    --j;
    if (j === -1) j = n / 2 - 1;
    t = k * i;
    ci = cos(t);
    si = sin(t);
    x[i] = xe[j] + ci * xo[j] - si * yo[j];
    y[i] = ye[j] + ci * yo[j] + si * xo[j];
  }
};
numeric._ifftpow2 = function _ifftpow2(x, y) {
  var n = x.length;
  if (n === 1) return;
  var cos = Math.cos,
    sin = Math.sin,
    i,
    j;
  var xe = Array(n / 2),
    ye = Array(n / 2),
    xo = Array(n / 2),
    yo = Array(n / 2);
  j = n / 2;
  for (i = n - 1; i !== -1; --i) {
    --j;
    xo[j] = x[i];
    yo[j] = y[i];
    --i;
    xe[j] = x[i];
    ye[j] = y[i];
  }
  _ifftpow2(xe, ye);
  _ifftpow2(xo, yo);
  j = n / 2;
  var t,
    k = 6.2831853071795864769252867665590057683943387987502116419 / n,
    ci,
    si;
  for (i = n - 1; i !== -1; --i) {
    --j;
    if (j === -1) j = n / 2 - 1;
    t = k * i;
    ci = cos(t);
    si = sin(t);
    x[i] = xe[j] + ci * xo[j] - si * yo[j];
    y[i] = ye[j] + ci * yo[j] + si * xo[j];
  }
};
numeric.ifftpow2 = function ifftpow2(x, y) {
  numeric._ifftpow2(x, y);
  numeric.diveq(x, x.length);
  numeric.diveq(y, y.length);
};
numeric.convpow2 = function convpow2(ax, ay, bx, by) {
  numeric.fftpow2(ax, ay);
  numeric.fftpow2(bx, by);
  var i,
    n = ax.length,
    axi,
    bxi,
    ayi,
    byi;
  for (i = n - 1; i !== -1; --i) {
    axi = ax[i];
    ayi = ay[i];
    bxi = bx[i];
    byi = by[i];
    ax[i] = axi * bxi - ayi * byi;
    ay[i] = axi * byi + ayi * bxi;
  }
  numeric.ifftpow2(ax, ay);
};
numeric.T.prototype.fft = function fft() {
  var x = this.x,
    y = this.y;
  var n = x.length,
    log = Math.log,
    log2 = log(2),
    p = Math.ceil(log(2 * n - 1) / log2),
    m = Math.pow(2, p);
  var cx = numeric.rep([m], 0),
    cy = numeric.rep([m], 0),
    cos = Math.cos,
    sin = Math.sin;
  var k,
    c = -3.14159265358979323846264338327950288419716939937510582 / n,
    t;
  var a = numeric.rep([m], 0),
    b = numeric.rep([m], 0),
    nhalf = Math.floor(n / 2);
  for (k = 0; k < n; k++) a[k] = x[k];
  if (typeof y !== "undefined") for (k = 0; k < n; k++) b[k] = y[k];
  cx[0] = 1;
  for (k = 1; k <= m / 2; k++) {
    t = c * k * k;
    cx[k] = cos(t);
    cy[k] = sin(t);
    cx[m - k] = cos(t);
    cy[m - k] = sin(t);
  }
  var X = new numeric.T(a, b),
    Y = new numeric.T(cx, cy);
  X = X.mul(Y);
  numeric.convpow2(X.x, X.y, numeric.clone(Y.x), numeric.neg(Y.y));
  X = X.mul(Y);
  X.x.length = n;
  X.y.length = n;
  return X;
};
numeric.T.prototype.ifft = function ifft() {
  var x = this.x,
    y = this.y;
  var n = x.length,
    log = Math.log,
    log2 = log(2),
    p = Math.ceil(log(2 * n - 1) / log2),
    m = Math.pow(2, p);
  var cx = numeric.rep([m], 0),
    cy = numeric.rep([m], 0),
    cos = Math.cos,
    sin = Math.sin;
  var k,
    c = 3.14159265358979323846264338327950288419716939937510582 / n,
    t;
  var a = numeric.rep([m], 0),
    b = numeric.rep([m], 0),
    nhalf = Math.floor(n / 2);
  for (k = 0; k < n; k++) a[k] = x[k];
  if (typeof y !== "undefined") for (k = 0; k < n; k++) b[k] = y[k];
  cx[0] = 1;
  for (k = 1; k <= m / 2; k++) {
    t = c * k * k;
    cx[k] = cos(t);
    cy[k] = sin(t);
    cx[m - k] = cos(t);
    cy[m - k] = sin(t);
  }
  var X = new numeric.T(a, b),
    Y = new numeric.T(cx, cy);
  X = X.mul(Y);
  numeric.convpow2(X.x, X.y, numeric.clone(Y.x), numeric.neg(Y.y));
  X = X.mul(Y);
  X.x.length = n;
  X.y.length = n;
  return X.div(n);
};

//9. Unconstrained optimization
numeric.gradient = function gradient(f, x) {
  var n = x.length;
  var f0 = f(x);
  if (isNaN(f0)) throw new Error("gradient: f(x) is a NaN!");
  var max = Math.max;
  var i,
    x0 = numeric.clone(x),
    f1,
    f2,
    J = Array(n);
  var div = numeric.div,
    sub = numeric.sub,
    errest,
    roundoff,
    max = Math.max,
    eps = 1e-3,
    abs = Math.abs,
    min = Math.min;
  var t0,
    t1,
    t2,
    it = 0,
    d1,
    d2,
    N;
  for (i = 0; i < n; i++) {
    var h = max(1e-6 * f0, 1e-8);
    while (1) {
      ++it;
      if (it > 20) {
        throw new Error("Numerical gradient fails");
      }
      x0[i] = x[i] + h;
      f1 = f(x0);
      x0[i] = x[i] - h;
      f2 = f(x0);
      x0[i] = x[i];
      if (isNaN(f1) || isNaN(f2)) {
        h /= 16;
        continue;
      }
      J[i] = (f1 - f2) / (2 * h);
      t0 = x[i] - h;
      t1 = x[i];
      t2 = x[i] + h;
      d1 = (f1 - f0) / h;
      d2 = (f0 - f2) / h;
      N = max(
        abs(J[i]),
        abs(f0),
        abs(f1),
        abs(f2),
        abs(t0),
        abs(t1),
        abs(t2),
        1e-8
      );
      errest = min(
        max(abs(d1 - J[i]), abs(d2 - J[i]), abs(d1 - d2)) / N,
        h / N
      );
      if (errest > eps) {
        h /= 16;
      } else break;
    }
  }
  return J;
};

numeric.uncmin = function uncmin(
  f,
  x0,
  tol,
  gradient,
  maxit,
  callback,
  options
) {
  var grad = numeric.gradient;
  if (typeof options === "undefined") {
    options = {};
  }
  if (typeof tol === "undefined") {
    tol = 1e-8;
  }
  if (typeof gradient === "undefined") {
    gradient = function(x) {
      return grad(f, x);
    };
  }
  if (typeof maxit === "undefined") maxit = 1000;
  x0 = numeric.clone(x0);
  var n = x0.length;
  var f0 = f(x0),
    f1,
    df0;
  if (isNaN(f0)) throw new Error("uncmin: f(x0) is a NaN!");
  var max = Math.max,
    norm2 = numeric.norm2;
  tol = max(tol, numeric.epsilon);
  var step,
    g0,
    g1,
    H1 = options.Hinv || numeric.identity(n);
  var dot = numeric.dot,
    inv = numeric.inv,
    sub = numeric.sub,
    add = numeric.add,
    ten = numeric.tensor,
    div = numeric.div,
    mul = numeric.mul;
  var all = numeric.all,
    isfinite = numeric.isFinite,
    neg = numeric.neg;
  var it = 0,
    i,
    s,
    x1,
    y,
    Hy,
    Hs,
    ys,
    i0,
    t,
    nstep,
    t1,
    t2;
  var msg = "";
  g0 = gradient(x0);
  while (it < maxit) {
    if (typeof callback === "function") {
      if (callback(it, x0, f0, g0, H1)) {
        msg = "Callback returned true";
        break;
      }
    }
    if (!all(isfinite(g0))) {
      msg = "Gradient has Infinity or NaN";
      break;
    }
    step = neg(dot(H1, g0));
    if (!all(isfinite(step))) {
      msg = "Search direction has Infinity or NaN";
      break;
    }
    nstep = norm2(step);
    if (nstep < tol) {
      msg = "Newton step smaller than tol";
      break;
    }
    t = 1;
    df0 = dot(g0, step);
    // line search
    x1 = x0;
    while (it < maxit) {
      if (t * nstep < tol) {
        break;
      }
      s = mul(step, t);
      x1 = add(x0, s);
      f1 = f(x1);
      if (f1 - f0 >= 0.1 * t * df0 || isNaN(f1)) {
        t *= 0.5;
        ++it;
        continue;
      }
      break;
    }
    if (t * nstep < tol) {
      msg = "Line search step size smaller than tol";
      break;
    }
    if (it === maxit) {
      msg = "maxit reached during line search";
      break;
    }
    g1 = gradient(x1);
    y = sub(g1, g0);
    ys = dot(y, s);
    Hy = dot(H1, y);
    H1 = sub(
      add(H1, mul((ys + dot(y, Hy)) / (ys * ys), ten(s, s))),
      div(add(ten(Hy, s), ten(s, Hy)), ys)
    );
    x0 = x1;
    f0 = f1;
    g0 = g1;
    ++it;
  }
  return {
    solution: x0,
    f: f0,
    gradient: g0,
    invHessian: H1,
    iterations: it,
    message: msg
  };
};

// 10. Ode solver (Dormand-Prince)
numeric.Dopri = function Dopri(x, y, f, ymid, iterations, msg, events) {
  this.x = x;
  this.y = y;
  this.f = f;
  this.ymid = ymid;
  this.iterations = iterations;
  this.events = events;
  this.message = msg;
};
numeric.Dopri.prototype._at = function _at(xi, j) {
  function sqr(x) {
    return x * x;
  }
  var sol = this;
  var xs = sol.x;
  var ys = sol.y;
  var k1 = sol.f;
  var ymid = sol.ymid;
  var n = xs.length;
  var x0, x1, xh, y0, y1, yh, xi;
  var floor = Math.floor,
    h;
  var c = 0.5;
  var add = numeric.add,
    mul = numeric.mul,
    sub = numeric.sub,
    p,
    q,
    w;
  x0 = xs[j];
  x1 = xs[j + 1];
  y0 = ys[j];
  y1 = ys[j + 1];
  h = x1 - x0;
  xh = x0 + c * h;
  yh = ymid[j];
  p = sub(k1[j], mul(y0, 1 / (x0 - xh) + 2 / (x0 - x1)));
  q = sub(k1[j + 1], mul(y1, 1 / (x1 - xh) + 2 / (x1 - x0)));
  w = [
    (sqr(xi - x1) * (xi - xh)) / sqr(x0 - x1) / (x0 - xh),
    (sqr(xi - x0) * sqr(xi - x1)) / sqr(x0 - xh) / sqr(x1 - xh),
    (sqr(xi - x0) * (xi - xh)) / sqr(x1 - x0) / (x1 - xh),
    ((xi - x0) * sqr(xi - x1) * (xi - xh)) / sqr(x0 - x1) / (x0 - xh),
    ((xi - x1) * sqr(xi - x0) * (xi - xh)) / sqr(x0 - x1) / (x1 - xh)
  ];
  return add(
    add(add(add(mul(y0, w[0]), mul(yh, w[1])), mul(y1, w[2])), mul(p, w[3])),
    mul(q, w[4])
  );
};
numeric.Dopri.prototype.at = function at(x) {
  var i,
    j,
    k,
    floor = Math.floor;
  if (typeof x !== "number") {
    var n = x.length,
      ret = Array(n);
    for (i = n - 1; i !== -1; --i) {
      ret[i] = this.at(x[i]);
    }
    return ret;
  }
  var x0 = this.x;
  i = 0;
  j = x0.length - 1;
  while (j - i > 1) {
    k = floor(0.5 * (i + j));
    if (x0[k] <= x) i = k;
    else j = k;
  }
  return this._at(x, i);
};

numeric.dopri = function dopri(x0, x1, y0, f, tol, maxit, event) {
  if (typeof tol === "undefined") {
    tol = 1e-6;
  }
  if (typeof maxit === "undefined") {
    maxit = 1000;
  }
  var xs = [x0],
    ys = [y0],
    k1 = [f(x0, y0)],
    k2,
    k3,
    k4,
    k5,
    k6,
    k7,
    ymid = [];
  var A2 = 1 / 5;
  var A3 = [3 / 40, 9 / 40];
  var A4 = [44 / 45, -56 / 15, 32 / 9];
  var A5 = [19372 / 6561, -25360 / 2187, 64448 / 6561, -212 / 729];
  var A6 = [9017 / 3168, -355 / 33, 46732 / 5247, 49 / 176, -5103 / 18656];
  var b = [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84];
  var bm = [
    (0.5 * 6025192743) / 30085553152,
    0,
    (0.5 * 51252292925) / 65400821598,
    (0.5 * -2691868925) / 45128329728,
    (0.5 * 187940372067) / 1594534317056,
    (0.5 * -1776094331) / 19743644256,
    (0.5 * 11237099) / 235043384
  ];
  var c = [1 / 5, 3 / 10, 4 / 5, 8 / 9, 1, 1];
  var e = [
    -71 / 57600,
    0,
    71 / 16695,
    -71 / 1920,
    17253 / 339200,
    -22 / 525,
    1 / 40
  ];
  var i = 0,
    er,
    j;
  var h = (x1 - x0) / 10;
  var it = 0;
  var add = numeric.add,
    mul = numeric.mul,
    y1,
    erinf;
  var max = Math.max,
    min = Math.min,
    abs = Math.abs,
    norminf = numeric.norminf,
    pow = Math.pow;
  var any = numeric.any,
    lt = numeric.lt,
    and = numeric.and,
    sub = numeric.sub;
  var e0, e1, ev;
  var ret = new numeric.Dopri(xs, ys, k1, ymid, -1, "");
  if (typeof event === "function") e0 = event(x0, y0);
  while (x0 < x1 && it < maxit) {
    ++it;
    if (x0 + h > x1) h = x1 - x0;
    k2 = f(x0 + c[0] * h, add(y0, mul(A2 * h, k1[i])));
    k3 = f(
      x0 + c[1] * h,
      add(add(y0, mul(A3[0] * h, k1[i])), mul(A3[1] * h, k2))
    );
    k4 = f(
      x0 + c[2] * h,
      add(
        add(add(y0, mul(A4[0] * h, k1[i])), mul(A4[1] * h, k2)),
        mul(A4[2] * h, k3)
      )
    );
    k5 = f(
      x0 + c[3] * h,
      add(
        add(
          add(add(y0, mul(A5[0] * h, k1[i])), mul(A5[1] * h, k2)),
          mul(A5[2] * h, k3)
        ),
        mul(A5[3] * h, k4)
      )
    );
    k6 = f(
      x0 + c[4] * h,
      add(
        add(
          add(
            add(add(y0, mul(A6[0] * h, k1[i])), mul(A6[1] * h, k2)),
            mul(A6[2] * h, k3)
          ),
          mul(A6[3] * h, k4)
        ),
        mul(A6[4] * h, k5)
      )
    );
    y1 = add(
      add(
        add(
          add(add(y0, mul(k1[i], h * b[0])), mul(k3, h * b[2])),
          mul(k4, h * b[3])
        ),
        mul(k5, h * b[4])
      ),
      mul(k6, h * b[5])
    );
    k7 = f(x0 + h, y1);
    er = add(
      add(
        add(
          add(add(mul(k1[i], h * e[0]), mul(k3, h * e[2])), mul(k4, h * e[3])),
          mul(k5, h * e[4])
        ),
        mul(k6, h * e[5])
      ),
      mul(k7, h * e[6])
    );
    if (typeof er === "number") erinf = abs(er);
    else erinf = norminf(er);
    if (erinf > tol) {
      // reject
      h = 0.2 * h * pow(tol / erinf, 0.25);
      if (x0 + h === x0) {
        ret.msg = "Step size became too small";
        break;
      }
      continue;
    }
    ymid[i] = add(
      add(
        add(
          add(
            add(add(y0, mul(k1[i], h * bm[0])), mul(k3, h * bm[2])),
            mul(k4, h * bm[3])
          ),
          mul(k5, h * bm[4])
        ),
        mul(k6, h * bm[5])
      ),
      mul(k7, h * bm[6])
    );
    ++i;
    xs[i] = x0 + h;
    ys[i] = y1;
    k1[i] = k7;
    if (typeof event === "function") {
      var yi,
        xl = x0,
        xr = x0 + 0.5 * h,
        xi;
      e1 = event(xr, ymid[i - 1]);
      ev = and(lt(e0, 0), lt(0, e1));
      if (!any(ev)) {
        xl = xr;
        xr = x0 + h;
        e0 = e1;
        e1 = event(xr, y1);
        ev = and(lt(e0, 0), lt(0, e1));
      }
      if (any(ev)) {
        var xc, yc, en, ei;
        var side = 0,
          sl = 1.0,
          sr = 1.0;
        while (1) {
          if (typeof e0 === "number")
            xi = (sr * e1 * xl - sl * e0 * xr) / (sr * e1 - sl * e0);
          else {
            xi = xr;
            for (j = e0.length - 1; j !== -1; --j) {
              if (e0[j] < 0 && e1[j] > 0)
                xi = min(
                  xi,
                  (sr * e1[j] * xl - sl * e0[j] * xr) /
                    (sr * e1[j] - sl * e0[j])
                );
            }
          }
          if (xi <= xl || xi >= xr) break;
          yi = ret._at(xi, i - 1);
          ei = event(xi, yi);
          en = and(lt(e0, 0), lt(0, ei));
          if (any(en)) {
            xr = xi;
            e1 = ei;
            ev = en;
            sr = 1.0;
            if (side === -1) sl *= 0.5;
            else sl = 1.0;
            side = -1;
          } else {
            xl = xi;
            e0 = ei;
            sl = 1.0;
            if (side === 1) sr *= 0.5;
            else sr = 1.0;
            side = 1;
          }
        }
        y1 = ret._at(0.5 * (x0 + xi), i - 1);
        ret.f[i] = f(xi, yi);
        ret.x[i] = xi;
        ret.y[i] = yi;
        ret.ymid[i - 1] = y1;
        ret.events = ev;
        ret.iterations = it;
        return ret;
      }
    }
    x0 += h;
    y0 = y1;
    e0 = e1;
    h = min(0.8 * h * pow(tol / erinf, 0.25), 4 * h);
  }
  ret.iterations = it;
  return ret;
};

// 11. Ax = b
numeric.LU = function(A, fast) {
  fast = fast || false;

  var abs = Math.abs;
  var i, j, k, absAjk, Akk, Ak, Pk, Ai;
  var max;
  var n = A.length,
    n1 = n - 1;
  var P = new Array(n);
  if (!fast) A = numeric.clone(A);

  for (k = 0; k < n; ++k) {
    Pk = k;
    Ak = A[k];
    max = abs(Ak[k]);
    for (j = k + 1; j < n; ++j) {
      absAjk = abs(A[j][k]);
      if (max < absAjk) {
        max = absAjk;
        Pk = j;
      }
    }
    P[k] = Pk;

    if (Pk != k) {
      A[k] = A[Pk];
      A[Pk] = Ak;
      Ak = A[k];
    }

    Akk = Ak[k];

    for (i = k + 1; i < n; ++i) {
      A[i][k] /= Akk;
    }

    for (i = k + 1; i < n; ++i) {
      Ai = A[i];
      for (j = k + 1; j < n1; ++j) {
        Ai[j] -= Ai[k] * Ak[j];
        ++j;
        Ai[j] -= Ai[k] * Ak[j];
      }
      if (j === n1) Ai[j] -= Ai[k] * Ak[j];
    }
  }

  return {
    LU: A,
    P: P
  };
};

numeric.LUsolve = function LUsolve(LUP, b) {
  var i, j;
  var LU = LUP.LU;
  var n = LU.length;
  var x = numeric.clone(b);
  var P = LUP.P;
  var Pi, LUi, LUii, tmp;

  for (i = n - 1; i !== -1; --i) x[i] = b[i];
  for (i = 0; i < n; ++i) {
    Pi = P[i];
    if (P[i] !== i) {
      tmp = x[i];
      x[i] = x[Pi];
      x[Pi] = tmp;
    }

    LUi = LU[i];
    for (j = 0; j < i; ++j) {
      x[i] -= x[j] * LUi[j];
    }
  }

  for (i = n - 1; i >= 0; --i) {
    LUi = LU[i];
    for (j = i + 1; j < n; ++j) {
      x[i] -= x[j] * LUi[j];
    }

    x[i] /= LUi[i];
  }

  return x;
};

numeric.solve = function solve(A, b, fast) {
  return numeric.LUsolve(numeric.LU(A, fast), b);
};

// 12. Linear programming
numeric.echelonize = function echelonize(A) {
  var s = numeric.dim(A),
    m = s[0],
    n = s[1];
  var I = numeric.identity(m);
  var P = Array(m);
  var i, j, k, l, Ai, Ii, Z, a;
  var abs = Math.abs;
  var diveq = numeric.diveq;
  A = numeric.clone(A);
  for (i = 0; i < m; ++i) {
    k = 0;
    Ai = A[i];
    Ii = I[i];
    for (j = 1; j < n; ++j) if (abs(Ai[k]) < abs(Ai[j])) k = j;
    P[i] = k;
    diveq(Ii, Ai[k]);
    diveq(Ai, Ai[k]);
    for (j = 0; j < m; ++j)
      if (j !== i) {
        Z = A[j];
        a = Z[k];
        for (l = n - 1; l !== -1; --l) Z[l] -= Ai[l] * a;
        Z = I[j];
        for (l = m - 1; l !== -1; --l) Z[l] -= Ii[l] * a;
      }
  }
  return { I: I, A: A, P: P };
};

numeric.__solveLP = function __solveLP(c, A, b, tol, maxit, x, flag) {
  var sum = numeric.sum,
    log = numeric.log,
    mul = numeric.mul,
    sub = numeric.sub,
    dot = numeric.dot,
    div = numeric.div,
    add = numeric.add;
  var m = c.length,
    n = b.length,
    y;
  var unbounded = false,
    cb,
    i0 = 0;
  var alpha = 1.0;
  var f0,
    df0,
    AT = numeric.transpose(A),
    svd = numeric.svd,
    transpose = numeric.transpose,
    leq = numeric.leq,
    sqrt = Math.sqrt,
    abs = Math.abs;
  var muleq = numeric.muleq;
  var norm = numeric.norminf,
    any = numeric.any,
    min = Math.min;
  var all = numeric.all,
    gt = numeric.gt;
  var p = Array(m),
    A0 = Array(n),
    e = numeric.rep([n], 1),
    H;
  var solve = numeric.solve,
    z = sub(b, dot(A, x)),
    count;
  var dotcc = dot(c, c);
  var g;
  for (count = i0; count < maxit; ++count) {
    var i, j, d;
    for (i = n - 1; i !== -1; --i) A0[i] = div(A[i], z[i]);
    var A1 = transpose(A0);
    for (i = m - 1; i !== -1; --i) p[i] = /*x[i]+*/ sum(A1[i]);
    alpha = 0.25 * abs(dotcc / dot(c, p));
    var a1 = 100 * sqrt(dotcc / dot(p, p));
    if (!isFinite(alpha) || alpha > a1) alpha = a1;
    g = add(c, mul(alpha, p));
    H = dot(A1, A0);
    for (i = m - 1; i !== -1; --i) H[i][i] += 1;
    d = solve(H, div(g, alpha), true);
    var t0 = div(z, dot(A, d));
    var t = 1.0;
    for (i = n - 1; i !== -1; --i) if (t0[i] < 0) t = min(t, -0.999 * t0[i]);
    y = sub(x, mul(d, t));
    z = sub(b, dot(A, y));
    if (!all(gt(z, 0))) return { solution: x, message: "", iterations: count };
    x = y;
    if (alpha < tol) return { solution: y, message: "", iterations: count };
    if (flag) {
      var s = dot(c, g),
        Ag = dot(A, g);
      unbounded = true;
      for (i = n - 1; i !== -1; --i)
        if (s * Ag[i] < 0) {
          unbounded = false;
          break;
        }
    } else {
      if (x[m - 1] >= 0) unbounded = false;
      else unbounded = true;
    }
    if (unbounded)
      return { solution: y, message: "Unbounded", iterations: count };
  }
  return {
    solution: x,
    message: "maximum iteration count exceeded",
    iterations: count
  };
};

numeric._solveLP = function _solveLP(c, A, b, tol, maxit) {
  var m = c.length,
    n = b.length,
    y;
  var sum = numeric.sum,
    log = numeric.log,
    mul = numeric.mul,
    sub = numeric.sub,
    dot = numeric.dot,
    div = numeric.div,
    add = numeric.add;
  var c0 = numeric.rep([m], 0).concat([1]);
  var J = numeric.rep([n, 1], -1);
  var A0 = numeric.blockMatrix([[A, J]]);
  var b0 = b;
  var y = numeric
    .rep([m], 0)
    .concat(Math.max(0, numeric.sup(numeric.neg(b))) + 1);
  var x0 = numeric.__solveLP(c0, A0, b0, tol, maxit, y, false);
  var x = numeric.clone(x0.solution);
  x.length = m;
  var foo = numeric.inf(sub(b, dot(A, x)));
  if (foo < 0) {
    return { solution: NaN, message: "Infeasible", iterations: x0.iterations };
  }
  var ret = numeric.__solveLP(c, A, b, tol, maxit - x0.iterations, x, true);
  ret.iterations += x0.iterations;
  return ret;
};

numeric.solveLP = function solveLP(c, A, b, Aeq, beq, tol, maxit) {
  if (typeof maxit === "undefined") maxit = 1000;
  if (typeof tol === "undefined") tol = numeric.epsilon;
  if (typeof Aeq === "undefined") return numeric._solveLP(c, A, b, tol, maxit);
  var m = Aeq.length,
    n = Aeq[0].length,
    o = A.length;
  var B = numeric.echelonize(Aeq);
  var flags = numeric.rep([n], 0);
  var P = B.P;
  var Q = [];
  var i;
  for (i = P.length - 1; i !== -1; --i) flags[P[i]] = 1;
  for (i = n - 1; i !== -1; --i) if (flags[i] === 0) Q.push(i);
  var g = numeric.getRange;
  var I = numeric.linspace(0, m - 1),
    J = numeric.linspace(0, o - 1);
  var Aeq2 = g(Aeq, I, Q),
    A1 = g(A, J, P),
    A2 = g(A, J, Q),
    dot = numeric.dot,
    sub = numeric.sub;
  var A3 = dot(A1, B.I);
  var A4 = sub(A2, dot(A3, Aeq2)),
    b4 = sub(b, dot(A3, beq));
  var c1 = Array(P.length),
    c2 = Array(Q.length);
  for (i = P.length - 1; i !== -1; --i) c1[i] = c[P[i]];
  for (i = Q.length - 1; i !== -1; --i) c2[i] = c[Q[i]];
  var c4 = sub(c2, dot(c1, dot(B.I, Aeq2)));
  var S = numeric._solveLP(c4, A4, b4, tol, maxit);
  var x2 = S.solution;
  if (x2 !== x2) return S;
  var x1 = dot(B.I, sub(beq, dot(Aeq2, x2)));
  var x = Array(c.length);
  for (i = P.length - 1; i !== -1; --i) x[P[i]] = x1[i];
  for (i = Q.length - 1; i !== -1; --i) x[Q[i]] = x2[i];
  return { solution: x, message: S.message, iterations: S.iterations };
};

numeric.MPStoLP = function MPStoLP(MPS) {
  if (MPS instanceof String) {
    MPS.split("\n");
  }
  var state = 0;
  var states = [
    "Initial state",
    "NAME",
    "ROWS",
    "COLUMNS",
    "RHS",
    "BOUNDS",
    "ENDATA"
  ];
  var n = MPS.length;
  var i,
    j,
    z,
    N = 0,
    rows = {},
    sign = [],
    rl = 0,
    vars = {},
    nv = 0;
  var name;
  var c = [],
    A = [],
    b = [];
  function err(e) {
    throw new Error(
      "MPStoLP: " +
        e +
        "\nLine " +
        i +
        ": " +
        MPS[i] +
        "\nCurrent state: " +
        states[state] +
        "\n"
    );
  }
  for (i = 0; i < n; ++i) {
    z = MPS[i];
    var w0 = z.match(/\S*/g);
    var w = [];
    for (j = 0; j < w0.length; ++j) if (w0[j] !== "") w.push(w0[j]);
    if (w.length === 0) continue;
    for (j = 0; j < states.length; ++j)
      if (z.substr(0, states[j].length) === states[j]) break;
    if (j < states.length) {
      state = j;
      if (j === 1) {
        name = w[1];
      }
      if (j === 6)
        return {
          name: name,
          c: c,
          A: numeric.transpose(A),
          b: b,
          rows: rows,
          vars: vars
        };
      continue;
    }
    switch (state) {
      case 0:
      case 1:
        err("Unexpected line");
      case 2:
        switch (w[0]) {
          case "N":
            if (N === 0) N = w[1];
            else err("Two or more N rows");
            break;
          case "L":
            rows[w[1]] = rl;
            sign[rl] = 1;
            b[rl] = 0;
            ++rl;
            break;
          case "G":
            rows[w[1]] = rl;
            sign[rl] = -1;
            b[rl] = 0;
            ++rl;
            break;
          case "E":
            rows[w[1]] = rl;
            sign[rl] = 0;
            b[rl] = 0;
            ++rl;
            break;
          default:
            err("Parse error " + numeric.prettyPrint(w));
        }
        break;
      case 3:
        if (!vars.hasOwnProperty(w[0])) {
          vars[w[0]] = nv;
          c[nv] = 0;
          A[nv] = numeric.rep([rl], 0);
          ++nv;
        }
        var p = vars[w[0]];
        for (j = 1; j < w.length; j += 2) {
          if (w[j] === N) {
            c[p] = parseFloat(w[j + 1]);
            continue;
          }
          var q = rows[w[j]];
          A[p][q] = (sign[q] < 0 ? -1 : 1) * parseFloat(w[j + 1]);
        }
        break;
      case 4:
        for (j = 1; j < w.length; j += 2)
          b[rows[w[j]]] =
            (sign[rows[w[j]]] < 0 ? -1 : 1) * parseFloat(w[j + 1]);
        break;
      case 5:
        /*FIXME*/ break;
      case 6:
        err("Internal error");
    }
  }
  err("Reached end of file without ENDATA");
};
// seedrandom.js version 2.0.
// Author: David Bau 4/2/2011
//
// Defines a method Math.seedrandom() that, when called, substitutes
// an explicitly seeded RC4-based algorithm for Math.random().  Also
// supports automatic seeding from local or network sources of entropy.
//
// Usage:
//
//   <script src=http://davidbau.com/encode/seedrandom-min.js></script>
//
//   Math.seedrandom('yipee'); Sets Math.random to a function that is
//                             initialized using the given explicit seed.
//
//   Math.seedrandom();        Sets Math.random to a function that is
//                             seeded using the current time, dom state,
//                             and other accumulated local entropy.
//                             The generated seed string is returned.
//
//   Math.seedrandom('yowza', true);
//                             Seeds using the given explicit seed mixed
//                             together with accumulated entropy.
//
//   <script src="http://bit.ly/srandom-512"></script>
//                             Seeds using physical random bits downloaded
//                             from random.org.
//
//   <script src="https://jsonlib.appspot.com/urandom?callback=Math.seedrandom">
//   </script>                 Seeds using urandom bits from call.jsonlib.com,
//                             which is faster than random.org.
//
// Examples:
//
//   Math.seedrandom("hello");            // Use "hello" as the seed.
//   document.write(Math.random());       // Always 0.5463663768140734
//   document.write(Math.random());       // Always 0.43973793770592234
//   var rng1 = Math.random;              // Remember the current prng.
//
//   var autoseed = Math.seedrandom();    // New prng with an automatic seed.
//   document.write(Math.random());       // Pretty much unpredictable.
//
//   Math.random = rng1;                  // Continue "hello" prng sequence.
//   document.write(Math.random());       // Always 0.554769432473455
//
//   Math.seedrandom(autoseed);           // Restart at the previous seed.
//   document.write(Math.random());       // Repeat the 'unpredictable' value.
//
// Notes:
//
// Each time seedrandom('arg') is called, entropy from the passed seed
// is accumulated in a pool to help generate future seeds for the
// zero-argument form of Math.seedrandom, so entropy can be injected over
// time by calling seedrandom with explicit data repeatedly.
//
// On speed - This javascript implementation of Math.random() is about
// 3-10x slower than the built-in Math.random() because it is not native
// code, but this is typically fast enough anyway.  Seeding is more expensive,
// especially if you use auto-seeding.  Some details (timings on Chrome 4):
//
// Our Math.random()            - avg less than 0.002 milliseconds per call
// seedrandom('explicit')       - avg less than 0.5 milliseconds per call
// seedrandom('explicit', true) - avg less than 2 milliseconds per call
// seedrandom()                 - avg about 38 milliseconds per call
//
// LICENSE (BSD):
//
// Copyright 2010 David Bau, all rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//   1. Redistributions of source code must retain the above copyright
//      notice, this list of conditions and the following disclaimer.
//
//   2. Redistributions in binary form must reproduce the above copyright
//      notice, this list of conditions and the following disclaimer in the
//      documentation and/or other materials provided with the distribution.
//
//   3. Neither the name of this module nor the names of its contributors may
//      be used to endorse or promote products derived from this software
//      without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
/**
 * All code is in an anonymous closure to keep the global namespace clean.
 *
 * @param {number=} overflow
 * @param {number=} startdenom
 */

// Patched by Seb so that seedrandom.js does not pollute the Math object.
// My tests suggest that doing Math.trouble = 1 makes Math lookups about 5%
// slower.
numeric.seedrandom = { pow: Math.pow, random: Math.random };

(function(pool, math, width, chunks, significance, overflow, startdenom) {
  //
  // seedrandom()
  // This is the seedrandom function described above.
  //
  math["seedrandom"] = function seedrandom(seed, use_entropy) {
    var key = [];
    var arc4;

    // Flatten the seed string or build one from local entropy if needed.
    seed = mixkey(
      flatten(
        use_entropy
          ? [seed, pool]
          : arguments.length
          ? seed
          : [new Date().getTime(), pool, window],
        3
      ),
      key
    );

    // Use the seed to initialize an ARC4 generator.
    arc4 = new ARC4(key);

    // Mix the randomness into accumulated entropy.
    mixkey(arc4.S, pool);

    // Override Math.random

    // This function returns a random double in [0, 1) that contains
    // randomness in every bit of the mantissa of the IEEE 754 value.

    math["random"] = function random() {
      // Closure to return a random double:
      var n = arc4.g(chunks); // Start with a numerator n < 2 ^ 48
      var d = startdenom; //   and denominator d = 2 ^ 48.
      var x = 0; //   and no 'extra last byte'.
      while (n < significance) {
        // Fill up all significant digits by
        n = (n + x) * width; //   shifting numerator and
        d *= width; //   denominator and generating a
        x = arc4.g(1); //   new least-significant-byte.
      }
      while (n >= overflow) {
        // To avoid rounding up, before adding
        n /= 2; //   last byte, shift everything
        d /= 2; //   right using integer math until
        x >>>= 1; //   we have exactly the desired bits.
      }
      return (n + x) / d; // Form the number within [0, 1).
    };

    // Return the seed that was used
    return seed;
  };

  //
  // ARC4
  //
  // An ARC4 implementation.  The constructor takes a key in the form of
  // an array of at most (width) integers that should be 0 <= x < (width).
  //
  // The g(count) method returns a pseudorandom integer that concatenates
  // the next (count) outputs from ARC4.  Its return value is a number x
  // that is in the range 0 <= x < (width ^ count).
  //
  /** @constructor */
  function ARC4(key) {
    var t,
      u,
      me = this,
      keylen = key.length;
    var i = 0,
      j = (me.i = me.j = me.m = 0);
    me.S = [];
    me.c = [];

    // The empty key [] is treated as [0].
    if (!keylen) {
      key = [keylen++];
    }

    // Set up S using the standard key scheduling algorithm.
    while (i < width) {
      me.S[i] = i++;
    }
    for (i = 0; i < width; i++) {
      t = me.S[i];
      j = lowbits(j + t + key[i % keylen]);
      u = me.S[j];
      me.S[i] = u;
      me.S[j] = t;
    }

    // The "g" method returns the next (count) outputs as one number.
    me.g = function getnext(count) {
      var s = me.S;
      var i = lowbits(me.i + 1);
      var t = s[i];
      var j = lowbits(me.j + t);
      var u = s[j];
      s[i] = u;
      s[j] = t;
      var r = s[lowbits(t + u)];
      while (--count) {
        i = lowbits(i + 1);
        t = s[i];
        j = lowbits(j + t);
        u = s[j];
        s[i] = u;
        s[j] = t;
        r = r * width + s[lowbits(t + u)];
      }
      me.i = i;
      me.j = j;
      return r;
    };
    // For robust unpredictability discard an initial batch of values.
    // See http://www.rsa.com/rsalabs/node.asp?id=2009
    me.g(width);
  }

  //
  // flatten()
  // Converts an object tree to nested arrays of strings.
  //
  /** @param {Object=} result
   * @param {string=} prop
   * @param {string=} typ */
  function flatten(obj, depth, result, prop, typ) {
    result = [];
    typ = typeof obj;
    if (depth && typ == "object") {
      for (prop in obj) {
        if (prop.indexOf("S") < 5) {
          // Avoid FF3 bug (local/sessionStorage)
          try {
            result.push(flatten(obj[prop], depth - 1));
          } catch (e) {}
        }
      }
    }
    return result.length ? result : obj + (typ != "string" ? "\0" : "");
  }

  //
  // mixkey()
  // Mixes a string seed into a key that is an array of integers, and
  // returns a shortened string seed that is equivalent to the result key.
  //
  /** @param {number=} smear
   * @param {number=} j */
  function mixkey(seed, key, smear, j) {
    seed += ""; // Ensure the seed is a string
    smear = 0;
    for (j = 0; j < seed.length; j++) {
      key[lowbits(j)] = lowbits(
        (smear ^= key[lowbits(j)] * 19) + seed.charCodeAt(j)
      );
    }
    seed = "";
    for (j in key) {
      seed += String.fromCharCode(key[j]);
    }
    return seed;
  }

  //
  // lowbits()
  // A quick "n mod width" for width a power of 2.
  //
  function lowbits(n) {
    return n & (width - 1);
  }

  //
  // The following constants are related to IEEE 754 limits.
  //
  startdenom = math.pow(width, chunks);
  significance = math.pow(2, significance);
  overflow = significance * 2;

  //
  // When seedrandom.js is loaded, we immediately mix a few bits
  // from the built-in RNG into the entropy pool.  Because we do
  // not want to intefere with determinstic PRNG state later,
  // seedrandom will not call math.random on its own again after
  // initialization.
  //
  mixkey(math.random(), pool);

  // End anonymous scope, and pass initial values.
})(
  [], // pool: entropy pool starts empty
  numeric.seedrandom, // math: package containing random, pow, and seedrandom
  256, // width: each RC4 output is 0 <= x < 256
  6, // chunks: at least six RC4 outputs for each double
  52 // significance: there are 52 significant digits in a double
);
/* This file is a slightly modified version of quadprog.js from Alberto Santini.
 * It has been slightly modified by Sébastien Loisel to make sure that it handles
 * 0-based Arrays instead of 1-based Arrays.
 * License is in resources/LICENSE.quadprog */
(function(exports) {
  function base0to1(A) {
    if (typeof A !== "object") {
      return A;
    }
    var ret = [],
      i,
      n = A.length;
    for (i = 0; i < n; i++) ret[i + 1] = base0to1(A[i]);
    return ret;
  }
  function base1to0(A) {
    if (typeof A !== "object") {
      return A;
    }
    var ret = [],
      i,
      n = A.length;
    for (i = 1; i < n; i++) ret[i - 1] = base1to0(A[i]);
    return ret;
  }

  function dpori(a, lda, n) {
    var i, j, k, kp1, t;

    for (k = 1; k <= n; k = k + 1) {
      a[k][k] = 1 / a[k][k];
      t = -a[k][k];
      //~ dscal(k - 1, t, a[1][k], 1);
      for (i = 1; i < k; i = i + 1) {
        a[i][k] = t * a[i][k];
      }

      kp1 = k + 1;
      if (n < kp1) {
        break;
      }
      for (j = kp1; j <= n; j = j + 1) {
        t = a[k][j];
        a[k][j] = 0;
        //~ daxpy(k, t, a[1][k], 1, a[1][j], 1);
        for (i = 1; i <= k; i = i + 1) {
          a[i][j] = a[i][j] + t * a[i][k];
        }
      }
    }
  }

  function dposl(a, lda, n, b) {
    var i, k, kb, t;

    for (k = 1; k <= n; k = k + 1) {
      //~ t = ddot(k - 1, a[1][k], 1, b[1], 1);
      t = 0;
      for (i = 1; i < k; i = i + 1) {
        t = t + a[i][k] * b[i];
      }

      b[k] = (b[k] - t) / a[k][k];
    }

    for (kb = 1; kb <= n; kb = kb + 1) {
      k = n + 1 - kb;
      b[k] = b[k] / a[k][k];
      t = -b[k];
      //~ daxpy(k - 1, t, a[1][k], 1, b[1], 1);
      for (i = 1; i < k; i = i + 1) {
        b[i] = b[i] + t * a[i][k];
      }
    }
  }

  function dpofa(a, lda, n, info) {
    var i, j, jm1, k, t, s;

    for (j = 1; j <= n; j = j + 1) {
      info[1] = j;
      s = 0;
      jm1 = j - 1;
      if (jm1 < 1) {
        s = a[j][j] - s;
        if (s <= 0) {
          break;
        }
        a[j][j] = Math.sqrt(s);
      } else {
        for (k = 1; k <= jm1; k = k + 1) {
          //~ t = a[k][j] - ddot(k - 1, a[1][k], 1, a[1][j], 1);
          t = a[k][j];
          for (i = 1; i < k; i = i + 1) {
            t = t - a[i][j] * a[i][k];
          }
          t = t / a[k][k];
          a[k][j] = t;
          s = s + t * t;
        }
        s = a[j][j] - s;
        if (s <= 0) {
          break;
        }
        a[j][j] = Math.sqrt(s);
      }
      info[1] = 0;
    }
  }

  function qpgen2(
    dmat,
    dvec,
    fddmat,
    n,
    sol,
    crval,
    amat,
    bvec,
    fdamat,
    q,
    meq,
    iact,
    nact,
    iter,
    work,
    ierr
  ) {
    var i,
      j,
      l,
      l1,
      info,
      it1,
      iwzv,
      iwrv,
      iwrm,
      iwsv,
      iwuv,
      nvl,
      r,
      iwnbv,
      temp,
      sum,
      t1,
      tt,
      gc,
      gs,
      nu,
      t1inf,
      t2min,
      vsmall,
      tmpa,
      tmpb,
      go;

    r = Math.min(n, q);
    l = 2 * n + (r * (r + 5)) / 2 + 2 * q + 1;

    vsmall = 1.0e-60;
    do {
      vsmall = vsmall + vsmall;
      tmpa = 1 + 0.1 * vsmall;
      tmpb = 1 + 0.2 * vsmall;
    } while (tmpa <= 1 || tmpb <= 1);

    for (i = 1; i <= n; i = i + 1) {
      work[i] = dvec[i];
    }
    for (i = n + 1; i <= l; i = i + 1) {
      work[i] = 0;
    }
    for (i = 1; i <= q; i = i + 1) {
      iact[i] = 0;
    }

    info = [];

    if (ierr[1] === 0) {
      dpofa(dmat, fddmat, n, info);
      if (info[1] !== 0) {
        ierr[1] = 2;
        return;
      }
      dposl(dmat, fddmat, n, dvec);
      dpori(dmat, fddmat, n);
    } else {
      for (j = 1; j <= n; j = j + 1) {
        sol[j] = 0;
        for (i = 1; i <= j; i = i + 1) {
          sol[j] = sol[j] + dmat[i][j] * dvec[i];
        }
      }
      for (j = 1; j <= n; j = j + 1) {
        dvec[j] = 0;
        for (i = j; i <= n; i = i + 1) {
          dvec[j] = dvec[j] + dmat[j][i] * sol[i];
        }
      }
    }

    crval[1] = 0;
    for (j = 1; j <= n; j = j + 1) {
      sol[j] = dvec[j];
      crval[1] = crval[1] + work[j] * sol[j];
      work[j] = 0;
      for (i = j + 1; i <= n; i = i + 1) {
        dmat[i][j] = 0;
      }
    }
    crval[1] = -crval[1] / 2;
    ierr[1] = 0;

    iwzv = n;
    iwrv = iwzv + n;
    iwuv = iwrv + r;
    iwrm = iwuv + r + 1;
    iwsv = iwrm + (r * (r + 1)) / 2;
    iwnbv = iwsv + q;

    for (i = 1; i <= q; i = i + 1) {
      sum = 0;
      for (j = 1; j <= n; j = j + 1) {
        sum = sum + amat[j][i] * amat[j][i];
      }
      work[iwnbv + i] = Math.sqrt(sum);
    }
    nact = 0;
    iter[1] = 0;
    iter[2] = 0;

    function fn_goto_50() {
      iter[1] = iter[1] + 1;

      l = iwsv;
      for (i = 1; i <= q; i = i + 1) {
        l = l + 1;
        sum = -bvec[i];
        for (j = 1; j <= n; j = j + 1) {
          sum = sum + amat[j][i] * sol[j];
        }
        if (Math.abs(sum) < vsmall) {
          sum = 0;
        }
        if (i > meq) {
          work[l] = sum;
        } else {
          work[l] = -Math.abs(sum);
          if (sum > 0) {
            for (j = 1; j <= n; j = j + 1) {
              amat[j][i] = -amat[j][i];
            }
            bvec[i] = -bvec[i];
          }
        }
      }

      for (i = 1; i <= nact; i = i + 1) {
        work[iwsv + iact[i]] = 0;
      }

      nvl = 0;
      temp = 0;
      for (i = 1; i <= q; i = i + 1) {
        if (work[iwsv + i] < temp * work[iwnbv + i]) {
          nvl = i;
          temp = work[iwsv + i] / work[iwnbv + i];
        }
      }
      if (nvl === 0) {
        return 999;
      }

      return 0;
    }

    function fn_goto_55() {
      for (i = 1; i <= n; i = i + 1) {
        sum = 0;
        for (j = 1; j <= n; j = j + 1) {
          sum = sum + dmat[j][i] * amat[j][nvl];
        }
        work[i] = sum;
      }

      l1 = iwzv;
      for (i = 1; i <= n; i = i + 1) {
        work[l1 + i] = 0;
      }
      for (j = nact + 1; j <= n; j = j + 1) {
        for (i = 1; i <= n; i = i + 1) {
          work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j];
        }
      }

      t1inf = true;
      for (i = nact; i >= 1; i = i - 1) {
        sum = work[i];
        l = iwrm + (i * (i + 3)) / 2;
        l1 = l - i;
        for (j = i + 1; j <= nact; j = j + 1) {
          sum = sum - work[l] * work[iwrv + j];
          l = l + j;
        }
        sum = sum / work[l1];
        work[iwrv + i] = sum;
        if (iact[i] < meq) {
          // continue;
          break;
        }
        if (sum < 0) {
          // continue;
          break;
        }
        t1inf = false;
        it1 = i;
      }

      if (!t1inf) {
        t1 = work[iwuv + it1] / work[iwrv + it1];
        for (i = 1; i <= nact; i = i + 1) {
          if (iact[i] < meq) {
            // continue;
            break;
          }
          if (work[iwrv + i] < 0) {
            // continue;
            break;
          }
          temp = work[iwuv + i] / work[iwrv + i];
          if (temp < t1) {
            t1 = temp;
            it1 = i;
          }
        }
      }

      sum = 0;
      for (i = iwzv + 1; i <= iwzv + n; i = i + 1) {
        sum = sum + work[i] * work[i];
      }
      if (Math.abs(sum) <= vsmall) {
        if (t1inf) {
          ierr[1] = 1;
          // GOTO 999
          return 999;
        } else {
          for (i = 1; i <= nact; i = i + 1) {
            work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i];
          }
          work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1;
          // GOTO 700
          return 700;
        }
      } else {
        sum = 0;
        for (i = 1; i <= n; i = i + 1) {
          sum = sum + work[iwzv + i] * amat[i][nvl];
        }
        tt = -work[iwsv + nvl] / sum;
        t2min = true;
        if (!t1inf) {
          if (t1 < tt) {
            tt = t1;
            t2min = false;
          }
        }

        for (i = 1; i <= n; i = i + 1) {
          sol[i] = sol[i] + tt * work[iwzv + i];
          if (Math.abs(sol[i]) < vsmall) {
            sol[i] = 0;
          }
        }

        crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]);
        for (i = 1; i <= nact; i = i + 1) {
          work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i];
        }
        work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt;

        if (t2min) {
          nact = nact + 1;
          iact[nact] = nvl;

          l = iwrm + ((nact - 1) * nact) / 2 + 1;
          for (i = 1; i <= nact - 1; i = i + 1) {
            work[l] = work[i];
            l = l + 1;
          }

          if (nact === n) {
            work[l] = work[n];
          } else {
            for (i = n; i >= nact + 1; i = i - 1) {
              if (work[i] === 0) {
                // continue;
                break;
              }
              gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i]));
              gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i]));
              if (work[i - 1] >= 0) {
                temp = Math.abs(gc * Math.sqrt(1 + (gs * gs) / (gc * gc)));
              } else {
                temp = -Math.abs(gc * Math.sqrt(1 + (gs * gs) / (gc * gc)));
              }
              gc = work[i - 1] / temp;
              gs = work[i] / temp;

              if (gc === 1) {
                // continue;
                break;
              }
              if (gc === 0) {
                work[i - 1] = gs * temp;
                for (j = 1; j <= n; j = j + 1) {
                  temp = dmat[j][i - 1];
                  dmat[j][i - 1] = dmat[j][i];
                  dmat[j][i] = temp;
                }
              } else {
                work[i - 1] = temp;
                nu = gs / (1 + gc);
                for (j = 1; j <= n; j = j + 1) {
                  temp = gc * dmat[j][i - 1] + gs * dmat[j][i];
                  dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i];
                  dmat[j][i - 1] = temp;
                }
              }
            }
            work[l] = work[nact];
          }
        } else {
          sum = -bvec[nvl];
          for (j = 1; j <= n; j = j + 1) {
            sum = sum + sol[j] * amat[j][nvl];
          }
          if (nvl > meq) {
            work[iwsv + nvl] = sum;
          } else {
            work[iwsv + nvl] = -Math.abs(sum);
            if (sum > 0) {
              for (j = 1; j <= n; j = j + 1) {
                amat[j][nvl] = -amat[j][nvl];
              }
              bvec[nvl] = -bvec[nvl];
            }
          }
          // GOTO 700
          return 700;
        }
      }

      return 0;
    }

    function fn_goto_797() {
      l = iwrm + (it1 * (it1 + 1)) / 2 + 1;
      l1 = l + it1;
      if (work[l1] === 0) {
        // GOTO 798
        return 798;
      }
      gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
      gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
      if (work[l1 - 1] >= 0) {
        temp = Math.abs(gc * Math.sqrt(1 + (gs * gs) / (gc * gc)));
      } else {
        temp = -Math.abs(gc * Math.sqrt(1 + (gs * gs) / (gc * gc)));
      }
      gc = work[l1 - 1] / temp;
      gs = work[l1] / temp;

      if (gc === 1) {
        // GOTO 798
        return 798;
      }
      if (gc === 0) {
        for (i = it1 + 1; i <= nact; i = i + 1) {
          temp = work[l1 - 1];
          work[l1 - 1] = work[l1];
          work[l1] = temp;
          l1 = l1 + i;
        }
        for (i = 1; i <= n; i = i + 1) {
          temp = dmat[i][it1];
          dmat[i][it1] = dmat[i][it1 + 1];
          dmat[i][it1 + 1] = temp;
        }
      } else {
        nu = gs / (1 + gc);
        for (i = it1 + 1; i <= nact; i = i + 1) {
          temp = gc * work[l1 - 1] + gs * work[l1];
          work[l1] = nu * (work[l1 - 1] + temp) - work[l1];
          work[l1 - 1] = temp;
          l1 = l1 + i;
        }
        for (i = 1; i <= n; i = i + 1) {
          temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1];
          dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1];
          dmat[i][it1] = temp;
        }
      }

      return 0;
    }

    function fn_goto_798() {
      l1 = l - it1;
      for (i = 1; i <= it1; i = i + 1) {
        work[l1] = work[l];
        l = l + 1;
        l1 = l1 + 1;
      }

      work[iwuv + it1] = work[iwuv + it1 + 1];
      iact[it1] = iact[it1 + 1];
      it1 = it1 + 1;
      if (it1 < nact) {
        // GOTO 797
        return 797;
      }

      return 0;
    }

    function fn_goto_799() {
      work[iwuv + nact] = work[iwuv + nact + 1];
      work[iwuv + nact + 1] = 0;
      iact[nact] = 0;
      nact = nact - 1;
      iter[2] = iter[2] + 1;

      return 0;
    }

    go = 0;
    while (true) {
      go = fn_goto_50();
      if (go === 999) {
        return;
      }
      while (true) {
        go = fn_goto_55();
        if (go === 0) {
          break;
        }
        if (go === 999) {
          return;
        }
        if (go === 700) {
          if (it1 === nact) {
            fn_goto_799();
          } else {
            while (true) {
              fn_goto_797();
              go = fn_goto_798();
              if (go !== 797) {
                break;
              }
            }
            fn_goto_799();
          }
        }
      }
    }
  }

  function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) {
    Dmat = base0to1(Dmat);
    dvec = base0to1(dvec);
    Amat = base0to1(Amat);
    var i,
      n,
      q,
      nact,
      r,
      crval = [],
      iact = [],
      sol = [],
      work = [],
      iter = [],
      message;

    meq = meq || 0;
    factorized = factorized ? base0to1(factorized) : [undefined, 0];
    bvec = bvec ? base0to1(bvec) : [];

    // In Fortran the array index starts from 1
    n = Dmat.length - 1;
    q = Amat[1].length - 1;

    if (!bvec) {
      for (i = 1; i <= q; i = i + 1) {
        bvec[i] = 0;
      }
    }
    for (i = 1; i <= q; i = i + 1) {
      iact[i] = 0;
    }
    nact = 0;
    r = Math.min(n, q);
    for (i = 1; i <= n; i = i + 1) {
      sol[i] = 0;
    }
    crval[1] = 0;
    for (i = 1; i <= 2 * n + (r * (r + 5)) / 2 + 2 * q + 1; i = i + 1) {
      work[i] = 0;
    }
    for (i = 1; i <= 2; i = i + 1) {
      iter[i] = 0;
    }

    qpgen2(
      Dmat,
      dvec,
      n,
      n,
      sol,
      crval,
      Amat,
      bvec,
      n,
      q,
      meq,
      iact,
      nact,
      iter,
      work,
      factorized
    );

    message = "";
    if (factorized[1] === 1) {
      message = "constraints are inconsistent, no solution!";
    }
    if (factorized[1] === 2) {
      message = "matrix D in quadratic function is not positive definite!";
    }

    return {
      solution: base1to0(sol),
      value: base1to0(crval),
      unconstrained_solution: base1to0(dvec),
      iterations: base1to0(iter),
      iact: base1to0(iact),
      message: message
    };
  }
  exports.solveQP = solveQP;
})(numeric);
/*
Shanti Rao sent me this routine by private email. I had to modify it
slightly to work on Arrays instead of using a Matrix object.
It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py
*/

numeric.svd = function svd(A) {
  var temp;
  //Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970)
  var prec = numeric.epsilon; //Math.pow(2,-52) // assumes double prec
  var tolerance = 1e-64 / prec;
  var itmax = 50;
  var c = 0;
  var i = 0;
  var j = 0;
  var k = 0;
  var l = 0;

  var u = numeric.clone(A);
  var m = u.length;

  var n = u[0].length;

  if (m < n) throw "Need more rows than columns";

  var e = new Array(n);
  var q = new Array(n);
  for (i = 0; i < n; i++) e[i] = q[i] = 0.0;
  var v = numeric.rep([n, n], 0);
  //	v.zero();

  function pythag(a, b) {
    a = Math.abs(a);
    b = Math.abs(b);
    if (a > b) return a * Math.sqrt(1.0 + (b * b) / a / a);
    else if (b == 0.0) return a;
    return b * Math.sqrt(1.0 + (a * a) / b / b);
  }

  //Householder's reduction to bidiagonal form

  var f = 0.0;
  var g = 0.0;
  var h = 0.0;
  var x = 0.0;
  var y = 0.0;
  var z = 0.0;
  var s = 0.0;

  for (i = 0; i < n; i++) {
    e[i] = g;
    s = 0.0;
    l = i + 1;
    for (j = i; j < m; j++) s += u[j][i] * u[j][i];
    if (s <= tolerance) g = 0.0;
    else {
      f = u[i][i];
      g = Math.sqrt(s);
      if (f >= 0.0) g = -g;
      h = f * g - s;
      u[i][i] = f - g;
      for (j = l; j < n; j++) {
        s = 0.0;
        for (k = i; k < m; k++) s += u[k][i] * u[k][j];
        f = s / h;
        for (k = i; k < m; k++) u[k][j] += f * u[k][i];
      }
    }
    q[i] = g;
    s = 0.0;
    for (j = l; j < n; j++) s = s + u[i][j] * u[i][j];
    if (s <= tolerance) g = 0.0;
    else {
      f = u[i][i + 1];
      g = Math.sqrt(s);
      if (f >= 0.0) g = -g;
      h = f * g - s;
      u[i][i + 1] = f - g;
      for (j = l; j < n; j++) e[j] = u[i][j] / h;
      for (j = l; j < m; j++) {
        s = 0.0;
        for (k = l; k < n; k++) s += u[j][k] * u[i][k];
        for (k = l; k < n; k++) u[j][k] += s * e[k];
      }
    }
    y = Math.abs(q[i]) + Math.abs(e[i]);
    if (y > x) x = y;
  }

  // accumulation of right hand gtransformations
  for (i = n - 1; i != -1; i += -1) {
    if (g != 0.0) {
      h = g * u[i][i + 1];
      for (j = l; j < n; j++) v[j][i] = u[i][j] / h;
      for (j = l; j < n; j++) {
        s = 0.0;
        for (k = l; k < n; k++) s += u[i][k] * v[k][j];
        for (k = l; k < n; k++) v[k][j] += s * v[k][i];
      }
    }
    for (j = l; j < n; j++) {
      v[i][j] = 0;
      v[j][i] = 0;
    }
    v[i][i] = 1;
    g = e[i];
    l = i;
  }

  // accumulation of left hand transformations
  for (i = n - 1; i != -1; i += -1) {
    l = i + 1;
    g = q[i];
    for (j = l; j < n; j++) u[i][j] = 0;
    if (g != 0.0) {
      h = u[i][i] * g;
      for (j = l; j < n; j++) {
        s = 0.0;
        for (k = l; k < m; k++) s += u[k][i] * u[k][j];
        f = s / h;
        for (k = i; k < m; k++) u[k][j] += f * u[k][i];
      }
      for (j = i; j < m; j++) u[j][i] = u[j][i] / g;
    } else for (j = i; j < m; j++) u[j][i] = 0;
    u[i][i] += 1;
  }

  // diagonalization of the bidiagonal form
  prec = prec * x;
  for (k = n - 1; k != -1; k += -1) {
    for (var iteration = 0; iteration < itmax; iteration++) {
      // test f splitting
      var test_convergence = false;
      for (l = k; l != -1; l += -1) {
        if (Math.abs(e[l]) <= prec) {
          test_convergence = true;
          break;
        }
        if (Math.abs(q[l - 1]) <= prec) break;
      }
      if (!test_convergence) {
        // cancellation of e[l] if l>0
        c = 0.0;
        s = 1.0;
        var l1 = l - 1;
        for (i = l; i < k + 1; i++) {
          f = s * e[i];
          e[i] = c * e[i];
          if (Math.abs(f) <= prec) break;
          g = q[i];
          h = pythag(f, g);
          q[i] = h;
          c = g / h;
          s = -f / h;
          for (j = 0; j < m; j++) {
            y = u[j][l1];
            z = u[j][i];
            u[j][l1] = y * c + z * s;
            u[j][i] = -y * s + z * c;
          }
        }
      }
      // test f convergence
      z = q[k];
      if (l == k) {
        //convergence
        if (z < 0.0) {
          //q[k] is made non-negative
          q[k] = -z;
          for (j = 0; j < n; j++) v[j][k] = -v[j][k];
        }
        break; //break out of iteration loop and move on to next k value
      }
      if (iteration >= itmax - 1) throw "Error: no convergence.";
      // shift from bottom 2x2 minor
      x = q[l];
      y = q[k - 1];
      g = e[k - 1];
      h = e[k];
      f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
      g = pythag(f, 1.0);
      if (f < 0.0) f = ((x - z) * (x + z) + h * (y / (f - g) - h)) / x;
      else f = ((x - z) * (x + z) + h * (y / (f + g) - h)) / x;
      // next QR transformation
      c = 1.0;
      s = 1.0;
      for (i = l + 1; i < k + 1; i++) {
        g = e[i];
        y = q[i];
        h = s * g;
        g = c * g;
        z = pythag(f, h);
        e[i - 1] = z;
        c = f / z;
        s = h / z;
        f = x * c + g * s;
        g = -x * s + g * c;
        h = y * s;
        y = y * c;
        for (j = 0; j < n; j++) {
          x = v[j][i - 1];
          z = v[j][i];
          v[j][i - 1] = x * c + z * s;
          v[j][i] = -x * s + z * c;
        }
        z = pythag(f, h);
        q[i - 1] = z;
        c = f / z;
        s = h / z;
        f = c * g + s * y;
        x = -s * g + c * y;
        for (j = 0; j < m; j++) {
          y = u[j][i - 1];
          z = u[j][i];
          u[j][i - 1] = y * c + z * s;
          u[j][i] = -y * s + z * c;
        }
      }
      e[l] = 0.0;
      e[k] = f;
      q[k] = x;
    }
  }

  //vt= transpose(v)
  //return (u,q,vt)
  for (i = 0; i < q.length; i++) if (q[i] < prec) q[i] = 0;

  //sort eigenvalues
  for (i = 0; i < n; i++) {
    //writeln(q)
    for (j = i - 1; j >= 0; j--) {
      if (q[j] < q[i]) {
        //  writeln(i,'-',j)
        c = q[j];
        q[j] = q[i];
        q[i] = c;
        for (k = 0; k < u.length; k++) {
          temp = u[k][i];
          u[k][i] = u[k][j];
          u[k][j] = temp;
        }
        for (k = 0; k < v.length; k++) {
          temp = v[k][i];
          v[k][i] = v[k][j];
          v[k][j] = temp;
        }
        //	   u.swapCols(i,j)
        //	   v.swapCols(i,j)
        i = j;
      }
    }
  }

  return { U: u, S: q, V: v };
};