      1 /*
      2     Copyright 2008-2022
      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 /*global JXG: true, define: true*/
     33 /*jslint nomen: true, plusplus: true*/
     34 /*eslint no-loss-of-precision: off */
     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(['jxg', 'utils/type', 'utils/env', 'math/math'], function (JXG, Type, Env, 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      * @name JXG.Math.Numerics
     85      * @exports Mat.Numerics as JXG.Math.Numerics
     86      * @namespace
     87      */
     88     Mat.Numerics = {
     89
     90     //JXG.extend(Mat.Numerics, /** @lends JXG.Math.Numerics */ {
     91         /**
     92          * Solves a system of linear equations given by A and b using the Gauss-Jordan-elimination.
     93          * The algorithm runs in-place. I.e. the entries of A and b are changed.
     94          * @param {Array} A Square matrix represented by an array of rows, containing the coefficients of the lineare equation system.
     95          * @param {Array} b A vector containing the linear equation system's right hand side.
     96          * @throws {Error} If a non-square-matrix is given or if b has not the right length or A's rank is not full.
     97          * @returns {Array} A vector that solves the linear equation system.
     98          * @memberof JXG.Math.Numerics
     99          */
    100         Gauss: function (A, b) {
    101             var i, j, k,
    102                 // copy the matrix to prevent changes in the original
    103                 Acopy,
    104                 // solution vector, to prevent changing b
    105                 x,
    106                 eps = Mat.eps,
    107                 // number of columns of A
    108                 n = A.length > 0 ? A[0].length : 0;
    109
    110             if ((n !== b.length) || (n !== A.length)) {
    111                 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.");
    112             }
    113
    114             // initialize solution vector
    115             Acopy = [];
    116             x = b.slice(0, n);
    117
    118             for (i = 0; i < n; i++) {
    119                 Acopy[i] = A[i].slice(0, n);
    120             }
    121
    122             // Gauss-Jordan-elimination
    123             for (j = 0; j < n; j++) {
    124                 for (i = n - 1; i > j; i--) {
    125                     // Is the element which is to eliminate greater than zero?
    126                     if (Math.abs(Acopy[i][j]) > eps) {
    127                         // Equals pivot element zero?
    128                         if (Math.abs(Acopy[j][j]) < eps) {
    129                             // At least numerically, so we have to exchange the rows
    130                             Type.swap(Acopy, i, j);
    131                             Type.swap(x, i, j);
    132                         } else {
    133                             // Saves the L matrix of the LR-decomposition. unnecessary.
    134                             Acopy[i][j] /= Acopy[j][j];
    135                             // Transform right-hand-side b
    136                             x[i] -= Acopy[i][j] * x[j];
    137
    138                             // subtract the multiple of A[i][j] / A[j][j] of the j-th row from the i-th.
    139                             for (k = j + 1; k < n; k++) {
    140                                 Acopy[i][k] -= Acopy[i][j] * Acopy[j][k];
    141                             }
    142                         }
    143                     }
    144                 }
    145
    146                 // The absolute values of all coefficients below the j-th row in the j-th column are smaller than JXG.Math.eps.
    147                 if (Math.abs(Acopy[j][j]) < eps) {
    148                     throw new Error("JXG.Math.Numerics.Gauss(): The given matrix seems to be singular.");
    149                 }
    150             }
    151
    152             this.backwardSolve(Acopy, x, true);
    153
    154             return x;
    155         },
    156
    157         /**
    158          * Solves a system of linear equations given by the right triangular matrix R and vector b.
    159          * @param {Array} R Right triangular matrix represented by an array of rows. All entries a_(i,j) with i < j are ignored.
    160          * @param {Array} b Right hand side of the linear equation system.
    161          * @param {Boolean} [canModify=false] If true, the right hand side vector is allowed to be changed by this method.
    162          * @returns {Array} An array representing a vector that solves the system of linear equations.
    163          * @memberof JXG.Math.Numerics
    164          */
    165         backwardSolve: function (R, b, canModify) {
    166             var x, m, n, i, j;
    167
    168             if (canModify) {
    169                 x = b;
    170             } else {
    171                 x = b.slice(0, b.length);
    172             }
    173
    174             // m: number of rows of R
    175             // n: number of columns of R
    176             m = R.length;
    177             n = R.length > 0 ? R[0].length : 0;
    178
    179             for (i = m - 1; i >= 0; i--) {
    180                 for (j = n - 1; j > i; j--) {
    181                     x[i] -= R[i][j] * x[j];
    182                 }
    183                 x[i] /= R[i][i];
    184             }
    185
    186             return x;
    187         },
    188
    189         /**
    190          * @private
    191          * Gauss-Bareiss algorithm to compute the
    192          * determinant of matrix without fractions.
    193          * See Henri Cohen, "A Course in Computational
    194          * Algebraic Number Theory (Graduate texts
    195          * in mathematics; 138)", Springer-Verlag,
    196          * ISBN 3-540-55640-0 / 0-387-55640-0
    197          * Third, Corrected Printing 1996
    198          * "Algorithm 2.2.6", pg. 52-53
    199          * @memberof JXG.Math.Numerics
    200          */
    201         gaussBareiss: function (mat) {
    202             var k, c, s, i, j, p, n, M, t,
    203                 eps = Mat.eps;
    204
    205             n = mat.length;
    206
    207             if (n <= 0) {
    208                 return 0;
    209             }
    210
    211             if (mat[0].length < n) {
    212                 n = mat[0].length;
    213             }
    214
    215             // Copy the input matrix to M
    216             M = [];
    217
    218             for (i = 0; i < n; i++) {
    219                 M[i] = mat[i].slice(0, n);
    220             }
    221
    222             c = 1;
    223             s = 1;
    224
    225             for (k = 0; k < n - 1; k++) {
    226                 p = M[k][k];
    227
    228                 // Pivot step
    229                 if (Math.abs(p) < eps) {
    230                     for (i = k + 1; i < n; i++) {
    231                         if (Math.abs(M[i][k]) >= eps) {
    232                             break;
    233                         }
    234                     }
    235
    236                     // No nonzero entry found in column k -> det(M) = 0
    237                     if (i === n) {
    238                         return 0.0;
    239                     }
    240
    241                     // swap row i and k partially
    242                     for (j = k; j < n; j++) {
    243                         t = M[i][j];
    244                         M[i][j] = M[k][j];
    245                         M[k][j] = t;
    246                     }
    247                     s = -s;
    248                     p = M[k][k];
    249                 }
    250
    251                 // Main step
    252                 for (i = k + 1; i < n; i++) {
    253                     for (j = k + 1; j < n; j++) {
    254                         t = p * M[i][j] - M[i][k] * M[k][j];
    255                         M[i][j] = t / c;
    256                     }
    257                 }
    258
    259                 c = p;
    260             }
    261
    262             return s * M[n - 1][n - 1];
    263         },
    264
    265         /**
    266          * Computes the determinant of a square nxn matrix with the
    267          * Gauss-Bareiss algorithm.
    268          * @param {Array} mat Matrix.
    269          * @returns {Number} The determinant pf the matrix mat.
    270          *                   The empty matrix returns 0.
    271          * @memberof JXG.Math.Numerics
    272          */
    273         det: function (mat) {
    274             var n = mat.length;
    275
    276             if (n === 2 && mat[0].length === 2) {
    277                 return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
    278             }
    279
    280             return this.gaussBareiss(mat);
    281         },
    282
    283         /**
    284          * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method
    285          * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990
    286          * @param {Array} Ain A symmetric 3x3 matrix.
    287          * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors.
    288          * @memberof JXG.Math.Numerics
    289          */
    290         Jacobi: function (Ain) {
    291             var i, j, k, aa, si, co, tt, ssum, amax,
    292                 eps = Mat.eps * Mat.eps,
    293                 sum = 0.0,
    294                 n = Ain.length,
    295                 V = [
    296                     [0, 0, 0],
    297                     [0, 0, 0],
    298                     [0, 0, 0]
    299                 ],
    300                 A = [
    301                     [0, 0, 0],
    302                     [0, 0, 0],
    303                     [0, 0, 0]
    304                 ],
    305                 nloops = 0;
    306
    307             // Initialization. Set initial Eigenvectors.
    308             for (i = 0; i < n; i++) {
    309                 for (j = 0; j < n; j++) {
    310                     V[i][j] = 0.0;
    311                     A[i][j] = Ain[i][j];
    312                     sum += Math.abs(A[i][j]);
    313                 }
    314                 V[i][i] = 1.0;
    315             }
    316
    317             // Trivial problems
    318             if (n === 1) {
    319                 return [A, V];
    320             }
    321
    322             if (sum <= 0.0) {
    323                 return [A, V];
    324             }
    325
    326             sum /= (n * n);
    327
    328             // Reduce matrix to diagonal
    329             do {
    330                 ssum = 0.0;
    331                 amax = 0.0;
    332                 for (j = 1; j < n; j++) {
    333                     for (i = 0; i < j; i++) {
    334                         // Check if A[i][j] is to be reduced
    335                         aa = Math.abs(A[i][j]);
    336
    337                         if (aa > amax) {
    338                             amax = aa;
    339                         }
    340
    341                         ssum += aa;
    342
    343                         if (aa >= eps) {
    344                             // calculate rotation angle
    345                             aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5;
    346                             si = Math.sin(aa);
    347                             co = Math.cos(aa);
    348
    349                             // Modify 'i' and 'j' columns
    350                             for (k = 0; k < n; k++) {
    351                                 tt = A[k][i];
    352                                 A[k][i] = co * tt + si * A[k][j];
    353                                 A[k][j] = -si * tt + co * A[k][j];
    354                                 tt = V[k][i];
    355                                 V[k][i] = co * tt + si * V[k][j];
    356                                 V[k][j] = -si * tt + co * V[k][j];
    357                             }
    358
    359                             // Modify diagonal terms
    360                             A[i][i] = co * A[i][i] + si * A[j][i];
    361                             A[j][j] = -si * A[i][j] + co * A[j][j];
    362                             A[i][j] = 0.0;
    363
    364                             // Make 'A' matrix symmetrical
    365                             for (k = 0; k < n; k++) {
    366                                 A[i][k] = A[k][i];
    367                                 A[j][k] = A[k][j];
    368                             }
    369                             // A[i][j] made zero by rotation
    370                         }
    371                     }
    372                 }
    373                 nloops += 1;
    374             } while (Math.abs(ssum) / sum > eps && nloops < 2000);
    375
    376             return [A, V];
    377         },
    378
    379         /**
    380          * Calculates the integral of function f over interval using Newton-Cotes-algorithm.
    381          * @param {Array} interval The integration interval, e.g. [0, 3].
    382          * @param {function} f A function which takes one argument of type number and returns a number.
    383          * @param {Object} [config] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type
    384          * with value being either 'trapez', 'simpson', or 'milne'.
    385          * @param {Number} [config.number_of_nodes=28]
    386          * @param {String} [config.integration_type='milne'] Possible values are 'milne', 'simpson', 'trapez'
    387          * @returns {Number} Integral value of f over interval
    388          * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use
    389          * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4.
    390          * @example
    391          * function f(x) {
    392          *   return x*x;
    393          * }
    394          *
    395          * // calculates integral of <tt>f</tt> from 0 to 2.
    396          * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f);
    397          *
    398          * // the same with an anonymous function
    399          * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; });
    400          *
    401          * // use trapez rule with 16 nodes
    402          * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f,
    403          *                                   {number_of_nodes: 16, integration_type: 'trapez'});
    404          * @memberof JXG.Math.Numerics
    405          */
    406         NewtonCotes: function (interval, f, config) {
    407             var evaluation_point, i, number_of_intervals,
    408                 integral_value = 0.0,
    409                 number_of_nodes = config && Type.isNumber(config.number_of_nodes) ? config.number_of_nodes : 28,
    410                 available_types = {trapez: true, simpson: true, milne: true},
    411                 integration_type = config && config.integration_type && available_types.hasOwnProperty(config.integration_type) && available_types[config.integration_type] ? config.integration_type : 'milne',
    412                 step_size = (interval[1] - interval[0]) / number_of_nodes;
    413
    414             switch (integration_type) {
    415             case 'trapez':
    416                 integral_value = (f(interval[0]) + f(interval[1])) * 0.5;
    417                 evaluation_point = interval[0];
    418
    419                 for (i = 0; i < number_of_nodes - 1; i++) {
    420                     evaluation_point += step_size;
    421                     integral_value += f(evaluation_point);
    422                 }
    423
    424                 integral_value *= step_size;
    425                 break;
    426             case 'simpson':
    427                 if (number_of_nodes % 2 > 0) {
    428                     throw new Error("JSXGraph:  INT_SIMPSON requires config.number_of_nodes dividable by 2.");
    429                 }
    430
    431                 number_of_intervals = number_of_nodes / 2.0;
    432                 integral_value = f(interval[0]) + f(interval[1]);
    433                 evaluation_point = interval[0];
    434
    435                 for (i = 0; i < number_of_intervals - 1; i++) {
    436                     evaluation_point += 2.0 * step_size;
    437                     integral_value += 2.0 * f(evaluation_point);
    438                 }
    439
    440                 evaluation_point = interval[0] - step_size;
    441
    442                 for (i = 0; i < number_of_intervals; i++) {
    443                     evaluation_point += 2.0 * step_size;
    444                     integral_value += 4.0 * f(evaluation_point);
    445                 }
    446
    447                 integral_value *= step_size / 3.0;
    448                 break;
    449             default:
    450                 if (number_of_nodes % 4 > 0) {
    451                     throw new Error("JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4");
    452                 }
    453
    454                 number_of_intervals = number_of_nodes * 0.25;
    455                 integral_value = 7.0 * (f(interval[0]) + f(interval[1]));
    456                 evaluation_point = interval[0];
    457
    458                 for (i = 0; i < number_of_intervals - 1; i++) {
    459                     evaluation_point += 4.0 * step_size;
    460                     integral_value += 14.0 * f(evaluation_point);
    461                 }
    462
    463                 evaluation_point = interval[0] - 3.0 * step_size;
    464
    465                 for (i = 0; i < number_of_intervals; i++) {
    466                     evaluation_point += 4.0 * step_size;
    467                     integral_value += 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size));
    468                 }
    469
    470                 evaluation_point = interval[0] - 2.0 * step_size;
    471
    472                 for (i = 0; i < number_of_intervals; i++) {
    473                     evaluation_point += 4.0 * step_size;
    474                     integral_value += 12.0 * f(evaluation_point);
    475                 }
    476
    477                 integral_value *= 2.0 * step_size / 45.0;
    478             }
    479             return integral_value;
    480         },
    481
    482        /**
    483          * Calculates the integral of function f over interval using Romberg iteration.
    484          * @param {Array} interval The integration interval, e.g. [0, 3].
    485          * @param {function} f A function which takes one argument of type number and returns a number.
    486          * @param {Object} [config] The algorithm setup. Accepted properties are max_iterations of type number and precision eps.
    487          * @param {Number} [config.max_iterations=20]
    488          * @param {Number} [config.eps=0.0000001]
    489          * @returns {Number} Integral value of f over interval
    490          * @example
    491          * function f(x) {
    492          *   return x*x;
    493          * }
    494          *
    495          * // calculates integral of <tt>f</tt> from 0 to 2.
    496          * var area1 = JXG.Math.Numerics.Romberg([0, 2], f);
    497          *
    498          * // the same with an anonymous function
    499          * var area2 = JXG.Math.Numerics.Romberg([0, 2], function (x) { return x*x; });
    500          *
    501          * // use trapez rule with maximum of 16 iterations or stop if the precision 0.0001 has been reached.
    502          * var area3 = JXG.Math.Numerics.Romberg([0, 2], f,
    503          *                                   {max_iterations: 16, eps: 0.0001});
    504          * @memberof JXG.Math.Numerics
    505          */
    506         Romberg: function (interval, f, config) {
    507             var a, b, h, s, n,
    508                 k, i, q,
    509                 p = [],
    510                 integral = 0.0,
    511                 last = Infinity,
    512                 m = config && Type.isNumber(config.max_iterations) ? config.max_iterations : 20,
    513                 eps = config && Type.isNumber(config.eps) ? config.eps : config.eps || 0.0000001;
    514
    515             a = interval[0];
    516             b = interval[1];
    517             h = b - a;
    518             n = 1;
    519
    520             p[0] = 0.5 * h * (f(a) + f(b));
    521
    522             for (k = 0; k < m; ++k) {
    523                 s = 0;
    524                 h *= 0.5;
    525                 n *= 2;
    526                 q = 1;
    527
    528                 for (i = 1; i < n; i += 2) {
    529                     s += f(a + i * h);
    530                 }
    531
    532                 p[k + 1] = 0.5 * p[k] + s * h;
    533
    534                 integral = p[k + 1];
    535                 for (i = k - 1; i >= 0; --i) {
    536                     q *= 4;
    537                     p[i] = p[i + 1] + (p[i + 1] - p[i]) / (q - 1.0);
    538                     integral = p[i];
    539                 }
    540
    541                 if (Math.abs(integral - last) < eps * Math.abs(integral)) {
    542                     break;
    543                 }
    544                 last = integral;
    545             }
    546
    547             return integral;
    548         },
    549
    550        /**
    551          * Calculates the integral of function f over interval using Gauss-Legendre quadrature.
    552          * @param {Array} interval The integration interval, e.g. [0, 3].
    553          * @param {function} f A function which takes one argument of type number and returns a number.
    554          * @param {Object} [config] The algorithm setup. Accepted property is the order n of type number. n is allowed to take
    555          * values between 2 and 18, default value is 12.
    556          * @param {Number} [config.n=16]
    557          * @returns {Number} Integral value of f over interval
    558          * @example
    559          * function f(x) {
    560          *   return x*x;
    561          * }
    562          *
    563          * // calculates integral of <tt>f</tt> from 0 to 2.
    564          * var area1 = JXG.Math.Numerics.GaussLegendre([0, 2], f);
    565          *
    566          * // the same with an anonymous function
    567          * var area2 = JXG.Math.Numerics.GaussLegendre([0, 2], function (x) { return x*x; });
    568          *
    569          * // use 16 point Gauss-Legendre rule.
    570          * var area3 = JXG.Math.Numerics.GaussLegendre([0, 2], f,
    571          *                                   {n: 16});
    572          * @memberof JXG.Math.Numerics
    573          */
    574         GaussLegendre: function (interval, f, config) {
    575             var a, b,
    576                 i, m,
    577                 xp, xm,
    578                 result = 0.0,
    579                 table_xi = [],
    580                 table_w = [],
    581                 xi, w,
    582                 n = config && Type.isNumber(config.n) ? config.n : 12;
    583
    584             if (n > 18) {
    585                 n = 18;
    586             }
    587
    588             /* n = 2 */
    589             table_xi[2] = [0.5773502691896257645091488];
    590             table_w[2] = [1.0000000000000000000000000];
    591
    592             /* n = 4 */
    593             table_xi[4] = [0.3399810435848562648026658, 0.8611363115940525752239465];
    594             table_w[4] = [0.6521451548625461426269361, 0.3478548451374538573730639];
    595
    596             /* n = 6 */
    597             table_xi[6] = [0.2386191860831969086305017, 0.6612093864662645136613996, 0.9324695142031520278123016];
    598             table_w[6] = [0.4679139345726910473898703, 0.3607615730481386075698335, 0.1713244923791703450402961];
    599
    600             /* n = 8 */
    601             table_xi[8] = [0.1834346424956498049394761, 0.5255324099163289858177390, 0.7966664774136267395915539, 0.9602898564975362316835609];
    602             table_w[8] = [0.3626837833783619829651504, 0.3137066458778872873379622, 0.2223810344533744705443560, 0.1012285362903762591525314];
    603
    604             /* n = 10 */
    605             table_xi[10] = [0.1488743389816312108848260, 0.4333953941292471907992659, 0.6794095682990244062343274, 0.8650633666889845107320967, 0.9739065285171717200779640];
    606             table_w[10] = [0.2955242247147528701738930, 0.2692667193099963550912269, 0.2190863625159820439955349, 0.1494513491505805931457763, 0.0666713443086881375935688];
    607
    608             /* n = 12 */
    609             table_xi[12] = [0.1252334085114689154724414, 0.3678314989981801937526915, 0.5873179542866174472967024, 0.7699026741943046870368938, 0.9041172563704748566784659, 0.9815606342467192506905491];
    610             table_w[12] = [0.2491470458134027850005624, 0.2334925365383548087608499, 0.2031674267230659217490645, 0.1600783285433462263346525, 0.1069393259953184309602547, 0.0471753363865118271946160];
    611
    612             /* n = 14 */
    613             table_xi[14] = [0.1080549487073436620662447, 0.3191123689278897604356718, 0.5152486363581540919652907, 0.6872929048116854701480198, 0.8272013150697649931897947, 0.9284348836635735173363911, 0.9862838086968123388415973];
    614             table_w[14] = [0.2152638534631577901958764, 0.2051984637212956039659241, 0.1855383974779378137417166, 0.1572031671581935345696019, 0.1215185706879031846894148, 0.0801580871597602098056333, 0.0351194603317518630318329];
    615
    616             /* n = 16 */
    617             table_xi[16] = [0.0950125098376374401853193, 0.2816035507792589132304605, 0.4580167776572273863424194, 0.6178762444026437484466718, 0.7554044083550030338951012, 0.8656312023878317438804679, 0.9445750230732325760779884, 0.9894009349916499325961542];
    618             table_w[16] = [0.1894506104550684962853967, 0.1826034150449235888667637, 0.1691565193950025381893121, 0.1495959888165767320815017, 0.1246289712555338720524763, 0.0951585116824927848099251, 0.0622535239386478928628438, 0.0271524594117540948517806];
    619
    620             /* n = 18 */
    621             table_xi[18] = [0.0847750130417353012422619, 0.2518862256915055095889729, 0.4117511614628426460359318, 0.5597708310739475346078715, 0.6916870430603532078748911, 0.8037049589725231156824175, 0.8926024664975557392060606, 0.9558239495713977551811959, 0.9915651684209309467300160];
    622             table_w[18] = [0.1691423829631435918406565, 0.1642764837458327229860538, 0.1546846751262652449254180, 0.1406429146706506512047313, 0.1225552067114784601845191, 0.1009420441062871655628140, 0.0764257302548890565291297, 0.0497145488949697964533349, 0.0216160135264833103133427];
    623
    624             /* n = 3 */
    625             table_xi[3] = [0.0000000000000000000000000, 0.7745966692414833770358531];
    626             table_w[3] = [0.8888888888888888888888889, 0.5555555555555555555555556];
    627
    628             /* n = 5 */
    629             table_xi[5] = [0.0000000000000000000000000, 0.5384693101056830910363144, 0.9061798459386639927976269];
    630             table_w[5] = [0.5688888888888888888888889, 0.4786286704993664680412915, 0.2369268850561890875142640];
    631
    632             /* n = 7 */
    633             table_xi[7] = [0.0000000000000000000000000, 0.4058451513773971669066064, 0.7415311855993944398638648, 0.9491079123427585245261897];
    634             table_w[7] = [0.4179591836734693877551020, 0.3818300505051189449503698, 0.2797053914892766679014678, 0.1294849661688696932706114];
    635
    636             /* n = 9 */
    637             table_xi[9] = [0.0000000000000000000000000, 0.3242534234038089290385380, 0.6133714327005903973087020, 0.8360311073266357942994298, 0.9681602395076260898355762];
    638             table_w[9] = [0.3302393550012597631645251, 0.3123470770400028400686304, 0.2606106964029354623187429, 0.1806481606948574040584720, 0.0812743883615744119718922];
    639
    640             /* n = 11 */
    641             table_xi[11] = [0.0000000000000000000000000, 0.2695431559523449723315320, 0.5190961292068118159257257, 0.7301520055740493240934163, 0.8870625997680952990751578, 0.9782286581460569928039380];
    642             table_w[11] = [0.2729250867779006307144835, 0.2628045445102466621806889, 0.2331937645919904799185237, 0.1862902109277342514260976, 0.1255803694649046246346943, 0.0556685671161736664827537];
    643
    644             /* n = 13 */
    645             table_xi[13] = [0.0000000000000000000000000, 0.2304583159551347940655281, 0.4484927510364468528779129, 0.6423493394403402206439846, 0.8015780907333099127942065, 0.9175983992229779652065478, 0.9841830547185881494728294];
    646             table_w[13] = [0.2325515532308739101945895, 0.2262831802628972384120902, 0.2078160475368885023125232, 0.1781459807619457382800467, 0.1388735102197872384636018, 0.0921214998377284479144218, 0.0404840047653158795200216];
    647
    648             /* n = 15 */
    649             table_xi[15] = [0.0000000000000000000000000, 0.2011940939974345223006283, 0.3941513470775633698972074, 0.5709721726085388475372267, 0.7244177313601700474161861, 0.8482065834104272162006483, 0.9372733924007059043077589, 0.9879925180204854284895657];
    650             table_w[15] = [0.2025782419255612728806202, 0.1984314853271115764561183, 0.1861610000155622110268006, 0.1662692058169939335532009, 0.1395706779261543144478048, 0.1071592204671719350118695, 0.0703660474881081247092674, 0.0307532419961172683546284];
    651
    652             /* n = 17 */
    653             table_xi[17] = [0.0000000000000000000000000, 0.1784841814958478558506775, 0.3512317634538763152971855, 0.5126905370864769678862466, 0.6576711592166907658503022, 0.7815140038968014069252301, 0.8802391537269859021229557, 0.9506755217687677612227170, 0.9905754753144173356754340];
    654             table_w[17] = [0.1794464703562065254582656, 0.1765627053669926463252710, 0.1680041021564500445099707, 0.1540457610768102880814316, 0.1351363684685254732863200, 0.1118838471934039710947884, 0.0850361483171791808835354, 0.0554595293739872011294402, 0.0241483028685479319601100];
    655
    656             a = interval[0];
    657             b = interval[1];
    658
    659             //m = Math.ceil(n * 0.5);
    660             m = (n + 1) >> 1;
    661
    662             xi = table_xi[n];
    663             w = table_w[n];
    664
    665             xm = 0.5 * (b - a);
    666             xp = 0.5 * (b + a);
    667
    668             if (n & 1 === 1) { // n odd
    669                 result = w[0] * f(xp);
    670                 for (i = 1; i < m; ++i) {
    671                     result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i]));
    672                 }
    673             } else { // n even
    674                 result = 0.0;
    675                 for (i = 0; i < m; ++i) {
    676                     result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i]));
    677                 }
    678             }
    679
    680             return xm * result;
    681         },
    682
    683         /**
    684          * Scale error in Gauss Kronrod quadrature.
    685          * Internal method used in {@link JXG.Math.Numerics._gaussKronrod}.
    686          * @private
    687          */
    688         _rescale_error: function (err, result_abs, result_asc) {
    689             var scale, min_err,
    690                 DBL_MIN = 2.2250738585072014e-308,
    691                 DBL_EPS = 2.2204460492503131e-16;
    692
    693             err = Math.abs(err);
    694             if (result_asc !== 0 && err !== 0) {
    695                 scale = Math.pow((200 * err / result_asc), 1.5);
    696
    697                 if (scale < 1.0) {
    698                     err = result_asc * scale;
    699                 } else {
    700                     err = result_asc;
    701                 }
    702             }
    703             if (result_abs > DBL_MIN / (50 * DBL_EPS)) {
    704                 min_err = 50 * DBL_EPS * result_abs;
    705
    706                 if (min_err > err) {
    707                     err = min_err;
    708                 }
    709             }
    710
    711             return err;
    712         },
    713
    714         /**
    715          * Generic Gauss-Kronrod quadrature algorithm.
    716          * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15},
    717          * {@link JXG.Math.Numerics.GaussKronrod21},
    718          * {@link JXG.Math.Numerics.GaussKronrod31}.
    719          * Taken from QUADPACK.
    720          *
    721          * @param {Array} interval The integration interval, e.g. [0, 3].
    722          * @param {function} f A function which takes one argument of type number and returns a number.
    723          * @param {Number} n order
    724          * @param {Array} xgk Kronrod quadrature abscissae
    725          * @param {Array} wg Weights of the Gauss rule
    726          * @param {Array} wgk Weights of the Kronrod rule
    727          * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc.
    728          * See the library QUADPACK for an explanation.
    729          *
    730          * @returns {Number} Integral value of f over interval
    731          *
    732          * @private
    733          */
    734         _gaussKronrod: function (interval, f, n, xgk, wg, wgk, resultObj) {
    735             var a = interval[0],
    736                 b = interval[1],
    737                 up,
    738                 result,
    739
    740                 center = 0.5 * (a + b),
    741                 half_length = 0.5 * (b - a),
    742                 abs_half_length = Math.abs(half_length),
    743                 f_center = f(center),
    744
    745                 result_gauss = 0.0,
    746                 result_kronrod = f_center * wgk[n - 1],
    747
    748                 result_abs = Math.abs(result_kronrod),
    749                 result_asc = 0.0,
    750                 mean = 0.0,
    751                 err = 0.0,
    752
    753                 j, jtw, abscissa, fval1, fval2, fsum,
    754                 jtwm1,
    755                 fv1 = [], fv2 = [];
    756
    757             if (n % 2 === 0) {
    758                 result_gauss = f_center * wg[n / 2 - 1];
    759             }
    760
    761             up = Math.floor((n - 1) / 2);
    762             for (j = 0; j < up; j++) {
    763                 jtw = j * 2 + 1;  // in original fortran j=1,2,3 jtw=2,4,6
    764                 abscissa = half_length * xgk[jtw];
    765                 fval1 = f(center - abscissa);
    766                 fval2 = f(center + abscissa);
    767                 fsum = fval1 + fval2;
    768                 fv1[jtw] = fval1;
    769                 fv2[jtw] = fval2;
    770                 result_gauss += wg[j] * fsum;
    771                 result_kronrod += wgk[jtw] * fsum;
    772                 result_abs += wgk[jtw] * (Math.abs(fval1) + Math.abs(fval2));
    773             }
    774
    775             up = Math.floor(n / 2);
    776             for (j = 0; j < up; j++) {
    777                 jtwm1 = j * 2;
    778                 abscissa = half_length * xgk[jtwm1];
    779                 fval1 = f(center - abscissa);
    780                 fval2 = f(center + abscissa);
    781                 fv1[jtwm1] = fval1;
    782                 fv2[jtwm1] = fval2;
    783                 result_kronrod += wgk[jtwm1] * (fval1 + fval2);
    784                 result_abs += wgk[jtwm1] * (Math.abs(fval1) + Math.abs(fval2));
    785             }
    786
    787             mean = result_kronrod * 0.5;
    788             result_asc = wgk[n - 1] * Math.abs(f_center - mean);
    789
    790             for (j = 0; j < n - 1; j++) {
    791                 result_asc += wgk[j] * (Math.abs(fv1[j] - mean) + Math.abs(fv2[j] - mean));
    792             }
    793
    794             // scale by the width of the integration region
    795             err = (result_kronrod - result_gauss) * half_length;
    796
    797             result_kronrod *= half_length;
    798             result_abs *= abs_half_length;
    799             result_asc *= abs_half_length;
    800             result = result_kronrod;
    801
    802             resultObj.abserr = this._rescale_error(err, result_abs, result_asc);
    803             resultObj.resabs = result_abs;
    804             resultObj.resasc = result_asc;
    805
    806             return result;
    807         },
    808
    809         /**
    810          * 15 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK
    811          * @param {Array} interval The integration interval, e.g. [0, 3].
    812          * @param {function} f A function which takes one argument of type number and returns a number.
    813          * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library
    814          *  QUADPACK for an explanation.
    815          *
    816          * @returns {Number} Integral value of f over interval
    817          *
    818          * @memberof JXG.Math.Numerics
    819          */
    820         GaussKronrod15: function (interval, f, resultObj) {
    821             /* Gauss quadrature weights and kronrod quadrature abscissae and
    822                 weights as evaluated with 80 decimal digit arithmetic by
    823                 L. W. Fullerton, Bell Labs, Nov. 1981. */
    824
    825             var xgk =    /* abscissae of the 15-point kronrod rule */
    826                     [
    827                         0.991455371120812639206854697526329,
    828                         0.949107912342758524526189684047851,
    829                         0.864864423359769072789712788640926,
    830                         0.741531185599394439863864773280788,
    831                         0.586087235467691130294144838258730,
    832                         0.405845151377397166906606412076961,
    833                         0.207784955007898467600689403773245,
    834                         0.000000000000000000000000000000000
    835                     ],
    836
    837             /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
    838                 xgk[0], xgk[2], ... abscissae to optimally extend the 7-point gauss rule */
    839
    840                 wg =     /* weights of the 7-point gauss rule */
    841                     [
    842                         0.129484966168869693270611432679082,
    843                         0.279705391489276667901467771423780,
    844                         0.381830050505118944950369775488975,
    845                         0.417959183673469387755102040816327
    846                     ],
    847
    848                 wgk =    /* weights of the 15-point kronrod rule */
    849                     [
    850                         0.022935322010529224963732008058970,
    851                         0.063092092629978553290700663189204,
    852                         0.104790010322250183839876322541518,
    853                         0.140653259715525918745189590510238,
    854                         0.169004726639267902826583426598550,
    855                         0.190350578064785409913256402421014,
    856                         0.204432940075298892414161999234649,
    857                         0.209482141084727828012999174891714
    858                     ];
    859
    860             return this._gaussKronrod(interval, f, 8, xgk, wg, wgk, resultObj);
    861         },
    862
    863         /**
    864          * 21 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK
    865          * @param {Array} interval The integration interval, e.g. [0, 3].
    866          * @param {function} f A function which takes one argument of type number and returns a number.
    867          * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library
    868          *  QUADPACK for an explanation.
    869          *
    870          * @returns {Number} Integral value of f over interval
    871          *
    872          * @memberof JXG.Math.Numerics
    873          */
    874         GaussKronrod21: function (interval, f, resultObj) {
    875             /* Gauss quadrature weights and kronrod quadrature abscissae and
    876                 weights as evaluated with 80 decimal digit arithmetic by
    877                 L. W. Fullerton, Bell Labs, Nov. 1981. */
    878
    879             var xgk =   /* abscissae of the 21-point kronrod rule */
    880                     [
    881                         0.995657163025808080735527280689003,
    882                         0.973906528517171720077964012084452,
    883                         0.930157491355708226001207180059508,
    884                         0.865063366688984510732096688423493,
    885                         0.780817726586416897063717578345042,
    886                         0.679409568299024406234327365114874,
    887                         0.562757134668604683339000099272694,
    888                         0.433395394129247190799265943165784,
    889                         0.294392862701460198131126603103866,
    890                         0.148874338981631210884826001129720,
    891                         0.000000000000000000000000000000000
    892                     ],
    893
    894                 /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule.
    895                 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */
    896                 wg =     /* weights of the 10-point gauss rule */
    897                     [
    898                         0.066671344308688137593568809893332,
    899                         0.149451349150580593145776339657697,
    900                         0.219086362515982043995534934228163,
    901                         0.269266719309996355091226921569469,
    902                         0.295524224714752870173892994651338
    903                     ],
    904
    905                 wgk =   /* weights of the 21-point kronrod rule */
    906                     [
    907                         0.011694638867371874278064396062192,
    908                         0.032558162307964727478818972459390,
    909                         0.054755896574351996031381300244580,
    910                         0.075039674810919952767043140916190,
    911                         0.093125454583697605535065465083366,
    912                         0.109387158802297641899210590325805,
    913                         0.123491976262065851077958109831074,
    914                         0.134709217311473325928054001771707,
    915                         0.142775938577060080797094273138717,
    916                         0.147739104901338491374841515972068,
    917                         0.149445554002916905664936468389821
    918                     ];
    919
    920             return this._gaussKronrod(interval, f, 11, xgk, wg, wgk, resultObj);
    921         },
    922
    923         /**
    924          * 31 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK
    925          * @param {Array} interval The integration interval, e.g. [0, 3].
    926          * @param {function} f A function which takes one argument of type number and returns a number.
    927          * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library
    928          *  QUADPACK for an explanation.
    929          *
    930          * @returns {Number} Integral value of f over interval
    931          *
    932          * @memberof JXG.Math.Numerics
    933          */
    934         GaussKronrod31: function (interval, f, resultObj) {
    935             /* Gauss quadrature weights and kronrod quadrature abscissae and
    936                 weights as evaluated with 80 decimal digit arithmetic by
    937                 L. W. Fullerton, Bell Labs, Nov. 1981. */
    938
    939             var xgk =   /* abscissae of the 21-point kronrod rule */
    940                     [
    941                         0.998002298693397060285172840152271,
    942                         0.987992518020485428489565718586613,
    943                         0.967739075679139134257347978784337,
    944                         0.937273392400705904307758947710209,
    945                         0.897264532344081900882509656454496,
    946                         0.848206583410427216200648320774217,
    947                         0.790418501442465932967649294817947,
    948                         0.724417731360170047416186054613938,
    949                         0.650996741297416970533735895313275,
    950                         0.570972172608538847537226737253911,
    951                         0.485081863640239680693655740232351,
    952                         0.394151347077563369897207370981045,
    953                         0.299180007153168812166780024266389,
    954                         0.201194093997434522300628303394596,
    955                         0.101142066918717499027074231447392,
    956                         0.000000000000000000000000000000000
    957                     ],
    958
    959                 /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule.
    960                 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */
    961                 wg =     /* weights of the 10-point gauss rule */
    962                     [
    963                         0.030753241996117268354628393577204,
    964                         0.070366047488108124709267416450667,
    965                         0.107159220467171935011869546685869,
    966                         0.139570677926154314447804794511028,
    967                         0.166269205816993933553200860481209,
    968                         0.186161000015562211026800561866423,
    969                         0.198431485327111576456118326443839,
    970                         0.202578241925561272880620199967519
    971                     ],
    972
    973                 wgk =   /* weights of the 21-point kronrod rule */
    974                     [
    975                         0.005377479872923348987792051430128,
    976                         0.015007947329316122538374763075807,
    977                         0.025460847326715320186874001019653,
    978                         0.035346360791375846222037948478360,
    979                         0.044589751324764876608227299373280,
    980                         0.053481524690928087265343147239430,
    981                         0.062009567800670640285139230960803,
    982                         0.069854121318728258709520077099147,
    983                         0.076849680757720378894432777482659,
    984                         0.083080502823133021038289247286104,
    985                         0.088564443056211770647275443693774,
    986                         0.093126598170825321225486872747346,
    987                         0.096642726983623678505179907627589,
    988                         0.099173598721791959332393173484603,
    989                         0.100769845523875595044946662617570,
    990                         0.101330007014791549017374792767493
    991                     ];
    992
    993             return this._gaussKronrod(interval, f, 16, xgk, wg, wgk, resultObj);
    994         },
    995
    996         /**
    997          * Generate workspace object for {@link JXG.Math.Numerics.Qag}.
    998          * @param {Array} interval The integration interval, e.g. [0, 3].
    999          * @param {Number} n Max. limit
    1000          * @returns {Object} Workspace object
    1001          *
    1002          * @private
    1003          * @memberof JXG.Math.Numerics
    1004          */
    1005         _workspace: function (interval, n) {
    1006             return {
    1007                 limit: n,
    1008                 size: 0,
    1009                 nrmax: 0,
    1010                 i: 0,
    1011                 alist: [interval[0]],
    1012                 blist: [interval[1]],
    1013                 rlist: [0.0],
    1014                 elist: [0.0],
    1015                 order: [0],
    1016                 level: [0],
    1017
    1018                 qpsrt: function () {
    1019                     var last = this.size - 1,
    1020                         limit = this.limit,
    1021                         errmax, errmin, i, k, top,
    1022                         i_nrmax = this.nrmax,
    1023                         i_maxerr = this.order[i_nrmax];
    1024
    1025                     /* Check whether the list contains more than two error estimates */
    1026                     if (last < 2) {
    1027                         this.order[0] = 0;
    1028                         this.order[1] = 1;
    1029                         this.i = i_maxerr;
    1030                         return;
    1031                     }
    1032
    1033                     errmax = this.elist[i_maxerr];
    1034
    1035                     /* This part of the routine is only executed if, due to a difficult
    1036                         integrand, subdivision increased the error estimate. In the normal
    1037                         case the insert procedure should start after the nrmax-th largest
    1038                         error estimate. */
    1039                     while (i_nrmax > 0 && errmax > this.elist[this.order[i_nrmax - 1]]) {
    1040                         this.order[i_nrmax] = this.order[i_nrmax - 1];
    1041                         i_nrmax--;
    1042                     }
    1043
    1044                     /* Compute the number of elements in the list to be maintained in
    1045                         descending order. This number depends on the number of
    1046                         subdivisions still allowed. */
    1047                     if (last < (limit / 2 + 2)) {
    1048                         top = last;
    1049                     } else {
    1050                         top = limit - last + 1;
    1051                     }
    1052
    1053                     /* Insert errmax by traversing the list top-down, starting
    1054                         comparison from the element elist(order(i_nrmax+1)). */
    1055                     i = i_nrmax + 1;
    1056
    1057                     /* The order of the tests in the following line is important to
    1058                         prevent a segmentation fault */
    1059                     while (i < top && errmax < this.elist[this.order[i]]) {
    1060                         this.order[i - 1] = this.order[i];
    1061                         i++;
    1062                     }
    1063
    1064                     this.order[i - 1] = i_maxerr;
    1065
    1066                     /* Insert errmin by traversing the list bottom-up */
    1067                     errmin = this.elist[last];
    1068                     k = top - 1;
    1069
    1070                     while (k > i - 2 && errmin >= this.elist[this.order[k]]) {
    1071                         this.order[k + 1] = this.order[k];
    1072                         k--;
    1073                     }
    1074
    1075                     this.order[k + 1] = last;
    1076
    1077                     /* Set i_max and e_max */
    1078                     i_maxerr = this.order[i_nrmax];
    1079                     this.i = i_maxerr;
    1080                     this.nrmax = i_nrmax;
    1081                 },
    1082
    1083                 set_initial_result: function (result, error) {
    1084                     this.size = 1;
    1085                     this.rlist[0] = result;
    1086                     this.elist[0] = error;
    1087                 },
    1088
    1089                 update: function (a1, b1, area1, error1, a2, b2, area2, error2) {
    1090                     var i_max = this.i,
    1091                         i_new = this.size,
    1092                         new_level = this.level[this.i] + 1;
    1093
    1094                     /* append the newly-created intervals to the list */
    1095
    1096                     if (error2 > error1) {
    1097                         this.alist[i_max] = a2;        /* blist[maxerr] is already == b2 */
    1098                         this.rlist[i_max] = area2;
    1099                         this.elist[i_max] = error2;
    1100                         this.level[i_max] = new_level;
    1101
    1102                         this.alist[i_new] = a1;
    1103                         this.blist[i_new] = b1;
    1104                         this.rlist[i_new] = area1;
    1105                         this.elist[i_new] = error1;
    1106                         this.level[i_new] = new_level;
    1107                     } else {
    1108                         this.blist[i_max] = b1;        /* alist[maxerr] is already == a1 */
    1109                         this.rlist[i_max] = area1;
    1110                         this.elist[i_max] = error1;
    1111                         this.level[i_max] = new_level;
    1112
    1113                         this.alist[i_new] = a2;
    1114                         this.blist[i_new] = b2;
    1115                         this.rlist[i_new] = area2;
    1116                         this.elist[i_new] = error2;
    1117                         this.level[i_new] = new_level;
    1118                     }
    1119
    1120                     this.size++;
    1121
    1122                     if (new_level > this.maximum_level) {
    1123                         this.maximum_level = new_level;
    1124                     }
    1125
    1126                     this.qpsrt();
    1127                 },
    1128
    1129                 retrieve: function() {
    1130                     var i = this.i;
    1131                     return {
    1132                         a: this.alist[i],
    1133                         b: this.blist[i],
    1134                         r: this.rlist[i],
    1135                         e: this.elist[i]
    1136                     };
    1137                 },
    1138
    1139                 sum_results: function () {
    1140                     var nn = this.size,
    1141                         k,
    1142                         result_sum = 0.0;
    1143
    1144                     for (k = 0; k < nn; k++) {
    1145                         result_sum += this.rlist[k];
    1146                     }
    1147
    1148                     return result_sum;
    1149                 },
    1150
    1151                 subinterval_too_small: function (a1, a2,  b2) {
    1152                     var e = 2.2204460492503131e-16,
    1153                         u = 2.2250738585072014e-308,
    1154                         tmp = (1 + 100 * e) * (Math.abs(a2) + 1000 * u);
    1155
    1156                     return Math.abs(a1) <= tmp && Math.abs(b2) <= tmp;
    1157                 }
    1158
    1159             };
    1160         },
    1161
    1162         /**
    1163          * Quadrature algorithm qag from QUADPACK.
    1164          * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15},
    1165          * {@link JXG.Math.Numerics.GaussKronrod21},
    1166          * {@link JXG.Math.Numerics.GaussKronrod31}.
    1167          *
    1168          * @param {Array} interval The integration interval, e.g. [0, 3].
    1169          * @param {function} f A function which takes one argument of type number and returns a number.
    1170          * @param {Object} [config] The algorithm setup. Accepted propert are max. recursion limit of type number,
    1171          * and epsrel and epsabs, the relative and absolute required precision of type number. Further,
    1172          * q the internal quadrature sub-algorithm of type function.
    1173          * @param {Number} [config.limit=15]
    1174          * @param {Number} [config.epsrel=0.0000001]
    1175          * @param {Number} [config.epsabs=0.0000001]
    1176          * @param {Number} [config.q=JXG.Math.Numerics.GaussKronrod15]
    1177          * @returns {Number} Integral value of f over interval
    1178          *
    1179          * @example
    1180          * function f(x) {
    1181          *   return x*x;
    1182          * }
    1183          *
    1184          * // calculates integral of <tt>f</tt> from 0 to 2.
    1185          * var area1 = JXG.Math.Numerics.Qag([0, 2], f);
    1186          *
    1187          * // the same with an anonymous function
    1188          * var area2 = JXG.Math.Numerics.Qag([0, 2], function (x) { return x*x; });
    1189          *
    1190          * // use JXG.Math.Numerics.GaussKronrod31 rule as sub-algorithm.
    1191          * var area3 = JXG.Math.Numerics.Quag([0, 2], f,
    1192          *                                   {q: JXG.Math.Numerics.GaussKronrod31});
    1193          * @memberof JXG.Math.Numerics
    1194          */
    1195         Qag: function (interval, f, config) {
    1196             var DBL_EPS = 2.2204460492503131e-16,
    1197                 ws = this._workspace(interval, 1000),
    1198
    1199                 limit = config && Type.isNumber(config.limit) ? config.limit : 15,
    1200                 epsrel = config && Type.isNumber(config.epsrel) ? config.epsrel : 0.0000001,
    1201                 epsabs = config && Type.isNumber(config.epsabs) ? config.epsabs : 0.0000001,
    1202                 q = config && Type.isFunction(config.q) ? config.q : this.GaussKronrod15,
    1203
    1204                 resultObj = {},
    1205                 area, errsum,
    1206                 result0, abserr0, resabs0, resasc0,
    1207                 result,
    1208                 tolerance,
    1209                 iteration = 0,
    1210                 roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0,
    1211                 round_off,
    1212
    1213                 a1, b1, a2, b2,
    1214                 a_i, b_i, r_i, e_i,
    1215                 area1 = 0, area2 = 0, area12 = 0,
    1216                 error1 = 0, error2 = 0, error12 = 0,
    1217                 resasc1, resasc2,
    1218                 // resabs1, resabs2,
    1219                 wsObj,
    1220                 delta;
    1221
    1222
    1223             if (limit > ws.limit) {
    1224                 JXG.warn('iteration limit exceeds available workspace');
    1225             }
    1226             if (epsabs <= 0 && (epsrel < 50 * Mat.eps || epsrel < 0.5e-28)) {
    1227                 JXG.warn('tolerance cannot be acheived with given epsabs and epsrel');
    1228             }
    1229
    1230             result0 = q.apply(this, [interval, f, resultObj]);
    1231             abserr0 = resultObj.abserr;
    1232             resabs0 = resultObj.resabs;
    1233             resasc0 = resultObj.resasc;
    1234
    1235             ws.set_initial_result(result0, abserr0);
    1236             tolerance = Math.max(epsabs, epsrel * Math.abs(result0));
    1237             round_off = 50 * DBL_EPS * resabs0;
    1238
    1239             if (abserr0 <= round_off && abserr0 > tolerance) {
    1240                 result = result0;
    1241                 // abserr = abserr0;
    1242
    1243                 JXG.warn('cannot reach tolerance because of roundoff error on first attempt');
    1244                 return -Infinity;
    1245             }
    1246
    1247             if ((abserr0 <= tolerance && abserr0 !== resasc0) || abserr0 === 0.0) {
    1248                 result = result0;
    1249                 // abserr = abserr0;
    1250
    1251                 return result;
    1252             }
    1253
    1254             if (limit === 1) {
    1255                 result = result0;
    1256                 // abserr = abserr0;
    1257
    1258                 JXG.warn('a maximum of one iteration was insufficient');
    1259                 return -Infinity;
    1260             }
    1261
    1262             area = result0;
    1263             errsum = abserr0;
    1264             iteration = 1;
    1265
    1266             do {
    1267                 area1 = 0;
    1268                 area2 = 0;
    1269                 area12 = 0;
    1270                 error1 = 0;
    1271                 error2 = 0;
    1272                 error12 = 0;
    1273
    1274                 /* Bisect the subinterval with the largest error estimate */
    1275                 wsObj = ws.retrieve();
    1276                 a_i = wsObj.a;
    1277                 b_i = wsObj.b;
    1278                 r_i = wsObj.r;
    1279                 e_i = wsObj.e;
    1280
    1281                 a1 = a_i;
    1282                 b1 = 0.5 * (a_i + b_i);
    1283                 a2 = b1;
    1284                 b2 = b_i;
    1285
    1286                 area1 = q.apply(this, [[a1, b1], f, resultObj]);
    1287                 error1 = resultObj.abserr;
    1288                 // resabs1 = resultObj.resabs;
    1289                 resasc1 = resultObj.resasc;
    1290
    1291                 area2 = q.apply(this, [[a2, b2], f, resultObj]);
    1292                 error2 = resultObj.abserr;
    1293                 // resabs2 = resultObj.resabs;
    1294                 resasc2 = resultObj.resasc;
    1295
    1296                 area12 = area1 + area2;
    1297                 error12 = error1 + error2;
    1298
    1299                 errsum += (error12 - e_i);
    1300                 area += area12 - r_i;
    1301
    1302                 if (resasc1 !== error1 && resasc2 !== error2) {
    1303                     delta = r_i - area12;
    1304                     if (Math.abs(delta) <= 1.0e-5 * Math.abs(area12) && error12 >= 0.99 * e_i) {
    1305                         roundoff_type1++;
    1306                     }
    1307                     if (iteration >= 10 && error12 > e_i) {
    1308                         roundoff_type2++;
    1309                     }
    1310                 }
    1311
    1312                 tolerance = Math.max(epsabs, epsrel * Math.abs(area));
    1313
    1314                 if (errsum > tolerance) {
    1315                     if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
    1316                         error_type = 2;   /* round off error */
    1317                     }
    1318
    1319                 /* set error flag in the case of bad integrand behaviour at
    1320                     a point of the integration range */
    1321
    1322                     if (ws.subinterval_too_small(a1, a2, b2)) {
    1323                         error_type = 3;
    1324                     }
    1325                 }
    1326
    1327                 ws.update(a1, b1, area1, error1, a2, b2, area2, error2);
    1328                 wsObj = ws.retrieve();
    1329                 a_i = wsObj.a_i;
    1330                 b_i = wsObj.b_i;
    1331                 r_i = wsObj.r_i;
    1332                 e_i = wsObj.e_i;
    1333
    1334                 iteration++;
    1335
    1336             } while (iteration < limit && !error_type && errsum > tolerance);
    1337
    1338             result = ws.sum_results();
    1339             // abserr = errsum;
    1340 /*
    1341   if (errsum <= tolerance)
    1342     {
    1343       return GSL_SUCCESS;
    1344     }
    1345   else if (error_type == 2)
    1346     {
    1347       GSL_ERROR ("roundoff error prevents tolerance from being achieved",
    1348                  GSL_EROUND);
    1349     }
    1350   else if (error_type == 3)
    1351     {
    1352       GSL_ERROR ("bad integrand behavior found in the integration interval",
    1353                  GSL_ESING);
    1354     }
    1355   else if (iteration == limit)
    1356     {
    1357       GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
    1358     }
    1359   else
    1360     {
    1361       GSL_ERROR ("could not integrate function", GSL_EFAILED);
    1362     }
    1363 */
    1364
    1365             return result;
    1366         },
    1367
    1368         /**
    1369          * Integral of function f over interval.
    1370          * @param {Array} interval The integration interval, e.g. [0, 3].
    1371          * @param {function} f A function which takes one argument of type number and returns a number.
    1372          * @returns {Number} The value of the integral of f over interval
    1373          * @see JXG.Math.Numerics.NewtonCotes
    1374          * @see JXG.Math.Numerics.Romberg
    1375          * @see JXG.Math.Numerics.Qag
    1376          * @memberof JXG.Math.Numerics
    1377          */
    1378         I: function (interval, f) {
    1379             // return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'});
    1380             // return this.Romberg(interval, f, {max_iterations: 20, eps: 0.0000001});
    1381             return this.Qag(interval, f, {q: this.GaussKronrod15, limit: 15, epsrel: 0.0000001, epsabs: 0.0000001});
    1382         },
    1383
    1384         /**
    1385          * Newton's method to find roots of a funtion in one variable.
    1386          * @param {function} f We search for a solution of f(x)=0.
    1387          * @param {Number} x initial guess for the root, i.e. start value.
    1388          * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
    1389          * the function is a method of an object and contains a reference to its parent object via "this".
    1390          * @returns {Number} A root of the function f.
    1391          * @memberof JXG.Math.Numerics
    1392          */
    1393         Newton: function (f, x, context) {
    1394             var df,
    1395                 i = 0,
    1396                 h = Mat.eps,
    1397                 newf = f.apply(context, [x]);
    1398                 // nfev = 1;
    1399
    1400             // For compatibility
    1401             if (Type.isArray(x)) {
    1402                 x = x[0];
    1403             }
    1404
    1405             while (i < 50 && Math.abs(newf) > h) {
    1406                 df = this.D(f, context)(x);
    1407                 // nfev += 2;
    1408
    1409                 if (Math.abs(df) > h) {
    1410                     x -= newf / df;
    1411                 } else {
    1412                     x += (Math.random() * 0.2 - 1.0);
    1413                 }
    1414
    1415                 newf = f.apply(context, [x]);
    1416                 // nfev += 1;
    1417                 i += 1;
    1418             }
    1419
    1420             return x;
    1421         },
    1422
    1423         /**
    1424          * Abstract method to find roots of univariate functions, which - for the time being -
    1425          * is an alias for {@link JXG.Math.Numerics.chandrupatla}.
    1426          * @param {function} f We search for a solution of f(x)=0.
    1427          * @param {Number|Array} x initial guess for the root, i.e. starting value, or start interval enclosing the root.
    1428          * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
    1429          * the function is a method of an object and contains a reference to its parent object via "this".
    1430          * @returns {Number} A root of the function f.
    1431          *
    1432          * @see JXG.Math.Numerics.chandrupatla
    1433          * @see JXG.Math.Numerics.fzero
    1434          * @memberof JXG.Math.Numerics
    1435          */
    1436         root: function (f, x, context) {
    1437             //return this.fzero(f, x, context);
    1438             return this.chandrupatla(f, x, context);
    1439         },
    1440
    1441         /**
    1442          * Compute an intersection of the curves c1 and c2
    1443          * with a generalized Newton method.
    1444          * We want to find values t1, t2 such that
    1445          * c1(t1) = c2(t2), i.e.
    1446          * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0).
    1447          * We set
    1448          * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2))
    1449          *
    1450          * The Jacobian J is defined by
    1451          * J = (a, b)
    1452          *     (c, d)
    1453          * where
    1454          * a = c1_x'(t1)
    1455          * b = -c2_x'(t2)
    1456          * c = c1_y'(t1)
    1457          * d = -c2_y'(t2)
    1458          *
    1459          * The inverse J^(-1) of J is equal to
    1460          *  (d, -b)/
    1461          *  (-c, a) / (ad-bc)
    1462          *
    1463          * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f).
    1464          * If the function meetCurveCurve possesses the properties
    1465          * t1memo and t2memo then these are taken as start values
    1466          * for the Newton algorithm.
    1467          * After stopping of the Newton algorithm the values of t1 and t2 are stored in
    1468          * t1memo and t2memo.
    1469          *
    1470          * @param {JXG.Curve} c1 Curve, Line or Circle
    1471          * @param {JXG.Curve} c2 Curve, Line or Circle
    1472          * @param {Number} t1ini start value for t1
    1473          * @param {Number} t2ini start value for t2
    1474          * @returns {JXG.Coords} intersection point
    1475          * @memberof JXG.Math.Numerics
    1476          */
    1477         generalizedNewton: function (c1, c2, t1ini, t2ini) {
    1478             var t1, t2,
    1479                 a, b, c, d, disc,
    1480                 e, f, F,
    1481                 D00, D01,
    1482                 D10, D11,
    1483                 count = 0;
    1484
    1485             if (this.generalizedNewton.t1memo) {
    1486                 t1 = this.generalizedNewton.t1memo;
    1487                 t2 = this.generalizedNewton.t2memo;
    1488             } else {
    1489                 t1 = t1ini;
    1490                 t2 = t2ini;
    1491             }
    1492
    1493             e = c1.X(t1) - c2.X(t2);
    1494             f = c1.Y(t1) - c2.Y(t2);
    1495             F = e * e + f * f;
    1496
    1497             D00 = this.D(c1.X, c1);
    1498             D01 = this.D(c2.X, c2);
    1499             D10 = this.D(c1.Y, c1);
    1500             D11 = this.D(c2.Y, c2);
    1501
    1502             while (F > Mat.eps && count < 10) {
    1503                 a = D00(t1);
    1504                 b = -D01(t2);
    1505                 c = D10(t1);
    1506                 d = -D11(t2);
    1507                 disc = a * d - b * c;
    1508                 t1 -= (d * e - b * f) / disc;
    1509                 t2 -= (a * f - c * e) / disc;
    1510                 e = c1.X(t1) - c2.X(t2);
    1511                 f = c1.Y(t1) - c2.Y(t2);
    1512                 F = e * e + f * f;
    1513                 count += 1;
    1514             }
    1515
    1516             this.generalizedNewton.t1memo = t1;
    1517             this.generalizedNewton.t2memo = t2;
    1518
    1519             if (Math.abs(t1) < Math.abs(t2)) {
    1520                 return [c1.X(t1), c1.Y(t1)];
    1521             }
    1522
    1523             return [c2.X(t2), c2.Y(t2)];
    1524         },
    1525
    1526         /**
    1527          * Returns the Lagrange polynomials for curves with equidistant nodes, see
    1528          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
    1529          * SIAM Review, Vol 46, No 3, (2004) 501-517.
    1530          * The graph of the parametric curve [x(t),y(t)] runs through the given points.
    1531          * @param {Array} p Array of JXG.Points
    1532          * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve
    1533          * f(t) = (x(t), y(t)), a number x1 (which equals 0) and a function x2 defining the curve's domain.
    1534          * That means the curve is defined between x1 and x2(). x2 returns the (length of array p minus one).
    1535          * @memberof JXG.Math.Numerics
    1536          *
    1537          * @example
    1538          * var p = [];
    1539          *
    1540          * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'});
    1541          * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''});
    1542          * p[2] = board.create('point', [1, 4], {size:2, name: ''});
    1543          * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'});
    1544          *
    1545          * // Curve
    1546          * var fg = JXG.Math.Numerics.Neville(p);
    1547          * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5});
    1548          *
    1549          * </pre><div id="JXG88a8b3a8-6561-44f5-a678-76bca13fd484" class="jxgbox" style="width: 300px; height: 300px;"></div>
    1550          * <script type="text/javascript">
    1551          *     (function() {
    1552          *         var board = JXG.JSXGraph.initBoard('JXG88a8b3a8-6561-44f5-a678-76bca13fd484',
    1553          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
    1554          *     var p = [];
    1555          *
    1556          *     p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'});
    1557          *     p[1] = board.create('point', [-1.5, 5], {size:2, name: ''});
    1558          *     p[2] = board.create('point', [1, 4], {size:2, name: ''});
    1559          *     p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'});
    1560          *
    1561          *     // Curve
    1562          *     var fg = JXG.Math.Numerics.Neville(p);
    1563          *     var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5});
    1564          *
    1565          *     })();
    1566          *
    1567          * </script><pre>
    1568          *
    1569          */
    1570         Neville: function (p) {
    1571             var w = [],
    1572                 /** @ignore */
    1573                 makeFct = function (fun) {
    1574                     return function (t, suspendedUpdate) {
    1575                         var i, d, s,
    1576                             bin = Mat.binomial,
    1577                             len = p.length,
    1578                             len1 = len - 1,
    1579                             num = 0.0,
    1580                             denom = 0.0;
    1581
    1582                         if (!suspendedUpdate) {
    1583                             s = 1;
    1584                             for (i = 0; i < len; i++) {
    1585                                 w[i] = bin(len1, i) * s;
    1586                                 s *= (-1);
    1587                             }
    1588                         }
    1589
    1590                         d = t;
    1591
    1592                         for (i = 0; i < len; i++) {
    1593                             if (d === 0) {
    1594                                 return p[i][fun]();
    1595                             }
    1596                             s = w[i] / d;
    1597                             d -= 1;
    1598                             num += p[i][fun]() * s;
    1599                             denom += s;
    1600                         }
    1601                         return num / denom;
    1602                     };
    1603                 },
    1604
    1605                 xfct = makeFct('X'),
    1606                 yfct = makeFct('Y');
    1607
    1608             return [xfct, yfct, 0, function () {
    1609                 return p.length - 1;
    1610             }];
    1611         },
    1612
    1613         /**
    1614          * Calculates second derivatives at the given knots.
    1615          * @param {Array} x x values of knots
    1616          * @param {Array} y y values of knots
    1617          * @returns {Array} Second derivatives of the interpolated function at the knots.
    1618          * @see #splineEval
    1619          * @memberof JXG.Math.Numerics
    1620          */
    1621         splineDef: function (x, y) {
    1622             var pair, i, l,
    1623                 n = Math.min(x.length, y.length),
    1624                 diag = [],
    1625                 z = [],
    1626                 data = [],
    1627                 dx = [],
    1628                 delta = [],
    1629                 F = [];
    1630
    1631             if (n === 2) {
    1632                 return [0, 0];
    1633             }
    1634
    1635             for (i = 0; i < n; i++) {
    1636                 pair = {X: x[i], Y: y[i]};
    1637                 data.push(pair);
    1638             }
    1639             data.sort(function (a, b) {
    1640                 return a.X - b.X;
    1641             });
    1642             for (i = 0; i < n; i++) {
    1643                 x[i] = data[i].X;
    1644                 y[i] = data[i].Y;
    1645             }
    1646
    1647             for (i = 0; i < n - 1; i++) {
    1648                 dx.push(x[i + 1] - x[i]);
    1649             }
    1650             for (i = 0; i < n - 2; i++) {
    1651                 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i]));
    1652             }
    1653
    1654             // ForwardSolve
    1655             diag.push(2 * (dx[0] + dx[1]));
    1656             z.push(delta[0]);
    1657
    1658             for (i = 0; i < n - 3; i++) {
    1659                 l = dx[i + 1] / diag[i];
    1660                 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]);
    1661                 z.push(delta[i + 1] - l * z[i]);
    1662             }
    1663
    1664             // BackwardSolve
    1665             F[n - 3] = z[n - 3] / diag[n - 3];
    1666             for (i = n - 4; i >= 0; i--) {
    1667                 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i];
    1668             }
    1669
    1670             // Generate f''-Vector
    1671             for (i = n - 3; i >= 0; i--) {
    1672                 F[i + 1] = F[i];
    1673             }
    1674
    1675             // natural cubic spline
    1676             F[0] = 0;
    1677             F[n - 1] = 0;
    1678
    1679             return F;
    1680         },
    1681
    1682         /**
    1683          * Evaluate points on spline.
    1684          * @param {Number,Array} x0 A single float value or an array of values to evaluate
    1685          * @param {Array} x x values of knots
    1686          * @param {Array} y y values of knots
    1687          * @param {Array} F Second derivatives at knots, calculated by {@link JXG.Math.Numerics.splineDef}
    1688          * @see #splineDef
    1689          * @returns {Number,Array} A single value or an array, depending on what is given as x0.
    1690          * @memberof JXG.Math.Numerics
    1691          */
    1692         splineEval: function (x0, x, y, F) {
    1693             var i, j, a, b, c, d, x_,
    1694                 n = Math.min(x.length, y.length),
    1695                 l = 1,
    1696                 asArray = false,
    1697                 y0 = [];
    1698
    1699             // number of points to be evaluated
    1700             if (Type.isArray(x0)) {
    1701                 l = x0.length;
    1702                 asArray = true;
    1703             } else {
    1704                 x0 = [x0];
    1705             }
    1706
    1707             for (i = 0; i < l; i++) {
    1708                 // is x0 in defining interval?
    1709                 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) {
    1710                     return NaN;
    1711                 }
    1712
    1713                 // determine part of spline in which x0 lies
    1714                 for (j = 1; j < n; j++) {
    1715                     if (x0[i] <= x[j]) {
    1716                         break;
    1717                     }
    1718                 }
    1719
    1720                 j -= 1;
    1721
    1722                 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1];
    1723                 // determine the coefficients of the polynomial in this interval
    1724                 a = y[j];
    1725                 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]);
    1726                 c = F[j] / 2;
    1727                 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j]));
    1728                 // evaluate x0[i]
    1729                 x_ = x0[i] - x[j];
    1730                 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_);
    1731                 y0.push(a + (b + (c + d * x_) * x_) * x_);
    1732             }
    1733
    1734             if (asArray) {
    1735                 return y0;
    1736             }
    1737
    1738             return y0[0];
    1739         },
    1740
    1741         /**
    1742          * Generate a string containing the function term of a polynomial.
    1743          * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i.
    1744          * @param {Number} deg Degree of the polynomial
    1745          * @param {String} varname Name of the variable (usually 'x')
    1746          * @param {Number} prec Precision
    1747          * @returns {String} A string containg the function term of the polynomial.
    1748          * @memberof JXG.Math.Numerics
    1749          */
    1750         generatePolynomialTerm: function (coeffs, deg, varname, prec) {
    1751             var i, t = [];
    1752
    1753             for (i = deg; i >= 0; i--) {
    1754                 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']);
    1755
    1756                 if (i > 1) {
    1757                     t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']);
    1758                 } else if (i === 1) {
    1759                     t = t.concat(['*', varname, ' + ']);
    1760                 }
    1761             }
    1762
    1763             return t.join('');
    1764         },
    1765
    1766         /**
    1767          * Computes the polynomial through a given set of coordinates in Lagrange form.
    1768          * Returns the Lagrange polynomials, see
    1769          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
    1770          * SIAM Review, Vol 46, No 3, (2004) 501-517.
    1771          * <p>
    1772          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
    1773          * @param {Array} p Array of JXG.Points
    1774          * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points.
    1775          * @memberof JXG.Math.Numerics
    1776          *
    1777          * @example
    1778          * var p = [];
    1779          * p[0] = board.create('point', [-1,2], {size:4});
    1780          * p[1] = board.create('point', [0,3], {size:4});
    1781          * p[2] = board.create('point', [1,1], {size:4});
    1782          * p[3] = board.create('point', [3,-1], {size:4});
    1783          * var f = JXG.Math.Numerics.lagrangePolynomial(p);
    1784          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1785          *
    1786          * </pre><div id="JXGc058aa6b-74d4-41e1-af94-df06169a2d89" class="jxgbox" style="width: 300px; height: 300px;"></div>
    1787          * <script type="text/javascript">
    1788          *     (function() {
    1789          *         var board = JXG.JSXGraph.initBoard('JXGc058aa6b-74d4-41e1-af94-df06169a2d89',
    1790          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
    1791          *     var p = [];
    1792          *     p[0] = board.create('point', [-1,2], {size:4});
    1793          *     p[1] = board.create('point', [0,3], {size:4});
    1794          *     p[2] = board.create('point', [1,1], {size:4});
    1795          *     p[3] = board.create('point', [3,-1], {size:4});
    1796          *     var f = JXG.Math.Numerics.lagrangePolynomial(p);
    1797          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1798          *
    1799          *     })();
    1800          *
    1801          * </script><pre>
    1802          *
    1803          * @example
    1804          * var points = [];
    1805          * points[0] = board.create('point', [-1,2], {size:4});
    1806          * points[1] = board.create('point', [0, 0], {size:4});
    1807          * points[2] = board.create('point', [2, 1], {size:4});
    1808          *
    1809          * var f = JXG.Math.Numerics.lagrangePolynomial(points);
    1810          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1811          * var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
    1812          *
    1813          * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eecbb" class="jxgbox" style="width: 300px; height: 300px;"></div>
    1814          * <script type="text/javascript">
    1815          *     (function() {
    1816          *         var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eecbb',
    1817          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
    1818          *     var points = [];
    1819          *     points[0] = board.create('point', [-1,2], {size:4});
    1820          *     points[1] = board.create('point', [0, 0], {size:4});
    1821          *     points[2] = board.create('point', [2, 1], {size:4});
    1822          *
    1823          *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
    1824          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1825          *     var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
    1826          *
    1827          *     })();
    1828          *
    1829          * </script><pre>
    1830          *
    1831          */
    1832         lagrangePolynomial: function (p) {
    1833             var w = [],
    1834                 that = this,
    1835                 /** @ignore */
    1836                 fct = function (x, suspendedUpdate) {
    1837                     var i, // j,
    1838                         k, xi, s, //M,
    1839                         len = p.length,
    1840                         num = 0,
    1841                         denom = 0;
    1842
    1843                     if (!suspendedUpdate) {
    1844                         for (i = 0; i < len; i++) {
    1845                             w[i] = 1.0;
    1846                             xi = p[i].X();
    1847
    1848                             for (k = 0; k < len; k++) {
    1849                                 if (k !== i) {
    1850                                     w[i] *= (xi - p[k].X());
    1851                                 }
    1852                             }
    1853
    1854                             w[i] = 1 / w[i];
    1855                         }
    1856
    1857                         // M = [];
    1858                         // for (k = 0; k < len; k++) {
    1859                         //     M.push([1]);
    1860                         // }
    1861                     }
    1862
    1863                     for (i = 0; i < len; i++) {
    1864                         xi = p[i].X();
    1865
    1866                         if (x === xi) {
    1867                             return p[i].Y();
    1868                         }
    1869
    1870                         s = w[i] / (x - xi);
    1871                         denom += s;
    1872                         num += s * p[i].Y();
    1873                     }
    1874
    1875                     return num / denom;
    1876                 };
    1877
    1878             /**
    1879              * Get the term of the Lagrange polynomial as string.
    1880              * Calls {@link JXG.Math.Numerics#lagrangePolynomialTerm}.
    1881              *
    1882              * @name JXG.Math.Numerics#lagrangePolynomial.getTerm
    1883              * @param {Number} digits Number of digits of the coefficients
    1884              * @param {String} param Variable name
    1885              * @param {String} dot Dot symbol
    1886              * @returns {String} containing the term of Lagrange polynomial as string.
    1887              * @see JXG.Math.Numerics#lagrangePolynomialTerm
    1888              * @example
    1889              * var points = [];
    1890              * points[0] = board.create('point', [-1,2], {size:4});
    1891              * points[1] = board.create('point', [0, 0], {size:4});
    1892              * points[2] = board.create('point', [2, 1], {size:4});
    1893              *
    1894              * var f = JXG.Math.Numerics.lagrangePolynomial(points);
    1895              * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1896              * var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
    1897              *
    1898              * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eeccf" class="jxgbox" style="width: 300px; height: 300px;"></div>
    1899              * <script type="text/javascript">
    1900              *     (function() {
    1901              *         var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eeccf',
    1902              *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
    1903              *     var points = [];
    1904              *     points[0] = board.create('point', [-1,2], {size:4});
    1905              *     points[1] = board.create('point', [0, 0], {size:4});
    1906              *     points[2] = board.create('point', [2, 1], {size:4});
    1907              *
    1908              *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
    1909              *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1910              *     var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
    1911              *
    1912              *     })();
    1913              *
    1914              * </script><pre>
    1915              *
    1916              */
    1917             fct.getTerm = function(digits, param, dot) {
    1918                 return that.lagrangePolynomialTerm(p, digits, param, dot)();
    1919             };
    1920
    1921             return fct;
    1922         },
    1923         // fct.getTerm = that.lagrangePolynomialTerm(p, 2, 'x');
    1924
    1925         /**
    1926          * Determine the Lagrange polynomial through an array of points and
    1927          * return the term of the polynomial as string.
    1928          *
    1929          * @param {Array} points Array of JXG.Points
    1930          * @param {Number} digits Number of decimal digits of the coefficients
    1931          * @param {String} param Name of the parameter. Default: 'x'.
    1932          * @param {String} dot Multiplication symbol. Default: ' * '.
    1933          * @returns {String} containing the Lagrange polynomial through
    1934          *    the supplied points.
    1935          * @memberof JXG.Math.Numerics
    1936          *
    1937          * @example
    1938          * var points = [];
    1939          * points[0] = board.create('point', [-1,2], {size:4});
    1940          * points[1] = board.create('point', [0, 0], {size:4});
    1941          * points[2] = board.create('point', [2, 1], {size:4});
    1942          *
    1943          * var f = JXG.Math.Numerics.lagrangePolynomial(points);
    1944          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1945          *
    1946          * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * ');
    1947          * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16});
    1948          *
    1949          * </pre><div id="JXGd45e9e96-7526-486d-aa43-e1178d5f2baa" class="jxgbox" style="width: 300px; height: 300px;"></div>
    1950          * <script type="text/javascript">
    1951          *     (function() {
    1952          *         var board = JXG.JSXGraph.initBoard('JXGd45e9e96-7526-486d-aa43-e1178d5f2baa',
    1953          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
    1954          *     var points = [];
    1955          *     points[0] = board.create('point', [-1,2], {size:4});
    1956          *     points[1] = board.create('point', [0, 0], {size:4});
    1957          *     points[2] = board.create('point', [2, 1], {size:4});
    1958          *
    1959          *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
    1960          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
    1961          *
    1962          *     var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * ');
    1963          *     var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16});
    1964          *
    1965          *     })();
    1966          *
    1967          * </script><pre>
    1968          *
    1969          */
    1970         lagrangePolynomialTerm: function(points, digits, param, dot) {
    1971             return function() {
    1972                 var len = points.length,
    1973                     zeroes = [],
    1974                     coeffs = [],
    1975                     coeffs_sum = [],
    1976                     isLeading = true,
    1977                     n, t,
    1978                     i, j, c, p;
    1979
    1980                 param = param || 'x';
    1981                 if (dot === undefined) {
    1982                     dot = ' * ';
    1983                 }
    1984
    1985                 n = len - 1;  // (Max) degree of the polynomial
    1986                 for (j = 0; j < len; j++) {
    1987                     coeffs_sum[j] = 0;
    1988                 }
    1989
    1990                 for (i = 0; i < len; i++) {
    1991                     c = points[i].Y();
    1992                     p = points[i].X();
    1993                     zeroes = [];
    1994                     for (j = 0; j < len; j++) {
    1995                         if (j !== i) {
    1996                             c /= p - points[j].X();
    1997                             zeroes.push(points[j].X());
    1998                         }
    1999                     }
    2000                     coeffs = [1].concat(Mat.Vieta(zeroes));
    2001                     for (j = 0; j < coeffs.length; j++) {
    2002                         coeffs_sum[j] += (j%2===1?(-1):1) * coeffs[j] * c;
    2003                     }
    2004                 }
    2005
    2006                 t = '';
    2007                 for (j = 0; j < coeffs_sum.length; j++) {
    2008                     c = coeffs_sum[j];
    2009                     if (Math.abs(c) < Mat.eps) {
    2010                         continue;
    2011                     }
    2012                     if (JXG.exists(digits)) {
    2013                         c = Env._round10(c, -digits);
    2014                     }
    2015                     if (isLeading) {
    2016                         t += (c > 0) ? (c) : ('-' + (-c));
    2017                         isLeading = false;
    2018                     } else {
    2019                         t += (c > 0) ? (' + ' + c) : (' - ' + (-c));
    2020                     }
    2021
    2022                     if (n - j > 1) {
    2023                         t += dot + param + '^' + (n - j);
    2024                     } else if (n - j === 1) {
    2025                         t += dot + param;
    2026                     }
    2027                 }
    2028                 return t; // board.jc.manipulate('f = map(x) -> ' + t + ';');
    2029             };
    2030         },
    2031
    2032         /**
    2033          * Determine the coefficients of a cardinal spline polynom, See
    2034          * http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections
    2035          * @param  {Number} x1 point 1
    2036          * @param  {Number} x2 point 2
    2037          * @param  {Number} t1 tangent slope 1
    2038          * @param  {Number} t2 tangent slope 2
    2039          * @return {Array}    coefficents array c for the polynomial t maps to
    2040          * c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t
    2041          */
    2042         _initCubicPoly: function(x1, x2, t1, t2) {
    2043             return [
    2044                 x1,
    2045                 t1,
    2046                 -3 * x1 + 3 * x2 - 2 * t1 - t2,
    2047                 2 * x1 -  2 * x2 + t1 + t2
    2048             ];
    2049         },
    2050
    2051         /**
    2052          * Computes the cubic cardinal spline curve through a given set of points. The curve
    2053          * is uniformly parametrized.
    2054          * Two artificial control points at the beginning and the end are added.
    2055          *
    2056          * The implementation (especially  the centripetal parametrization) is from
    2057          * http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections .
    2058          * @param {Array} points Array consisting of JXG.Points.
    2059          * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1.
    2060          * tau=1/2 give Catmull-Rom splines.
    2061          * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and
    2062          * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal".
    2063          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
    2064          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value,
    2065          * and a function simply returning the length of the points array
    2066          * minus three.
    2067          * @memberof JXG.Math.Numerics
    2068         */
    2069         CardinalSpline: function (points, tau_param, type) {
    2070             var p,
    2071                 coeffs = [],
    2072                 makeFct,
    2073                 tau, _tau,
    2074                 that = this;
    2075
    2076             if (Type.isFunction(tau_param)) {
    2077                 _tau = tau_param;
    2078             } else {
    2079                 _tau = function () { return tau_param; };
    2080             }
    2081
    2082             if (type === undefined) {
    2083                 type = 'uniform';
    2084             }
    2085
    2086             /** @ignore */
    2087             makeFct = function (which) {
    2088                 return function (t, suspendedUpdate) {
    2089                     var s, c,
    2090                         // control point at the beginning and at the end
    2091                         first, last,
    2092                         t1, t2, dt0, dt1, dt2,
    2093                         // dx, dy,
    2094                         len;
    2095
    2096                     if (points.length < 2) {
    2097                         return NaN;
    2098                     }
    2099
    2100                     if (!suspendedUpdate) {
    2101                         tau = _tau();
    2102
    2103                         // New point list p: [first, points ..., last]
    2104                         first = {
    2105                             X: function () { return 2 * points[0].X() - points[1].X(); },
    2106                             Y: function () { return 2 * points[0].Y() - points[1].Y(); },
    2107                             Dist: function(p) {
    2108                                 var dx = this.X() - p.X(),
    2109                                     dy = this.Y() - p.Y();
    2110                                 return Math.sqrt(dx * dx + dy * dy);
    2111                             }
    2112                         };
    2113
    2114                         last = {
    2115                             X: function () { return 2 * points[points.length - 1].X() - points[points.length - 2].X(); },
    2116                             Y: function () { return 2 * points[points.length - 1].Y() - points[points.length - 2].Y(); },
    2117                             Dist: function(p) {
    2118                                 var dx = this.X() - p.X(),
    2119                                     dy = this.Y() - p.Y();
    2120                                 return Math.sqrt(dx * dx + dy * dy);
    2121                             }
    2122                         };
    2123
    2124                         p = [first].concat(points, [last]);
    2125                         len = p.length;
    2126
    2127                         coeffs[which] = [];
    2128
    2129                         for (s = 0; s < len - 3; s++) {
    2130                             if (type === 'centripetal') {
    2131                                 // The order is important, since p[0].coords === undefined
    2132                                 dt0 = p[s].Dist(p[s + 1]);
    2133                                 dt1 = p[s + 2].Dist(p[s + 1]);
    2134                                 dt2 = p[s + 3].Dist(p[s + 2]);
    2135
    2136                                 dt0 = Math.sqrt(dt0);
    2137                                 dt1 = Math.sqrt(dt1);
    2138                                 dt2 = Math.sqrt(dt2);
    2139
    2140                                 if (dt1 < Mat.eps) { dt1 = 1.0; }
    2141                                 if (dt0 < Mat.eps) { dt0 = dt1; }
    2142                                 if (dt2 < Mat.eps) { dt2 = dt1; }
    2143
    2144                                 t1 = (p[s + 1][which]() - p[s][which]()) / dt0 -
    2145                                      (p[s + 2][which]() - p[s][which]()) / (dt1 + dt0) +
    2146                                      (p[s + 2][which]() - p[s + 1][which]()) / dt1;
    2147
    2148                                 t2 = (p[s + 2][which]() - p[s + 1][which]()) / dt1 -
    2149                                      (p[s + 3][which]() - p[s + 1][which]()) / (dt2 + dt1) +
    2150                                      (p[s + 3][which]() - p[s + 2][which]()) / dt2;
    2151
    2152                                 t1 *= dt1;
    2153                                 t2 *= dt1;
    2154
    2155                                 coeffs[which][s] = that._initCubicPoly(
    2156                                      p[s + 1][which](),
    2157                                      p[s + 2][which](),
    2158                                      tau * t1,
    2159                                      tau * t2
    2160                                 );
    2161                             } else {
    2162                                 coeffs[which][s] = that._initCubicPoly(
    2163                                     p[s + 1][which](),
    2164                                     p[s + 2][which](),
    2165                                     tau * (p[s + 2][which]() - p[s][which]()),
    2166                                     tau * (p[s + 3][which]() - p[s + 1][which]())
    2167                                 );
    2168                             }
    2169                         }
    2170                     }
    2171
    2172                     if (isNaN(t)) {
    2173                         return NaN;
    2174                     }
    2175
    2176                     len = points.length;
    2177                     // This is necessary for our advanced plotting algorithm:
    2178                     if (t <= 0.0) {
    2179                         return points[0][which]();
    2180                     }
    2181                     if (t >= len) {
    2182                         return points[len - 1][which]();
    2183                     }
    2184
    2185                     s = Math.floor(t);
    2186                     if (s === t) {
    2187                         return points[s][which]();
    2188                     }
    2189
    2190                     t -= s;
    2191                     c = coeffs[which][s];
    2192                     if (c === undefined) {
    2193                         return NaN;
    2194                     }
    2195
    2196                     return (((c[3] * t + c[2]) * t + c[1]) * t + c[0]);
    2197                 };
    2198             };
    2199
    2200             return [makeFct('X'), makeFct('Y'), 0,
    2201                 function () {
    2202                     return points.length - 1;
    2203                 }];
    2204         },
    2205
    2206         /**
    2207          * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve
    2208          * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5.
    2209          * Two artificial control points at the beginning and the end are added.
    2210          * @param {Array} points Array consisting of JXG.Points.
    2211          * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and
    2212          * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal".
    2213          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
    2214          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
    2215          * returning the length of the points array minus three.
    2216          * @memberof JXG.Math.Numerics
    2217         */
    2218         CatmullRomSpline: function (points, type) {
    2219             return this.CardinalSpline(points, 0.5, type);
    2220         },
    2221
    2222         /**
    2223          * Computes the regression polynomial of a given degree through a given set of coordinates.
    2224          * Returns the regression polynomial function.
    2225          * @param {Number,function,Slider} degree number, function or slider.
    2226          * Either
    2227          * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in
    2228          * an array of {@link JXG.Point}s or {@link JXG.Coords}.
    2229          * In the latter case, the <tt>dataY</tt> parameter will be ignored.
    2230          * @param {Array} dataY Array containing the y-coordinates of the data set,
    2231          * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree.
    2232          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
    2233          * @memberof JXG.Math.Numerics
    2234          */
    2235         regressionPolynomial: function (degree, dataX, dataY) {
    2236             var coeffs, deg, dX, dY, inputType, fct,
    2237                 term = '';
    2238
    2239             // Slider
    2240             if (Type.isPoint(degree) && Type.isFunction(degree.Value)) {
    2241                 /** @ignore */
    2242                 deg = function () {
    2243                     return degree.Value();
    2244                 };
    2245             // function
    2246             } else if (Type.isFunction(degree)) {
    2247                 deg = degree;
    2248             // number
    2249             } else if (Type.isNumber(degree)) {
    2250                 /** @ignore */
    2251                 deg = function () {
    2252                     return degree;
    2253                 };
    2254             } else {
    2255                 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'.");
    2256             }
    2257
    2258             // Parameters degree, dataX, dataY
    2259             if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) {
    2260                 inputType = 0;
    2261             // Parameters degree, point array
    2262             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) {
    2263                 inputType = 1;
    2264             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) {
    2265                 inputType = 2;
    2266             } else {
    2267                 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters.");
    2268             }
    2269
    2270             /** @ignore */
    2271             fct = function (x, suspendedUpdate) {
    2272                 var i, j, M, MT, y, B, c, s, d,
    2273                     // input data
    2274                     len = dataX.length;
    2275
    2276                 d = Math.floor(deg());
    2277
    2278                 if (!suspendedUpdate) {
    2279                     // point list as input
    2280                     if (inputType === 1) {
    2281                         dX = [];
    2282                         dY = [];
    2283
    2284                         for (i = 0; i < len; i++) {
    2285                             dX[i] = dataX[i].X();
    2286                             dY[i] = dataX[i].Y();
    2287                         }
    2288                     }
    2289
    2290                     if (inputType === 2) {
    2291                         dX = [];
    2292                         dY = [];
    2293
    2294                         for (i = 0; i < len; i++) {
    2295                             dX[i] = dataX[i].usrCoords[1];
    2296                             dY[i] = dataX[i].usrCoords[2];
    2297                         }
    2298                     }
    2299
    2300                     // check for functions
    2301                     if (inputType === 0) {
    2302                         dX = [];
    2303                         dY = [];
    2304
    2305                         for (i = 0; i < len; i++) {
    2306                             if (Type.isFunction(dataX[i])) {
    2307                                 dX.push(dataX[i]());
    2308                             } else {
    2309                                 dX.push(dataX[i]);
    2310                             }
    2311
    2312                             if (Type.isFunction(dataY[i])) {
    2313                                 dY.push(dataY[i]());
    2314                             } else {
    2315                                 dY.push(dataY[i]);
    2316                             }
    2317                         }
    2318                     }
    2319
    2320                     M = [];
    2321
    2322                     for (j = 0; j < len; j++) {
    2323                         M.push([1]);
    2324                     }
    2325
    2326                     for (i = 1; i <= d; i++) {
    2327                         for (j = 0; j < len; j++) {
    2328                             M[j][i] = M[j][i - 1] * dX[j];
    2329                         }
    2330                     }
    2331
    2332                     y = dY;
    2333                     MT = Mat.transpose(M);
    2334                     B = Mat.matMatMult(MT, M);
    2335                     c = Mat.matVecMult(MT, y);
    2336                     coeffs = Mat.Numerics.Gauss(B, c);
    2337                     term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3);
    2338                 }
    2339
    2340                 // Horner's scheme to evaluate polynomial
    2341                 s = coeffs[d];
    2342
    2343                 for (i = d - 1; i >= 0; i--) {
    2344                     s = (s * x + coeffs[i]);
    2345                 }
    2346
    2347                 return s;
    2348             };
    2349
    2350             fct.getTerm = function () {
    2351                 return term;
    2352             };
    2353
    2354             return fct;
    2355         },
    2356
    2357         /**
    2358          * Computes the cubic Bezier curve through a given set of points.
    2359          * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}.
    2360          * The points at position k with k mod 3 = 0 are the data points,
    2361          * points at position k with k mod 3 = 1 or 2 are the control points.
    2362          * @returns {Array} An array consisting of two functions of one parameter t which return the
    2363          * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting
    2364          * no parameters and returning one third of the length of the points.
    2365          * @memberof JXG.Math.Numerics
    2366          */
    2367         bezier: function (points) {
    2368             var len, flen,
    2369                 /** @ignore */
    2370                 makeFct = function (which) {
    2371                     return function (t, suspendedUpdate) {
    2372                         var z = Math.floor(t) * 3,
    2373                             t0 = t % 1,
    2374                             t1 = 1 - t0;
    2375
    2376                         if (!suspendedUpdate) {
    2377                             flen = 3 * Math.floor((points.length - 1) / 3);
    2378                             len = Math.floor(flen / 3);
    2379                         }
    2380
    2381                         if (t < 0) {
    2382                             return points[0][which]();
    2383                         }
    2384
    2385                         if (t >= len) {
    2386                             return points[flen][which]();
    2387                         }
    2388
    2389                         if (isNaN(t)) {
    2390                             return NaN;
    2391                         }
    2392
    2393                         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;
    2394                     };
    2395                 };
    2396
    2397             return [makeFct('X'), makeFct('Y'), 0,
    2398                 function () {
    2399                     return Math.floor(points.length / 3);
    2400                 }];
    2401         },
    2402
    2403         /**
    2404          * Computes the B-spline curve of order k (order = degree+1) through a given set of points.
    2405          * @param {Array} points Array consisting of JXG.Points.
    2406          * @param {Number} order Order of the B-spline curve.
    2407          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
    2408          * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply
    2409          * returning the length of the points array minus one.
    2410          * @memberof JXG.Math.Numerics
    2411          */
    2412         bspline: function (points, order) {
    2413             var knots,
    2414                 _knotVector = function (n, k) {
    2415                     var j,
    2416                         kn = [];
    2417
    2418                     for (j = 0; j < n + k + 1; j++) {
    2419                         if (j < k) {
    2420                             kn[j] = 0.0;
    2421                         } else if (j <= n) {
    2422                             kn[j] = j - k + 1;
    2423                         } else {
    2424                             kn[j] = n - k + 2;
    2425                         }
    2426                     }
    2427
    2428                     return kn;
    2429                 },
    2430
    2431                 _evalBasisFuncs = function (t, kn, k, s) {
    2432                     var i, j, a, b, den,
    2433                         N = [];
    2434
    2435                     if (kn[s] <= t && t < kn[s + 1]) {
    2436                         N[s] = 1;
    2437                     } else {
    2438                         N[s] = 0;
    2439                     }
    2440
    2441                     for (i = 2; i <= k; i++) {
    2442                         for (j = s - i + 1; j <= s; j++) {
    2443                             if (j <= s - i + 1 || j < 0) {
    2444                                 a = 0.0;
    2445                             } else {
    2446                                 a = N[j];
    2447                             }
    2448
    2449                             if (j >= s) {
    2450                                 b = 0.0;
    2451                             } else {
    2452                                 b = N[j + 1];
    2453                             }
    2454
    2455                             den = kn[j + i - 1] - kn[j];
    2456
    2457                             if (den === 0) {
    2458                                 N[j] = 0;
    2459                             } else {
    2460                                 N[j] = (t - kn[j]) / den * a;
    2461                             }
    2462
    2463                             den = kn[j + i] - kn[j + 1];
    2464
    2465                             if (den !== 0) {
    2466                                 N[j] += (kn[j + i] - t) / den * b;
    2467                             }
    2468                         }
    2469                     }
    2470                     return N;
    2471                 },
    2472                 /** @ignore */
    2473                 makeFct = function (which) {
    2474                     return function (t, suspendedUpdate) {
    2475                         var y, j, s, N = [],
    2476                             len = points.length,
    2477                             n = len - 1,
    2478                             k = order;
    2479
    2480                         if (n <= 0) {
    2481                             return NaN;
    2482                         }
    2483
    2484                         if (n + 2 <= k) {
    2485                             k = n + 1;
    2486                         }
    2487
    2488                         if (t <= 0) {
    2489                             return points[0][which]();
    2490                         }
    2491
    2492                         if (t >= n - k + 2) {
    2493                             return points[n][which]();
    2494                         }
    2495
    2496                         s = Math.floor(t) + k - 1;
    2497                         knots = _knotVector(n, k);
    2498                         N = _evalBasisFuncs(t, knots, k, s);
    2499
    2500                         y = 0.0;
    2501                         for (j = s - k + 1; j <= s; j++) {
    2502                             if (j < len && j >= 0) {
    2503                                 y += points[j][which]() * N[j];
    2504                             }
    2505                         }
    2506
    2507                         return y;
    2508                     };
    2509                 };
    2510
    2511             return [makeFct('X'), makeFct('Y'), 0,
    2512                 function () {
    2513                     return points.length - 1;
    2514                 }];
    2515         },
    2516
    2517         /**
    2518          * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through,
    2519          * see {@link JXG.Curve#updateCurve}
    2520          * and {@link JXG.Curve#hasPoint}.
    2521          * @param {function} f Function in one variable to be differentiated.
    2522          * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a
    2523          * method of an object and contains a reference to its parent object via "this".
    2524          * @returns {function} Derivative function of a given function f.
    2525          * @memberof JXG.Math.Numerics
    2526          */
    2527         D: function (f, obj) {
    2528             if (!Type.exists(obj)) {
    2529                 return function (x, suspendedUpdate) {
    2530                     var h = 0.00001,
    2531                         h2 = (h * 2.0);
    2532
    2533                     // Experiments with Richardsons rule
    2534                     /*
    2535                     var phi = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
    2536                     var phi2;
    2537                     h *= 0.5;
    2538                     h2 *= 0.5;
    2539                     phi2 = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
    2540
    2541                     return phi2 + (phi2 - phi) / 3.0;
    2542                     */
    2543                     return (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
    2544                 };
    2545             }
    2546
    2547             return function (x, suspendedUpdate) {
    2548                 var h = 0.00001,
    2549                     h2 = (h * 2.0);
    2550
    2551                 return (f.apply(obj, [x + h, suspendedUpdate]) - f.apply(obj, [x - h, suspendedUpdate])) / h2;
    2552             };
    2553         },
    2554
    2555         /**
    2556          * Evaluate the function term for {@see #riemann}.
    2557          * @private
    2558          * @param {Number} x function argument
    2559          * @param {function} f JavaScript function returning a number
    2560          * @param {String} type Name of the Riemann sum type, e.g. 'lower', see {@see #riemann}.
    2561          * @param {Number} delta Width of the bars in user coordinates
    2562          * @returns {Number} Upper (delta > 0) or lower (delta < 0) value of the bar containing x of the Riemann sum.
    2563          *
    2564          * @memberof JXG.Math.Numerics
    2565          */
    2566         _riemannValue: function (x, f, type, delta) {
    2567             var y, y1, x1, delta1;
    2568
    2569             if (delta < 0) { // delta is negative if the lower function term is evaluated
    2570                 if (type !== 'trapezoidal') {
    2571                     x = x + delta;
    2572                 }
    2573                 delta *= -1;
    2574                 if (type === 'lower') {
    2575                     type = 'upper';
    2576                 } else if (type === 'upper') {
    2577                     type = 'lower';
    2578                 }
    2579             }
    2580
    2581             delta1 = delta * 0.01; // for 'lower' and 'upper'
    2582
    2583             if (type === 'right') {
    2584                 y = f(x + delta);
    2585             } else if (type === 'middle') {
    2586                 y = f(x + delta * 0.5);
    2587             } else if (type === 'left' || type === 'trapezoidal') {
    2588                 y = f(x);
    2589             } else if (type === 'lower') {
    2590                 y = f(x);
    2591
    2592                 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
    2593                     y1 = f(x1);
    2594
    2595                     if (y1 < y) {
    2596                         y = y1;
    2597                     }
    2598                 }
    2599
    2600                 y1 = f(x + delta);
    2601                 if (y1 < y) {
    2602                     y = y1;
    2603                 }
    2604             } else if (type === 'upper') {
    2605                 y = f(x);
    2606
    2607                 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
    2608                     y1 = f(x1);
    2609                     if (y1 > y) {
    2610                         y = y1;
    2611                     }
    2612                 }
    2613
    2614                 y1 = f(x + delta);
    2615                 if (y1 > y) {
    2616                     y = y1;
    2617                 }
    2618             } else if (type === 'random') {
    2619                 y = f(x + delta * Math.random());
    2620             } else if (type === 'simpson') {
    2621                 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0;
    2622             } else {
    2623                 y = f(x);  // default is lower
    2624             }
    2625
    2626             return y;
    2627         },
    2628
    2629         /**
    2630          * Helper function to create curve which displays Riemann sums.
    2631          * Compute coordinates for the rectangles showing the Riemann sum.
    2632          * @param {Function,Array} f Function or array of two functions.
    2633          * If f is a function the integral of this function is approximated by the Riemann sum.
    2634          * If f is an array consisting of two functions the area between the two functions is filled
    2635          * by the Riemann sum bars.
    2636          * @param {Number} n number of rectangles.
    2637          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'.
    2638          * @param {Number} start Left border of the approximation interval
    2639          * @param {Number} end Right border of the approximation interval
    2640          * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This
    2641          * 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
    2642          * rectangles.
    2643          * @memberof JXG.Math.Numerics
    2644          */
    2645         riemann: function (gf, n, type, start, end) {
    2646             var i, delta,
    2647                 xarr = [],
    2648                 yarr = [],
    2649                 j = 0,
    2650                 x = start, y,
    2651                 sum = 0,
    2652                 f, g,
    2653                 ylow, yup;
    2654
    2655             if (Type.isArray(gf)) {
    2656                 g = gf[0];
    2657                 f = gf[1];
    2658             } else {
    2659                 f = gf;
    2660             }
    2661
    2662             n = Math.floor(n);
    2663
    2664             if (n <= 0) {
    2665                 return [xarr, yarr, sum];
    2666             }
    2667
    2668             delta = (end - start) / n;
    2669
    2670             // Upper bar ends
    2671             for (i = 0; i < n; i++) {
    2672                 y = this._riemannValue(x, f, type, delta);
    2673                 xarr[j] = x;
    2674                 yarr[j] = y;
    2675
    2676                 j += 1;
    2677                 x += delta;
    2678                 if (type === 'trapezoidal') {
    2679                     y = f(x);
    2680                 }
    2681                 xarr[j] = x;
    2682                 yarr[j] = y;
    2683
    2684                 j += 1;
    2685             }
    2686
    2687             // Lower bar ends
    2688             for (i = 0; i < n; i++) {
    2689                 if (g) {
    2690                     y = this._riemannValue(x, g, type, -delta);
    2691                 } else {
    2692                     y = 0.0;
    2693                 }
    2694                 xarr[j] = x;
    2695                 yarr[j] = y;
    2696
    2697                 j += 1;
    2698                 x -= delta;
    2699                 if (type === 'trapezoidal' && g) {
    2700                     y = g(x);
    2701                 }
    2702                 xarr[j] = x;
    2703                 yarr[j] = y;
    2704
    2705                 // Add the area of the bar to 'sum'
    2706                 if (type !== 'trapezoidal') {
    2707                     ylow = y;
    2708                     yup = yarr[2 * (n - 1) - 2 * i];
    2709                 } else {
    2710                     yup = 0.5 * (f(x + delta) + f(x));
    2711                     if (g) {
    2712                         ylow = 0.5 * (g(x + delta) + g(x));
    2713                     } else {
    2714                         ylow = 0.0;
    2715                     }
    2716                 }
    2717                 sum += (yup - ylow) * delta;
    2718
    2719                 // Draw the vertical lines
    2720                 j += 1;
    2721                 xarr[j] = x;
    2722                 yarr[j] = yarr[2 * (n - 1) - 2 * i];
    2723
    2724                 j += 1;
    2725             }
    2726
    2727             return [xarr, yarr, sum];
    2728         },
    2729
    2730         /**
    2731          * Approximate the integral by Riemann sums.
    2732          * Compute the area described by the riemann sum rectangles.
    2733          *
    2734          * If there is an element of type {@link Riemannsum}, then it is more efficient
    2735          * to use the method JXG.Curve.Value() of this element instead.
    2736          *
    2737          * @param {Function_Array} f Function or array of two functions.
    2738          * If f is a function the integral of this function is approximated by the Riemann sum.
    2739          * If f is an array consisting of two functions the area between the two functions is approximated
    2740          * by the Riemann sum.
    2741          * @param {Number} n number of rectangles.
    2742          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'.
    2743          *
    2744          * @param {Number} start Left border of the approximation interval
    2745          * @param {Number} end Right border of the approximation interval
    2746          * @returns {Number} The sum of the areas of the rectangles.
    2747          * @memberof JXG.Math.Numerics
    2748          */
    2749         riemannsum: function (f, n, type, start, end) {
    2750             JXG.deprecated('Numerics.riemannsum()', 'Numerics.riemann()');
    2751             return this.riemann(f, n, type, start, end)[2];
    2752         },
    2753
    2754         /**
    2755          * Solve initial value problems numerically using Runge-Kutta-methods.
    2756          * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
    2757          * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
    2758          * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
    2759          * <pre>
    2760          * {
    2761          *     s: <Number>,
    2762          *     A: <matrix>,
    2763          *     b: <Array>,
    2764          *     c: <Array>
    2765          * }
    2766          * </pre>
    2767          * 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
    2768          * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array.
    2769          * @param {Array} I Interval on which to integrate.
    2770          * @param {Number} N Number of evaluation points.
    2771          * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
    2772          * 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
    2773          * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has.
    2774          * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
    2775          * @example
    2776          * // A very simple autonomous system dx(t)/dt = x(t);
    2777          * function f(t, x) {
    2778          *     return x;
    2779          * }
    2780          *
    2781          * // Solve it with initial value x(0) = 1 on the interval [0, 2]
    2782          * // with 20 evaluation points.
    2783          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
    2784          *
    2785          * // Prepare data for plotting the solution of the ode using a curve.
    2786          * var dataX = [];
    2787          * var dataY = [];
    2788          * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
    2789          * for(var i=0; i<data.length; i++) {
    2790          *     dataX[i] = i*h;
    2791          *     dataY[i] = data[i][0];
    2792          * }
    2793          * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
    2794          * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
    2795          * <script type="text/javascript">
    2796          * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
    2797          * function f(t, x) {
    2798          *     // we have to copy the value.
    2799          *     // return x; would just return the reference.
    2800          *     return [x[0]];
    2801          * }
    2802          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
    2803          * var dataX = [];
    2804          * var dataY = [];
    2805          * var h = 0.1;
    2806          * for(var i=0; i<data.length; i++) {
    2807          *     dataX[i] = i*h;
    2808          *     dataY[i] = data[i][0];
    2809          * }
    2810          * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
    2811          * </script><pre>
    2812          * @memberof JXG.Math.Numerics
    2813          */
    2814         rungeKutta: function (butcher, x0, I, N, f) {
    2815             var e, i, j, k, l, s,
    2816                 x = [],
    2817                 y = [],
    2818                 h = (I[1] - I[0]) / N,
    2819                 t = I[0],
    2820                 dim = x0.length,
    2821                 result = [],
    2822                 r = 0;
    2823
    2824             if (Type.isString(butcher)) {
    2825                 butcher = predefinedButcher[butcher] || predefinedButcher.euler;
    2826             }
    2827             s = butcher.s;
    2828
    2829             // don't change x0, so copy it
    2830             for (e = 0; e < dim; e++) {
    2831                 x[e] = x0[e];
    2832             }
    2833
    2834             for (i = 0; i < N; i++) {
    2835                 // Optimization doesn't work for ODEs plotted using time
    2836                 //        if((i % quotient == 0) || (i == N-1)) {
    2837                 result[r] = [];
    2838                 for (e = 0; e < dim; e++) {
    2839                     result[r][e] = x[e];
    2840                 }
    2841
    2842                 r += 1;
    2843                 k = [];
    2844
    2845                 for (j = 0; j < s; j++) {
    2846                     // init y = 0
    2847                     for (e = 0; e < dim; e++) {
    2848                         y[e] = 0.0;
    2849                     }
    2850
    2851
    2852                     // Calculate linear combination of former k's and save it in y
    2853                     for (l = 0; l < j; l++) {
    2854                         for (e = 0; e < dim; e++) {
    2855                             y[e] += (butcher.A[j][l]) * h * k[l][e];
    2856                         }
    2857                     }
    2858
    2859                     // add x(t) to y
    2860                     for (e = 0; e < dim; e++) {
    2861                         y[e] += x[e];
    2862                     }
    2863
    2864                     // calculate new k and add it to the k matrix
    2865                     k.push(f(t + butcher.c[j] * h, y));
    2866                 }
    2867
    2868                 // init y = 0
    2869                 for (e = 0; e < dim; e++) {
    2870                     y[e] = 0.0;
    2871                 }
    2872
    2873                 for (l = 0; l < s; l++) {
    2874                     for (e = 0; e < dim; e++) {
    2875                         y[e] += butcher.b[l] * k[l][e];
    2876                     }
    2877                 }
    2878
    2879                 for (e = 0; e < dim; e++) {
    2880                     x[e] = x[e] + h * y[e];
    2881                 }
    2882
    2883                 t += h;
    2884             }
    2885
    2886             return result;
    2887         },
    2888
    2889         /**
    2890          * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and
    2891          * {@link JXG.Math.Numerics.chandrupatla}
    2892          * @type Number
    2893          * @default 80
    2894          * @memberof JXG.Math.Numerics
    2895          */
    2896         maxIterationsRoot: 80,
    2897
    2898         /**
    2899          * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr}
    2900          * @type Number
    2901          * @default 500
    2902          * @memberof JXG.Math.Numerics
    2903          */
    2904         maxIterationsMinimize: 500,
    2905
    2906         /**
    2907          * Given a value x_0, this function tries to find a second value x_1 such that
    2908          * the function f has opposite signs at x_0 and x_1.
    2909          * The return values have to be tested if the method succeeded.
    2910          *
    2911          * @param {Function} f Function, whose root is to be found
    2912          * @param {Number} x0 Start value
    2913          * @param {Object} object Parent object in case f is method of it
    2914          * @returns {Array} [x_0, f(x_0), x_1, f(x_1)]
    2915          *
    2916          * @see JXG.Math.Numerics.fzero
    2917          * @see JXG.Math.Numerics.chandrupatla
    2918          *
    2919          * @memberof JXG.Math.Numerics
    2920          */
    2921         findBracket: function(f, x0, object) {
    2922             var a, aa, fa,
    2923                 blist, b, fb,
    2924                 u, fu,
    2925                 i, len;
    2926
    2927             if (Type.isArray(x0)) {
    2928                 return x0;
    2929             }
    2930
    2931             a = x0;
    2932             fa = f.call(object, a);
    2933             // nfev += 1;
    2934
    2935             // Try to get b, by trying several values related to a
    2936             aa = (a === 0) ? 1 : a;
    2937             blist = [
    2938                 a - 0.1 * aa, a + 0.1 * aa,
    2939                 a - 1, a + 1,
    2940                 a - 0.5 * aa, a + 0.5 * aa,
    2941                 a - 0.6 * aa, a + 0.6 * aa,
    2942                 a - 1 * aa, a + 1 * aa,
    2943                 a - 2 * aa, a + 2 * aa,
    2944                 a - 5 * aa, a + 5 * aa,
    2945                 a - 10 * aa, a + 10 * aa,
    2946                 a - 50 * aa, a + 50 * aa,
    2947                 a - 100 * aa, a + 100 * aa
    2948             ];
    2949             len = blist.length;
    2950
    2951             for (i = 0; i < len; i++) {
    2952                 b = blist[i];
    2953                 fb = f.call(object, b);
    2954                 // nfev += 1;
    2955
    2956                 if (fa * fb <= 0) {
    2957                     break;
    2958                 }
    2959             }
    2960             if (b < a) {
    2961                 u = a;
    2962                 a = b;
    2963                 b = u;
    2964
    2965                 fu = fa;
    2966                 fa = fb;
    2967                 fb = fu;
    2968             }
    2969             return [a, fa, b, fb];
    2970         },
    2971
    2972         /**
    2973          *
    2974          * Find zero of an univariate function f.
    2975          * @param {function} f Function, whose root is to be found
    2976          * @param {Array,Number} x0  Start value or start interval enclosing the root
    2977          * @param {Object} object Parent object in case f is method of it
    2978          * @returns {Number} the approximation of the root
    2979          * Algorithm:
    2980          *  Brent's root finder from
    2981          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
    2982          *  computations. M., Mir, 1980, p.180 of the Russian edition
    2983          *  http://www.netlib.org/c/brent.shar
    2984          *
    2985          * If x0 is an array containing lower and upper bound for the zero
    2986          * algorithm 748 is applied. Otherwise, if x0 is a number,
    2987          * the algorithm tries to bracket a zero of f starting from x0.
    2988          * If this fails, we fall back to Newton's method.
    2989          *
    2990          * @see JXG.Math.Numerics.chandrupatla
    2991          * @see JXG.Math.Numerics.root
    2992          * @memberof JXG.Math.Numerics
    2993          */
    2994         fzero: function (f, x0, object) {
    2995             var a, b, c,
    2996                 d, e,
    2997                 fa, fb, fc,
    2998                 res,
    2999                 prev_step, t1, cb, t2,
    3000                 // Actual tolerance
    3001                 tol_act,
    3002                 // Interpolation step is calculated in the form p/q; division
    3003                 // operations is delayed until the last moment
    3004                 p, q,
    3005                 // Step at this iteration
    3006                 new_step,
    3007                 eps = Mat.eps,
    3008                 maxiter = this.maxIterationsRoot,
    3009                 niter = 0;
    3010                 // nfev = 0;
    3011
    3012             if (Type.isArray(x0)) {
    3013                 if (x0.length < 2) {
    3014                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
    3015                 }
    3016
    3017                 a = x0[0];
    3018                 fa = f.call(object, a);
    3019                 // nfev += 1;
    3020                 b = x0[1];
    3021                 fb = f.call(object, b);
    3022                 // nfev += 1;
    3023             } else {
    3024                 res = this.findBracket(f, x0, object);
    3025                 a = res[0];
    3026                 fa = res[1];
    3027                 b = res[2];
    3028                 fb = res[3];
    3029             }
    3030
    3031             if (Math.abs(fa) <= eps) {
    3032                 return a;
    3033             }
    3034             if (Math.abs(fb) <= eps) {
    3035                 return b;
    3036             }
    3037
    3038             if (fa * fb > 0) {
    3039                 // Bracketing not successful, fall back to Newton's method or to fminbr
    3040                 if (Type.isArray(x0)) {
    3041                     return this.fminbr(f, [a, b], object);
    3042                 }
    3043
    3044                 return this.Newton(f, a, object);
    3045             }
    3046
    3047             // OK, we have enclosed a zero of f.
    3048             // Now we can start Brent's method
    3049             c = a;
    3050             fc = fa;
    3051
    3052             // Main iteration loop
    3053             while (niter < maxiter) {
    3054                 // Distance from the last but one to the last approximation
    3055                 prev_step = b - a;
    3056
    3057                 // Swap data for b to be the best approximation
    3058                 if (Math.abs(fc) < Math.abs(fb)) {
    3059                     a = b;
    3060                     b = c;
    3061                     c = a;
    3062
    3063                     fa = fb;
    3064                     fb = fc;
    3065                     fc = fa;
    3066                 }
    3067                 tol_act = 2 * eps * Math.abs(b) + eps * 0.5;
    3068                 new_step = (c - b) * 0.5;
    3069
    3070                 if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) {
    3071                     //  Acceptable approx. is found
    3072                     return b;
    3073                 }
    3074
    3075                 // Decide if the interpolation can be tried
    3076                 // If prev_step was large enough and was in true direction Interpolatiom may be tried
    3077                 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
    3078                     cb = c - b;
    3079
    3080                     // If we have only two distinct points linear interpolation can only be applied
    3081                     if (a === c) {
    3082                         t1 = fb / fa;
    3083                         p = cb * t1;
    3084                         q = 1.0 - t1;
    3085                     // Quadric inverse interpolation
    3086                     } else {
    3087                         q = fa / fc;
    3088                         t1 = fb / fc;
    3089                         t2 = fb / fa;
    3090
    3091                         p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
    3092                         q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
    3093                     }
    3094
    3095                     // p was calculated with the opposite sign; make p positive
    3096                     if (p > 0) {
    3097                         q = -q;
    3098                     // and assign possible minus to q
    3099                     } else {
    3100                         p = -p;
    3101                     }
    3102
    3103                     // If b+p/q falls in [b,c] and isn't too large it is accepted
    3104                     // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
    3105                     if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) &&
    3106                             p < Math.abs(prev_step * q * 0.5)) {
    3107                         new_step = p / q;
    3108                     }
    3109                 }
    3110
    3111                 // Adjust the step to be not less than tolerance
    3112                 if (Math.abs(new_step) < tol_act) {
    3113                     new_step = (new_step > 0) ? tol_act : -tol_act;
    3114                 }
    3115
    3116                 // Save the previous approx.
    3117                 a  = b;
    3118                 fa = fb;
    3119                 b += new_step;
    3120                 fb = f.call(object, b);
    3121                 // Do step to a new approxim.
    3122                 // nfev += 1;
    3123
    3124                 // Adjust c for it to have a sign opposite to that of b
    3125                 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
    3126                     c = a;
    3127                     fc = fa;
    3128                 }
    3129                 niter++;
    3130             } // End while
    3131
    3132             return b;
    3133         },
    3134
    3135         /**
    3136          * Find zero of an univariate function f.
    3137          * @param {function} f Function, whose root is to be found
    3138          * @param {Array,Number} x0  Start value or start interval enclosing the root
    3139          * @param {Object} object Parent object in case f is method of it
    3140          * @returns {Number} the approximation of the root
    3141          * Algorithm:
    3142          * Chandrupatla's method, see
    3143          * Tirupathi R. Chandrupatla,
    3144          * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives",
    3145          * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149.
    3146          *
    3147          * If x0 is an array containing lower and upper bound for the zero
    3148          * algorithm 748 is applied. Otherwise, if x0 is a number,
    3149          * the algorithm tries to bracket a zero of f starting from x0.
    3150          * If this fails, we fall back to Newton's method.
    3151          *
    3152          * @see JXG.Math.Numerics.root
    3153          * @see JXG.Math.Numerics.fzero
    3154          * @memberof JXG.Math.Numerics
    3155          */
    3156         chandrupatla: function (f, x0, object) {
    3157             var a, fa, b, fb, res,
    3158                 niter = 0,
    3159                 maxiter = this.maxIterationsRoot,
    3160                 rand = (1 + Math.random() * 0.001),
    3161                 t = 0.5 * rand,
    3162                 eps = Mat.eps, // 1.e-10,
    3163                 dlt = 0.00001,
    3164                 x1, x2, x3, x,
    3165                 f1, f2, f3, y,
    3166                 xm, fm, tol, tl,
    3167                 xi, ph, fl, fh,
    3168                 AL, A, B, C, D;
    3169
    3170             if (Type.isArray(x0)) {
    3171                 if (x0.length < 2) {
    3172                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
    3173                 }
    3174
    3175                 a = x0[0];
    3176                 fa = f.call(object, a);
    3177                 // nfev += 1;
    3178                 b = x0[1];
    3179                 fb = f.call(object, b);
    3180                 // nfev += 1;
    3181             } else {
    3182                 res = this.findBracket(f, x0, object);
    3183                 a = res[0];
    3184                 fa = res[1];
    3185                 b = res[2];
    3186                 fb = res[3];
    3187             }
    3188
    3189             if (fa * fb > 0) {
    3190                 // Bracketing not successful, fall back to Newton's method or to fminbr
    3191                 if (Type.isArray(x0)) {
    3192                     return this.fminbr(f, [a, b], object);
    3193                 }
    3194
    3195                 return this.Newton(f, a, object);
    3196             }
    3197
    3198             x1 = a;  x2 = b;
    3199             f1 = fa; f2 = fb;
    3200             do {
    3201                 x = x1 + t * (x2 - x1);
    3202                 y = f.call(object, x);
    3203
    3204                 // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point
    3205                 if (Math.sign(y) === Math.sign(f1)) {
    3206                     x3 = x1; x1 = x;
    3207                     f3 = f1; f1 = y;
    3208                 } else {
    3209                     x3 = x2; x2 = x1;
    3210                     f3 = f2; f2 = f1;
    3211                 }
    3212                 x1 = x; f1 = y;
    3213
    3214                 xm = x1; fm = f1;
    3215                 if (Math.abs(f2) < Math.abs(f1)) {
    3216                     xm = x2; fm = f2;
    3217                 }
    3218                 tol = 2 * eps * Math.abs(xm) + 0.5 * dlt;
    3219                 tl = tol / Math.abs(x2 - x1);
    3220                 if (tl > 0.5 || fm === 0) {
    3221                     break;
    3222                 }
    3223                 // If inverse quadratic interpolation holds, use it
    3224                 xi = (x1 - x2) / (x3 - x2);
    3225                 ph = (f1 - f2) / (f3 - f2);
    3226                 fl = 1 - Math.sqrt(1 - xi);
    3227                 fh = Math.sqrt(xi);
    3228                 if (fl < ph && ph < fh) {
    3229                     AL = (x3 - x1) / (x2 - x1);
    3230                     A = f1 / (f2 - f1);
    3231                     B = f3 / (f2 - f3);
    3232                     C = f1 / (f3 - f1);
    3233                     D = f2 / (f3 - f2);
    3234                     t = A * B + C * D * AL;
    3235                 } else {
    3236                     t = 0.5 * rand;
    3237                 }
    3238                 // Adjust t away from the interval boundary
    3239                 if (t < tl) {
    3240                     t = tl;
    3241                 }
    3242                 if (t > 1 - tl) {
    3243                     t = 1 - tl;
    3244                 }
    3245                 niter++;
    3246             } while (niter <= maxiter);
    3247             // console.log(niter);
    3248
    3249             return xm;
    3250         },
    3251
    3252         /**
    3253          *
    3254          * Find minimum of an univariate function f.
    3255          * <p>
    3256          * Algorithm:
    3257          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
    3258          *  computations. M., Mir, 1980, p.180 of the Russian edition
    3259          *
    3260          * @param {function} f Function, whose minimum is to be found
    3261          * @param {Array} x0  Start interval enclosing the minimum
    3262          * @param {Object} context Parent object in case f is method of it
    3263          * @returns {Number} the approximation of the minimum value position
    3264          * @memberof JXG.Math.Numerics
    3265          **/
    3266         fminbr: function (f, x0, context) {
    3267             var a, b, x, v, w,
    3268                 fx, fv, fw,
    3269                 range, middle_range, tol_act, new_step,
    3270                 p, q, t, ft,
    3271                 // Golden section ratio
    3272                 r = (3.0 - Math.sqrt(5.0)) * 0.5,
    3273                 tol = Mat.eps,
    3274                 sqrteps = Mat.eps, //Math.sqrt(Mat.eps),
    3275                 maxiter = this.maxIterationsMinimize,
    3276                 niter = 0;
    3277                 // nfev = 0;
    3278
    3279             if (!Type.isArray(x0) || x0.length < 2) {
    3280                 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two.");
    3281             }
    3282
    3283             a = x0[0];
    3284             b = x0[1];
    3285             v = a + r * (b - a);
    3286             fv = f.call(context, v);
    3287
    3288             // First step - always gold section
    3289             // nfev += 1;
    3290             x = v;
    3291             w = v;
    3292             fx = fv;
    3293             fw = fv;
    3294
    3295             while (niter < maxiter) {
    3296                 // Range over the interval in which we are looking for the minimum
    3297                 range = b - a;
    3298                 middle_range = (a + b) * 0.5;
    3299
    3300                 // Actual tolerance
    3301                 tol_act = sqrteps * Math.abs(x) + tol / 3.0;
    3302
    3303                 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) {
    3304                     // Acceptable approx. is found
    3305                     return x;
    3306                 }
    3307
    3308                 // Obtain the golden section step
    3309                 new_step = r * (x < middle_range ? b - x : a - x);
    3310
    3311                 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried
    3312                 if (Math.abs(x - w) >= tol_act) {
    3313                     // Interpolation step is calculated as p/q;
    3314                     // division operation is delayed until last moment
    3315                     t = (x - w) * (fx - fv);
    3316                     q = (x - v) * (fx - fw);
    3317                     p = (x - v) * q - (x - w) * t;
    3318                     q = 2 * (q - t);
    3319
    3320                     if (q > 0) {                        // q was calculated with the op-
    3321                         p = -p;                         // posite sign; make q positive
    3322                     } else {                            // and assign possible minus to
    3323                         q = -q;                         // p
    3324                     }
    3325                     if (Math.abs(p) < Math.abs(new_step * q) &&     // If x+p/q falls in [a,b]
    3326                             p > q * (a - x + 2 * tol_act) &&        //  not too close to a and
    3327                             p < q * (b - x - 2 * tol_act)) {        // b, and isn't too large
    3328                         new_step = p / q;                          // it is accepted
    3329                     }
    3330                     // If p/q is too large then the
    3331                     // golden section procedure can
    3332                     // reduce [a,b] range to more
    3333                     // extent
    3334                 }
    3335
    3336                 // Adjust the step to be not less than tolerance
    3337                 if (Math.abs(new_step) < tol_act) {
    3338                     if (new_step > 0) {
    3339                         new_step = tol_act;
    3340                     } else {
    3341                         new_step = -tol_act;
    3342                     }
    3343                 }
    3344
    3345                 // Obtain the next approximation to min
    3346                 // and reduce the enveloping range
    3347
    3348                 // Tentative point for the min
    3349                 t = x + new_step;
    3350                 ft = f.call(context, t);
    3351                 // nfev += 1;
    3352
    3353                 // t is a better approximation
    3354                 if (ft <= fx) {
    3355                     // Reduce the range so that t would fall within it
    3356                     if (t < x) {
    3357                         b = x;
    3358                     } else {
    3359                         a = x;
    3360                     }
    3361
    3362                     // Assign the best approx to x
    3363                     v = w;
    3364                     w = x;
    3365                     x = t;
    3366
    3367                     fv = fw;
    3368                     fw = fx;
    3369                     fx = ft;
    3370                 // x remains the better approx
    3371                 } else {
    3372                     // Reduce the range enclosing x
    3373                     if (t < x) {
    3374                         a = t;
    3375                     } else {
    3376                         b = t;
    3377                     }
    3378
    3379                     if (ft <= fw || w === x) {
    3380                         v = w;
    3381                         w = t;
    3382                         fv = fw;
    3383                         fw = ft;
    3384                     } else if (ft <= fv || v === x || v === w) {
    3385                         v = t;
    3386                         fv = ft;
    3387                     }
    3388                 }
    3389                 niter += 1;
    3390             }
    3391
    3392             return x;
    3393         },
    3394
    3395         /**
    3396          * Implements the Ramer-Douglas-Peucker algorithm.
    3397          * It discards points which are not necessary from the polygonal line defined by the point array
    3398          * pts. The computation is done in screen coordinates.
    3399          * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points.
    3400          * @param {Array} pts Array of {@link JXG.Coords}
    3401          * @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>.
    3402          * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points.
    3403          * @memberof JXG.Math.Numerics
    3404          */
    3405         RamerDouglasPeucker: function (pts, eps) {
    3406             var allPts = [], newPts = [], i, k, len,
    3407
    3408                 /**
    3409                  * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
    3410                  * It searches for the point between index i and j which
    3411                  * has the largest distance from the line between the points i and j.
    3412                  * @param {Array} pts Array of {@link JXG.Coords}
    3413                  * @param {Number} i Index of a point in pts
    3414                  * @param {Number} j Index of a point in pts
    3415                  * @ignore
    3416                  * @private
    3417                  */
    3418                 findSplit = function (pts, i, j) {
    3419                     var d, k, ci, cj, ck,
    3420                         x0, y0, x1, y1,
    3421                         den, lbda,
    3422                         huge = 10000,
    3423                         dist = 0,
    3424                         f = i;
    3425
    3426                     if (j - i < 2) {
    3427                         return [-1.0, 0];
    3428                     }
    3429
    3430                     ci = pts[i].scrCoords;
    3431                     cj = pts[j].scrCoords;
    3432
    3433                     if (isNaN(ci[1]) || isNaN(ci[2])) {
    3434                         return [NaN, i];
    3435                     }
    3436                     if (isNaN(cj[1]) || isNaN(cj[2])) {
    3437                         return [NaN, j];
    3438                     }
    3439
    3440                     for (k = i + 1; k < j; k++) {
    3441                         ck = pts[k].scrCoords;
    3442                         if (isNaN(ck[1]) || isNaN(ck[2])) {
    3443                             return [NaN, k];
    3444                         }
    3445
    3446                         x0 = ck[1] - ci[1];
    3447                         y0 = ck[2] - ci[2];
    3448                         x1 = cj[1] - ci[1];
    3449                         y1 = cj[2] - ci[2];
    3450                         x0 = x0 === Infinity ? huge : x0;
    3451                         y0 = y0 === Infinity ? huge : y0;
    3452                         x1 = x1 === Infinity ? huge : x1;
    3453                         y1 = y1 === Infinity ? huge : y1;
    3454                         x0 = x0 === -Infinity ? -huge : x0;
    3455                         y0 = y0 === -Infinity ? -huge : y0;
    3456                         x1 = x1 === -Infinity ? -huge : x1;
    3457                         y1 = y1 === -Infinity ? -huge : y1;
    3458                         den = x1 * x1 + y1 * y1;
    3459
    3460                         if (den >= Mat.eps) {
    3461                             lbda = (x0 * x1 + y0 * y1) / den;
    3462
    3463                             if (lbda < 0.0) {
    3464                                 lbda = 0.0;
    3465                             } else if (lbda > 1.0) {
    3466                                 lbda = 1.0;
    3467                             }
    3468
    3469                             x0 = x0 - lbda * x1;
    3470                             y0 = y0 - lbda * y1;
    3471                             d = x0 * x0 + y0 * y0;
    3472                         } else {
    3473                             lbda = 0.0;
    3474                             d = x0 * x0 + y0 * y0;
    3475                         }
    3476
    3477                         if (d > dist) {
    3478                             dist = d;
    3479                             f = k;
    3480                         }
    3481                     }
    3482                     return [Math.sqrt(dist), f];
    3483                 },
    3484
    3485                 /**
    3486                  * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
    3487                  * It runs recursively through the point set and searches the
    3488                  * point which has the largest distance from the line between the first point and
    3489                  * the last point. If the distance from the line is greater than eps, this point is
    3490                  * included in our new point set otherwise it is discarded.
    3491                  * If it is taken, we recursively apply the subroutine to the point set before
    3492                  * and after the chosen point.
    3493                  * @param {Array} pts Array of {@link JXG.Coords}
    3494                  * @param {Number} i Index of an element of pts
    3495                  * @param {Number} j Index of an element of pts
    3496                  * @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>.
    3497                  * @param {Array} newPts Array of {@link JXG.Coords}
    3498                  * @ignore
    3499                  * @private
    3500                  */
    3501                 RDP = function (pts, i, j, eps, newPts) {
    3502                     var result = findSplit(pts, i, j),
    3503                         k = result[1];
    3504
    3505                     if (isNaN(result[0])) {
    3506                         RDP(pts, i, k - 1, eps, newPts);
    3507                         newPts.push(pts[k]);
    3508                         do {
    3509                             ++k;
    3510                         } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2]));
    3511                         if (k <= j) {
    3512                             newPts.push(pts[k]);
    3513                         }
    3514                         RDP(pts, k + 1, j, eps, newPts);
    3515                     } else if (result[0] > eps) {
    3516                         RDP(pts, i, k, eps, newPts);
    3517                         RDP(pts, k, j, eps, newPts);
    3518                     } else {
    3519                         newPts.push(pts[j]);
    3520                     }
    3521                 };
    3522
    3523             len = pts.length;
    3524
    3525             i = 0;
    3526             while (true) {
    3527                 // Search for the next point without NaN coordinates
    3528                 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) {
    3529                     i += 1;
    3530                 }
    3531                 // Search for the next position of a NaN point
    3532                 k = i + 1;
    3533                 while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) {
    3534                     k += 1;
    3535                 }
    3536                 k--;
    3537
    3538                 // Only proceed if something is left
    3539                 if (i < len && k > i) {
    3540                     newPts = [];
    3541                     newPts[0] = pts[i];
    3542                     RDP(pts, i, k, eps, newPts);
    3543                     allPts = allPts.concat(newPts);
    3544                 }
    3545                 if (i >= len) {
    3546                     break;
    3547                 }
    3548                 // Push the NaN point
    3549                 if (k < len - 1) {
    3550                     allPts.push(pts[k + 1]);
    3551                 }
    3552                 i = k + 1;
    3553             }
    3554
    3555             return allPts;
    3556         },
    3557
    3558         /**
    3559          * Old name for the implementation of the Ramer-Douglas-Peucker algorithm.
    3560          * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker}
    3561          * @memberof JXG.Math.Numerics
    3562          */
    3563         RamerDouglasPeuker: function (pts, eps) {
    3564             JXG.deprecated('Numerics.RamerDouglasPeuker()', 'Numerics.RamerDouglasPeucker()');
    3565             return this.RamerDouglasPeucker(pts, eps);
    3566         },
    3567
    3568         /**
    3569          * Implements the Visvalingam-Whyatt algorithm.
    3570          * See M. Visvalingam, J. D. Whyatt:
    3571          * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992
    3572          *
    3573          * The algorithm discards points which are not necessary from the polygonal line defined by the point array
    3574          * pts (consisting of type JXG.Coords).
    3575          * @param {Array} pts Array of {@link JXG.Coords}
    3576          * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will
    3577          *    be taken in any case.
    3578          * @returns {Array} An array containing points which approximates the curve defined by pts.
    3579          * @memberof JXG.Math.Numerics
    3580          *
    3581          * @example
    3582          *     var i, p = [];
    3583          *     for (i = 0; i < 5; ++i) {
    3584          *         p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
    3585          *     }
    3586          *
    3587          *     // Plot a cardinal spline curve
    3588          *     var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
    3589          *     var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
    3590          *
    3591          *     var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
    3592          *     c.updateDataArray = function() {
    3593          *         var i, len, points;
    3594          *
    3595          *         // Reduce number of intermediate points with Visvakingam-Whyatt to 6
    3596          *         points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
    3597          *         // Plot the remaining points
    3598          *         len = points.length;
    3599          *         this.dataX = [];
    3600          *         this.dataY = [];
    3601          *         for (i = 0; i < len; i++) {
    3602          *             this.dataX.push(points[i].usrCoords[1]);
    3603          *             this.dataY.push(points[i].usrCoords[2]);
    3604          *         }
    3605          *     };
    3606          *     board.update();
    3607          *
    3608          * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div>
    3609          * <script type="text/javascript">
    3610          *     (function() {
    3611          *         var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb',
    3612          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
    3613          *
    3614          *         var i, p = [];
    3615          *         for (i = 0; i < 5; ++i) {
    3616          *             p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
    3617          *         }
    3618          *
    3619          *         // Plot a cardinal spline curve
    3620          *         var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
    3621          *         var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
    3622          *
    3623          *         var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
    3624          *         c.updateDataArray = function() {
    3625          *             var i, len, points;
    3626          *
    3627          *             // Reduce number of intermediate points with Visvakingam-Whyatt to 6
    3628          *             points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
    3629          *             // Plot the remaining points
    3630          *             len = points.length;
    3631          *             this.dataX = [];
    3632          *             this.dataY = [];
    3633          *             for (i = 0; i < len; i++) {
    3634          *                 this.dataX.push(points[i].usrCoords[1]);
    3635          *                 this.dataY.push(points[i].usrCoords[2]);
    3636          *             }
    3637          *         };
    3638          *         board.update();
    3639          *
    3640          *     })();
    3641          *
    3642          * </script><pre>
    3643          *
    3644          */
    3645         Visvalingam: function(pts, numPoints) {
    3646             var i, len, vol, lastVol,
    3647                 linkedList = [],
    3648                 heap = [],
    3649                 points = [],
    3650                 lft, rt, lft2, rt2,
    3651                 obj;
    3652
    3653             len = pts.length;
    3654
    3655             // At least one intermediate point is needed
    3656             if (len <= 2) {
    3657                 return pts;
    3658             }
    3659
    3660             // Fill the linked list
    3661             // Add first point to the linked list
    3662             linkedList[0] = {
    3663                     used: true,
    3664                     lft: null,
    3665                     node: null
    3666                 };
    3667
    3668             // Add all intermediate points to the linked list,
    3669             // whose triangle area is nonzero.
    3670             lft = 0;
    3671             for (i = 1; i < len - 1; i++) {
    3672                 vol = Math.abs(JXG.Math.Numerics.det([pts[i - 1].usrCoords,
    3673                                               pts[i].usrCoords,
    3674                                               pts[i + 1].usrCoords]));
    3675                 if (!isNaN(vol)) {
    3676                     obj = {
    3677                         v: vol,
    3678                         idx: i
    3679                     };
    3680                     heap.push(obj);
    3681                     linkedList[i] = {
    3682                             used: true,
    3683                             lft: lft,
    3684                             node: obj
    3685                         };
    3686                     linkedList[lft].rt = i;
    3687                     lft = i;
    3688                 }
    3689             }
    3690
    3691             // Add last point to the linked list
    3692             linkedList[len - 1] = {
    3693                     used: true,
    3694                     rt: null,
    3695                     lft: lft,
    3696                     node: null
    3697                 };
    3698             linkedList[lft].rt = len - 1;
    3699
    3700             // Remove points until only numPoints intermediate points remain
    3701             lastVol = -Infinity;
    3702             while (heap.length > numPoints) {
    3703                 // Sort the heap with the updated volume values
    3704                 heap.sort(function(a, b) {
    3705                     // descending sort
    3706                     return b.v - a.v;
    3707                 });
    3708
    3709                 // Remove the point with the smallest triangle
    3710                 i = heap.pop().idx;
    3711                 linkedList[i].used = false;
    3712                 lastVol = linkedList[i].node.v;
    3713
    3714                 // Update the pointers of the linked list
    3715                 lft = linkedList[i].lft;
    3716                 rt = linkedList[i].rt;
    3717                 linkedList[lft].rt = rt;
    3718                 linkedList[rt].lft = lft;
    3719
    3720                 // Update the values for the volumes in the linked list
    3721                 lft2 = linkedList[lft].lft;
    3722                 if (lft2 !== null) {
    3723                     vol = Math.abs(JXG.Math.Numerics.det(
    3724                                 [pts[lft2].usrCoords,
    3725                                  pts[lft].usrCoords,
    3726                                  pts[rt].usrCoords]));
    3727
    3728                     linkedList[lft].node.v = (vol >= lastVol) ? vol : lastVol;
    3729                 }
    3730                 rt2 = linkedList[rt].rt;
    3731                 if (rt2 !== null) {
    3732                     vol = Math.abs(JXG.Math.Numerics.det(
    3733                                 [pts[lft].usrCoords,
    3734                                  pts[rt].usrCoords,
    3735                                  pts[rt2].usrCoords]));
    3736
    3737                     linkedList[rt].node.v = (vol >= lastVol) ? vol : lastVol;
    3738                 }
    3739             }
    3740
    3741             // Return an array with the remaining points
    3742             i = 0;
    3743             points = [pts[i]];
    3744             do {
    3745                 i = linkedList[i].rt;
    3746                 points.push(pts[i]);
    3747             } while (linkedList[i].rt !== null);
    3748
    3749             return points;
    3750         }
    3751     };
    3752
    3753     return Mat.Numerics;
    3754 });
    3755