      1 /*
      2     Copyright 2008-2017
      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          * Determine the relative difference between two numbers.
     97          * @param  {Number} a First number
     98          * @param  {Number} b Second number
     99          * @returns {Number}  Relative difference between a and b: |a-b| / max(|a|, |b|)
    100          */
    101         relDif: function(a, b) {
    102             var c = Math.abs(a),
    103                 d = Math.abs(b);
    104
    105             d = Math.max(c, d);
    106
    107             return (d === 0.0) ? 0.0 : Math.abs(a - b) / d;
    108         },
    109
    110         /**
    111          * The JavaScript implementation of the % operator returns the symmetric modulo.
    112          * They are both identical if a >= 0 and m >= 0 but the results differ if a or m < 0.
    113          * @param {Number} a
    114          * @param {Number} m
    115          * @returns {Number} Mathematical modulo <tt>a mod m</tt>
    116          */
    117         mod: function (a, m) {
    118             return a - Math.floor(a / m) * m;
    119         },
    120
    121         /**
    122          * Initializes a vector as an array with the coefficients set to the given value resp. zero.
    123          * @param {Number} n Length of the vector
    124          * @param {Number} [init=0] Initial value for each coefficient
    125          * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a
    126          * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows.
    127          */
    128         vector: function (n, init) {
    129             var r, i;
    130
    131             init = init || 0;
    132             r = [];
    133
    134             for (i = 0; i < n; i++) {
    135                 r[i] = init;
    136             }
    137
    138             return r;
    139         },
    140
    141         /**
    142          * Initializes a matrix as an array of rows with the given value.
    143          * @param {Number} n Number of rows
    144          * @param {Number} [m=n] Number of columns
    145          * @param {Number} [init=0] Initial value for each coefficient
    146          * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a
    147          * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows.
    148          */
    149         matrix: function (n, m, init) {
    150             var r, i, j;
    151
    152             init = init || 0;
    153             m = m || n;
    154             r = [];
    155
    156             for (i = 0; i < n; i++) {
    157                 r[i] = [];
    158
    159                 for (j = 0; j < m; j++) {
    160                     r[i][j] = init;
    161                 }
    162             }
    163
    164             return r;
    165         },
    166
    167         /**
    168          * Generates an identity matrix. If n is a number and m is undefined or not a number, a square matrix is generated,
    169          * if n and m are both numbers, an nxm matrix is generated.
    170          * @param {Number} n Number of rows
    171          * @param {Number} [m=n] Number of columns
    172          * @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
    173          * 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.
    174          */
    175         identity: function (n, m) {
    176             var r, i;
    177
    178             if ((m === undef) && (typeof m !== 'number')) {
    179                 m = n;
    180             }
    181
    182             r = this.matrix(n, m);
    183
    184             for (i = 0; i < Math.min(n, m); i++) {
    185                 r[i][i] = 1;
    186             }
    187
    188             return r;
    189         },
    190
    191         /**
    192          * Generates a 4x4 matrix for 3D to 2D projections.
    193          * @param {Number} l Left
    194          * @param {Number} r Right
    195          * @param {Number} t Top
    196          * @param {Number} b Bottom
    197          * @param {Number} n Near
    198          * @param {Number} f Far
    199          * @returns {Array} 4x4 Matrix
    200          */
    201         frustum: function (l, r, b, t, n, f) {
    202             var ret = this.matrix(4, 4);
    203
    204             ret[0][0] = (n * 2) / (r - l);
    205             ret[0][1] = 0;
    206             ret[0][2] = (r + l) / (r - l);
    207             ret[0][3] = 0;
    208
    209             ret[1][0] = 0;
    210             ret[1][1] = (n * 2) / (t - b);
    211             ret[1][2] = (t + b) / (t - b);
    212             ret[1][3] = 0;
    213
    214             ret[2][0] = 0;
    215             ret[2][1] = 0;
    216             ret[2][2] = -(f + n) / (f - n);
    217             ret[2][3] = -(f * n * 2) / (f - n);
    218
    219             ret[3][0] = 0;
    220             ret[3][1] = 0;
    221             ret[3][2] = -1;
    222             ret[3][3] = 0;
    223
    224             return ret;
    225         },
    226
    227         /**
    228          * Generates a 4x4 matrix for 3D to 2D projections.
    229          * @param {Number} fov Field of view in vertical direction, given in rad.
    230          * @param {Number} ratio Aspect ratio of the projection plane.
    231          * @param {Number} n Near
    232          * @param {Number} f Far
    233          * @returns {Array} 4x4 Projection Matrix
    234          */
    235         projection: function (fov, ratio, n, f) {
    236             var t = n * Math.tan(fov / 2),
    237                 r = t * ratio;
    238
    239             return this.frustum(-r, r, -t, t, n, f);
    240         },
    241
    242         /**
    243          * 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
    244          * function does not check if the dimensions match.
    245          * @param {Array} mat Two dimensional array of numbers. The inner arrays describe the columns, the outer ones the matrix' rows.
    246          * @param {Array} vec Array of numbers
    247          * @returns {Array} Array of numbers containing the result
    248          * @example
    249          * var A = [[2, 1],
    250          *          [1, 3]],
    251          *     b = [4, 5],
    252          *     c;
    253          * c = JXG.Math.matVecMult(A, b)
    254          * // c === [13, 19];
    255          */
    256         matVecMult: function (mat, vec) {
    257             var i, s, k,
    258                 m = mat.length,
    259                 n = vec.length,
    260                 res = [];
    261
    262             if (n === 3) {
    263                 for (i = 0; i < m; i++) {
    264                     res[i] = mat[i][0] * vec[0] + mat[i][1] * vec[1] + mat[i][2] * vec[2];
    265                 }
    266             } else {
    267                 for (i = 0; i < m; i++) {
    268                     s = 0;
    269                     for (k = 0; k < n; k++) {
    270                         s += mat[i][k] * vec[k];
    271                     }
    272                     res[i] = s;
    273                 }
    274             }
    275             return res;
    276         },
    277
    278         /**
    279          * Computes the product of the two matrices mat1*mat2.
    280          * @param {Array} mat1 Two dimensional array of numbers
    281          * @param {Array} mat2 Two dimensional array of numbers
    282          * @returns {Array} Two dimensional Array of numbers containing result
    283          */
    284         matMatMult: function (mat1, mat2) {
    285             var i, j, s, k,
    286                 m = mat1.length,
    287                 n = m > 0 ? mat2[0].length : 0,
    288                 m2 = mat2.length,
    289                 res = this.matrix(m, n);
    290
    291             for (i = 0; i < m; i++) {
    292                 for (j = 0; j < n; j++) {
    293                     s = 0;
    294                     for (k = 0; k < m2; k++) {
    295                         s += mat1[i][k] * mat2[k][j];
    296                     }
    297                     res[i][j] = s;
    298                 }
    299             }
    300             return res;
    301         },
    302
    303         /**
    304          * Transposes a matrix given as a two dimensional array.
    305          * @param {Array} M The matrix to be transposed
    306          * @returns {Array} The transpose of M
    307          */
    308         transpose: function (M) {
    309             var MT, i, j,
    310                 m, n;
    311
    312             // number of rows of M
    313             m = M.length;
    314             // number of columns of M
    315             n = M.length > 0 ? M[0].length : 0;
    316             MT = this.matrix(n, m);
    317
    318             for (i = 0; i < n; i++) {
    319                 for (j = 0; j < m; j++) {
    320                     MT[i][j] = M[j][i];
    321                 }
    322             }
    323
    324             return MT;
    325         },
    326
    327         /**
    328          * Compute the inverse of an nxn matrix with Gauss elimination.
    329          * @param {Array} Ain
    330          * @returns {Array} Inverse matrix of Ain
    331          */
    332         inverse: function (Ain) {
    333             var i, j, k, s, ma, r, swp,
    334                 n = Ain.length,
    335                 A = [],
    336                 p = [],
    337                 hv = [];
    338
    339             for (i = 0; i < n; i++) {
    340                 A[i] = [];
    341                 for (j = 0; j < n; j++) {
    342                     A[i][j] = Ain[i][j];
    343                 }
    344                 p[i] = i;
    345             }
    346
    347             for (j = 0; j < n; j++) {
    348                 // pivot search:
    349                 ma = Math.abs(A[j][j]);
    350                 r = j;
    351
    352                 for (i = j + 1; i < n; i++) {
    353                     if (Math.abs(A[i][j]) > ma) {
    354                         ma = Math.abs(A[i][j]);
    355                         r = i;
    356                     }
    357                 }
    358
    359                 // Singular matrix
    360                 if (ma <= this.eps) {
    361                     return [];
    362                 }
    363
    364                 // swap rows:
    365                 if (r > j) {
    366                     for (k = 0; k < n; k++) {
    367                         swp = A[j][k];
    368                         A[j][k] = A[r][k];
    369                         A[r][k] = swp;
    370                     }
    371
    372                     swp = p[j];
    373                     p[j] = p[r];
    374                     p[r] = swp;
    375                 }
    376
    377                 // transformation:
    378                 s = 1.0 / A[j][j];
    379                 for (i = 0; i < n; i++) {
    380                     A[i][j] *= s;
    381                 }
    382                 A[j][j] = s;
    383
    384                 for (k = 0; k < n; k++) {
    385                     if (k !== j) {
    386                         for (i = 0; i < n; i++) {
    387                             if (i !== j) {
    388                                 A[i][k] -= A[i][j] * A[j][k];
    389                             }
    390                         }
    391                         A[j][k] = -s * A[j][k];
    392                     }
    393                 }
    394             }
    395
    396             // swap columns:
    397             for (i = 0; i < n; i++) {
    398                 for (k = 0; k < n; k++) {
    399                     hv[p[k]] = A[i][k];
    400                 }
    401                 for (k = 0; k < n; k++) {
    402                     A[i][k] = hv[k];
    403                 }
    404             }
    405
    406             return A;
    407         },
    408
    409         /**
    410          * Inner product of two vectors a and b. n is the length of the vectors.
    411          * @param {Array} a Vector
    412          * @param {Array} b Vector
    413          * @param {Number} [n] Length of the Vectors. If not given the length of the first vector is taken.
    414          * @returns {Number} The inner product of a and b.
    415          */
    416         innerProduct: function (a, b, n) {
    417             var i,
    418                 s = 0;
    419
    420             if (n === undef || !Type.isNumber(n)) {
    421                 n = a.length;
    422             }
    423
    424             for (i = 0; i < n; i++) {
    425                 s += a[i] * b[i];
    426             }
    427
    428             return s;
    429         },
    430
    431         /**
    432          * Calculates the cross product of two vectors both of length three.
    433          * In case of homogeneous coordinates this is either
    434          * <ul>
    435          * <li>the intersection of two lines</li>
    436          * <li>the line through two points</li>
    437          * </ul>
    438          * @param {Array} c1 Homogeneous coordinates of line or point 1
    439          * @param {Array} c2 Homogeneous coordinates of line or point 2
    440          * @returns {Array} vector of length 3: homogeneous coordinates of the resulting point / line.
    441          */
    442         crossProduct: function (c1, c2) {
    443             return [c1[1] * c2[2] - c1[2] * c2[1],
    444                 c1[2] * c2[0] - c1[0] * c2[2],
    445                 c1[0] * c2[1] - c1[1] * c2[0]];
    446         },
    447
    448         /**
    449          * Compute the factorial of a positive integer. If a non-integer value
    450          * is given, the fraction will be ignored.
    451          * @function
    452          * @param {Number} n
    453          * @returns {Number} n! = n*(n-1)*...*2*1
    454          */
    455         factorial: memoizer(function (n) {
    456             if (n < 0) {
    457                 return NaN;
    458             }
    459
    460             n = Math.floor(n);
    461
    462             if (n === 0 || n === 1) {
    463                 return 1;
    464             }
    465
    466             return n * this.factorial(n - 1);
    467         }),
    468
    469         /**
    470          * Computes the binomial coefficient n over k.
    471          * @function
    472          * @param {Number} n Fraction will be ignored
    473          * @param {Number} k Fraction will be ignored
    474          * @returns {Number} The binomial coefficient n over k
    475          */
    476         binomial: memoizer(function (n, k) {
    477             var b, i;
    478
    479             if (k > n || k < 0) {
    480                 return NaN;
    481             }
    482
    483             k = Math.round(k);
    484             n = Math.round(n);
    485
    486             if (k === 0 || k === n) {
    487                 return 1;
    488             }
    489
    490             b = 1;
    491
    492             for (i = 0; i < k; i++) {
    493                 b *= (n - i);
    494                 b /= (i + 1);
    495             }
    496
    497             return b;
    498         }),
    499
    500         /**
    501          * Calculates the cosine hyperbolicus of x.
    502          * @param {Number} x The number the cosine hyperbolicus will be calculated of.
    503          * @returns {Number} Cosine hyperbolicus of the given value.
    504          */
    505         cosh: function (x) {
    506             return (Math.exp(x) + Math.exp(-x)) * 0.5;
    507         },
    508
    509         /**
    510          * Sine hyperbolicus of x.
    511          * @param {Number} x The number the sine hyperbolicus will be calculated of.
    512          * @returns {Number} Sine hyperbolicus of the given value.
    513          */
    514         sinh: function (x) {
    515             return (Math.exp(x) - Math.exp(-x)) * 0.5;
    516         },
    517
    518         /**
    519          * Compute base to the power of exponent.
    520          * @param {Number} base
    521          * @param {Number} exponent
    522          * @returns {Number} base to the power of exponent.
    523          */
    524         pow: function (base, exponent) {
    525             if (base === 0) {
    526                 if (exponent === 0) {
    527                     return 1;
    528                 }
    529
    530                 return 0;
    531             }
    532
    533             if (Math.floor(exponent) === exponent) {
    534                 // a is an integer
    535                 return Math.pow(base, exponent);
    536             }
    537
    538             // a is not an integer
    539             if (base > 0) {
    540                 return Math.exp(exponent * Math.log(Math.abs(base)));
    541             }
    542
    543             return NaN;
    544         },
    545
    546         /**
    547          * Logarithm to base 10.
    548          * @param {Number} x
    549          * @returns {Number} log10(x) Logarithm of x to base 10.
    550          */
    551         log10: function (x) {
    552             return Math.log(x) / Math.log(10.0);
    553         },
    554
    555         /**
    556          * Logarithm to base 2.
    557          * @param {Number} x
    558          * @returns {Number} log2(x) Logarithm of x to base 2.
    559          */
    560         log2: function (x) {
    561             return Math.log(x) / Math.log(2.0);
    562         },
    563
    564         /**
    565          * Logarithm to arbitrary base b. If b is not given, natural log is taken, i.e. b = e.
    566          * @param {Number} x
    567          * @param {Number} b base
    568          * @returns {Number} log(x, b) Logarithm of x to base b, that is log(x)/log(b).
    569          */
    570         log: function (x, b) {
    571             if (b !== undefined && Type.isNumber(b)) {
    572                 return Math.log(x) / Math.log(b);
    573             }
    574
    575             return Math.log(x);
    576         },
    577
    578         /**
    579          * The sign() function returns the sign of a number, indicating whether the number is positive, negative or zero.
    580          * @param  {Number} x A Number
    581          * @returns {[type]}  This function has 5 kinds of return values,
    582          *    1, -1, 0, -0, NaN, which represent "positive number", "negative number", "positive zero", "negative zero"
    583          *    and NaN respectively.
    584          */
    585         sign: Math.sign || function(x) {
    586             x = +x; // convert to a number
    587             if (x === 0 || isNaN(x)) {
    588                 return x;
    589             }
    590             return x > 0 ? 1 : -1;
    591         },
    592
    593         /**
    594          * A square & multiply algorithm to compute base to the power of exponent.
    595          * Implementated by Wolfgang Riedl.
    596          * @param {Number} base
    597          * @param {Number} exponent
    598          * @returns {Number} Base to the power of exponent
    599          */
    600         squampow: function (base, exponent) {
    601             var result;
    602
    603             if (Math.floor(exponent) === exponent) {
    604                 // exponent is integer (could be zero)
    605                 result = 1;
    606
    607                 if (exponent < 0) {
    608                     // invert: base
    609                     base = 1.0 / base;
    610                     exponent *= -1;
    611                 }
    612
    613                 while (exponent !== 0) {
    614                     if (exponent & 1) {
    615                         result *= base;
    616                     }
    617
    618                     exponent >>= 1;
    619                     base *= base;
    620                 }
    621                 return result;
    622             }
    623
    624             return this.pow(base, exponent);
    625         },
    626
    627         /**
    628          * Greatest common divisor (gcd) of two numbers.
    629          * @see http://rosettacode.org/wiki/Greatest_common_divisor#JavaScript
    630          *
    631          * @param  {Number} a First number
    632          * @param  {Number} b Second number
    633          * @returns {Number}   gcd(a, b) if a and b are numbers, NaN else.
    634          */
    635         gcd: function (a,b) {
    636             a = Math.abs(a);
    637             b = Math.abs(b);
    638
    639             if (!(Type.isNumber(a) && Type.isNumber(b))) {
    640                 return NaN;
    641             }
    642             if (b > a) {
    643                 var temp = a;
    644                 a = b;
    645                 b = temp;
    646             }
    647
    648             while (true) {
    649                 a %= b;
    650                 if (a === 0) { return b; }
    651                 b %= a;
    652                 if (b === 0) { return a; }
    653             }
    654         },
    655
    656         /**
    657          * Normalize the standard form [c, b0, b1, a, k, r, q0, q1].
    658          * @private
    659          * @param {Array} stdform The standard form to be normalized.
    660          * @returns {Array} The normalized standard form.
    661          */
    662         normalize: function (stdform) {
    663             var n, signr,
    664                 a2 = 2 * stdform[3],
    665                 r = stdform[4] / a2;
    666
    667             stdform[5] = r;
    668             stdform[6] = -stdform[1] / a2;
    669             stdform[7] = -stdform[2] / a2;
    670
    671             if (!isFinite(r)) {
    672                 n = Math.sqrt(stdform[1] * stdform[1] + stdform[2] * stdform[2]);
    673
    674                 stdform[0] /= n;
    675                 stdform[1] /= n;
    676                 stdform[2] /= n;
    677                 stdform[3] = 0;
    678                 stdform[4] = 1;
    679             } else if (Math.abs(r) >= 1) {
    680                 stdform[0] = (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) / (2 * r);
    681                 stdform[1] = -stdform[6] / r;
    682                 stdform[2] = -stdform[7] / r;
    683                 stdform[3] = 1 / (2 * r);
    684                 stdform[4] = 1;
    685             } else {
    686                 signr = (r <= 0 ? -1 : 1);
    687                 stdform[0] = signr * (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) * 0.5;
    688                 stdform[1] = -signr * stdform[6];
    689                 stdform[2] = -signr * stdform[7];
    690                 stdform[3] = signr / 2;
    691                 stdform[4] = signr * r;
    692             }
    693
    694             return stdform;
    695         },
    696
    697         /**
    698          * Converts a two dimensional array to a one dimensional Float32Array that can be processed by WebGL.
    699          * @param {Array} m A matrix in a two dimensional array.
    700          * @returns {Float32Array} A one dimensional array containing the matrix in column wise notation. Provides a fall
    701          * back to the default JavaScript Array if Float32Array is not available.
    702          */
    703         toGL: function (m) {
    704             var v, i, j;
    705
    706             if (typeof Float32Array === 'function') {
    707                 v = new Float32Array(16);
    708             } else {
    709                 v = new Array(16);
    710             }
    711
    712             if (m.length !== 4 && m[0].length !== 4) {
    713                 return v;
    714             }
    715
    716             for (i = 0; i < 4; i++) {
    717                 for (j = 0; j < 4; j++) {
    718                     v[i + 4 * j] = m[i][j];
    719                 }
    720             }
    721
    722             return v;
    723         }
    724     };
    725
    726     return JXG.Math;
    727 });
    728