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