1 /*
  2     Copyright 2008-2026
  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 <https://www.gnu.org/licenses/>
 29     and <https://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 /**
 37  * @fileoverview In this file the namespace Math.Numerics is defined, which holds numerical
 38  * algorithms for solving linear equations etc.
 39  */
 40 
 41 import JXG from "../jxg.js";
 42 import Type from "../utils/type.js";
 43 import Env from "../utils/env.js";
 44 import Mat from "./math.js";
 45 
 46 // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order).
 47 var predefinedButcher = {
 48     rk4: {
 49         s: 4,
 50         A: [
 51             [0, 0, 0, 0],
 52             [0.5, 0, 0, 0],
 53             [0, 0.5, 0, 0],
 54             [0, 0, 1, 0]
 55         ],
 56         b: [1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0],
 57         c: [0, 0.5, 0.5, 1]
 58     },
 59     heun: {
 60         s: 2,
 61         A: [
 62             [0, 0],
 63             [1, 0]
 64         ],
 65         b: [0.5, 0.5],
 66         c: [0, 1]
 67     },
 68     euler: {
 69         s: 1,
 70         A: [[0]],
 71         b: [1],
 72         c: [0]
 73     }
 74 };
 75 
 76 /**
 77  * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables.
 78  * @name JXG.Math.Numerics
 79  * @exports Mat.Numerics as JXG.Math.Numerics
 80  * @namespace
 81  */
 82 Mat.Numerics = {
 83     //JXG.extend(Mat.Numerics, /** @lends JXG.Math.Numerics */ {
 84     /**
 85      * Solves a system of linear equations given by A and b using the Gauss-Jordan-elimination.
 86      * The algorithm runs in-place. I.e. the entries of A and b are changed.
 87      * @param {Array} A Square matrix represented by an array of rows, containing the coefficients of the lineare equation system.
 88      * @param {Array} b A vector containing the linear equation system's right hand side.
 89      * @throws {Error} If a non-square-matrix is given or if b has not the right length or A's rank is not full.
 90      * @returns {Array} A vector that solves the linear equation system.
 91      * @memberof JXG.Math.Numerics
 92      */
 93     Gauss: function (A, b) {
 94         var i,
 95             j,
 96             k,
 97             // copy the matrix to prevent changes in the original
 98             Acopy,
 99             // solution vector, to prevent changing b
100             x,
101             eps = Mat.eps,
102             // number of columns of A
103             n = A.length > 0 ? A[0].length : 0;
104 
105         if (n !== b.length || n !== A.length) {
106             throw new Error(
107                 "JXG.Math.Numerics.Gauss: Dimensions don't match. A must be a square matrix and b must be of the same length as A."
108             );
109         }
110 
111         // initialize solution vector
112         Acopy = [];
113         x = b.slice(0, n);
114 
115         for (i = 0; i < n; i++) {
116             Acopy[i] = A[i].slice(0, n);
117         }
118 
119         // Gauss-Jordan-elimination
120         for (j = 0; j < n; j++) {
121             for (i = n - 1; i > j; i--) {
122                 // Is the element which is to eliminate greater than zero?
123                 if (Math.abs(Acopy[i][j]) > eps) {
124                     // Equals pivot element zero?
125                     if (Math.abs(Acopy[j][j]) < eps) {
126                         // At least numerically, so we have to exchange the rows
127                         Type.swap(Acopy, i, j);
128                         Type.swap(x, i, j);
129                     } else {
130                         // Saves the L matrix of the LR-decomposition. unnecessary.
131                         Acopy[i][j] /= Acopy[j][j];
132                         // Transform right-hand-side b
133                         x[i] -= Acopy[i][j] * x[j];
134 
135                         // subtract the multiple of A[i][j] / A[j][j] of the j-th row from the i-th.
136                         for (k = j + 1; k < n; k++) {
137                             Acopy[i][k] -= Acopy[i][j] * Acopy[j][k];
138                         }
139                     }
140                 }
141             }
142 
143             // The absolute values of all coefficients below the j-th row in the j-th column are smaller than JXG.Math.eps.
144             if (Math.abs(Acopy[j][j]) < eps) {
145                 throw new Error(
146                     "JXG.Math.Numerics.Gauss(): The given matrix seems to be singular."
147                 );
148             }
149         }
150 
151         this.backwardSolve(Acopy, x, true);
152 
153         return x;
154     },
155 
156     /**
157      * Solves a system of linear equations given by the right triangular matrix R and vector b.
158      * @param {Array} R Right triangular matrix represented by an array of rows. All entries a_(i,j) with i < j are ignored.
159      * @param {Array} b Right hand side of the linear equation system.
160      * @param {Boolean} [canModify=false] If true, the right hand side vector is allowed to be changed by this method.
161      * @returns {Array} An array representing a vector that solves the system of linear equations.
162      * @memberof JXG.Math.Numerics
163      */
164     backwardSolve: function (R, b, canModify) {
165         var x, m, n, i, j;
166 
167         if (canModify) {
168             x = b;
169         } else {
170             x = b.slice(0, b.length);
171         }
172 
173         // m: number of rows of R
174         // n: number of columns of R
175         m = R.length;
176         n = R.length > 0 ? R[0].length : 0;
177 
178         for (i = m - 1; i >= 0; i--) {
179             for (j = n - 1; j > i; j--) {
180                 x[i] -= R[i][j] * x[j];
181             }
182             x[i] /= R[i][i];
183         }
184 
185         return x;
186     },
187 
188     /**
189      *  Gauss-Bareiss algorithm to compute the
190      *  determinant of matrix without fractions.
191      *  See Henri Cohen, "A Course in Computational
192      *  Algebraic Number Theory (Graduate texts
193      *  in mathematics; 138)", Springer-Verlag,
194      *  ISBN 3-540-55640-0 / 0-387-55640-0
195      *  Third, Corrected Printing 1996
196      *  "Algorithm 2.2.6", pg. 52-53
197      *
198      * @param {Array} mat Matrix
199      * @returns Number
200      * @private
201      * @memberof JXG.Math.Numerics
202      */
203     gaussBareiss: function (mat) {
204         var k, c, s,
205             i, j, p,
206             n, M, t,
207             eps = Mat.eps;
208 
209         n = mat.length;
210 
211         if (n <= 0) {
212             return 0;
213         }
214 
215         if (mat[0].length < n) {
216             n = mat[0].length;
217         }
218 
219         // Copy the input matrix to M
220         M = [];
221 
222         for (i = 0; i < n; i++) {
223             M[i] = mat[i].slice(0, n);
224         }
225 
226         c = 1;
227         s = 1;
228 
229         for (k = 0; k < n - 1; k++) {
230             p = M[k][k];
231 
232             // Pivot step
233             if (Math.abs(p) < eps) {
234                 for (i = k + 1; i < n; i++) {
235                     if (Math.abs(M[i][k]) >= eps) {
236                         break;
237                     }
238                 }
239 
240                 // No nonzero entry found in column k -> det(M) = 0
241                 if (i === n) {
242                     return 0.0;
243                 }
244 
245                 // swap row i and k partially
246                 for (j = k; j < n; j++) {
247                     t = M[i][j];
248                     M[i][j] = M[k][j];
249                     M[k][j] = t;
250                 }
251                 s = -s;
252                 p = M[k][k];
253             }
254 
255             // Main step
256             for (i = k + 1; i < n; i++) {
257                 for (j = k + 1; j < n; j++) {
258                     t = p * M[i][j] - M[i][k] * M[k][j];
259                     M[i][j] = t / c;
260                 }
261             }
262 
263             c = p;
264         }
265 
266         return s * M[n - 1][n - 1];
267     },
268 
269     /**
270      * Computes the determinant of a square nxn matrix with the
271      * Gauss-Bareiss algorithm.
272      * @param {Array} mat Matrix.
273      * @returns {Number} The determinant pf the matrix mat.
274      *                   The empty matrix returns 0.
275      * @memberof JXG.Math.Numerics
276      */
277     det: function (mat) {
278         var n = mat.length;
279 
280         if (n === 2 && mat[0].length === 2) {
281             return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
282         }
283 
284         return this.gaussBareiss(mat);
285     },
286 
287     /**
288      * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method
289      * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990
290      * @param {Array} Ain A symmetric 3x3 matrix.
291      * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors.
292      * @memberof JXG.Math.Numerics
293      */
294     Jacobi: function (Ain) {
295         var i,
296             j,
297             k,
298             aa,
299             si,
300             co,
301             tt,
302             ssum,
303             amax,
304             eps = Mat.eps * Mat.eps,
305             sum = 0.0,
306             n = Ain.length,
307             V = [
308                 [0, 0, 0],
309                 [0, 0, 0],
310                 [0, 0, 0]
311             ],
312             A = [
313                 [0, 0, 0],
314                 [0, 0, 0],
315                 [0, 0, 0]
316             ],
317             nloops = 0;
318 
319         // Initialization. Set initial Eigenvectors.
320         for (i = 0; i < n; i++) {
321             for (j = 0; j < n; j++) {
322                 V[i][j] = 0.0;
323                 A[i][j] = Ain[i][j];
324                 sum += Math.abs(A[i][j]);
325             }
326             V[i][i] = 1.0;
327         }
328 
329         // Trivial problems
330         if (n === 1) {
331             return [A, V];
332         }
333 
334         if (sum <= 0.0) {
335             return [A, V];
336         }
337 
338         sum /= n * n;
339 
340         // Reduce matrix to diagonal
341         do {
342             ssum = 0.0;
343             amax = 0.0;
344             for (j = 1; j < n; j++) {
345                 for (i = 0; i < j; i++) {
346                     // Check if A[i][j] is to be reduced
347                     aa = Math.abs(A[i][j]);
348 
349                     if (aa > amax) {
350                         amax = aa;
351                     }
352 
353                     ssum += aa;
354 
355                     if (aa >= eps) {
356                         // calculate rotation angle
357                         aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5;
358                         si = Math.sin(aa);
359                         co = Math.cos(aa);
360 
361                         // Modify 'i' and 'j' columns
362                         for (k = 0; k < n; k++) {
363                             tt = A[k][i];
364                             A[k][i] = co * tt + si * A[k][j];
365                             A[k][j] = -si * tt + co * A[k][j];
366                             tt = V[k][i];
367                             V[k][i] = co * tt + si * V[k][j];
368                             V[k][j] = -si * tt + co * V[k][j];
369                         }
370 
371                         // Modify diagonal terms
372                         A[i][i] = co * A[i][i] + si * A[j][i];
373                         A[j][j] = -si * A[i][j] + co * A[j][j];
374                         A[i][j] = 0.0;
375 
376                         // Make 'A' matrix symmetrical
377                         for (k = 0; k < n; k++) {
378                             A[i][k] = A[k][i];
379                             A[j][k] = A[k][j];
380                         }
381                         // A[i][j] made zero by rotation
382                     }
383                 }
384             }
385             nloops += 1;
386         } while (Math.abs(ssum) / sum > eps && nloops < 2000);
387 
388         return [A, V];
389     },
390 
391     /**
392      * Calculates the integral of function f over interval using Newton-Cotes-algorithm.
393      * @param {Array} interval The integration interval, e.g. [0, 3].
394      * @param {function} f A function which takes one argument of type number and returns a number.
395      * @param {Object} [config] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type
396      * with value being either 'trapez', 'simpson', or 'milne'.
397      * @param {Number} [config.number_of_nodes=28]
398      * @param {String} [config.integration_type='milne'] Possible values are 'milne', 'simpson', 'trapez'
399      * @returns {Number} Integral value of f over interval
400      * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use
401      * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4.
402      * @example
403      * function f(x) {
404      *   return x*x;
405      * }
406      *
407      * // calculates integral of <tt>f</tt> from 0 to 2.
408      * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f);
409      *
410      * // the same with an anonymous function
411      * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; });
412      *
413      * // use trapez rule with 16 nodes
414      * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f,
415      *                                   {number_of_nodes: 16, integration_type: 'trapez'});
416      * @memberof JXG.Math.Numerics
417      */
418     NewtonCotes: function (interval, f, config) {
419         var evaluation_point,
420             i,
421             number_of_intervals,
422             integral_value = 0.0,
423             number_of_nodes =
424                 config && Type.isNumber(config.number_of_nodes) ? config.number_of_nodes : 28,
425             available_types = { trapez: true, simpson: true, milne: true },
426             integration_type =
427                 config &&
428                     config.integration_type &&
429                     available_types.hasOwnProperty(config.integration_type) &&
430                     available_types[config.integration_type]
431                     ? config.integration_type
432                     : "milne",
433             step_size = (interval[1] - interval[0]) / number_of_nodes;
434 
435         switch (integration_type) {
436             case "trapez":
437                 integral_value = (f(interval[0]) + f(interval[1])) * 0.5;
438                 evaluation_point = interval[0];
439 
440                 for (i = 0; i < number_of_nodes - 1; i++) {
441                     evaluation_point += step_size;
442                     integral_value += f(evaluation_point);
443                 }
444 
445                 integral_value *= step_size;
446                 break;
447             case "simpson":
448                 if (number_of_nodes % 2 > 0) {
449                     throw new Error(
450                         "JSXGraph:  INT_SIMPSON requires config.number_of_nodes dividable by 2."
451                     );
452                 }
453 
454                 number_of_intervals = number_of_nodes / 2.0;
455                 integral_value = 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 += 2.0 * step_size;
460                     integral_value += 2.0 * f(evaluation_point);
461                 }
462 
463                 evaluation_point = interval[0] - step_size;
464 
465                 for (i = 0; i < number_of_intervals; i++) {
466                     evaluation_point += 2.0 * step_size;
467                     integral_value += 4.0 * f(evaluation_point);
468                 }
469 
470                 integral_value *= step_size / 3.0;
471                 break;
472             default:
473                 if (number_of_nodes % 4 > 0) {
474                     throw new Error(
475                         "JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4"
476                     );
477                 }
478 
479                 number_of_intervals = number_of_nodes * 0.25;
480                 integral_value = 7.0 * (f(interval[0]) + f(interval[1]));
481                 evaluation_point = interval[0];
482 
483                 for (i = 0; i < number_of_intervals - 1; i++) {
484                     evaluation_point += 4.0 * step_size;
485                     integral_value += 14.0 * f(evaluation_point);
486                 }
487 
488                 evaluation_point = interval[0] - 3.0 * step_size;
489 
490                 for (i = 0; i < number_of_intervals; i++) {
491                     evaluation_point += 4.0 * step_size;
492                     integral_value +=
493                         32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size));
494                 }
495 
496                 evaluation_point = interval[0] - 2.0 * step_size;
497 
498                 for (i = 0; i < number_of_intervals; i++) {
499                     evaluation_point += 4.0 * step_size;
500                     integral_value += 12.0 * f(evaluation_point);
501                 }
502 
503                 integral_value *= (2.0 * step_size) / 45.0;
504         }
505         return integral_value;
506     },
507 
508     /**
509      * Calculates the integral of function f over interval using Romberg iteration.
510      * @param {Array} interval The integration interval, e.g. [0, 3].
511      * @param {function} f A function which takes one argument of type number and returns a number.
512      * @param {Object} [config] The algorithm setup. Accepted properties are max_iterations of type number and precision eps.
513      * @param {Number} [config.max_iterations=20]
514      * @param {Number} [config.eps=0.0000001]
515      * @returns {Number} Integral value of f over interval
516      * @example
517      * function f(x) {
518      *   return x*x;
519      * }
520      *
521      * // calculates integral of <tt>f</tt> from 0 to 2.
522      * var area1 = JXG.Math.Numerics.Romberg([0, 2], f);
523      *
524      * // the same with an anonymous function
525      * var area2 = JXG.Math.Numerics.Romberg([0, 2], function (x) { return x*x; });
526      *
527      * // use trapez rule with maximum of 16 iterations or stop if the precision 0.0001 has been reached.
528      * var area3 = JXG.Math.Numerics.Romberg([0, 2], f,
529      *                                   {max_iterations: 16, eps: 0.0001});
530      * @memberof JXG.Math.Numerics
531      */
532     Romberg: function (interval, f, config) {
533         var a,
534             b,
535             h,
536             s,
537             n,
538             k,
539             i,
540             q,
541             p = [],
542             integral = 0.0,
543             last = Infinity,
544             m = config && Type.isNumber(config.max_iterations) ? config.max_iterations : 20,
545             eps = config && Type.isNumber(config.eps) ? config.eps : config.eps || 0.0000001;
546 
547         a = interval[0];
548         b = interval[1];
549         h = b - a;
550         n = 1;
551 
552         p[0] = 0.5 * h * (f(a) + f(b));
553 
554         for (k = 0; k < m; ++k) {
555             s = 0;
556             h *= 0.5;
557             n *= 2;
558             q = 1;
559 
560             for (i = 1; i < n; i += 2) {
561                 s += f(a + i * h);
562             }
563 
564             p[k + 1] = 0.5 * p[k] + s * h;
565 
566             integral = p[k + 1];
567             for (i = k - 1; i >= 0; --i) {
568                 q *= 4;
569                 p[i] = p[i + 1] + (p[i + 1] - p[i]) / (q - 1.0);
570                 integral = p[i];
571             }
572 
573             if (Math.abs(integral - last) < eps * Math.abs(integral)) {
574                 break;
575             }
576             last = integral;
577         }
578 
579         return integral;
580     },
581 
582     /**
583      * Calculates the integral of function f over interval using Gauss-Legendre quadrature.
584      * @param {Array} interval The integration interval, e.g. [0, 3].
585      * @param {function} f A function which takes one argument of type number and returns a number.
586      * @param {Object} [config] The algorithm setup. Accepted property is the order n of type number. n is allowed to take
587      * values between 2 and 18, default value is 12.
588      * @param {Number} [config.n=16]
589      * @returns {Number} Integral value of f over interval
590      * @example
591      * function f(x) {
592      *   return x*x;
593      * }
594      *
595      * // calculates integral of <tt>f</tt> from 0 to 2.
596      * var area1 = JXG.Math.Numerics.GaussLegendre([0, 2], f);
597      *
598      * // the same with an anonymous function
599      * var area2 = JXG.Math.Numerics.GaussLegendre([0, 2], function (x) { return x*x; });
600      *
601      * // use 16 point Gauss-Legendre rule.
602      * var area3 = JXG.Math.Numerics.GaussLegendre([0, 2], f,
603      *                                   {n: 16});
604      * @memberof JXG.Math.Numerics
605      */
606     GaussLegendre: function (interval, f, config) {
607         var a,
608             b,
609             i,
610             m,
611             xp,
612             xm,
613             result = 0.0,
614             table_xi = [],
615             table_w = [],
616             xi,
617             w,
618             n = config && Type.isNumber(config.n) ? config.n : 12;
619 
620         if (n > 18) {
621             n = 18;
622         }
623 
624         /* n = 2 */
625         table_xi[2] = [0.5773502691896257645091488];
626         table_w[2] = [1.0];
627 
628         /* n = 4 */
629         table_xi[4] = [0.3399810435848562648026658, 0.8611363115940525752239465];
630         table_w[4] = [0.6521451548625461426269361, 0.3478548451374538573730639];
631 
632         /* n = 6 */
633         table_xi[6] = [
634             0.2386191860831969086305017, 0.6612093864662645136613996,
635             0.9324695142031520278123016
636         ];
637         table_w[6] = [
638             0.4679139345726910473898703, 0.3607615730481386075698335,
639             0.1713244923791703450402961
640         ];
641 
642         /* n = 8 */
643         table_xi[8] = [
644             0.1834346424956498049394761, 0.525532409916328985817739,
645             0.7966664774136267395915539, 0.9602898564975362316835609
646         ];
647         table_w[8] = [
648             0.3626837833783619829651504, 0.3137066458778872873379622,
649             0.222381034453374470544356, 0.1012285362903762591525314
650         ];
651 
652         /* n = 10 */
653         table_xi[10] = [
654             0.148874338981631210884826, 0.4333953941292471907992659,
655             0.6794095682990244062343274, 0.8650633666889845107320967, 0.973906528517171720077964
656         ];
657         table_w[10] = [
658             0.295524224714752870173893, 0.2692667193099963550912269,
659             0.2190863625159820439955349, 0.1494513491505805931457763,
660             0.0666713443086881375935688
661         ];
662 
663         /* n = 12 */
664         table_xi[12] = [
665             0.1252334085114689154724414, 0.3678314989981801937526915,
666             0.5873179542866174472967024, 0.7699026741943046870368938,
667             0.9041172563704748566784659, 0.9815606342467192506905491
668         ];
669         table_w[12] = [
670             0.2491470458134027850005624, 0.2334925365383548087608499,
671             0.2031674267230659217490645, 0.1600783285433462263346525,
672             0.1069393259953184309602547, 0.047175336386511827194616
673         ];
674 
675         /* n = 14 */
676         table_xi[14] = [
677             0.1080549487073436620662447, 0.3191123689278897604356718,
678             0.5152486363581540919652907, 0.6872929048116854701480198,
679             0.8272013150697649931897947, 0.9284348836635735173363911,
680             0.9862838086968123388415973
681         ];
682         table_w[14] = [
683             0.2152638534631577901958764, 0.2051984637212956039659241,
684             0.1855383974779378137417166, 0.1572031671581935345696019,
685             0.1215185706879031846894148, 0.0801580871597602098056333,
686             0.0351194603317518630318329
687         ];
688 
689         /* n = 16 */
690         table_xi[16] = [
691             0.0950125098376374401853193, 0.2816035507792589132304605,
692             0.4580167776572273863424194, 0.6178762444026437484466718,
693             0.7554044083550030338951012, 0.8656312023878317438804679,
694             0.9445750230732325760779884, 0.9894009349916499325961542
695         ];
696         table_w[16] = [
697             0.1894506104550684962853967, 0.1826034150449235888667637,
698             0.1691565193950025381893121, 0.1495959888165767320815017,
699             0.1246289712555338720524763, 0.0951585116824927848099251,
700             0.0622535239386478928628438, 0.0271524594117540948517806
701         ];
702 
703         /* n = 18 */
704         table_xi[18] = [
705             0.0847750130417353012422619, 0.2518862256915055095889729,
706             0.4117511614628426460359318, 0.5597708310739475346078715,
707             0.6916870430603532078748911, 0.8037049589725231156824175,
708             0.8926024664975557392060606, 0.9558239495713977551811959, 0.991565168420930946730016
709         ];
710         table_w[18] = [
711             0.1691423829631435918406565, 0.1642764837458327229860538,
712             0.154684675126265244925418, 0.1406429146706506512047313,
713             0.1225552067114784601845191, 0.100942044106287165562814,
714             0.0764257302548890565291297, 0.0497145488949697964533349,
715             0.0216160135264833103133427
716         ];
717 
718         /* n = 3 */
719         table_xi[3] = [0.0, 0.7745966692414833770358531];
720         table_w[3] = [0.8888888888888888888888889, 0.5555555555555555555555556];
721 
722         /* n = 5 */
723         table_xi[5] = [0.0, 0.5384693101056830910363144, 0.9061798459386639927976269];
724         table_w[5] = [
725             0.5688888888888888888888889, 0.4786286704993664680412915, 0.236926885056189087514264
726         ];
727 
728         /* n = 7 */
729         table_xi[7] = [
730             0.0, 0.4058451513773971669066064, 0.7415311855993944398638648,
731             0.9491079123427585245261897
732         ];
733         table_w[7] = [
734             0.417959183673469387755102, 0.3818300505051189449503698,
735             0.2797053914892766679014678, 0.1294849661688696932706114
736         ];
737 
738         /* n = 9 */
739         table_xi[9] = [
740             0.0, 0.324253423403808929038538, 0.613371432700590397308702,
741             0.8360311073266357942994298, 0.9681602395076260898355762
742         ];
743         table_w[9] = [
744             0.3302393550012597631645251, 0.3123470770400028400686304,
745             0.2606106964029354623187429, 0.180648160694857404058472, 0.0812743883615744119718922
746         ];
747 
748         /* n = 11 */
749         table_xi[11] = [
750             0.0, 0.269543155952344972331532, 0.5190961292068118159257257,
751             0.7301520055740493240934163, 0.8870625997680952990751578, 0.978228658146056992803938
752         ];
753         table_w[11] = [
754             0.2729250867779006307144835, 0.2628045445102466621806889,
755             0.2331937645919904799185237, 0.1862902109277342514260976,
756             0.1255803694649046246346943, 0.0556685671161736664827537
757         ];
758 
759         /* n = 13 */
760         table_xi[13] = [
761             0.0, 0.2304583159551347940655281, 0.4484927510364468528779129,
762             0.6423493394403402206439846, 0.8015780907333099127942065,
763             0.9175983992229779652065478, 0.9841830547185881494728294
764         ];
765         table_w[13] = [
766             0.2325515532308739101945895, 0.2262831802628972384120902,
767             0.2078160475368885023125232, 0.1781459807619457382800467,
768             0.1388735102197872384636018, 0.0921214998377284479144218,
769             0.0404840047653158795200216
770         ];
771 
772         /* n = 15 */
773         table_xi[15] = [
774             0.0, 0.2011940939974345223006283, 0.3941513470775633698972074,
775             0.5709721726085388475372267, 0.7244177313601700474161861,
776             0.8482065834104272162006483, 0.9372733924007059043077589,
777             0.9879925180204854284895657
778         ];
779         table_w[15] = [
780             0.2025782419255612728806202, 0.1984314853271115764561183,
781             0.1861610000155622110268006, 0.1662692058169939335532009,
782             0.1395706779261543144478048, 0.1071592204671719350118695,
783             0.0703660474881081247092674, 0.0307532419961172683546284
784         ];
785 
786         /* n = 17 */
787         table_xi[17] = [
788             0.0, 0.1784841814958478558506775, 0.3512317634538763152971855,
789             0.5126905370864769678862466, 0.6576711592166907658503022,
790             0.7815140038968014069252301, 0.8802391537269859021229557,
791             0.950675521768767761222717, 0.990575475314417335675434
792         ];
793         table_w[17] = [
794             0.1794464703562065254582656, 0.176562705366992646325271,
795             0.1680041021564500445099707, 0.1540457610768102880814316, 0.13513636846852547328632,
796             0.1118838471934039710947884, 0.0850361483171791808835354,
797             0.0554595293739872011294402, 0.02414830286854793196011
798         ];
799 
800         a = interval[0];
801         b = interval[1];
802 
803         //m = Math.ceil(n * 0.5);
804         m = (n + 1) >> 1;
805 
806         xi = table_xi[n];
807         w = table_w[n];
808 
809         xm = 0.5 * (b - a);
810         xp = 0.5 * (b + a);
811 
812         if (n & (1 === 1)) {
813             // n odd
814             result = w[0] * f(xp);
815             for (i = 1; i < m; ++i) {
816                 result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i]));
817             }
818         } else {
819             // n even
820             result = 0.0;
821             for (i = 0; i < m; ++i) {
822                 result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i]));
823             }
824         }
825 
826         return xm * result;
827     },
828 
829     /**
830      * Scale error in Gauss Kronrod quadrature.
831      * Internal method used in {@link JXG.Math.Numerics._gaussKronrod}.
832      * @private
833      */
834     _rescale_error: function (err, result_abs, result_asc) {
835         var scale,
836             min_err,
837             DBL_MIN = 2.2250738585072014e-308,
838             DBL_EPS = 2.2204460492503131e-16;
839 
840         err = Math.abs(err);
841         if (result_asc !== 0 && err !== 0) {
842             scale = Math.pow((200 * err) / result_asc, 1.5);
843 
844             if (scale < 1.0) {
845                 err = result_asc * scale;
846             } else {
847                 err = result_asc;
848             }
849         }
850         if (result_abs > DBL_MIN / (50 * DBL_EPS)) {
851             min_err = 50 * DBL_EPS * result_abs;
852 
853             if (min_err > err) {
854                 err = min_err;
855             }
856         }
857 
858         return err;
859     },
860 
861     /**
862      * Generic Gauss-Kronrod quadrature algorithm.
863      * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15},
864      * {@link JXG.Math.Numerics.GaussKronrod21},
865      * {@link JXG.Math.Numerics.GaussKronrod31}.
866      * Taken from QUADPACK.
867      *
868      * @param {Array} interval The integration interval, e.g. [0, 3].
869      * @param {function} f A function which takes one argument of type number and returns a number.
870      * @param {Number} n order of approximation. Actually, n is the length of the array xgk. For example, for the 15-point Kronrod rule, n is equal to 8.
871      * @param {Array} xgk Kronrod quadrature abscissae
872      * @param {Array} wg Weights of the Gauss rule
873      * @param {Array} wgk Weights of the Kronrod rule
874      * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc.
875      * See the library QUADPACK for an explanation.
876      *
877      * @returns {Number} Integral value of f over interval
878      *
879      * @private
880      */
881     _gaussKronrod: function (interval, f, n, xgk, wg, wgk, resultObj) {
882         var a = interval[0],
883             b = interval[1],
884             up,
885             result,
886             center = 0.5 * (a + b),
887             half_length = 0.5 * (b - a),
888             abs_half_length = Math.abs(half_length),
889             f_center = f(center),
890             result_gauss = 0.0,
891             result_kronrod = f_center * wgk[n - 1],
892             result_abs = Math.abs(result_kronrod),
893             result_asc = 0.0,
894             mean = 0.0,
895             err = 0.0,
896             j, jtw, jtwm1,
897             abscissa,
898             fval1, fval2, fsum,
899             fv1 = [],
900             fv2 = [];
901 
902         if (n % 2 === 0) {
903             result_gauss = f_center * wg[n / 2 - 1];
904         }
905 
906         up = Math.floor((n - 1) / 2);
907         for (j = 0; j < up; j++) {
908             jtw = j * 2 + 1; // in original fortran j=1,2,3 jtw=2,4,6
909             abscissa = half_length * xgk[jtw];
910             fval1 = f(center - abscissa);
911             fval2 = f(center + abscissa);
912             fsum = fval1 + fval2;
913             fv1[jtw] = fval1;
914             fv2[jtw] = fval2;
915             result_gauss += wg[j] * fsum;
916             result_kronrod += wgk[jtw] * fsum;
917             result_abs += wgk[jtw] * (Math.abs(fval1) + Math.abs(fval2));
918         }
919 
920         up = Math.floor(n / 2);
921         for (j = 0; j < up; j++) {
922             jtwm1 = j * 2;
923             abscissa = half_length * xgk[jtwm1];
924             fval1 = f(center - abscissa);
925             fval2 = f(center + abscissa);
926             fv1[jtwm1] = fval1;
927             fv2[jtwm1] = fval2;
928             result_kronrod += wgk[jtwm1] * (fval1 + fval2);
929             result_abs += wgk[jtwm1] * (Math.abs(fval1) + Math.abs(fval2));
930         }
931 
932         mean = result_kronrod * 0.5;
933         result_asc = wgk[n - 1] * Math.abs(f_center - mean);
934 
935         for (j = 0; j < n - 1; j++) {
936             result_asc += wgk[j] * (Math.abs(fv1[j] - mean) + Math.abs(fv2[j] - mean));
937         }
938 
939         // scale by the width of the integration region
940         err = (result_kronrod - result_gauss) * half_length;
941 
942         result_kronrod *= half_length;
943         result_abs *= abs_half_length;
944         result_asc *= abs_half_length;
945         result = result_kronrod;
946 
947         resultObj.abserr = this._rescale_error(err, result_abs, result_asc);
948         resultObj.resabs = result_abs;
949         resultObj.resasc = result_asc;
950 
951         return result;
952     },
953 
954     /**
955      * 15-point Gauss-Kronrod quadrature algorithm, see the library QUADPACK
956      * @param {Array} interval The integration interval, e.g. [0, 3].
957      * @param {function} f A function which takes one argument of type number and returns a number.
958      * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library
959      *  QUADPACK for an explanation.
960      *
961      * @returns {Number} Integral value of f over interval
962      *
963      * @memberof JXG.Math.Numerics
964      */
965     GaussKronrod15: function (interval, f, resultObj) {
966         /* Gauss quadrature weights and kronrod quadrature abscissae and
967                 weights as evaluated with 80 decimal digit arithmetic by
968                 L. W. Fullerton, Bell Labs, Nov. 1981. */
969 
970         var xgk =
971             /* abscissae of the 15-point kronrod rule */
972             [
973                 0.991455371120812639206854697526329, 0.949107912342758524526189684047851,
974                 0.864864423359769072789712788640926, 0.741531185599394439863864773280788,
975                 0.58608723546769113029414483825873, 0.405845151377397166906606412076961,
976                 0.207784955007898467600689403773245, 0.0
977             ],
978             /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
979                 xgk[0], xgk[2], ... abscissae to optimally extend the 7-point gauss rule */
980 
981             wg =
982                 /* weights of the 7-point gauss rule */
983                 [
984                     0.129484966168869693270611432679082, 0.27970539148927666790146777142378,
985                     0.381830050505118944950369775488975, 0.417959183673469387755102040816327
986                 ],
987             wgk =
988                 /* weights of the 15-point kronrod rule */
989                 [
990                     0.02293532201052922496373200805897, 0.063092092629978553290700663189204,
991                     0.104790010322250183839876322541518, 0.140653259715525918745189590510238,
992                     0.16900472663926790282658342659855, 0.190350578064785409913256402421014,
993                     0.204432940075298892414161999234649, 0.209482141084727828012999174891714
994                 ];
995 
996         return this._gaussKronrod(interval, f, 8, xgk, wg, wgk, resultObj);
997     },
998 
999     /**
1000      * 21 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK
1001      * @param {Array} interval The integration interval, e.g. [0, 3].
1002      * @param {function} f A function which takes one argument of type number and returns a number.
1003      * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library
1004      *  QUADPACK for an explanation.
1005      *
1006      * @returns {Number} Integral value of f over interval
1007      *
1008      * @memberof JXG.Math.Numerics
1009      */
1010     GaussKronrod21: function (interval, f, resultObj) {
1011         /* Gauss quadrature weights and kronrod quadrature abscissae and
1012                 weights as evaluated with 80 decimal digit arithmetic by
1013                 L. W. Fullerton, Bell Labs, Nov. 1981. */
1014 
1015         var xgk =
1016             /* abscissae of the 21-point kronrod rule */
1017             [
1018                 0.995657163025808080735527280689003, 0.973906528517171720077964012084452,
1019                 0.930157491355708226001207180059508, 0.865063366688984510732096688423493,
1020                 0.780817726586416897063717578345042, 0.679409568299024406234327365114874,
1021                 0.562757134668604683339000099272694, 0.433395394129247190799265943165784,
1022                 0.294392862701460198131126603103866, 0.14887433898163121088482600112972, 0.0
1023             ],
1024             /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule.
1025                 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */
1026             wg =
1027                 /* weights of the 10-point gauss rule */
1028                 [
1029                     0.066671344308688137593568809893332, 0.149451349150580593145776339657697,
1030                     0.219086362515982043995534934228163, 0.269266719309996355091226921569469,
1031                     0.295524224714752870173892994651338
1032                 ],
1033             wgk =
1034                 /* weights of the 21-point kronrod rule */
1035                 [
1036                     0.011694638867371874278064396062192, 0.03255816230796472747881897245939,
1037                     0.05475589657435199603138130024458, 0.07503967481091995276704314091619,
1038                     0.093125454583697605535065465083366, 0.109387158802297641899210590325805,
1039                     0.123491976262065851077958109831074, 0.134709217311473325928054001771707,
1040                     0.142775938577060080797094273138717, 0.147739104901338491374841515972068,
1041                     0.149445554002916905664936468389821
1042                 ];
1043 
1044         return this._gaussKronrod(interval, f, 11, xgk, wg, wgk, resultObj);
1045     },
1046 
1047     /**
1048      * 31 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK
1049      * @param {Array} interval The integration interval, e.g. [0, 3].
1050      * @param {function} f A function which takes one argument of type number and returns a number.
1051      * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library
1052      *  QUADPACK for an explanation.
1053      *
1054      * @returns {Number} Integral value of f over interval
1055      *
1056      * @memberof JXG.Math.Numerics
1057      */
1058     GaussKronrod31: function (interval, f, resultObj) {
1059         /* Gauss quadrature weights and kronrod quadrature abscissae and
1060                 weights as evaluated with 80 decimal digit arithmetic by
1061                 L. W. Fullerton, Bell Labs, Nov. 1981. */
1062 
1063         var xgk =
1064             /* abscissae of the 21-point kronrod rule */
1065             [
1066                 0.998002298693397060285172840152271, 0.987992518020485428489565718586613,
1067                 0.967739075679139134257347978784337, 0.937273392400705904307758947710209,
1068                 0.897264532344081900882509656454496, 0.848206583410427216200648320774217,
1069                 0.790418501442465932967649294817947, 0.724417731360170047416186054613938,
1070                 0.650996741297416970533735895313275, 0.570972172608538847537226737253911,
1071                 0.485081863640239680693655740232351, 0.394151347077563369897207370981045,
1072                 0.299180007153168812166780024266389, 0.201194093997434522300628303394596,
1073                 0.101142066918717499027074231447392, 0.0
1074             ],
1075             /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule.
1076                 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */
1077             wg =
1078                 /* weights of the 10-point gauss rule */
1079                 [
1080                     0.030753241996117268354628393577204, 0.070366047488108124709267416450667,
1081                     0.107159220467171935011869546685869, 0.139570677926154314447804794511028,
1082                     0.166269205816993933553200860481209, 0.186161000015562211026800561866423,
1083                     0.198431485327111576456118326443839, 0.202578241925561272880620199967519
1084                 ],
1085             wgk =
1086                 /* weights of the 21-point kronrod rule */
1087                 [
1088                     0.005377479872923348987792051430128, 0.015007947329316122538374763075807,
1089                     0.025460847326715320186874001019653, 0.03534636079137584622203794847836,
1090                     0.04458975132476487660822729937328, 0.05348152469092808726534314723943,
1091                     0.062009567800670640285139230960803, 0.069854121318728258709520077099147,
1092                     0.076849680757720378894432777482659, 0.083080502823133021038289247286104,
1093                     0.088564443056211770647275443693774, 0.093126598170825321225486872747346,
1094                     0.096642726983623678505179907627589, 0.099173598721791959332393173484603,
1095                     0.10076984552387559504494666261757, 0.101330007014791549017374792767493
1096                 ];
1097 
1098         return this._gaussKronrod(interval, f, 16, xgk, wg, wgk, resultObj);
1099     },
1100 
1101     /**
1102      * Generate workspace object for {@link JXG.Math.Numerics.Qag}.
1103      * @param {Array} interval The integration interval, e.g. [0, 3].
1104      * @param {Number} n Max. limit
1105      * @returns {Object} Workspace object
1106      *
1107      * @private
1108      * @memberof JXG.Math.Numerics
1109      */
1110     _workspace: function (interval, n) {
1111         return {
1112             limit: n,
1113             size: 0,
1114             nrmax: 0,
1115             i: 0,
1116             alist: [interval[0]],
1117             blist: [interval[1]],
1118             rlist: [0.0],
1119             elist: [0.0],
1120             order: [0],
1121             level: [0],
1122 
1123             qpsrt: function () {
1124                 var last = this.size - 1,
1125                     limit = this.limit,
1126                     errmax,
1127                     errmin,
1128                     i,
1129                     k,
1130                     top,
1131                     i_nrmax = this.nrmax,
1132                     i_maxerr = this.order[i_nrmax];
1133 
1134                 /* Check whether the list contains more than two error estimates */
1135                 if (last < 2) {
1136                     this.order[0] = 0;
1137                     this.order[1] = 1;
1138                     this.i = i_maxerr;
1139                     return;
1140                 }
1141 
1142                 errmax = this.elist[i_maxerr];
1143 
1144                 /* This part of the routine is only executed if, due to a difficult
1145                         integrand, subdivision increased the error estimate. In the normal
1146                         case the insert procedure should start after the nrmax-th largest
1147                         error estimate. */
1148                 while (i_nrmax > 0 && errmax > this.elist[this.order[i_nrmax - 1]]) {
1149                     this.order[i_nrmax] = this.order[i_nrmax - 1];
1150                     i_nrmax--;
1151                 }
1152 
1153                 /* Compute the number of elements in the list to be maintained in
1154                         descending order. This number depends on the number of
1155                         subdivisions still allowed. */
1156                 if (last < limit / 2 + 2) {
1157                     top = last;
1158                 } else {
1159                     top = limit - last + 1;
1160                 }
1161 
1162                 /* Insert errmax by traversing the list top-down, starting
1163                         comparison from the element elist(order(i_nrmax+1)). */
1164                 i = i_nrmax + 1;
1165 
1166                 /* The order of the tests in the following line is important to
1167                         prevent a segmentation fault */
1168                 while (i < top && errmax < this.elist[this.order[i]]) {
1169                     this.order[i - 1] = this.order[i];
1170                     i++;
1171                 }
1172 
1173                 this.order[i - 1] = i_maxerr;
1174 
1175                 /* Insert errmin by traversing the list bottom-up */
1176                 errmin = this.elist[last];
1177                 k = top - 1;
1178 
1179                 while (k > i - 2 && errmin >= this.elist[this.order[k]]) {
1180                     this.order[k + 1] = this.order[k];
1181                     k--;
1182                 }
1183 
1184                 this.order[k + 1] = last;
1185 
1186                 /* Set i_max and e_max */
1187                 i_maxerr = this.order[i_nrmax];
1188                 this.i = i_maxerr;
1189                 this.nrmax = i_nrmax;
1190             },
1191 
1192             set_initial_result: function (result, error) {
1193                 this.size = 1;
1194                 this.rlist[0] = result;
1195                 this.elist[0] = error;
1196             },
1197 
1198             update: function (a1, b1, area1, error1, a2, b2, area2, error2) {
1199                 var i_max = this.i,
1200                     i_new = this.size,
1201                     new_level = this.level[this.i] + 1;
1202 
1203                 /* append the newly-created intervals to the list */
1204 
1205                 if (error2 > error1) {
1206                     this.alist[i_max] = a2; /* blist[maxerr] is already == b2 */
1207                     this.rlist[i_max] = area2;
1208                     this.elist[i_max] = error2;
1209                     this.level[i_max] = new_level;
1210 
1211                     this.alist[i_new] = a1;
1212                     this.blist[i_new] = b1;
1213                     this.rlist[i_new] = area1;
1214                     this.elist[i_new] = error1;
1215                     this.level[i_new] = new_level;
1216                 } else {
1217                     this.blist[i_max] = b1; /* alist[maxerr] is already == a1 */
1218                     this.rlist[i_max] = area1;
1219                     this.elist[i_max] = error1;
1220                     this.level[i_max] = new_level;
1221 
1222                     this.alist[i_new] = a2;
1223                     this.blist[i_new] = b2;
1224                     this.rlist[i_new] = area2;
1225                     this.elist[i_new] = error2;
1226                     this.level[i_new] = new_level;
1227                 }
1228 
1229                 this.size++;
1230 
1231                 if (new_level > this.maximum_level) {
1232                     this.maximum_level = new_level;
1233                 }
1234 
1235                 this.qpsrt();
1236             },
1237 
1238             retrieve: function () {
1239                 var i = this.i;
1240                 return {
1241                     a: this.alist[i],
1242                     b: this.blist[i],
1243                     r: this.rlist[i],
1244                     e: this.elist[i]
1245                 };
1246             },
1247 
1248             sum_results: function () {
1249                 var nn = this.size,
1250                     k,
1251                     result_sum = 0.0;
1252 
1253                 for (k = 0; k < nn; k++) {
1254                     result_sum += this.rlist[k];
1255                 }
1256 
1257                 return result_sum;
1258             },
1259 
1260             subinterval_too_small: function (a1, a2, b2) {
1261                 var e = 2.2204460492503131e-16,
1262                     u = 2.2250738585072014e-308,
1263                     tmp = (1 + 100 * e) * (Math.abs(a2) + 1000 * u);
1264 
1265                 return Math.abs(a1) <= tmp && Math.abs(b2) <= tmp;
1266             }
1267         };
1268     },
1269 
1270     /**
1271      * Quadrature algorithm qag from QUADPACK.
1272      * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15},
1273      * {@link JXG.Math.Numerics.GaussKronrod21},
1274      * {@link JXG.Math.Numerics.GaussKronrod31}.
1275      *
1276      * @param {Array} interval The integration interval, e.g. [0, 3].
1277      * @param {function} f A function which takes one argument of type number and returns a number.
1278      * @param {Object} [config] The algorithm setup. Accepted propert are max. recursion limit of type number,
1279      * and epsrel and epsabs, the relative and absolute required precision of type number. Further,
1280      * q the internal quadrature sub-algorithm of type function.
1281      * @param {Number} [config.limit=15]
1282      * @param {Number} [config.epsrel=0.0000001]
1283      * @param {Number} [config.epsabs=0.0000001]
1284      * @param {Number} [config.q=JXG.Math.Numerics.GaussKronrod15]
1285      * @returns {Number} Integral value of f over interval
1286      *
1287      * @example
1288      * function f(x) {
1289      *   return x*x;
1290      * }
1291      *
1292      * // calculates integral of <tt>f</tt> from 0 to 2.
1293      * var area1 = JXG.Math.Numerics.Qag([0, 2], f);
1294      *
1295      * // the same with an anonymous function
1296      * var area2 = JXG.Math.Numerics.Qag([0, 2], function (x) { return x*x; });
1297      *
1298      * // use JXG.Math.Numerics.GaussKronrod31 rule as sub-algorithm.
1299      * var area3 = JXG.Math.Numerics.Quag([0, 2], f,
1300      *                                   {q: JXG.Math.Numerics.GaussKronrod31});
1301      * @memberof JXG.Math.Numerics
1302      */
1303     Qag: function (interval, f, config) {
1304         var DBL_EPS = 2.2204460492503131e-16,
1305             ws = this._workspace(interval, 1000),
1306             limit = config && Type.isNumber(config.limit) ? config.limit : 15,
1307             epsrel = config && Type.isNumber(config.epsrel) ? config.epsrel : 0.0000001,
1308             epsabs = config && Type.isNumber(config.epsabs) ? config.epsabs : 0.0000001,
1309             q = config && Type.isFunction(config.q) ? config.q : this.GaussKronrod15,
1310             resultObj = {},
1311             area,
1312             errsum,
1313             result0,
1314             abserr0,
1315             resabs0,
1316             resasc0,
1317             result,
1318             tolerance,
1319             iteration = 0,
1320             roundoff_type1 = 0,
1321             roundoff_type2 = 0,
1322             error_type = 0,
1323             round_off,
1324             a1,
1325             b1,
1326             a2,
1327             b2,
1328             a_i,
1329             b_i,
1330             r_i,
1331             e_i,
1332             area1 = 0,
1333             area2 = 0,
1334             area12 = 0,
1335             error1 = 0,
1336             error2 = 0,
1337             error12 = 0,
1338             resasc1,
1339             resasc2,
1340             // resabs1, resabs2,
1341             wsObj,
1342             delta;
1343 
1344         if (limit > ws.limit) {
1345             JXG.warn("iteration limit exceeds available workspace");
1346         }
1347         if (epsabs <= 0 && (epsrel < 50 * Mat.eps || epsrel < 0.5e-28)) {
1348             JXG.warn("tolerance cannot be acheived with given epsabs and epsrel");
1349         }
1350 
1351         result0 = q.apply(this, [interval, f, resultObj]);
1352         abserr0 = resultObj.abserr;
1353         resabs0 = resultObj.resabs;
1354         resasc0 = resultObj.resasc;
1355 
1356         ws.set_initial_result(result0, abserr0);
1357         tolerance = Math.max(epsabs, epsrel * Math.abs(result0));
1358         round_off = 50 * DBL_EPS * resabs0;
1359 
1360         if (abserr0 <= round_off && abserr0 > tolerance) {
1361             result = result0;
1362             // abserr = abserr0;
1363 
1364             JXG.warn("cannot reach tolerance because of roundoff error on first attempt");
1365             return -Infinity;
1366         }
1367 
1368         if ((abserr0 <= tolerance && abserr0 !== resasc0) || abserr0 === 0.0) {
1369             result = result0;
1370             // abserr = abserr0;
1371 
1372             return result;
1373         }
1374 
1375         if (limit === 1) {
1376             result = result0;
1377             // abserr = abserr0;
1378 
1379             JXG.warn("a maximum of one iteration was insufficient");
1380             return -Infinity;
1381         }
1382 
1383         area = result0;
1384         errsum = abserr0;
1385         iteration = 1;
1386 
1387         do {
1388             area1 = 0;
1389             area2 = 0;
1390             area12 = 0;
1391             error1 = 0;
1392             error2 = 0;
1393             error12 = 0;
1394 
1395             /* Bisect the subinterval with the largest error estimate */
1396             wsObj = ws.retrieve();
1397             a_i = wsObj.a;
1398             b_i = wsObj.b;
1399             r_i = wsObj.r;
1400             e_i = wsObj.e;
1401 
1402             a1 = a_i;
1403             b1 = 0.5 * (a_i + b_i);
1404             a2 = b1;
1405             b2 = b_i;
1406 
1407             area1 = q.apply(this, [[a1, b1], f, resultObj]);
1408             error1 = resultObj.abserr;
1409             // resabs1 = resultObj.resabs;
1410             resasc1 = resultObj.resasc;
1411 
1412             area2 = q.apply(this, [[a2, b2], f, resultObj]);
1413             error2 = resultObj.abserr;
1414             // resabs2 = resultObj.resabs;
1415             resasc2 = resultObj.resasc;
1416 
1417             area12 = area1 + area2;
1418             error12 = error1 + error2;
1419 
1420             errsum += error12 - e_i;
1421             area += area12 - r_i;
1422 
1423             if (resasc1 !== error1 && resasc2 !== error2) {
1424                 delta = r_i - area12;
1425                 if (Math.abs(delta) <= 1.0e-5 * Math.abs(area12) && error12 >= 0.99 * e_i) {
1426                     roundoff_type1++;
1427                 }
1428                 if (iteration >= 10 && error12 > e_i) {
1429                     roundoff_type2++;
1430                 }
1431             }
1432 
1433             tolerance = Math.max(epsabs, epsrel * Math.abs(area));
1434 
1435             if (errsum > tolerance) {
1436                 if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
1437                     error_type = 2; /* round off error */
1438                 }
1439 
1440                 /* set error flag in the case of bad integrand behaviour at
1441                     a point of the integration range */
1442 
1443                 if (ws.subinterval_too_small(a1, a2, b2)) {
1444                     error_type = 3;
1445                 }
1446             }
1447 
1448             ws.update(a1, b1, area1, error1, a2, b2, area2, error2);
1449             wsObj = ws.retrieve();
1450             a_i = wsObj.a_i;
1451             b_i = wsObj.b_i;
1452             r_i = wsObj.r_i;
1453             e_i = wsObj.e_i;
1454 
1455             iteration++;
1456         } while (iteration < limit && !error_type && errsum > tolerance);
1457 
1458         result = ws.sum_results();
1459         // abserr = errsum;
1460         /*
1461   if (errsum <= tolerance)
1462     {
1463       return GSL_SUCCESS;
1464     }
1465   else if (error_type == 2)
1466     {
1467       GSL_ERROR ("roundoff error prevents tolerance from being achieved",
1468                  GSL_EROUND);
1469     }
1470   else if (error_type == 3)
1471     {
1472       GSL_ERROR ("bad integrand behavior found in the integration interval",
1473                  GSL_ESING);
1474     }
1475   else if (iteration == limit)
1476     {
1477       GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
1478     }
1479   else
1480     {
1481       GSL_ERROR ("could not integrate function", GSL_EFAILED);
1482     }
1483 */
1484 
1485         return result;
1486     },
1487 
1488     /**
1489      * Integral of function f over interval.
1490      * @param {Array} interval The integration interval, e.g. [0, 3].
1491      * @param {function} f A function which takes one argument of type number and returns a number.
1492      * @returns {Number} The value of the integral of f over interval
1493      * @see JXG.Math.Numerics.NewtonCotes
1494      * @see JXG.Math.Numerics.Romberg
1495      * @see JXG.Math.Numerics.Qag
1496      * @memberof JXG.Math.Numerics
1497      */
1498     I: function (interval, f) {
1499         // return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'});
1500         // return this.Romberg(interval, f, {max_iterations: 20, eps: 0.0000001});
1501         return this.Qag(interval, f, {
1502             q: this.GaussKronrod15,
1503             limit: 15,
1504             epsrel: 0.0000001,
1505             epsabs: 0.0000001
1506         });
1507     },
1508 
1509     /**
1510      * Newton's method to find roots of a funtion in one variable.
1511      * @param {function} f We search for a solution of f(x)=0.
1512      * @param {Number} x initial guess for the root, i.e. start value.
1513      * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
1514      * the function is a method of an object and contains a reference to its parent object via "this".
1515      * @returns {Number} A root of the function f.
1516      * @memberof JXG.Math.Numerics
1517      */
1518     Newton: function (f, x, context) {
1519         var df,
1520             i = 0,
1521             h = Mat.eps,
1522             newf = f.apply(context, [x]);
1523         // nfev = 1;
1524 
1525         // For compatibility
1526         if (Type.isArray(x)) {
1527             x = x[0];
1528         }
1529 
1530         while (i < 50 && Math.abs(newf) > h) {
1531             df = this.D(f, context)(x);
1532             // nfev += 2;
1533 
1534             if (Math.abs(df) > h) {
1535                 x -= newf / df;
1536             } else {
1537                 x += Math.random() * 0.2 - 1.0;
1538             }
1539 
1540             newf = f.apply(context, [x]);
1541             // nfev += 1;
1542             i += 1;
1543         }
1544 
1545         return x;
1546     },
1547 
1548     /**
1549      * Abstract method to find roots of univariate functions, which - for the time being -
1550      * is an alias for {@link JXG.Math.Numerics.chandrupatla}.
1551      * @param {function} f We search for a solution of f(x)=0.
1552      * @param {Number|Array} x initial guess for the root, i.e. starting value, or start interval enclosing the root.
1553      * If x is an interval [a,b], it is required that f(a)f(b) <= 0, otherwise the minimum of f in [a, b] will be returned.
1554      * If x is a number, the algorithms tries to enclose the root by an interval [a, b] containing x and the root and
1555      * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method.
1556      *
1557      * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
1558      * the function is a method of an object and contains a reference to its parent object via "this".
1559      * @returns {Number} A root of the function f.
1560      *
1561      * @see JXG.Math.Numerics.chandrupatla
1562      * @see JXG.Math.Numerics.fzero
1563      * @see JXG.Math.Numerics.polzeros
1564      * @see JXG.Math.Numerics.Newton
1565      * @memberof JXG.Math.Numerics
1566      */
1567     root: function (f, x, context) {
1568         //return this.fzero(f, x, context);
1569         return this.chandrupatla(f, x, context);
1570     },
1571 
1572     /**
1573      * Compute an intersection of the curves c1 and c2
1574      * with a generalized Newton method (Newton-Raphson).
1575      * We want to find values t1, t2 such that
1576      * c1(t1) = c2(t2), i.e.
1577      * <br>
1578      * (c1_x(t1) - c2_x(t2), c1_y(t1) - c2_y(t2)) = (0, 0).
1579      * <p>
1580      * We set
1581      * (e, f) := (c1_x(t1) - c2_x(t2), c1_y(t1) - c2_y(t2))
1582      * <p>
1583      * The Jacobian J is defined by
1584      * <pre>
1585      * J = (a, b)
1586      *     (c, d)
1587      * </pre>
1588      * where
1589      * <ul>
1590      * <li> a = c1_x'(t1)
1591      * <li> b = -c2_x'(t2)
1592      * <li> c = c1_y'(t1)
1593      * <li> d = -c2_y'(t2)
1594      * </ul>
1595      * The inverse J^(-1) of J is equal to
1596      * <pre>
1597      *  (d, -b) / (ad - bc)
1598      *  (-c, a) / (ad - bc)
1599      * </pre>
1600      *
1601      * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f).
1602      * <p>
1603      *
1604      * @param {JXG.Curve} c1 Curve, Line or Circle
1605      * @param {JXG.Curve} c2 Curve, Line or Circle
1606      * @param {Number} t1ini start value for t1
1607      * @param {Number} t2ini start value for t2
1608      * @returns {JXG.Coords} intersection point
1609      * @memberof JXG.Math.Numerics
1610      */
1611     generalizedNewton: function (c1, c2, t1ini, t2ini) {
1612         var t1, t2,
1613             a, b, c, d, e, f,
1614             disc,
1615             F,
1616             D00, D01, D10, D11,
1617             count = 0;
1618 
1619         // if (this.generalizedNewton.t1memo) {
1620         //     t1 = this.generalizedNewton.t1memo;
1621         //     t2 = this.generalizedNewton.t2memo;
1622         // } else {
1623         t1 = t1ini;
1624         t2 = t2ini;
1625         // }
1626 
1627         e = c1.X(t1) - c2.X(t2);
1628         f = c1.Y(t1) - c2.Y(t2);
1629         F = e * e + f * f;
1630 
1631         D00 = this.D(c1.X, c1);
1632         D01 = this.D(c2.X, c2);
1633         D10 = this.D(c1.Y, c1);
1634         D11 = this.D(c2.Y, c2);
1635 
1636         while (F > Mat.eps && count < 10) {
1637             a = D00(t1);
1638             b = -D01(t2);
1639             c = D10(t1);
1640             d = -D11(t2);
1641             disc = a * d - b * c;
1642             t1 -= (d * e - b * f) / disc;
1643             t2 -= (a * f - c * e) / disc;
1644             e = c1.X(t1) - c2.X(t2);
1645             f = c1.Y(t1) - c2.Y(t2);
1646             F = e * e + f * f;
1647             count += 1;
1648         }
1649 
1650         // this.generalizedNewton.t1memo = t1;
1651         // this.generalizedNewton.t2memo = t2;
1652 
1653         if (Math.abs(t1) < Math.abs(t2)) {
1654             return [c1.X(t1), c1.Y(t1)];
1655         }
1656 
1657         return [c2.X(t2), c2.Y(t2)];
1658     },
1659 
1660     /**
1661      * Apply damped Newton-Raphson algorithm to determine the intersection
1662      * between the curve elements c1 and c2. Transformations of the curves
1663      * are already taken into regard.
1664      * <p>
1665      * We use a very high accuracy: Mat.eps**3
1666      *
1667      * @deprecated
1668      * @param {JXG.Curve} c1 Curve, Line or Circle
1669      * @param {JXG.Curve} c2 Curve, Line or Circle
1670      * @param {Number} t1ini Start value for curve c1
1671      * @param {Number} t2ini Start value for curve c2
1672      * @param {Number} gamma Damping factor, should be in the open interval (0, 1)
1673      * @param {Number} eps Stop if function value is smaller than eps
1674      * @returns {Array} [t1, t2, F2], where t1 and t2 are the parameters of the intersection for both curves, F2 is ||c1[t1]-c2[t2]||**2.
1675      */
1676     generalizedDampedNewtonCurves: function (c1, c2, t1ini, t2ini, gamma, eps) {
1677         var t1, t2,
1678             a, b, c, d, e, f,
1679             disc,
1680             F2,
1681             f1, f2,
1682             D, Dt,
1683             max_it = 40,
1684             count = 0;
1685 
1686         t1 = t1ini;
1687         t2 = t2ini;
1688 
1689         f1 = c1.Ft(t1);
1690         f2 = c2.Ft(t2);
1691         e = f1[1] - f2[1];
1692         f = f1[2] - f2[2];
1693         F2 = e * e + f * f;
1694 
1695         D = function(t1, t2) {
1696             var h = Mat.eps,
1697                 f1_1 = c1.Ft(t1 - h),
1698                 f1_2 = c1.Ft(t1 + h),
1699                 f2_1 = c2.Ft(t2 - h),
1700                 f2_2 = c2.Ft(t2 + h);
1701             return [
1702                 [ (f1_2[1] - f1_1[1]) / (2 * h),
1703                  -(f2_2[1] - f2_1[1]) / (2 * h)],
1704                 [ (f1_2[2] - f1_1[2]) / (2 * h),
1705                  -(f2_2[2] - f2_1[2]) / (2 * h)]
1706             ];
1707         };
1708 
1709         while (F2 > eps && count < max_it) {
1710             Dt = D(t1, t2);
1711             a = Dt[0][0];
1712             b = Dt[0][1];
1713             c = Dt[1][0];
1714             d = Dt[1][1];
1715 
1716             disc = a * d - b * c;
1717             t1 -= gamma * (d * e - b * f) / disc;
1718             t2 -= gamma * (a * f - c * e) / disc;
1719             f1 = c1.Ft(t1);
1720             f2 = c2.Ft(t2);
1721 
1722             e = f1[1] - f2[1];
1723             f = f1[2] - f2[2];
1724             F2 = e * e + f * f;
1725             count += 1;
1726         }
1727 
1728         return [t1, t2, F2];
1729     },
1730 
1731     /**
1732      * Apply the damped Newton-Raphson algorithm to determine to find a root of a
1733      * function F: R^n to R^n.
1734      *
1735      * @param {Function} F Function with n parameters, returns a vactor of length n.
1736      * @param {Function} D Function returning the Jacobian matrix (n \times n) of F
1737      * @param {Number} n
1738      * @param {Array} t_ini Array of length n, containing start values
1739      * @param {Number} gamma Damping factor should be between 0 and 1. If equal to 1,
1740      * the algorithm is Newton-Raphson.
1741      * @param {Number} eps The algorithm stops if the square norm of the root is less than this eps
1742      * or if the maximum number of steps is reached.
1743      * @param {Number} [max_steps=40] maximum number of steps
1744      * @returns {Array} [t, F2] array of length, containing t, the approximation of the root (array of length n),
1745      * and the square norm of F(t).
1746      */
1747     generalizedDampedNewton: function (F, D, n, t_ini, gamma, eps, max_steps) {
1748         var i,
1749             t = [],
1750             a, b, c, d, e, f,
1751             disc,
1752             Ft, Dt,
1753             F2, vec,
1754             count = 0;
1755 
1756         max_steps = max_steps || 40;
1757 
1758         t = t_ini.slice(0, n);
1759         Ft = F(t, n);
1760 
1761         if (n === 2) {
1762             // Special case n = 2
1763             Ft = F(t, n);
1764             e = Ft[0];
1765             f = Ft[1];
1766             F2 = e * e + f * f;
1767 
1768             while (F2 > eps && count < max_steps) {
1769                 Dt = D(t, n);
1770 
1771                 a = Dt[0][0];
1772                 b = Dt[0][1];
1773                 c = Dt[1][0];
1774                 d = Dt[1][1];
1775 
1776                 disc = a * d - b * c;
1777                 t[0] -= gamma * (d * e - b * f) / disc;
1778                 t[1] -= gamma * (a * f - c * e) / disc;
1779 
1780                 Ft = F(t, n);
1781                 e = Ft[0];
1782                 f = Ft[1];
1783                 F2 = e * e + f * f;
1784 
1785                 count += 1;
1786             }
1787 
1788             return [t, F2];
1789         } else {
1790             // General case, arbitrary n
1791             Ft = F(t, n);
1792             F2 = Mat.innerProduct(Ft, Ft, n);
1793 
1794             while (F2 > eps && count < max_steps) {
1795                 Dt = Mat.inverse(D(t, n));
1796 
1797                 vec = Mat.matVecMult(Dt, Ft);
1798                 for (i = 0; i < n; i++) {
1799                     t[i] -= gamma * vec[i];
1800                 }
1801 
1802                 Ft = F(t, n);
1803                 F2 = Mat.innerProduct(Ft, Ft, n);
1804 
1805                 count += 1;
1806             }
1807 
1808             return [t, F2];
1809         }
1810     },
1811 
1812     /**
1813      * Returns the Lagrange polynomials for curves with equidistant nodes, see
1814      * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
1815      * SIAM Review, Vol 46, No 3, (2004) 501-517.
1816      * The graph of the parametric curve [x(t),y(t)] runs through the given points.
1817      * @param {Array} p Array of JXG.Points
1818      * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve
1819      * f(t) = (x(t), y(t)), a number x1 (which equals 0) and a function x2 defining the curve's domain.
1820      * That means the curve is defined between x1 and x2(). x2 returns the (length of array p minus one).
1821      * @memberof JXG.Math.Numerics
1822      *
1823      * @example
1824      * var p = [];
1825      *
1826      * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'});
1827      * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''});
1828      * p[2] = board.create('point', [1, 4], {size:2, name: ''});
1829      * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'});
1830      *
1831      * // Curve
1832      * var fg = JXG.Math.Numerics.Neville(p);
1833      * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5});
1834      *
1835      * </pre><div id="JXG88a8b3a8-6561-44f5-a678-76bca13fd484" class="jxgbox" style="width: 300px; height: 300px;"></div>
1836      * <script type="text/javascript">
1837      *     (function() {
1838      *         var board = JXG.JSXGraph.initBoard('JXG88a8b3a8-6561-44f5-a678-76bca13fd484',
1839      *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
1840      *     var p = [];
1841      *
1842      *     p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'});
1843      *     p[1] = board.create('point', [-1.5, 5], {size:2, name: ''});
1844      *     p[2] = board.create('point', [1, 4], {size:2, name: ''});
1845      *     p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'});
1846      *
1847      *     // Curve
1848      *     var fg = JXG.Math.Numerics.Neville(p);
1849      *     var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5});
1850      *
1851      *     })();
1852      *
1853      * </script><pre>
1854      *
1855      */
1856     Neville: function (p) {
1857         var w = [],
1858             /** @ignore */
1859             makeFct = function (fun) {
1860                 return function (t, suspendedUpdate) {
1861                     var i,
1862                         d,
1863                         s,
1864                         bin = Mat.binomial,
1865                         len = p.length,
1866                         len1 = len - 1,
1867                         num = 0.0,
1868                         denom = 0.0;
1869 
1870                     if (!suspendedUpdate) {
1871                         s = 1;
1872                         for (i = 0; i < len; i++) {
1873                             w[i] = bin(len1, i) * s;
1874                             s *= -1;
1875                         }
1876                     }
1877 
1878                     d = t;
1879 
1880                     for (i = 0; i < len; i++) {
1881                         if (d === 0) {
1882                             return p[i][fun]();
1883                         }
1884                         s = w[i] / d;
1885                         d -= 1;
1886                         num += p[i][fun]() * s;
1887                         denom += s;
1888                     }
1889                     return num / denom;
1890                 };
1891             },
1892             xfct = makeFct('X'),
1893             yfct = makeFct('Y');
1894 
1895         return [
1896             xfct,
1897             yfct,
1898             0,
1899             function () {
1900                 return p.length - 1;
1901             }
1902         ];
1903     },
1904 
1905     /**
1906      * Calculates second derivatives at the given knots.
1907      * @param {Array} x x values of knots
1908      * @param {Array} y y values of knots
1909      * @returns {Array} Second derivatives of the interpolated function at the knots.
1910      * @see JXG.Math.Numerics.splineEval
1911      * @memberof JXG.Math.Numerics
1912      */
1913     splineDef: function (x, y) {
1914         var pair,
1915             i,
1916             l,
1917             n = Math.min(x.length, y.length),
1918             diag = [],
1919             z = [],
1920             data = [],
1921             dx = [],
1922             delta = [],
1923             F = [];
1924 
1925         if (n === 2) {
1926             return [0, 0];
1927         }
1928 
1929         for (i = 0; i < n; i++) {
1930             pair = { X: x[i], Y: y[i] };
1931             data.push(pair);
1932         }
1933         data.sort(function (a, b) {
1934             return a.X - b.X;
1935         });
1936         for (i = 0; i < n; i++) {
1937             x[i] = data[i].X;
1938             y[i] = data[i].Y;
1939         }
1940 
1941         for (i = 0; i < n - 1; i++) {
1942             dx.push(x[i + 1] - x[i]);
1943         }
1944         for (i = 0; i < n - 2; i++) {
1945             delta.push(
1946                 (6 * (y[i + 2] - y[i + 1])) / dx[i + 1] - (6 * (y[i + 1] - y[i])) / dx[i]
1947             );
1948         }
1949 
1950         // ForwardSolve
1951         diag.push(2 * (dx[0] + dx[1]));
1952         z.push(delta[0]);
1953 
1954         for (i = 0; i < n - 3; i++) {
1955             l = dx[i + 1] / diag[i];
1956             diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]);
1957             z.push(delta[i + 1] - l * z[i]);
1958         }
1959 
1960         // BackwardSolve
1961         F[n - 3] = z[n - 3] / diag[n - 3];
1962         for (i = n - 4; i >= 0; i--) {
1963             F[i] = (z[i] - dx[i + 1] * F[i + 1]) / diag[i];
1964         }
1965 
1966         // Generate f''-Vector
1967         for (i = n - 3; i >= 0; i--) {
1968             F[i + 1] = F[i];
1969         }
1970 
1971         // natural cubic spline
1972         F[0] = 0;
1973         F[n - 1] = 0;
1974 
1975         return F;
1976     },
1977 
1978     /**
1979      * Evaluate points on spline.
1980      * @param {Number|Array} x0 A single float value or an array of values to evaluate
1981      * @param {Array} x x values of knots
1982      * @param {Array} y y values of knots
1983      * @param {Array} F Second derivatives at knots, calculated by {@link JXG.Math.Numerics.splineDef}
1984      * @see JXG.Math.Numerics.splineDef
1985      * @returns {Number|Array} A single value or an array, depending on what is given as x0.
1986      * @memberof JXG.Math.Numerics
1987      */
1988     splineEval: function (x0, x, y, F) {
1989         var i,
1990             j,
1991             a,
1992             b,
1993             c,
1994             d,
1995             x_,
1996             n = Math.min(x.length, y.length),
1997             l = 1,
1998             asArray = false,
1999             y0 = [];
2000 
2001         // number of points to be evaluated
2002         if (Type.isArray(x0)) {
2003             l = x0.length;
2004             asArray = true;
2005         } else {
2006             x0 = [x0];
2007         }
2008 
2009         for (i = 0; i < l; i++) {
2010             // is x0 in defining interval?
2011             if (x0[i] < x[0] || x[i] > x[n - 1]) {
2012                 return NaN;
2013             }
2014 
2015             // determine part of spline in which x0 lies
2016             for (j = 1; j < n; j++) {
2017                 if (x0[i] <= x[j]) {
2018                     break;
2019                 }
2020             }
2021 
2022             j -= 1;
2023 
2024             // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1];
2025             // determine the coefficients of the polynomial in this interval
2026             a = y[j];
2027             b =
2028                 (y[j + 1] - y[j]) / (x[j + 1] - x[j]) -
2029                 ((x[j + 1] - x[j]) / 6) * (F[j + 1] + 2 * F[j]);
2030             c = F[j] / 2;
2031             d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j]));
2032             // evaluate x0[i]
2033             x_ = x0[i] - x[j];
2034             //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_);
2035             y0.push(a + (b + (c + d * x_) * x_) * x_);
2036         }
2037 
2038         if (asArray) {
2039             return y0;
2040         }
2041 
2042         return y0[0];
2043     },
2044 
2045     /**
2046      * Generate a string containing the function term of a polynomial.
2047      * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i.
2048      * @param {Number} deg Degree of the polynomial
2049      * @param {String} varname Name of the variable (usually 'x')
2050      * @param {Number} prec Precision
2051      * @returns {String} A string containing the function term of the polynomial.
2052      * @memberof JXG.Math.Numerics
2053      */
2054     generatePolynomialTerm: function (coeffs, deg, varname, prec) {
2055         var i,
2056             t = [];
2057 
2058         for (i = deg; i >= 0; i--) {
2059             Type.concat(t, ["(", coeffs[i].toPrecision(prec), ")"]);
2060 
2061             if (i > 1) {
2062                 Type.concat(t, ["*", varname, "<sup>", i, "<", "/sup> + "]);
2063             } else if (i === 1) {
2064                 Type.concat(t, ["*", varname, " + "]);
2065             }
2066         }
2067 
2068         return t.join("");
2069     },
2070 
2071     /**
2072      * Computes the polynomial through a given set of coordinates in Lagrange form.
2073      * Returns the Lagrange polynomials, see
2074      * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
2075      * SIAM Review, Vol 46, No 3, (2004) 501-517.
2076      * <p>
2077      * It possesses the method getTerm() which returns the string containing the function term of the polynomial and
2078      * the method getCoefficients() which returns an array containing the coefficients of the polynomial.
2079      * @param {Array} p Array of JXG.Points
2080      * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points.
2081      * @memberof JXG.Math.Numerics
2082      *
2083      * @example
2084      * var p = [];
2085      * p[0] = board.create('point', [-1,2], {size:4});
2086      * p[1] = board.create('point', [0,3], {size:4});
2087      * p[2] = board.create('point', [1,1], {size:4});
2088      * p[3] = board.create('point', [3,-1], {size:4});
2089      * var f = JXG.Math.Numerics.lagrangePolynomial(p);
2090      * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2091      *
2092      * </pre><div id="JXGc058aa6b-74d4-41e1-af94-df06169a2d89" class="jxgbox" style="width: 300px; height: 300px;"></div>
2093      * <script type="text/javascript">
2094      *     (function() {
2095      *         var board = JXG.JSXGraph.initBoard('JXGc058aa6b-74d4-41e1-af94-df06169a2d89',
2096      *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
2097      *     var p = [];
2098      *     p[0] = board.create('point', [-1,2], {size:4});
2099      *     p[1] = board.create('point', [0,3], {size:4});
2100      *     p[2] = board.create('point', [1,1], {size:4});
2101      *     p[3] = board.create('point', [3,-1], {size:4});
2102      *     var f = JXG.Math.Numerics.lagrangePolynomial(p);
2103      *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2104      *
2105      *     })();
2106      *
2107      * </script><pre>
2108      *
2109      * @example
2110      * var points = [];
2111      * points[0] = board.create('point', [-1,2], {size:4});
2112      * points[1] = board.create('point', [0, 0], {size:4});
2113      * points[2] = board.create('point', [2, 1], {size:4});
2114      *
2115      * var f = JXG.Math.Numerics.lagrangePolynomial(points);
2116      * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2117      * var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
2118      * var txt2 = board.create('text', [-3, -6,  () => f.getCoefficients()], {fontSize: 12});
2119      *
2120      * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eecbb" class="jxgbox" style="width: 300px; height: 300px;"></div>
2121      * <script type="text/javascript">
2122      *     (function() {
2123      *         var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eecbb',
2124      *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
2125      *     var points = [];
2126      *     points[0] = board.create('point', [-1,2], {size:4});
2127      *     points[1] = board.create('point', [0, 0], {size:4});
2128      *     points[2] = board.create('point', [2, 1], {size:4});
2129      *
2130      *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
2131      *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2132      *     var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
2133      *     var txt2 = board.create('text', [-3, -6,  () => f.getCoefficients()], {fontSize: 12});
2134      *
2135      *     })();
2136      *
2137      * </script><pre>
2138      *
2139      */
2140     lagrangePolynomial: function (p) {
2141         var w = [],
2142             that = this,
2143             /** @ignore */
2144             fct = function (x, suspendedUpdate) {
2145                 var i, // j,
2146                     k,
2147                     xi,
2148                     s, //M,
2149                     len = p.length,
2150                     num = 0,
2151                     denom = 0;
2152 
2153                 if (!suspendedUpdate) {
2154                     for (i = 0; i < len; i++) {
2155                         w[i] = 1.0;
2156                         xi = p[i].X();
2157 
2158                         for (k = 0; k < len; k++) {
2159                             if (k !== i) {
2160                                 w[i] *= xi - p[k].X();
2161                             }
2162                         }
2163 
2164                         w[i] = 1 / w[i];
2165                     }
2166 
2167                     // M = [];
2168                     // for (k = 0; k < len; k++) {
2169                     //     M.push([1]);
2170                     // }
2171                 }
2172 
2173                 for (i = 0; i < len; i++) {
2174                     xi = p[i].X();
2175 
2176                     if (x === xi) {
2177                         return p[i].Y();
2178                     }
2179 
2180                     s = w[i] / (x - xi);
2181                     denom += s;
2182                     num += s * p[i].Y();
2183                 }
2184 
2185                 return num / denom;
2186             };
2187 
2188         /**
2189          * Get the term of the Lagrange polynomial as string.
2190          * Calls {@link JXG.Math.Numerics#lagrangePolynomialTerm}.
2191          *
2192          * @name JXG.Math.Numerics.lagrangePolynomial#getTerm
2193          * @param {Number} digits Number of digits of the coefficients
2194          * @param {String} param Variable name
2195          * @param {String} dot Dot symbol
2196          * @returns {String} containing the term of Lagrange polynomial as string.
2197          * @see JXG.Math.Numerics.lagrangePolynomialTerm
2198          * @example
2199          * var points = [];
2200          * points[0] = board.create('point', [-1,2], {size:4});
2201          * points[1] = board.create('point', [0, 0], {size:4});
2202          * points[2] = board.create('point', [2, 1], {size:4});
2203          *
2204          * var f = JXG.Math.Numerics.lagrangePolynomial(points);
2205          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2206          * var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
2207          *
2208          * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eeccf" class="jxgbox" style="width: 300px; height: 300px;"></div>
2209          * <script type="text/javascript">
2210          *     (function() {
2211          *         var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eeccf',
2212          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
2213          *     var points = [];
2214          *     points[0] = board.create('point', [-1,2], {size:4});
2215          *     points[1] = board.create('point', [0, 0], {size:4});
2216          *     points[2] = board.create('point', [2, 1], {size:4});
2217          *
2218          *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
2219          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2220          *     var txt = board.create('text', [-3, -4,  () => f.getTerm(2, 't', ' * ')], {fontSize: 16});
2221          *
2222          *     })();
2223          *
2224          * </script><pre>
2225          *
2226          */
2227         fct.getTerm = function (digits, param, dot) {
2228             return that.lagrangePolynomialTerm(p, digits, param, dot)();
2229         };
2230 
2231         /**
2232          * Get the coefficients of the Lagrange polynomial as array. The leading
2233          * coefficient is at position 0.
2234          * Calls {@link JXG.Math.Numerics#lagrangePolynomialCoefficients}.
2235          *
2236          * @name JXG.Math.Numerics.lagrangePolynomial#getCoefficients
2237          * @returns {Array} containing the coefficients of the Lagrange polynomial.
2238          * @see JXG.Math.Numerics.lagrangePolynomial.getTerm
2239          * @see JXG.Math.Numerics.lagrangePolynomialTerm
2240          * @see JXG.Math.Numerics.lagrangePolynomialCoefficients
2241          * @example
2242          * var points = [];
2243          * points[0] = board.create('point', [-1,2], {size:4});
2244          * points[1] = board.create('point', [0, 0], {size:4});
2245          * points[2] = board.create('point', [2, 1], {size:4});
2246          *
2247          * var f = JXG.Math.Numerics.lagrangePolynomial(points);
2248          * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2249          * var txt = board.create('text', [1, -4,  () => f.getCoefficients()], {fontSize: 10});
2250          *
2251          * </pre><div id="JXG52a883a5-2e0c-4caf-8f84-8650c173c365" class="jxgbox" style="width: 300px; height: 300px;"></div>
2252          * <script type="text/javascript">
2253          *     (function() {
2254          *         var board = JXG.JSXGraph.initBoard('JXG52a883a5-2e0c-4caf-8f84-8650c173c365',
2255          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
2256          *     var points = [];
2257          *     points[0] = board.create('point', [-1,2], {size:4});
2258          *     points[1] = board.create('point', [0, 0], {size:4});
2259          *     points[2] = board.create('point', [2, 1], {size:4});
2260          *
2261          *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
2262          *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2263          *     var txt = board.create('text', [1, -4,  () => f.getCoefficients()], {fontSize: 10});
2264          *
2265          *     })();
2266          *
2267          * </script><pre>
2268          *
2269          */
2270         fct.getCoefficients = function () {
2271             return that.lagrangePolynomialCoefficients(p)();
2272         };
2273 
2274         return fct;
2275     },
2276 
2277     /**
2278      * Determine the Lagrange polynomial through an array of points and
2279      * return the term of the polynomial as string.
2280      *
2281      * @param {Array} points Array of JXG.Points
2282      * @param {Number} digits Number of decimal digits of the coefficients
2283      * @param {String} param Name of the parameter. Default: 'x'.
2284      * @param {String} dot Multiplication symbol. Default: ' * '.
2285      * @returns {Function} returning the Lagrange polynomial term through
2286      *    the supplied points as string
2287      * @memberof JXG.Math.Numerics
2288      *
2289      * @example
2290      * var points = [];
2291      * points[0] = board.create('point', [-1,2], {size:4});
2292      * points[1] = board.create('point', [0, 0], {size:4});
2293      * points[2] = board.create('point', [2, 1], {size:4});
2294      *
2295      * var f = JXG.Math.Numerics.lagrangePolynomial(points);
2296      * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2297      *
2298      * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * ');
2299      * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16});
2300      *
2301      * </pre><div id="JXGd45e9e96-7526-486d-aa43-e1178d5f2baa" class="jxgbox" style="width: 300px; height: 300px;"></div>
2302      * <script type="text/javascript">
2303      *     (function() {
2304      *         var board = JXG.JSXGraph.initBoard('JXGd45e9e96-7526-486d-aa43-e1178d5f2baa',
2305      *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
2306      *     var points = [];
2307      *     points[0] = board.create('point', [-1,2], {size:4});
2308      *     points[1] = board.create('point', [0, 0], {size:4});
2309      *     points[2] = board.create('point', [2, 1], {size:4});
2310      *
2311      *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
2312      *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2313      *
2314      *     var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * ');
2315      *     var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16});
2316      *
2317      *     })();
2318      *
2319      * </script><pre>
2320      *
2321      */
2322     lagrangePolynomialTerm: function (points, digits, param, dot) {
2323         var that = this;
2324 
2325         return function () {
2326             var len = points.length,
2327                 coeffs = [],
2328                 isLeading = true,
2329                 n, t, j, c;
2330 
2331             param = param || 'x';
2332             if (dot === undefined) {
2333                 dot = " * ";
2334             }
2335 
2336             n = len - 1; // (Max) degree of the polynomial
2337             coeffs = that.lagrangePolynomialCoefficients(points)();
2338 
2339             t = "";
2340             for (j = 0; j < coeffs.length; j++) {
2341                 c = coeffs[j];
2342                 if (Math.abs(c) < Mat.eps) {
2343                     continue;
2344                 }
2345                 if (JXG.exists(digits)) {
2346                     c = Env._round10(c, -digits);
2347                 }
2348                 if (isLeading) {
2349                     t += c > 0 ? c : "-" + -c;
2350                     isLeading = false;
2351                 } else {
2352                     t += c > 0 ? " + " + c : " - " + -c;
2353                 }
2354 
2355                 if (n - j > 1) {
2356                     t += dot + param + "^" + (n - j);
2357                 } else if (n - j === 1) {
2358                     t += dot + param;
2359                 }
2360             }
2361             return t; // board.jc.manipulate('f = map(x) -> ' + t + ';');
2362         };
2363     },
2364 
2365     /**
2366      * Determine the Lagrange polynomial through an array of points and
2367      * return the coefficients of the polynomial as array.
2368      * The leading coefficient is at position 0.
2369      *
2370      * @param {Array} points Array of JXG.Points
2371      * @returns {Function} returning the coefficients of the Lagrange polynomial through
2372      *    the supplied points.
2373      * @memberof JXG.Math.Numerics
2374      *
2375      * @example
2376      * var points = [];
2377      * points[0] = board.create('point', [-1,2], {size:4});
2378      * points[1] = board.create('point', [0, 0], {size:4});
2379      * points[2] = board.create('point', [2, 1], {size:4});
2380      *
2381      * var f = JXG.Math.Numerics.lagrangePolynomial(points);
2382      * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2383      *
2384      * var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points);
2385      * var txt = board.create('text', [1, -4, f_arr], {fontSize: 10});
2386      *
2387      * </pre><div id="JXG1778f0d1-a420-473f-99e8-1755ef4be97e" class="jxgbox" style="width: 300px; height: 300px;"></div>
2388      * <script type="text/javascript">
2389      *     (function() {
2390      *         var board = JXG.JSXGraph.initBoard('JXG1778f0d1-a420-473f-99e8-1755ef4be97e',
2391      *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
2392      *     var points = [];
2393      *     points[0] = board.create('point', [-1,2], {size:4});
2394      *     points[1] = board.create('point', [0, 0], {size:4});
2395      *     points[2] = board.create('point', [2, 1], {size:4});
2396      *
2397      *     var f = JXG.Math.Numerics.lagrangePolynomial(points);
2398      *     var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3});
2399      *
2400      *     var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points);
2401      *     var txt = board.create('text', [1, -4, f_arr], {fontSize: 10});
2402      *
2403      *     })();
2404      *
2405      * </script><pre>
2406      *
2407      */
2408     lagrangePolynomialCoefficients: function (points) {
2409         return function () {
2410             var len = points.length,
2411                 zeroes = [],
2412                 coeffs = [],
2413                 coeffs_sum = [],
2414                 i, j, c, p;
2415 
2416             // n = len - 1; // (Max) degree of the polynomial
2417             for (j = 0; j < len; j++) {
2418                 coeffs_sum[j] = 0;
2419             }
2420 
2421             for (i = 0; i < len; i++) {
2422                 c = points[i].Y();
2423                 p = points[i].X();
2424                 zeroes = [];
2425                 for (j = 0; j < len; j++) {
2426                     if (j !== i) {
2427                         c /= p - points[j].X();
2428                         zeroes.push(points[j].X());
2429                     }
2430                 }
2431                 coeffs = [1].concat(Mat.Vieta(zeroes));
2432                 for (j = 0; j < coeffs.length; j++) {
2433                     coeffs_sum[j] += (j % 2 === 1 ? -1 : 1) * coeffs[j] * c;
2434                 }
2435             }
2436 
2437             return coeffs_sum;
2438         };
2439     },
2440 
2441     /**
2442      * Determine the coefficients of a cardinal spline polynom, See
2443      * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections
2444      * @param  {Number} x1 point 1
2445      * @param  {Number} x2 point 2
2446      * @param  {Number} t1 tangent slope 1
2447      * @param  {Number} t2 tangent slope 2
2448      * @return {Array}    coefficents array c for the polynomial t maps to
2449      * c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t
2450      */
2451     _initCubicPoly: function (x1, x2, t1, t2) {
2452         return [x1, t1, -3 * x1 + 3 * x2 - 2 * t1 - t2, 2 * x1 - 2 * x2 + t1 + t2];
2453     },
2454 
2455     /**
2456      * Computes the cubic cardinal spline curve through a given set of points. The curve
2457      * is uniformly parametrized.
2458      * Two artificial control points at the beginning and the end are added.
2459      *
2460      * The implementation (especially the centripetal parametrization) is from
2461      * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections .
2462      * @param {Array} points Array consisting of JXG.Points.
2463      * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1.
2464      * tau=1/2 give Catmull-Rom splines.
2465      * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and
2466      * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal".
2467      * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
2468      * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value,
2469      * and a function simply returning the length of the points array
2470      * minus three.
2471      * @memberof JXG.Math.Numerics
2472      */
2473     CardinalSpline: function (points, tau_param, type) {
2474         var p,
2475             coeffs = [],
2476             makeFct,
2477             tau, _tau,
2478             that = this;
2479 
2480         if (Type.isFunction(tau_param)) {
2481             _tau = tau_param;
2482         } else {
2483             _tau = function () {
2484                 return tau_param;
2485             };
2486         }
2487 
2488         if (type === undefined) {
2489             type = 'uniform';
2490         }
2491 
2492         /** @ignore */
2493         makeFct = function (which) {
2494             return function (t, suspendedUpdate) {
2495                 var s,
2496                     c,
2497                     // control point at the beginning and at the end
2498                     first,
2499                     last,
2500                     t1,
2501                     t2,
2502                     dt0,
2503                     dt1,
2504                     dt2,
2505                     // dx, dy,
2506                     len;
2507 
2508                 if (points.length < 2) {
2509                     return NaN;
2510                 }
2511 
2512                 if (!suspendedUpdate) {
2513                     tau = _tau();
2514 
2515                     // New point list p: [first, points ..., last]
2516                     first = {
2517                         X: function () {
2518                             return 2 * points[0].X() - points[1].X();
2519                         },
2520                         Y: function () {
2521                             return 2 * points[0].Y() - points[1].Y();
2522                         },
2523                         Dist: function (p) {
2524                             var dx = this.X() - p.X(),
2525                                 dy = this.Y() - p.Y();
2526                             return Mat.hypot(dx, dy);
2527                         }
2528                     };
2529 
2530                     last = {
2531                         X: function () {
2532                             return (
2533                                 2 * points[points.length - 1].X() -
2534                                 points[points.length - 2].X()
2535                             );
2536                         },
2537                         Y: function () {
2538                             return (
2539                                 2 * points[points.length - 1].Y() -
2540                                 points[points.length - 2].Y()
2541                             );
2542                         },
2543                         Dist: function (p) {
2544                             var dx = this.X() - p.X(),
2545                                 dy = this.Y() - p.Y();
2546                             return Mat.hypot(dx, dy);
2547                         }
2548                     };
2549 
2550                     p = [first].concat(points, [last]);
2551                     len = p.length;
2552 
2553                     coeffs[which] = [];
2554 
2555                     for (s = 0; s < len - 3; s++) {
2556                         if (type === 'centripetal') {
2557                             // The order is important, since p[0].coords === undefined
2558                             dt0 = p[s].Dist(p[s + 1]);
2559                             dt1 = p[s + 2].Dist(p[s + 1]);
2560                             dt2 = p[s + 3].Dist(p[s + 2]);
2561 
2562                             dt0 = Math.sqrt(dt0);
2563                             dt1 = Math.sqrt(dt1);
2564                             dt2 = Math.sqrt(dt2);
2565 
2566                             if (dt1 < Mat.eps) {
2567                                 dt1 = 1.0;
2568                             }
2569                             if (dt0 < Mat.eps) {
2570                                 dt0 = dt1;
2571                             }
2572                             if (dt2 < Mat.eps) {
2573                                 dt2 = dt1;
2574                             }
2575 
2576                             t1 =
2577                                 (p[s + 1][which]() - p[s][which]()) / dt0 -
2578                                 (p[s + 2][which]() - p[s][which]()) / (dt1 + dt0) +
2579                                 (p[s + 2][which]() - p[s + 1][which]()) / dt1;
2580 
2581                             t2 =
2582                                 (p[s + 2][which]() - p[s + 1][which]()) / dt1 -
2583                                 (p[s + 3][which]() - p[s + 1][which]()) / (dt2 + dt1) +
2584                                 (p[s + 3][which]() - p[s + 2][which]()) / dt2;
2585 
2586                             t1 *= dt1;
2587                             t2 *= dt1;
2588 
2589                             coeffs[which][s] = that._initCubicPoly(
2590                                 p[s + 1][which](),
2591                                 p[s + 2][which](),
2592                                 tau * t1,
2593                                 tau * t2
2594                             );
2595                         } else {
2596                             coeffs[which][s] = that._initCubicPoly(
2597                                 p[s + 1][which](),
2598                                 p[s + 2][which](),
2599                                 tau * (p[s + 2][which]() - p[s][which]()),
2600                                 tau * (p[s + 3][which]() - p[s + 1][which]())
2601                             );
2602                         }
2603                     }
2604                 }
2605 
2606                 if (isNaN(t)) {
2607                     return NaN;
2608                 }
2609 
2610                 len = points.length;
2611                 // This is necessary for our advanced plotting algorithm:
2612                 if (t <= 0.0) {
2613                     return points[0][which]();
2614                 }
2615                 if (t >= len) {
2616                     return points[len - 1][which]();
2617                 }
2618 
2619                 s = Math.floor(t);
2620                 if (s === t) {
2621                     return points[s][which]();
2622                 }
2623 
2624                 t -= s;
2625                 c = coeffs[which][s];
2626                 if (c === undefined) {
2627                     return NaN;
2628                 }
2629 
2630                 return ((c[3] * t + c[2]) * t + c[1]) * t + c[0];
2631             };
2632         };
2633 
2634         return [
2635             makeFct('X'),
2636             makeFct('Y'),
2637             0,
2638             function () {
2639                 return points.length - 1;
2640             }
2641         ];
2642     },
2643 
2644     /**
2645      * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve
2646      * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5.
2647      * Two artificial control points at the beginning and the end are added.
2648      * @param {Array} points Array consisting of JXG.Points.
2649      * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and
2650      * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal".
2651      * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
2652      * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
2653      * returning the length of the points array minus three.
2654      * @memberof JXG.Math.Numerics
2655      */
2656     CatmullRomSpline: function (points, type) {
2657         return this.CardinalSpline(points, 0.5, type);
2658     },
2659 
2660     /**
2661      * Computes the regression polynomial of a given degree through a given set of coordinates.
2662      * Returns the regression polynomial function.
2663      * @param {Number|function|Slider} degree number, function or slider.
2664      * Either
2665      * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in
2666      * an array of {@link JXG.Point}s or {@link JXG.Coords}.
2667      * In the latter case, the <tt>dataY</tt> parameter will be ignored.
2668      * @param {Array} dataY Array containing the y-coordinates of the data set,
2669      * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree.
2670      * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
2671      * The function returned will throw an exception, if the data set is malformed.
2672      * @memberof JXG.Math.Numerics
2673      */
2674     regressionPolynomial: function (degree, dataX, dataY) {
2675         var coeffs, deg, dX, dY, inputType, fct,
2676             term = "";
2677 
2678         // Slider
2679         if (Type.isPoint(degree) && Type.isFunction(degree.Value)) {
2680             /** @ignore */
2681             deg = function () {
2682                 return degree.Value();
2683             };
2684             // function
2685         } else if (Type.isFunction(degree)) {
2686             deg = degree;
2687             // number
2688         } else if (Type.isNumber(degree)) {
2689             /** @ignore */
2690             deg = function () {
2691                 return degree;
2692             };
2693         } else {
2694             throw new Error(
2695                 "JSXGraph: Can't create regressionPolynomial from degree of type'" +
2696                 typeof degree +
2697                 "'."
2698             );
2699         }
2700 
2701         // Parameters degree, dataX, dataY
2702         if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) {
2703             inputType = 0;
2704             // Parameters degree, point array
2705         } else if (
2706             arguments.length === 2 &&
2707             Type.isArray(dataX) &&
2708             dataX.length > 0 &&
2709             Type.isPoint(dataX[0])
2710         ) {
2711             inputType = 1;
2712         } else if (
2713             arguments.length === 2 &&
2714             Type.isArray(dataX) &&
2715             dataX.length > 0 &&
2716             dataX[0].usrCoords &&
2717             dataX[0].scrCoords
2718         ) {
2719             inputType = 2;
2720         } else {
2721             throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters.");
2722         }
2723 
2724         /** @ignore */
2725         fct = function (x, suspendedUpdate) {
2726             var i, j,
2727                 M, MT, y, B, c, s, d,
2728                 // input data
2729                 len = dataX.length;
2730 
2731             d = Math.floor(deg());
2732 
2733             if (!suspendedUpdate) {
2734                 // point list as input
2735                 if (inputType === 1) {
2736                     dX = [];
2737                     dY = [];
2738 
2739                     for (i = 0; i < len; i++) {
2740                         dX[i] = dataX[i].X();
2741                         dY[i] = dataX[i].Y();
2742                     }
2743                 }
2744 
2745                 if (inputType === 2) {
2746                     dX = [];
2747                     dY = [];
2748 
2749                     for (i = 0; i < len; i++) {
2750                         dX[i] = dataX[i].usrCoords[1];
2751                         dY[i] = dataX[i].usrCoords[2];
2752                     }
2753                 }
2754 
2755                 // check for functions
2756                 if (inputType === 0) {
2757                     dX = [];
2758                     dY = [];
2759 
2760                     for (i = 0; i < len; i++) {
2761                         if (Type.isFunction(dataX[i])) {
2762                             dX.push(dataX[i]());
2763                         } else {
2764                             dX.push(dataX[i]);
2765                         }
2766 
2767                         if (Type.isFunction(dataY[i])) {
2768                             dY.push(dataY[i]());
2769                         } else {
2770                             dY.push(dataY[i]);
2771                         }
2772                     }
2773                 }
2774 
2775                 M = [];
2776                 for (j = 0; j < len; j++) {
2777                     M.push([1]);
2778                 }
2779                 for (i = 1; i <= d; i++) {
2780                     for (j = 0; j < len; j++) {
2781                         M[j][i] = M[j][i - 1] * dX[j];
2782                     }
2783                 }
2784 
2785                 y = dY;
2786                 MT = Mat.transpose(M);
2787                 B = Mat.matMatMult(MT, M);
2788                 c = Mat.matVecMult(MT, y);
2789                 coeffs = Mat.Numerics.Gauss(B, c);
2790                 term = Mat.Numerics.generatePolynomialTerm(coeffs, d, "x", 3);
2791             }
2792 
2793             // Horner's scheme to evaluate polynomial
2794             s = coeffs[d];
2795             for (i = d - 1; i >= 0; i--) {
2796                 s = s * x + coeffs[i];
2797             }
2798 
2799             return s;
2800         };
2801 
2802         /** @ignore */
2803         fct.getTerm = function () {
2804             return term;
2805         };
2806 
2807         return fct;
2808     },
2809 
2810     /**
2811      * Computes the cubic Bezier curve through a given set of points.
2812      * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}.
2813      * The points at position k with k mod 3 = 0 are the data points,
2814      * points at position k with k mod 3 = 1 or 2 are the control points.
2815      * @returns {Array} An array consisting of two functions of one parameter t which return the
2816      * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting
2817      * no parameters and returning one third of the length of the points.
2818      * @memberof JXG.Math.Numerics
2819      */
2820     bezier: function (points) {
2821         var len,
2822             flen,
2823             /** @ignore */
2824             makeFct = function (which) {
2825                 return function (t, suspendedUpdate) {
2826                     var z = Math.floor(t) * 3,
2827                         t0 = t % 1,
2828                         t1 = 1 - t0;
2829 
2830                     if (!suspendedUpdate) {
2831                         flen = 3 * Math.floor((points.length - 1) / 3);
2832                         len = Math.floor(flen / 3);
2833                     }
2834 
2835                     if (t < 0) {
2836                         return points[0][which]();
2837                     }
2838 
2839                     if (t >= len) {
2840                         return points[flen][which]();
2841                     }
2842 
2843                     if (isNaN(t)) {
2844                         return NaN;
2845                     }
2846 
2847                     return (
2848                         t1 * t1 * (t1 * points[z][which]() + 3 * t0 * points[z + 1][which]()) +
2849                         (3 * t1 * points[z + 2][which]() + t0 * points[z + 3][which]()) *
2850                         t0 *
2851                         t0
2852                     );
2853                 };
2854             };
2855 
2856         return [
2857             makeFct('X'),
2858             makeFct('Y'),
2859             0,
2860             function () {
2861                 return Math.floor(points.length / 3);
2862             }
2863         ];
2864     },
2865 
2866     /**
2867      * Computes the B-spline curve of order k (order = degree+1) through a given set of points.
2868      * @param {Array} points Array consisting of JXG.Points.
2869      * @param {Number} order Order of the B-spline curve.
2870      * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
2871      * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply
2872      * returning the length of the points array minus one.
2873      * @memberof JXG.Math.Numerics
2874      */
2875     bspline: function (points, order) {
2876         var knots,
2877             _knotVector = function (n, k) {
2878                 var j,
2879                     kn = [];
2880 
2881                 for (j = 0; j < n + k + 1; j++) {
2882                     if (j < k) {
2883                         kn[j] = 0.0;
2884                     } else if (j <= n) {
2885                         kn[j] = j - k + 1;
2886                     } else {
2887                         kn[j] = n - k + 2;
2888                     }
2889                 }
2890 
2891                 return kn;
2892             },
2893             _evalBasisFuncs = function (t, kn, k, s) {
2894                 var i,
2895                     j,
2896                     a,
2897                     b,
2898                     den,
2899                     N = [];
2900 
2901                 if (kn[s] <= t && t < kn[s + 1]) {
2902                     N[s] = 1;
2903                 } else {
2904                     N[s] = 0;
2905                 }
2906 
2907                 for (i = 2; i <= k; i++) {
2908                     for (j = s - i + 1; j <= s; j++) {
2909                         if (j <= s - i + 1 || j < 0) {
2910                             a = 0.0;
2911                         } else {
2912                             a = N[j];
2913                         }
2914 
2915                         if (j >= s) {
2916                             b = 0.0;
2917                         } else {
2918                             b = N[j + 1];
2919                         }
2920 
2921                         den = kn[j + i - 1] - kn[j];
2922 
2923                         if (den === 0) {
2924                             N[j] = 0;
2925                         } else {
2926                             N[j] = ((t - kn[j]) / den) * a;
2927                         }
2928 
2929                         den = kn[j + i] - kn[j + 1];
2930 
2931                         if (den !== 0) {
2932                             N[j] += ((kn[j + i] - t) / den) * b;
2933                         }
2934                     }
2935                 }
2936                 return N;
2937             },
2938             /** @ignore */
2939             makeFct = function (which) {
2940                 return function (t, suspendedUpdate) {
2941                     var y,
2942                         j,
2943                         s,
2944                         N = [],
2945                         len = points.length,
2946                         n = len - 1,
2947                         k = order;
2948 
2949                     if (n <= 0) {
2950                         return NaN;
2951                     }
2952 
2953                     if (n + 2 <= k) {
2954                         k = n + 1;
2955                     }
2956 
2957                     if (t <= 0) {
2958                         return points[0][which]();
2959                     }
2960 
2961                     if (t >= n - k + 2) {
2962                         return points[n][which]();
2963                     }
2964 
2965                     s = Math.floor(t) + k - 1;
2966                     knots = _knotVector(n, k);
2967                     N = _evalBasisFuncs(t, knots, k, s);
2968 
2969                     y = 0.0;
2970                     for (j = s - k + 1; j <= s; j++) {
2971                         if (j < len && j >= 0) {
2972                             y += points[j][which]() * N[j];
2973                         }
2974                     }
2975 
2976                     return y;
2977                 };
2978             };
2979 
2980         return [
2981             makeFct('X'),
2982             makeFct('Y'),
2983             0,
2984             function () {
2985                 return points.length - 1;
2986             }
2987         ];
2988     },
2989 
2990     /**
2991      * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through,
2992      * see {@link JXG.Curve#updateCurve}
2993      * and {@link JXG.Curve#hasPoint}.
2994      * @param {function} f Function in one variable to be differentiated.
2995      * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a
2996      * method of an object and contains a reference to its parent object via "this".
2997      * @returns {function} Derivative function of a given function f.
2998      * @memberof JXG.Math.Numerics
2999      */
3000     D: function (f, obj) {
3001         if (!Type.exists(obj)) {
3002             return function (x, suspendedUpdate) {
3003                 var h = 0.00001,
3004                     h2 = h * 2.0;
3005 
3006                 // Experiments with Richardsons rule
3007                 /*
3008                     var phi = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
3009                     var phi2;
3010                     h *= 0.5;
3011                     h2 *= 0.5;
3012                     phi2 = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
3013 
3014                     return phi2 + (phi2 - phi) / 3.0;
3015                     */
3016                 return (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2;
3017             };
3018         }
3019 
3020         return function (x, suspendedUpdate) {
3021             var h = 0.00001,
3022                 h2 = h * 2.0;
3023 
3024             return (
3025                 (f.apply(obj, [x + h, suspendedUpdate]) -
3026                     f.apply(obj, [x - h, suspendedUpdate])) /
3027                 h2
3028             );
3029         };
3030     },
3031 
3032     /**
3033      * Evaluate the function term for {@link JXG.Math.Numerics.riemann}.
3034      * @private
3035      * @param {Number} x function argument
3036      * @param {function} f JavaScript function returning a number
3037      * @param {String} type Name of the Riemann sum type, e.g. 'lower'.
3038      * @param {Number} delta Width of the bars in user coordinates
3039      * @returns {Number} Upper (delta > 0) or lower (delta < 0) value of the bar containing x of the Riemann sum.
3040      * @see JXG.Math.Numerics.riemann
3041      * @private
3042      * @memberof JXG.Math.Numerics
3043      */
3044     _riemannValue: function (x, f, type, delta) {
3045         var y, y1, x1, delta1;
3046 
3047         if (delta < 0) {
3048             // delta is negative if the lower function term is evaluated
3049             if (type !== 'trapezoidal') {
3050                 x = x + delta;
3051             }
3052             delta *= -1;
3053             if (type === 'lower') {
3054                 type = 'upper';
3055             } else if (type === 'upper') {
3056                 type = 'lower';
3057             }
3058         }
3059 
3060         delta1 = delta * 0.01; // for 'lower' and 'upper'
3061 
3062         if (type === 'right') {
3063             y = f(x + delta);
3064         } else if (type === 'middle') {
3065             y = f(x + delta * 0.5);
3066         } else if (type === "left" || type === 'trapezoidal') {
3067             y = f(x);
3068         } else if (type === 'lower') {
3069             y = f(x);
3070 
3071             for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
3072                 y1 = f(x1);
3073 
3074                 if (y1 < y) {
3075                     y = y1;
3076                 }
3077             }
3078 
3079             y1 = f(x + delta);
3080             if (y1 < y) {
3081                 y = y1;
3082             }
3083         } else if (type === 'upper') {
3084             y = f(x);
3085 
3086             for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
3087                 y1 = f(x1);
3088                 if (y1 > y) {
3089                     y = y1;
3090                 }
3091             }
3092 
3093             y1 = f(x + delta);
3094             if (y1 > y) {
3095                 y = y1;
3096             }
3097         } else if (type === 'random') {
3098             y = f(x + delta * Math.random());
3099         } else if (type === 'simpson') {
3100             y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0;
3101         } else {
3102             y = f(x); // default is lower
3103         }
3104 
3105         return y;
3106     },
3107 
3108     /**
3109      * Helper function to create curve which displays Riemann sums.
3110      * Compute coordinates for the rectangles showing the Riemann sum.
3111      * <p>
3112      * In case of type "simpson" and "trapezoidal", the horizontal line approximating the function value
3113      * is replaced by a parabola or a secant. IN case of "simpson",
3114      * the parabola is approximated visually by a polygonal chain of fixed step width.
3115      *
3116      * @param {Function|Array} f Function or array of two functions.
3117      * If f is a function the integral of this function is approximated by the Riemann sum.
3118      * If f is an array consisting of two functions the area between the two functions is filled
3119      * by the Riemann sum bars.
3120      * @param {Number} n number of rectangles.
3121      * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'.
3122      * "simpson" is Simpson's 1/3 rule.
3123      * @param {Number} start Left border of the approximation interval
3124      * @param {Number} end Right border of the approximation interval
3125      * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This
3126      * 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
3127      * rectangles.
3128      * @memberof JXG.Math.Numerics
3129      */
3130     riemann: function (gf, n, type, start, end) {
3131         var i, delta,
3132             k, a, b, c, f0, f1, f2, xx, h,
3133             steps = 30, // Constant step width for Simpson's rule
3134             xarr = [],
3135             yarr = [],
3136             x = start,
3137             sum = 0,
3138             y, f, g;
3139 
3140         if (Type.isArray(gf)) {
3141             g = gf[0];
3142             f = gf[1];
3143         } else {
3144             f = gf;
3145         }
3146 
3147         n = Math.floor(n);
3148 
3149         if (n <= 0) {
3150             return [xarr, yarr, sum];
3151         }
3152 
3153         delta = (end - start) / n;
3154 
3155         // Horizontal lines defined by function
3156         // Go forward
3157         for (i = 0; i < n; i++) {
3158             if (type === 'simpson') {
3159                 sum += this._riemannValue(x, f, type, delta) * delta;
3160 
3161                 h = delta * 0.5;
3162                 f0 = f(x);
3163                 f1 = f(x + h);
3164                 f2 = f(x + 2 * h);
3165 
3166                 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5;
3167                 b = (f2 - f0) / (2 * h);
3168                 c = f1;
3169                 for (k = 0; k < steps; k++) {
3170                     xx = k * delta / steps - h;
3171                     xarr.push(x + xx + h);
3172                     yarr.push(a * xx * xx + b * xx + c);
3173                 }
3174                 x += delta;
3175                 y = f2;
3176             } else {
3177                 y = this._riemannValue(x, f, type, delta);
3178                 xarr.push(x);
3179                 yarr.push(y);
3180 
3181                 x += delta;
3182                 if (type === 'trapezoidal') {
3183                     f2 = f(x);
3184                     sum += (y + f2) * 0.5 * delta;
3185                     y = f2;
3186                 } else {
3187                     sum += y * delta;
3188                 }
3189             }
3190             xarr.push(x);
3191             yarr.push(y);
3192         }
3193 
3194         // Horizontal lines defined by x-axis or second function
3195         // Go backwards
3196         for (i = 0; i < n; i++) {
3197             if (type === "simpson" && g) {
3198                 sum -= this._riemannValue(x, g, type, -delta) * delta;
3199 
3200                 h = delta * 0.5;
3201                 f0 = g(x);
3202                 f1 = g(x - h);
3203                 f2 = g(x - 2 * h);
3204 
3205                 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5;
3206                 b = (f2 - f0) / (2 * h);
3207                 c = f1;
3208                 for (k = 0; k < steps; k++) {
3209                     xx = k * delta / steps - h;
3210                     xarr.push(x - xx - h);
3211                     yarr.push(a * xx * xx + b * xx + c);
3212                 }
3213                 x -= delta;
3214                 y = f2;
3215             } else {
3216                 if (g) {
3217                     y = this._riemannValue(x, g, type, -delta);
3218                 } else {
3219                     y = 0.0;
3220                 }
3221                 xarr.push(x);
3222                 yarr.push(y);
3223 
3224                 x -= delta;
3225                 if (g) {
3226                     if (type === 'trapezoidal') {
3227                         f2 = g(x);
3228                         sum -= (y + f2) * 0.5 * delta;
3229                         y = f2;
3230                     } else {
3231                         sum -= y * delta;
3232                     }
3233                 }
3234             }
3235             xarr.push(x);
3236             yarr.push(y);
3237 
3238             // Close rectangle
3239             xarr.push(x);
3240             if (type === 'simpson') {
3241                 yarr.push(f(x));
3242             } else {
3243                 y = this._riemannValue(x, f, type, delta);
3244                 yarr.push(y);
3245             }
3246         }
3247 
3248         return [xarr, yarr, sum];
3249     },
3250 
3251     /**
3252      * Approximate the integral by Riemann sums.
3253      * Compute the area described by the riemann sum rectangles.
3254      *
3255      * If there is an element of type {@link Riemannsum}, then it is more efficient
3256      * to use the method JXG.Curve.Value() of this element instead.
3257      *
3258      * @param {Function_Array} f Function or array of two functions.
3259      * If f is a function the integral of this function is approximated by the Riemann sum.
3260      * If f is an array consisting of two functions the area between the two functions is approximated
3261      * by the Riemann sum.
3262      * @param {Number} n number of rectangles.
3263      * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'.
3264      *
3265      * @param {Number} start Left border of the approximation interval
3266      * @param {Number} end Right border of the approximation interval
3267      * @returns {Number} The sum of the areas of the rectangles.
3268      * @memberof JXG.Math.Numerics
3269      */
3270     riemannsum: function (f, n, type, start, end) {
3271         JXG.deprecated("Numerics.riemannsum()", "Numerics.riemann()[2]");
3272         return this.riemann(f, n, type, start, end)[2];
3273     },
3274 
3275     /**
3276      * Solve initial value problems numerically using <i>explicit</i> Runge-Kutta methods.
3277      * See {@link https://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
3278      * @param {object|String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
3279      * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
3280      * <pre>
3281      * {
3282      *     s: <Number>,
3283      *     A: <matrix>,
3284      *     b: <Array>,
3285      *     c: <Array>
3286      * }
3287      * </pre>
3288      * which corresponds to the Butcher tableau structure
3289      * shown here: https://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696 .
3290      * <i>Default</i> is 'euler'.
3291      * @param {Array} x0 Initial value vector. Even if the problem is one-dimensional, the initial value has to be given in an array.
3292      * @param {Array} I Interval on which to integrate.
3293      * @param {Number} N Number of integration intervals, i.e. there are <i>N+1</i> evaluation points.
3294      * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
3295      * 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
3296      * vector <tt>x</tt>, and has to return a vector of the same length as <tt>x</tt> has.
3297      * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
3298      * @example
3299      * // A very simple autonomous system dx(t)/dt = x(t);
3300      * var f = function(t, x) {
3301      *     return [x[0]];
3302      * }
3303      *
3304      * // Solve it with initial value x(0) = 1 on the interval [0, 2]
3305      * // with 20 evaluation points.
3306      * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
3307      *
3308      * // Prepare data for plotting the solution of the ode using a curve.
3309      * var dataX = [];
3310      * var dataY = [];
3311      * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
3312      * var i;
3313      * for(i=0; i<data.length; i++) {
3314      *     dataX[i] = i*h;
3315      *     dataY[i] = data[i][0];
3316      * }
3317      * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
3318      * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
3319      * <script type="text/javascript">
3320      * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
3321      * var f = function(t, x) {
3322      *     // we have to copy the value.
3323      *     // return x; would just return the reference.
3324      *     return [x[0]];
3325      * }
3326      * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
3327      * var dataX = [];
3328      * var dataY = [];
3329      * var h = 0.1;
3330      * for(var i=0; i<data.length; i++) {
3331      *     dataX[i] = i*h;
3332      *     dataY[i] = data[i][0];
3333      * }
3334      * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
3335      * </script><pre>
3336      * @memberof JXG.Math.Numerics
3337      */
3338     rungeKutta: function (butcher, x0, I, N, f) {
3339         var e,
3340             i, j, k, l, s,
3341             x = [],
3342             y = [],
3343             h = (I[1] - I[0]) / N,
3344             t = I[0],
3345             dim = x0.length,
3346             result = [],
3347             r = 0;
3348 
3349         if (Type.isString(butcher)) {
3350             butcher = predefinedButcher[butcher] || predefinedButcher.euler;
3351         }
3352         s = butcher.s;
3353 
3354         // Don't change x0, so copy it
3355         x = x0.slice();
3356 
3357         for (i = 0; i <= N; i++) {
3358             result[r] = x.slice();
3359 
3360             r++;
3361             k = [];
3362 
3363             for (j = 0; j < s; j++) {
3364                 // Init y = 0
3365                 for (e = 0; e < dim; e++) {
3366                     y[e] = 0.0;
3367                 }
3368 
3369                 // Calculate linear combination of former k's and save it in y
3370                 for (l = 0; l < j; l++) {
3371                     for (e = 0; e < dim; e++) {
3372                         y[e] += butcher.A[j][l] * h * k[l][e];
3373                     }
3374                 }
3375 
3376                 // Add x(t) to y
3377                 for (e = 0; e < dim; e++) {
3378                     y[e] += x[e];
3379                 }
3380 
3381                 // Calculate new k and add it to the k matrix
3382                 k.push(f(t + butcher.c[j] * h, y));
3383             }
3384 
3385             // Init y = 0
3386             for (e = 0; e < dim; e++) {
3387                 y[e] = 0.0;
3388             }
3389 
3390             for (l = 0; l < s; l++) {
3391                 for (e = 0; e < dim; e++) {
3392                     y[e] += butcher.b[l] * k[l][e];
3393                 }
3394             }
3395 
3396             for (e = 0; e < dim; e++) {
3397                 x[e] = x[e] + h * y[e];
3398             }
3399 
3400             t += h;
3401         }
3402 
3403         return result;
3404     },
3405 
3406     /**
3407      * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and
3408      * {@link JXG.Math.Numerics.chandrupatla}
3409      * @type Number
3410      * @default 80
3411      * @memberof JXG.Math.Numerics
3412      */
3413     maxIterationsRoot: 80,
3414 
3415     /**
3416      * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr}
3417      * @type Number
3418      * @default 500
3419      * @memberof JXG.Math.Numerics
3420      */
3421     maxIterationsMinimize: 500,
3422 
3423     /**
3424      * Given a number x_0, this function tries to find a second number x_1 such that
3425      * the function f has opposite signs at x_0 and x_1.
3426      * The return values have to be tested if the method succeeded.
3427      *
3428      * @param {Function} f Function, whose root is to be found
3429      * @param {Number} x0 Start value
3430      * @param {Object} [context] Parent object in case f is method of it
3431      * @returns {Array} [x_0, f(x_0), x_1, f(x_1)] in case that x_0 <= x_1
3432      *   or [x_1, f(x_1), x_0, f(x_0)] in case that x_1 < x_0.
3433      *
3434      * @see JXG.Math.Numerics.fzero
3435      * @see JXG.Math.Numerics.chandrupatla
3436      *
3437      * @memberof JXG.Math.Numerics
3438      */
3439     findBracket: function (f, x0, context) {
3440         var a, aa, fa, blist, b, fb, u, fu, i, len;
3441 
3442         if (Type.isArray(x0)) {
3443             return x0;
3444         }
3445 
3446         a = x0;
3447         fa = f.call(context, a);
3448         // nfev += 1;
3449 
3450         // Try to get b, by trying several values related to a
3451         aa = a === 0 ? 1 : a;
3452         blist = [
3453             a - 0.1 * aa,
3454             a + 0.1 * aa,
3455             a - 1,
3456             a + 1,
3457             a - 0.5 * aa,
3458             a + 0.5 * aa,
3459             a - 0.6 * aa,
3460             a + 0.6 * aa,
3461             a - 1 * aa,
3462             a + 1 * aa,
3463             a - 2 * aa,
3464             a + 2 * aa,
3465             a - 5 * aa,
3466             a + 5 * aa,
3467             a - 10 * aa,
3468             a + 10 * aa,
3469             a - 50 * aa,
3470             a + 50 * aa,
3471             a - 100 * aa,
3472             a + 100 * aa
3473         ];
3474         len = blist.length;
3475 
3476         for (i = 0; i < len; i++) {
3477             b = blist[i];
3478             fb = f.call(context, b);
3479             // nfev += 1;
3480 
3481             if (fa * fb <= 0) {
3482                 break;
3483             }
3484         }
3485         if (b < a) {
3486             u = a;
3487             a = b;
3488             b = u;
3489 
3490             fu = fa;
3491             fa = fb;
3492             fb = fu;
3493         }
3494         return [a, fa, b, fb];
3495     },
3496 
3497     /**
3498      *
3499      * Find zero of an univariate function f.
3500      * @param {function} f Function, whose root is to be found
3501      * @param {Array|Number} x0  Start value or start interval enclosing the root.
3502      * If x0 is an interval [a,b], it is required that f(a)f(b) <= 0, otherwise the minimum of f in [a, b] will be returned.
3503      * If x0 is a number, the algorithms tries to enclose the root by an interval [a, b] containing x0 and the root and
3504      * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method.
3505      * @param {Object} [context] Parent object in case f is method of it
3506      * @returns {Number} the approximation of the root
3507      * Algorithm:
3508      *  Brent's root finder from
3509      *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
3510      *  computations. M., Mir, 1980, p.180 of the Russian edition
3511      *  https://www.netlib.org/c/brent.shar
3512      *
3513      * If x0 is an array containing lower and upper bound for the zero
3514      * algorithm 748 is applied. Otherwise, if x0 is a number,
3515      * the algorithm tries to bracket a zero of f starting from x0.
3516      * If this fails, we fall back to Newton's method.
3517      *
3518      * @see JXG.Math.Numerics.chandrupatla
3519      * @see JXG.Math.Numerics.root
3520      * @see JXG.Math.Numerics.findBracket
3521      * @see JXG.Math.Numerics.Newton
3522      * @see JXG.Math.Numerics.fminbr
3523      * @memberof JXG.Math.Numerics
3524      */
3525     fzero: function (f, x0, context) {
3526         var a, b, c,
3527             fa, fb, fc,
3528             res, x00,
3529             prev_step,
3530             t1, t2,
3531             cb,
3532             tol_act,   // Actual tolerance
3533             p,         // Interpolation step is calculated in the form p/q; division
3534             q,         // operations is delayed until the last moment
3535             new_step,  // Step at this iteration
3536             eps = Mat.eps,
3537             maxiter = this.maxIterationsRoot,
3538             niter = 0;
3539         // nfev = 0;
3540 
3541         if (Type.isArray(x0)) {
3542             if (x0.length < 2) {
3543                 throw new Error(
3544                     "JXG.Math.Numerics.fzero: length of array x0 has to be at least two."
3545                 );
3546             }
3547 
3548             x00 = this.findDomain(f, x0, context);
3549             a = x00[0];
3550             b = x00[1];
3551             // a = x0[0];
3552             // b = x0[1];
3553 
3554             fa = f.call(context, a);
3555             // nfev += 1;
3556             fb = f.call(context, b);
3557             // nfev += 1;
3558         } else {
3559             res = this.findBracket(f, x0, context);
3560             a = res[0];
3561             fa = res[1];
3562             b = res[2];
3563             fb = res[3];
3564         }
3565 
3566         if (Math.abs(fa) <= eps) {
3567             return a;
3568         }
3569         if (Math.abs(fb) <= eps) {
3570             return b;
3571         }
3572 
3573         if (fa * fb > 0) {
3574             // Bracketing not successful, fall back to Newton's method or to fminbr
3575             if (Type.isArray(x0)) {
3576                 return this.fminbr(f, [a, b], context);
3577             }
3578 
3579             return this.Newton(f, a, context);
3580         }
3581 
3582         // OK, we have enclosed a zero of f.
3583         // Now we can start Brent's method
3584         c = a;
3585         fc = fa;
3586 
3587         // Main iteration loop
3588         while (niter < maxiter) {
3589             // Distance from the last but one to the last approximation
3590             prev_step = b - a;
3591 
3592             // Swap data for b to be the best approximation
3593             if (Math.abs(fc) < Math.abs(fb)) {
3594                 a = b;
3595                 b = c;
3596                 c = a;
3597 
3598                 fa = fb;
3599                 fb = fc;
3600                 fc = fa;
3601             }
3602             tol_act = 2 * eps * Math.abs(b) + eps * 0.5;
3603             new_step = (c - b) * 0.5;
3604 
3605             if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) {
3606                 //  Acceptable approx. is found
3607                 return b;
3608             }
3609 
3610             // Decide if the interpolation can be tried
3611             // If prev_step was large enough and was in true direction Interpolatiom may be tried
3612             if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
3613                 cb = c - b;
3614 
3615                 // If we have only two distinct points linear interpolation can only be applied
3616                 if (a === c) {
3617                     t1 = fb / fa;
3618                     p = cb * t1;
3619                     q = 1.0 - t1;
3620                     // Quadric inverse interpolation
3621                 } else {
3622                     q = fa / fc;
3623                     t1 = fb / fc;
3624                     t2 = fb / fa;
3625 
3626                     p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
3627                     q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
3628                 }
3629 
3630                 // p was calculated with the opposite sign; make p positive
3631                 if (p > 0) {
3632                     q = -q;
3633                     // and assign possible minus to q
3634                 } else {
3635                     p = -p;
3636                 }
3637 
3638                 // If b+p/q falls in [b,c] and isn't too large it is accepted
3639                 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
3640                 if (
3641                     p < 0.75 * cb * q - Math.abs(tol_act * q) * 0.5 &&
3642                     p < Math.abs(prev_step * q * 0.5)
3643                 ) {
3644                     new_step = p / q;
3645                 }
3646             }
3647 
3648             // Adjust the step to be not less than tolerance
3649             if (Math.abs(new_step) < tol_act) {
3650                 new_step = new_step > 0 ? tol_act : -tol_act;
3651             }
3652 
3653             // Save the previous approx.
3654             a = b;
3655             fa = fb;
3656             b += new_step;
3657             fb = f.call(context, b);
3658             // Do step to a new approxim.
3659             // nfev += 1;
3660 
3661             // Adjust c for it to have a sign opposite to that of b
3662             if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
3663                 c = a;
3664                 fc = fa;
3665             }
3666             niter++;
3667         } // End while
3668 
3669         return b;
3670     },
3671 
3672     /**
3673      * Find zero of an univariate function f.
3674      * @param {function} f Function, whose root is to be found
3675      * @param {Array|Number} x0  Start value or start interval enclosing the root.
3676      * If x0 is an interval [a,b], it is required that f(a)f(b) <= 0, otherwise the minimum of f in [a, b] will be returned.
3677      * If x0 is a number, the algorithms tries to enclose the root by an interval [a, b] containing x0 and the root and
3678      * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method.
3679      * @param {Object} [context] Parent object in case f is method of it
3680      * @returns {Number} the approximation of the root
3681      * Algorithm:
3682      * Chandrupatla's method, see
3683      * Tirupathi R. Chandrupatla,
3684      * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives",
3685      * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149.
3686      *
3687      * If x0 is an array containing lower and upper bound for the zero
3688      * algorithm 748 is applied. Otherwise, if x0 is a number,
3689      * the algorithm tries to bracket a zero of f starting from x0.
3690      * If this fails, we fall back to Newton's method.
3691      *
3692      * @see JXG.Math.Numerics.root
3693      * @see JXG.Math.Numerics.fzero
3694      * @see JXG.Math.Numerics.findBracket
3695      * @see JXG.Math.Numerics.Newton
3696      * @see JXG.Math.Numerics.fminbr
3697      * @memberof JXG.Math.Numerics
3698      */
3699     chandrupatla: function (f, x0, context) {
3700         var a, b, fa, fb,
3701             res,
3702             niter = 0,
3703             maxiter = this.maxIterationsRoot,
3704             rand = 1 + Math.random() * 0.001,
3705             t = 0.5 * rand,
3706             eps = Mat.eps, // 1.e-10,
3707             dlt = 0.00001,
3708             x1, x2, x3, x,
3709             f1, f2, f3, y,
3710             xm, fm,
3711             tol, tl,
3712             xi, ph, fl, fh,
3713             AL, A, B, C, D;
3714 
3715         if (Type.isArray(x0)) {
3716             if (x0.length < 2) {
3717                 throw new Error(
3718                     "JXG.Math.Numerics.fzero: length of array x0 has to be at least two."
3719                 );
3720             }
3721 
3722             a = x0[0];
3723             fa = f.call(context, a);
3724             // nfev += 1;
3725             b = x0[1];
3726             fb = f.call(context, b);
3727             // nfev += 1;
3728         } else {
3729             res = this.findBracket(f, x0, context);
3730             a = res[0];
3731             fa = res[1];
3732             b = res[2];
3733             fb = res[3];
3734         }
3735 
3736         if (fa * fb > 0) {
3737             // Bracketing not successful, fall back to Newton's method or to fminbr
3738             if (Type.isArray(x0)) {
3739                 return this.fminbr(f, [a, b], context);
3740             }
3741 
3742             return this.Newton(f, a, context);
3743         }
3744 
3745         x1 = a;
3746         x2 = b;
3747         f1 = fa;
3748         f2 = fb;
3749         do {
3750             x = x1 + t * (x2 - x1);
3751             y = f.call(context, x);
3752 
3753             // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point
3754             if (Math.sign(y) === Math.sign(f1)) {
3755                 x3 = x1;
3756                 x1 = x;
3757                 f3 = f1;
3758                 f1 = y;
3759             } else {
3760                 x3 = x2;
3761                 x2 = x1;
3762                 f3 = f2;
3763                 f2 = f1;
3764             }
3765             x1 = x;
3766             f1 = y;
3767 
3768             xm = x1;
3769             fm = f1;
3770             if (Math.abs(f2) < Math.abs(f1)) {
3771                 xm = x2;
3772                 fm = f2;
3773             }
3774             tol = 2 * eps * Math.abs(xm) + 0.5 * dlt;
3775             tl = tol / Math.abs(x2 - x1);
3776             if (tl > 0.5 || fm === 0) {
3777                 break;
3778             }
3779             // If inverse quadratic interpolation holds, use it
3780             xi = (x1 - x2) / (x3 - x2);
3781             ph = (f1 - f2) / (f3 - f2);
3782             fl = 1 - Math.sqrt(1 - xi);
3783             fh = Math.sqrt(xi);
3784             if (fl < ph && ph < fh) {
3785                 AL = (x3 - x1) / (x2 - x1);
3786                 A = f1 / (f2 - f1);
3787                 B = f3 / (f2 - f3);
3788                 C = f1 / (f3 - f1);
3789                 D = f2 / (f3 - f2);
3790                 t = A * B + C * D * AL;
3791             } else {
3792                 t = 0.5 * rand;
3793             }
3794             // Adjust t away from the interval boundary
3795             if (t < tl) {
3796                 t = tl;
3797             }
3798             if (t > 1 - tl) {
3799                 t = 1 - tl;
3800             }
3801             niter++;
3802         } while (niter <= maxiter);
3803         // console.log(niter);
3804 
3805         return xm;
3806     },
3807 
3808     /**
3809      * Find a small enclosing interval of the domain of a function by
3810      * tightening the input interval x0.
3811      * <p>
3812      * This is a helper function which is used in {@link JXG.Math.Numerics.fminbr},
3813      * {@link JXG.Math.Numerics.fzero}, and  {@link JXG.Curve.getLabelPosition}
3814      * to avoid search in an interval where the function is mostly undefined.
3815      *
3816      * @param {function} f
3817      * @param {Array} x0 Start interval
3818      * @param {Object} context Parent object in case f is method of it
3819      * @param {Boolean} [outer=true] if true take a proper enclosing array. If false return the domain such that the function is defined
3820      * at its  borders.
3821      * @returns Array
3822      *
3823      * @example
3824      * var f = (x) => Math.sqrt(x);
3825      * console.log(JXG.Math.Numerics.findDomain(f, [-5, 5]));
3826      *
3827      * // Output: [ -0.00020428174445492973, 5 ]
3828      *
3829      * @example
3830      * var f = (x) => Math.sqrt(x);
3831      * console.log(JXG.Math.Numerics.findDomain(f, [-5, 5], null, false));
3832      *
3833      * // Output: [ 0.00020428174562965915, 5 ]
3834      */
3835     findDomain: function (f, x0, context, outer) {
3836         var a, b, c, fc,
3837             x,
3838             gr = 1 - 1 / 1.61803398875,
3839             eps = 0.001,
3840             cnt,
3841             max_cnt = 20;
3842 
3843         if (outer === undefined) {
3844             outer = true;
3845         }
3846 
3847         if (!Type.isArray(x0) || x0.length < 2) {
3848             throw new Error(
3849                 "JXG.Math.Numerics.findDomain: length of array x0 has to be at least two."
3850             );
3851         }
3852 
3853         x = x0.slice();
3854         a = x[0];
3855         b = x[1];
3856         fc = f.call(context, a);
3857         if (isNaN(fc)) {
3858             // Divide the interval with the golden ratio
3859             // and keep a such that f(a) = NaN
3860             cnt = 0;
3861             while (b - a > eps && cnt < max_cnt) {
3862                 c = (b - a) * gr + a;
3863                 fc = f.call(context, c);
3864                 if (isNaN(fc)) {
3865                     a = c;
3866                 } else {
3867                     b = c;
3868                 }
3869                 cnt++;
3870             }
3871             if (outer) {
3872                 x[0] = a;
3873             } else {
3874                 x[0] = b;
3875             }
3876             // x[0] = a;
3877         }
3878 
3879         a = x[0];
3880         b = x[1];
3881         fc = f.call(context, b);
3882         if (isNaN(fc)) {
3883             // Divide the interval with the golden ratio
3884             // and keep b such that f(b) = NaN
3885             cnt = 0;
3886             while (b - a > eps && cnt < max_cnt) {
3887                 c = b - (b - a) * gr;
3888                 fc = f.call(context, c);
3889                 if (isNaN(fc)) {
3890                     b = c;
3891                 } else {
3892                     a = c;
3893                 }
3894                 cnt++;
3895             }
3896             if (outer) {
3897                 x[1] = b;
3898             } else {
3899                 x[1] = a;
3900             }
3901             // x[1] = b;
3902         }
3903 
3904         return x;
3905     },
3906 
3907     /**
3908      *
3909      * Find minimum of an univariate function f.
3910      * <p>
3911      * Algorithm:
3912      *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
3913      *  computations. M., Mir, 1980, p.180 of the Russian edition
3914      *
3915      * @param {function} f Function, whose minimum is to be found
3916      * @param {Array} x0  Start interval enclosing the minimum
3917      * @param {Object} [context] Parent object in case f is method of it
3918      * @returns {Number} the approximation of the minimum value position
3919      * @memberof JXG.Math.Numerics
3920      **/
3921     fminbr: function (f, x0, context) {
3922         var a, b, x, v, w,
3923             fx, fv, fw,
3924             x00,
3925             range, middle_range, tol_act, new_step,
3926             p, q, t, ft,
3927             r = (3.0 - Math.sqrt(5.0)) * 0.5,      // Golden section ratio
3928             tol = Mat.eps,
3929             sqrteps = Mat.eps, // Math.sqrt(Mat.eps),
3930             maxiter = this.maxIterationsMinimize,
3931             niter = 0;
3932         // nfev = 0;
3933 
3934         if (!Type.isArray(x0) || x0.length < 2) {
3935             throw new Error(
3936                 "JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."
3937             );
3938         }
3939 
3940         x00 = this.findDomain(f, x0, context);
3941         a = x00[0];
3942         b = x00[1];
3943         v = a + r * (b - a);
3944         fv = f.call(context, v);
3945 
3946         // First step - always gold section
3947         // nfev += 1;
3948         x = v;
3949         w = v;
3950         fx = fv;
3951         fw = fv;
3952 
3953         while (niter < maxiter) {
3954             // Range over the interval in which we are looking for the minimum
3955             range = b - a;
3956             middle_range = (a + b) * 0.5;
3957 
3958             // Actual tolerance
3959             tol_act = sqrteps * Math.abs(x) + tol / 3.0;
3960 
3961             if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) {
3962                 // Acceptable approx. is found
3963                 return x;
3964             }
3965 
3966             // Obtain the golden section step
3967             new_step = r * (x < middle_range ? b - x : a - x);
3968 
3969             // Decide if the interpolation can be tried. If x and w are distinct, interpolatiom may be tried
3970             if (Math.abs(x - w) >= tol_act) {
3971                 // Interpolation step is calculated as p/q;
3972                 // division operation is delayed until last moment
3973                 t = (x - w) * (fx - fv);
3974                 q = (x - v) * (fx - fw);
3975                 p = (x - v) * q - (x - w) * t;
3976                 q = 2 * (q - t);
3977 
3978                 if (q > 0) {
3979                     p = -p; // q was calculated with the opposite sign; make q positive
3980                 } else {
3981                     q = -q; // // and assign possible minus to p
3982                 }
3983                 if (
3984                     Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b]
3985                     p > q * (a - x + 2 * tol_act) &&        //  not too close to a and
3986                     p < q * (b - x - 2 * tol_act)
3987                 ) {
3988                     // b, and isn't too large
3989                     new_step = p / q; // it is accepted
3990                 }
3991                 // If p / q is too large then the
3992                 // golden section procedure can
3993                 // reduce [a,b] range to more
3994                 // extent
3995             }
3996 
3997             // Adjust the step to be not less than tolerance
3998             if (Math.abs(new_step) < tol_act) {
3999                 if (new_step > 0) {
4000                     new_step = tol_act;
4001                 } else {
4002                     new_step = -tol_act;
4003                 }
4004             }
4005 
4006             // Obtain the next approximation to min
4007             // and reduce the enveloping range
4008 
4009             // Tentative point for the min
4010             t = x + new_step;
4011             ft = f.call(context, t);
4012             // nfev += 1;
4013 
4014             // t is a better approximation
4015             if (ft <= fx) {
4016                 // Reduce the range so that t would fall within it
4017                 if (t < x) {
4018                     b = x;
4019                 } else {
4020                     a = x;
4021                 }
4022 
4023                 // Assign the best approx to x
4024                 v = w;
4025                 w = x;
4026                 x = t;
4027 
4028                 fv = fw;
4029                 fw = fx;
4030                 fx = ft;
4031                 // x remains the better approx
4032             } else {
4033                 // Reduce the range enclosing x
4034                 if (t < x) {
4035                     a = t;
4036                 } else {
4037                     b = t;
4038                 }
4039 
4040                 if (ft <= fw || w === x) {
4041                     v = w;
4042                     w = t;
4043                     fv = fw;
4044                     fw = ft;
4045                 } else if (ft <= fv || v === x || v === w) {
4046                     v = t;
4047                     fv = ft;
4048                 }
4049             }
4050             niter += 1;
4051         }
4052 
4053         return x;
4054     },
4055 
4056     /**
4057      * GLOMIN seeks a global minimum of a function F(X) in an interval [A,B]
4058      * and is the adaption of the algorithm GLOMIN by Richard Brent.
4059      *
4060      * Here is the original documentation:
4061      * <pre>
4062      *
4063      * Discussion:
4064      *
4065      * This function assumes that F(X) is twice continuously differentiable over [A,B]
4066      * and that |F''(X)| <= M for all X in [A,B].
4067      *
4068      * Licensing:
4069      *   This code is distributed under the GNU LGPL license.
4070      *
4071      * Modified:
4072      *
4073      *   17 April 2008
4074      *
4075      * Author:
4076      *
4077      *   Original FORTRAN77 version by Richard Brent.
4078      *   C version by John Burkardt.
4079      *   https://people.math.sc.edu/Burkardt/c_src/brent/brent.c
4080      *
4081      * Reference:
4082      *
4083      *   Richard Brent,
4084      *  Algorithms for Minimization Without Derivatives,
4085      *   Dover, 2002,
4086      *  ISBN: 0-486-41998-3,
4087      *   LC: QA402.5.B74.
4088      *
4089      * Parameters:
4090      *
4091      *   Input, double A, B, the endpoints of the interval.
4092      *  It must be the case that A < B.
4093      *
4094      *   Input, double C, an initial guess for the global
4095      *  minimizer.  If no good guess is known, C = A or B is acceptable.
4096      *
4097      *  Input, double M, the bound on the second derivative.
4098      *
4099      *   Input, double MACHEP, an estimate for the relative machine
4100      *  precision.
4101      *
4102      *   Input, double E, a positive tolerance, a bound for the
4103      *  absolute error in the evaluation of F(X) for any X in [A,B].
4104      *
4105      *   Input, double T, a positive error tolerance.
4106      *
4107      *    Input, double F (double x ), a user-supplied
4108      *  function whose global minimum is being sought.
4109      *
4110      *   Output, double *X, the estimated value of the abscissa
4111      *  for which F attains its global minimum value in [A,B].
4112      *
4113      *   Output, double GLOMIN, the value F(X).
4114      * </pre>
4115      *
4116      * In JSXGraph, some parameters of the original algorithm are set to fixed values:
4117      * <ul>
4118      *  <li> M = 10000000.0
4119      *  <li> C = A or B, depending if f(A) <= f(B)
4120      *  <li> T = JXG.Math.eps
4121      *  <li> E = JXG.Math.eps * JXG.Math.eps
4122      *  <li> MACHEP = JXG.Math.eps * JXG.Math.eps * JXG.Math.eps
4123      * </ul>
4124      * @param {function} f Function, whose global minimum is to be found
4125      * @param {Array} x0 Array of length 2 determining the interval [A, B] for which the global minimum is to be found
4126      * @returns {Array} [x, y] x is the position of the global minimum and y = f(x).
4127      */
4128     glomin: function (f, x0) {
4129         var a0, a2, a3, d0, d1, d2, h,
4130             k, m2,
4131             p, q, qs,
4132             r, s, sc,
4133             y, y0, y1, y2, y3, yb,
4134             z0, z1, z2,
4135             a, b, c, x,
4136             m = 10000000.0,
4137             t = Mat.eps, // * Mat.eps,
4138             e = Mat.eps * Mat.eps,
4139             machep = Mat.eps * Mat.eps * Mat.eps;
4140 
4141         a = x0[0];
4142         b = x0[1];
4143         c = (f(a) < f(b)) ? a : b;
4144 
4145         a0 = b;
4146         x = a0;
4147         a2 = a;
4148         y0 = f(b);
4149         yb = y0;
4150         y2 = f(a);
4151         y = y2;
4152 
4153         if (y0 < y) {
4154             y = y0;
4155         } else {
4156             x = a;
4157         }
4158 
4159         if (m <= 0.0 || b <= a) {
4160             return y;
4161         }
4162 
4163         m2 = 0.5 * (1.0 + 16.0 * machep) * m;
4164 
4165         if (c <= a || b <= c) {
4166             sc = 0.5 * (a + b);
4167         } else {
4168             sc = c;
4169         }
4170 
4171         y1 = f(sc);
4172         k = 3;
4173         d0 = a2 - sc;
4174         h = 9.0 / 11.0;
4175 
4176         if (y1 < y) {
4177             x = sc;
4178             y = y1;
4179         }
4180 
4181         for (; ;) {
4182             d1 = a2 - a0;
4183             d2 = sc - a0;
4184             z2 = b - a2;
4185             z0 = y2 - y1;
4186             z1 = y2 - y0;
4187             r = d1 * d1 * z0 - d0 * d0 * z1;
4188             p = r;
4189             qs = 2.0 * (d0 * z1 - d1 * z0);
4190             q = qs;
4191 
4192             if (k < 1000000 || y2 <= y) {
4193                 for (; ;) {
4194                     if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) <
4195                         z2 * m2 * r * (z2 * q - r)) {
4196 
4197                         a3 = a2 + r / q;
4198                         y3 = f(a3);
4199 
4200                         if (y3 < y) {
4201                             x = a3;
4202                             y = y3;
4203                         }
4204                     }
4205                     k = ((1611 * k) % 1048576);
4206                     q = 1.0;
4207                     r = (b - a) * 0.00001 * k;
4208 
4209                     if (z2 <= r) {
4210                         break;
4211                     }
4212                 }
4213             } else {
4214                 k = ((1611 * k) % 1048576);
4215                 q = 1.0;
4216                 r = (b - a) * 0.00001 * k;
4217 
4218                 while (r < z2) {
4219                     if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) <
4220                         z2 * m2 * r * (z2 * q - r)) {
4221 
4222                         a3 = a2 + r / q;
4223                         y3 = f(a3);
4224 
4225                         if (y3 < y) {
4226                             x = a3;
4227                             y = y3;
4228                         }
4229                     }
4230                     k = ((1611 * k) % 1048576);
4231                     q = 1.0;
4232                     r = (b - a) * 0.00001 * k;
4233                 }
4234             }
4235 
4236             r = m2 * d0 * d1 * d2;
4237             s = Math.sqrt(((y2 - y) + t) / m2);
4238             h = 0.5 * (1.0 + h);
4239             p = h * (p + 2.0 * r * s);
4240             q = q + 0.5 * qs;
4241             r = - 0.5 * (d0 + (z0 + 2.01 * e) / (d0 * m2));
4242 
4243             if (r < s || d0 < 0.0) {
4244                 r = a2 + s;
4245             } else {
4246                 r = a2 + r;
4247             }
4248 
4249             if (0.0 < p * q) {
4250                 a3 = a2 + p / q;
4251             } else {
4252                 a3 = r;
4253             }
4254 
4255             for (; ;) {
4256                 a3 = Math.max(a3, r);
4257 
4258                 if (b <= a3) {
4259                     a3 = b;
4260                     y3 = yb;
4261                 } else {
4262                     y3 = f(a3);
4263                 }
4264 
4265                 if (y3 < y) {
4266                     x = a3;
4267                     y = y3;
4268                 }
4269 
4270                 d0 = a3 - a2;
4271 
4272                 if (a3 <= r) {
4273                     break;
4274                 }
4275 
4276                 p = 2.0 * (y2 - y3) / (m * d0);
4277 
4278                 if ((1.0 + 9.0 * machep) * d0 <= Math.abs(p)) {
4279                     break;
4280                 }
4281 
4282                 if (0.5 * m2 * (d0 * d0 + p * p) <= (y2 - y) + (y3 - y) + 2.0 * t) {
4283                     break;
4284                 }
4285                 a3 = 0.5 * (a2 + a3);
4286                 h = 0.9 * h;
4287             }
4288 
4289             if (b <= a3) {
4290                 break;
4291             }
4292 
4293             a0 = sc;
4294             sc = a2;
4295             a2 = a3;
4296             y0 = y1;
4297             y1 = y2;
4298             y2 = y3;
4299         }
4300 
4301         return [x, y];
4302     },
4303 
4304     /**
4305      * Determine all roots of a polynomial with real or complex coefficients by using the
4306      * iterative method attributed to Weierstrass, Durand, Kerner, Aberth, and Ehrlich. In particular,
4307      * the iteration method with cubic convergence is used that is usually attributed to Ehrlich-Aberth.
4308      * <p>
4309      * The returned roots are sorted with respect to their real values.
4310      * <p> This method makes use of the JSXGraph classes {@link JXG.Complex} and {@link JXG.C} to handle
4311      * complex numbers.
4312      *
4313      * @param {Array} a Array of coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4314      * The coefficients are of type Number or JXG.Complex.
4315      * @param {Number} [deg] Optional degree of the polynomial. Otherwise all entries are taken, with
4316      * leading zeros removed.
4317      * @param {Number} [tol=Number.EPSILON] Approximation tolerance
4318      * @param {Number} [max_it=30] Maximum number of iterations
4319      * @param {Array} [initial_values=null] Array of initial values for the roots. If not given,
4320      * starting values are determined by the method of Ozawa.
4321      * @returns {Array} Array of complex numbers (of JXG.Complex) approximating the roots of the polynomial.
4322      * @memberof JXG.Math.Numerics
4323      * @see JXG.Complex
4324      * @see JXG.C
4325      *
4326      * @example
4327      * // Polynomial p(z) = -1 + 1z^2
4328      * var i, roots,
4329      *     p = [-1, 0, 1];
4330      *
4331      * roots = JXG.Math.Numerics.polzeros(p);
4332      * for (i = 0; i < roots.length; i++) {
4333      *     console.log(i, roots[i].toString());
4334      * }
4335      * // Output:
4336      *   0 -1 + -3.308722450212111e-24i
4337      *   1 1 + 0i
4338      *
4339      * @example
4340      * // Polynomial p(z) = -1 + 3z - 9z^2 + z^3 - 8z^6 + 9z^7 - 9z^8 + z^9
4341      * var i, roots,
4342      *     p = [-1, 3, -9, 1, 0, 0, -8, 9, -9, 1];
4343      *
4344      * roots = JXG.Math.Numerics.polzeros(p);
4345      * for (i = 0; i < roots.length; i++) {
4346      *     console.log(i, roots[i].toString());
4347      * }
4348      * // Output:
4349      * 0 -0.7424155888401961 + 0.4950476539211721i
4350      * 1 -0.7424155888401961 + -0.4950476539211721i
4351      * 2 0.16674869833354108 + 0.2980502714610669i
4352      * 3 0.16674869833354108 + -0.29805027146106694i
4353      * 4 0.21429002063640837 + 1.0682775088132996i
4354      * 5 0.21429002063640842 + -1.0682775088132999i
4355      * 6 0.861375497926218 + -0.6259177003583295i
4356      * 7 0.8613754979262181 + 0.6259177003583295i
4357      * 8 8.000002743888055 + -1.8367099231598242e-40i
4358      *
4359      */
4360     polzeros: function (coeffs, deg, tol, max_it, initial_values) {
4361         var i, le, off, it,
4362             debug = false,
4363             cc = [],
4364             obvious = [],
4365             roots = [],
4366 
4367             /**
4368              * Horner method to evaluate polynomial or the derivative thereof for complex numbers,
4369              * i.e. coefficients and variable are complex.
4370              * @function
4371              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4372              * @param {JXG.Complex} z Value for which the polynomial will be evaluated.
4373              * @param {Boolean} [derivative=false] If true the derivative will be evaluated.
4374              * @ignore
4375              */
4376             hornerComplex = function (a, z, derivative) {
4377                 var i, s,
4378                     n = a.length - 1;
4379 
4380                 derivative = derivative || false;
4381                 if (derivative) {
4382                     // s = n * a_n
4383                     s = JXG.C.mult(n, a[n]);
4384                     for (i = n - 1; i > 0; i--) {
4385                         // s = s * z + i * a_i
4386                         s.mult(z);
4387                         s.add(JXG.C.mult(a[i], i));
4388                     }
4389                 } else {
4390                     // s = a_n
4391                     s = JXG.C.copy(a[n]);
4392                     for (i = n - 1; i >= 0; i--) {
4393                         // s = s * z + a_i
4394                         s.mult(z);
4395                         s.add(a[i]);
4396                     }
4397                 }
4398                 return s;
4399             },
4400 
4401             /**
4402              * Horner method to evaluate reciprocal polynomial or the derivative thereof for complex numbers,
4403              * i.e. coefficients and variable are complex.
4404              * @function
4405              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4406              * @param {JXG.Complex} z Value for which the reciprocal polynomial will be evaluated.
4407              * @param {Boolean} [derivative=false] If true the derivative will be evaluated.
4408              * @ignore
4409              */
4410             hornerRec = function (a, x, derivative) {
4411                 var i, s,
4412                     n = a.length - 1;
4413 
4414                 derivative = derivative || false;
4415                 if (derivative) {
4416                     // s = n * a_0
4417                     s = JXG.C.mult(n, a[0]);
4418                     for (i = n - 1; i > 0; i--) {
4419                         // s = s * x + i * a_{n-i}
4420                         s.mult(x);
4421                         s.add(JXG.C.mult(a[n - i], i));
4422                     }
4423                 } else {
4424                     // s = a_0
4425                     s = JXG.C.copy(a[0]);
4426                     for (i = n - 1; i >= 0; i--) {
4427                         // s = s * x + a_{n-i}
4428                         s.mult(x);
4429                         s.add(a[n - i]);
4430                     }
4431                 }
4432                 return s;
4433             },
4434 
4435             /**
4436              * Horner method to evaluate real polynomial at a real value.
4437              * @function
4438              * @param {Array} a Array of real coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4439              * @param {Number} z Value for which the polynomial will be evaluated.
4440              * @ignore
4441              */
4442             horner = function (a, x) {
4443                 var i, s,
4444                     n = a.length - 1;
4445 
4446                 s = a[n];
4447                 for (i = n - 1; i >= 0; i--) {
4448                     s = s * x + a[i];
4449                 }
4450                 return s;
4451             },
4452 
4453             /**
4454              * Determine start values for the Aberth iteration, see
4455              * Ozawa, "An experimental study of the starting values
4456              * of the Durand-Kerner-Aberth iteration" (1995).
4457              *
4458              * @function
4459              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4460              * @returns {Array} Array Initial values for the roots.
4461              * @ignore
4462              */
4463             initial_guess = function (a) {
4464                 var i, r,
4465                     n = a.length - 1, // degree
4466                     alpha1 = Math.PI * 2 / n,
4467                     alpha0 = Math.PI / n * 0.5,
4468                     b, z,
4469                     init = [];
4470 
4471 
4472                 // From Ozawa, "An experimental study of the starting values
4473                 // of the Durand-Kerner-Aberth iteration" (1995)
4474 
4475                 // b is the arithmetic mean of the roots.
4476                 // With is Vieta's formula <https://en.wikipedia.org/wiki/Vieta%27s_formulas>
4477                 //   b = -a_{n-1} / (n * a_n)
4478                 b = JXG.C.mult(-1, a[n - 1]);
4479                 b.div(JXG.C.mult(n, a[n]));
4480 
4481                 // r is the geometric mean of the deviations |b - root_i|.
4482                 // Using
4483                 //   p(z) = a_n prod(z - root_i)
4484                 // and therefore
4485                 //   |p(b)| = |a_n| prod(|b - root_i|)
4486                 // we arrive at:
4487                 //   r = |p(b)/a_n|^(1/n)
4488                 z = JXG.C.div(hornerComplex(a, b), a[n]);
4489                 r = Math.pow(JXG.C.abs(z), 1 / n);
4490                 if (r === 0) { r = 1; }
4491 
4492                 for (i = 0; i < n; i++) {
4493                     a = new JXG.Complex(r * Math.cos(alpha1 * i + alpha0), r * Math.sin(alpha1 * i + alpha0));
4494                     init[i] = JXG.C.add(b, a);
4495                 }
4496 
4497                 return init;
4498             },
4499 
4500             /**
4501              * Ehrlich-Aberth iteration. The stopping criterion is from
4502              * D.A. Bini, "Numerical computation of polynomial zeros
4503              * by means of Aberths's method", Numerical Algorithms (1996).
4504              *
4505              * @function
4506              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4507              * @param {Number} mu Machine precision
4508              * @param {Number} max_it Maximum number of iterations
4509              * @param {Array} z Initial guess for the roots. Will be changed in place.
4510              * @returns {Number} Number of iterations
4511              * @ignore
4512              */
4513             aberthIteration = function (cc, mu, max_it, z) {
4514                 var k, i, j,
4515                     done = [],
4516                     cr = [],
4517                     gamma, x,
4518                     done_sum = 0,
4519                     num, denom, s, pp,
4520                     n = z.length;
4521 
4522                 for (i = 0; i < n; i++) {
4523                     done.push(false);
4524                 }
4525                 for (i = 0; i < cc.length; i++) {
4526                     cr.push(JXG.C.abs(cc[i]) * (4 * i + 1));
4527                 }
4528                 for (k = 0; k < max_it && done_sum < n; k++) {
4529                     for (i = 0; i < n; i++) {
4530                         if (done[i]) {
4531                             continue;
4532                         }
4533                         num = hornerComplex(cc, z[i]);
4534                         x = JXG.C.abs(z[i]);
4535 
4536                         // Stopping criterion by D.A. Bini
4537                         // "Numerical computation of polynomial zeros
4538                         // by means of Aberths's method", Numerical Algorithms (1996).
4539                         //
4540                         if (JXG.C.abs(num) < mu * horner(cr, x)) {
4541                             done[i] = true;
4542                             done_sum++;
4543                             if (done_sum === n) {
4544                                 break;
4545                             }
4546                             continue;
4547                         }
4548 
4549                         // num = P(z_i) / P'(z_i)
4550                         if (x > 1) {
4551                             gamma = JXG.C.div(1, z[i]);
4552                             pp = hornerRec(cc, gamma, true);
4553                             pp.div(hornerRec(cc, gamma));
4554                             pp.mult(gamma);
4555                             num = JXG.C.sub(n, pp);
4556                             num = JXG.C.div(z[i], num);
4557                         } else {
4558                             num.div(hornerComplex(cc, z[i], true));
4559                         }
4560 
4561                         // denom = sum_{i\neq j} 1 / (z_i  - z_j)
4562                         denom = new JXG.Complex(0);
4563                         for (j = 0; j < n; j++) {
4564                             if (j === i) {
4565                                 continue;
4566                             }
4567                             s = JXG.C.sub(z[i], z[j]);
4568                             s = JXG.C.div(1, s);
4569                             denom.add(s);
4570                         }
4571 
4572                         // num = num / 1 - num * sum_{i\neq j} 1 / (z_i - z_j)
4573                         denom.mult(num);
4574                         denom = JXG.C.sub(1, denom);
4575                         num.div(denom);
4576                         // z_i = z_i - num
4577                         z[i].sub(num);
4578                     }
4579                 }
4580 
4581                 return k;
4582             };
4583 
4584 
4585         tol = tol || Number.EPSILON;
4586         max_it = max_it || 30;
4587 
4588         le = coeffs.length;
4589         if (JXG.isNumber(deg) && deg >= 0 && deg < le - 1) {
4590             le = deg + 1;
4591         }
4592 
4593         // Convert coefficient array to complex numbers
4594         for (i = 0; i < le; i++) {
4595             cc.push(new JXG.Complex(coeffs[i]));
4596         }
4597 
4598         // Search for (multiple) roots at x=0
4599         for (i = 0; i < le; i++) {
4600             if (cc[i].real !== 0 || cc[i].imaginary !== 0) {
4601                 off = i;
4602                 break;
4603             }
4604         }
4605 
4606         // Deflate root x=0, store roots at x=0 in obvious
4607         for (i = 0; i < off; i++) {
4608             obvious.push(new JXG.Complex(0));
4609         }
4610         cc = cc.slice(off);
4611         le = cc.length;
4612 
4613         // Remove leading zeros from the coefficient array
4614         for (i = le - 1; i >= 0; i--) {
4615             if (cc[i].real !== 0 || cc[i].imaginary !== 0) {
4616                 break;
4617             }
4618             cc.pop();
4619         }
4620         le = cc.length;
4621         if (le === 0) {
4622             return [];
4623         }
4624 
4625         // From now on we can assume that the
4626         // constant coefficient and the leading coefficient
4627         // are not zero.
4628         if (initial_values) {
4629             for (i = 0; i < le - 1; i++) {
4630                 roots.push(new JXG.Complex(initial_values[i]));
4631             }
4632         } else {
4633             roots = initial_guess(cc);
4634         }
4635         it = aberthIteration(cc, tol, max_it, roots);
4636 
4637         // Append the roots at x=0
4638         roots = obvious.concat(roots);
4639 
4640         if (debug) {
4641             console.log("Iterations:", it);
4642             console.log('Roots:');
4643             for (i = 0; i < roots.length; i++) {
4644                 console.log(i, roots[i].toString(), JXG.C.abs(hornerComplex(cc, roots[i])));
4645             }
4646         }
4647 
4648         // Sort roots according to their real part
4649         roots.sort(function (a, b) {
4650             if (a.real < b.real) {
4651                 return -1;
4652             }
4653             if (a.real > b.real) {
4654                 return 1;
4655             }
4656             return 0;
4657         });
4658 
4659         return roots;
4660     },
4661 
4662     /**
4663      * Implements the Ramer-Douglas-Peucker algorithm.
4664      * It discards points which are not necessary from the polygonal line defined by the point array
4665      * pts. The computation is done in screen coordinates.
4666      * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points.
4667      * @param {Array} pts Array of {@link JXG.Coords}
4668      * @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>.
4669      * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points.
4670      * @memberof JXG.Math.Numerics
4671      */
4672     RamerDouglasPeucker: function (pts, eps) {
4673         var allPts = [],
4674             newPts = [],
4675             i, k, len,
4676             endless = true,
4677 
4678             /**
4679              * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
4680              * It searches for the point between index i and j which
4681              * has the largest distance from the line between the points i and j.
4682              * @param {Array} pts Array of {@link JXG.Coords}
4683              * @param {Number} i Index of a point in pts
4684              * @param {Number} j Index of a point in pts
4685              * @ignore
4686              * @private
4687              */
4688             findSplit = function (pts, i, j) {
4689                 var d, k, ci, cj, ck,
4690                     x0, y0, x1, y1,
4691                     den, lbda,
4692                     eps = Mat.eps * Mat.eps,
4693                     huge = 10000,
4694                     dist = 0,
4695                     f = i;
4696 
4697                 if (j - i < 2) {
4698                     return [-1.0, 0];
4699                 }
4700 
4701                 ci = pts[i].scrCoords;
4702                 cj = pts[j].scrCoords;
4703 
4704                 if (isNaN(ci[1]) || isNaN(ci[2])) {
4705                     return [NaN, i];
4706                 }
4707                 if (isNaN(cj[1]) || isNaN(cj[2])) {
4708                     return [NaN, j];
4709                 }
4710 
4711                 for (k = i + 1; k < j; k++) {
4712                     ck = pts[k].scrCoords;
4713                     if (isNaN(ck[1]) || isNaN(ck[2])) {
4714                         return [NaN, k];
4715                     }
4716 
4717                     x0 = ck[1] - ci[1];
4718                     y0 = ck[2] - ci[2];
4719                     x1 = cj[1] - ci[1];
4720                     y1 = cj[2] - ci[2];
4721                     x0 = x0 === Infinity ? huge : x0;
4722                     y0 = y0 === Infinity ? huge : y0;
4723                     x1 = x1 === Infinity ? huge : x1;
4724                     y1 = y1 === Infinity ? huge : y1;
4725                     x0 = x0 === -Infinity ? -huge : x0;
4726                     y0 = y0 === -Infinity ? -huge : y0;
4727                     x1 = x1 === -Infinity ? -huge : x1;
4728                     y1 = y1 === -Infinity ? -huge : y1;
4729                     den = x1 * x1 + y1 * y1;
4730 
4731                     if (den > eps) {
4732                         lbda = (x0 * x1 + y0 * y1) / den;
4733 
4734                         if (lbda < 0.0) {
4735                             lbda = 0.0;
4736                         } else if (lbda > 1.0) {
4737                             lbda = 1.0;
4738                         }
4739 
4740                         x0 = x0 - lbda * x1;
4741                         y0 = y0 - lbda * y1;
4742                         d = x0 * x0 + y0 * y0;
4743                     } else {
4744                         lbda = 0.0;
4745                         d = x0 * x0 + y0 * y0;
4746                     }
4747 
4748                     if (d > dist) {
4749                         dist = d;
4750                         f = k;
4751                     }
4752                 }
4753                 return [Math.sqrt(dist), f];
4754             },
4755             /**
4756              * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
4757              * It runs recursively through the point set and searches the
4758              * point which has the largest distance from the line between the first point and
4759              * the last point. If the distance from the line is greater than eps, this point is
4760              * included in our new point set otherwise it is discarded.
4761              * If it is taken, we recursively apply the subroutine to the point set before
4762              * and after the chosen point.
4763              * @param {Array} pts Array of {@link JXG.Coords}
4764              * @param {Number} i Index of an element of pts
4765              * @param {Number} j Index of an element of pts
4766              * @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>.
4767              * @param {Array} newPts Array of {@link JXG.Coords}
4768              * @ignore
4769              * @private
4770              */
4771             RDP = function (pts, i, j, eps, newPts) {
4772                 var result = findSplit(pts, i, j),
4773                     k = result[1];
4774 
4775                 if (isNaN(result[0])) {
4776                     RDP(pts, i, k - 1, eps, newPts);
4777                     newPts.push(pts[k]);
4778                     do {
4779                         ++k;
4780                     } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2]));
4781                     if (k <= j) {
4782                         newPts.push(pts[k]);
4783                     }
4784                     RDP(pts, k + 1, j, eps, newPts);
4785                 } else if (result[0] > eps) {
4786                     RDP(pts, i, k, eps, newPts);
4787                     RDP(pts, k, j, eps, newPts);
4788                 } else {
4789                     newPts.push(pts[j]);
4790                 }
4791             };
4792 
4793         len = pts.length;
4794 
4795         i = 0;
4796         while (endless) {
4797             // Search for the next point without NaN coordinates
4798             while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) {
4799                 i += 1;
4800             }
4801             // Search for the next position of a NaN point
4802             k = i + 1;
4803             while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) {
4804                 k += 1;
4805             }
4806             k--;
4807 
4808             // Only proceed if something is left
4809             if (i < len && k > i) {
4810                 newPts = [];
4811                 newPts[0] = pts[i];
4812                 RDP(pts, i, k, eps, newPts);
4813                 allPts = allPts.concat(newPts);
4814             }
4815             if (i >= len) {
4816                 break;
4817             }
4818             // Push the NaN point
4819             if (k < len - 1) {
4820                 allPts.push(pts[k + 1]);
4821             }
4822             i = k + 1;
4823         }
4824 
4825         return allPts;
4826     },
4827 
4828     /**
4829      * Old name for the implementation of the Ramer-Douglas-Peucker algorithm.
4830      * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker}
4831      * @memberof JXG.Math.Numerics
4832      */
4833     RamerDouglasPeuker: function (pts, eps) {
4834         JXG.deprecated("Numerics.RamerDouglasPeuker()", "Numerics.RamerDouglasPeucker()");
4835         return this.RamerDouglasPeucker(pts, eps);
4836     },
4837 
4838     /**
4839      * Implements the Visvalingam-Whyatt algorithm.
4840      * See M. Visvalingam, J. D. Whyatt:
4841      * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992
4842      *
4843      * The algorithm discards points which are not necessary from the polygonal line defined by the point array
4844      * pts (consisting of type JXG.Coords).
4845      * @param {Array} pts Array of {@link JXG.Coords}
4846      * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will
4847      *    be taken in any case.
4848      * @returns {Array} An array containing points which approximates the curve defined by pts.
4849      * @memberof JXG.Math.Numerics
4850      *
4851      * @example
4852      *     var i, p = [];
4853      *     for (i = 0; i < 5; ++i) {
4854      *         p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
4855      *     }
4856      *
4857      *     // Plot a cardinal spline curve
4858      *     var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
4859      *     var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
4860      *
4861      *     var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
4862      *     c.updateDataArray = function() {
4863      *         var i, len, points;
4864      *
4865      *         // Reduce number of intermediate points with Visvakingam-Whyatt to 6
4866      *         points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
4867      *         // Plot the remaining points
4868      *         len = points.length;
4869      *         this.dataX = [];
4870      *         this.dataY = [];
4871      *         for (i = 0; i < len; i++) {
4872      *             this.dataX.push(points[i].usrCoords[1]);
4873      *             this.dataY.push(points[i].usrCoords[2]);
4874      *         }
4875      *     };
4876      *     board.update();
4877      *
4878      * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div>
4879      * <script type="text/javascript">
4880      *     (function() {
4881      *         var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb',
4882      *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
4883      *
4884      *         var i, p = [];
4885      *         for (i = 0; i < 5; ++i) {
4886      *             p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
4887      *         }
4888      *
4889      *         // Plot a cardinal spline curve
4890      *         var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
4891      *         var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
4892      *
4893      *         var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
4894      *         c.updateDataArray = function() {
4895      *             var i, len, points;
4896      *
4897      *             // Reduce number of intermediate points with Visvakingam-Whyatt to 6
4898      *             points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
4899      *             // Plot the remaining points
4900      *             len = points.length;
4901      *             this.dataX = [];
4902      *             this.dataY = [];
4903      *             for (i = 0; i < len; i++) {
4904      *                 this.dataX.push(points[i].usrCoords[1]);
4905      *                 this.dataY.push(points[i].usrCoords[2]);
4906      *             }
4907      *         };
4908      *         board.update();
4909      *
4910      *     })();
4911      *
4912      * </script><pre>
4913      *
4914      */
4915     Visvalingam: function (pts, numPoints) {
4916         var i,
4917             len,
4918             vol,
4919             lastVol,
4920             linkedList = [],
4921             heap = [],
4922             points = [],
4923             lft,
4924             rt,
4925             lft2,
4926             rt2,
4927             obj;
4928 
4929         len = pts.length;
4930 
4931         // At least one intermediate point is needed
4932         if (len <= 2) {
4933             return pts;
4934         }
4935 
4936         // Fill the linked list
4937         // Add first point to the linked list
4938         linkedList[0] = {
4939             used: true,
4940             lft: null,
4941             node: null
4942         };
4943 
4944         // Add all intermediate points to the linked list,
4945         // whose triangle area is nonzero.
4946         lft = 0;
4947         for (i = 1; i < len - 1; i++) {
4948             vol = Math.abs(
4949                 JXG.Math.Numerics.det([
4950                     pts[i - 1].usrCoords,
4951                     pts[i].usrCoords,
4952                     pts[i + 1].usrCoords
4953                 ])
4954             );
4955             if (!isNaN(vol)) {
4956                 obj = {
4957                     v: vol,
4958                     idx: i
4959                 };
4960                 heap.push(obj);
4961                 linkedList[i] = {
4962                     used: true,
4963                     lft: lft,
4964                     node: obj
4965                 };
4966                 linkedList[lft].rt = i;
4967                 lft = i;
4968             }
4969         }
4970 
4971         // Add last point to the linked list
4972         linkedList[len - 1] = {
4973             used: true,
4974             rt: null,
4975             lft: lft,
4976             node: null
4977         };
4978         linkedList[lft].rt = len - 1;
4979 
4980         // Remove points until only numPoints intermediate points remain
4981         lastVol = -Infinity;
4982         while (heap.length > numPoints) {
4983             // Sort the heap with the updated volume values
4984             heap.sort(function (a, b) {
4985                 // descending sort
4986                 return b.v - a.v;
4987             });
4988 
4989             // Remove the point with the smallest triangle
4990             i = heap.pop().idx;
4991             linkedList[i].used = false;
4992             lastVol = linkedList[i].node.v;
4993 
4994             // Update the pointers of the linked list
4995             lft = linkedList[i].lft;
4996             rt = linkedList[i].rt;
4997             linkedList[lft].rt = rt;
4998             linkedList[rt].lft = lft;
4999 
5000             // Update the values for the volumes in the linked list
5001             lft2 = linkedList[lft].lft;
5002             if (lft2 !== null) {
5003                 vol = Math.abs(
5004                     JXG.Math.Numerics.det([
5005                         pts[lft2].usrCoords,
5006                         pts[lft].usrCoords,
5007                         pts[rt].usrCoords
5008                     ])
5009                 );
5010 
5011                 linkedList[lft].node.v = vol >= lastVol ? vol : lastVol;
5012             }
5013             rt2 = linkedList[rt].rt;
5014             if (rt2 !== null) {
5015                 vol = Math.abs(
5016                     JXG.Math.Numerics.det([
5017                         pts[lft].usrCoords,
5018                         pts[rt].usrCoords,
5019                         pts[rt2].usrCoords
5020                     ])
5021                 );
5022 
5023                 linkedList[rt].node.v = vol >= lastVol ? vol : lastVol;
5024             }
5025         }
5026 
5027         // Return an array with the remaining points
5028         i = 0;
5029         points = [pts[i]];
5030         do {
5031             i = linkedList[i].rt;
5032             points.push(pts[i]);
5033         } while (linkedList[i].rt !== null);
5034 
5035         return points;
5036     }
5037 };
5038 
5039 export default Mat.Numerics;
5040