1 /*
  2     Copyright 2008-2015
  3         Matthias Ehmann,
  4         Michael Gerhaeuser,
  5         Carsten Miller,
  6         Bianca Valentin,
  7         Alfred Wassermann,
  8         Peter Wilfahrt
  9 
 10     This file is part of JSXGraph.
 11 
 12     JSXGraph is free software dual licensed under the GNU LGPL or MIT License.
 13 
 14     You can redistribute it and/or modify it under the terms of the
 15 
 16       * GNU Lesser General Public License as published by
 17         the Free Software Foundation, either version 3 of the License, or
 18         (at your option) any later version
 19       OR
 20       * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT
 21 
 22     JSXGraph is distributed in the hope that it will be useful,
 23     but WITHOUT ANY WARRANTY; without even the implied warranty of
 24     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 25     GNU Lesser General Public License for more details.
 26 
 27     You should have received a copy of the GNU Lesser General Public License and
 28     the MIT License along with JSXGraph. If not, see <http://www.gnu.org/licenses/>
 29     and <http://opensource.org/licenses/MIT/>.
 30  */
 31 
 32 
 33 /*global JXG: true, define: true, Float32Array: true */
 34 /*jslint nomen: true, plusplus: true, bitwise: true*/
 35 
 36 /* depends:
 37  jxg
 38  */
 39 
 40 /**
 41  * @fileoverview In this file the namespace JXG.Math is defined, which is the base namespace
 42  * for namespaces like Math.Numerics, Math.Algebra, Math.Statistics etc.
 43  */
 44 
 45 define(['jxg', 'utils/type'], function (JXG, Type) {
 46 
 47     "use strict";
 48 
 49     var undef,
 50 
 51         /*
 52          * Dynamic programming approach for recursive functions.
 53          * From "Speed up your JavaScript, Part 3" by Nicholas C. Zakas.
 54          * @see JXG.Math.factorial
 55          * @see JXG.Math.binomial
 56          * http://blog.thejit.org/2008/09/05/memoization-in-javascript/
 57          *
 58          * This method is hidden, because it is only used in JXG.Math. If someone wants
 59          * to use it in JSXGraph outside of JXG.Math, it should be moved to jsxgraph.js
 60          */
 61         memoizer = function (f) {
 62             var cache, join;
 63 
 64             if (f.memo) {
 65                 return f.memo;
 66             }
 67 
 68             cache = {};
 69             join = Array.prototype.join;
 70 
 71             f.memo = function () {
 72                 var key = join.call(arguments);
 73 
 74                 // Seems to be a bit faster than "if (a in b)"
 75                 return (cache[key] !== undef) ?
 76                         cache[key] :
 77                         cache[key] = f.apply(this, arguments);
 78             };
 79 
 80             return f.memo;
 81         };
 82 
 83     /**
 84      * Math namespace.
 85      * @namespace
 86      */
 87     JXG.Math = {
 88         /**
 89          * eps defines the closeness to zero. If the absolute value of a given number is smaller
 90          * than eps, it is considered to be equal to zero.
 91          * @type number
 92          */
 93         eps: 0.000001,
 94 
 95         /**
 96          * The JavaScript implementation of the % operator returns the symmetric modulo.
 97          * They are both identical if a >= 0 and m >= 0 but the results differ if a or m < 0.
 98          * @param {Number} a
 99          * @param {Number} m
100          * @returns {Number} Mathematical modulo <tt>a mod m</tt>
101          */
102         mod: function (a, m) {
103             return a - Math.floor(a / m) * m;
104         },
105 
106         /**
107          * Initializes a vector as an array with the coefficients set to the given value resp. zero.
108          * @param {Number} n Length of the vector
109          * @param {Number} [init=0] Initial value for each coefficient
110          * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a
111          * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows.
112          */
113         vector: function (n, init) {
114             var r, i;
115 
116             init = init || 0;
117             r = [];
118 
119             for (i = 0; i < n; i++) {
120                 r[i] = init;
121             }
122 
123             return r;
124         },
125 
126         /**
127          * Initializes a matrix as an array of rows with the given value.
128          * @param {Number} n Number of rows
129          * @param {Number} [m=n] Number of columns
130          * @param {Number} [init=0] Initial value for each coefficient
131          * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a
132          * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows.
133          */
134         matrix: function (n, m, init) {
135             var r, i, j;
136 
137             init = init || 0;
138             m = m || n;
139             r = [];
140 
141             for (i = 0; i < n; i++) {
142                 r[i] = [];
143 
144                 for (j = 0; j < m; j++) {
145                     r[i][j] = init;
146                 }
147             }
148 
149             return r;
150         },
151 
152         /**
153          * Generates an identity matrix. If n is a number and m is undefined or not a number, a square matrix is generated,
154          * if n and m are both numbers, an nxm matrix is generated.
155          * @param {Number} n Number of rows
156          * @param {Number} [m=n] Number of columns
157          * @returns {Array} A square matrix of length <tt>n</tt> with all coefficients equal to 0 except a_(i,i), i out of (1, ..., n), if <tt>m</tt> is undefined or not a number
158          * or a <tt>n</tt> times <tt>m</tt>-matrix with a_(i,j) = 0 and a_(i,i) = 1 if m is a number.
159          */
160         identity: function (n, m) {
161             var r, i;
162 
163             if ((m === undef) && (typeof m !== 'number')) {
164                 m = n;
165             }
166 
167             r = this.matrix(n, m);
168 
169             for (i = 0; i < Math.min(n, m); i++) {
170                 r[i][i] = 1;
171             }
172 
173             return r;
174         },
175 
176         /**
177          * Generates a 4x4 matrix for 3D to 2D projections.
178          * @param {Number} l Left
179          * @param {Number} r Right
180          * @param {Number} t Top
181          * @param {Number} b Bottom
182          * @param {Number} n Near
183          * @param {Number} f Far
184          * @returns {Array} 4x4 Matrix
185          */
186         frustum: function (l, r, b, t, n, f) {
187             var ret = this.matrix(4, 4);
188 
189             ret[0][0] = (n * 2) / (r - l);
190             ret[0][1] = 0;
191             ret[0][2] = (r + l) / (r - l);
192             ret[0][3] = 0;
193 
194             ret[1][0] = 0;
195             ret[1][1] = (n * 2) / (t - b);
196             ret[1][2] = (t + b) / (t - b);
197             ret[1][3] = 0;
198 
199             ret[2][0] = 0;
200             ret[2][1] = 0;
201             ret[2][2] = -(f + n) / (f - n);
202             ret[2][3] = -(f * n * 2) / (f - n);
203 
204             ret[3][0] = 0;
205             ret[3][1] = 0;
206             ret[3][2] = -1;
207             ret[3][3] = 0;
208 
209             return ret;
210         },
211 
212         /**
213          * Generates a 4x4 matrix for 3D to 2D projections.
214          * @param {Number} fov Field of view in vertical direction, given in rad.
215          * @param {Number} ratio Aspect ratio of the projection plane.
216          * @param {Number} n Near
217          * @param {Number} f Far
218          * @returns {Array} 4x4 Projection Matrix
219          */
220         projection: function (fov, ratio, n, f) {
221             var t = n * Math.tan(fov / 2),
222                 r = t * ratio;
223 
224             return this.frustum(-r, r, -t, t, n, f);
225         },
226 
227         /**
228          * Multiplies a vector vec to a matrix mat: mat * vec. The matrix is interpreted by this function as an array of rows. Please note: This
229          * function does not check if the dimensions match.
230          * @param {Array} mat Two dimensional array of numbers. The inner arrays describe the columns, the outer ones the matrix' rows.
231          * @param {Array} vec Array of numbers
232          * @returns {Array} Array of numbers containing the result
233          * @example
234          * var A = [[2, 1],
235          *          [1, 3]],
236          *     b = [4, 5],
237          *     c;
238          * c = JXG.Math.matVecMult(A, b)
239          * // c === [13, 19];
240          */
241         matVecMult: function (mat, vec) {
242             var i, s, k,
243                 m = mat.length,
244                 n = vec.length,
245                 res = [];
246 
247             if (n === 3) {
248                 for (i = 0; i < m; i++) {
249                     res[i] = mat[i][0] * vec[0] + mat[i][1] * vec[1] + mat[i][2] * vec[2];
250                 }
251             } else {
252                 for (i = 0; i < m; i++) {
253                     s = 0;
254                     for (k = 0; k < n; k++) {
255                         s += mat[i][k] * vec[k];
256                     }
257                     res[i] = s;
258                 }
259             }
260             return res;
261         },
262 
263         /**
264          * Computes the product of the two matrices mat1*mat2.
265          * @param {Array} mat1 Two dimensional array of numbers
266          * @param {Array} mat2 Two dimensional array of numbers
267          * @returns {Array} Two dimensional Array of numbers containing result
268          */
269         matMatMult: function (mat1, mat2) {
270             var i, j, s, k,
271                 m = mat1.length,
272                 n = m > 0 ? mat2[0].length : 0,
273                 m2 = mat2.length,
274                 res = this.matrix(m, n);
275 
276             for (i = 0; i < m; i++) {
277                 for (j = 0; j < n; j++) {
278                     s = 0;
279                     for (k = 0; k < m2; k++) {
280                         s += mat1[i][k] * mat2[k][j];
281                     }
282                     res[i][j] = s;
283                 }
284             }
285             return res;
286         },
287 
288         /**
289          * Transposes a matrix given as a two dimensional array.
290          * @param {Array} M The matrix to be transposed
291          * @returns {Array} The transpose of M
292          */
293         transpose: function (M) {
294             var MT, i, j,
295                 m, n;
296 
297             // number of rows of M
298             m = M.length;
299             // number of columns of M
300             n = M.length > 0 ? M[0].length : 0;
301             MT = this.matrix(n, m);
302 
303             for (i = 0; i < n; i++) {
304                 for (j = 0; j < m; j++) {
305                     MT[i][j] = M[j][i];
306                 }
307             }
308 
309             return MT;
310         },
311 
312         /**
313          * Compute the inverse of an nxn matrix with Gauss elimination.
314          * @param {Array} Ain
315          * @returns {Array} Inverse matrix of Ain
316          */
317         inverse: function (Ain) {
318             var i, j, k, s, ma, r, swp,
319                 n = Ain.length,
320                 A = [],
321                 p = [],
322                 hv = [];
323 
324             for (i = 0; i < n; i++) {
325                 A[i] = [];
326                 for (j = 0; j < n; j++) {
327                     A[i][j] = Ain[i][j];
328                 }
329                 p[i] = i;
330             }
331 
332             for (j = 0; j < n; j++) {
333                 // pivot search:
334                 ma = Math.abs(A[j][j]);
335                 r = j;
336 
337                 for (i = j + 1; i < n; i++) {
338                     if (Math.abs(A[i][j]) > ma) {
339                         ma = Math.abs(A[i][j]);
340                         r = i;
341                     }
342                 }
343 
344                 // Singular matrix
345                 if (ma <= this.eps) {
346                     return [];
347                 }
348 
349                 // swap rows:
350                 if (r > j) {
351                     for (k = 0; k < n; k++) {
352                         swp = A[j][k];
353                         A[j][k] = A[r][k];
354                         A[r][k] = swp;
355                     }
356 
357                     swp = p[j];
358                     p[j] = p[r];
359                     p[r] = swp;
360                 }
361 
362                 // transformation:
363                 s = 1.0 / A[j][j];
364                 for (i = 0; i < n; i++) {
365                     A[i][j] *= s;
366                 }
367                 A[j][j] = s;
368 
369                 for (k = 0; k < n; k++) {
370                     if (k !== j) {
371                         for (i = 0; i < n; i++) {
372                             if (i !== j) {
373                                 A[i][k] -= A[i][j] * A[j][k];
374                             }
375                         }
376                         A[j][k] = -s * A[j][k];
377                     }
378                 }
379             }
380 
381             // swap columns:
382             for (i = 0; i < n; i++) {
383                 for (k = 0; k < n; k++) {
384                     hv[p[k]] = A[i][k];
385                 }
386                 for (k = 0; k < n; k++) {
387                     A[i][k] = hv[k];
388                 }
389             }
390 
391             return A;
392         },
393 
394         /**
395          * Inner product of two vectors a and b. n is the length of the vectors.
396          * @param {Array} a Vector
397          * @param {Array} b Vector
398          * @param {Number} [n] Length of the Vectors. If not given the length of the first vector is taken.
399          * @returns {Number} The inner product of a and b.
400          */
401         innerProduct: function (a, b, n) {
402             var i,
403                 s = 0;
404 
405             if ((n === undef) || (typeof n !== 'number')) {
406                 n = a.length;
407             }
408 
409             for (i = 0; i < n; i++) {
410                 s += a[i] * b[i];
411             }
412 
413             return s;
414         },
415 
416         /**
417          * Calculates the cross product of two vectors both of length three.
418          * In case of homogeneous coordinates this is either
419          * <ul>
420          * <li>the intersection of two lines</li>
421          * <li>the line through two points</li>
422          * </ul>
423          * @param {Array} c1 Homogeneous coordinates of line or point 1
424          * @param {Array} c2 Homogeneous coordinates of line or point 2
425          * @returns {Array} vector of length 3: homogeneous coordinates of the resulting point / line.
426          */
427         crossProduct: function (c1, c2) {
428             return [c1[1] * c2[2] - c1[2] * c2[1],
429                 c1[2] * c2[0] - c1[0] * c2[2],
430                 c1[0] * c2[1] - c1[1] * c2[0]];
431         },
432 
433         /**
434          * Compute the factorial of a positive integer. If a non-integer value
435          * is given, the fraction will be ignored.
436          * @function
437          * @param {Number} n
438          * @returns {Number} n! = n*(n-1)*...*2*1
439          */
440         factorial: memoizer(function (n) {
441             if (n < 0) {
442                 return NaN;
443             }
444 
445             n = Math.floor(n);
446 
447             if (n === 0 || n === 1) {
448                 return 1;
449             }
450 
451             return n * this.factorial(n - 1);
452         }),
453 
454         /**
455          * Computes the binomial coefficient n over k.
456          * @function
457          * @param {Number} n Fraction will be ignored
458          * @param {Number} k Fraction will be ignored
459          * @returns {Number} The binomial coefficient n over k
460          */
461         binomial: memoizer(function (n, k) {
462             var b, i;
463 
464             if (k > n || k < 0) {
465                 return NaN;
466             }
467 
468             k = Math.round(k);
469             n = Math.round(n);
470 
471             if (k === 0 || k === n) {
472                 return 1;
473             }
474 
475             b = 1;
476 
477             for (i = 0; i < k; i++) {
478                 b *= (n - i);
479                 b /= (i + 1);
480             }
481 
482             return b;
483         }),
484 
485         /**
486          * Calculates the cosine hyperbolicus of x.
487          * @param {Number} x The number the cosine hyperbolicus will be calculated of.
488          * @returns {Number} Cosine hyperbolicus of the given value.
489          */
490         cosh: function (x) {
491             return (Math.exp(x) + Math.exp(-x)) * 0.5;
492         },
493 
494         /**
495          * Sine hyperbolicus of x.
496          * @param {Number} x The number the sine hyperbolicus will be calculated of.
497          * @returns {Number} Sine hyperbolicus of the given value.
498          */
499         sinh: function (x) {
500             return (Math.exp(x) - Math.exp(-x)) * 0.5;
501         },
502 
503         /**
504          * Compute base to the power of exponent.
505          * @param {Number} base
506          * @param {Number} exponent
507          * @returns {Number} base to the power of exponent.
508          */
509         pow: function (base, exponent) {
510             if (base === 0) {
511                 if (exponent === 0) {
512                     return 1;
513                 }
514 
515                 return 0;
516             }
517 
518             if (Math.floor(exponent) === exponent) {
519                 // a is an integer
520                 return Math.pow(base, exponent);
521             }
522 
523             // a is not an integer
524             if (base > 0) {
525                 return Math.exp(exponent * Math.log(Math.abs(base)));
526             }
527 
528             return NaN;
529         },
530 
531         /**
532          * Logarithm to base 10.
533          * @param {Number} x
534          * @returns {Number} log10(x) Logarithm of x to base 10.
535          */
536         log10: function (x) {
537             return Math.log(x) / Math.log(10.0);
538         },
539 
540         /**
541          * Logarithm to base 2.
542          * @param {Number} x
543          * @returns {Number} log2(x) Logarithm of x to base 2.
544          */
545         log2: function (x) {
546             return Math.log(x) / Math.log(2.0);
547         },
548 
549         /**
550          * Logarithm to arbitrary base b. If b is not given, natural log is taken, i.e. b = e.
551          * @param {Number} x
552          * @param {Number} b base
553          * @returns {Number} log(x, b) Logarithm of x to base b, that is log(x)/log(b).
554          */
555         log: function (x, b) {
556             if (typeof b !== 'undefined' && Type.isNumber(b)) {
557                 return Math.log(x) / Math.log(b);
558             } else {
559                 return Math.log(x);
560             }
561         },
562 
563         /**
564          * A square & multiply algorithm to compute base to the power of exponent.
565          * Implementated by Wolfgang Riedl.
566          * @param {Number} base
567          * @param {Number} exponent
568          * @returns {Number} Base to the power of exponent
569          */
570         squampow: function (base, exponent) {
571             var result;
572 
573             if (Math.floor(exponent) === exponent) {
574                 // exponent is integer (could be zero)
575                 result = 1;
576 
577                 if (exponent < 0) {
578                     // invert: base
579                     base = 1.0 / base;
580                     exponent *= -1;
581                 }
582 
583                 while (exponent !== 0) {
584                     if (exponent & 1) {
585                         result *= base;
586                     }
587 
588                     exponent >>= 1;
589                     base *= base;
590                 }
591                 return result;
592             }
593 
594             return this.pow(base, exponent);
595         },
596 
597         /**
598          * Normalize the standard form [c, b0, b1, a, k, r, q0, q1].
599          * @private
600          * @param {Array} stdform The standard form to be normalized.
601          * @returns {Array} The normalized standard form.
602          */
603         normalize: function (stdform) {
604             var n, signr,
605                 a2 = 2 * stdform[3],
606                 r = stdform[4] / a2;
607 
608             stdform[5] = r;
609             stdform[6] = -stdform[1] / a2;
610             stdform[7] = -stdform[2] / a2;
611 
612             if (!isFinite(r)) {
613                 n = Math.sqrt(stdform[1] * stdform[1] + stdform[2] * stdform[2]);
614 
615                 stdform[0] /= n;
616                 stdform[1] /= n;
617                 stdform[2] /= n;
618                 stdform[3] = 0;
619                 stdform[4] = 1;
620             } else if (Math.abs(r) >= 1) {
621                 stdform[0] = (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) / (2 * r);
622                 stdform[1] = -stdform[6] / r;
623                 stdform[2] = -stdform[7] / r;
624                 stdform[3] = 1 / (2 * r);
625                 stdform[4] = 1;
626             } else {
627                 signr = (r <= 0 ? -1 : 1);
628                 stdform[0] = signr * (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) * 0.5;
629                 stdform[1] = -signr * stdform[6];
630                 stdform[2] = -signr * stdform[7];
631                 stdform[3] = signr / 2;
632                 stdform[4] = signr * r;
633             }
634 
635             return stdform;
636         },
637 
638         /**
639          * Converts a two dimensional array to a one dimensional Float32Array that can be processed by WebGL.
640          * @param {Array} m A matrix in a two dimensional array.
641          * @returns {Float32Array} A one dimensional array containing the matrix in column wise notation. Provides a fall
642          * back to the default JavaScript Array if Float32Array is not available.
643          */
644         toGL: function (m) {
645             var v, i, j;
646 
647             if (typeof Float32Array === 'function') {
648                 v = new Float32Array(16);
649             } else {
650                 v = new Array(16);
651             }
652 
653             if (m.length !== 4 && m[0].length !== 4) {
654                 return v;
655             }
656 
657             for (i = 0; i < 4; i++) {
658                 for (j = 0; j < 4; j++) {
659                     v[i + 4 * j] = m[i][j];
660                 }
661             }
662 
663             return v;
664         }
665     };
666 
667     return JXG.Math;
668 });
669