1 /*
  2     Copyright 2008-2015
  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 
 33 /*global JXG: true, define: true*/
 34 /*jslint nomen: true, plusplus: true*/
 35 
 36 /* depends:
 37  utils/type
 38  math/math
 39  */
 40 
 41 /**
 42  * @fileoverview In this file the namespace Math.Numerics is defined, which holds numerical
 43  * algorithms for solving linear equations etc.
 44  */
 45 
 46 define(['utils/type', 'math/math'], function (Type, Mat) {
 47 
 48     "use strict";
 49 
 50     // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order).
 51     var predefinedButcher = {
 52         rk4: {
 53             s: 4,
 54             A: [
 55                 [ 0,  0,  0, 0],
 56                 [0.5, 0,  0, 0],
 57                 [ 0, 0.5, 0, 0],
 58                 [ 0,  0,  1, 0]
 59             ],
 60             b: [1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0],
 61             c: [0, 0.5, 0.5, 1]
 62         },
 63         heun: {
 64             s: 2,
 65             A: [
 66                 [0, 0],
 67                 [1, 0]
 68             ],
 69             b: [0.5, 0.5],
 70             c: [0, 1]
 71         },
 72         euler: {
 73             s: 1,
 74             A: [
 75                 [0]
 76             ],
 77             b: [1],
 78             c: [0]
 79         }
 80     };
 81 
 82     /**
 83      * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables.
 84      * @namespace
 85      * @name JXG.Math.Numerics
 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, intergration_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 && typeof config.number_of_nodes === 'number' ? 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          * Integral of function f over interval.
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          * @returns {Number} The value of the integral of f over interval
486          * @see JXG.Math.Numerics.NewtonCotes
487          * @memberof JXG.Math.Numerics
488          */
489         I: function (interval, f) {
490             return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'});
491         },
492 
493         /**
494          * Newton's method to find roots of a funtion in one variable.
495          * @param {function} f We search for a solution of f(x)=0.
496          * @param {Number} x initial guess for the root, i.e. start value.
497          * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
498          * the function is a method of an object and contains a reference to its parent object via "this".
499          * @returns {Number} A root of the function f.
500          * @memberof JXG.Math.Numerics
501          */
502         Newton: function (f, x, context) {
503             var df,
504                 i = 0,
505                 h = Mat.eps,
506                 newf = f.apply(context, [x]),
507                 nfev = 1;
508 
509             // For compatibility
510             if (Type.isArray(x)) {
511                 x = x[0];
512             }
513 
514             while (i < 50 && Math.abs(newf) > h) {
515                 df = this.D(f, context)(x);
516                 nfev += 2;
517 
518                 if (Math.abs(df) > h) {
519                     x -= newf / df;
520                 } else {
521                     x += (Math.random() * 0.2 - 1.0);
522                 }
523 
524                 newf = f.apply(context, [x]);
525                 nfev += 1;
526                 i += 1;
527             }
528 
529             return x;
530         },
531 
532         /**
533          * Abstract method to find roots of univariate functions.
534          * @param {function} f We search for a solution of f(x)=0.
535          * @param {Number} x initial guess for the root, i.e. starting value.
536          * @param {Object} context optional object that is treated as "this" in the function body. This is useful if
537          * the function is a method of an object and contains a reference to its parent object via "this".
538          * @returns {Number} A root of the function f.
539          * @memberof JXG.Math.Numerics
540          */
541         root: function (f, x, context) {
542             return this.fzero(f, x, context);
543         },
544 
545         /**
546          * Compute an intersection of the curves c1 and c2
547          * with a generalized Newton method.
548          * We want to find values t1, t2 such that
549          * c1(t1) = c2(t2), i.e.
550          * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0).
551          * We set
552          * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2))
553          *
554          * The Jacobian J is defined by
555          * J = (a, b)
556          *     (c, d)
557          * where
558          * a = c1_x'(t1)
559          * b = -c2_x'(t2)
560          * c = c1_y'(t1)
561          * d = -c2_y'(t2)
562          *
563          * The inverse J^(-1) of J is equal to
564          *  (d, -b)/
565          *  (-c, a) / (ad-bc)
566          *
567          * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f).
568          * If the function meetCurveCurve possesses the properties
569          * t1memo and t2memo then these are taken as start values
570          * for the Newton algorithm.
571          * After stopping of the Newton algorithm the values of t1 and t2 are stored in
572          * t1memo and t2memo.
573          *
574          * @param {JXG.Curve} c1 Curve, Line or Circle
575          * @param {JXG.Curve} c2 Curve, Line or Circle
576          * @param {Number} t1ini start value for t1
577          * @param {Number} t2ini start value for t2
578          * @returns {JXG.Coords} intersection point
579          * @memberof JXG.Math.Numerics
580          */
581         generalizedNewton: function (c1, c2, t1ini, t2ini) {
582             var t1, t2,
583                 a, b, c, d, disc,
584                 e, f, F,
585                 D00, D01,
586                 D10, D11,
587                 count = 0;
588 
589             if (this.generalizedNewton.t1memo) {
590                 t1 = this.generalizedNewton.t1memo;
591                 t2 = this.generalizedNewton.t2memo;
592             } else {
593                 t1 = t1ini;
594                 t2 = t2ini;
595             }
596 
597             e = c1.X(t1) - c2.X(t2);
598             f = c1.Y(t1) - c2.Y(t2);
599             F = e * e + f * f;
600 
601             D00 = this.D(c1.X, c1);
602             D01 = this.D(c2.X, c2);
603             D10 = this.D(c1.Y, c1);
604             D11 = this.D(c2.Y, c2);
605 
606             while (F > Mat.eps && count < 10) {
607                 a = D00(t1);
608                 b = -D01(t2);
609                 c = D10(t1);
610                 d = -D11(t2);
611                 disc = a * d - b * c;
612                 t1 -= (d * e - b * f) / disc;
613                 t2 -= (a * f - c * e) / disc;
614                 e = c1.X(t1) - c2.X(t2);
615                 f = c1.Y(t1) - c2.Y(t2);
616                 F = e * e + f * f;
617                 count += 1;
618             }
619 
620             this.generalizedNewton.t1memo = t1;
621             this.generalizedNewton.t2memo = t2;
622 
623             if (Math.abs(t1) < Math.abs(t2)) {
624                 return [c1.X(t1), c1.Y(t1)];
625             }
626 
627             return [c2.X(t2), c2.Y(t2)];
628         },
629 
630         /**
631          * Returns the Lagrange polynomials for curves with equidistant nodes, see
632          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
633          * SIAM Review, Vol 46, No 3, (2004) 501-517.
634          * The graph of the parametric curve [x(t),y(t)] runs through the given points.
635          * @param {Array} p Array of JXG.Points
636          * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve
637          * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero.
638          * @memberof JXG.Math.Numerics
639          */
640         Neville: function (p) {
641             var w = [],
642                 /** @ignore */
643                 makeFct = function (fun) {
644                     return function (t, suspendedUpdate) {
645                         var i, d, s,
646                             bin = Mat.binomial,
647                             len = p.length,
648                             len1 = len - 1,
649                             num = 0.0,
650                             denom = 0.0;
651 
652                         if (!suspendedUpdate) {
653                             s = 1;
654                             for (i = 0; i < len; i++) {
655                                 w[i] = bin(len1, i) * s;
656                                 s *= (-1);
657                             }
658                         }
659 
660                         d = t;
661 
662                         for (i = 0; i < len; i++) {
663                             if (d === 0) {
664                                 return p[i][fun]();
665                             }
666                             s = w[i] / d;
667                             d -= 1;
668                             num += p[i][fun]() * s;
669                             denom += s;
670                         }
671                         return num / denom;
672                     };
673                 },
674 
675                 xfct = makeFct('X'),
676                 yfct = makeFct('Y');
677 
678             return [xfct, yfct, 0, function () {
679                 return p.length - 1;
680             }];
681         },
682 
683         /**
684          * Calculates second derivatives at the given knots.
685          * @param {Array} x x values of knots
686          * @param {Array} y y values of knots
687          * @returns {Array} Second derivatives of the interpolated function at the knots.
688          * @see #splineEval
689          * @memberof JXG.Math.Numerics
690          */
691         splineDef: function (x, y) {
692             var pair, i, l,
693                 n = Math.min(x.length, y.length),
694                 diag = [],
695                 z = [],
696                 data = [],
697                 dx = [],
698                 delta = [],
699                 F = [];
700 
701             if (n === 2) {
702                 return [0, 0];
703             }
704 
705             for (i = 0; i < n; i++) {
706                 pair = {X: x[i], Y: y[i]};
707                 data.push(pair);
708             }
709             data.sort(function (a, b) {
710                 return a.X - b.X;
711             });
712             for (i = 0; i < n; i++) {
713                 x[i] = data[i].X;
714                 y[i] = data[i].Y;
715             }
716 
717             for (i = 0; i < n - 1; i++) {
718                 dx.push(x[i + 1] - x[i]);
719             }
720             for (i = 0; i < n - 2; i++) {
721                 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i]));
722             }
723 
724             // ForwardSolve
725             diag.push(2 * (dx[0] + dx[1]));
726             z.push(delta[0]);
727 
728             for (i = 0; i < n - 3; i++) {
729                 l = dx[i + 1] / diag[i];
730                 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]);
731                 z.push(delta[i + 1] - l * z[i]);
732             }
733 
734             // BackwardSolve
735             F[n - 3] = z[n - 3] / diag[n - 3];
736             for (i = n - 4; i >= 0; i--) {
737                 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i];
738             }
739 
740             // Generate f''-Vector
741             for (i = n - 3; i >= 0; i--) {
742                 F[i + 1] = F[i];
743             }
744 
745             // natural cubic spline
746             F[0] = 0;
747             F[n - 1] = 0;
748 
749             return F;
750         },
751 
752         /**
753          * Evaluate points on spline.
754          * @param {Number,Array} x0 A single float value or an array of values to evaluate
755          * @param {Array} x x values of knots
756          * @param {Array} y y values of knots
757          * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef}
758          * @see #splineDef
759          * @returns {Number,Array} A single value or an array, depending on what is given as x0.
760          * @memberof JXG.Math.Numerics
761          */
762         splineEval: function (x0, x, y, F) {
763             var i, j, a, b, c, d, x_,
764                 n = Math.min(x.length, y.length),
765                 l = 1,
766                 asArray = false,
767                 y0 = [];
768 
769             // number of points to be evaluated
770             if (Type.isArray(x0)) {
771                 l = x0.length;
772                 asArray = true;
773             } else {
774                 x0 = [x0];
775             }
776 
777             for (i = 0; i < l; i++) {
778                 // is x0 in defining interval?
779                 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) {
780                     return NaN;
781                 }
782 
783                 // determine part of spline in which x0 lies
784                 for (j = 1; j < n; j++) {
785                     if (x0[i] <= x[j]) {
786                         break;
787                     }
788                 }
789 
790                 j -= 1;
791 
792                 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1];
793                 // determine the coefficients of the polynomial in this interval
794                 a = y[j];
795                 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]);
796                 c = F[j] / 2;
797                 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j]));
798                 // evaluate x0[i]
799                 x_ = x0[i] - x[j];
800                 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_);
801                 y0.push(a + (b + (c + d * x_) * x_) * x_);
802             }
803 
804             if (asArray) {
805                 return y0;
806             }
807 
808             return y0[0];
809         },
810 
811         /**
812          * Generate a string containing the function term of a polynomial.
813          * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i.
814          * @param {Number} deg Degree of the polynomial
815          * @param {String} varname Name of the variable (usually 'x')
816          * @param {Number} prec Precision
817          * @returns {String} A string containg the function term of the polynomial.
818          * @memberof JXG.Math.Numerics
819          */
820         generatePolynomialTerm: function (coeffs, deg, varname, prec) {
821             var i, t = [];
822 
823             for (i = deg; i >= 0; i--) {
824                 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']);
825 
826                 if (i > 1) {
827                     t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']);
828                 } else if (i === 1) {
829                     t = t.concat(['*', varname, ' + ']);
830                 }
831             }
832 
833             return t.join('');
834         },
835 
836         /**
837          * Computes the polynomial through a given set of coordinates in Lagrange form.
838          * Returns the Lagrange polynomials, see
839          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
840          * SIAM Review, Vol 46, No 3, (2004) 501-517.
841          * @param {Array} p Array of JXG.Points
842          * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points.
843          * @memberof JXG.Math.Numerics
844          */
845         lagrangePolynomial: function (p) {
846             var w = [],
847                 /** @ignore */
848                 fct = function (x, suspendedUpdate) {
849                     var i, j, k, xi, s, M,
850                         len = p.length,
851                         num = 0,
852                         denom = 0;
853 
854                     if (!suspendedUpdate) {
855                         for (i = 0; i < len; i++) {
856                             w[i] = 1.0;
857                             xi = p[i].X();
858 
859                             for (k = 0; k < len; k++) {
860                                 if (k !== i) {
861                                     w[i] *= (xi - p[k].X());
862                                 }
863                             }
864 
865                             w[i] = 1 / w[i];
866                         }
867 
868                         M = [];
869 
870                         for (j = 0; j < len; j++) {
871                             M.push([1]);
872                         }
873                     }
874 
875                     for (i = 0; i < len; i++) {
876                         xi = p[i].X();
877 
878                         if (x === xi) {
879                             return p[i].Y();
880                         }
881 
882                         s = w[i] / (x - xi);
883                         denom += s;
884                         num += s * p[i].Y();
885                     }
886 
887                     return num / denom;
888                 };
889 
890             fct.getTerm = function () {
891                 return '';
892             };
893 
894             return fct;
895         },
896 
897         /**
898          * Computes the cubic cardinal spline curve through a given set of points. The curve
899          * is uniformly parametrized.
900          * Two artificial control points at the beginning and the end are added.
901          * @param {Array} points Array consisting of JXG.Points.
902          * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1.
903          * tau=1/2 give Catmull-Rom splines.
904          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
905          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
906          * returning the length of the points array minus three.
907          * @memberof JXG.Math.Numerics
908         */
909         CardinalSpline: function (points, tau) {
910             var p,
911                 coeffs = [],
912                 // control point at the beginning and at the end
913                 first = {},
914                 last = {},
915                 makeFct,
916                 _tau;
917 
918             if (Type.isFunction(tau)) {
919                 _tau = tau;
920             } else {
921                 _tau = function () { return tau; };
922             }
923 
924             /** @ignore */
925             makeFct = function (which) {
926                 return function (t, suspendedUpdate) {
927                     var s, c,
928                         len = points.length,
929                         tau = _tau();
930 
931                     if (len < 2) {
932                         return NaN;
933                     }
934 
935                     if (!suspendedUpdate) {
936                         first[which] = function () {
937                             return 2 * points[0][which]() - points[1][which]();
938                         };
939 
940                         last[which] = function () {
941                             return 2 * points[len - 1][which]() - points[len - 2][which]();
942                         };
943 
944                         p = [first].concat(points, [last]);
945                         coeffs[which] = [];
946 
947                         for (s = 0; s < len - 1; s++) {
948                             coeffs[which][s] = [
949                                 1 / tau * p[s + 1][which](),
950                                 -p[s][which]() +   p[s + 2][which](),
951                                 2 * p[s][which]() + (-3 / tau + 1) * p[s + 1][which]() + (3 / tau - 2) * p[s + 2][which]() - p[s + 3][which](),
952                                 -p[s][which]() + (2 / tau - 1) * p[s + 1][which]() + (-2 / tau + 1) * p[s + 2][which]() + p[s + 3][which]()
953                             ];
954                         }
955                     }
956 
957                     len += 2;  // add the two control points
958 
959                     if (isNaN(t)) {
960                         return NaN;
961                     }
962 
963                     // This is necessary for our advanced plotting algorithm:
964                     if (t <= 0.0) {
965                         return p[1][which]();
966                     }
967 
968                     if (t >= len - 3) {
969                         return p[len - 2][which]();
970                     }
971 
972                     s = Math.floor(t);
973 
974                     if (s === t) {
975                         return p[s][which]();
976                     }
977 
978                     t -= s;
979                     c = coeffs[which][s];
980 
981                     return tau * (((c[3] * t + c[2]) * t + c[1]) * t + c[0]);
982                 };
983             };
984 
985             return [makeFct('X'), makeFct('Y'), 0,
986                 function () {
987                     return points.length - 1;
988                 }];
989         },
990 
991         /**
992          * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve
993          * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5.
994          * Two artificial control points at the beginning and the end are added.
995          * @param {Array} points Array consisting of JXG.Points.
996          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
997          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
998          * returning the length of the points array minus three.
999          * @memberof JXG.Math.Numerics
1000         */
1001         CatmullRomSpline: function (points) {
1002             return this.CardinalSpline(points, 0.5);
1003         },
1004 
1005         /**
1006          * Computes the regression polynomial of a given degree through a given set of coordinates.
1007          * Returns the regression polynomial function.
1008          * @param {Number,function,Slider} degree number, function or slider.
1009          * Either
1010          * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in
1011          * an array of {@link JXG.Point}s or {@link JXG.Coords}. In the latter case, the <tt>dataY</tt> parameter will be ignored.
1012          * @param {Array} dataY Array containing the y-coordinates of the data set,
1013          * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree.
1014          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
1015          * @memberof JXG.Math.Numerics
1016          */
1017         regressionPolynomial: function (degree, dataX, dataY) {
1018             var coeffs, deg, dX, dY, inputType, fct,
1019                 term = '';
1020 
1021             // Slider
1022             if (Type.isPoint(degree) && typeof degree.Value === 'function') {
1023                 /** @ignore */
1024                 deg = function () {
1025                     return degree.Value();
1026                 };
1027             // function
1028             } else if (Type.isFunction(degree)) {
1029                 deg = degree;
1030             // number
1031             } else if (Type.isNumber(degree)) {
1032                 /** @ignore */
1033                 deg = function () {
1034                     return degree;
1035                 };
1036             } else {
1037                 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'.");
1038             }
1039 
1040             // Parameters degree, dataX, dataY
1041             if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) {
1042                 inputType = 0;
1043             // Parameters degree, point array
1044             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) {
1045                 inputType = 1;
1046             } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) {
1047                 inputType = 2;
1048             } else {
1049                 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters.");
1050             }
1051 
1052             /** @ignore */
1053             fct = function (x, suspendedUpdate) {
1054                 var i, j, M, MT, y, B, c, s, d,
1055                     // input data
1056                     len = dataX.length;
1057 
1058                 d = Math.floor(deg());
1059 
1060                 if (!suspendedUpdate) {
1061                     // point list as input
1062                     if (inputType === 1) {
1063                         dX = [];
1064                         dY = [];
1065 
1066                         for (i = 0; i < len; i++) {
1067                             dX[i] = dataX[i].X();
1068                             dY[i] = dataX[i].Y();
1069                         }
1070                     }
1071 
1072                     if (inputType === 2) {
1073                         dX = [];
1074                         dY = [];
1075 
1076                         for (i = 0; i < len; i++) {
1077                             dX[i] = dataX[i].usrCoords[1];
1078                             dY[i] = dataX[i].usrCoords[2];
1079                         }
1080                     }
1081 
1082                     // check for functions
1083                     if (inputType === 0) {
1084                         dX = [];
1085                         dY = [];
1086 
1087                         for (i = 0; i < len; i++) {
1088                             if (Type.isFunction(dataX[i])) {
1089                                 dX.push(dataX[i]());
1090                             } else {
1091                                 dX.push(dataX[i]);
1092                             }
1093 
1094                             if (Type.isFunction(dataY[i])) {
1095                                 dY.push(dataY[i]());
1096                             } else {
1097                                 dY.push(dataY[i]);
1098                             }
1099                         }
1100                     }
1101 
1102                     M = [];
1103 
1104                     for (j = 0; j < len; j++) {
1105                         M.push([1]);
1106                     }
1107 
1108                     for (i = 1; i <= d; i++) {
1109                         for (j = 0; j < len; j++) {
1110                             M[j][i] = M[j][i - 1] * dX[j];
1111                         }
1112                     }
1113 
1114                     y = dY;
1115                     MT = Mat.transpose(M);
1116                     B = Mat.matMatMult(MT, M);
1117                     c = Mat.matVecMult(MT, y);
1118                     coeffs = Mat.Numerics.Gauss(B, c);
1119                     term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3);
1120                 }
1121 
1122                 // Horner's scheme to evaluate polynomial
1123                 s = coeffs[d];
1124 
1125                 for (i = d - 1; i >= 0; i--) {
1126                     s = (s * x + coeffs[i]);
1127                 }
1128 
1129                 return s;
1130             };
1131 
1132             fct.getTerm = function () {
1133                 return term;
1134             };
1135 
1136             return fct;
1137         },
1138 
1139         /**
1140          * Computes the cubic Bezier curve through a given set of points.
1141          * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}.
1142          * The points at position k with k mod 3 = 0 are the data points,
1143          * points at position k with k mod 3 = 1 or 2 are the control points.
1144          * @returns {Array} An array consisting of two functions of one parameter t which return the
1145          * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting
1146          * no parameters and returning one third of the length of the points.
1147          * @memberof JXG.Math.Numerics
1148          */
1149         bezier: function (points) {
1150             var len, flen,
1151                 /** @ignore */
1152                 makeFct = function (which) {
1153                     return function (t, suspendedUpdate) {
1154                         var z = Math.floor(t) * 3,
1155                             t0 = t % 1,
1156                             t1 = 1 - t0;
1157 
1158                         if (!suspendedUpdate) {
1159                             flen = 3 * Math.floor((points.length - 1) / 3);
1160                             len = Math.floor(flen / 3);
1161                         }
1162 
1163                         if (t < 0) {
1164                             return points[0][which]();
1165                         }
1166 
1167                         if (t >= len) {
1168                             return points[flen][which]();
1169                         }
1170 
1171                         if (isNaN(t)) {
1172                             return NaN;
1173                         }
1174 
1175                         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;
1176                     };
1177                 };
1178 
1179             return [makeFct('X'), makeFct('Y'), 0,
1180                 function () {
1181                     return Math.floor(points.length / 3);
1182                 }];
1183         },
1184 
1185         /**
1186          * Computes the B-spline curve of order k (order = degree+1) through a given set of points.
1187          * @param {Array} points Array consisting of JXG.Points.
1188          * @param {Number} order Order of the B-spline curve.
1189          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
1190          * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply
1191          * returning the length of the points array minus one.
1192          * @memberof JXG.Math.Numerics
1193          */
1194         bspline: function (points, order) {
1195             var knots, N = [],
1196                 _knotVector = function (n, k) {
1197                     var j,
1198                         kn = [];
1199 
1200                     for (j = 0; j < n + k + 1; j++) {
1201                         if (j < k) {
1202                             kn[j] = 0.0;
1203                         } else if (j <= n) {
1204                             kn[j] = j - k + 1;
1205                         } else {
1206                             kn[j] = n - k + 2;
1207                         }
1208                     }
1209 
1210                     return kn;
1211                 },
1212 
1213                 _evalBasisFuncs = function (t, kn, n, k, s) {
1214                     var i, j, a, b, den,
1215                         N = [];
1216 
1217                     if (kn[s] <= t && t < kn[s + 1]) {
1218                         N[s] = 1;
1219                     } else {
1220                         N[s] = 0;
1221                     }
1222 
1223                     for (i = 2; i <= k; i++) {
1224                         for (j = s - i + 1; j <= s; j++) {
1225                             if (j <= s - i + 1 || j < 0) {
1226                                 a = 0.0;
1227                             } else {
1228                                 a = N[j];
1229                             }
1230 
1231                             if (j >= s) {
1232                                 b = 0.0;
1233                             } else {
1234                                 b = N[j + 1];
1235                             }
1236 
1237                             den = kn[j + i - 1] - kn[j];
1238 
1239                             if (den === 0) {
1240                                 N[j] = 0;
1241                             } else {
1242                                 N[j] = (t - kn[j]) / den * a;
1243                             }
1244 
1245                             den = kn[j + i] - kn[j + 1];
1246 
1247                             if (den !== 0) {
1248                                 N[j] += (kn[j + i] - t) / den * b;
1249                             }
1250                         }
1251                     }
1252                     return N;
1253                 },
1254                 /** @ignore */
1255                 makeFct = function (which) {
1256                     return function (t, suspendedUpdate) {
1257                         var y, j, s,
1258                             len = points.length,
1259                             n = len - 1,
1260                             k = order;
1261 
1262                         if (n <= 0) {
1263                             return NaN;
1264                         }
1265 
1266                         if (n + 2 <= k) {
1267                             k = n + 1;
1268                         }
1269 
1270                         if (t <= 0) {
1271                             return points[0][which]();
1272                         }
1273 
1274                         if (t >= n - k + 2) {
1275                             return points[n][which]();
1276                         }
1277 
1278                         s = Math.floor(t) + k - 1;
1279                         knots = _knotVector(n, k);
1280                         N = _evalBasisFuncs(t, knots, n, k, s);
1281 
1282                         y = 0.0;
1283                         for (j = s - k + 1; j <= s; j++) {
1284                             if (j < len && j >= 0) {
1285                                 y += points[j][which]() * N[j];
1286                             }
1287                         }
1288 
1289                         return y;
1290                     };
1291                 };
1292 
1293             return [makeFct('X'), makeFct('Y'), 0,
1294                 function () {
1295                     return points.length - 1;
1296                 }];
1297         },
1298 
1299         /**
1300          * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve}
1301          * and {@link JXG.Curve#hasPoint}.
1302          * @param {function} f Function in one variable to be differentiated.
1303          * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a
1304          * method of an object and contains a reference to its parent object via "this".
1305          * @returns {function} Derivative function of a given function f.
1306          * @memberof JXG.Math.Numerics
1307          */
1308         D: function (f, obj) {
1309             var h = 0.00001,
1310                 h2 = 1.0 / (h * 2.0);
1311 
1312             if (!Type.exists(obj)) {
1313                 return function (x, suspendUpdate) {
1314                     return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2;
1315                 };
1316             }
1317 
1318             return function (x, suspendUpdate) {
1319                 return (f.apply(obj, [x + h, suspendUpdate]) - f.apply(obj, [x - h, suspendUpdate])) * h2;
1320             };
1321         },
1322 
1323         /**
1324          * Helper function to create curve which displays Riemann sums.
1325          * Compute coordinates for the rectangles showing the Riemann sum.
1326          * @param {function} f Function, whose integral is approximated by the Riemann sum.
1327          * @param {Number} n number of rectangles.
1328          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'.
1329          * @param {Number} start Left border of the approximation interval
1330          * @param {Number} end Right border of the approximation interval
1331          * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This
1332          * 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
1333          * rectangles.
1334          * @memberof JXG.Math.Numerics
1335          */
1336         riemann: function (f, n, type, start, end) {
1337             var i, x1, y1, delta1, delta,
1338                 xarr = [],
1339                 yarr = [],
1340                 j = 0,
1341                 x = start, y,
1342                 sum = 0;
1343 
1344             n = Math.round(n);
1345 
1346             xarr[j] = x;
1347             yarr[j] = 0.0;
1348 
1349             if (n > 0) {
1350                 delta = (end - start) / n;
1351                 // for 'lower' and 'upper'
1352                 delta1 = delta * 0.01;
1353 
1354                 for (i = 0; i < n; i++) {
1355                     if (type === 'right') {
1356                         y = f(x + delta);
1357                     } else if (type === 'middle') {
1358                         y = f(x + delta * 0.5);
1359                     } else if (type === 'left' || type === 'trapezoidal') {
1360                         y = f(x);
1361                     } else if (type === 'lower') {
1362                         y = f(x);
1363 
1364                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1365                             y1 = f(x1);
1366 
1367                             if (y1 < y) {
1368                                 y = y1;
1369                             }
1370                         }
1371 
1372                         y1 = f(x + delta);
1373                         if (y1 < y) {
1374                             y = y1;
1375                         }
1376                     } else if (type === 'upper') {
1377                         y = f(x);
1378 
1379                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1380                             y1 = f(x1);
1381 
1382                             if (y1 > y) {
1383                                 y = y1;
1384                             }
1385                         }
1386 
1387                         y1 = f(x + delta);
1388                         if (y1 > y) {
1389                             y = y1;
1390                         }
1391                     } else if (type === 'random') {
1392                         y = f(x + delta * Math.random());
1393                     } else if (type === 'simpson') {
1394                         y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0;
1395                     } else {
1396                         y = f(x);  // default is lower
1397                     }
1398 
1399                     j += 1;
1400                     xarr[j] = x;
1401                     yarr[j] = y;
1402                     j += 1;
1403                     x += delta;
1404 
1405                     if (type === 'trapezoidal') {
1406                         y = f(x);
1407                     }
1408 
1409                     xarr[j] = x;
1410                     yarr[j] = y;
1411                     j += 1;
1412                     xarr[j] = x;
1413                     yarr[j] = 0.0;
1414                     sum += y * delta;
1415                 }
1416             }
1417             return [xarr, yarr, sum];
1418         },
1419 
1420         /**
1421          * Approximate the integral by Riemann sums.
1422          * Compute the area described by the riemann sum rectangles.
1423          * @deprecated Replaced by JXG.Curve.Value(), see {@link JXG.Curve#riemannsum}
1424          * @param {function} f Function, whose integral is approximated by the Riemann sum.
1425          * @param {Number} n number of rectangles.
1426          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'.
1427          *
1428          * @param {Number} start Left border of the approximation interval
1429          * @param {Number} end Right border of the approximation interval
1430          * @returns {Number} The sum of the areas of the rectangles.
1431          * @memberof JXG.Math.Numerics
1432          */
1433         riemannsum: function (f, n, type, start, end) {
1434             var i, x1, y1, delta1, delta, y,
1435                 sum = 0.0,
1436                 x = start;
1437 
1438             // this function looks very similar to this.riemann... maybe there is some merge potential?
1439 
1440             n = Math.floor(n);
1441 
1442             if (n > 0) {
1443                 delta = (end - start) / n;
1444                 // for 'lower' and 'upper'
1445                 delta1 = delta * 0.01;
1446 
1447                 for (i = 0; i < n; i++) {
1448                     if (type === 'right') {
1449                         y = f(x + delta);
1450                     } else if (type === 'middle') {
1451                         y = f(x + delta * 0.5);
1452                     } else if (type === 'trapezoidal') {
1453                         y = 0.5 * (f(x + delta) + f(x));
1454                     } else if (type === 'left') {
1455                         y = f(x);
1456                     } else if (type === 'lower') {
1457                         y = f(x);
1458                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1459                             y1 = f(x1);
1460 
1461                             if (y1 < y) {
1462                                 y = y1;
1463                             }
1464                         }
1465 
1466                         y1 = f(x + delta);
1467                         if (y1 < y) {
1468                             y = y1;
1469                         }
1470                     } else if (type === 'upper') {
1471                         y = f(x);
1472 
1473                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1474                             y1 = f(x1);
1475 
1476                             if (y1 > y) {
1477                                 y = y1;
1478                             }
1479                         }
1480 
1481                         y1 = f(x + delta);
1482                         if (y1 > y) {
1483                             y = y1;
1484                         }
1485                     } else if (type === 'random') {
1486                         y = f(x + delta * Math.random());
1487                     } else if (type === 'simpson') {
1488                         y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0;
1489                     } else {
1490                         y = f(x);  // default is lower
1491                     }
1492 
1493                     sum += delta * y;
1494                     x += delta;
1495                 }
1496             }
1497 
1498             return sum;
1499         },
1500 
1501         /**
1502          * Solve initial value problems numerically using Runge-Kutta-methods.
1503          * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
1504          * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
1505          * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
1506          * <pre>
1507          * {
1508          *     s: <Number>,
1509          *     A: <matrix>,
1510          *     b: <Array>,
1511          *     c: <Array>
1512          * }
1513          * </pre>
1514          * 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
1515          * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array.
1516          * @param {Array} I Interval on which to integrate.
1517          * @param {Number} N Number of evaluation points.
1518          * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
1519          * 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
1520          * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has.
1521          * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
1522          * @example
1523          * // A very simple autonomous system dx(t)/dt = x(t);
1524          * function f(t, x) {
1525          *     return x;
1526          * }
1527          *
1528          * // Solve it with initial value x(0) = 1 on the interval [0, 2]
1529          * // with 20 evaluation points.
1530          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
1531          *
1532          * // Prepare data for plotting the solution of the ode using a curve.
1533          * var dataX = [];
1534          * var dataY = [];
1535          * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
1536          * for(var i=0; i<data.length; i++) {
1537          *     dataX[i] = i*h;
1538          *     dataY[i] = data[i][0];
1539          * }
1540          * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
1541          * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
1542          * <script type="text/javascript">
1543          * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
1544          * function f(t, x) {
1545          *     // we have to copy the value.
1546          *     // return x; would just return the reference.
1547          *     return [x[0]];
1548          * }
1549          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
1550          * var dataX = [];
1551          * var dataY = [];
1552          * var h = 0.1;
1553          * for(var i=0; i<data.length; i++) {
1554          *     dataX[i] = i*h;
1555          *     dataY[i] = data[i][0];
1556          * }
1557          * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
1558          * </script><pre>
1559          * @memberof JXG.Math.Numerics
1560          */
1561         rungeKutta: function (butcher, x0, I, N, f) {
1562             var e, i, j, k, l, s,
1563                 x = [],
1564                 y = [],
1565                 h = (I[1] - I[0]) / N,
1566                 t = I[0],
1567                 dim = x0.length,
1568                 result = [],
1569                 r = 0;
1570 
1571             if (Type.isString(butcher)) {
1572                 butcher = predefinedButcher[butcher] || predefinedButcher.euler;
1573             }
1574             s = butcher.s;
1575 
1576             // don't change x0, so copy it
1577             for (e = 0; e < dim; e++) {
1578                 x[e] = x0[e];
1579             }
1580 
1581             for (i = 0; i < N; i++) {
1582                 // Optimization doesn't work for ODEs plotted using time
1583                 //        if((i % quotient == 0) || (i == N-1)) {
1584                 result[r] = [];
1585                 for (e = 0; e < dim; e++) {
1586                     result[r][e] = x[e];
1587                 }
1588 
1589                 r += 1;
1590                 k = [];
1591 
1592                 for (j = 0; j < s; j++) {
1593                     // init y = 0
1594                     for (e = 0; e < dim; e++) {
1595                         y[e] = 0.0;
1596                     }
1597 
1598 
1599                     // Calculate linear combination of former k's and save it in y
1600                     for (l = 0; l < j; l++) {
1601                         for (e = 0; e < dim; e++) {
1602                             y[e] += (butcher.A[j][l]) * h * k[l][e];
1603                         }
1604                     }
1605 
1606                     // add x(t) to y
1607                     for (e = 0; e < dim; e++) {
1608                         y[e] += x[e];
1609                     }
1610 
1611                     // calculate new k and add it to the k matrix
1612                     k.push(f(t + butcher.c[j] * h, y));
1613                 }
1614 
1615                 // init y = 0
1616                 for (e = 0; e < dim; e++) {
1617                     y[e] = 0.0;
1618                 }
1619 
1620                 for (l = 0; l < s; l++) {
1621                     for (e = 0; e < dim; e++) {
1622                         y[e] += butcher.b[l] * k[l][e];
1623                     }
1624                 }
1625 
1626                 for (e = 0; e < dim; e++) {
1627                     x[e] = x[e] + h * y[e];
1628                 }
1629 
1630                 t += h;
1631             }
1632 
1633             return result;
1634         },
1635 
1636 
1637         /**
1638          * Maximum number of iterations in {@link JXG.Math.Numerics#fzero}
1639          * @type Number
1640          * @default 80
1641          * @memberof JXG.Math.Numerics
1642          */
1643         maxIterationsRoot: 80,
1644 
1645         /**
1646          * Maximum number of iterations in {@link JXG.Math.Numerics#fminbr}
1647          * @type Number
1648          * @default 500
1649          * @memberof JXG.Math.Numerics
1650          */
1651         maxIterationsMinimize: 500,
1652 
1653         /**
1654          *
1655          * Find zero of an univariate function f.
1656          * @param {function} f Function, whose root is to be found
1657          * @param {Array,Number} x0  Start value or start interval enclosing the root
1658          * @param {Object} object Parent object in case f is method of it
1659          * @returns {Number} the approximation of the root
1660          * Algorithm:
1661          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
1662          *  computations. M., Mir, 1980, p.180 of the Russian edition
1663          *
1664          * If x0 is an array containing lower and upper bound for the zero
1665          * algorithm 748 is applied. Otherwise, if x0 is a number,
1666          * the algorithm tries to bracket a zero of f starting from x0.
1667          * If this fails, we fall back to Newton's method.
1668          * @memberof JXG.Math.Numerics
1669          */
1670         fzero: function (f, x0, object) {
1671             var a, b, c,
1672                 fa, fb, fc,
1673                 aa, blist, i, len, u, fu,
1674                 prev_step, t1, cb, t2,
1675                 // Actual tolerance
1676                 tol_act,
1677                 // Interpolation step is calculated in the form p/q; division
1678                 // operations is delayed until the last moment
1679                 p, q,
1680                 // Step at this iteration
1681                 new_step,
1682                 eps = Mat.eps,
1683                 maxiter = this.maxIterationsRoot,
1684                 niter = 0,
1685                 nfev = 0;
1686 
1687             if (Type.isArray(x0)) {
1688                 if (x0.length < 2) {
1689                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
1690                 }
1691 
1692                 a = x0[0];
1693                 fa = f.call(object, a);
1694                 nfev += 1;
1695                 b = x0[1];
1696                 fb = f.call(object, b);
1697                 nfev += 1;
1698             } else {
1699                 a = x0;
1700                 fa = f.call(object, a);
1701                 nfev += 1;
1702 
1703                 // Try to get b.
1704                 if (a === 0) {
1705                     aa = 1;
1706                 } else {
1707                     aa = a;
1708                 }
1709 
1710                 blist = [0.9 * aa, 1.1 * aa, aa - 1, aa + 1, 0.5 * aa, 1.5 * aa, -aa, 2 * aa, -10 * aa, 10 * aa];
1711                 len = blist.length;
1712 
1713                 for (i = 0; i < len; i++) {
1714                     b = blist[i];
1715                     fb = f.call(object, b);
1716                     nfev += 1;
1717 
1718                     if (fa * fb <= 0) {
1719                         break;
1720                     }
1721                 }
1722                 if (b < a) {
1723                     u = a;
1724                     a = b;
1725                     b = u;
1726 
1727                     fu = fa;
1728                     fa = fb;
1729                     fb = fu;
1730                 }
1731             }
1732 
1733             if (fa * fb > 0) {
1734                 // Bracketing not successful, fall back to Newton's method or to fminbr
1735                 if (Type.isArray(x0)) {
1736                     return this.fminbr(f, [a, b], object);
1737                 }
1738 
1739                 return this.Newton(f, a, object);
1740             }
1741 
1742             // OK, we have enclosed a zero of f.
1743             // Now we can start Brent's method
1744 
1745             c = a;
1746             fc = fa;
1747 
1748             // Main iteration loop
1749             while (niter < maxiter) {
1750                 // Distance from the last but one to the last approximation
1751                 prev_step = b - a;
1752 
1753                 // Swap data for b to be the best approximation
1754                 if (Math.abs(fc) < Math.abs(fb)) {
1755                     a = b;
1756                     b = c;
1757                     c = a;
1758 
1759                     fa = fb;
1760                     fb = fc;
1761                     fc = fa;
1762                 }
1763                 tol_act = 2 * eps * Math.abs(b) + eps * 0.5;
1764                 new_step = (c - b) * 0.5;
1765 
1766                 if (Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps) {
1767                     //  Acceptable approx. is found
1768                     return b;
1769                 }
1770 
1771                 // Decide if the interpolation can be tried
1772                 // If prev_step was large enough and was in true direction Interpolatiom may be tried
1773                 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) {
1774                     cb = c - b;
1775 
1776                     // If we have only two distinct points linear interpolation can only be applied
1777                     if (a === c) {
1778                         t1 = fb / fa;
1779                         p = cb * t1;
1780                         q = 1.0 - t1;
1781                     // Quadric inverse interpolation
1782                     } else {
1783                         q = fa / fc;
1784                         t1 = fb / fc;
1785                         t2 = fb / fa;
1786 
1787                         p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0));
1788                         q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
1789                     }
1790 
1791                     // p was calculated with the opposite sign; make p positive
1792                     if (p > 0) {
1793                         q = -q;
1794                     // and assign possible minus to q
1795                     } else {
1796                         p = -p;
1797                     }
1798 
1799                     // If b+p/q falls in [b,c] and isn't too large it is accepted
1800                     // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
1801                     if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) &&
1802                             p < Math.abs(prev_step * q * 0.5)) {
1803                         new_step = p / q;
1804                     }
1805                 }
1806 
1807                 // Adjust the step to be not less than tolerance
1808                 if (Math.abs(new_step) < tol_act) {
1809                     if (new_step > 0) {
1810                         new_step = tol_act;
1811                     } else {
1812                         new_step = -tol_act;
1813                     }
1814                 }
1815 
1816                 // Save the previous approx.
1817                 a = b;
1818                 fa = fb;
1819                 b += new_step;
1820                 fb = f.call(object, b);
1821                 // Do step to a new approxim.
1822                 nfev += 1;
1823 
1824                 // Adjust c for it to have a sign opposite to that of b
1825                 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
1826                     c = a;
1827                     fc = fa;
1828                 }
1829                 niter++;
1830             } // End while
1831 
1832             return b;
1833         },
1834 
1835 
1836         /**
1837          *
1838          * Find minimum of an univariate function f.
1839          * @param {function} f Function, whose minimum is to be found
1840          * @param {Array} x0  Start interval enclosing the minimum
1841          * @param {Object} context Parent object in case f is method of it
1842          * @returns {Number} the approximation of the minimum value position
1843          * Algorithm:
1844          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
1845          *  computations. M., Mir, 1980, p.180 of the Russian edition
1846          * x0
1847          * @memberof JXG.Math.Numerics
1848          **/
1849         fminbr: function (f, x0, context) {
1850             var a, b, x, v, w,
1851                 fx, fv, fw,
1852                 range, middle_range, tol_act, new_step,
1853                 p, q, t, ft,
1854                 // Golden section ratio
1855                 r = (3.0 - Math.sqrt(5.0)) * 0.5,
1856                 tol = Mat.eps,
1857                 sqrteps = Mat.eps, //Math.sqrt(Mat.eps),
1858                 maxiter = this.maxIterationsMinimize,
1859                 niter = 0,
1860                 nfev = 0;
1861 
1862             if (!Type.isArray(x0) || x0.length < 2) {
1863                 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two.");
1864             }
1865 
1866             a = x0[0];
1867             b = x0[1];
1868             v = a + r * (b - a);
1869             fv = f.call(context, v);
1870 
1871             // First step - always gold section
1872             nfev += 1;
1873             x = v;
1874             w = v;
1875             fx = fv;
1876             fw = fv;
1877 
1878             while (niter < maxiter) {
1879                 // Range over the interval in which we are looking for the minimum
1880                 range = b - a;
1881                 middle_range = (a + b) * 0.5;
1882 
1883                 // Actual tolerance
1884                 tol_act = sqrteps * Math.abs(x) + tol / 3.0;
1885 
1886                 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) {
1887                     // Acceptable approx. is found
1888                     return x;
1889                 }
1890 
1891                 // Obtain the golden section step
1892                 new_step = r * (x < middle_range ? b - x : a - x);
1893 
1894                 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried
1895                 if (Math.abs(x - w) >= tol_act) {
1896                     // Interpolation step is calculated as p/q;
1897                     // division operation is delayed until last moment
1898                     t = (x - w) * (fx - fv);
1899                     q = (x - v) * (fx - fw);
1900                     p = (x - v) * q - (x - w) * t;
1901                     q = 2 * (q - t);
1902 
1903                     if (q > 0) {                        // q was calculated with the op-
1904                         p = -p;                         // posite sign; make q positive
1905                     } else {                            // and assign possible minus to
1906                         q = -q;                         // p
1907                     }
1908                     if (Math.abs(p) < Math.abs(new_step * q) &&     // If x+p/q falls in [a,b]
1909                             p > q * (a - x + 2 * tol_act) &&        //  not too close to a and
1910                             p < q * (b - x - 2 * tol_act)) {        // b, and isn't too large
1911                         new_step = p / q;                          // it is accepted
1912                     }
1913                     // If p/q is too large then the
1914                     // golden section procedure can
1915                     // reduce [a,b] range to more
1916                     // extent
1917                 }
1918 
1919                 // Adjust the step to be not less than tolerance
1920                 if (Math.abs(new_step) < tol_act) {
1921                     if (new_step > 0) {
1922                         new_step = tol_act;
1923                     } else {
1924                         new_step = -tol_act;
1925                     }
1926                 }
1927 
1928                 // Obtain the next approximation to min
1929                 // and reduce the enveloping range
1930 
1931                 // Tentative point for the min
1932                 t = x + new_step;
1933                 ft = f.call(context, t);
1934                 nfev += 1;
1935 
1936                 // t is a better approximation
1937                 if (ft <= fx) {
1938                     // Reduce the range so that t would fall within it
1939                     if (t < x) {
1940                         b = x;
1941                     } else {
1942                         a = x;
1943                     }
1944 
1945                     // Assign the best approx to x
1946                     v = w;
1947                     w = x;
1948                     x = t;
1949 
1950                     fv = fw;
1951                     fw = fx;
1952                     fx = ft;
1953                 // x remains the better approx
1954                 } else {
1955                     // Reduce the range enclosing x
1956                     if (t < x) {
1957                         a = t;
1958                     } else {
1959                         b = t;
1960                     }
1961 
1962                     if (ft <= fw || w === x) {
1963                         v = w;
1964                         w = t;
1965                         fv = fw;
1966                         fw = ft;
1967                     } else if (ft <= fv || v === x || v === w) {
1968                         v = t;
1969                         fv = ft;
1970                     }
1971                 }
1972                 niter += 1;
1973             }
1974 
1975             return x;
1976         },
1977 
1978         /**
1979          * Implements the Ramer-Douglas-Peucker algorithm.
1980          * It discards points which are not necessary from the polygonal line defined by the point array
1981          * pts. The computation is done in screen coordinates.
1982          * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points.
1983          * @param {Array} pts Array of {@link JXG.Coords}
1984          * @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>.
1985          * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points.
1986          * @memberof JXG.Math.Numerics
1987          */
1988         RamerDouglasPeucker: function (pts, eps) {
1989             var newPts = [], i, k, len,
1990 
1991                 /**
1992                  * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeucker}.
1993                  * It searches for the point between index i and j which
1994                  * has the largest distance from the line between the points i and j.
1995                  * @param {Array} pts Array of {@link JXG.Coords}
1996                  * @param {Number} i Index of a point in pts
1997                  * @param {Number} j Index of a point in pts
1998                  * @ignore
1999                  * @private
2000                  */
2001                 findSplit = function (pts, i, j) {
2002                     var d, k, ci, cj, ck,
2003                         x0, y0, x1, y1,
2004                         den, lbda,
2005                         dist = 0,
2006                         f = i;
2007 
2008                     if (j - i < 2) {
2009                         return [-1.0, 0];
2010                     }
2011 
2012                     ci = pts[i].scrCoords;
2013                     cj = pts[j].scrCoords;
2014 
2015                     if (isNaN(ci[1] + ci[2])) {
2016                         return [NaN, i];
2017                     }
2018                     if (isNaN(cj[1] + cj[2])) {
2019                         return [NaN, j];
2020                     }
2021 
2022                     for (k = i + 1; k < j; k++) {
2023                         ck = pts[k].scrCoords;
2024                         if (isNaN(ck[1] + ck[2])) {
2025                             return [NaN, k];
2026                         }
2027 
2028                         x0 = ck[1] - ci[1];
2029                         y0 = ck[2] - ci[2];
2030                         x1 = cj[1] - ci[1];
2031                         y1 = cj[2] - ci[2];
2032                         den = x1 * x1 + y1 * y1;
2033 
2034                         if (den >= Mat.eps) {
2035                             lbda = (x0 * x1 + y0 * y1) / den;
2036 
2037                             if (lbda < 0.0) {
2038                                 lbda = 0.0;
2039                             } else if (lbda > 1.0) {
2040                                 lbda = 1.0;
2041                             }
2042 
2043                             x0 = x0 - lbda * x1;
2044                             y0 = y0 - lbda * y1;
2045                             d = x0 * x0 + y0 * y0;
2046                         } else {
2047                             lbda = 0.0;
2048                             d = x0 * x0 + y0 * y0;
2049                         }
2050 
2051                         if (d > dist) {
2052                             dist = d;
2053                             f = k;
2054                         }
2055                     }
2056                     return [Math.sqrt(dist), f];
2057                 },
2058 
2059                 /**
2060                  * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeucker}.
2061                  * It runs recursively through the point set and searches the
2062                  * point which has the largest distance from the line between the first point and
2063                  * the last point. If the distance from the line is greater than eps, this point is
2064                  * included in our new point set otherwise it is discarded.
2065                  * If it is taken, we recursively apply the subroutine to the point set before
2066                  * and after the chosen point.
2067                  * @param {Array} pts Array of {@link JXG.Coords}
2068                  * @param {Number} i Index of an element of pts
2069                  * @param {Number} j Index of an element of pts
2070                  * @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>.
2071                  * @param {Array} newPts Array of {@link JXG.Coords}
2072                  * @ignore
2073                  * @private
2074                  */
2075                 RDP = function (pts, i, j, eps, newPts) {
2076                     var result = findSplit(pts, i, j),
2077                         k = result[1];
2078 
2079                     if (isNaN(result[0])) {
2080                         RDP(pts, i, k - 1, eps, newPts);
2081                         newPts.push(pts[k]);
2082                         do {
2083                             ++k;
2084                         } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2]));
2085                         if (k <= j) {
2086                             newPts.push(pts[k]);
2087                         }
2088                         RDP(pts, k + 1, j, eps, newPts);
2089                     } else if (result[0] > eps) {
2090                         RDP(pts, i, k, eps, newPts);
2091                         RDP(pts, k, j, eps, newPts);
2092                     } else {
2093                         newPts.push(pts[j]);
2094                     }
2095                 };
2096 
2097             len = pts.length;
2098 
2099             // Search for the left most point woithout NaN coordinates
2100             i = 0;
2101             while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) {
2102                 i += 1;
2103             }
2104             // Search for the right most point woithout NaN coordinates
2105             k = len - 1;
2106             while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) {
2107                 k -= 1;
2108             }
2109 
2110             // Only proceed if something is left
2111             if (!(i > k || i === len)) {
2112                 newPts[0] = pts[i];
2113                 RDP(pts, i, k, eps, newPts);
2114             }
2115 
2116             return newPts;
2117         },
2118 
2119         /**
2120          * Old name for the implementation of the Ramer-Douglas-Peucker algorithm.
2121          * @deprecated Use {@link JXG.Math.Numerics#RamerDouglasPeucker}
2122          * @memberof JXG.Math.Numerics
2123          */
2124         RamerDouglasPeuker: function (pts, eps) {
2125             return this.RamerDouglasPeucker(pts, eps);
2126         }
2127     };
2128 
2129     return Mat.Numerics;
2130 });
2131