|
1 (function(){science = {version: "1.7.0"}; // semver |
|
2 science.ascending = function(a, b) { |
|
3 return a - b; |
|
4 }; |
|
5 // Euler's constant. |
|
6 science.EULER = .5772156649015329; |
|
7 // Compute exp(x) - 1 accurately for small x. |
|
8 science.expm1 = function(x) { |
|
9 return (x < 1e-5 && x > -1e-5) ? x + .5 * x * x : Math.exp(x) - 1; |
|
10 }; |
|
11 science.functor = function(v) { |
|
12 return typeof v === "function" ? v : function() { return v; }; |
|
13 }; |
|
14 // Based on: |
|
15 // http://www.johndcook.com/blog/2010/06/02/whats-so-hard-about-finding-a-hypotenuse/ |
|
16 science.hypot = function(x, y) { |
|
17 x = Math.abs(x); |
|
18 y = Math.abs(y); |
|
19 var max, |
|
20 min; |
|
21 if (x > y) { max = x; min = y; } |
|
22 else { max = y; min = x; } |
|
23 var r = min / max; |
|
24 return max * Math.sqrt(1 + r * r); |
|
25 }; |
|
26 science.quadratic = function() { |
|
27 var complex = false; |
|
28 |
|
29 function quadratic(a, b, c) { |
|
30 var d = b * b - 4 * a * c; |
|
31 if (d > 0) { |
|
32 d = Math.sqrt(d) / (2 * a); |
|
33 return complex |
|
34 ? [{r: -b - d, i: 0}, {r: -b + d, i: 0}] |
|
35 : [-b - d, -b + d]; |
|
36 } else if (d === 0) { |
|
37 d = -b / (2 * a); |
|
38 return complex ? [{r: d, i: 0}] : [d]; |
|
39 } else { |
|
40 if (complex) { |
|
41 d = Math.sqrt(-d) / (2 * a); |
|
42 return [ |
|
43 {r: -b, i: -d}, |
|
44 {r: -b, i: d} |
|
45 ]; |
|
46 } |
|
47 return []; |
|
48 } |
|
49 } |
|
50 |
|
51 quadratic.complex = function(x) { |
|
52 if (!arguments.length) return complex; |
|
53 complex = x; |
|
54 return quadratic; |
|
55 }; |
|
56 |
|
57 return quadratic; |
|
58 }; |
|
59 // Constructs a multi-dimensional array filled with zeroes. |
|
60 science.zeroes = function(n) { |
|
61 var i = -1, |
|
62 a = []; |
|
63 if (arguments.length === 1) |
|
64 while (++i < n) |
|
65 a[i] = 0; |
|
66 else |
|
67 while (++i < n) |
|
68 a[i] = science.zeroes.apply( |
|
69 this, Array.prototype.slice.call(arguments, 1)); |
|
70 return a; |
|
71 }; |
|
72 science.vector = {}; |
|
73 science.vector.cross = function(a, b) { |
|
74 // TODO how to handle non-3D vectors? |
|
75 // TODO handle 7D vectors? |
|
76 return [ |
|
77 a[1] * b[2] - a[2] * b[1], |
|
78 a[2] * b[0] - a[0] * b[2], |
|
79 a[0] * b[1] - a[1] * b[0] |
|
80 ]; |
|
81 }; |
|
82 science.vector.dot = function(a, b) { |
|
83 var s = 0, |
|
84 i = -1, |
|
85 n = Math.min(a.length, b.length); |
|
86 while (++i < n) s += a[i] * b[i]; |
|
87 return s; |
|
88 }; |
|
89 science.vector.length = function(p) { |
|
90 return Math.sqrt(science.vector.dot(p, p)); |
|
91 }; |
|
92 science.vector.normalize = function(p) { |
|
93 var length = science.vector.length(p); |
|
94 return p.map(function(d) { return d / length; }); |
|
95 }; |
|
96 // 4x4 matrix determinant. |
|
97 science.vector.determinant = function(matrix) { |
|
98 var m = matrix[0].concat(matrix[1]).concat(matrix[2]).concat(matrix[3]); |
|
99 return ( |
|
100 m[12] * m[9] * m[6] * m[3] - m[8] * m[13] * m[6] * m[3] - |
|
101 m[12] * m[5] * m[10] * m[3] + m[4] * m[13] * m[10] * m[3] + |
|
102 m[8] * m[5] * m[14] * m[3] - m[4] * m[9] * m[14] * m[3] - |
|
103 m[12] * m[9] * m[2] * m[7] + m[8] * m[13] * m[2] * m[7] + |
|
104 m[12] * m[1] * m[10] * m[7] - m[0] * m[13] * m[10] * m[7] - |
|
105 m[8] * m[1] * m[14] * m[7] + m[0] * m[9] * m[14] * m[7] + |
|
106 m[12] * m[5] * m[2] * m[11] - m[4] * m[13] * m[2] * m[11] - |
|
107 m[12] * m[1] * m[6] * m[11] + m[0] * m[13] * m[6] * m[11] + |
|
108 m[4] * m[1] * m[14] * m[11] - m[0] * m[5] * m[14] * m[11] - |
|
109 m[8] * m[5] * m[2] * m[15] + m[4] * m[9] * m[2] * m[15] + |
|
110 m[8] * m[1] * m[6] * m[15] - m[0] * m[9] * m[6] * m[15] - |
|
111 m[4] * m[1] * m[10] * m[15] + m[0] * m[5] * m[10] * m[15]); |
|
112 }; |
|
113 // Performs in-place Gauss-Jordan elimination. |
|
114 // |
|
115 // Based on Jarno Elonen's Python version (public domain): |
|
116 // http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html |
|
117 science.vector.gaussjordan = function(m, eps) { |
|
118 if (!eps) eps = 1e-10; |
|
119 |
|
120 var h = m.length, |
|
121 w = m[0].length, |
|
122 y = -1, |
|
123 y2, |
|
124 x; |
|
125 |
|
126 while (++y < h) { |
|
127 var maxrow = y; |
|
128 |
|
129 // Find max pivot. |
|
130 y2 = y; while (++y2 < h) { |
|
131 if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y])) |
|
132 maxrow = y2; |
|
133 } |
|
134 |
|
135 // Swap. |
|
136 var tmp = m[y]; |
|
137 m[y] = m[maxrow]; |
|
138 m[maxrow] = tmp; |
|
139 |
|
140 // Singular? |
|
141 if (Math.abs(m[y][y]) <= eps) return false; |
|
142 |
|
143 // Eliminate column y. |
|
144 y2 = y; while (++y2 < h) { |
|
145 var c = m[y2][y] / m[y][y]; |
|
146 x = y - 1; while (++x < w) { |
|
147 m[y2][x] -= m[y][x] * c; |
|
148 } |
|
149 } |
|
150 } |
|
151 |
|
152 // Backsubstitute. |
|
153 y = h; while (--y >= 0) { |
|
154 var c = m[y][y]; |
|
155 y2 = -1; while (++y2 < y) { |
|
156 x = w; while (--x >= y) { |
|
157 m[y2][x] -= m[y][x] * m[y2][y] / c; |
|
158 } |
|
159 } |
|
160 m[y][y] /= c; |
|
161 // Normalize row y. |
|
162 x = h - 1; while (++x < w) { |
|
163 m[y][x] /= c; |
|
164 } |
|
165 } |
|
166 return true; |
|
167 }; |
|
168 // Find matrix inverse using Gauss-Jordan. |
|
169 science.vector.inverse = function(m) { |
|
170 var n = m.length |
|
171 i = -1; |
|
172 |
|
173 // Check if the matrix is square. |
|
174 if (n !== m[0].length) return; |
|
175 |
|
176 // Augment with identity matrix I to get AI. |
|
177 m = m.map(function(row, i) { |
|
178 var identity = new Array(n), |
|
179 j = -1; |
|
180 while (++j < n) identity[j] = i === j ? 1 : 0; |
|
181 return row.concat(identity); |
|
182 }); |
|
183 |
|
184 // Compute IA^-1. |
|
185 science.vector.gaussjordan(m); |
|
186 |
|
187 // Remove identity matrix I to get A^-1. |
|
188 while (++i < n) { |
|
189 m[i] = m[i].slice(n); |
|
190 } |
|
191 |
|
192 return m; |
|
193 }; |
|
194 science.vector.multiply = function(a, b) { |
|
195 var m = a.length, |
|
196 n = b[0].length, |
|
197 p = b.length, |
|
198 i = -1, |
|
199 j, |
|
200 k; |
|
201 if (p !== a[0].length) throw {"error": "columns(a) != rows(b); " + a[0].length + " != " + p}; |
|
202 var ab = new Array(m); |
|
203 while (++i < m) { |
|
204 ab[i] = new Array(n); |
|
205 j = -1; while(++j < n) { |
|
206 var s = 0; |
|
207 k = -1; while (++k < p) s += a[i][k] * b[k][j]; |
|
208 ab[i][j] = s; |
|
209 } |
|
210 } |
|
211 return ab; |
|
212 }; |
|
213 science.vector.transpose = function(a) { |
|
214 var m = a.length, |
|
215 n = a[0].length, |
|
216 i = -1, |
|
217 j, |
|
218 b = new Array(n); |
|
219 while (++i < n) { |
|
220 b[i] = new Array(m); |
|
221 j = -1; while (++j < m) b[i][j] = a[j][i]; |
|
222 } |
|
223 return b; |
|
224 }; |
|
225 })() |