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          * Initializes a vector as an array with the coefficients set to the given value resp. zero.
 77          * @param {Number} n Length of the vector
 78          * @param {Number} [init=0] Initial value for each coefficient
 79          * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a
 80          * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows.
 81          */
 82         vector: function(n, init) {
 83             var r, i;
 84 
 85             init = init || 0;
 86 
 87             r = new Array(Math.ceil(n));
 88             for(i=0; i<n; i++) {
 89                 r[i] = init;
 90             }
 91 
 92             return r;
 93         },
 94 
 95         /**
 96          * Initializes a matrix as an array of rows with the given value.
 97          * @param {Number} n Number of rows
 98          * @param {Number} [m=n] Number of columns
 99          * @param {Number} [init=0] Initial value for each coefficient
100          * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a
101          * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows.
102          */
103         matrix: function(n, m, init) {
104             var r, i, j;
105 
106             init = init || 0;
107             m = m || n;
108 
109             r = new Array(Math.ceil(n));
110             for(i=0; i<n; i++) {
111                 r[i] = new Array(Math.ceil(m));
112                 for(j=0; j<m; j++) {
113                     r[i][j] = init;
114                 }
115             }
116 
117             return r;
118         },
119 
120         /**
121          * Generates an identity matrix. If n is a number and m is undefined or not a number, a square matrix is generated,
122          * if n and m are both numbers, an nxm matrix is generated.
123          * @param {Number} n Number of rows
124          * @param {Number} [m=n] Number of columns
125          * @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
126          * 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.
127          */
128         identity: function(n, m) {
129             var r, i;
130 
131             if((m === undefined) && (typeof m !== 'number')) {
132                 m = n;
133             }
134 
135             r = this.matrix(n, m);
136             for(i=0; i<Math.min(n, m); i++) {
137                 r[i][i] = 1;
138             }
139 
140             return r;
141         },
142 
143         /**
144          * 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
145          * function does not check if the dimensions match.
146          * @param {Array} mat Two dimensional array of numbers. The inner arrays describe the columns, the outer ones the matrix' rows.
147          * @param {Array} vec Array of numbers
148          * @returns {Array} Array of numbers containing the result
149          * @example
150          * var A = [[2, 1],
151          *          [1, 3]],
152          *     b = [4, 5],
153          *     c;
154          * c = JXG.Math.matVecMult(A, b)
155          * // c === [13, 19];
156          */
157         matVecMult: function(mat, vec) {
158             var m = mat.length,
159                     n = vec.length,
160                     res = [],
161                     i, s, k;
162 
163             if (n===3) {
164                 for (i=0;i<m;i++) {
165                     res[i] = mat[i][0]*vec[0] + mat[i][1]*vec[1] + mat[i][2]*vec[2];
166                 }
167             } else {
168                 for (i=0;i<m;i++) {
169                     s = 0;
170                     for (k=0;k<n;k++) { s += mat[i][k]*vec[k]; }
171                     res[i] = s;
172                 }
173             }
174             return res;
175         },
176 
177         /**
178          * Computes the product of the two matrices mat1*mat2.
179          * @param {Array} mat1 Two dimensional array of numbers
180          * @param {Array} mat2 Two dimensional array of numbers
181          * @returns {Array} Two dimensional Array of numbers containing result
182          */
183         matMatMult: function(mat1, mat2) {
184             var m = mat1.length,
185                     n = m>0 ? mat2[0].length : 0,
186                     m2 = mat2.length,
187                     res = this.matrix(m,n),
188                     i, j, s, k;
189   
190             for (i=0;i<m;i++) {
191                 for (j=0;j<n;j++) {
192                     s = 0;
193                     for (k=0;k<m2;k++) {
194                         s += mat1[i][k]*mat2[k][j];
195                     }
196                     res[i][j] = s;
197                 }
198             }
199             return res;
200         },
201 
202         /**
203          * Transposes a matrix given as a two dimensional array.
204          * @param {Array} M The matrix to be transposed
205          * @returns {Array} The transpose of M
206          */
207         transpose: function(M) {
208             var MT, i, j,
209                     m, n;
210 
211             m = M.length;                     // number of rows of M
212             n = M.length>0 ? M[0].length : 0; // number of columns of M
213             MT = this.matrix(n,m);
214 
215             for (i=0; i<n; i++) {
216                 for (j=0;j<m;j++) {
217                     MT[i][j] = M[j][i];
218                 }
219             }
220             return MT;
221         },
222 
223         /**
224          * Compute the inverse of an nxn matrix with Gauss elimination.
225          * @param {Array} Ain
226          * @returns {Array} Inverse matrix of Ain
227          */
228         inverse: function(Ain) {
229             var i,j,k,s,ma,r,swp,
230                 n = Ain.length,
231                 A = [],
232                 p = [],
233                 hv = [];
234 
235             for (i=0;i<n;i++) {
236                 A[i] = [];
237                 for (j=0;j<n;j++) { A[i][j] = Ain[i][j]; }
238                 p[i] = i;
239             }
240 
241             for (j=0;j<n;j++) {
242                 // pivot search:
243                 ma = Math.abs(A[j][j]);
244                 r = j;
245                 for (i=j+1;i<n;i++) {
246                     if (Math.abs(A[i][j])>ma) {
247                         ma = Math.abs(A[i][j]);
248                         r = i;
249                     }
250                 }
251                 if (ma<=JXG.Math.eps) { // Singular matrix
252                     return false;
253                 }
254                 // swap rows:
255                 if (r>j) {
256                     for (k=0;k<n;k++) {
257                         swp = A[j][k]; A[j][k] = A[r][k]; A[r][k] = swp;
258                     }
259                     swp = p[j]; p[j] = p[r]; p[r] = swp;
260                 }
261                 // transformation:
262                 s = 1.0/A[j][j];
263                 for (i=0;i<n;i++) {
264                     A[i][j] *= s;
265                 }
266                 A[j][j] = s;
267                 for (k=0;k<n;k++) if (k!=j) {
268                     for (i=0;i<n;i++) if (i!=j) {
269                         A[i][k] -= A[i][j]*A[j][k];
270                     }
271                     A[j][k] = -s*A[j][k];
272                 }
273             }
274             // swap columns:
275             for (i=0;i<n;i++) {
276                 for (k=0;k<n;k++) { hv[p[k]] = A[i][k]; }
277                 for (k=0;k<n;k++) { A[i][k] = hv[k]; }
278             }
279             return A;
280         },
281 
282         /**
283          * Inner product of two vectors a and b. n is the length of the vectors.
284          * @param {Array} a Vector
285          * @param {Array} b Vector
286          * @param {Number} [n] Length of the Vectors. If not given the length of the first vector is taken.
287          * @return The inner product of a and b.
288          */
289         innerProduct: function(a, b, n) {
290             var i, s = 0;
291 
292             if((n === undefined) || (typeof n !== 'number')) {
293                 n = a.length;
294             }
295 
296             for (i=0; i<n; i++) {
297                 s += a[i]*b[i];
298             }
299 
300             return s;
301         },
302 
303         /**
304          * Calculates the cross product of two vectors both of length three.
305          * In case of homogeneous coordinates this is either
306          * <ul>
307          * <li>the intersection of two lines</li>
308          * <li>the line through two points</li>
309          * </ul>
310          * @param {Array} c1 Homogeneous coordinates of line or point 1
311          * @param {Array} c2 Homogeneous coordinates of line or point 2
312          * @returns {Array} vector of length 3: homogeneous coordinates of the resulting point / line.
313          */
314         crossProduct: function(c1,c2) {
315             return [c1[1]*c2[2]-c1[2]*c2[1],
316                 c1[2]*c2[0]-c1[0]*c2[2],
317                 c1[0]*c2[1]-c1[1]*c2[0]];
318         },
319 
320         /**
321          * Compute the factorial of a positive integer. If a non-integer value
322          * is given, the fraction will be ignored.
323          * @function
324          * @param {Number} n
325          * @returns {Number} n! = n*(n-1)*...*2*1
326          */
327         factorial: memoizer(function (n) {
328             if (n<0) return NaN;
329             n = Math.floor(n);
330             if (n===0 || n===1) return 1;
331             return n*arguments.callee(n-1);
332         }),
333 
334         /**
335          * Computes the binomial coefficient n over k.
336          * @function
337          * @param {Number} n Fraction will be ignored
338          * @param {Number} k Fraction will be ignored
339          * @returns {Number} The binomial coefficient n over k
340          */
341         binomial: memoizer(function(n,k) {
342             var b, i;
343 
344             if (k>n || k<0) return NaN;
345 
346             k = Math.floor(k);
347             n = Math.floor(n);
348 
349             if (k===0 || k===n) return 1;
350 
351             b = 1;
352             for (i=0; i<k; i++) {
353                 b *= (n-i);
354                 b /= (i+1);
355             }
356             return b;
357         }),
358 
359         /**
360          * Calculates the cosine hyperbolicus of x.
361          * @param {Number} x The number the cosine hyperbolicus will be calculated of.
362          * @returns {Number} Cosine hyperbolicus of the given value.
363          */
364         cosh: function(x) {
365             return (Math.exp(x)+Math.exp(-x))*0.5;
366         },
367 
368         /**
369          * Sine hyperbolicus of x.
370          * @param {Number} x The number the sine hyperbolicus will be calculated of.
371          * @returns {Number} Sine hyperbolicus of the given value.
372          */
373         sinh: function(x) {
374             return (Math.exp(x)-Math.exp(-x))*0.5;
375         },
376 
377         /**
378          * Compute base to the power of exponent.
379          * @param {Number} base
380          * @param {Number} exponent
381          * @returns {Number} base to the power of exponent.
382          */
383         pow: function(base, exponent) {
384             if (base===0) {
385                 if (exponent===0) {
386                     return 1;
387                 } else {
388                     return 0;
389                 }
390             }
391 
392             if (Math.floor(exponent)===exponent) {// a is an integer
393                 return Math.pow(base,exponent);
394             } else { // a is not an integer
395                 if (base>0) {
396                     return Math.exp(exponent*Math.log(Math.abs(base)));
397                 } else {
398                     return NaN;
399                 }
400             }
401         },
402 
403         /**
404          * A square & multiply algorithm to compute base to the power of exponent.
405          * Implementated by Wolfgang Riedl.
406          * @param {Number} base
407          * @param {Number} exponent
408          * @returns {Number} Base to the power of exponent
409          */
410         squampow: function(base, exponent) {
411             var result;
412 
413             if (Math.floor(exponent)===exponent) { // exponent is integer (could be zero)
414                 result = 1;
415                 if(exponent < 0) {
416                     // invert: base
417                     base = 1.0 / base;
418                     exponent *= -1;
419                 }
420                 while(exponent != 0) {
421                     if(exponent & 1)
422                         result *= base;
423                     exponent >>= 1;
424                     base *= base;
425                 }
426                 return result;
427             } else {
428                 return this.pow(base, exponent);
429             }
430         },
431 
432         /**
433          * Normalize the standard form [c, b0, b1, a, k, r, q0, q1].
434          * @private
435          * @param {Array} stdform The standard form to be normalized.
436          * @returns {Array} The normalized standard form.
437          */
438         normalize: function(stdform) {
439             var a2 = 2*stdform[3],
440                     r = stdform[4]/(a2),  // k/(2a)
441                     n, signr;
442             stdform[5] = r;
443             stdform[6] = -stdform[1]/a2;
444             stdform[7] = -stdform[2]/a2;
445             if (r===Infinity || isNaN(r)) {
446                 n = Math.sqrt(stdform[1]*stdform[1]+stdform[2]*stdform[2]);
447                 stdform[0] /= n;
448                 stdform[1] /= n;
449                 stdform[2] /= n;
450                 stdform[3] = 0;
451                 stdform[4] = 1;
452             } else if (Math.abs(r)>=1) {
453                 stdform[0] = (stdform[6]*stdform[6]+stdform[7]*stdform[7]-r*r)/(2*r);
454                 stdform[1] = -stdform[6]/r;
455                 stdform[2] = -stdform[7]/r;
456                 stdform[3] = 1/(2*r);
457                 stdform[4] = 1;
458             } else {
459                 signr = (r<=0)?(-1):(1/*(r==0)?0:1*/);
460                 stdform[0] = signr*(stdform[6]*stdform[6]+stdform[7]*stdform[7]-r*r)*0.5;
461                 stdform[1] = -signr*stdform[6];
462                 stdform[2] = -signr*stdform[7];
463                 stdform[3] = signr/2;
464                 stdform[4] = signr*r;
465             }
466             return stdform;
467         }
468     }; // JXG.Math
469 })(JXG, Math);
470