(function(){science = {version: "1.7.0"}; // semver

science.ascending = function(a, b) {

  return a - b;

};

// Euler's constant.

science.EULER = .5772156649015329;

// Compute exp(x) - 1 accurately for small x.

science.expm1 = function(x) {

  return (x < 1e-5 && x > -1e-5) ? x + .5 * x * x : Math.exp(x) - 1;

};

science.functor = function(v) {

  return typeof v === "function" ? v : function() { return v; };

};

// Based on:

// http://www.johndcook.com/blog/2010/06/02/whats-so-hard-about-finding-a-hypotenuse/

science.hypot = function(x, y) {

  x = Math.abs(x);

  y = Math.abs(y);

  var max,

      min;

  if (x > y) { max = x; min = y; }

  else       { max = y; min = x; }

  var r = min / max;

  return max * Math.sqrt(1 + r * r);

};

science.quadratic = function() {

  var complex = false;

  function quadratic(a, b, c) {

    var d = b * b - 4 * a * c;

    if (d > 0) {

      d = Math.sqrt(d) / (2 * a);

      return complex

        ? [{r: -b - d, i: 0}, {r: -b + d, i: 0}]

        : [-b - d, -b + d];

    } else if (d === 0) {

      d = -b / (2 * a);

      return complex ? [{r: d, i: 0}] : [d];

    } else {

      if (complex) {

        d = Math.sqrt(-d) / (2 * a);

        return [

          {r: -b, i: -d},

          {r: -b, i: d}

        ];

      }

      return [];

    }

  }

  quadratic.complex = function(x) {

    if (!arguments.length) return complex;

    complex = x;

    return quadratic;

  };

  return quadratic;

};

// Constructs a multi-dimensional array filled with zeroes.

science.zeroes = function(n) {

  var i = -1,

      a = [];

  if (arguments.length === 1)

    while (++i < n)

      a[i] = 0;

  else

    while (++i < n)

      a[i] = science.zeroes.apply(

        this, Array.prototype.slice.call(arguments, 1));

  return a;

};

science.vector = {};

science.vector.cross = function(a, b) {

  // TODO how to handle non-3D vectors?

  // TODO handle 7D vectors?

  return [

    a[1] * b[2] - a[2] * b[1],

    a[2] * b[0] - a[0] * b[2],

    a[0] * b[1] - a[1] * b[0]

  ];

};

science.vector.dot = function(a, b) {

  var s = 0,

      i = -1,

      n = Math.min(a.length, b.length);

  while (++i < n) s += a[i] * b[i];

  return s;

};

science.vector.length = function(p) {

  return Math.sqrt(science.vector.dot(p, p));

};

science.vector.normalize = function(p) {

  var length = science.vector.length(p);

  return p.map(function(d) { return d / length; });

};

// 4x4 matrix determinant.

science.vector.determinant = function(matrix) {

  var m = matrix[0].concat(matrix[1]).concat(matrix[2]).concat(matrix[3]);

  return (

    m[12] * m[9]  * m[6]  * m[3]  - m[8] * m[13] * m[6]  * m[3]  -

    m[12] * m[5]  * m[10] * m[3]  + m[4] * m[13] * m[10] * m[3]  +

    m[8]  * m[5]  * m[14] * m[3]  - m[4] * m[9]  * m[14] * m[3]  -

    m[12] * m[9]  * m[2]  * m[7]  + m[8] * m[13] * m[2]  * m[7]  +

    m[12] * m[1]  * m[10] * m[7]  - m[0] * m[13] * m[10] * m[7]  -

    m[8]  * m[1]  * m[14] * m[7]  + m[0] * m[9]  * m[14] * m[7]  +

    m[12] * m[5]  * m[2]  * m[11] - m[4] * m[13] * m[2]  * m[11] -

    m[12] * m[1]  * m[6]  * m[11] + m[0] * m[13] * m[6]  * m[11] +

    m[4]  * m[1]  * m[14] * m[11] - m[0] * m[5]  * m[14] * m[11] -

    m[8]  * m[5]  * m[2]  * m[15] + m[4] * m[9]  * m[2]  * m[15] +

    m[8]  * m[1]  * m[6]  * m[15] - m[0] * m[9]  * m[6]  * m[15] -

    m[4]  * m[1]  * m[10] * m[15] + m[0] * m[5]  * m[10] * m[15]);

};

// Performs in-place Gauss-Jordan elimination.

//

// Based on Jarno Elonen's Python version (public domain):

// http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html

science.vector.gaussjordan = function(m, eps) {

  if (!eps) eps = 1e-10;

  var h = m.length,

      w = m[0].length,

      y = -1,

      y2,

      x;

  while (++y < h) {

    var maxrow = y;

    // Find max pivot.

    y2 = y; while (++y2 < h) {

      if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y]))

        maxrow = y2;

    }

    // Swap.

    var tmp = m[y];

    m[y] = m[maxrow];

    m[maxrow] = tmp;

    // Singular?

    if (Math.abs(m[y][y]) <= eps) return false;

    // Eliminate column y.

    y2 = y; while (++y2 < h) {

      var c = m[y2][y] / m[y][y];

      x = y - 1; while (++x < w) {

        m[y2][x] -= m[y][x] * c;

      }

    }

  }

  // Backsubstitute.

  y = h; while (--y >= 0) {

    var c = m[y][y];

    y2 = -1; while (++y2 < y) {

      x = w; while (--x >= y) {

        m[y2][x] -=  m[y][x] * m[y2][y] / c;

      }

    }

    m[y][y] /= c;

    // Normalize row y.

    x = h - 1; while (++x < w) {

      m[y][x] /= c;

    }

  }

  return true;

};

// Find matrix inverse using Gauss-Jordan.

science.vector.inverse = function(m) {

  var n = m.length

      i = -1;

  // Check if the matrix is square.

  if (n !== m[0].length) return;

  // Augment with identity matrix I to get AI.

  m = m.map(function(row, i) {

    var identity = new Array(n),

        j = -1;

    while (++j < n) identity[j] = i === j ? 1 : 0;

    return row.concat(identity);

  });

  // Compute IA^-1.

  science.vector.gaussjordan(m);

  // Remove identity matrix I to get A^-1.

  while (++i < n) {

    m[i] = m[i].slice(n);

  }

  return m;

};

science.vector.multiply = function(a, b) {

  var m = a.length,

      n = b[0].length,

      p = b.length,

      i = -1,

      j,

      k;

  if (p !== a[0].length) throw {"error": "columns(a) != rows(b); " + a[0].length + " != " + p};

  var ab = new Array(m);

  while (++i < m) {

    ab[i] = new Array(n);

    j = -1; while(++j < n) {

      var s = 0;

      k = -1; while (++k < p) s += a[i][k] * b[k][j];

      ab[i][j] = s;

    }

  }

  return ab;

};

science.vector.transpose = function(a) {

  var m = a.length,

      n = a[0].length,

      i = -1,

      j,

      b = new Array(n);

  while (++i < n) {

    b[i] = new Array(m);

    j = -1; while (++j < m) b[i][j] = a[j][i];

  }

  return b;

};

})()