1 /*
  2     Copyright 2008-2025
  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, // Fixed 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         // "Upper" horizontal line defined by function
3156         for (i = 0; i < n; i++) {
3157             if (type === 'simpson') {
3158                 sum += this._riemannValue(x, f, type, delta) * delta;
3159 
3160                 h = delta * 0.5;
3161                 f0 = f(x);
3162                 f1 = f(x + h);
3163                 f2 = f(x + 2 * h);
3164 
3165                 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5;
3166                 b = (f2 - f0) / (2 * h);
3167                 c = f1;
3168                 for (k = 0; k < steps; k++) {
3169                     xx = k * delta / steps - h;
3170                     xarr.push(x + xx + h);
3171                     yarr.push(a * xx * xx + b * xx + c);
3172                 }
3173                 x += delta;
3174                 y = f2;
3175             } else {
3176                 y = this._riemannValue(x, f, type, delta);
3177                 xarr.push(x);
3178                 yarr.push(y);
3179 
3180                 x += delta;
3181                 if (type === 'trapezoidal') {
3182                     f2 = f(x);
3183                     sum += (y + f2) * 0.5 * delta;
3184                     y = f2;
3185                 } else {
3186                     sum += y * delta;
3187                 }
3188 
3189                 xarr.push(x);
3190                 yarr.push(y);
3191             }
3192             xarr.push(x);
3193             yarr.push(y);
3194         }
3195 
3196         // "Lower" horizontal line
3197         // Go backwards
3198         for (i = 0; i < n; i++) {
3199             if (type === "simpson" && g) {
3200                 sum -= this._riemannValue(x, g, type, -delta) * delta;
3201 
3202                 h = delta * 0.5;
3203                 f0 = g(x);
3204                 f1 = g(x - h);
3205                 f2 = g(x - 2 * h);
3206 
3207                 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5;
3208                 b = (f2 - f0) / (2 * h);
3209                 c = f1;
3210                 for (k = 0; k < steps; k++) {
3211                     xx = k * delta / steps - h;
3212                     xarr.push(x - xx - h);
3213                     yarr.push(a * xx * xx + b * xx + c);
3214                 }
3215                 x -= delta;
3216                 y = f2;
3217             } else {
3218                 if (g) {
3219                     y = this._riemannValue(x, g, type, -delta);
3220                 } else {
3221                     y = 0.0;
3222                 }
3223                 xarr.push(x);
3224                 yarr.push(y);
3225 
3226                 x -= delta;
3227                 if (g) {
3228                     if (type === 'trapezoidal') {
3229                         f2 = g(x);
3230                         sum -= (y + f2) * 0.5 * delta;
3231                         y = f2;
3232                     } else {
3233                         sum -= y * delta;
3234                     }
3235                 }
3236             }
3237             xarr.push(x);
3238             yarr.push(y);
3239 
3240             // Draw the vertical lines
3241             xarr.push(x);
3242             yarr.push(f(x));
3243         }
3244 
3245         return [xarr, yarr, sum];
3246     },
3247 
3248     /**
3249      * Approximate the integral by Riemann sums.
3250      * Compute the area described by the riemann sum rectangles.
3251      *
3252      * If there is an element of type {@link Riemannsum}, then it is more efficient
3253      * to use the method JXG.Curve.Value() of this element instead.
3254      *
3255      * @param {Function_Array} f Function or array of two functions.
3256      * If f is a function the integral of this function is approximated by the Riemann sum.
3257      * If f is an array consisting of two functions the area between the two functions is approximated
3258      * by the Riemann sum.
3259      * @param {Number} n number of rectangles.
3260      * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'.
3261      *
3262      * @param {Number} start Left border of the approximation interval
3263      * @param {Number} end Right border of the approximation interval
3264      * @returns {Number} The sum of the areas of the rectangles.
3265      * @memberof JXG.Math.Numerics
3266      */
3267     riemannsum: function (f, n, type, start, end) {
3268         JXG.deprecated("Numerics.riemannsum()", "Numerics.riemann()[2]");
3269         return this.riemann(f, n, type, start, end)[2];
3270     },
3271 
3272     /**
3273      * Solve initial value problems numerically using <i>explicit</i> Runge-Kutta methods.
3274      * See {@link https://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
3275      * @param {object|String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
3276      * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
3277      * <pre>
3278      * {
3279      *     s: <Number>,
3280      *     A: <matrix>,
3281      *     b: <Array>,
3282      *     c: <Array>
3283      * }
3284      * </pre>
3285      * which corresponds to the Butcher tableau structure
3286      * shown here: https://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696 .
3287      * <i>Default</i> is 'euler'.
3288      * @param {Array} x0 Initial value vector. Even if the problem is one-dimensional, the initial value has to be given in an array.
3289      * @param {Array} I Interval on which to integrate.
3290      * @param {Number} N Number of integration intervals, i.e. there are <i>N+1</i> evaluation points.
3291      * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
3292      * 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
3293      * vector <tt>x</tt>, and has to return a vector of the same length as <tt>x</tt> has.
3294      * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
3295      * @example
3296      * // A very simple autonomous system dx(t)/dt = x(t);
3297      * var f = function(t, x) {
3298      *     return [x[0]];
3299      * }
3300      *
3301      * // Solve it with initial value x(0) = 1 on the interval [0, 2]
3302      * // with 20 evaluation points.
3303      * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
3304      *
3305      * // Prepare data for plotting the solution of the ode using a curve.
3306      * var dataX = [];
3307      * var dataY = [];
3308      * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
3309      * var i;
3310      * for(i=0; i<data.length; i++) {
3311      *     dataX[i] = i*h;
3312      *     dataY[i] = data[i][0];
3313      * }
3314      * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
3315      * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
3316      * <script type="text/javascript">
3317      * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
3318      * var f = function(t, x) {
3319      *     // we have to copy the value.
3320      *     // return x; would just return the reference.
3321      *     return [x[0]];
3322      * }
3323      * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
3324      * var dataX = [];
3325      * var dataY = [];
3326      * var h = 0.1;
3327      * for(var i=0; i<data.length; i++) {
3328      *     dataX[i] = i*h;
3329      *     dataY[i] = data[i][0];
3330      * }
3331      * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
3332      * </script><pre>
3333      * @memberof JXG.Math.Numerics
3334      */
3335     rungeKutta: function (butcher, x0, I, N, f) {
3336         var e,
3337             i, j, k, l, s,
3338             x = [],
3339             y = [],
3340             h = (I[1] - I[0]) / N,
3341             t = I[0],
3342             dim = x0.length,
3343             result = [],
3344             r = 0;
3345 
3346         if (Type.isString(butcher)) {
3347             butcher = predefinedButcher[butcher] || predefinedButcher.euler;
3348         }
3349         s = butcher.s;
3350 
3351         // Don't change x0, so copy it
3352         x = x0.slice();
3353 
3354         for (i = 0; i <= N; i++) {
3355             result[r] = x.slice();
3356 
3357             r++;
3358             k = [];
3359 
3360             for (j = 0; j < s; j++) {
3361                 // Init y = 0
3362                 for (e = 0; e < dim; e++) {
3363                     y[e] = 0.0;
3364                 }
3365 
3366                 // Calculate linear combination of former k's and save it in y
3367                 for (l = 0; l < j; l++) {
3368                     for (e = 0; e < dim; e++) {
3369                         y[e] += butcher.A[j][l] * h * k[l][e];
3370                     }
3371                 }
3372 
3373                 // Add x(t) to y
3374                 for (e = 0; e < dim; e++) {
3375                     y[e] += x[e];
3376                 }
3377 
3378                 // Calculate new k and add it to the k matrix
3379                 k.push(f(t + butcher.c[j] * h, y));
3380             }
3381 
3382             // Init y = 0
3383             for (e = 0; e < dim; e++) {
3384                 y[e] = 0.0;
3385             }
3386 
3387             for (l = 0; l < s; l++) {
3388                 for (e = 0; e < dim; e++) {
3389                     y[e] += butcher.b[l] * k[l][e];
3390                 }
3391             }
3392 
3393             for (e = 0; e < dim; e++) {
3394                 x[e] = x[e] + h * y[e];
3395             }
3396 
3397             t += h;
3398         }
3399 
3400         return result;
3401     },
3402 
3403     /**
3404      * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and
3405      * {@link JXG.Math.Numerics.chandrupatla}
3406      * @type Number
3407      * @default 80
3408      * @memberof JXG.Math.Numerics
3409      */
3410     maxIterationsRoot: 80,
3411 
3412     /**
3413      * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr}
3414      * @type Number
3415      * @default 500
3416      * @memberof JXG.Math.Numerics
3417      */
3418     maxIterationsMinimize: 500,
3419 
3420     /**
3421      * Given a number x_0, this function tries to find a second number x_1 such that
3422      * the function f has opposite signs at x_0 and x_1.
3423      * The return values have to be tested if the method succeeded.
3424      *
3425      * @param {Function} f Function, whose root is to be found
3426      * @param {Number} x0 Start value
3427      * @param {Object} [context] Parent object in case f is method of it
3428      * @returns {Array} [x_0, f(x_0), x_1, f(x_1)] in case that x_0 <= x_1
3429      *   or [x_1, f(x_1), x_0, f(x_0)] in case that x_1 < x_0.
3430      *
3431      * @see JXG.Math.Numerics.fzero
3432      * @see JXG.Math.Numerics.chandrupatla
3433      *
3434      * @memberof JXG.Math.Numerics
3435      */
3436     findBracket: function (f, x0, context) {
3437         var a, aa, fa, blist, b, fb, u, fu, i, len;
3438 
3439         if (Type.isArray(x0)) {
3440             return x0;
3441         }
3442 
3443         a = x0;
3444         fa = f.call(context, a);
3445         // nfev += 1;
3446 
3447         // Try to get b, by trying several values related to a
3448         aa = a === 0 ? 1 : a;
3449         blist = [
3450             a - 0.1 * aa,
3451             a + 0.1 * aa,
3452             a - 1,
3453             a + 1,
3454             a - 0.5 * aa,
3455             a + 0.5 * aa,
3456             a - 0.6 * aa,
3457             a + 0.6 * aa,
3458             a - 1 * aa,
3459             a + 1 * aa,
3460             a - 2 * aa,
3461             a + 2 * aa,
3462             a - 5 * aa,
3463             a + 5 * aa,
3464             a - 10 * aa,
3465             a + 10 * aa,
3466             a - 50 * aa,
3467             a + 50 * aa,
3468             a - 100 * aa,
3469             a + 100 * aa
3470         ];
3471         len = blist.length;
3472 
3473         for (i = 0; i < len; i++) {
3474             b = blist[i];
3475             fb = f.call(context, b);
3476             // nfev += 1;
3477 
3478             if (fa * fb <= 0) {
3479                 break;
3480             }
3481         }
3482         if (b < a) {
3483             u = a;
3484             a = b;
3485             b = u;
3486 
3487             fu = fa;
3488             fa = fb;
3489             fb = fu;
3490         }
3491         return [a, fa, b, fb];
3492     },
3493 
3494     /**
3495      *
3496      * Find zero of an univariate function f.
3497      * @param {function} f Function, whose root is to be found
3498      * @param {Array|Number} x0  Start value or start interval enclosing the root.
3499      * 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.
3500      * If x0 is a number, the algorithms tries to enclose the root by an interval [a, b] containing x0 and the root and
3501      * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method.
3502      * @param {Object} [context] Parent object in case f is method of it
3503      * @returns {Number} the approximation of the root
3504      * Algorithm:
3505      *  Brent's root finder from
3506      *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
3507      *  computations. M., Mir, 1980, p.180 of the Russian edition
3508      *  https://www.netlib.org/c/brent.shar
3509      *
3510      * If x0 is an array containing lower and upper bound for the zero
3511      * algorithm 748 is applied. Otherwise, if x0 is a number,
3512      * the algorithm tries to bracket a zero of f starting from x0.
3513      * If this fails, we fall back to Newton's method.
3514      *
3515      * @see JXG.Math.Numerics.chandrupatla
3516      * @see JXG.Math.Numerics.root
3517      * @see JXG.Math.Numerics.findBracket
3518      * @see JXG.Math.Numerics.Newton
3519      * @see JXG.Math.Numerics.fminbr
3520      * @memberof JXG.Math.Numerics
3521      */
3522     fzero: function (f, x0, context) {
3523         var a, b, c,
3524             fa, fb, fc,
3525             res, x00,
3526             prev_step,
3527             t1, t2,
3528             cb,
3529             tol_act,   // Actual tolerance
3530             p,         // Interpolation step is calculated in the form p/q; division
3531             q,         // operations is delayed until the last moment
3532             new_step,  // Step at this iteration
3533             eps = Mat.eps,
3534             maxiter = this.maxIterationsRoot,
3535             niter = 0;
3536         // nfev = 0;
3537 
3538         if (Type.isArray(x0)) {
3539             if (x0.length < 2) {
3540                 throw new Error(
3541                     "JXG.Math.Numerics.fzero: length of array x0 has to be at least two."
3542                 );
3543             }
3544 
3545             x00 = this.findDomain(f, x0, context);
3546             a = x00[0];
3547             b = x00[1];
3548             // a = x0[0];
3549             // b = x0[1];
3550 
3551             fa = f.call(context, a);
3552             // nfev += 1;
3553             fb = f.call(context, b);
3554             // nfev += 1;
3555         } else {
3556             res = this.findBracket(f, x0, context);
3557             a = res[0];
3558             fa = res[1];
3559             b = res[2];
3560             fb = res[3];
3561         }
3562 
3563         if (Math.abs(fa) <= eps) {
3564             return a;
3565         }
3566         if (Math.abs(fb) <= eps) {
3567             return b;
3568         }
3569 
3570         if (fa * fb > 0) {
3571             // Bracketing not successful, fall back to Newton's method or to fminbr
3572             if (Type.isArray(x0)) {
3573                 return this.fminbr(f, [a, b], context);
3574             }
3575 
3576             return this.Newton(f, a, context);
3577         }
3578 
3579         // OK, we have enclosed a zero of f.
3580         // Now we can start Brent's method
3581         c = a;
3582         fc = fa;
3583 
3584         // Main iteration loop
3585         while (niter < maxiter) {
3586             // Distance from the last but one to the last approximation
3587             prev_step = b - a;
3588 
3589             // Swap data for b to be the best approximation
3590             if (Math.abs(fc) < Math.abs(fb)) {
3591                 a = b;
3592                 b = c;
3593                 c = a;
3594 
3595                 fa = fb;
3596                 fb = fc;
3597                 fc = fa;
3598             }
3599             tol_act = 2 * eps * Math.abs(b) + eps * 0.5;
3600             new_step = (c - b) * 0.5;
3601 
3602             if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) {
3603                 //  Acceptable approx. is found
3604                 return b;
3605             }
3606 
3607             // Decide if the interpolation can be tried
3608             // If prev_step was large enough and was in true direction Interpolatiom may be tried
3609             if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
3610                 cb = c - b;
3611 
3612                 // If we have only two distinct points linear interpolation can only be applied
3613                 if (a === c) {
3614                     t1 = fb / fa;
3615                     p = cb * t1;
3616                     q = 1.0 - t1;
3617                     // Quadric inverse interpolation
3618                 } else {
3619                     q = fa / fc;
3620                     t1 = fb / fc;
3621                     t2 = fb / fa;
3622 
3623                     p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
3624                     q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
3625                 }
3626 
3627                 // p was calculated with the opposite sign; make p positive
3628                 if (p > 0) {
3629                     q = -q;
3630                     // and assign possible minus to q
3631                 } else {
3632                     p = -p;
3633                 }
3634 
3635                 // If b+p/q falls in [b,c] and isn't too large it is accepted
3636                 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
3637                 if (
3638                     p < 0.75 * cb * q - Math.abs(tol_act * q) * 0.5 &&
3639                     p < Math.abs(prev_step * q * 0.5)
3640                 ) {
3641                     new_step = p / q;
3642                 }
3643             }
3644 
3645             // Adjust the step to be not less than tolerance
3646             if (Math.abs(new_step) < tol_act) {
3647                 new_step = new_step > 0 ? tol_act : -tol_act;
3648             }
3649 
3650             // Save the previous approx.
3651             a = b;
3652             fa = fb;
3653             b += new_step;
3654             fb = f.call(context, b);
3655             // Do step to a new approxim.
3656             // nfev += 1;
3657 
3658             // Adjust c for it to have a sign opposite to that of b
3659             if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
3660                 c = a;
3661                 fc = fa;
3662             }
3663             niter++;
3664         } // End while
3665 
3666         return b;
3667     },
3668 
3669     /**
3670      * Find zero of an univariate function f.
3671      * @param {function} f Function, whose root is to be found
3672      * @param {Array|Number} x0  Start value or start interval enclosing the root.
3673      * 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.
3674      * If x0 is a number, the algorithms tries to enclose the root by an interval [a, b] containing x0 and the root and
3675      * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method.
3676      * @param {Object} [context] Parent object in case f is method of it
3677      * @returns {Number} the approximation of the root
3678      * Algorithm:
3679      * Chandrupatla's method, see
3680      * Tirupathi R. Chandrupatla,
3681      * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives",
3682      * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149.
3683      *
3684      * If x0 is an array containing lower and upper bound for the zero
3685      * algorithm 748 is applied. Otherwise, if x0 is a number,
3686      * the algorithm tries to bracket a zero of f starting from x0.
3687      * If this fails, we fall back to Newton's method.
3688      *
3689      * @see JXG.Math.Numerics.root
3690      * @see JXG.Math.Numerics.fzero
3691      * @see JXG.Math.Numerics.findBracket
3692      * @see JXG.Math.Numerics.Newton
3693      * @see JXG.Math.Numerics.fminbr
3694      * @memberof JXG.Math.Numerics
3695      */
3696     chandrupatla: function (f, x0, context) {
3697         var a, b, fa, fb,
3698             res,
3699             niter = 0,
3700             maxiter = this.maxIterationsRoot,
3701             rand = 1 + Math.random() * 0.001,
3702             t = 0.5 * rand,
3703             eps = Mat.eps, // 1.e-10,
3704             dlt = 0.00001,
3705             x1, x2, x3, x,
3706             f1, f2, f3, y,
3707             xm, fm,
3708             tol, tl,
3709             xi, ph, fl, fh,
3710             AL, A, B, C, D;
3711 
3712         if (Type.isArray(x0)) {
3713             if (x0.length < 2) {
3714                 throw new Error(
3715                     "JXG.Math.Numerics.fzero: length of array x0 has to be at least two."
3716                 );
3717             }
3718 
3719             a = x0[0];
3720             fa = f.call(context, a);
3721             // nfev += 1;
3722             b = x0[1];
3723             fb = f.call(context, b);
3724             // nfev += 1;
3725         } else {
3726             res = this.findBracket(f, x0, context);
3727             a = res[0];
3728             fa = res[1];
3729             b = res[2];
3730             fb = res[3];
3731         }
3732 
3733         if (fa * fb > 0) {
3734             // Bracketing not successful, fall back to Newton's method or to fminbr
3735             if (Type.isArray(x0)) {
3736                 return this.fminbr(f, [a, b], context);
3737             }
3738 
3739             return this.Newton(f, a, context);
3740         }
3741 
3742         x1 = a;
3743         x2 = b;
3744         f1 = fa;
3745         f2 = fb;
3746         do {
3747             x = x1 + t * (x2 - x1);
3748             y = f.call(context, x);
3749 
3750             // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point
3751             if (Math.sign(y) === Math.sign(f1)) {
3752                 x3 = x1;
3753                 x1 = x;
3754                 f3 = f1;
3755                 f1 = y;
3756             } else {
3757                 x3 = x2;
3758                 x2 = x1;
3759                 f3 = f2;
3760                 f2 = f1;
3761             }
3762             x1 = x;
3763             f1 = y;
3764 
3765             xm = x1;
3766             fm = f1;
3767             if (Math.abs(f2) < Math.abs(f1)) {
3768                 xm = x2;
3769                 fm = f2;
3770             }
3771             tol = 2 * eps * Math.abs(xm) + 0.5 * dlt;
3772             tl = tol / Math.abs(x2 - x1);
3773             if (tl > 0.5 || fm === 0) {
3774                 break;
3775             }
3776             // If inverse quadratic interpolation holds, use it
3777             xi = (x1 - x2) / (x3 - x2);
3778             ph = (f1 - f2) / (f3 - f2);
3779             fl = 1 - Math.sqrt(1 - xi);
3780             fh = Math.sqrt(xi);
3781             if (fl < ph && ph < fh) {
3782                 AL = (x3 - x1) / (x2 - x1);
3783                 A = f1 / (f2 - f1);
3784                 B = f3 / (f2 - f3);
3785                 C = f1 / (f3 - f1);
3786                 D = f2 / (f3 - f2);
3787                 t = A * B + C * D * AL;
3788             } else {
3789                 t = 0.5 * rand;
3790             }
3791             // Adjust t away from the interval boundary
3792             if (t < tl) {
3793                 t = tl;
3794             }
3795             if (t > 1 - tl) {
3796                 t = 1 - tl;
3797             }
3798             niter++;
3799         } while (niter <= maxiter);
3800         // console.log(niter);
3801 
3802         return xm;
3803     },
3804 
3805     /**
3806      * Find a small enclosing interval of the domain of a function by
3807      * tightening the input interval x0.
3808      * <p>
3809      * This is a helper function which is used in {@link JXG.Math.Numerics.fminbr},
3810      * {@link JXG.Math.Numerics.fzero}, and  {@link JXG.Curve.getLabelPosition}
3811      * to avoid search in an interval where the function is mostly undefined.
3812      *
3813      * @param {function} f
3814      * @param {Array} x0 Start interval
3815      * @param {Object} context Parent object in case f is method of it
3816      * @param {Boolean} [outer=true] if true take a proper enclosing array. If false return the domain such that the function is defined
3817      * at its  borders.
3818      * @returns Array
3819      *
3820      * @example
3821      * var f = (x) => Math.sqrt(x);
3822      * console.log(JXG.Math.Numerics.findDomain(f, [-5, 5]));
3823      *
3824      * // Output: [ -0.00020428174445492973, 5 ]
3825      *
3826      * @example
3827      * var f = (x) => Math.sqrt(x);
3828      * console.log(JXG.Math.Numerics.findDomain(f, [-5, 5], null, false));
3829      *
3830      * // Output: [ 0.00020428174562965915, 5 ]
3831      */
3832     findDomain: function (f, x0, context, outer) {
3833         var a, b, c, fc,
3834             x,
3835             gr = 1 - 1 / 1.61803398875,
3836             eps = 0.001,
3837             cnt,
3838             max_cnt = 20;
3839 
3840         if (outer === undefined) {
3841             outer = true;
3842         }
3843 
3844         if (!Type.isArray(x0) || x0.length < 2) {
3845             throw new Error(
3846                 "JXG.Math.Numerics.findDomain: length of array x0 has to be at least two."
3847             );
3848         }
3849 
3850         x = x0.slice();
3851         a = x[0];
3852         b = x[1];
3853         fc = f.call(context, a);
3854         if (isNaN(fc)) {
3855             // Divide the interval with the golden ratio
3856             // and keep a such that f(a) = NaN
3857             cnt = 0;
3858             while (b - a > eps && cnt < max_cnt) {
3859                 c = (b - a) * gr + a;
3860                 fc = f.call(context, c);
3861                 if (isNaN(fc)) {
3862                     a = c;
3863                 } else {
3864                     b = c;
3865                 }
3866                 cnt++;
3867             }
3868             if (outer) {
3869                 x[0] = a;
3870             } else {
3871                 x[0] = b;
3872             }
3873             // x[0] = a;
3874         }
3875 
3876         a = x[0];
3877         b = x[1];
3878         fc = f.call(context, b);
3879         if (isNaN(fc)) {
3880             // Divide the interval with the golden ratio
3881             // and keep b such that f(b) = NaN
3882             cnt = 0;
3883             while (b - a > eps && cnt < max_cnt) {
3884                 c = b - (b - a) * gr;
3885                 fc = f.call(context, c);
3886                 if (isNaN(fc)) {
3887                     b = c;
3888                 } else {
3889                     a = c;
3890                 }
3891                 cnt++;
3892             }
3893             if (outer) {
3894                 x[1] = b;
3895             } else {
3896                 x[1] = a;
3897             }
3898             // x[1] = b;
3899         }
3900 
3901         return x;
3902     },
3903 
3904     /**
3905      *
3906      * Find minimum of an univariate function f.
3907      * <p>
3908      * Algorithm:
3909      *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
3910      *  computations. M., Mir, 1980, p.180 of the Russian edition
3911      *
3912      * @param {function} f Function, whose minimum is to be found
3913      * @param {Array} x0  Start interval enclosing the minimum
3914      * @param {Object} [context] Parent object in case f is method of it
3915      * @returns {Number} the approximation of the minimum value position
3916      * @memberof JXG.Math.Numerics
3917      **/
3918     fminbr: function (f, x0, context) {
3919         var a, b, x, v, w,
3920             fx, fv, fw,
3921             x00,
3922             range, middle_range, tol_act, new_step,
3923             p, q, t, ft,
3924             r = (3.0 - Math.sqrt(5.0)) * 0.5,      // Golden section ratio
3925             tol = Mat.eps,
3926             sqrteps = Mat.eps, // Math.sqrt(Mat.eps),
3927             maxiter = this.maxIterationsMinimize,
3928             niter = 0;
3929         // nfev = 0;
3930 
3931         if (!Type.isArray(x0) || x0.length < 2) {
3932             throw new Error(
3933                 "JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."
3934             );
3935         }
3936 
3937         x00 = this.findDomain(f, x0, context);
3938         a = x00[0];
3939         b = x00[1];
3940         v = a + r * (b - a);
3941         fv = f.call(context, v);
3942 
3943         // First step - always gold section
3944         // nfev += 1;
3945         x = v;
3946         w = v;
3947         fx = fv;
3948         fw = fv;
3949 
3950         while (niter < maxiter) {
3951             // Range over the interval in which we are looking for the minimum
3952             range = b - a;
3953             middle_range = (a + b) * 0.5;
3954 
3955             // Actual tolerance
3956             tol_act = sqrteps * Math.abs(x) + tol / 3.0;
3957 
3958             if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) {
3959                 // Acceptable approx. is found
3960                 return x;
3961             }
3962 
3963             // Obtain the golden section step
3964             new_step = r * (x < middle_range ? b - x : a - x);
3965 
3966             // Decide if the interpolation can be tried. If x and w are distinct, interpolatiom may be tried
3967             if (Math.abs(x - w) >= tol_act) {
3968                 // Interpolation step is calculated as p/q;
3969                 // division operation is delayed until last moment
3970                 t = (x - w) * (fx - fv);
3971                 q = (x - v) * (fx - fw);
3972                 p = (x - v) * q - (x - w) * t;
3973                 q = 2 * (q - t);
3974 
3975                 if (q > 0) {
3976                     p = -p; // q was calculated with the opposite sign; make q positive
3977                 } else {
3978                     q = -q; // // and assign possible minus to p
3979                 }
3980                 if (
3981                     Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b]
3982                     p > q * (a - x + 2 * tol_act) &&        //  not too close to a and
3983                     p < q * (b - x - 2 * tol_act)
3984                 ) {
3985                     // b, and isn't too large
3986                     new_step = p / q; // it is accepted
3987                 }
3988                 // If p / q is too large then the
3989                 // golden section procedure can
3990                 // reduce [a,b] range to more
3991                 // extent
3992             }
3993 
3994             // Adjust the step to be not less than tolerance
3995             if (Math.abs(new_step) < tol_act) {
3996                 if (new_step > 0) {
3997                     new_step = tol_act;
3998                 } else {
3999                     new_step = -tol_act;
4000                 }
4001             }
4002 
4003             // Obtain the next approximation to min
4004             // and reduce the enveloping range
4005 
4006             // Tentative point for the min
4007             t = x + new_step;
4008             ft = f.call(context, t);
4009             // nfev += 1;
4010 
4011             // t is a better approximation
4012             if (ft <= fx) {
4013                 // Reduce the range so that t would fall within it
4014                 if (t < x) {
4015                     b = x;
4016                 } else {
4017                     a = x;
4018                 }
4019 
4020                 // Assign the best approx to x
4021                 v = w;
4022                 w = x;
4023                 x = t;
4024 
4025                 fv = fw;
4026                 fw = fx;
4027                 fx = ft;
4028                 // x remains the better approx
4029             } else {
4030                 // Reduce the range enclosing x
4031                 if (t < x) {
4032                     a = t;
4033                 } else {
4034                     b = t;
4035                 }
4036 
4037                 if (ft <= fw || w === x) {
4038                     v = w;
4039                     w = t;
4040                     fv = fw;
4041                     fw = ft;
4042                 } else if (ft <= fv || v === x || v === w) {
4043                     v = t;
4044                     fv = ft;
4045                 }
4046             }
4047             niter += 1;
4048         }
4049 
4050         return x;
4051     },
4052 
4053     /**
4054      * GLOMIN seeks a global minimum of a function F(X) in an interval [A,B]
4055      * and is the adaption of the algorithm GLOMIN by Richard Brent.
4056      *
4057      * Here is the original documentation:
4058      * <pre>
4059      *
4060      * Discussion:
4061      *
4062      * This function assumes that F(X) is twice continuously differentiable over [A,B]
4063      * and that |F''(X)| <= M for all X in [A,B].
4064      *
4065      * Licensing:
4066      *   This code is distributed under the GNU LGPL license.
4067      *
4068      * Modified:
4069      *
4070      *   17 April 2008
4071      *
4072      * Author:
4073      *
4074      *   Original FORTRAN77 version by Richard Brent.
4075      *   C version by John Burkardt.
4076      *   https://people.math.sc.edu/Burkardt/c_src/brent/brent.c
4077      *
4078      * Reference:
4079      *
4080      *   Richard Brent,
4081      *  Algorithms for Minimization Without Derivatives,
4082      *   Dover, 2002,
4083      *  ISBN: 0-486-41998-3,
4084      *   LC: QA402.5.B74.
4085      *
4086      * Parameters:
4087      *
4088      *   Input, double A, B, the endpoints of the interval.
4089      *  It must be the case that A < B.
4090      *
4091      *   Input, double C, an initial guess for the global
4092      *  minimizer.  If no good guess is known, C = A or B is acceptable.
4093      *
4094      *  Input, double M, the bound on the second derivative.
4095      *
4096      *   Input, double MACHEP, an estimate for the relative machine
4097      *  precision.
4098      *
4099      *   Input, double E, a positive tolerance, a bound for the
4100      *  absolute error in the evaluation of F(X) for any X in [A,B].
4101      *
4102      *   Input, double T, a positive error tolerance.
4103      *
4104      *    Input, double F (double x ), a user-supplied
4105      *  function whose global minimum is being sought.
4106      *
4107      *   Output, double *X, the estimated value of the abscissa
4108      *  for which F attains its global minimum value in [A,B].
4109      *
4110      *   Output, double GLOMIN, the value F(X).
4111      * </pre>
4112      *
4113      * In JSXGraph, some parameters of the original algorithm are set to fixed values:
4114      * <ul>
4115      *  <li> M = 10000000.0
4116      *  <li> C = A or B, depending if f(A) <= f(B)
4117      *  <li> T = JXG.Math.eps
4118      *  <li> E = JXG.Math.eps * JXG.Math.eps
4119      *  <li> MACHEP = JXG.Math.eps * JXG.Math.eps * JXG.Math.eps
4120      * </ul>
4121      * @param {function} f Function, whose global minimum is to be found
4122      * @param {Array} x0 Array of length 2 determining the interval [A, B] for which the global minimum is to be found
4123      * @returns {Array} [x, y] x is the position of the global minimum and y = f(x).
4124      */
4125     glomin: function (f, x0) {
4126         var a0, a2, a3, d0, d1, d2, h,
4127             k, m2,
4128             p, q, qs,
4129             r, s, sc,
4130             y, y0, y1, y2, y3, yb,
4131             z0, z1, z2,
4132             a, b, c, x,
4133             m = 10000000.0,
4134             t = Mat.eps, // * Mat.eps,
4135             e = Mat.eps * Mat.eps,
4136             machep = Mat.eps * Mat.eps * Mat.eps;
4137 
4138         a = x0[0];
4139         b = x0[1];
4140         c = (f(a) < f(b)) ? a : b;
4141 
4142         a0 = b;
4143         x = a0;
4144         a2 = a;
4145         y0 = f(b);
4146         yb = y0;
4147         y2 = f(a);
4148         y = y2;
4149 
4150         if (y0 < y) {
4151             y = y0;
4152         } else {
4153             x = a;
4154         }
4155 
4156         if (m <= 0.0 || b <= a) {
4157             return y;
4158         }
4159 
4160         m2 = 0.5 * (1.0 + 16.0 * machep) * m;
4161 
4162         if (c <= a || b <= c) {
4163             sc = 0.5 * (a + b);
4164         } else {
4165             sc = c;
4166         }
4167 
4168         y1 = f(sc);
4169         k = 3;
4170         d0 = a2 - sc;
4171         h = 9.0 / 11.0;
4172 
4173         if (y1 < y) {
4174             x = sc;
4175             y = y1;
4176         }
4177 
4178         for (; ;) {
4179             d1 = a2 - a0;
4180             d2 = sc - a0;
4181             z2 = b - a2;
4182             z0 = y2 - y1;
4183             z1 = y2 - y0;
4184             r = d1 * d1 * z0 - d0 * d0 * z1;
4185             p = r;
4186             qs = 2.0 * (d0 * z1 - d1 * z0);
4187             q = qs;
4188 
4189             if (k < 1000000 || y2 <= y) {
4190                 for (; ;) {
4191                     if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) <
4192                         z2 * m2 * r * (z2 * q - r)) {
4193 
4194                         a3 = a2 + r / q;
4195                         y3 = f(a3);
4196 
4197                         if (y3 < y) {
4198                             x = a3;
4199                             y = y3;
4200                         }
4201                     }
4202                     k = ((1611 * k) % 1048576);
4203                     q = 1.0;
4204                     r = (b - a) * 0.00001 * k;
4205 
4206                     if (z2 <= r) {
4207                         break;
4208                     }
4209                 }
4210             } else {
4211                 k = ((1611 * k) % 1048576);
4212                 q = 1.0;
4213                 r = (b - a) * 0.00001 * k;
4214 
4215                 while (r < z2) {
4216                     if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) <
4217                         z2 * m2 * r * (z2 * q - r)) {
4218 
4219                         a3 = a2 + r / q;
4220                         y3 = f(a3);
4221 
4222                         if (y3 < y) {
4223                             x = a3;
4224                             y = y3;
4225                         }
4226                     }
4227                     k = ((1611 * k) % 1048576);
4228                     q = 1.0;
4229                     r = (b - a) * 0.00001 * k;
4230                 }
4231             }
4232 
4233             r = m2 * d0 * d1 * d2;
4234             s = Math.sqrt(((y2 - y) + t) / m2);
4235             h = 0.5 * (1.0 + h);
4236             p = h * (p + 2.0 * r * s);
4237             q = q + 0.5 * qs;
4238             r = - 0.5 * (d0 + (z0 + 2.01 * e) / (d0 * m2));
4239 
4240             if (r < s || d0 < 0.0) {
4241                 r = a2 + s;
4242             } else {
4243                 r = a2 + r;
4244             }
4245 
4246             if (0.0 < p * q) {
4247                 a3 = a2 + p / q;
4248             } else {
4249                 a3 = r;
4250             }
4251 
4252             for (; ;) {
4253                 a3 = Math.max(a3, r);
4254 
4255                 if (b <= a3) {
4256                     a3 = b;
4257                     y3 = yb;
4258                 } else {
4259                     y3 = f(a3);
4260                 }
4261 
4262                 if (y3 < y) {
4263                     x = a3;
4264                     y = y3;
4265                 }
4266 
4267                 d0 = a3 - a2;
4268 
4269                 if (a3 <= r) {
4270                     break;
4271                 }
4272 
4273                 p = 2.0 * (y2 - y3) / (m * d0);
4274 
4275                 if ((1.0 + 9.0 * machep) * d0 <= Math.abs(p)) {
4276                     break;
4277                 }
4278 
4279                 if (0.5 * m2 * (d0 * d0 + p * p) <= (y2 - y) + (y3 - y) + 2.0 * t) {
4280                     break;
4281                 }
4282                 a3 = 0.5 * (a2 + a3);
4283                 h = 0.9 * h;
4284             }
4285 
4286             if (b <= a3) {
4287                 break;
4288             }
4289 
4290             a0 = sc;
4291             sc = a2;
4292             a2 = a3;
4293             y0 = y1;
4294             y1 = y2;
4295             y2 = y3;
4296         }
4297 
4298         return [x, y];
4299     },
4300 
4301     /**
4302      * Determine all roots of a polynomial with real or complex coefficients by using the
4303      * iterative method attributed to Weierstrass, Durand, Kerner, Aberth, and Ehrlich. In particular,
4304      * the iteration method with cubic convergence is used that is usually attributed to Ehrlich-Aberth.
4305      * <p>
4306      * The returned roots are sorted with respect to their real values.
4307      * <p> This method makes use of the JSXGraph classes {@link JXG.Complex} and {@link JXG.C} to handle
4308      * complex numbers.
4309      *
4310      * @param {Array} a Array of coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4311      * The coefficients are of type Number or JXG.Complex.
4312      * @param {Number} [deg] Optional degree of the polynomial. Otherwise all entries are taken, with
4313      * leading zeros removed.
4314      * @param {Number} [tol=Number.EPSILON] Approximation tolerance
4315      * @param {Number} [max_it=30] Maximum number of iterations
4316      * @param {Array} [initial_values=null] Array of initial values for the roots. If not given,
4317      * starting values are determined by the method of Ozawa.
4318      * @returns {Array} Array of complex numbers (of JXG.Complex) approximating the roots of the polynomial.
4319      * @memberof JXG.Math.Numerics
4320      * @see JXG.Complex
4321      * @see JXG.C
4322      *
4323      * @example
4324      * // Polynomial p(z) = -1 + 1z^2
4325      * var i, roots,
4326      *     p = [-1, 0, 1];
4327      *
4328      * roots = JXG.Math.Numerics.polzeros(p);
4329      * for (i = 0; i < roots.length; i++) {
4330      *     console.log(i, roots[i].toString());
4331      * }
4332      * // Output:
4333      *   0 -1 + -3.308722450212111e-24i
4334      *   1 1 + 0i
4335      *
4336      * @example
4337      * // Polynomial p(z) = -1 + 3z - 9z^2 + z^3 - 8z^6 + 9z^7 - 9z^8 + z^9
4338      * var i, roots,
4339      *     p = [-1, 3, -9, 1, 0, 0, -8, 9, -9, 1];
4340      *
4341      * roots = JXG.Math.Numerics.polzeros(p);
4342      * for (i = 0; i < roots.length; i++) {
4343      *     console.log(i, roots[i].toString());
4344      * }
4345      * // Output:
4346      * 0 -0.7424155888401961 + 0.4950476539211721i
4347      * 1 -0.7424155888401961 + -0.4950476539211721i
4348      * 2 0.16674869833354108 + 0.2980502714610669i
4349      * 3 0.16674869833354108 + -0.29805027146106694i
4350      * 4 0.21429002063640837 + 1.0682775088132996i
4351      * 5 0.21429002063640842 + -1.0682775088132999i
4352      * 6 0.861375497926218 + -0.6259177003583295i
4353      * 7 0.8613754979262181 + 0.6259177003583295i
4354      * 8 8.000002743888055 + -1.8367099231598242e-40i
4355      *
4356      */
4357     polzeros: function (coeffs, deg, tol, max_it, initial_values) {
4358         var i, le, off, it,
4359             debug = false,
4360             cc = [],
4361             obvious = [],
4362             roots = [],
4363 
4364             /**
4365              * Horner method to evaluate polynomial or the derivative thereof for complex numbers,
4366              * i.e. coefficients and variable are complex.
4367              * @function
4368              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4369              * @param {JXG.Complex} z Value for which the polynomial will be evaluated.
4370              * @param {Boolean} [derivative=false] If true the derivative will be evaluated.
4371              * @ignore
4372              */
4373             hornerComplex = function (a, z, derivative) {
4374                 var i, s,
4375                     n = a.length - 1;
4376 
4377                 derivative = derivative || false;
4378                 if (derivative) {
4379                     // s = n * a_n
4380                     s = JXG.C.mult(n, a[n]);
4381                     for (i = n - 1; i > 0; i--) {
4382                         // s = s * z + i * a_i
4383                         s.mult(z);
4384                         s.add(JXG.C.mult(a[i], i));
4385                     }
4386                 } else {
4387                     // s = a_n
4388                     s = JXG.C.copy(a[n]);
4389                     for (i = n - 1; i >= 0; i--) {
4390                         // s = s * z + a_i
4391                         s.mult(z);
4392                         s.add(a[i]);
4393                     }
4394                 }
4395                 return s;
4396             },
4397 
4398             /**
4399              * Horner method to evaluate reciprocal polynomial or the derivative thereof for complex numbers,
4400              * i.e. coefficients and variable are complex.
4401              * @function
4402              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4403              * @param {JXG.Complex} z Value for which the reciprocal polynomial will be evaluated.
4404              * @param {Boolean} [derivative=false] If true the derivative will be evaluated.
4405              * @ignore
4406              */
4407             hornerRec = function (a, x, derivative) {
4408                 var i, s,
4409                     n = a.length - 1;
4410 
4411                 derivative = derivative || false;
4412                 if (derivative) {
4413                     // s = n * a_0
4414                     s = JXG.C.mult(n, a[0]);
4415                     for (i = n - 1; i > 0; i--) {
4416                         // s = s * x + i * a_{n-i}
4417                         s.mult(x);
4418                         s.add(JXG.C.mult(a[n - i], i));
4419                     }
4420                 } else {
4421                     // s = a_0
4422                     s = JXG.C.copy(a[0]);
4423                     for (i = n - 1; i >= 0; i--) {
4424                         // s = s * x + a_{n-i}
4425                         s.mult(x);
4426                         s.add(a[n - i]);
4427                     }
4428                 }
4429                 return s;
4430             },
4431 
4432             /**
4433              * Horner method to evaluate real polynomial at a real value.
4434              * @function
4435              * @param {Array} a Array of real coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4436              * @param {Number} z Value for which the polynomial will be evaluated.
4437              * @ignore
4438              */
4439             horner = function (a, x) {
4440                 var i, s,
4441                     n = a.length - 1;
4442 
4443                 s = a[n];
4444                 for (i = n - 1; i >= 0; i--) {
4445                     s = s * x + a[i];
4446                 }
4447                 return s;
4448             },
4449 
4450             /**
4451              * Determine start values for the Aberth iteration, see
4452              * Ozawa, "An experimental study of the starting values
4453              * of the Durand-Kerner-Aberth iteration" (1995).
4454              *
4455              * @function
4456              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4457              * @returns {Array} Array Initial values for the roots.
4458              * @ignore
4459              */
4460             initial_guess = function (a) {
4461                 var i, r,
4462                     n = a.length - 1, // degree
4463                     alpha1 = Math.PI * 2 / n,
4464                     alpha0 = Math.PI / n * 0.5,
4465                     b, z,
4466                     init = [];
4467 
4468 
4469                 // From Ozawa, "An experimental study of the starting values
4470                 // of the Durand-Kerner-Aberth iteration" (1995)
4471 
4472                 // b is the arithmetic mean of the roots.
4473                 // With is Vieta's formula <https://en.wikipedia.org/wiki/Vieta%27s_formulas>
4474                 //   b = -a_{n-1} / (n * a_n)
4475                 b = JXG.C.mult(-1, a[n - 1]);
4476                 b.div(JXG.C.mult(n, a[n]));
4477 
4478                 // r is the geometric mean of the deviations |b - root_i|.
4479                 // Using
4480                 //   p(z) = a_n prod(z - root_i)
4481                 // and therefore
4482                 //   |p(b)| = |a_n| prod(|b - root_i|)
4483                 // we arrive at:
4484                 //   r = |p(b)/a_n|^(1/n)
4485                 z = JXG.C.div(hornerComplex(a, b), a[n]);
4486                 r = Math.pow(JXG.C.abs(z), 1 / n);
4487                 if (r === 0) { r = 1; }
4488 
4489                 for (i = 0; i < n; i++) {
4490                     a = new JXG.Complex(r * Math.cos(alpha1 * i + alpha0), r * Math.sin(alpha1 * i + alpha0));
4491                     init[i] = JXG.C.add(b, a);
4492                 }
4493 
4494                 return init;
4495             },
4496 
4497             /**
4498              * Ehrlich-Aberth iteration. The stopping criterion is from
4499              * D.A. Bini, "Numerical computation of polynomial zeros
4500              * by means of Aberths's method", Numerical Algorithms (1996).
4501              *
4502              * @function
4503              * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2...
4504              * @param {Number} mu Machine precision
4505              * @param {Number} max_it Maximum number of iterations
4506              * @param {Array} z Initial guess for the roots. Will be changed in place.
4507              * @returns {Number} Number of iterations
4508              * @ignore
4509              */
4510             aberthIteration = function (cc, mu, max_it, z) {
4511                 var k, i, j,
4512                     done = [],
4513                     cr = [],
4514                     gamma, x,
4515                     done_sum = 0,
4516                     num, denom, s, pp,
4517                     n = z.length;
4518 
4519                 for (i = 0; i < n; i++) {
4520                     done.push(false);
4521                 }
4522                 for (i = 0; i < cc.length; i++) {
4523                     cr.push(JXG.C.abs(cc[i]) * (4 * i + 1));
4524                 }
4525                 for (k = 0; k < max_it && done_sum < n; k++) {
4526                     for (i = 0; i < n; i++) {
4527                         if (done[i]) {
4528                             continue;
4529                         }
4530                         num = hornerComplex(cc, z[i]);
4531                         x = JXG.C.abs(z[i]);
4532 
4533                         // Stopping criterion by D.A. Bini
4534                         // "Numerical computation of polynomial zeros
4535                         // by means of Aberths's method", Numerical Algorithms (1996).
4536                         //
4537                         if (JXG.C.abs(num) < mu * horner(cr, x)) {
4538                             done[i] = true;
4539                             done_sum++;
4540                             if (done_sum === n) {
4541                                 break;
4542                             }
4543                             continue;
4544                         }
4545 
4546                         // num = P(z_i) / P'(z_i)
4547                         if (x > 1) {
4548                             gamma = JXG.C.div(1, z[i]);
4549                             pp = hornerRec(cc, gamma, true);
4550                             pp.div(hornerRec(cc, gamma));
4551                             pp.mult(gamma);
4552                             num = JXG.C.sub(n, pp);
4553                             num = JXG.C.div(z[i], num);
4554                         } else {
4555                             num.div(hornerComplex(cc, z[i], true));
4556                         }
4557 
4558                         // denom = sum_{i\neq j} 1 / (z_i  - z_j)
4559                         denom = new JXG.Complex(0);
4560                         for (j = 0; j < n; j++) {
4561                             if (j === i) {
4562                                 continue;
4563                             }
4564                             s = JXG.C.sub(z[i], z[j]);
4565                             s = JXG.C.div(1, s);
4566                             denom.add(s);
4567                         }
4568 
4569                         // num = num / 1 - num * sum_{i\neq j} 1 / (z_i - z_j)
4570                         denom.mult(num);
4571                         denom = JXG.C.sub(1, denom);
4572                         num.div(denom);
4573                         // z_i = z_i - num
4574                         z[i].sub(num);
4575                     }
4576                 }
4577 
4578                 return k;
4579             };
4580 
4581 
4582         tol = tol || Number.EPSILON;
4583         max_it = max_it || 30;
4584 
4585         le = coeffs.length;
4586         if (JXG.isNumber(deg) && deg >= 0 && deg < le - 1) {
4587             le = deg + 1;
4588         }
4589 
4590         // Convert coefficient array to complex numbers
4591         for (i = 0; i < le; i++) {
4592             cc.push(new JXG.Complex(coeffs[i]));
4593         }
4594 
4595         // Search for (multiple) roots at x=0
4596         for (i = 0; i < le; i++) {
4597             if (cc[i].real !== 0 || cc[i].imaginary !== 0) {
4598                 off = i;
4599                 break;
4600             }
4601         }
4602 
4603         // Deflate root x=0, store roots at x=0 in obvious
4604         for (i = 0; i < off; i++) {
4605             obvious.push(new JXG.Complex(0));
4606         }
4607         cc = cc.slice(off);
4608         le = cc.length;
4609 
4610         // Remove leading zeros from the coefficient array
4611         for (i = le - 1; i >= 0; i--) {
4612             if (cc[i].real !== 0 || cc[i].imaginary !== 0) {
4613                 break;
4614             }
4615             cc.pop();
4616         }
4617         le = cc.length;
4618         if (le === 0) {
4619             return [];
4620         }
4621 
4622         // From now on we can assume that the
4623         // constant coefficient and the leading coefficient
4624         // are not zero.
4625         if (initial_values) {
4626             for (i = 0; i < le - 1; i++) {
4627                 roots.push(new JXG.Complex(initial_values[i]));
4628             }
4629         } else {
4630             roots = initial_guess(cc);
4631         }
4632         it = aberthIteration(cc, tol, max_it, roots);
4633 
4634         // Append the roots at x=0
4635         roots = obvious.concat(roots);
4636 
4637         if (debug) {
4638             console.log("Iterations:", it);
4639             console.log('Roots:');
4640             for (i = 0; i < roots.length; i++) {
4641                 console.log(i, roots[i].toString(), JXG.C.abs(hornerComplex(cc, roots[i])));
4642             }
4643         }
4644 
4645         // Sort roots according to their real part
4646         roots.sort(function (a, b) {
4647             if (a.real < b.real) {
4648                 return -1;
4649             }
4650             if (a.real > b.real) {
4651                 return 1;
4652             }
4653             return 0;
4654         });
4655 
4656         return roots;
4657     },
4658 
4659     /**
4660      * Implements the Ramer-Douglas-Peucker algorithm.
4661      * It discards points which are not necessary from the polygonal line defined by the point array
4662      * pts. The computation is done in screen coordinates.
4663      * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points.
4664      * @param {Array} pts Array of {@link JXG.Coords}
4665      * @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>.
4666      * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points.
4667      * @memberof JXG.Math.Numerics
4668      */
4669     RamerDouglasPeucker: function (pts, eps) {
4670         var allPts = [],
4671             newPts = [],
4672             i, k, len,
4673             endless = true,
4674 
4675             /**
4676              * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
4677              * It searches for the point between index i and j which
4678              * has the largest distance from the line between the points i and j.
4679              * @param {Array} pts Array of {@link JXG.Coords}
4680              * @param {Number} i Index of a point in pts
4681              * @param {Number} j Index of a point in pts
4682              * @ignore
4683              * @private
4684              */
4685             findSplit = function (pts, i, j) {
4686                 var d, k, ci, cj, ck,
4687                     x0, y0, x1, y1,
4688                     den, lbda,
4689                     eps = Mat.eps * Mat.eps,
4690                     huge = 10000,
4691                     dist = 0,
4692                     f = i;
4693 
4694                 if (j - i < 2) {
4695                     return [-1.0, 0];
4696                 }
4697 
4698                 ci = pts[i].scrCoords;
4699                 cj = pts[j].scrCoords;
4700 
4701                 if (isNaN(ci[1]) || isNaN(ci[2])) {
4702                     return [NaN, i];
4703                 }
4704                 if (isNaN(cj[1]) || isNaN(cj[2])) {
4705                     return [NaN, j];
4706                 }
4707 
4708                 for (k = i + 1; k < j; k++) {
4709                     ck = pts[k].scrCoords;
4710                     if (isNaN(ck[1]) || isNaN(ck[2])) {
4711                         return [NaN, k];
4712                     }
4713 
4714                     x0 = ck[1] - ci[1];
4715                     y0 = ck[2] - ci[2];
4716                     x1 = cj[1] - ci[1];
4717                     y1 = cj[2] - ci[2];
4718                     x0 = x0 === Infinity ? huge : x0;
4719                     y0 = y0 === Infinity ? huge : y0;
4720                     x1 = x1 === Infinity ? huge : x1;
4721                     y1 = y1 === Infinity ? huge : y1;
4722                     x0 = x0 === -Infinity ? -huge : x0;
4723                     y0 = y0 === -Infinity ? -huge : y0;
4724                     x1 = x1 === -Infinity ? -huge : x1;
4725                     y1 = y1 === -Infinity ? -huge : y1;
4726                     den = x1 * x1 + y1 * y1;
4727 
4728                     if (den > eps) {
4729                         lbda = (x0 * x1 + y0 * y1) / den;
4730 
4731                         if (lbda < 0.0) {
4732                             lbda = 0.0;
4733                         } else if (lbda > 1.0) {
4734                             lbda = 1.0;
4735                         }
4736 
4737                         x0 = x0 - lbda * x1;
4738                         y0 = y0 - lbda * y1;
4739                         d = x0 * x0 + y0 * y0;
4740                     } else {
4741                         lbda = 0.0;
4742                         d = x0 * x0 + y0 * y0;
4743                     }
4744 
4745                     if (d > dist) {
4746                         dist = d;
4747                         f = k;
4748                     }
4749                 }
4750                 return [Math.sqrt(dist), f];
4751             },
4752             /**
4753              * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}.
4754              * It runs recursively through the point set and searches the
4755              * point which has the largest distance from the line between the first point and
4756              * the last point. If the distance from the line is greater than eps, this point is
4757              * included in our new point set otherwise it is discarded.
4758              * If it is taken, we recursively apply the subroutine to the point set before
4759              * and after the chosen point.
4760              * @param {Array} pts Array of {@link JXG.Coords}
4761              * @param {Number} i Index of an element of pts
4762              * @param {Number} j Index of an element of pts
4763              * @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>.
4764              * @param {Array} newPts Array of {@link JXG.Coords}
4765              * @ignore
4766              * @private
4767              */
4768             RDP = function (pts, i, j, eps, newPts) {
4769                 var result = findSplit(pts, i, j),
4770                     k = result[1];
4771 
4772                 if (isNaN(result[0])) {
4773                     RDP(pts, i, k - 1, eps, newPts);
4774                     newPts.push(pts[k]);
4775                     do {
4776                         ++k;
4777                     } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2]));
4778                     if (k <= j) {
4779                         newPts.push(pts[k]);
4780                     }
4781                     RDP(pts, k + 1, j, eps, newPts);
4782                 } else if (result[0] > eps) {
4783                     RDP(pts, i, k, eps, newPts);
4784                     RDP(pts, k, j, eps, newPts);
4785                 } else {
4786                     newPts.push(pts[j]);
4787                 }
4788             };
4789 
4790         len = pts.length;
4791 
4792         i = 0;
4793         while (endless) {
4794             // Search for the next point without NaN coordinates
4795             while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) {
4796                 i += 1;
4797             }
4798             // Search for the next position of a NaN point
4799             k = i + 1;
4800             while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) {
4801                 k += 1;
4802             }
4803             k--;
4804 
4805             // Only proceed if something is left
4806             if (i < len && k > i) {
4807                 newPts = [];
4808                 newPts[0] = pts[i];
4809                 RDP(pts, i, k, eps, newPts);
4810                 allPts = allPts.concat(newPts);
4811             }
4812             if (i >= len) {
4813                 break;
4814             }
4815             // Push the NaN point
4816             if (k < len - 1) {
4817                 allPts.push(pts[k + 1]);
4818             }
4819             i = k + 1;
4820         }
4821 
4822         return allPts;
4823     },
4824 
4825     /**
4826      * Old name for the implementation of the Ramer-Douglas-Peucker algorithm.
4827      * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker}
4828      * @memberof JXG.Math.Numerics
4829      */
4830     RamerDouglasPeuker: function (pts, eps) {
4831         JXG.deprecated("Numerics.RamerDouglasPeuker()", "Numerics.RamerDouglasPeucker()");
4832         return this.RamerDouglasPeucker(pts, eps);
4833     },
4834 
4835     /**
4836      * Implements the Visvalingam-Whyatt algorithm.
4837      * See M. Visvalingam, J. D. Whyatt:
4838      * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992
4839      *
4840      * The algorithm discards points which are not necessary from the polygonal line defined by the point array
4841      * pts (consisting of type JXG.Coords).
4842      * @param {Array} pts Array of {@link JXG.Coords}
4843      * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will
4844      *    be taken in any case.
4845      * @returns {Array} An array containing points which approximates the curve defined by pts.
4846      * @memberof JXG.Math.Numerics
4847      *
4848      * @example
4849      *     var i, p = [];
4850      *     for (i = 0; i < 5; ++i) {
4851      *         p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
4852      *     }
4853      *
4854      *     // Plot a cardinal spline curve
4855      *     var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
4856      *     var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
4857      *
4858      *     var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
4859      *     c.updateDataArray = function() {
4860      *         var i, len, points;
4861      *
4862      *         // Reduce number of intermediate points with Visvakingam-Whyatt to 6
4863      *         points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
4864      *         // Plot the remaining points
4865      *         len = points.length;
4866      *         this.dataX = [];
4867      *         this.dataY = [];
4868      *         for (i = 0; i < len; i++) {
4869      *             this.dataX.push(points[i].usrCoords[1]);
4870      *             this.dataY.push(points[i].usrCoords[2]);
4871      *         }
4872      *     };
4873      *     board.update();
4874      *
4875      * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div>
4876      * <script type="text/javascript">
4877      *     (function() {
4878      *         var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb',
4879      *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
4880      *
4881      *         var i, p = [];
4882      *         for (i = 0; i < 5; ++i) {
4883      *             p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6]));
4884      *         }
4885      *
4886      *         // Plot a cardinal spline curve
4887      *         var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5);
4888      *         var cu1 = board.create('curve', splineArr, {strokeColor: 'green'});
4889      *
4890      *         var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'});
4891      *         c.updateDataArray = function() {
4892      *             var i, len, points;
4893      *
4894      *             // Reduce number of intermediate points with Visvakingam-Whyatt to 6
4895      *             points = JXG.Math.Numerics.Visvalingam(cu1.points, 6);
4896      *             // Plot the remaining points
4897      *             len = points.length;
4898      *             this.dataX = [];
4899      *             this.dataY = [];
4900      *             for (i = 0; i < len; i++) {
4901      *                 this.dataX.push(points[i].usrCoords[1]);
4902      *                 this.dataY.push(points[i].usrCoords[2]);
4903      *             }
4904      *         };
4905      *         board.update();
4906      *
4907      *     })();
4908      *
4909      * </script><pre>
4910      *
4911      */
4912     Visvalingam: function (pts, numPoints) {
4913         var i,
4914             len,
4915             vol,
4916             lastVol,
4917             linkedList = [],
4918             heap = [],
4919             points = [],
4920             lft,
4921             rt,
4922             lft2,
4923             rt2,
4924             obj;
4925 
4926         len = pts.length;
4927 
4928         // At least one intermediate point is needed
4929         if (len <= 2) {
4930             return pts;
4931         }
4932 
4933         // Fill the linked list
4934         // Add first point to the linked list
4935         linkedList[0] = {
4936             used: true,
4937             lft: null,
4938             node: null
4939         };
4940 
4941         // Add all intermediate points to the linked list,
4942         // whose triangle area is nonzero.
4943         lft = 0;
4944         for (i = 1; i < len - 1; i++) {
4945             vol = Math.abs(
4946                 JXG.Math.Numerics.det([
4947                     pts[i - 1].usrCoords,
4948                     pts[i].usrCoords,
4949                     pts[i + 1].usrCoords
4950                 ])
4951             );
4952             if (!isNaN(vol)) {
4953                 obj = {
4954                     v: vol,
4955                     idx: i
4956                 };
4957                 heap.push(obj);
4958                 linkedList[i] = {
4959                     used: true,
4960                     lft: lft,
4961                     node: obj
4962                 };
4963                 linkedList[lft].rt = i;
4964                 lft = i;
4965             }
4966         }
4967 
4968         // Add last point to the linked list
4969         linkedList[len - 1] = {
4970             used: true,
4971             rt: null,
4972             lft: lft,
4973             node: null
4974         };
4975         linkedList[lft].rt = len - 1;
4976 
4977         // Remove points until only numPoints intermediate points remain
4978         lastVol = -Infinity;
4979         while (heap.length > numPoints) {
4980             // Sort the heap with the updated volume values
4981             heap.sort(function (a, b) {
4982                 // descending sort
4983                 return b.v - a.v;
4984             });
4985 
4986             // Remove the point with the smallest triangle
4987             i = heap.pop().idx;
4988             linkedList[i].used = false;
4989             lastVol = linkedList[i].node.v;
4990 
4991             // Update the pointers of the linked list
4992             lft = linkedList[i].lft;
4993             rt = linkedList[i].rt;
4994             linkedList[lft].rt = rt;
4995             linkedList[rt].lft = lft;
4996 
4997             // Update the values for the volumes in the linked list
4998             lft2 = linkedList[lft].lft;
4999             if (lft2 !== null) {
5000                 vol = Math.abs(
5001                     JXG.Math.Numerics.det([
5002                         pts[lft2].usrCoords,
5003                         pts[lft].usrCoords,
5004                         pts[rt].usrCoords
5005                     ])
5006                 );
5007 
5008                 linkedList[lft].node.v = vol >= lastVol ? vol : lastVol;
5009             }
5010             rt2 = linkedList[rt].rt;
5011             if (rt2 !== null) {
5012                 vol = Math.abs(
5013                     JXG.Math.Numerics.det([
5014                         pts[lft].usrCoords,
5015                         pts[rt].usrCoords,
5016                         pts[rt2].usrCoords
5017                     ])
5018                 );
5019 
5020                 linkedList[rt].node.v = vol >= lastVol ? vol : lastVol;
5021             }
5022         }
5023 
5024         // Return an array with the remaining points
5025         i = 0;
5026         points = [pts[i]];
5027         do {
5028             i = linkedList[i].rt;
5029             points.push(pts[i]);
5030         } while (linkedList[i].rt !== null);
5031 
5032         return points;
5033     }
5034 };
5035 
5036 export default Mat.Numerics;
5037