2011-09-10 02:37:59 +04:00
|
|
|
(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;
|
|
|
|
};
|
2011-07-06 21:57:14 +04:00
|
|
|
science.functor = function(v) {
|
|
|
|
return typeof v === "function" ? v : function() { return v; };
|
|
|
|
};
|
2011-09-10 02:37:59 +04:00
|
|
|
// 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;
|
|
|
|
};
|
2011-07-06 21:57:14 +04:00
|
|
|
})()
|