      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