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