      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*/
     34 /*jslint nomen: true, plusplus: true*/
     35
     36 /* depends:
     37  utils/type
     38  math/math
     39  */
     40
     41 /**
     42  * @fileoverview In this file the namespace Math.Numerics is defined, which holds numerical
     43  * algorithms for solving linear equations etc.
     44  */
     45
     46 define(['utils/type', 'math/math'], function (Type, Mat) {
     47
     48     "use strict";
     49
     50     // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order).
     51     var predefinedButcher = {
     52         rk4: {
     53             s: 4,
     54             A: [
     55                 [ 0,  0,  0, 0],
     56                 [0.5, 0,  0, 0],
     57                 [ 0, 0.5, 0, 0],
     58                 [ 0,  0,  1, 0]
     59             ],
     60             b: [1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0],
     61             c: [0, 0.5, 0.5, 1]
     62         },
     63         heun: {
     64             s: 2,
     65             A: [
     66                 [0, 0],
     67                 [1, 0]
     68             ],
     69             b: [0.5, 0.5],
     70             c: [0, 1]
     71         },
     72         euler: {
     73             s: 1,
     74             A: [
     75                 [0]
     76             ],
     77             b: [1],
     78             c: [0]
     79         }
     80     };
     81
     82     /**
     83      * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables.
     84      * @namespace
     85      * @name JXG.Math.Numerics
     86      */
     87     Mat.Numerics = {
     88
     89     //JXG.extend(Mat.Numerics, /** @lends JXG.Math.Numerics */ {
     90         /**
     91          * Solves a system of linear equations given by A and b using the Gauss-Jordan-elimination.
     92          * The algorithm runs in-place. I.e. the entries of A and b are changed.
     93          * @param {Array} A Square matrix represented by an array of rows, containing the coefficients of the lineare equation system.
     94          * @param {Array} b A vector containing the linear equation system's right hand side.
     95          * @throws {Error} If a non-square-matrix is given or if b has not the right length or A's rank is not full.
     96          * @returns {Array} A vector that solves the linear equation system.
     97          * @memberof JXG.Math.Numerics
     98          */
     99         Gauss: function (A, b) {
    100             var i, j, k,
    101                 // copy the matrix to prevent changes in the original
    102                 Acopy,
    103                 // solution vector, to prevent changing b
    104                 x,
    105                 eps = Mat.eps,
    106                 // number of columns of A
    107                 n = A.length > 0 ? A[0].length : 0;
    108
    109             if ((n !== b.length) || (n !== A.length)) {
    110                 throw new Error("JXG.Math.Numerics.Gauss: Dimensions don't match. A must be a square matrix and b must be of the same length as A.");
    111             }
    112
    113             // initialize solution vector
    114             Acopy = [];
    115             x = b.slice(0, n);
    116
    117             for (i = 0; i < n; i++) {
    118                 Acopy[i] = A[i].slice(0, n);
    119             }
    120
    121             // Gauss-Jordan-elimination
    122             for (j = 0; j < n; j++) {
    123                 for (i = n - 1; i > j; i--) {
    124                     // Is the element which is to eliminate greater than zero?
    125                     if (Math.abs(Acopy[i][j]) > eps) {
    126                         // Equals pivot element zero?
    127                         if (Math.abs(Acopy[j][j]) < eps) {
    128                             // At least numerically, so we have to exchange the rows
    129                             Type.swap(Acopy, i, j);
    130                             Type.swap(x, i, j);
    131                         } else {
    132                             // Saves the L matrix of the LR-decomposition. unnecessary.
    133                             Acopy[i][j] /= Acopy[j][j];
    134                             // Transform right-hand-side b
    135                             x[i] -= Acopy[i][j] * x[j];
    136
    137                             // subtract the multiple of A[i][j] / A[j][j] of the j-th row from the i-th.
    138                             for (k = j + 1; k < n; k++) {
    139                                 Acopy[i][k] -= Acopy[i][j] * Acopy[j][k];
    140                             }
    141                         }
    142                     }
    143                 }
    144
    145                 // The absolute values of all coefficients below the j-th row in the j-th column are smaller than JXG.Math.eps.
    146                 if (Math.abs(Acopy[j][j]) < eps) {
    147                     throw new Error("JXG.Math.Numerics.Gauss(): The given matrix seems to be singular.");
    148                 }
    149             }
    150
    151             this.backwardSolve(Acopy, x, true);
    152
    153             return x;
    154         },
    155
    156         /**
    157          * Solves a system of linear equations given by the right triangular matrix R and vector b.
    158          * @param {Array} R Right triangular matrix represented by an array of rows. All entries a_(i,j) with i < j are ignored.
    159          * @param {Array} b Right hand side of the linear equation system.
    160          * @param {Boolean} [canModify=false] If true, the right hand side vector is allowed to be changed by this method.
    161          * @returns {Array} An array representing a vector that solves the system of linear equations.
    162          * @memberof JXG.Math.Numerics
    163          */
    164         backwardSolve: function (R, b, canModify) {
    165             var x, m, n, i, j;
    166
    167             if (canModify) {
    168                 x = b;
    169             } else {
    170                 x = b.slice(0, b.length);
    171             }
    172
    173             // m: number of rows of R
    174             // n: number of columns of R
    175             m = R.length;
    176             n = R.length > 0 ? R[0].length : 0;
    177
    178             for (i = m - 1; i >= 0; i--) {
    179                 for (j = n - 1; j > i; j--) {
    180                     x[i] -= R[i][j] * x[j];
    181                 }
    182                 x[i] /= R[i][i];
    183             }
    184
    185             return x;
    186         },
    187
    188         /**
    189          * @private
    190          * Gauss-Bareiss algorithm to compute the
    191          * determinant of matrix without fractions.
    192          * See Henri Cohen, "A Course in Computational
    193          * Algebraic Number Theory (Graduate texts
    194          * in mathematics; 138)", Springer-Verlag,
    195          * ISBN 3-540-55640-0 / 0-387-55640-0
    196          * Third, Corrected Printing 1996
    197          * "Algorithm 2.2.6", pg. 52-53
    198          * @memberof JXG.Math.Numerics
    199          */
    200         gaussBareiss: function (mat) {
    201             var k, c, s, i, j, p, n, M, t,
    202                 eps = Mat.eps;
    203
    204             n = mat.length;
    205
    206             if (n <= 0) {
    207                 return 0;
    208             }
    209
    210             if (mat[0].length < n) {
    211                 n = mat[0].length;
    212             }
    213
    214             // Copy the input matrix to M
    215             M = [];
    216
    217             for (i = 0; i < n; i++) {
    218                 M[i] = mat[i].slice(0, n);
    219             }
    220
    221             c = 1;
    222             s = 1;
    223
    224             for (k = 0; k < n - 1; k++) {
    225                 p = M[k][k];
    226
    227                 // Pivot step
    228                 if (Math.abs(p) < eps) {
    229                     for (i = k + 1; i < n; i++) {
    230                         if (Math.abs(M[i][k]) >= eps) {
    231                             break;
    232                         }
    233                     }
    234
    235                     // No nonzero entry found in column k -> det(M) = 0
    236                     if (i === n) {
    237                         return 0.0;
    238                     }
    239
    240                     // swap row i and k partially
    241                     for (j = k; j < n; j++) {
    242                         t = M[i][j];
    243                         M[i][j] = M[k][j];
    244                         M[k][j] = t;
    245                     }
    246                     s = -s;
    247                     p = M[k][k];
    248                 }
    249
    250                 // Main step
    251                 for (i = k + 1; i < n; i++) {
    252                     for (j = k + 1; j < n; j++) {
    253                         t = p * M[i][j] - M[i][k] * M[k][j];
    254                         M[i][j] = t / c;
    255                     }
    256                 }
    257
    258                 c = p;
    259             }
    260
    261             return s * M[n - 1][n - 1];
    262         },
    263
    264         /**
    265          * Computes the determinant of a square nxn matrix with the
    266          * Gauss-Bareiss algorithm.
    267          * @param {Array} mat Matrix.
    268          * @returns {Number} The determinant pf the matrix mat.
    269          *                   The empty matrix returns 0.
    270          * @memberof JXG.Math.Numerics
    271          */
    272         det: function (mat) {
    273             var n = mat.length;
    274
    275             if (n === 2 && mat[0].length === 2) {
    276                 return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
    277             }
    278
    279             return this.gaussBareiss(mat);
    280         },
    281
    282         /**
    283          * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method
    284          * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990
    285          * @param {Array} Ain A symmetric 3x3 matrix.
    286          * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors.
    287          * @memberof JXG.Math.Numerics
    288          */
    289         Jacobi: function (Ain) {
    290             var i, j, k, aa, si, co, tt, ssum, amax,
    291                 eps = Mat.eps,
    292                 sum = 0.0,
    293                 n = Ain.length,
    294                 V = [
    295                     [0, 0, 0],
    296                     [0, 0, 0],
    297                     [0, 0, 0]
    298                 ],
    299                 A = [
    300                     [0, 0, 0],
    301                     [0, 0, 0],
    302                     [0, 0, 0]
    303                 ],
    304                 nloops = 0;
    305
    306             // Initialization. Set initial Eigenvectors.
    307             for (i = 0; i < n; i++) {
    308                 for (j = 0; j < n; j++) {
    309                     V[i][j] = 0.0;
    310                     A[i][j] = Ain[i][j];
    311                     sum += Math.abs(A[i][j]);
    312                 }
    313                 V[i][i] = 1.0;
    314             }
    315
    316             // Trivial problems
    317             if (n === 1) {
    318                 return [A, V];
    319             }
    320
    321             if (sum <= 0.0) {
    322                 return [A, V];
    323             }
    324
    325             sum /= (n * n);
    326
    327             // Reduce matrix to diagonal
    328             do {
    329                 ssum = 0.0;
    330                 amax = 0.0;
    331                 for (j = 1; j < n; j++) {
    332                     for (i = 0; i < j; i++) {
    333                         // Check if A[i][j] is to be reduced
    334                         aa = Math.abs(A[i][j]);
    335
    336                         if (aa > amax) {
    337                             amax = aa;
    338                         }
    339
    340                         ssum += aa;
    341
    342                         if (aa >= eps) {
    343                             // calculate rotation angle
    344                             aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5;
    345                             si = Math.sin(aa);
    346                             co = Math.cos(aa);
    347
    348                             // Modify 'i' and 'j' columns
    349                             for (k = 0; k < n; k++) {
    350                                 tt = A[k][i];
    351                                 A[k][i] = co * tt + si * A[k][j];
    352                                 A[k][j] = -si * tt + co * A[k][j];
    353                                 tt = V[k][i];
    354                                 V[k][i] = co * tt + si * V[k][j];
    355                                 V[k][j] = -si * tt + co * V[k][j];
    356                             }
    357
    358                             // Modify diagonal terms
    359                             A[i][i] = co * A[i][i] + si * A[j][i];
    360                             A[j][j] = -si * A[i][j] + co * A[j][j];
    361                             A[i][j] = 0.0;
    362
    363                             // Make 'A' matrix symmetrical
    364                             for (k = 0; k < n; k++) {
    365                                 A[i][k] = A[k][i];
    366                                 A[j][k] = A[k][j];
    367                             }
    368                             // A[i][j] made zero by rotation
    369                         }
    370                     }
    371                 }
    372                 nloops += 1;
    373             } while (Math.abs(ssum) / sum > eps && nloops < 2000);
    374
    375             return [A, V];
    376         },
    377
    378         /**
    379          * Calculates the integral of function f over interval using Newton-Cotes-algorithm.
    380          * @param {Array} interval The integration interval, e.g. [0, 3].
    381          * @param {function} f A function which takes one argument of type number and returns a number.
    382          * @param {Object} [config] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type
    383          * with value being either 'trapez', 'simpson', or 'milne'.
    384          * @param {Number} [config.number_of_nodes=28]
    385          * @param {String} [config.integration_type='milne'] Possible values are 'milne', 'simpson', 'trapez'
    386          * @returns {Number} Integral value of f over interval
    387          * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use
    388          * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4.
    389          * @example
    390          * function f(x) {
    391          *   return x*x;
    392          * }
    393          *
    394          * // calculates integral of <tt>f</tt> from 0 to 2.
    395          * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f);
    396          *
    397          * // the same with an anonymous function
    398          * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; });
    399          *
    400          * // use trapez rule with 16 nodes
    401          * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f,
    402          *                                   {number_of_nodes: 16, intergration_type: 'trapez'});
    403          * @memberof JXG.Math.Numerics
    404          */
    405         NewtonCotes: function (interval, f, config) {
    406             var evaluation_point, i, number_of_intervals,
    407                 integral_value = 0.0,
    408                 number_of_nodes = config && typeof config.number_of_nodes === 'number' ? config.number_of_nodes : 28,
    409                 available_types = {trapez: true, simpson: true, milne: true},
    410                 integration_type = config && config.integration_type && available_types.hasOwnProperty(config.integration_type) && available_types[config.integration_type] ? config.integration_type : 'milne',
    411                 step_size = (interval[1] - interval[0]) / number_of_nodes;
    412
    413             switch (integration_type) {
    414             case 'trapez':
    415                 integral_value = (f(interval[0]) + f(interval[1])) * 0.5;
    416                 evaluation_point = interval[0];
    417
    418                 for (i = 0; i < number_of_nodes - 1; i++) {
    419                     evaluation_point += step_size;
    420                     integral_value += f(evaluation_point);
    421                 }
    422
    423                 integral_value *= step_size;
    424                 break;
    425             case 'simpson':
    426                 if (number_of_nodes % 2 > 0) {
    427                     throw new Error("JSXGraph:  INT_SIMPSON requires config.number_of_nodes dividable by 2.");
    428                 }
    429
    430                 number_of_intervals = number_of_nodes / 2.0;
    431                 integral_value = f(interval[0]) + f(interval[1]);
    432                 evaluation_point = interval[0];
    433
    434                 for (i = 0; i < number_of_intervals - 1; i++) {
    435                     evaluation_point += 2.0 * step_size;
    436                     integral_value += 2.0 * f(evaluation_point);
    437                 }
    438
    439                 evaluation_point = interval[0] - step_size;
    440
    441                 for (i = 0; i < number_of_intervals; i++) {
    442                     evaluation_point += 2.0 * step_size;
    443                     integral_value += 4.0 * f(evaluation_point);
    444                 }
    445
    446                 integral_value *= step_size / 3.0;
    447                 break;
    448             default:
    449                 if (number_of_nodes % 4 > 0) {
    450                     throw new Error("JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4");
    451                 }
    452
    453                 number_of_intervals = number_of_nodes * 0.25;
    454                 integral_value = 7.0 * (f(interval[0]) + f(interval[1]));
    455                 evaluation_point = interval[0];
    456
    457                 for (i = 0; i < number_of_intervals - 1; i++) {
    458                     evaluation_point += 4.0 * step_size;
    459                     integral_value += 14.0 * f(evaluation_point);
    460                 }
    461
    462                 evaluation_point = interval[0] - 3.0 * step_size;
    463
    464                 for (i = 0; i < number_of_intervals; i++) {
    465                     evaluation_point += 4.0 * step_size;
    466                     integral_value += 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size));
    467                 }
    468
    469                 evaluation_point = interval[0] - 2.0 * step_size;
    470
    471                 for (i = 0; i < number_of_intervals; i++) {
    472                     evaluation_point += 4.0 * step_size;
    473                     integral_value += 12.0 * f(evaluation_point);
    474                 }
    475
    476                 integral_value *= 2.0 * step_size / 45.0;
    477             }
    478             return integral_value;
    479         },
    480
    481         /**
    482          * Integral of function f over interval.
    483          * @param {Array} interval The integration interval, e.g. [0, 3].
    484          * @param {function} f A function which takes one argument of type number and returns a number.
    485          * @returns {Number} The value of the integral of f over interval
    486          * @see JXG.Math.Numerics.NewtonCotes
    487          * @memberof JXG.Math.Numerics
    488          */
    489         I: function (interval, f) {
    490             return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'});
    491         },
    492
    493         /**
    494          * Newton's method to find roots of a funtion in one variable.
    495          * @param {function} f We search for a solution of f(x)=0.
    496          * @param {Number} x initial guess for the root, i.e. start value.
    497          * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
    498          * the function is a method of an object and contains a reference to its parent object via "this".
    499          * @returns {Number} A root of the function f.
    500          * @memberof JXG.Math.Numerics
    501          */
    502         Newton: function (f, x, context) {
    503             var df,
    504                 i = 0,
    505                 h = Mat.eps,
    506                 newf = f.apply(context, [x]),
    507                 nfev = 1;
    508
    509             // For compatibility
    510             if (Type.isArray(x)) {
    511                 x = x[0];
    512             }
    513
    514             while (i < 50 && Math.abs(newf) > h) {
    515                 df = this.D(f, context)(x);
    516                 nfev += 2;
    517
    518                 if (Math.abs(df) > h) {
    519                     x -= newf / df;
    520                 } else {
    521                     x += (Math.random() * 0.2 - 1.0);
    522                 }
    523
    524                 newf = f.apply(context, [x]);
    525                 nfev += 1;
    526                 i += 1;
    527             }
    528
    529             return x;
    530         },
    531
    532         /**
    533          * Abstract method to find roots of univariate functions.
    534          * @param {function} f We search for a solution of f(x)=0.
    535          * @param {Number} x initial guess for the root, i.e. starting value.
    536          * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
    537          * the function is a method of an object and contains a reference to its parent object via "this".
    538          * @returns {Number} A root of the function f.
    539          * @memberof JXG.Math.Numerics
    540          */
    541         root: function (f, x, context) {
    542             return this.fzero(f, x, context);
    543         },
    544
    545         /**
    546          * Compute an intersection of the curves c1 and c2
    547          * with a generalized Newton method.
    548          * We want to find values t1, t2 such that
    549          * c1(t1) = c2(t2), i.e.
    550          * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0).
    551          * We set
    552          * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2))
    553          *
    554          * The Jacobian J is defined by
    555          * J = (a, b)
    556          *     (c, d)
    557          * where
    558          * a = c1_x'(t1)
    559          * b = -c2_x'(t2)
    560          * c = c1_y'(t1)
    561          * d = -c2_y'(t2)
    562          *
    563          * The inverse J^(-1) of J is equal to
    564          *  (d, -b)/
    565          *  (-c, a) / (ad-bc)
    566          *
    567          * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f).
    568          * If the function meetCurveCurve possesses the properties
    569          * t1memo and t2memo then these are taken as start values
    570          * for the Newton algorithm.
    571          * After stopping of the Newton algorithm the values of t1 and t2 are stored in
    572          * t1memo and t2memo.
    573          *
    574          * @param {JXG.Curve} c1 Curve, Line or Circle
    575          * @param {JXG.Curve} c2 Curve, Line or Circle
    576          * @param {Number} t1ini start value for t1
    577          * @param {Number} t2ini start value for t2
    578          * @returns {JXG.Coords} intersection point
    579          * @memberof JXG.Math.Numerics
    580          */
    581         generalizedNewton: function (c1, c2, t1ini, t2ini) {
    582             var t1, t2,
    583                 a, b, c, d, disc,
    584                 e, f, F,
    585                 D00, D01,
    586                 D10, D11,
    587                 count = 0;
    588
    589             if (this.generalizedNewton.t1memo) {
    590                 t1 = this.generalizedNewton.t1memo;
    591                 t2 = this.generalizedNewton.t2memo;
    592             } else {
    593                 t1 = t1ini;
    594                 t2 = t2ini;
    595             }
    596
    597             e = c1.X(t1) - c2.X(t2);
    598             f = c1.Y(t1) - c2.Y(t2);
    599             F = e * e + f * f;
    600
    601             D00 = this.D(c1.X, c1);
    602             D01 = this.D(c2.X, c2);
    603             D10 = this.D(c1.Y, c1);
    604             D11 = this.D(c2.Y, c2);
    605
    606             while (F > Mat.eps && count < 10) {
    607                 a = D00(t1);
    608                 b = -D01(t2);
    609                 c = D10(t1);
    610                 d = -D11(t2);
    611                 disc = a * d - b * c;
    612                 t1 -= (d * e - b * f) / disc;
    613                 t2 -= (a * f - c * e) / disc;
    614                 e = c1.X(t1) - c2.X(t2);
    615                 f = c1.Y(t1) - c2.Y(t2);
    616                 F = e * e + f * f;
    617                 count += 1;
    618             }
    619
    620             this.generalizedNewton.t1memo = t1;
    621             this.generalizedNewton.t2memo = t2;
    622
    623             if (Math.abs(t1) < Math.abs(t2)) {
    624                 return [c1.X(t1), c1.Y(t1)];
    625             }
    626
    627             return [c2.X(t2), c2.Y(t2)];
    628         },
    629
    630         /**
    631          * Returns the Lagrange polynomials for curves with equidistant nodes, see
    632          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
    633          * SIAM Review, Vol 46, No 3, (2004) 501-517.
    634          * The graph of the parametric curve [x(t),y(t)] runs through the given points.
    635          * @param {Array} p Array of JXG.Points
    636          * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve
    637          * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero.
    638          * @memberof JXG.Math.Numerics
    639          */
    640         Neville: function (p) {
    641             var w = [],
    642                 /** @ignore */
    643                 makeFct = function (fun) {
    644                     return function (t, suspendedUpdate) {
    645                         var i, d, s,
    646                             bin = Mat.binomial,
    647                             len = p.length,
    648                             len1 = len - 1,
    649                             num = 0.0,
    650                             denom = 0.0;
    651
    652                         if (!suspendedUpdate) {
    653                             s = 1;
    654                             for (i = 0; i < len; i++) {
    655                                 w[i] = bin(len1, i) * s;
    656                                 s *= (-1);
    657                             }
    658                         }
    659
    660                         d = t;
    661
    662                         for (i = 0; i < len; i++) {
    663                             if (d === 0) {
    664                                 return p[i][fun]();
    665                             }
    666                             s = w[i] / d;
    667                             d -= 1;
    668                             num += p[i][fun]() * s;
    669                             denom += s;
    670                         }
    671                         return num / denom;
    672                     };
    673                 },
    674
    675                 xfct = makeFct('X'),
    676                 yfct = makeFct('Y');
    677
    678             return [xfct, yfct, 0, function () {
    679                 return p.length - 1;
    680             }];
    681         },
    682
    683         /**
    684          * Calculates second derivatives at the given knots.
    685          * @param {Array} x x values of knots
    686          * @param {Array} y y values of knots
    687          * @returns {Array} Second derivatives of the interpolated function at the knots.
    688          * @see #splineEval
    689          * @memberof JXG.Math.Numerics
    690          */
    691         splineDef: function (x, y) {
    692             var pair, i, l,
    693                 n = Math.min(x.length, y.length),
    694                 diag = [],
    695                 z = [],
    696                 data = [],
    697                 dx = [],
    698                 delta = [],
    699                 F = [];
    700
    701             if (n === 2) {
    702                 return [0, 0];
    703             }
    704
    705             for (i = 0; i < n; i++) {
    706                 pair = {X: x[i], Y: y[i]};
    707                 data.push(pair);
    708             }
    709             data.sort(function (a, b) {
    710                 return a.X - b.X;
    711             });
    712             for (i = 0; i < n; i++) {
    713                 x[i] = data[i].X;
    714                 y[i] = data[i].Y;
    715             }
    716
    717             for (i = 0; i < n - 1; i++) {
    718                 dx.push(x[i + 1] - x[i]);
    719             }
    720             for (i = 0; i < n - 2; i++) {
    721                 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i]));
    722             }
    723
    724             // ForwardSolve
    725             diag.push(2 * (dx[0] + dx[1]));
    726             z.push(delta[0]);
    727
    728             for (i = 0; i < n - 3; i++) {
    729                 l = dx[i + 1] / diag[i];
    730                 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]);
    731                 z.push(delta[i + 1] - l * z[i]);
    732             }
    733
    734             // BackwardSolve
    735             F[n - 3] = z[n - 3] / diag[n - 3];
    736             for (i = n - 4; i >= 0; i--) {
    737                 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i];
    738             }
    739
    740             // Generate f''-Vector
    741             for (i = n - 3; i >= 0; i--) {
    742                 F[i + 1] = F[i];
    743             }
    744
    745             // natural cubic spline
    746             F[0] = 0;
    747             F[n - 1] = 0;
    748
    749             return F;
    750         },
    751
    752         /**
    753          * Evaluate points on spline.
    754          * @param {Number,Array} x0 A single float value or an array of values to evaluate
    755          * @param {Array} x x values of knots
    756          * @param {Array} y y values of knots
    757          * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef}
    758          * @see #splineDef
    759          * @returns {Number,Array} A single value or an array, depending on what is given as x0.
    760          * @memberof JXG.Math.Numerics
    761          */
    762         splineEval: function (x0, x, y, F) {
    763             var i, j, a, b, c, d, x_,
    764                 n = Math.min(x.length, y.length),
    765                 l = 1,
    766                 asArray = false,
    767                 y0 = [];
    768
    769             // number of points to be evaluated
    770             if (Type.isArray(x0)) {
    771                 l = x0.length;
    772                 asArray = true;
    773             } else {
    774                 x0 = [x0];
    775             }
    776
    777             for (i = 0; i < l; i++) {
    778                 // is x0 in defining interval?
    779                 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) {
    780                     return NaN;
    781                 }
    782
    783                 // determine part of spline in which x0 lies
    784                 for (j = 1; j < n; j++) {
    785                     if (x0[i] <= x[j]) {
    786                         break;
    787                     }
    788                 }
    789
    790                 j -= 1;
    791
    792                 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1];
    793                 // determine the coefficients of the polynomial in this interval
    794                 a = y[j];
    795                 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]);
    796                 c = F[j] / 2;
    797                 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j]));
    798                 // evaluate x0[i]
    799                 x_ = x0[i] - x[j];
    800                 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_);
    801                 y0.push(a + (b + (c + d * x_) * x_) * x_);
    802             }
    803
    804             if (asArray) {
    805                 return y0;
    806             }
    807
    808             return y0[0];
    809         },
    810
    811         /**
    812          * Generate a string containing the function term of a polynomial.
    813          * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i.
    814          * @param {Number} deg Degree of the polynomial
    815          * @param {String} varname Name of the variable (usually 'x')
    816          * @param {Number} prec Precision
    817          * @returns {String} A string containg the function term of the polynomial.
    818          * @memberof JXG.Math.Numerics
    819          */
    820         generatePolynomialTerm: function (coeffs, deg, varname, prec) {
    821             var i, t = [];
    822
    823             for (i = deg; i >= 0; i--) {
    824                 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']);
    825
    826                 if (i > 1) {
    827                     t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']);
    828                 } else if (i === 1) {
    829                     t = t.concat(['*', varname, ' + ']);
    830                 }
    831             }
    832
    833             return t.join('');
    834         },
    835
    836         /**
    837          * Computes the polynomial through a given set of coordinates in Lagrange form.
    838          * Returns the Lagrange polynomials, see
    839          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
    840          * SIAM Review, Vol 46, No 3, (2004) 501-517.
    841          * @param {Array} p Array of JXG.Points
    842          * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points.
    843          * @memberof JXG.Math.Numerics
    844          */
    845         lagrangePolynomial: function (p) {
    846             var w = [],
    847                 /** @ignore */
    848                 fct = function (x, suspendedUpdate) {
    849                     var i, j, k, xi, s, M,
    850                         len = p.length,
    851                         num = 0,
    852                         denom = 0;
    853
    854                     if (!suspendedUpdate) {
    855                         for (i = 0; i < len; i++) {
    856                             w[i] = 1.0;
    857                             xi = p[i].X();
    858
    859                             for (k = 0; k < len; k++) {
    860                                 if (k !== i) {
    861                                     w[i] *= (xi - p[k].X());
    862                                 }
    863                             }
    864
    865                             w[i] = 1 / w[i];
    866                         }
    867
    868                         M = [];
    869
    870                         for (j = 0; j < len; j++) {
    871                             M.push([1]);
    872                         }
    873                     }
    874
    875                     for (i = 0; i < len; i++) {
    876                         xi = p[i].X();
    877
    878                         if (x === xi) {
    879                             return p[i].Y();
    880                         }
    881
    882                         s = w[i] / (x - xi);
    883                         denom += s;
    884                         num += s * p[i].Y();
    885                     }
    886
    887                     return num / denom;
    888                 };
    889
    890             fct.getTerm = function () {
    891                 return '';
    892             };
    893
    894             return fct;
    895         },
    896
    897         /**
    898          * Computes the cubic cardinal spline curve through a given set of points. The curve
    899          * is uniformly parametrized.
    900          * Two artificial control points at the beginning and the end are added.
    901          * @param {Array} points Array consisting of JXG.Points.
    902          * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1.
    903          * tau=1/2 give Catmull-Rom splines.
    904          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
    905          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
    906          * returning the length of the points array minus three.
    907          * @memberof JXG.Math.Numerics
    908         */
    909         CardinalSpline: function (points, tau) {
    910             var p,
    911                 coeffs = [],
    912                 // control point at the beginning and at the end
    913                 first = {},
    914                 last = {},
    915                 makeFct,
    916                 _tau;
    917
    918             if (Type.isFunction(tau)) {
    919                 _tau = tau;
    920             } else {
    921                 _tau = function () { return tau; };
    922             }
    923
    924             /** @ignore */
    925             makeFct = function (which) {
    926                 return function (t, suspendedUpdate) {
    927                     var s, c,
    928                         len = points.length,
    929                         tau = _tau();
    930
    931                     if (len < 2) {
    932                         return NaN;
    933                     }
    934
    935                     if (!suspendedUpdate) {
    936                         first[which] = function () {
    937                             return 2 * points[0][which]() - points[1][which]();
    938                         };
    939
    940                         last[which] = function () {
    941                             return 2 * points[len - 1][which]() - points[len - 2][which]();
    942                         };
    943
    944                         p = [first].concat(points, [last]);
    945                         coeffs[which] = [];
    946
    947                         for (s = 0; s < len - 1; s++) {
    948                             coeffs[which][s] = [
    949                                 1 / tau * p[s + 1][which](),
    950                                 -p[s][which]() +   p[s + 2][which](),
    951                                 2 * p[s][which]() + (-3 / tau + 1) * p[s + 1][which]() + (3 / tau - 2) * p[s + 2][which]() - p[s + 3][which](),
    952                                 -p[s][which]() + (2 / tau - 1) * p[s + 1][which]() + (-2 / tau + 1) * p[s + 2][which]() + p[s + 3][which]()
    953                             ];
    954                         }
    955                     }
    956
    957                     len += 2;  // add the two control points
    958
    959                     if (isNaN(t)) {
    960                         return NaN;
    961                     }
    962
    963                     // This is necessary for our advanced plotting algorithm:
    964                     if (t <= 0.0) {
    965                         return p[1][which]();
    966                     }
    967
    968                     if (t >= len - 3) {
    969                         return p[len - 2][which]();
    970                     }
    971
    972                     s = Math.floor(t);
    973
    974                     if (s === t) {
    975                         return p[s][which]();
    976                     }
    977
    978                     t -= s;
    979                     c = coeffs[which][s];
    980
    981                     return tau * (((c[3] * t + c[2]) * t + c[1]) * t + c[0]);
    982                 };
    983             };
    984
    985             return [makeFct('X'), makeFct('Y'), 0,
    986                 function () {
    987                     return points.length - 1;
    988                 }];
    989         },
    990
    991         /**
    992          * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve
    993          * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5.
    994          * Two artificial control points at the beginning and the end are added.
    995          * @param {Array} points Array consisting of JXG.Points.
    996          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
    997          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
    998          * returning the length of the points array minus three.
    999          * @memberof JXG.Math.Numerics
    1000         */
    1001         CatmullRomSpline: function (points) {
    1002             return this.CardinalSpline(points, 0.5);
    1003         },
    1004
    1005         /**
    1006          * Computes the regression polynomial of a given degree through a given set of coordinates.
    1007          * Returns the regression polynomial function.
    1008          * @param {Number,function,Slider} degree number, function or slider.
    1009          * Either
    1010          * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in
    1011          * an array of {@link JXG.Point}s or {@link JXG.Coords}. In the latter case, the <tt>dataY</tt> parameter will be ignored.
    1012          * @param {Array} dataY Array containing the y-coordinates of the data set,
    1013          * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree.
    1014          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
    1015          * @memberof JXG.Math.Numerics
    1016          */
    1017         regressionPolynomial: function (degree, dataX, dataY) {
    1018             var coeffs, deg, dX, dY, inputType, fct,
    1019                 term = '';
    1020
    1021             // Slider
    1022             if (Type.isPoint(degree) && typeof degree.Value === 'function') {
    1023                 /** @ignore */
    1024                 deg = function () {
    1025                     return degree.Value();
    1026                 };
    1027             // function
    1028             } else if (Type.isFunction(degree)) {
    1029                 deg = degree;
    1030             // number
    1031             } else if (Type.isNumber(degree)) {
    1032                 /** @ignore */
    1033                 deg = function () {
    1034                     return degree;
    1035                 };
    1036             } else {
    1037                 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'.");
    1038             }
    1039
    1040             // Parameters degree, dataX, dataY
    1041             if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) {
    1042                 inputType = 0;
    1043             // Parameters degree, point array
    1044             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) {
    1045                 inputType = 1;
    1046             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) {
    1047                 inputType = 2;
    1048             } else {
    1049                 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters.");
    1050             }
    1051
    1052             /** @ignore */
    1053             fct = function (x, suspendedUpdate) {
    1054                 var i, j, M, MT, y, B, c, s, d,
    1055                     // input data
    1056                     len = dataX.length;
    1057
    1058                 d = Math.floor(deg());
    1059
    1060                 if (!suspendedUpdate) {
    1061                     // point list as input
    1062                     if (inputType === 1) {
    1063                         dX = [];
    1064                         dY = [];
    1065
    1066                         for (i = 0; i < len; i++) {
    1067                             dX[i] = dataX[i].X();
    1068                             dY[i] = dataX[i].Y();
    1069                         }
    1070                     }
    1071
    1072                     if (inputType === 2) {
    1073                         dX = [];
    1074                         dY = [];
    1075
    1076                         for (i = 0; i < len; i++) {
    1077                             dX[i] = dataX[i].usrCoords[1];
    1078                             dY[i] = dataX[i].usrCoords[2];
    1079                         }
    1080                     }
    1081
    1082                     // check for functions
    1083                     if (inputType === 0) {
    1084                         dX = [];
    1085                         dY = [];
    1086
    1087                         for (i = 0; i < len; i++) {
    1088                             if (Type.isFunction(dataX[i])) {
    1089                                 dX.push(dataX[i]());
    1090                             } else {
    1091                                 dX.push(dataX[i]);
    1092                             }
    1093
    1094                             if (Type.isFunction(dataY[i])) {
    1095                                 dY.push(dataY[i]());
    1096                             } else {
    1097                                 dY.push(dataY[i]);
    1098                             }
    1099                         }
    1100                     }
    1101
    1102                     M = [];
    1103
    1104                     for (j = 0; j < len; j++) {
    1105                         M.push([1]);
    1106                     }
    1107
    1108                     for (i = 1; i <= d; i++) {
    1109                         for (j = 0; j < len; j++) {
    1110                             M[j][i] = M[j][i - 1] * dX[j];
    1111                         }
    1112                     }
    1113
    1114                     y = dY;
    1115                     MT = Mat.transpose(M);
    1116                     B = Mat.matMatMult(MT, M);
    1117                     c = Mat.matVecMult(MT, y);
    1118                     coeffs = Mat.Numerics.Gauss(B, c);
    1119                     term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3);
    1120                 }
    1121
    1122                 // Horner's scheme to evaluate polynomial
    1123                 s = coeffs[d];
    1124
    1125                 for (i = d - 1; i >= 0; i--) {
    1126                     s = (s * x + coeffs[i]);
    1127                 }
    1128
    1129                 return s;
    1130             };
    1131
    1132             fct.getTerm = function () {
    1133                 return term;
    1134             };
    1135
    1136             return fct;
    1137         },
    1138
    1139         /**
    1140          * Computes the cubic Bezier curve through a given set of points.
    1141          * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}.
    1142          * The points at position k with k mod 3 = 0 are the data points,
    1143          * points at position k with k mod 3 = 1 or 2 are the control points.
    1144          * @returns {Array} An array consisting of two functions of one parameter t which return the
    1145          * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting
    1146          * no parameters and returning one third of the length of the points.
    1147          * @memberof JXG.Math.Numerics
    1148          */
    1149         bezier: function (points) {
    1150             var len, flen,
    1151                 /** @ignore */
    1152                 makeFct = function (which) {
    1153                     return function (t, suspendedUpdate) {
    1154                         var z = Math.floor(t) * 3,
    1155                             t0 = t % 1,
    1156                             t1 = 1 - t0;
    1157
    1158                         if (!suspendedUpdate) {
    1159                             flen = 3 * Math.floor((points.length - 1) / 3);
    1160                             len = Math.floor(flen / 3);
    1161                         }
    1162
    1163                         if (t < 0) {
    1164                             return points[0][which]();
    1165                         }
    1166
    1167                         if (t >= len) {
    1168                             return points[flen][which]();
    1169                         }
    1170
    1171                         if (isNaN(t)) {
    1172                             return NaN;
    1173                         }
    1174
    1175                         return t1 * t1 * (t1 * points[z][which]() + 3 * t0 * points[z + 1][which]()) + (3 * t1 * points[z + 2][which]() + t0 * points[z + 3][which]()) * t0 * t0;
    1176                     };
    1177                 };
    1178
    1179             return [makeFct('X'), makeFct('Y'), 0,
    1180                 function () {
    1181                     return Math.floor(points.length / 3);
    1182                 }];
    1183         },
    1184
    1185         /**
    1186          * Computes the B-spline curve of order k (order = degree+1) through a given set of points.
    1187          * @param {Array} points Array consisting of JXG.Points.
    1188          * @param {Number} order Order of the B-spline curve.
    1189          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
    1190          * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply
    1191          * returning the length of the points array minus one.
    1192          * @memberof JXG.Math.Numerics
    1193          */
    1194         bspline: function (points, order) {
    1195             var knots, N = [],
    1196                 _knotVector = function (n, k) {
    1197                     var j,
    1198                         kn = [];
    1199
    1200                     for (j = 0; j < n + k + 1; j++) {
    1201                         if (j < k) {
    1202                             kn[j] = 0.0;
    1203                         } else if (j <= n) {
    1204                             kn[j] = j - k + 1;
    1205                         } else {
    1206                             kn[j] = n - k + 2;
    1207                         }
    1208                     }
    1209
    1210                     return kn;
    1211                 },
    1212
    1213                 _evalBasisFuncs = function (t, kn, n, k, s) {
    1214                     var i, j, a, b, den,
    1215                         N = [];
    1216
    1217                     if (kn[s] <= t && t < kn[s + 1]) {
    1218                         N[s] = 1;
    1219                     } else {
    1220                         N[s] = 0;
    1221                     }
    1222
    1223                     for (i = 2; i <= k; i++) {
    1224                         for (j = s - i + 1; j <= s; j++) {
    1225                             if (j <= s - i + 1 || j < 0) {
    1226                                 a = 0.0;
    1227                             } else {
    1228                                 a = N[j];
    1229                             }
    1230
    1231                             if (j >= s) {
    1232                                 b = 0.0;
    1233                             } else {
    1234                                 b = N[j + 1];
    1235                             }
    1236
    1237                             den = kn[j + i - 1] - kn[j];
    1238
    1239                             if (den === 0) {
    1240                                 N[j] = 0;
    1241                             } else {
    1242                                 N[j] = (t - kn[j]) / den * a;
    1243                             }
    1244
    1245                             den = kn[j + i] - kn[j + 1];
    1246
    1247                             if (den !== 0) {
    1248                                 N[j] += (kn[j + i] - t) / den * b;
    1249                             }
    1250                         }
    1251                     }
    1252                     return N;
    1253                 },
    1254                 /** @ignore */
    1255                 makeFct = function (which) {
    1256                     return function (t, suspendedUpdate) {
    1257                         var y, j, s,
    1258                             len = points.length,
    1259                             n = len - 1,
    1260                             k = order;
    1261
    1262                         if (n <= 0) {
    1263                             return NaN;
    1264                         }
    1265
    1266                         if (n + 2 <= k) {
    1267                             k = n + 1;
    1268                         }
    1269
    1270                         if (t <= 0) {
    1271                             return points[0][which]();
    1272                         }
    1273
    1274                         if (t >= n - k + 2) {
    1275                             return points[n][which]();
    1276                         }
    1277
    1278                         s = Math.floor(t) + k - 1;
    1279                         knots = _knotVector(n, k);
    1280                         N = _evalBasisFuncs(t, knots, n, k, s);
    1281
    1282                         y = 0.0;
    1283                         for (j = s - k + 1; j <= s; j++) {
    1284                             if (j < len && j >= 0) {
    1285                                 y += points[j][which]() * N[j];
    1286                             }
    1287                         }
    1288
    1289                         return y;
    1290                     };
    1291                 };
    1292
    1293             return [makeFct('X'), makeFct('Y'), 0,
    1294                 function () {
    1295                     return points.length - 1;
    1296                 }];
    1297         },
    1298
    1299         /**
    1300          * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve}
    1301          * and {@link JXG.Curve#hasPoint}.
    1302          * @param {function} f Function in one variable to be differentiated.
    1303          * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a
    1304          * method of an object and contains a reference to its parent object via "this".
    1305          * @returns {function} Derivative function of a given function f.
    1306          * @memberof JXG.Math.Numerics
    1307          */
    1308         D: function (f, obj) {
    1309             var h = 0.00001,
    1310                 h2 = 1.0 / (h * 2.0);
    1311
    1312             if (!Type.exists(obj)) {
    1313                 return function (x, suspendUpdate) {
    1314                     return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2;
    1315                 };
    1316             }
    1317
    1318             return function (x, suspendUpdate) {
    1319                 return (f.apply(obj, [x + h, suspendUpdate]) - f.apply(obj, [x - h, suspendUpdate])) * h2;
    1320             };
    1321         },
    1322
    1323         /**
    1324          * Helper function to create curve which displays Riemann sums.
    1325          * Compute coordinates for the rectangles showing the Riemann sum.
    1326          * @param {function} f Function, whose integral is approximated by the Riemann sum.
    1327          * @param {Number} n number of rectangles.
    1328          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'.
    1329          * @param {Number} start Left border of the approximation interval
    1330          * @param {Number} end Right border of the approximation interval
    1331          * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This
    1332          * array may be used as parent array of a {@link JXG.Curve}. The third parameteris the riemann sum, i.e. the sum of the volumes of all
    1333          * rectangles.
    1334          * @memberof JXG.Math.Numerics
    1335          */
    1336         riemann: function (f, n, type, start, end) {
    1337             var i, x1, y1, delta1, delta,
    1338                 xarr = [],
    1339                 yarr = [],
    1340                 j = 0,
    1341                 x = start, y,
    1342                 sum = 0;
    1343
    1344             n = Math.round(n);
    1345
    1346             xarr[j] = x;
    1347             yarr[j] = 0.0;
    1348
    1349             if (n > 0) {
    1350                 delta = (end - start) / n;
    1351                 // for 'lower' and 'upper'
    1352                 delta1 = delta * 0.01;
    1353
    1354                 for (i = 0; i < n; i++) {
    1355                     if (type === 'right') {
    1356                         y = f(x + delta);
    1357                     } else if (type === 'middle') {
    1358                         y = f(x + delta * 0.5);
    1359                     } else if (type === 'left' || type === 'trapezoidal') {
    1360                         y = f(x);
    1361                     } else if (type === 'lower') {
    1362                         y = f(x);
    1363
    1364                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
    1365                             y1 = f(x1);
    1366
    1367                             if (y1 < y) {
    1368                                 y = y1;
    1369                             }
    1370                         }
    1371
    1372                         y1 = f(x + delta);
    1373                         if (y1 < y) {
    1374                             y = y1;
    1375                         }
    1376                     } else if (type === 'upper') {
    1377                         y = f(x);
    1378
    1379                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
    1380                             y1 = f(x1);
    1381
    1382                             if (y1 > y) {
    1383                                 y = y1;
    1384                             }
    1385                         }
    1386
    1387                         y1 = f(x + delta);
    1388                         if (y1 > y) {
    1389                             y = y1;
    1390                         }
    1391                     } else if (type === 'random') {
    1392                         y = f(x + delta * Math.random());
    1393                     } else if (type === 'simpson') {
    1394                         y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0;
    1395                     } else {
    1396                         y = f(x);  // default is lower
    1397                     }
    1398
    1399                     j += 1;
    1400                     xarr[j] = x;
    1401                     yarr[j] = y;
    1402                     j += 1;
    1403                     x += delta;
    1404
    1405                     if (type === 'trapezoidal') {
    1406                         y = f(x);
    1407                     }
    1408
    1409                     xarr[j] = x;
    1410                     yarr[j] = y;
    1411                     j += 1;
    1412                     xarr[j] = x;
    1413                     yarr[j] = 0.0;
    1414                     sum += y * delta;
    1415                 }
    1416             }
    1417             return [xarr, yarr, sum];
    1418         },
    1419
    1420         /**
    1421          * Approximate the integral by Riemann sums.
    1422          * Compute the area described by the riemann sum rectangles.
    1423          * @deprecated Replaced by JXG.Curve.Value(), see {@link JXG.Curve#riemannsum}
    1424          * @param {function} f Function, whose integral is approximated by the Riemann sum.
    1425          * @param {Number} n number of rectangles.
    1426          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'.
    1427          *
    1428          * @param {Number} start Left border of the approximation interval
    1429          * @param {Number} end Right border of the approximation interval
    1430          * @returns {Number} The sum of the areas of the rectangles.
    1431          * @memberof JXG.Math.Numerics
    1432          */
    1433         riemannsum: function (f, n, type, start, end) {
    1434             var i, x1, y1, delta1, delta, y,
    1435                 sum = 0.0,
    1436                 x = start;
    1437
    1438             // this function looks very similar to this.riemann... maybe there is some merge potential?
    1439
    1440             n = Math.floor(n);
    1441
    1442             if (n > 0) {
    1443                 delta = (end - start) / n;
    1444                 // for 'lower' and 'upper'
    1445                 delta1 = delta * 0.01;
    1446
    1447                 for (i = 0; i < n; i++) {
    1448                     if (type === 'right') {
    1449                         y = f(x + delta);
    1450                     } else if (type === 'middle') {
    1451                         y = f(x + delta * 0.5);
    1452                     } else if (type === 'trapezoidal') {
    1453                         y = 0.5 * (f(x + delta) + f(x));
    1454                     } else if (type === 'left') {
    1455                         y = f(x);
    1456                     } else if (type === 'lower') {
    1457                         y = f(x);
    1458                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
    1459                             y1 = f(x1);
    1460
    1461                             if (y1 < y) {
    1462                                 y = y1;
    1463                             }
    1464                         }
    1465
    1466                         y1 = f(x + delta);
    1467                         if (y1 < y) {
    1468                             y = y1;
    1469                         }
    1470                     } else if (type === 'upper') {
    1471                         y = f(x);
    1472
    1473                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
    1474                             y1 = f(x1);
    1475
    1476                             if (y1 > y) {
    1477                                 y = y1;
    1478                             }
    1479                         }
    1480
    1481                         y1 = f(x + delta);
    1482                         if (y1 > y) {
    1483                             y = y1;
    1484                         }
    1485                     } else if (type === 'random') {
    1486                         y = f(x + delta * Math.random());
    1487                     } else if (type === 'simpson') {
    1488                         y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0;
    1489                     } else {
    1490                         y = f(x);  // default is lower
    1491                     }
    1492
    1493                     sum += delta * y;
    1494                     x += delta;
    1495                 }
    1496             }
    1497
    1498             return sum;
    1499         },
    1500
    1501         /**
    1502          * Solve initial value problems numerically using Runge-Kutta-methods.
    1503          * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
    1504          * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
    1505          * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
    1506          * <pre>
    1507          * {
    1508          *     s: <Number>,
    1509          *     A: <matrix>,
    1510          *     b: <Array>,
    1511          *     c: <Array>
    1512          * }
    1513          * </pre>
    1514          * which corresponds to the Butcher tableau structure shown here: http://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696
    1515          * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array.
    1516          * @param {Array} I Interval on which to integrate.
    1517          * @param {Number} N Number of evaluation points.
    1518          * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
    1519          * is given by the equation <pre>dx/dt = f(t, x(t)).</pre> So f has to take two parameters, a number <tt>t</tt> and a
    1520          * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has.
    1521          * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
    1522          * @example
    1523          * // A very simple autonomous system dx(t)/dt = x(t);
    1524          * function f(t, x) {
    1525          *     return x;
    1526          * }
    1527          *
    1528          * // Solve it with initial value x(0) = 1 on the interval [0, 2]
    1529          * // with 20 evaluation points.
    1530          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
    1531          *
    1532          * // Prepare data for plotting the solution of the ode using a curve.
    1533          * var dataX = [];
    1534          * var dataY = [];
    1535          * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
    1536          * for(var i=0; i<data.length; i++) {
    1537          *     dataX[i] = i*h;
    1538          *     dataY[i] = data[i][0];
    1539          * }
    1540          * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
    1541          * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
    1542          * <script type="text/javascript">
    1543          * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
    1544          * function f(t, x) {
    1545          *     // we have to copy the value.
    1546          *     // return x; would just return the reference.
    1547          *     return [x[0]];
    1548          * }
    1549          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
    1550          * var dataX = [];
    1551          * var dataY = [];
    1552          * var h = 0.1;
    1553          * for(var i=0; i<data.length; i++) {
    1554          *     dataX[i] = i*h;
    1555          *     dataY[i] = data[i][0];
    1556          * }
    1557          * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
    1558          * </script><pre>
    1559          * @memberof JXG.Math.Numerics
    1560          */
    1561         rungeKutta: function (butcher, x0, I, N, f) {
    1562             var e, i, j, k, l, s,
    1563                 x = [],
    1564                 y = [],
    1565                 h = (I[1] - I[0]) / N,
    1566                 t = I[0],
    1567                 dim = x0.length,
    1568                 result = [],
    1569                 r = 0;
    1570
    1571             if (Type.isString(butcher)) {
    1572                 butcher = predefinedButcher[butcher] || predefinedButcher.euler;
    1573             }
    1574             s = butcher.s;
    1575
    1576             // don't change x0, so copy it
    1577             for (e = 0; e < dim; e++) {
    1578                 x[e] = x0[e];
    1579             }
    1580
    1581             for (i = 0; i < N; i++) {
    1582                 // Optimization doesn't work for ODEs plotted using time
    1583                 //        if((i % quotient == 0) || (i == N-1)) {
    1584                 result[r] = [];
    1585                 for (e = 0; e < dim; e++) {
    1586                     result[r][e] = x[e];
    1587                 }
    1588
    1589                 r += 1;
    1590                 k = [];
    1591
    1592                 for (j = 0; j < s; j++) {
    1593                     // init y = 0
    1594                     for (e = 0; e < dim; e++) {
    1595                         y[e] = 0.0;
    1596                     }
    1597
    1598
    1599                     // Calculate linear combination of former k's and save it in y
    1600                     for (l = 0; l < j; l++) {
    1601                         for (e = 0; e < dim; e++) {
    1602                             y[e] += (butcher.A[j][l]) * h * k[l][e];
    1603                         }
    1604                     }
    1605
    1606                     // add x(t) to y
    1607                     for (e = 0; e < dim; e++) {
    1608                         y[e] += x[e];
    1609                     }
    1610
    1611                     // calculate new k and add it to the k matrix
    1612                     k.push(f(t + butcher.c[j] * h, y));
    1613                 }
    1614
    1615                 // init y = 0
    1616                 for (e = 0; e < dim; e++) {
    1617                     y[e] = 0.0;
    1618                 }
    1619
    1620                 for (l = 0; l < s; l++) {
    1621                     for (e = 0; e < dim; e++) {
    1622                         y[e] += butcher.b[l] * k[l][e];
    1623                     }
    1624                 }
    1625
    1626                 for (e = 0; e < dim; e++) {
    1627                     x[e] = x[e] + h * y[e];
    1628                 }
    1629
    1630                 t += h;
    1631             }
    1632
    1633             return result;
    1634         },
    1635
    1636
    1637         /**
    1638          * Maximum number of iterations in {@link JXG.Math.Numerics#fzero}
    1639          * @type Number
    1640          * @default 80
    1641          * @memberof JXG.Math.Numerics
    1642          */
    1643         maxIterationsRoot: 80,
    1644
    1645         /**
    1646          * Maximum number of iterations in {@link JXG.Math.Numerics#fminbr}
    1647          * @type Number
    1648          * @default 500
    1649          * @memberof JXG.Math.Numerics
    1650          */
    1651         maxIterationsMinimize: 500,
    1652
    1653         /**
    1654          *
    1655          * Find zero of an univariate function f.
    1656          * @param {function} f Function, whose root is to be found
    1657          * @param {Array,Number} x0  Start value or start interval enclosing the root
    1658          * @param {Object} object Parent object in case f is method of it
    1659          * @returns {Number} the approximation of the root
    1660          * Algorithm:
    1661          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
    1662          *  computations. M., Mir, 1980, p.180 of the Russian edition
    1663          *
    1664          * If x0 is an array containing lower and upper bound for the zero
    1665          * algorithm 748 is applied. Otherwise, if x0 is a number,
    1666          * the algorithm tries to bracket a zero of f starting from x0.
    1667          * If this fails, we fall back to Newton's method.
    1668          * @memberof JXG.Math.Numerics
    1669          */
    1670         fzero: function (f, x0, object) {
    1671             var a, b, c,
    1672                 fa, fb, fc,
    1673                 aa, blist, i, len, u, fu,
    1674                 prev_step, t1, cb, t2,
    1675                 // Actual tolerance
    1676                 tol_act,
    1677                 // Interpolation step is calculated in the form p/q; division
    1678                 // operations is delayed until the last moment
    1679                 p, q,
    1680                 // Step at this iteration
    1681                 new_step,
    1682                 eps = Mat.eps,
    1683                 maxiter = this.maxIterationsRoot,
    1684                 niter = 0,
    1685                 nfev = 0;
    1686
    1687             if (Type.isArray(x0)) {
    1688                 if (x0.length < 2) {
    1689                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
    1690                 }
    1691
    1692                 a = x0[0];
    1693                 fa = f.call(object, a);
    1694                 nfev += 1;
    1695                 b = x0[1];
    1696                 fb = f.call(object, b);
    1697                 nfev += 1;
    1698             } else {
    1699                 a = x0;
    1700                 fa = f.call(object, a);
    1701                 nfev += 1;
    1702
    1703                 // Try to get b.
    1704                 if (a === 0) {
    1705                     aa = 1;
    1706                 } else {
    1707                     aa = a;
    1708                 }
    1709
    1710                 blist = [0.9 * aa, 1.1 * aa, aa - 1, aa + 1, 0.5 * aa, 1.5 * aa, -aa, 2 * aa, -10 * aa, 10 * aa];
    1711                 len = blist.length;
    1712
    1713                 for (i = 0; i < len; i++) {
    1714                     b = blist[i];
    1715                     fb = f.call(object, b);
    1716                     nfev += 1;
    1717
    1718                     if (fa * fb <= 0) {
    1719                         break;
    1720                     }
    1721                 }
    1722                 if (b < a) {
    1723                     u = a;
    1724                     a = b;
    1725                     b = u;
    1726
    1727                     fu = fa;
    1728                     fa = fb;
    1729                     fb = fu;
    1730                 }
    1731             }
    1732
    1733             if (fa * fb > 0) {
    1734                 // Bracketing not successful, fall back to Newton's method or to fminbr
    1735                 if (Type.isArray(x0)) {
    1736                     return this.fminbr(f, [a, b], object);
    1737                 }
    1738
    1739                 return this.Newton(f, a, object);
    1740             }
    1741
    1742             // OK, we have enclosed a zero of f.
    1743             // Now we can start Brent's method
    1744
    1745             c = a;
    1746             fc = fa;
    1747
    1748             // Main iteration loop
    1749             while (niter < maxiter) {
    1750                 // Distance from the last but one to the last approximation
    1751                 prev_step = b - a;
    1752
    1753                 // Swap data for b to be the best approximation
    1754                 if (Math.abs(fc) < Math.abs(fb)) {
    1755                     a = b;
    1756                     b = c;
    1757                     c = a;
    1758
    1759                     fa = fb;
    1760                     fb = fc;
    1761                     fc = fa;
    1762                 }
    1763                 tol_act = 2 * eps * Math.abs(b) + eps * 0.5;
    1764                 new_step = (c - b) * 0.5;
    1765
    1766                 if (Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps) {
    1767                     //  Acceptable approx. is found
    1768                     return b;
    1769                 }
    1770
    1771                 // Decide if the interpolation can be tried
    1772                 // If prev_step was large enough and was in true direction Interpolatiom may be tried
    1773                 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
    1774                     cb = c - b;
    1775
    1776                     // If we have only two distinct points linear interpolation can only be applied
    1777                     if (a === c) {
    1778                         t1 = fb / fa;
    1779                         p = cb * t1;
    1780                         q = 1.0 - t1;
    1781                     // Quadric inverse interpolation
    1782                     } else {
    1783                         q = fa / fc;
    1784                         t1 = fb / fc;
    1785                         t2 = fb / fa;
    1786
    1787                         p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
    1788                         q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
    1789                     }
    1790
    1791                     // p was calculated with the opposite sign; make p positive
    1792                     if (p > 0) {
    1793                         q = -q;
    1794                     // and assign possible minus to q
    1795                     } else {
    1796                         p = -p;
    1797                     }
    1798
    1799                     // If b+p/q falls in [b,c] and isn't too large it is accepted
    1800                     // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
    1801                     if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) &&
    1802                             p < Math.abs(prev_step * q * 0.5)) {
    1803                         new_step = p / q;
    1804                     }
    1805                 }
    1806
    1807                 // Adjust the step to be not less than tolerance
    1808                 if (Math.abs(new_step) < tol_act) {
    1809                     if (new_step > 0) {
    1810                         new_step = tol_act;
    1811                     } else {
    1812                         new_step = -tol_act;
    1813                     }
    1814                 }
    1815
    1816                 // Save the previous approx.
    1817                 a = b;
    1818                 fa = fb;
    1819                 b += new_step;
    1820                 fb = f.call(object, b);
    1821                 // Do step to a new approxim.
    1822                 nfev += 1;
    1823
    1824                 // Adjust c for it to have a sign opposite to that of b
    1825                 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
    1826                     c = a;
    1827                     fc = fa;
    1828                 }
    1829                 niter++;
    1830             } // End while
    1831
    1832             return b;
    1833         },
    1834
    1835
    1836         /**
    1837          *
    1838          * Find minimum of an univariate function f.
    1839          * @param {function} f Function, whose minimum is to be found
    1840          * @param {Array} x0  Start interval enclosing the minimum
    1841          * @param {Object} context Parent object in case f is method of it
    1842          * @returns {Number} the approximation of the minimum value position
    1843          * Algorithm:
    1844          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
    1845          *  computations. M., Mir, 1980, p.180 of the Russian edition
    1846          * x0
    1847          * @memberof JXG.Math.Numerics
    1848          **/
    1849         fminbr: function (f, x0, context) {
    1850             var a, b, x, v, w,
    1851                 fx, fv, fw,
    1852                 range, middle_range, tol_act, new_step,
    1853                 p, q, t, ft,
    1854                 // Golden section ratio
    1855                 r = (3.0 - Math.sqrt(5.0)) * 0.5,
    1856                 tol = Mat.eps,
    1857                 sqrteps = Mat.eps, //Math.sqrt(Mat.eps),
    1858                 maxiter = this.maxIterationsMinimize,
    1859                 niter = 0,
    1860                 nfev = 0;
    1861
    1862             if (!Type.isArray(x0) || x0.length < 2) {
    1863                 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two.");
    1864             }
    1865
    1866             a = x0[0];
    1867             b = x0[1];
    1868             v = a + r * (b - a);
    1869             fv = f.call(context, v);
    1870
    1871             // First step - always gold section
    1872             nfev += 1;
    1873             x = v;
    1874             w = v;
    1875             fx = fv;
    1876             fw = fv;
    1877
    1878             while (niter < maxiter) {
    1879                 // Range over the interval in which we are looking for the minimum
    1880                 range = b - a;
    1881                 middle_range = (a + b) * 0.5;
    1882
    1883                 // Actual tolerance
    1884                 tol_act = sqrteps * Math.abs(x) + tol / 3.0;
    1885
    1886                 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) {
    1887                     // Acceptable approx. is found
    1888                     return x;
    1889                 }
    1890
    1891                 // Obtain the golden section step
    1892                 new_step = r * (x < middle_range ? b - x : a - x);
    1893
    1894                 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried
    1895                 if (Math.abs(x - w) >= tol_act) {
    1896                     // Interpolation step is calculated as p/q;
    1897                     // division operation is delayed until last moment
    1898                     t = (x - w) * (fx - fv);
    1899                     q = (x - v) * (fx - fw);
    1900                     p = (x - v) * q - (x - w) * t;
    1901                     q = 2 * (q - t);
    1902
    1903                     if (q > 0) {                        // q was calculated with the op-
    1904                         p = -p;                         // posite sign; make q positive
    1905                     } else {                            // and assign possible minus to
    1906                         q = -q;                         // p
    1907                     }
    1908                     if (Math.abs(p) < Math.abs(new_step * q) &&     // If x+p/q falls in [a,b]
    1909                             p > q * (a - x + 2 * tol_act) &&        //  not too close to a and
    1910                             p < q * (b - x - 2 * tol_act)) {        // b, and isn't too large
    1911                         new_step = p / q;                          // it is accepted
    1912                     }
    1913                     // If p/q is too large then the
    1914                     // golden section procedure can
    1915                     // reduce [a,b] range to more
    1916                     // extent
    1917                 }
    1918
    1919                 // Adjust the step to be not less than tolerance
    1920                 if (Math.abs(new_step) < tol_act) {
    1921                     if (new_step > 0) {
    1922                         new_step = tol_act;
    1923                     } else {
    1924                         new_step = -tol_act;
    1925                     }
    1926                 }
    1927
    1928                 // Obtain the next approximation to min
    1929                 // and reduce the enveloping range
    1930
    1931                 // Tentative point for the min
    1932                 t = x + new_step;
    1933                 ft = f.call(context, t);
    1934                 nfev += 1;
    1935
    1936                 // t is a better approximation
    1937                 if (ft <= fx) {
    1938                     // Reduce the range so that t would fall within it
    1939                     if (t < x) {
    1940                         b = x;
    1941                     } else {
    1942                         a = x;
    1943                     }
    1944
    1945                     // Assign the best approx to x
    1946                     v = w;
    1947                     w = x;
    1948                     x = t;
    1949
    1950                     fv = fw;
    1951                     fw = fx;
    1952                     fx = ft;
    1953                 // x remains the better approx
    1954                 } else {
    1955                     // Reduce the range enclosing x
    1956                     if (t < x) {
    1957                         a = t;
    1958                     } else {
    1959                         b = t;
    1960                     }
    1961
    1962                     if (ft <= fw || w === x) {
    1963                         v = w;
    1964                         w = t;
    1965                         fv = fw;
    1966                         fw = ft;
    1967                     } else if (ft <= fv || v === x || v === w) {
    1968                         v = t;
    1969                         fv = ft;
    1970                     }
    1971                 }
    1972                 niter += 1;
    1973             }
    1974
    1975             return x;
    1976         },
    1977
    1978         /**
    1979          * Implements the Ramer-Douglas-Peucker algorithm.
    1980          * It discards points which are not necessary from the polygonal line defined by the point array
    1981          * pts. The computation is done in screen coordinates.
    1982          * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points.
    1983          * @param {Array} pts Array of {@link JXG.Coords}
    1984          * @param {Number} eps If the absolute value of a given number <tt>x</tt> is smaller than <tt>eps</tt> it is considered to be equal <tt>0</tt>.
    1985          * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points.
    1986          * @memberof JXG.Math.Numerics
    1987          */
    1988         RamerDouglasPeucker: function (pts, eps) {
    1989             var newPts = [], i, k, len,
    1990
    1991                 /**
    1992                  * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeucker}.
    1993                  * It searches for the point between index i and j which
    1994                  * has the largest distance from the line between the points i and j.
    1995                  * @param {Array} pts Array of {@link JXG.Coords}
    1996                  * @param {Number} i Index of a point in pts
    1997                  * @param {Number} j Index of a point in pts
    1998                  * @ignore
    1999                  * @private
    2000                  */
    2001                 findSplit = function (pts, i, j) {
    2002                     var d, k, ci, cj, ck,
    2003                         x0, y0, x1, y1,
    2004                         den, lbda,
    2005                         dist = 0,
    2006                         f = i;
    2007
    2008                     if (j - i < 2) {
    2009                         return [-1.0, 0];
    2010                     }
    2011
    2012                     ci = pts[i].scrCoords;
    2013                     cj = pts[j].scrCoords;
    2014
    2015                     if (isNaN(ci[1] + ci[2])) {
    2016                         return [NaN, i];
    2017                     }
    2018                     if (isNaN(cj[1] + cj[2])) {
    2019                         return [NaN, j];
    2020                     }
    2021
    2022                     for (k = i + 1; k < j; k++) {
    2023                         ck = pts[k].scrCoords;
    2024                         if (isNaN(ck[1] + ck[2])) {
    2025                             return [NaN, k];
    2026                         }
    2027
    2028                         x0 = ck[1] - ci[1];
    2029                         y0 = ck[2] - ci[2];
    2030                         x1 = cj[1] - ci[1];
    2031                         y1 = cj[2] - ci[2];
    2032                         den = x1 * x1 + y1 * y1;
    2033
    2034                         if (den >= Mat.eps) {
    2035                             lbda = (x0 * x1 + y0 * y1) / den;
    2036
    2037                             if (lbda < 0.0) {
    2038                                 lbda = 0.0;
    2039                             } else if (lbda > 1.0) {
    2040                                 lbda = 1.0;
    2041                             }
    2042
    2043                             x0 = x0 - lbda * x1;
    2044                             y0 = y0 - lbda * y1;
    2045                             d = x0 * x0 + y0 * y0;
    2046                         } else {
    2047                             lbda = 0.0;
    2048                             d = x0 * x0 + y0 * y0;
    2049                         }
    2050
    2051                         if (d > dist) {
    2052                             dist = d;
    2053                             f = k;
    2054                         }
    2055                     }
    2056                     return [Math.sqrt(dist), f];
    2057                 },
    2058
    2059                 /**
    2060                  * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeucker}.
    2061                  * It runs recursively through the point set and searches the
    2062                  * point which has the largest distance from the line between the first point and
    2063                  * the last point. If the distance from the line is greater than eps, this point is
    2064                  * included in our new point set otherwise it is discarded.
    2065                  * If it is taken, we recursively apply the subroutine to the point set before
    2066                  * and after the chosen point.
    2067                  * @param {Array} pts Array of {@link JXG.Coords}
    2068                  * @param {Number} i Index of an element of pts
    2069                  * @param {Number} j Index of an element of pts
    2070                  * @param {Number} eps If the absolute value of a given number <tt>x</tt> is smaller than <tt>eps</tt> it is considered to be equal <tt>0</tt>.
    2071                  * @param {Array} newPts Array of {@link JXG.Coords}
    2072                  * @ignore
    2073                  * @private
    2074                  */
    2075                 RDP = function (pts, i, j, eps, newPts) {
    2076                     var result = findSplit(pts, i, j),
    2077                         k = result[1];
    2078
    2079                     if (isNaN(result[0])) {
    2080                         RDP(pts, i, k - 1, eps, newPts);
    2081                         newPts.push(pts[k]);
    2082                         do {
    2083                             ++k;
    2084                         } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2]));
    2085                         if (k <= j) {
    2086                             newPts.push(pts[k]);
    2087                         }
    2088                         RDP(pts, k + 1, j, eps, newPts);
    2089                     } else if (result[0] > eps) {
    2090                         RDP(pts, i, k, eps, newPts);
    2091                         RDP(pts, k, j, eps, newPts);
    2092                     } else {
    2093                         newPts.push(pts[j]);
    2094                     }
    2095                 };
    2096
    2097             len = pts.length;
    2098
    2099             // Search for the left most point woithout NaN coordinates
    2100             i = 0;
    2101             while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) {
    2102                 i += 1;
    2103             }
    2104             // Search for the right most point woithout NaN coordinates
    2105             k = len - 1;
    2106             while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) {
    2107                 k -= 1;
    2108             }
    2109
    2110             // Only proceed if something is left
    2111             if (!(i > k || i === len)) {
    2112                 newPts[0] = pts[i];
    2113                 RDP(pts, i, k, eps, newPts);
    2114             }
    2115
    2116             return newPts;
    2117         },
    2118
    2119         /**
    2120          * Old name for the implementation of the Ramer-Douglas-Peucker algorithm.
    2121          * @deprecated Use {@link JXG.Math.Numerics#RamerDouglasPeucker}
    2122          * @memberof JXG.Math.Numerics
    2123          */
    2124         RamerDouglasPeuker: function (pts, eps) {
    2125             return this.RamerDouglasPeucker(pts, eps);
    2126         }
    2127     };
    2128
    2129     return Mat.Numerics;
    2130 });
    2131