      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