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