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