1 /*
  2     Copyright 2008-2021
  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', 'utils/env', 'math/math'], function (JXG, Type, Env, 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 * 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,
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,
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, which - for the time being -
1424          * is an alias for {@link JXG.Math.Numerics.chandrupatla}.
1425          * @param {function} f We search for a solution of f(x)=0.
1426          * @param {Number|Array} x initial guess for the root, i.e. starting value, or start interval enclosing the root.
1427          * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
1428          * the function is a method of an object and contains a reference to its parent object via "this".
1429          * @returns {Number} A root of the function f.
1430          *
1431          * @see JXG.Math.Numerics.chandrupatla
1432          * @see JXG.Math.Numerics.fzero
1433          * @memberof JXG.Math.Numerics
1434          */
1435         root: function (f, x, context) {
1436             //return this.fzero(f, x, context);
1437             return this.chandrupatla(f, x, context);
1438         },
1439 
1440         /**
1441          * Compute an intersection of the curves c1 and c2
1442          * with a generalized Newton method.
1443          * We want to find values t1, t2 such that
1444          * c1(t1) = c2(t2), i.e.
1445          * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0).
1446          * We set
1447          * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2))
1448          *
1449          * The Jacobian J is defined by
1450          * J = (a, b)
1451          *     (c, d)
1452          * where
1453          * a = c1_x'(t1)
1454          * b = -c2_x'(t2)
1455          * c = c1_y'(t1)
1456          * d = -c2_y'(t2)
1457          *
1458          * The inverse J^(-1) of J is equal to
1459          *  (d, -b)/
1460          *  (-c, a) / (ad-bc)
1461          *
1462          * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f).
1463          * If the function meetCurveCurve possesses the properties
1464          * t1memo and t2memo then these are taken as start values
1465          * for the Newton algorithm.
1466          * After stopping of the Newton algorithm the values of t1 and t2 are stored in
1467          * t1memo and t2memo.
1468          *
1469          * @param {JXG.Curve} c1 Curve, Line or Circle
1470          * @param {JXG.Curve} c2 Curve, Line or Circle
1471          * @param {Number} t1ini start value for t1
1472          * @param {Number} t2ini start value for t2
1473          * @returns {JXG.Coords} intersection point
1474          * @memberof JXG.Math.Numerics
1475          */
1476         generalizedNewton: function (c1, c2, t1ini, t2ini) {
1477             var t1, t2,
1478                 a, b, c, d, disc,
1479                 e, f, F,
1480                 D00, D01,
1481                 D10, D11,
1482                 count = 0;
1483 
1484             if (this.generalizedNewton.t1memo) {
1485                 t1 = this.generalizedNewton.t1memo;
1486                 t2 = this.generalizedNewton.t2memo;
1487             } else {
1488                 t1 = t1ini;
1489                 t2 = t2ini;
1490             }
1491 
1492             e = c1.X(t1) - c2.X(t2);
1493             f = c1.Y(t1) - c2.Y(t2);
1494             F = e * e + f * f;
1495 
1496             D00 = this.D(c1.X, c1);
1497             D01 = this.D(c2.X, c2);
1498             D10 = this.D(c1.Y, c1);
1499             D11 = this.D(c2.Y, c2);
1500 
1501             while (F > Mat.eps && count < 10) {
1502                 a = D00(t1);
1503                 b = -D01(t2);
1504                 c = D10(t1);
1505                 d = -D11(t2);
1506                 disc = a * d - b * c;
1507                 t1 -= (d * e - b * f) / disc;
1508                 t2 -= (a * f - c * e) / disc;
1509                 e = c1.X(t1) - c2.X(t2);
1510                 f = c1.Y(t1) - c2.Y(t2);
1511                 F = e * e + f * f;
1512                 count += 1;
1513             }
1514 
1515             this.generalizedNewton.t1memo = t1;
1516             this.generalizedNewton.t2memo = t2;
1517 
1518             if (Math.abs(t1) < Math.abs(t2)) {
1519                 return [c1.X(t1), c1.Y(t1)];
1520             }
1521 
1522             return [c2.X(t2), c2.Y(t2)];
1523         },
1524 
1525         /**
1526          * Returns the Lagrange polynomials for curves with equidistant nodes, see
1527          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
1528          * SIAM Review, Vol 46, No 3, (2004) 501-517.
1529          * The graph of the parametric curve [x(t),y(t)] runs through the given points.
1530          * @param {Array} p Array of JXG.Points
1531          * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve
1532          * f(t) = (x(t), y(t)), a number x1 (which equals 0) and a function x2 defining the curve's domain.
1533          * That means the curve is defined between x1 and x2(). x2 returns the (length of array p minus one).
1534          * @memberof JXG.Math.Numerics
1535          *
1536          * @example
1537          * var p = [];
1538          *
1539          * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'});
1540          * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''});
1541          * p[2] = board.create('point', [1, 4], {size:2, name: ''});
1542          * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'});
1543          *
1544          * // Curve
1545          * var fg = JXG.Math.Numerics.Neville(p);
1546          * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5});
1547          *
1548          * </pre><div id="JXG88a8b3a8-6561-44f5-a678-76bca13fd484" class="jxgbox" style="width: 300px; height: 300px;"></div>
1549          * <script type="text/javascript">
1550          *     (function() {
1551          *         var board = JXG.JSXGraph.initBoard('JXG88a8b3a8-6561-44f5-a678-76bca13fd484',
1552          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
1553          *     var p = [];
1554          *
1555          *     p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'});
1556          *     p[1] = board.create('point', [-1.5, 5], {size:2, name: ''});
1557          *     p[2] = board.create('point', [1, 4], {size:2, name: ''});
1558          *     p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'});
1559          *
1560          *     // Curve
1561          *     var fg = JXG.Math.Numerics.Neville(p);
1562          *     var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5});
1563          *
1564          *     })();
1565          *
1566          * </script><pre>
1567          *
1568          */
1569         Neville: function (p) {
1570             var w = [],
1571                 /** @ignore */
1572                 makeFct = function (fun) {
1573                     return function (t, suspendedUpdate) {
1574                         var i, d, s,
1575                             bin = Mat.binomial,
1576                             len = p.length,
1577                             len1 = len - 1,
1578                             num = 0.0,
1579                             denom = 0.0;
1580 
1581                         if (!suspendedUpdate) {
1582                             s = 1;
1583                             for (i = 0; i < len; i++) {
1584                                 w[i] = bin(len1, i) * s;
1585                                 s *= (-1);
1586                             }
1587                         }
1588 
1589                         d = t;
1590 
1591                         for (i = 0; i < len; i++) {
1592                             if (d === 0) {
1593                                 return p[i][fun]();
1594                             }
1595                             s = w[i] / d;
1596                             d -= 1;
1597                             num += p[i][fun]() * s;
1598                             denom += s;
1599                         }
1600                         return num / denom;
1601                     };
1602                 },
1603 
1604                 xfct = makeFct('X'),
1605                 yfct = makeFct('Y');
1606 
1607             return [xfct, yfct, 0, function () {
1608                 return p.length - 1;
1609             }];
1610         },
1611 
1612         /**
1613          * Calculates second derivatives at the given knots.
1614          * @param {Array} x x values of knots
1615          * @param {Array} y y values of knots
1616          * @returns {Array} Second derivatives of the interpolated function at the knots.
1617          * @see #splineEval
1618          * @memberof JXG.Math.Numerics
1619          */
1620         splineDef: function (x, y) {
1621             var pair, i, l,
1622                 n = Math.min(x.length, y.length),
1623                 diag = [],
1624                 z = [],
1625                 data = [],
1626                 dx = [],
1627                 delta = [],
1628                 F = [];
1629 
1630             if (n === 2) {
1631                 return [0, 0];
1632             }
1633 
1634             for (i = 0; i < n; i++) {
1635                 pair = {X: x[i], Y: y[i]};
1636                 data.push(pair);
1637             }
1638             data.sort(function (a, b) {
1639                 return a.X - b.X;
1640             });
1641             for (i = 0; i < n; i++) {
1642                 x[i] = data[i].X;
1643                 y[i] = data[i].Y;
1644             }
1645 
1646             for (i = 0; i < n - 1; i++) {
1647                 dx.push(x[i + 1] - x[i]);
1648             }
1649             for (i = 0; i < n - 2; i++) {
1650                 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i]));
1651             }
1652 
1653             // ForwardSolve
1654             diag.push(2 * (dx[0] + dx[1]));
1655             z.push(delta[0]);
1656 
1657             for (i = 0; i < n - 3; i++) {
1658                 l = dx[i + 1] / diag[i];
1659                 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]);
1660                 z.push(delta[i + 1] - l * z[i]);
1661             }
1662 
1663             // BackwardSolve
1664             F[n - 3] = z[n - 3] / diag[n - 3];
1665             for (i = n - 4; i >= 0; i--) {
1666                 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i];
1667             }
1668 
1669             // Generate f''-Vector
1670             for (i = n - 3; i >= 0; i--) {
1671                 F[i + 1] = F[i];
1672             }
1673 
1674             // natural cubic spline
1675             F[0] = 0;
1676             F[n - 1] = 0;
1677 
1678             return F;
1679         },
1680 
1681         /**
1682          * Evaluate points on spline.
1683          * @param {Number,Array} x0 A single float value or an array of values to evaluate
1684          * @param {Array} x x values of knots
1685          * @param {Array} y y values of knots
1686          * @param {Array} F Second derivatives at knots, calculated by {@link JXG.Math.Numerics.splineDef}
1687          * @see #splineDef
1688          * @returns {Number,Array} A single value or an array, depending on what is given as x0.
1689          * @memberof JXG.Math.Numerics
1690          */
1691         splineEval: function (x0, x, y, F) {
1692             var i, j, a, b, c, d, x_,
1693                 n = Math.min(x.length, y.length),
1694                 l = 1,
1695                 asArray = false,
1696                 y0 = [];
1697 
1698             // number of points to be evaluated
1699             if (Type.isArray(x0)) {
1700                 l = x0.length;
1701                 asArray = true;
1702             } else {
1703                 x0 = [x0];
1704             }
1705 
1706             for (i = 0; i < l; i++) {
1707                 // is x0 in defining interval?
1708                 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) {
1709                     return NaN;
1710                 }
1711 
1712                 // determine part of spline in which x0 lies
1713                 for (j = 1; j < n; j++) {
1714                     if (x0[i] <= x[j]) {
1715                         break;
1716                     }
1717                 }
1718 
1719                 j -= 1;
1720 
1721                 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1];
1722                 // determine the coefficients of the polynomial in this interval
1723                 a = y[j];
1724                 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]);
1725                 c = F[j] / 2;
1726                 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j]));
1727                 // evaluate x0[i]
1728                 x_ = x0[i] - x[j];
1729                 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_);
1730                 y0.push(a + (b + (c + d * x_) * x_) * x_);
1731             }
1732 
1733             if (asArray) {
1734                 return y0;
1735             }
1736 
1737             return y0[0];
1738         },
1739 
1740         /**
1741          * Generate a string containing the function term of a polynomial.
1742          * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i.
1743          * @param {Number} deg Degree of the polynomial
1744          * @param {String} varname Name of the variable (usually 'x')
1745          * @param {Number} prec Precision
1746          * @returns {String} A string containg the function term of the polynomial.
1747          * @memberof JXG.Math.Numerics
1748          */
1749         generatePolynomialTerm: function (coeffs, deg, varname, prec) {
1750             var i, t = [];
1751 
1752             for (i = deg; i >= 0; i--) {
1753                 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']);
1754 
1755                 if (i > 1) {
1756                     t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']);
1757                 } else if (i === 1) {
1758                     t = t.concat(['*', varname, ' + ']);
1759                 }
1760             }
1761 
1762             return t.join('');
1763         },
1764 
1765         /**
1766          * Computes the polynomial through a given set of coordinates in Lagrange form.
1767          * Returns the Lagrange polynomials, see
1768          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
1769          * SIAM Review, Vol 46, No 3, (2004) 501-517.
1770          * <p>
1771          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
1772          * @param {Array} p Array of JXG.Points
1773          * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points.
1774          * @memberof JXG.Math.Numerics
1775          *
1776          * @example
1777          * var p = [];
1778          * p[0] = board.create('point', [-1,2], {size:4});
1779          * p[1] = board.create('point', [0,3], {size:4});
1780          * p[2] = board.create('point', [1,1], {size:4});
1781          * p[3] = board.create('point', [3,-1], {size:4});
1782          * var f = JXG.Math.Numerics.lagrangePolynomial(p);
1783          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1784          *
1785          * </pre><div id="JXGc058aa6b-74d4-41e1-af94-df06169a2d89" class="jxgbox" style="width: 300px; height: 300px;"></div>
1786          * <script type="text/javascript">
1787          *     (function() {
1788          *         var board = JXG.JSXGraph.initBoard('JXGc058aa6b-74d4-41e1-af94-df06169a2d89',
1789          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
1790          *     var p = [];
1791          *     p[0] = board.create('point', [-1,2], {size:4});
1792          *     p[1] = board.create('point', [0,3], {size:4});
1793          *     p[2] = board.create('point', [1,1], {size:4});
1794          *     p[3] = board.create('point', [3,-1], {size:4});
1795          *     var f = JXG.Math.Numerics.lagrangePolynomial(p);
1796          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1797          *
1798          *     })();
1799          *
1800          * </script><pre>
1801          *
1802          * @example
1803          * var points = [];
1804          * points[0] = board.create('point', [-1,2], {size:4});
1805          * points[1] = board.create('point', [0, 0], {size:4});
1806          * points[2] = board.create('point', [2, 1], {size:4});
1807          *
1808          * var f = JXG.Math.Numerics.lagrangePolynomial(points);
1809          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1810          * var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
1811          *
1812          * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eecbb" class="jxgbox" style="width: 300px; height: 300px;"></div>
1813          * <script type="text/javascript">
1814          *     (function() {
1815          *         var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eecbb',
1816          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
1817          *     var points = [];
1818          *     points[0] = board.create('point', [-1,2], {size:4});
1819          *     points[1] = board.create('point', [0, 0], {size:4});
1820          *     points[2] = board.create('point', [2, 1], {size:4});
1821          *
1822          *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
1823          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1824          *     var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
1825          *
1826          *     })();
1827          *
1828          * </script><pre>
1829          *
1830          */
1831         lagrangePolynomial: function (p) {
1832             var w = [],
1833                 that = this,
1834                 /** @ignore */
1835                 fct = function (x, suspendedUpdate) {
1836                     var i, // j,
1837                         k, xi, s, //M,
1838                         len = p.length,
1839                         num = 0,
1840                         denom = 0;
1841 
1842                     if (!suspendedUpdate) {
1843                         for (i = 0; i < len; i++) {
1844                             w[i] = 1.0;
1845                             xi = p[i].X();
1846 
1847                             for (k = 0; k < len; k++) {
1848                                 if (k !== i) {
1849                                     w[i] *= (xi - p[k].X());
1850                                 }
1851                             }
1852 
1853                             w[i] = 1 / w[i];
1854                         }
1855 
1856                         // M = [];
1857                         // for (k = 0; k < len; k++) {
1858                         //     M.push([1]);
1859                         // }
1860                     }
1861 
1862                     for (i = 0; i < len; i++) {
1863                         xi = p[i].X();
1864 
1865                         if (x === xi) {
1866                             return p[i].Y();
1867                         }
1868 
1869                         s = w[i] / (x - xi);
1870                         denom += s;
1871                         num += s * p[i].Y();
1872                     }
1873 
1874                     return num / denom;
1875                 };
1876 
1877             /**
1878              * Get the term of the Lagrange polynomial as string.
1879              * Calls {@link JXG.Math.Numerics#lagrangePolynomialTerm}.
1880              *
1881              * @name JXG.Math.Numerics#lagrangePolynomial.getTerm
1882              * @param {Number} digits Number of digits of the coefficients
1883              * @param {String} param Variable name
1884              * @param {String} dot Dot symbol
1885              * @returns {String} containing the term of Lagrange polynomial as string.
1886              * @see JXG.Math.Numerics#lagrangePolynomialTerm
1887              * @example
1888              * var points = [];
1889              * points[0] = board.create('point', [-1,2], {size:4});
1890              * points[1] = board.create('point', [0, 0], {size:4});
1891              * points[2] = board.create('point', [2, 1], {size:4});
1892              *
1893              * var f = JXG.Math.Numerics.lagrangePolynomial(points);
1894              * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1895              * var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
1896              *
1897              * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eeccf" class="jxgbox" style="width: 300px; height: 300px;"></div>
1898              * <script type="text/javascript">
1899              *     (function() {
1900              *         var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eeccf',
1901              *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
1902              *     var points = [];
1903              *     points[0] = board.create('point', [-1,2], {size:4});
1904              *     points[1] = board.create('point', [0, 0], {size:4});
1905              *     points[2] = board.create('point', [2, 1], {size:4});
1906              *
1907              *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
1908              *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1909              *     var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
1910              *
1911              *     })();
1912              *
1913              * </script><pre>
1914              *
1915              */
1916             fct.getTerm = function(digits, param, dot) {
1917                 return that.lagrangePolynomialTerm(p, digits, param, dot)();
1918             };
1919 
1920             return fct;
1921         },
1922         // fct.getTerm = that.lagrangePolynomialTerm(p, 2, 'x');
1923 
1924         /**
1925          * Determine the Lagrange polynomial through an array of points and
1926          * return the term of the polynomial as string.
1927          *
1928          * @param {Array} points Array of JXG.Points
1929          * @param {Number} digits Number of decimal digits of the coefficients
1930          * @param {String} param Name of the parameter. Default: 'x'.
1931          * @param {String} dot Multiplication symbol. Default: ' * '.
1932          * @returns {String} containing the Lagrange polynomial through
1933          *    the supplied points.
1934          * @memberof JXG.Math.Numerics
1935          *
1936          * @example
1937          * var points = [];
1938          * points[0] = board.create('point', [-1,2], {size:4});
1939          * points[1] = board.create('point', [0, 0], {size:4});
1940          * points[2] = board.create('point', [2, 1], {size:4});
1941          *
1942          * var f = JXG.Math.Numerics.lagrangePolynomial(points);
1943          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1944          *
1945          * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * ');
1946          * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16});
1947          *
1948          * </pre><div id="JXGd45e9e96-7526-486d-aa43-e1178d5f2baa" class="jxgbox" style="width: 300px; height: 300px;"></div>
1949          * <script type="text/javascript">
1950          *     (function() {
1951          *         var board = JXG.JSXGraph.initBoard('JXGd45e9e96-7526-486d-aa43-e1178d5f2baa',
1952          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
1953          *     var points = [];
1954          *     points[0] = board.create('point', [-1,2], {size:4});
1955          *     points[1] = board.create('point', [0, 0], {size:4});
1956          *     points[2] = board.create('point', [2, 1], {size:4});
1957          *
1958          *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
1959          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
1960          *
1961          *     var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * ');
1962          *     var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16});
1963          *
1964          *     })();
1965          *
1966          * </script><pre>
1967          *
1968          */
1969         lagrangePolynomialTerm: function(points, digits, param, dot) {
1970             return function() {
1971                 var len = points.length,
1972                     zeroes = [],
1973                     coeffs = [],
1974                     coeffs_sum = [],
1975                     isLeading = true,
1976                     n, t,
1977                     i, j, c, p;
1978 
1979                 param = param || 'x';
1980                 if (dot === undefined) {
1981                     dot = ' * ';
1982                 }
1983 
1984                 n = len - 1;  // (Max) degree of the polynomial
1985                 for (j = 0; j < len; j++) {
1986                     coeffs_sum[j] = 0;
1987                 }
1988 
1989                 for (i = 0; i < len; i++) {
1990                     c = points[i].Y();
1991                     p = points[i].X();
1992                     zeroes = [];
1993                     for (j = 0; j < len; j++) {
1994                         if (j !== i) {
1995                             c /= p - points[j].X();
1996                             zeroes.push(points[j].X());
1997                         }
1998                     }
1999                     coeffs = [1].concat(Mat.Vieta(zeroes));
2000                     for (j = 0; j < coeffs.length; j++) {
2001                         coeffs_sum[j] += (j%2===1?(-1):1) * coeffs[j] * c;
2002                     }
2003                 }
2004 
2005                 t = '';
2006                 for (j = 0; j < coeffs_sum.length; j++) {
2007                     c = coeffs_sum[j];
2008                     if (Math.abs(c) < Mat.eps) {
2009                         continue;
2010                     }
2011                     if (JXG.exists(digits)) {
2012                         c = Env._round10(c, -digits);
2013                     }
2014                     if (isLeading) {
2015                         t += (c > 0) ? (c) : ('-' + (-c));
2016                         isLeading = false;
2017                     } else {
2018                         t += (c > 0) ? (' + ' + c) : (' - ' + (-c));
2019                     }
2020 
2021                     if (n - j > 1) {
2022                         t += dot + param + '^' + (n - j);
2023                     } else if (n - j === 1) {
2024                         t += dot + param;
2025                     }
2026                 }
2027                 return t; // board.jc.manipulate('f = map(x) -> ' + t + ';');
2028             };
2029         },
2030 
2031         /**
2032          * Determine the coefficients of a cardinal spline polynom, See
2033          * http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections
2034          * @param  {Number} x1 point 1
2035          * @param  {Number} x2 point 2
2036          * @param  {Number} t1 tangent slope 1
2037          * @param  {Number} t2 tangent slope 2
2038          * @return {Array}    coefficents array c for the polynomial t maps to
2039          * c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t
2040          */
2041         _initCubicPoly: function(x1, x2, t1, t2) {
2042             return [
2043                 x1,
2044                 t1,
2045                 -3 * x1 + 3 * x2 - 2 * t1 - t2,
2046                 2 * x1 -  2 * x2 + t1 + t2
2047             ];
2048         },
2049 
2050         /**
2051          * Computes the cubic cardinal spline curve through a given set of points. The curve
2052          * is uniformly parametrized.
2053          * Two artificial control points at the beginning and the end are added.
2054          *
2055          * The implementation (especially  the centripetal parametrization) is from
2056          * http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections .
2057          * @param {Array} points Array consisting of JXG.Points.
2058          * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1.
2059          * tau=1/2 give Catmull-Rom splines.
2060          * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and
2061          * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal".
2062          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
2063          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value,
2064          * and a function simply returning the length of the points array
2065          * minus three.
2066          * @memberof JXG.Math.Numerics
2067         */
2068         CardinalSpline: function (points, tau_param, type) {
2069             var p,
2070                 coeffs = [],
2071                 makeFct,
2072                 tau, _tau,
2073                 that = this;
2074 
2075             if (Type.isFunction(tau_param)) {
2076                 _tau = tau_param;
2077             } else {
2078                 _tau = function () { return tau_param; };
2079             }
2080 
2081             if (type === undefined) {
2082                 type = 'uniform';
2083             }
2084 
2085             /** @ignore */
2086             makeFct = function (which) {
2087                 return function (t, suspendedUpdate) {
2088                     var s, c,
2089                         // control point at the beginning and at the end
2090                         first, last,
2091                         t1, t2, dt0, dt1, dt2,
2092                         // dx, dy,
2093                         len;
2094 
2095                     if (points.length < 2) {
2096                         return NaN;
2097                     }
2098 
2099                     if (!suspendedUpdate) {
2100                         tau = _tau();
2101 
2102                         // New point list p: [first, points ..., last]
2103                         first = {
2104                             X: function () { return 2 * points[0].X() - points[1].X(); },
2105                             Y: function () { return 2 * points[0].Y() - points[1].Y(); },
2106                             Dist: function(p) {
2107                                 var dx = this.X() - p.X(),
2108                                     dy = this.Y() - p.Y();
2109                                 return Math.sqrt(dx * dx + dy * dy);
2110                             }
2111                         };
2112 
2113                         last = {
2114                             X: function () { return 2 * points[points.length - 1].X() - points[points.length - 2].X(); },
2115                             Y: function () { return 2 * points[points.length - 1].Y() - points[points.length - 2].Y(); },
2116                             Dist: function(p) {
2117                                 var dx = this.X() - p.X(),
2118                                     dy = this.Y() - p.Y();
2119                                 return Math.sqrt(dx * dx + dy * dy);
2120                             }
2121                         };
2122 
2123                         p = [first].concat(points, [last]);
2124                         len = p.length;
2125 
2126                         coeffs[which] = [];
2127 
2128                         for (s = 0; s < len - 3; s++) {
2129                             if (type === 'centripetal') {
2130                                 // The order is important, since p[0].coords === undefined
2131                                 dt0 = p[s].Dist(p[s + 1]);
2132                                 dt1 = p[s + 2].Dist(p[s + 1]);
2133                                 dt2 = p[s + 3].Dist(p[s + 2]);
2134 
2135                                 dt0 = Math.sqrt(dt0);
2136                                 dt1 = Math.sqrt(dt1);
2137                                 dt2 = Math.sqrt(dt2);
2138 
2139                                 if (dt1 < Mat.eps) { dt1 = 1.0; }
2140                                 if (dt0 < Mat.eps) { dt0 = dt1; }
2141                                 if (dt2 < Mat.eps) { dt2 = dt1; }
2142 
2143                                 t1 = (p[s + 1][which]() - p[s][which]()) / dt0 -
2144                                      (p[s + 2][which]() - p[s][which]()) / (dt1 + dt0) +
2145                                      (p[s + 2][which]() - p[s + 1][which]()) / dt1;
2146 
2147                                 t2 = (p[s + 2][which]() - p[s + 1][which]()) / dt1 -
2148                                      (p[s + 3][which]() - p[s + 1][which]()) / (dt2 + dt1) +
2149                                      (p[s + 3][which]() - p[s + 2][which]()) / dt2;
2150 
2151                                 t1 *= dt1;
2152                                 t2 *= dt1;
2153 
2154                                 coeffs[which][s] = that._initCubicPoly(
2155                                      p[s + 1][which](),
2156                                      p[s + 2][which](),
2157                                      tau * t1,
2158                                      tau * t2
2159                                 );
2160                             } else {
2161                                 coeffs[which][s] = that._initCubicPoly(
2162                                     p[s + 1][which](),
2163                                     p[s + 2][which](),
2164                                     tau * (p[s + 2][which]() - p[s][which]()),
2165                                     tau * (p[s + 3][which]() - p[s + 1][which]())
2166                                 );
2167                             }
2168                         }
2169                     }
2170 
2171                     if (isNaN(t)) {
2172                         return NaN;
2173                     }
2174 
2175                     len = points.length;
2176                     // This is necessary for our advanced plotting algorithm:
2177                     if (t <= 0.0) {
2178                         return points[0][which]();
2179                     }
2180                     if (t >= len) {
2181                         return points[len - 1][which]();
2182                     }
2183 
2184                     s = Math.floor(t);
2185                     if (s === t) {
2186                         return points[s][which]();
2187                     }
2188 
2189                     t -= s;
2190                     c = coeffs[which][s];
2191                     if (c === undefined) {
2192                         return NaN;
2193                     }
2194 
2195                     return (((c[3] * t + c[2]) * t + c[1]) * t + c[0]);
2196                 };
2197             };
2198 
2199             return [makeFct('X'), makeFct('Y'), 0,
2200                 function () {
2201                     return points.length - 1;
2202                 }];
2203         },
2204 
2205         /**
2206          * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve
2207          * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5.
2208          * Two artificial control points at the beginning and the end are added.
2209          * @param {Array} points Array consisting of JXG.Points.
2210          * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and
2211          * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal".
2212          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
2213          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
2214          * returning the length of the points array minus three.
2215          * @memberof JXG.Math.Numerics
2216         */
2217         CatmullRomSpline: function (points, type) {
2218             return this.CardinalSpline(points, 0.5, type);
2219         },
2220 
2221         /**
2222          * Computes the regression polynomial of a given degree through a given set of coordinates.
2223          * Returns the regression polynomial function.
2224          * @param {Number,function,Slider} degree number, function or slider.
2225          * Either
2226          * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in
2227          * an array of {@link JXG.Point}s or {@link JXG.Coords}.
2228          * In the latter case, the <tt>dataY</tt> parameter will be ignored.
2229          * @param {Array} dataY Array containing the y-coordinates of the data set,
2230          * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree.
2231          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
2232          * @memberof JXG.Math.Numerics
2233          */
2234         regressionPolynomial: function (degree, dataX, dataY) {
2235             var coeffs, deg, dX, dY, inputType, fct,
2236                 term = '';
2237 
2238             // Slider
2239             if (Type.isPoint(degree) && Type.isFunction(degree.Value)) {
2240                 /** @ignore */
2241                 deg = function () {
2242                     return degree.Value();
2243                 };
2244             // function
2245             } else if (Type.isFunction(degree)) {
2246                 deg = degree;
2247             // number
2248             } else if (Type.isNumber(degree)) {
2249                 /** @ignore */
2250                 deg = function () {
2251                     return degree;
2252                 };
2253             } else {
2254                 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'.");
2255             }
2256 
2257             // Parameters degree, dataX, dataY
2258             if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) {
2259                 inputType = 0;
2260             // Parameters degree, point array
2261             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) {
2262                 inputType = 1;
2263             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) {
2264                 inputType = 2;
2265             } else {
2266                 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters.");
2267             }
2268 
2269             /** @ignore */
2270             fct = function (x, suspendedUpdate) {
2271                 var i, j, M, MT, y, B, c, s, d,
2272                     // input data
2273                     len = dataX.length;
2274 
2275                 d = Math.floor(deg());
2276 
2277                 if (!suspendedUpdate) {
2278                     // point list as input
2279                     if (inputType === 1) {
2280                         dX = [];
2281                         dY = [];
2282 
2283                         for (i = 0; i < len; i++) {
2284                             dX[i] = dataX[i].X();
2285                             dY[i] = dataX[i].Y();
2286                         }
2287                     }
2288 
2289                     if (inputType === 2) {
2290                         dX = [];
2291                         dY = [];
2292 
2293                         for (i = 0; i < len; i++) {
2294                             dX[i] = dataX[i].usrCoords[1];
2295                             dY[i] = dataX[i].usrCoords[2];
2296                         }
2297                     }
2298 
2299                     // check for functions
2300                     if (inputType === 0) {
2301                         dX = [];
2302                         dY = [];
2303 
2304                         for (i = 0; i < len; i++) {
2305                             if (Type.isFunction(dataX[i])) {
2306                                 dX.push(dataX[i]());
2307                             } else {
2308                                 dX.push(dataX[i]);
2309                             }
2310 
2311                             if (Type.isFunction(dataY[i])) {
2312                                 dY.push(dataY[i]());
2313                             } else {
2314                                 dY.push(dataY[i]);
2315                             }
2316                         }
2317                     }
2318 
2319                     M = [];
2320 
2321                     for (j = 0; j < len; j++) {
2322                         M.push([1]);
2323                     }
2324 
2325                     for (i = 1; i <= d; i++) {
2326                         for (j = 0; j < len; j++) {
2327                             M[j][i] = M[j][i - 1] * dX[j];
2328                         }
2329                     }
2330 
2331                     y = dY;
2332                     MT = Mat.transpose(M);
2333                     B = Mat.matMatMult(MT, M);
2334                     c = Mat.matVecMult(MT, y);
2335                     coeffs = Mat.Numerics.Gauss(B, c);
2336                     term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3);
2337                 }
2338 
2339                 // Horner's scheme to evaluate polynomial
2340                 s = coeffs[d];
2341 
2342                 for (i = d - 1; i >= 0; i--) {
2343                     s = (s * x + coeffs[i]);
2344                 }
2345 
2346                 return s;
2347             };
2348 
2349             fct.getTerm = function () {
2350                 return term;
2351             };
2352 
2353             return fct;
2354         },
2355 
2356         /**
2357          * Computes the cubic Bezier curve through a given set of points.
2358          * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}.
2359          * The points at position k with k mod 3 = 0 are the data points,
2360          * points at position k with k mod 3 = 1 or 2 are the control points.
2361          * @returns {Array} An array consisting of two functions of one parameter t which return the
2362          * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting
2363          * no parameters and returning one third of the length of the points.
2364          * @memberof JXG.Math.Numerics
2365          */
2366         bezier: function (points) {
2367             var len, flen,
2368                 /** @ignore */
2369                 makeFct = function (which) {
2370                     return function (t, suspendedUpdate) {
2371                         var z = Math.floor(t) * 3,
2372                             t0 = t % 1,
2373                             t1 = 1 - t0;
2374 
2375                         if (!suspendedUpdate) {
2376                             flen = 3 * Math.floor((points.length - 1) / 3);
2377                             len = Math.floor(flen / 3);
2378                         }
2379 
2380                         if (t < 0) {
2381                             return points[0][which]();
2382                         }
2383 
2384                         if (t >= len) {
2385                             return points[flen][which]();
2386                         }
2387 
2388                         if (isNaN(t)) {
2389                             return NaN;
2390                         }
2391 
2392                         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;
2393                     };
2394                 };
2395 
2396             return [makeFct('X'), makeFct('Y'), 0,
2397                 function () {
2398                     return Math.floor(points.length / 3);
2399                 }];
2400         },
2401 
2402         /**
2403          * Computes the B-spline curve of order k (order = degree+1) through a given set of points.
2404          * @param {Array} points Array consisting of JXG.Points.
2405          * @param {Number} order Order of the B-spline curve.
2406          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
2407          * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply
2408          * returning the length of the points array minus one.
2409          * @memberof JXG.Math.Numerics
2410          */
2411         bspline: function (points, order) {
2412             var knots,
2413                 _knotVector = function (n, k) {
2414                     var j,
2415                         kn = [];
2416 
2417                     for (j = 0; j < n + k + 1; j++) {
2418                         if (j < k) {
2419                             kn[j] = 0.0;
2420                         } else if (j <= n) {
2421                             kn[j] = j - k + 1;
2422                         } else {
2423                             kn[j] = n - k + 2;
2424                         }
2425                     }
2426 
2427                     return kn;
2428                 },
2429 
2430                 _evalBasisFuncs = function (t, kn, k, s) {
2431                     var i, j, a, b, den,
2432                         N = [];
2433 
2434                     if (kn[s] <= t && t < kn[s + 1]) {
2435                         N[s] = 1;
2436                     } else {
2437                         N[s] = 0;
2438                     }
2439 
2440                     for (i = 2; i <= k; i++) {
2441                         for (j = s - i + 1; j <= s; j++) {
2442                             if (j <= s - i + 1 || j < 0) {
2443                                 a = 0.0;
2444                             } else {
2445                                 a = N[j];
2446                             }
2447 
2448                             if (j >= s) {
2449                                 b = 0.0;
2450                             } else {
2451                                 b = N[j + 1];
2452                             }
2453 
2454                             den = kn[j + i - 1] - kn[j];
2455 
2456                             if (den === 0) {
2457                                 N[j] = 0;
2458                             } else {
2459                                 N[j] = (t - kn[j]) / den * a;
2460                             }
2461 
2462                             den = kn[j + i] - kn[j + 1];
2463 
2464                             if (den !== 0) {
2465                                 N[j] += (kn[j + i] - t) / den * b;
2466                             }
2467                         }
2468                     }
2469                     return N;
2470                 },
2471                 /** @ignore */
2472                 makeFct = function (which) {
2473                     return function (t, suspendedUpdate) {
2474                         var y, j, s, N = [],
2475                             len = points.length,
2476                             n = len - 1,
2477                             k = order;
2478 
2479                         if (n <= 0) {
2480                             return NaN;
2481                         }
2482 
2483                         if (n + 2 <= k) {
2484                             k = n + 1;
2485                         }
2486 
2487                         if (t <= 0) {
2488                             return points[0][which]();
2489                         }
2490 
2491                         if (t >= n - k + 2) {
2492                             return points[n][which]();
2493                         }
2494 
2495                         s = Math.floor(t) + k - 1;
2496                         knots = _knotVector(n, k);
2497                         N = _evalBasisFuncs(t, knots, k, s);
2498 
2499                         y = 0.0;
2500                         for (j = s - k + 1; j <= s; j++) {
2501                             if (j < len && j >= 0) {
2502                                 y += points[j][which]() * N[j];
2503                             }
2504                         }
2505 
2506                         return y;
2507                     };
2508                 };
2509 
2510             return [makeFct('X'), makeFct('Y'), 0,
2511                 function () {
2512                     return points.length - 1;
2513                 }];
2514         },
2515 
2516         /**
2517          * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through,
2518          * see {@link JXG.Curve#updateCurve}
2519          * and {@link JXG.Curve#hasPoint}.
2520          * @param {function} f Function in one variable to be differentiated.
2521          * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a
2522          * method of an object and contains a reference to its parent object via "this".
2523          * @returns {function} Derivative function of a given function f.
2524          * @memberof JXG.Math.Numerics
2525          */
2526         D: function (f, obj) {
2527             if (!Type.exists(obj)) {
2528                 return function (x, suspendedUpdate) {
2529                     var h = 0.00001,
2530                         h2 = (h * 2.0);
2531 
2532                     // Experiments with Richardsons rule
2533                     /*
2534                     var phi = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
2535                     var phi2;
2536                     h *= 0.5;
2537                     h2 *= 0.5;
2538                     phi2 = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
2539 
2540                     return phi2 + (phi2 - phi) / 3.0;
2541                     */
2542                     return (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
2543                 };
2544             }
2545 
2546             return function (x, suspendedUpdate) {
2547                 var h = 0.00001,
2548                     h2 = (h * 2.0);
2549 
2550                 return (f.apply(obj, [x + h, suspendedUpdate]) - f.apply(obj, [x - h, suspendedUpdate])) / h2;
2551             };
2552         },
2553 
2554         /**
2555          * Evaluate the function term for {@see #riemann}.
2556          * @private
2557          * @param {Number} x function argument
2558          * @param {function} f JavaScript function returning a number
2559          * @param {String} type Name of the Riemann sum type, e.g. 'lower', see {@see #riemann}.
2560          * @param {Number} delta Width of the bars in user coordinates
2561          * @returns {Number} Upper (delta > 0) or lower (delta < 0) value of the bar containing x of the Riemann sum.
2562          *
2563          * @memberof JXG.Math.Numerics
2564          */
2565         _riemannValue: function (x, f, type, delta) {
2566             var y, y1, x1, delta1;
2567 
2568             if (delta < 0) { // delta is negative if the lower function term is evaluated
2569                 if (type !== 'trapezoidal') {
2570                     x = x + delta;
2571                 }
2572                 delta *= -1;
2573                 if (type === 'lower') {
2574                     type = 'upper';
2575                 } else if (type === 'upper') {
2576                     type = 'lower';
2577                 }
2578             }
2579 
2580             delta1 = delta * 0.01; // for 'lower' and 'upper'
2581 
2582             if (type === 'right') {
2583                 y = f(x + delta);
2584             } else if (type === 'middle') {
2585                 y = f(x + delta * 0.5);
2586             } else if (type === 'left' || type === 'trapezoidal') {
2587                 y = f(x);
2588             } else if (type === 'lower') {
2589                 y = f(x);
2590 
2591                 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
2592                     y1 = f(x1);
2593 
2594                     if (y1 < y) {
2595                         y = y1;
2596                     }
2597                 }
2598 
2599                 y1 = f(x + delta);
2600                 if (y1 < y) {
2601                     y = y1;
2602                 }
2603             } else if (type === 'upper') {
2604                 y = f(x);
2605 
2606                 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
2607                     y1 = f(x1);
2608                     if (y1 > y) {
2609                         y = y1;
2610                     }
2611                 }
2612 
2613                 y1 = f(x + delta);
2614                 if (y1 > y) {
2615                     y = y1;
2616                 }
2617             } else if (type === 'random') {
2618                 y = f(x + delta * Math.random());
2619             } else if (type === 'simpson') {
2620                 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0;
2621             } else {
2622                 y = f(x);  // default is lower
2623             }
2624 
2625             return y;
2626         },
2627 
2628         /**
2629          * Helper function to create curve which displays Riemann sums.
2630          * Compute coordinates for the rectangles showing the Riemann sum.
2631          * @param {Function,Array} f Function or array of two functions.
2632          * If f is a function the integral of this function is approximated by the Riemann sum.
2633          * If f is an array consisting of two functions the area between the two functions is filled
2634          * by the Riemann sum bars.
2635          * @param {Number} n number of rectangles.
2636          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'.
2637          * @param {Number} start Left border of the approximation interval
2638          * @param {Number} end Right border of the approximation interval
2639          * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This
2640          * 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
2641          * rectangles.
2642          * @memberof JXG.Math.Numerics
2643          */
2644         riemann: function (gf, n, type, start, end) {
2645             var i, delta,
2646                 xarr = [],
2647                 yarr = [],
2648                 j = 0,
2649                 x = start, y,
2650                 sum = 0,
2651                 f, g,
2652                 ylow, yup;
2653 
2654             if (Type.isArray(gf)) {
2655                 g = gf[0];
2656                 f = gf[1];
2657             } else {
2658                 f = gf;
2659             }
2660 
2661             n = Math.floor(n);
2662 
2663             if (n <= 0) {
2664                 return [xarr, yarr, sum];
2665             }
2666 
2667             delta = (end - start) / n;
2668 
2669             // Upper bar ends
2670             for (i = 0; i < n; i++) {
2671                 y = this._riemannValue(x, f, type, delta);
2672                 xarr[j] = x;
2673                 yarr[j] = y;
2674 
2675                 j += 1;
2676                 x += delta;
2677                 if (type === 'trapezoidal') {
2678                     y = f(x);
2679                 }
2680                 xarr[j] = x;
2681                 yarr[j] = y;
2682 
2683                 j += 1;
2684             }
2685 
2686             // Lower bar ends
2687             for (i = 0; i < n; i++) {
2688                 if (g) {
2689                     y = this._riemannValue(x, g, type, -delta);
2690                 } else {
2691                     y = 0.0;
2692                 }
2693                 xarr[j] = x;
2694                 yarr[j] = y;
2695 
2696                 j += 1;
2697                 x -= delta;
2698                 if (type === 'trapezoidal' && g) {
2699                     y = g(x);
2700                 }
2701                 xarr[j] = x;
2702                 yarr[j] = y;
2703 
2704                 // Add the area of the bar to 'sum'
2705                 if (type !== 'trapezoidal') {
2706                     ylow = y;
2707                     yup = yarr[2 * (n - 1) - 2 * i];
2708                 } else {
2709                     yup = 0.5 * (f(x + delta) + f(x));
2710                     if (g) {
2711                         ylow = 0.5 * (g(x + delta) + g(x));
2712                     } else {
2713                         ylow = 0.0;
2714                     }
2715                 }
2716                 sum += (yup - ylow) * delta;
2717 
2718                 // Draw the vertical lines
2719                 j += 1;
2720                 xarr[j] = x;
2721                 yarr[j] = yarr[2 * (n - 1) - 2 * i];
2722 
2723                 j += 1;
2724             }
2725 
2726             return [xarr, yarr, sum];
2727         },
2728 
2729         /**
2730          * Approximate the integral by Riemann sums.
2731          * Compute the area described by the riemann sum rectangles.
2732          *
2733          * If there is an element of type {@link Riemannsum}, then it is more efficient
2734          * to use the method JXG.Curve.Value() of this element instead.
2735          *
2736          * @param {Function_Array} f Function or array of two functions.
2737          * If f is a function the integral of this function is approximated by the Riemann sum.
2738          * If f is an array consisting of two functions the area between the two functions is approximated
2739          * by the Riemann sum.
2740          * @param {Number} n number of rectangles.
2741          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'.
2742          *
2743          * @param {Number} start Left border of the approximation interval
2744          * @param {Number} end Right border of the approximation interval
2745          * @returns {Number} The sum of the areas of the rectangles.
2746          * @memberof JXG.Math.Numerics
2747          */
2748         riemannsum: function (f, n, type, start, end) {
2749             JXG.deprecated('Numerics.riemannsum()', 'Numerics.riemann()');
2750             return this.riemann(f, n, type, start, end)[2];
2751         },
2752 
2753         /**
2754          * Solve initial value problems numerically using Runge-Kutta-methods.
2755          * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
2756          * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
2757          * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
2758          * <pre>
2759          * {
2760          *     s: <Number>,
2761          *     A: <matrix>,
2762          *     b: <Array>,
2763          *     c: <Array>
2764          * }
2765          * </pre>
2766          * 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
2767          * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array.
2768          * @param {Array} I Interval on which to integrate.
2769          * @param {Number} N Number of evaluation points.
2770          * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
2771          * 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
2772          * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has.
2773          * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
2774          * @example
2775          * // A very simple autonomous system dx(t)/dt = x(t);
2776          * function f(t, x) {
2777          *     return x;
2778          * }
2779          *
2780          * // Solve it with initial value x(0) = 1 on the interval [0, 2]
2781          * // with 20 evaluation points.
2782          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
2783          *
2784          * // Prepare data for plotting the solution of the ode using a curve.
2785          * var dataX = [];
2786          * var dataY = [];
2787          * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
2788          * for(var i=0; i<data.length; i++) {
2789          *     dataX[i] = i*h;
2790          *     dataY[i] = data[i][0];
2791          * }
2792          * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
2793          * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
2794          * <script type="text/javascript">
2795          * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
2796          * function f(t, x) {
2797          *     // we have to copy the value.
2798          *     // return x; would just return the reference.
2799          *     return [x[0]];
2800          * }
2801          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
2802          * var dataX = [];
2803          * var dataY = [];
2804          * var h = 0.1;
2805          * for(var i=0; i<data.length; i++) {
2806          *     dataX[i] = i*h;
2807          *     dataY[i] = data[i][0];
2808          * }
2809          * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
2810          * </script><pre>
2811          * @memberof JXG.Math.Numerics
2812          */
2813         rungeKutta: function (butcher, x0, I, N, f) {
2814             var e, i, j, k, l, s,
2815                 x = [],
2816                 y = [],
2817                 h = (I[1] - I[0]) / N,
2818                 t = I[0],
2819                 dim = x0.length,
2820                 result = [],
2821                 r = 0;
2822 
2823             if (Type.isString(butcher)) {
2824                 butcher = predefinedButcher[butcher] || predefinedButcher.euler;
2825             }
2826             s = butcher.s;
2827 
2828             // don't change x0, so copy it
2829             for (e = 0; e < dim; e++) {
2830                 x[e] = x0[e];
2831             }
2832 
2833             for (i = 0; i < N; i++) {
2834                 // Optimization doesn't work for ODEs plotted using time
2835                 //        if((i % quotient == 0) || (i == N-1)) {
2836                 result[r] = [];
2837                 for (e = 0; e < dim; e++) {
2838                     result[r][e] = x[e];
2839                 }
2840 
2841                 r += 1;
2842                 k = [];
2843 
2844                 for (j = 0; j < s; j++) {
2845                     // init y = 0
2846                     for (e = 0; e < dim; e++) {
2847                         y[e] = 0.0;
2848                     }
2849 
2850 
2851                     // Calculate linear combination of former k's and save it in y
2852                     for (l = 0; l < j; l++) {
2853                         for (e = 0; e < dim; e++) {
2854                             y[e] += (butcher.A[j][l]) * h * k[l][e];
2855                         }
2856                     }
2857 
2858                     // add x(t) to y
2859                     for (e = 0; e < dim; e++) {
2860                         y[e] += x[e];
2861                     }
2862 
2863                     // calculate new k and add it to the k matrix
2864                     k.push(f(t + butcher.c[j] * h, y));
2865                 }
2866 
2867                 // init y = 0
2868                 for (e = 0; e < dim; e++) {
2869                     y[e] = 0.0;
2870                 }
2871 
2872                 for (l = 0; l < s; l++) {
2873                     for (e = 0; e < dim; e++) {
2874                         y[e] += butcher.b[l] * k[l][e];
2875                     }
2876                 }
2877 
2878                 for (e = 0; e < dim; e++) {
2879                     x[e] = x[e] + h * y[e];
2880                 }
2881 
2882                 t += h;
2883             }
2884 
2885             return result;
2886         },
2887 
2888         /**
2889          * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and
2890          * {@link JXG.Math.Numerics.chandrupatla}
2891          * @type Number
2892          * @default 80
2893          * @memberof JXG.Math.Numerics
2894          */
2895         maxIterationsRoot: 80,
2896 
2897         /**
2898          * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr}
2899          * @type Number
2900          * @default 500
2901          * @memberof JXG.Math.Numerics
2902          */
2903         maxIterationsMinimize: 500,
2904 
2905         /**
2906          * Given a value x_0, this function tries to find a second value x_1 such that
2907          * the function f has opposite signs at x_0 and x_1.
2908          * The return values have to be tested if the method succeeded.
2909          *
2910          * @param {Function} f Function, whose root is to be found
2911          * @param {Number} x0 Start value
2912          * @param {Object} object Parent object in case f is method of it
2913          * @returns {Array} [x_0, f(x_0), x_1, f(x_1)]
2914          *
2915          * @see JXG.Math.Numerics.fzero
2916          * @see JXG.Math.Numerics.chandrupatla
2917          *
2918          * @memberof JXG.Math.Numerics
2919          */
2920         findBracket: function(f, x0, object) {
2921             var a, aa, fa,
2922                 blist, b, fb,
2923                 u, fu,
2924                 i, len;
2925 
2926             if (Type.isArray(x0)) {
2927                 return x0;
2928             }
2929 
2930             a = x0;
2931             fa = f.call(object, a);
2932             // nfev += 1;
2933 
2934             // Try to get b, by trying several values related to a
2935             aa = (a === 0) ? 1 : a;
2936             blist = [
2937                 a - 0.1 * aa, a + 0.1 * aa,
2938                 a - 1, a + 1,
2939                 a - 0.5 * aa, a + 0.5 * aa,
2940                 a - 0.6 * aa, a + 0.6 * aa,
2941                 a - 1 * aa, a + 1 * aa,
2942                 a - 2 * aa, a + 2 * aa,
2943                 a - 5 * aa, a + 5 * aa,
2944                 a - 10 * aa, a + 10 * aa,
2945                 a - 50 * aa, a + 50 * aa,
2946                 a - 100 * aa, a + 100 * aa
2947             ];
2948             len = blist.length;
2949 
2950             for (i = 0; i < len; i++) {
2951                 b = blist[i];
2952                 fb = f.call(object, b);
2953                 // nfev += 1;
2954 
2955                 if (fa * fb <= 0) {
2956                     break;
2957                 }
2958             }
2959             if (b < a) {
2960                 u = a;
2961                 a = b;
2962                 b = u;
2963 
2964                 fu = fa;
2965                 fa = fb;
2966                 fb = fu;
2967             }
2968             return [a, fa, b, fb];
2969         },
2970 
2971         /**
2972          *
2973          * Find zero of an univariate function f.
2974          * @param {function} f Function, whose root is to be found
2975          * @param {Array,Number} x0  Start value or start interval enclosing the root
2976          * @param {Object} object Parent object in case f is method of it
2977          * @returns {Number} the approximation of the root
2978          * Algorithm:
2979          *  Brent's root finder from
2980          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
2981          *  computations. M., Mir, 1980, p.180 of the Russian edition
2982          *  http://www.netlib.org/c/brent.shar
2983          *
2984          * If x0 is an array containing lower and upper bound for the zero
2985          * algorithm 748 is applied. Otherwise, if x0 is a number,
2986          * the algorithm tries to bracket a zero of f starting from x0.
2987          * If this fails, we fall back to Newton's method.
2988          *
2989          * @see JXG.Math.Numerics.chandrupatla
2990          * @see JXG.Math.Numerics.root
2991          * @memberof JXG.Math.Numerics
2992          */
2993         fzero: function (f, x0, object) {
2994             var a, b, c,
2995                 d, e,
2996                 fa, fb, fc,
2997                 res,
2998                 prev_step, t1, cb, t2,
2999                 // Actual tolerance
3000                 tol_act,
3001                 // Interpolation step is calculated in the form p/q; division
3002                 // operations is delayed until the last moment
3003                 p, q,
3004                 // Step at this iteration
3005                 new_step,
3006                 eps = Mat.eps,
3007                 maxiter = this.maxIterationsRoot,
3008                 niter = 0;
3009                 // nfev = 0;
3010 
3011             if (Type.isArray(x0)) {
3012                 if (x0.length < 2) {
3013                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
3014                 }
3015 
3016                 a = x0[0];
3017                 fa = f.call(object, a);
3018                 // nfev += 1;
3019                 b = x0[1];
3020                 fb = f.call(object, b);
3021                 // nfev += 1;
3022             } else {
3023                 res = this.findBracket(f, x0, object);
3024                 a = res[0];
3025                 fa = res[1];
3026                 b = res[2];
3027                 fb = res[3];
3028             }
3029 
3030             if (Math.abs(fa) <= eps) {
3031                 return a;
3032             }
3033             if (Math.abs(fb) <= eps) {
3034                 return b;
3035             }
3036 
3037             if (fa * fb > 0) {
3038                 // Bracketing not successful, fall back to Newton's method or to fminbr
3039                 if (Type.isArray(x0)) {
3040                     return this.fminbr(f, [a, b], object);
3041                 }
3042 
3043                 return this.Newton(f, a, object);
3044             }
3045 
3046             // OK, we have enclosed a zero of f.
3047             // Now we can start Brent's method
3048             c = a;
3049             fc = fa;
3050 
3051             // Main iteration loop
3052             while (niter < maxiter) {
3053                 // Distance from the last but one to the last approximation
3054                 prev_step = b - a;
3055 
3056                 // Swap data for b to be the best approximation
3057                 if (Math.abs(fc) < Math.abs(fb)) {
3058                     a = b;
3059                     b = c;
3060                     c = a;
3061 
3062                     fa = fb;
3063                     fb = fc;
3064                     fc = fa;
3065                 }
3066                 tol_act = 2 * eps * Math.abs(b) + eps * 0.5;
3067                 new_step = (c - b) * 0.5;
3068 
3069                 if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) {
3070                     //  Acceptable approx. is found
3071                     return b;
3072                 }
3073 
3074                 // Decide if the interpolation can be tried
3075                 // If prev_step was large enough and was in true direction Interpolatiom may be tried
3076                 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
3077                     cb = c - b;
3078 
3079                     // If we have only two distinct points linear interpolation can only be applied
3080                     if (a === c) {
3081                         t1 = fb / fa;
3082                         p = cb * t1;
3083                         q = 1.0 - t1;
3084                     // Quadric inverse interpolation
3085                     } else {
3086                         q = fa / fc;
3087                         t1 = fb / fc;
3088                         t2 = fb / fa;
3089 
3090                         p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
3091                         q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
3092                     }
3093 
3094                     // p was calculated with the opposite sign; make p positive
3095                     if (p > 0) {
3096                         q = -q;
3097                     // and assign possible minus to q
3098                     } else {
3099                         p = -p;
3100                     }
3101 
3102                     // If b+p/q falls in [b,c] and isn't too large it is accepted
3103                     // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
3104                     if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) &&
3105                             p < Math.abs(prev_step * q * 0.5)) {
3106                         new_step = p / q;
3107                     }
3108                 }
3109 
3110                 // Adjust the step to be not less than tolerance
3111                 if (Math.abs(new_step) < tol_act) {
3112                     new_step = (new_step > 0) ? tol_act : -tol_act;
3113                 }
3114 
3115                 // Save the previous approx.
3116                 a  = b;
3117                 fa = fb;
3118                 b += new_step;
3119                 fb = f.call(object, b);
3120                 // Do step to a new approxim.
3121                 // nfev += 1;
3122 
3123                 // Adjust c for it to have a sign opposite to that of b
3124                 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
3125                     c = a;
3126                     fc = fa;
3127                 }
3128                 niter++;
3129             } // End while
3130 
3131             return b;
3132         },
3133 
3134         /**
3135          * Find zero of an univariate function f.
3136          * @param {function} f Function, whose root is to be found
3137          * @param {Array,Number} x0  Start value or start interval enclosing the root
3138          * @param {Object} object Parent object in case f is method of it
3139          * @returns {Number} the approximation of the root
3140          * Algorithm:
3141          * Chandrupatla's method, see
3142          * Tirupathi R. Chandrupatla,
3143          * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives",
3144          * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149.
3145          *
3146          * If x0 is an array containing lower and upper bound for the zero
3147          * algorithm 748 is applied. Otherwise, if x0 is a number,
3148          * the algorithm tries to bracket a zero of f starting from x0.
3149          * If this fails, we fall back to Newton's method.
3150          *
3151          * @see JXG.Math.Numerics.root
3152          * @see JXG.Math.Numerics.fzero
3153          * @memberof JXG.Math.Numerics
3154          */
3155         chandrupatla: function (f, x0, object) {
3156             var a, fa, b, fb, res,
3157                 niter = 0,
3158                 maxiter = this.maxIterationsRoot,
3159                 rand = (1 + Math.random() * 0.001),
3160                 t = 0.5 * rand,
3161                 eps = Mat.eps, // 1.e-10,
3162                 dlt = 0.00001,
3163                 x1, x2, x3, x,
3164                 f1, f2, f3, y,
3165                 xm, fm, tol, tl,
3166                 xi, ph, fl, fh,
3167                 AL, A, B, C, D;
3168 
3169             if (Type.isArray(x0)) {
3170                 if (x0.length < 2) {
3171                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
3172                 }
3173 
3174                 a = x0[0];
3175                 fa = f.call(object, a);
3176                 // nfev += 1;
3177                 b = x0[1];
3178                 fb = f.call(object, b);
3179                 // nfev += 1;
3180             } else {
3181                 res = this.findBracket(f, x0, object);
3182                 a = res[0];
3183                 fa = res[1];
3184                 b = res[2];
3185                 fb = res[3];
3186             }
3187 
3188             if (fa * fb > 0) {
3189                 // Bracketing not successful, fall back to Newton's method or to fminbr
3190                 if (Type.isArray(x0)) {
3191                     return this.fminbr(f, [a, b], object);
3192                 }
3193 
3194                 return this.Newton(f, a, object);
3195             }
3196 
3197             x1 = a;  x2 = b;
3198             f1 = fa; f2 = fb;
3199             do {
3200                 x = x1 + t * (x2 - x1);
3201                 y = f.call(object, x);
3202 
3203                 // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point
3204                 if (Math.sign(y) === Math.sign(f1)) {
3205                     x3 = x1; x1 = x;
3206                     f3 = f1; f1 = y;
3207                 } else {
3208                     x3 = x2; x2 = x1;
3209                     f3 = f2; f2 = f1;
3210                 }
3211                 x1 = x; f1 = y;
3212 
3213                 xm = x1; fm = f1;
3214                 if (Math.abs(f2) < Math.abs(f1)) {
3215                     xm = x2; fm = f2;
3216                 }
3217                 tol = 2 * eps * Math.abs(xm) + 0.5 * dlt;
3218                 tl = tol / Math.abs(x2 - x1);
3219                 if (tl > 0.5 || fm === 0) {
3220                     break;
3221                 }
3222                 // If inverse quadratic interpolation holds, use it
3223                 xi = (x1 - x2) / (x3 - x2);
3224                 ph = (f1 - f2) / (f3 - f2);
3225                 fl = 1 - Math.sqrt(1 - xi);
3226                 fh = Math.sqrt(xi);
3227                 if (fl < ph && ph < fh) {
3228                     AL = (x3 - x1) / (x2 - x1);
3229                     A = f1 / (f2 - f1);
3230                     B = f3 / (f2 - f3);
3231                     C = f1 / (f3 - f1);
3232                     D = f2 / (f3 - f2);
3233                     t = A * B + C * D * AL;
3234                 } else {
3235                     t = 0.5 * rand;
3236                 }
3237                 // Adjust t away from the interval boundary
3238                 if (t < tl) {
3239                     t = tl;
3240                 }
3241                 if (t > 1 - tl) {
3242                     t = 1 - tl;
3243                 }
3244                 niter++;
3245             } while (niter <= maxiter);
3246             // console.log(niter);
3247 
3248             return xm;
3249         },
3250 
3251         /**
3252          *
3253          * Find minimum of an univariate function f.
3254          * <p>
3255          * Algorithm:
3256          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
3257          *  computations. M., Mir, 1980, p.180 of the Russian edition
3258          *
3259          * @param {function} f Function, whose minimum is to be found
3260          * @param {Array} x0  Start interval enclosing the minimum
3261          * @param {Object} context Parent object in case f is method of it
3262          * @returns {Number} the approximation of the minimum value position
3263          * @memberof JXG.Math.Numerics
3264          **/
3265         fminbr: function (f, x0, context) {
3266             var a, b, x, v, w,
3267                 fx, fv, fw,
3268                 range, middle_range, tol_act, new_step,
3269                 p, q, t, ft,
3270                 // Golden section ratio
3271                 r = (3.0 - Math.sqrt(5.0)) * 0.5,
3272                 tol = Mat.eps,
3273                 sqrteps = Mat.eps, //Math.sqrt(Mat.eps),
3274                 maxiter = this.maxIterationsMinimize,
3275                 niter = 0;
3276                 // nfev = 0;
3277 
3278             if (!Type.isArray(x0) || x0.length < 2) {
3279                 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two.");
3280             }
3281 
3282             a = x0[0];
3283             b = x0[1];
3284             v = a + r * (b - a);
3285             fv = f.call(context, v);
3286 
3287             // First step - always gold section
3288             // nfev += 1;
3289             x = v;
3290             w = v;
3291             fx = fv;
3292             fw = fv;
3293 
3294             while (niter < maxiter) {
3295                 // Range over the interval in which we are looking for the minimum
3296                 range = b - a;
3297                 middle_range = (a + b) * 0.5;
3298 
3299                 // Actual tolerance
3300                 tol_act = sqrteps * Math.abs(x) + tol / 3.0;
3301 
3302                 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) {
3303                     // Acceptable approx. is found
3304                     return x;
3305                 }
3306 
3307                 // Obtain the golden section step
3308                 new_step = r * (x < middle_range ? b - x : a - x);
3309 
3310                 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried
3311                 if (Math.abs(x - w) >= tol_act) {
3312                     // Interpolation step is calculated as p/q;
3313                     // division operation is delayed until last moment
3314                     t = (x - w) * (fx - fv);
3315                     q = (x - v) * (fx - fw);
3316                     p = (x - v) * q - (x - w) * t;
3317                     q = 2 * (q - t);
3318 
3319                     if (q > 0) {                        // q was calculated with the op-
3320                         p = -p;                         // posite sign; make q positive
3321                     } else {                            // and assign possible minus to
3322                         q = -q;                         // p
3323                     }
3324                     if (Math.abs(p) < Math.abs(new_step * q) &&     // If x+p/q falls in [a,b]
3325                             p > q * (a - x + 2 * tol_act) &&        //  not too close to a and
3326                             p < q * (b - x - 2 * tol_act)) {        // b, and isn't too large
3327                         new_step = p / q;                          // it is accepted
3328                     }
3329                     // If p/q is too large then the
3330                     // golden section procedure can
3331                     // reduce [a,b] range to more
3332                     // extent
3333                 }
3334 
3335                 // Adjust the step to be not less than tolerance
3336                 if (Math.abs(new_step) < tol_act) {
3337                     if (new_step > 0) {
3338                         new_step = tol_act;
3339                     } else {
3340                         new_step = -tol_act;
3341                     }
3342                 }
3343 
3344                 // Obtain the next approximation to min
3345                 // and reduce the enveloping range
3346 
3347                 // Tentative point for the min
3348                 t = x + new_step;
3349                 ft = f.call(context, t);
3350                 // nfev += 1;
3351 
3352                 // t is a better approximation
3353                 if (ft <= fx) {
3354                     // Reduce the range so that t would fall within it
3355                     if (t < x) {
3356                         b = x;
3357                     } else {
3358                         a = x;
3359                     }
3360 
3361                     // Assign the best approx to x
3362                     v = w;
3363                     w = x;
3364                     x = t;
3365 
3366                     fv = fw;
3367                     fw = fx;
3368                     fx = ft;
3369                 // x remains the better approx
3370                 } else {
3371                     // Reduce the range enclosing x
3372                     if (t < x) {
3373                         a = t;
3374                     } else {
3375                         b = t;
3376                     }
3377 
3378                     if (ft <= fw || w === x) {
3379                         v = w;
3380                         w = t;
3381                         fv = fw;
3382                         fw = ft;
3383                     } else if (ft <= fv || v === x || v === w) {
3384                         v = t;
3385                         fv = ft;
3386                     }
3387                 }
3388                 niter += 1;
3389             }
3390 
3391             return x;
3392         },
3393 
3394         /**
3395          * Implements the Ramer-Douglas-Peucker algorithm.
3396          * It discards points which are not necessary from the polygonal line defined by the point array
3397          * pts. The computation is done in screen coordinates.
3398          * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points.
3399          * @param {Array} pts Array of {@link JXG.Coords}
3400          * @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>.
3401          * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points.
3402          * @memberof JXG.Math.Numerics
3403          */
3404         RamerDouglasPeucker: function (pts, eps) {
3405             var allPts = [], newPts = [], i, k, len,
3406 
3407                 /**
3408                  * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
3409                  * It searches for the point between index i and j which
3410                  * has the largest distance from the line between the points i and j.
3411                  * @param {Array} pts Array of {@link JXG.Coords}
3412                  * @param {Number} i Index of a point in pts
3413                  * @param {Number} j Index of a point in pts
3414                  * @ignore
3415                  * @private
3416                  */
3417                 findSplit = function (pts, i, j) {
3418                     var d, k, ci, cj, ck,
3419                         x0, y0, x1, y1,
3420                         den, lbda,
3421                         huge = 10000,
3422                         dist = 0,
3423                         f = i;
3424 
3425                     if (j - i < 2) {
3426                         return [-1.0, 0];
3427                     }
3428 
3429                     ci = pts[i].scrCoords;
3430                     cj = pts[j].scrCoords;
3431 
3432                     if (isNaN(ci[1]) || isNaN(ci[2])) {
3433                         return [NaN, i];
3434                     }
3435                     if (isNaN(cj[1]) || isNaN(cj[2])) {
3436                         return [NaN, j];
3437                     }
3438 
3439                     for (k = i + 1; k < j; k++) {
3440                         ck = pts[k].scrCoords;
3441                         if (isNaN(ck[1]) || isNaN(ck[2])) {
3442                             return [NaN, k];
3443                         }
3444 
3445                         x0 = ck[1] - ci[1];
3446                         y0 = ck[2] - ci[2];
3447                         x1 = cj[1] - ci[1];
3448                         y1 = cj[2] - ci[2];
3449                         x0 = x0 === Infinity ? huge : x0;
3450                         y0 = y0 === Infinity ? huge : y0;
3451                         x1 = x1 === Infinity ? huge : x1;
3452                         y1 = y1 === Infinity ? huge : y1;
3453                         x0 = x0 === -Infinity ? -huge : x0;
3454                         y0 = y0 === -Infinity ? -huge : y0;
3455                         x1 = x1 === -Infinity ? -huge : x1;
3456                         y1 = y1 === -Infinity ? -huge : y1;
3457                         den = x1 * x1 + y1 * y1;
3458 
3459                         if (den >= Mat.eps) {
3460                             lbda = (x0 * x1 + y0 * y1) / den;
3461 
3462                             if (lbda < 0.0) {
3463                                 lbda = 0.0;
3464                             } else if (lbda > 1.0) {
3465                                 lbda = 1.0;
3466                             }
3467 
3468                             x0 = x0 - lbda * x1;
3469                             y0 = y0 - lbda * y1;
3470                             d = x0 * x0 + y0 * y0;
3471                         } else {
3472                             lbda = 0.0;
3473                             d = x0 * x0 + y0 * y0;
3474                         }
3475 
3476                         if (d > dist) {
3477                             dist = d;
3478                             f = k;
3479                         }
3480                     }
3481                     return [Math.sqrt(dist), f];
3482                 },
3483 
3484                 /**
3485                  * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
3486                  * It runs recursively through the point set and searches the
3487                  * point which has the largest distance from the line between the first point and
3488                  * the last point. If the distance from the line is greater than eps, this point is
3489                  * included in our new point set otherwise it is discarded.
3490                  * If it is taken, we recursively apply the subroutine to the point set before
3491                  * and after the chosen point.
3492                  * @param {Array} pts Array of {@link JXG.Coords}
3493                  * @param {Number} i Index of an element of pts
3494                  * @param {Number} j Index of an element of pts
3495                  * @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>.
3496                  * @param {Array} newPts Array of {@link JXG.Coords}
3497                  * @ignore
3498                  * @private
3499                  */
3500                 RDP = function (pts, i, j, eps, newPts) {
3501                     var result = findSplit(pts, i, j),
3502                         k = result[1];
3503 
3504                     if (isNaN(result[0])) {
3505                         RDP(pts, i, k - 1, eps, newPts);
3506                         newPts.push(pts[k]);
3507                         do {
3508                             ++k;
3509                         } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2]));
3510                         if (k <= j) {
3511                             newPts.push(pts[k]);
3512                         }
3513                         RDP(pts, k + 1, j, eps, newPts);
3514                     } else if (result[0] > eps) {
3515                         RDP(pts, i, k, eps, newPts);
3516                         RDP(pts, k, j, eps, newPts);
3517                     } else {
3518                         newPts.push(pts[j]);
3519                     }
3520                 };
3521 
3522             len = pts.length;
3523 
3524             i = 0;
3525             while (true) {
3526                 // Search for the next point without NaN coordinates
3527                 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) {
3528                     i += 1;
3529                 }
3530                 // Search for the next position of a NaN point
3531                 k = i + 1;
3532                 while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) {
3533                     k += 1;
3534                 }
3535                 k--;
3536 
3537                 // Only proceed if something is left
3538                 if (i < len && k > i) {
3539                     newPts = [];
3540                     newPts[0] = pts[i];
3541                     RDP(pts, i, k, eps, newPts);
3542                     allPts = allPts.concat(newPts);
3543                 }
3544                 if (i >= len) {
3545                     break;
3546                 }
3547                 // Push the NaN point
3548                 if (k < len - 1) {
3549                     allPts.push(pts[k + 1]);
3550                 }
3551                 i = k + 1;
3552             }
3553 
3554             return allPts;
3555         },
3556 
3557         /**
3558          * Old name for the implementation of the Ramer-Douglas-Peucker algorithm.
3559          * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker}
3560          * @memberof JXG.Math.Numerics
3561          */
3562         RamerDouglasPeuker: function (pts, eps) {
3563             JXG.deprecated('Numerics.RamerDouglasPeuker()', 'Numerics.RamerDouglasPeucker()');
3564             return this.RamerDouglasPeucker(pts, eps);
3565         },
3566 
3567         /**
3568          * Implements the Visvalingam-Whyatt algorithm.
3569          * See M. Visvalingam, J. D. Whyatt:
3570          * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992
3571          *
3572          * The algorithm discards points which are not necessary from the polygonal line defined by the point array
3573          * pts (consisting of type JXG.Coords).
3574          * @param {Array} pts Array of {@link JXG.Coords}
3575          * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will
3576          *    be taken in any case.
3577          * @returns {Array} An array containing points which approximates the curve defined by pts.
3578          * @memberof JXG.Math.Numerics
3579          *
3580          * @example
3581          *     var i, p = [];
3582          *     for (i = 0; i < 5; ++i) {
3583          *         p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
3584          *     }
3585          *
3586          *     // Plot a cardinal spline curve
3587          *     var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
3588          *     var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
3589          *
3590          *     var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
3591          *     c.updateDataArray = function() {
3592          *         var i, len, points;
3593          *
3594          *         // Reduce number of intermediate points with Visvakingam-Whyatt to 6
3595          *         points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
3596          *         // Plot the remaining points
3597          *         len = points.length;
3598          *         this.dataX = [];
3599          *         this.dataY = [];
3600          *         for (i = 0; i < len; i++) {
3601          *             this.dataX.push(points[i].usrCoords[1]);
3602          *             this.dataY.push(points[i].usrCoords[2]);
3603          *         }
3604          *     };
3605          *     board.update();
3606          *
3607          * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div>
3608          * <script type="text/javascript">
3609          *     (function() {
3610          *         var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb',
3611          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
3612          *
3613          *         var i, p = [];
3614          *         for (i = 0; i < 5; ++i) {
3615          *             p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
3616          *         }
3617          *
3618          *         // Plot a cardinal spline curve
3619          *         var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
3620          *         var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
3621          *
3622          *         var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
3623          *         c.updateDataArray = function() {
3624          *             var i, len, points;
3625          *
3626          *             // Reduce number of intermediate points with Visvakingam-Whyatt to 6
3627          *             points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
3628          *             // Plot the remaining points
3629          *             len = points.length;
3630          *             this.dataX = [];
3631          *             this.dataY = [];
3632          *             for (i = 0; i < len; i++) {
3633          *                 this.dataX.push(points[i].usrCoords[1]);
3634          *                 this.dataY.push(points[i].usrCoords[2]);
3635          *             }
3636          *         };
3637          *         board.update();
3638          *
3639          *     })();
3640          *
3641          * </script><pre>
3642          *
3643          */
3644         Visvalingam: function(pts, numPoints) {
3645             var i, len, vol, lastVol,
3646                 linkedList = [],
3647                 heap = [],
3648                 points = [],
3649                 lft, rt, lft2, rt2,
3650                 obj;
3651 
3652             len = pts.length;
3653 
3654             // At least one intermediate point is needed
3655             if (len <= 2) {
3656                 return pts;
3657             }
3658 
3659             // Fill the linked list
3660             // Add first point to the linked list
3661             linkedList[0] = {
3662                     used: true,
3663                     lft: null,
3664                     node: null
3665                 };
3666 
3667             // Add all intermediate points to the linked list,
3668             // whose triangle area is nonzero.
3669             lft = 0;
3670             for (i = 1; i < len - 1; i++) {
3671                 vol = Math.abs(JXG.Math.Numerics.det([pts[i - 1].usrCoords,
3672                                               pts[i].usrCoords,
3673                                               pts[i + 1].usrCoords]));
3674                 if (!isNaN(vol)) {
3675                     obj = {
3676                         v: vol,
3677                         idx: i
3678                     };
3679                     heap.push(obj);
3680                     linkedList[i] = {
3681                             used: true,
3682                             lft: lft,
3683                             node: obj
3684                         };
3685                     linkedList[lft].rt = i;
3686                     lft = i;
3687                 }
3688             }
3689 
3690             // Add last point to the linked list
3691             linkedList[len - 1] = {
3692                     used: true,
3693                     rt: null,
3694                     lft: lft,
3695                     node: null
3696                 };
3697             linkedList[lft].rt = len - 1;
3698 
3699             // Remove points until only numPoints intermediate points remain
3700             lastVol = -Infinity;
3701             while (heap.length > numPoints) {
3702                 // Sort the heap with the updated volume values
3703                 heap.sort(function(a, b) {
3704                     // descending sort
3705                     return b.v - a.v;
3706                 });
3707 
3708                 // Remove the point with the smallest triangle
3709                 i = heap.pop().idx;
3710                 linkedList[i].used = false;
3711                 lastVol = linkedList[i].node.v;
3712 
3713                 // Update the pointers of the linked list
3714                 lft = linkedList[i].lft;
3715                 rt = linkedList[i].rt;
3716                 linkedList[lft].rt = rt;
3717                 linkedList[rt].lft = lft;
3718 
3719                 // Update the values for the volumes in the linked list
3720                 lft2 = linkedList[lft].lft;
3721                 if (lft2 !== null) {
3722                     vol = Math.abs(JXG.Math.Numerics.det(
3723                                 [pts[lft2].usrCoords,
3724                                  pts[lft].usrCoords,
3725                                  pts[rt].usrCoords]));
3726 
3727                     linkedList[lft].node.v = (vol >= lastVol) ? vol : lastVol;
3728                 }
3729                 rt2 = linkedList[rt].rt;
3730                 if (rt2 !== null) {
3731                     vol = Math.abs(JXG.Math.Numerics.det(
3732                                 [pts[lft].usrCoords,
3733                                  pts[rt].usrCoords,
3734                                  pts[rt2].usrCoords]));
3735 
3736                     linkedList[rt].node.v = (vol >= lastVol) ? vol : lastVol;
3737                 }
3738             }
3739 
3740             // Return an array with the remaining points
3741             i = 0;
3742             points = [pts[i]];
3743             do {
3744                 i = linkedList[i].rt;
3745                 points.push(pts[i]);
3746             } while (linkedList[i].rt !== null);
3747 
3748             return points;
3749         }
3750     };
3751 
3752     return Mat.Numerics;
3753 });
3754