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