1 /*
  2  Copyright 2008,2009
  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: you can redistribute it and/or modify
 13  it under the terms of the GNU Lesser General Public License as published by
 14  the Free Software Foundation, either version 3 of the License, or
 15  (at your option) any later version.
 16 
 17  JSXGraph is distributed in the hope that it will be useful,
 18  but WITHOUT ANY WARRANTY; without even the implied warranty of
 19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 20  GNU Lesser General Public License for more details.
 21 
 22  You should have received a copy of the GNU Lesser General Public License
 23  along with JSXGraph. If not, see <http://www.gnu.org/licenses/>.
 24  */
 25 
 26 /**
 27  * @fileoverview In this file the namespace Math.Numerics is defined, which holds numerical
 28  * algorithms for solving linear equations etc.
 29  * @author graphjs
 30  */
 31 
 32 
 33 /**
 34  * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables.
 35  * @namespace
 36  */
 37 JXG.Math.Numerics = (function(JXG, Math) {
 38 
 39     // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order).
 40     var predefinedButcher = {
 41         rk4: {
 42             s: 4,
 43             A: [
 44                 [ 0,  0,  0, 0],
 45                 [0.5, 0,  0, 0],
 46                 [ 0, 0.5, 0, 0],
 47                 [ 0,  0,  1, 0]
 48             ],
 49             b: [1. / 6., 1. / 3., 1. / 3., 1. / 6.],
 50             c: [0, 0.5, 0.5, 1]
 51         },
 52         heun: {
 53             s: 2,
 54             A: [
 55                 [0, 0],
 56                 [1, 0]
 57             ],
 58             b: [0.5, 0.5],
 59             c: [0, 1]
 60         },
 61         euler: {
 62             s: 1,
 63             A: [
 64                 [0]
 65             ],
 66             b: [1],
 67             c: [0]
 68         }
 69     };
 70 
 71 
 72     /** @lends JXG.Math.Numerics */
 73     return {
 74 
 75         /**
 76          * Solves a system of linear equations given by A and b using the Gauss-Jordan-elimination.
 77          * The algorithm runs in-place. I.e. the entries of A and b are changed.
 78          * @param {Array} A Square matrix represented by an array of rows, containing the coefficients of the lineare equation system.
 79          * @param {Array} b A vector containing the linear equation system's right hand side.
 80          * @throws {Error} If a non-square-matrix is given or if b has not the right length or A's rank is not full.
 81          * @returns {Array} A vector that solves the linear equation system.
 82          */
 83         Gauss: function(A, b) {
 84             var eps = JXG.Math.eps,
 85                 // number of columns of A
 86                 n = A.length > 0 ? A[0].length : 0,
 87                 // copy the matrix to prevent changes in the original
 88                 Acopy,
 89                 // solution vector, to prevent changing b
 90                 x,
 91                 i, j, k,
 92                 // little helper to swap array elements
 93                 swap = function(i, j) {
 94                     var temp = this[i];
 95 
 96                     this[i] = this[j];
 97                     this[j] = temp;
 98                 };
 99 
100             if ((n !== b.length) || (n !== A.length))
101                 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.");
102 
103             // initialize solution vector
104             Acopy = new Array(n);
105             x = b.slice(0, n);
106             for (i = 0; i < n; i++) {
107                 Acopy[i] = A[i].slice(0, n);
108             }
109 
110             // Gauss-Jordan-elimination
111             for (j = 0; j < n; j++) {
112                 for (i = n - 1; i > j; i--) {
113                     // Is the element which is to eliminate greater than zero?
114                     if (Math.abs(Acopy[i][j]) > eps) {
115                         // Equals pivot element zero?
116                         if (Math.abs(Acopy[j][j]) < eps) {
117                             // At least numerically, so we have to exchange the rows
118                             swap.apply(Acopy, [i, j]);
119                             swap.apply(x, [i, j]);
120                         } else {
121                             // Saves the L matrix of the LR-decomposition. unnecessary.
122                             Acopy[i][j] /= Acopy[j][j];
123                             // Transform right-hand-side b
124                             x[i] -= Acopy[i][j] * x[j];
125                             // subtract the multiple of A[i][j] / A[j][j] of the j-th row from the i-th.
126                             for (k = j + 1; k < n; k ++) {
127                                 Acopy[i][k] -= Acopy[i][j] * Acopy[j][k];
128                             }
129                         }
130                     }
131                 }
132                 if (Math.abs(Acopy[j][j]) < eps) { // The absolute values of all coefficients below the j-th row in the j-th column are smaller than JXG.Math.eps.
133                     throw new Error("JXG.Math.Numerics.Gauss(): The given matrix seems to be singular.");
134                 }
135             }
136             this.backwardSolve(Acopy, x, true); // return Array
137 
138             return x;
139         },
140 
141         /**
142          * Solves a system of linear equations given by the right triangular matrix R and vector b.
143          * @param {Array} R Right triangular matrix represented by an array of rows. All entries a_(i,j) with i < j are ignored.
144          * @param {Array} b Right hand side of the linear equation system.
145          * @param {Boolean} [canModify=false] If true, the right hand side vector is allowed to be changed by this method.
146          * @returns {Array} An array representing a vector that solves the system of linear equations.
147          */
148         backwardSolve: function(R, b, canModify) {
149             var x, m, n, i, j;
150 
151             if (canModify) {
152                 x = b;
153             } else {
154                 x = b.slice(0, b.length);
155             }
156 
157             // m: number of rows of R
158             // n: number of columns of R
159             m = R.length;
160             n = R.length > 0 ? R[0].length : 0;
161             for (i = m - 1; i >= 0; i--) {
162                 for (j = n - 1; j > i; j--) {
163                     x[i] -= R[i][j] * x[j];
164                 }
165                 x[i] /= R[i][i];
166             }
167 
168             return x;
169         },
170         
171         /**
172          * @private
173          * Gauss-Bareiss algorithm to compute the 
174          * determinant of matrix without fractions.
175          * See H. Cohen, "A course in computational
176          * algebraic number thoery".
177          */
178         gaussBareiss: function(mat) {
179             var k, c, s, i, j, p, n, M, t, 
180                 eps = JXG.Math.eps;
181             
182             n = mat.length;
183             if (n<=0) { return 0; }
184             if (mat[0].length<n) { n = mat[0].length; }
185             
186             // Copy the input matrix to M
187             M = new Array(n);
188             for (i = 0; i < n; i++) {
189                 M[i] = mat[i].slice(0, n);
190             }
191             
192             c = 1;
193             s = 1;
194             for (k=0;k<n-1;k++) {
195                 p = M[k][k];
196                 if (Math.abs(p)<eps) {    // Pivot step
197                     for (i=0;i<n;i++) {
198                         if (Math.abs(M[i][k])>=eps) break
199                     }
200                     if (i==n) {      // No nonzero entry found in column k -> det(M) = 0
201                         return 0.0;
202                     } 
203                     for (j=k;j<n;j++) {   // swap row i and k partially
204                         t = M[i][j]; M[i][j] = M[k][j];  M[k][j] = t;
205                     }
206                     s = -s;
207                     p = M[k][k];
208                 }
209                 
210                 // Main step
211                 for (i=k+1;i<n;i++) {
212                     for (j=k+1;j<n;j++) {
213                         t = p*M[i][j] - M[i][k]*M[k][j];
214                         M[i][j] = t/c;
215                     }
216                 }
217                 c = p;
218             }
219             return s*M[n-1][n-1];
220         
221         },
222 
223         /** 
224          * Computes the determinant of a square nxn matrix with the
225          * Gauss-Bareiss algorithm.
226          * @param {Array} mat Matrix. 
227          * @returns {Number} The determinant pf the matrix mat. 
228                              The empty matrix returns 0.
229          */ 
230         det: function(mat) {
231             var n = mat.length;
232             if (n==2 && mat[0].length==2) {
233                 return mat[0][0]*mat[1][1] - mat[1][0]*mat[0][1];
234             } else {
235                 return this.gaussBareiss(mat);
236             }
237         },
238         
239         /**
240          * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method
241          * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990
242          * @param {Array} Ain A symmetric 3x3 matrix.
243          * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors.
244          */
245         Jacobi: function(Ain) {
246             var i,j,k,aa,si,co,tt,eps = JXG.Math.eps,
247                 sum = 0.0,
248                 ssum, amax,
249                 n = Ain.length,
250                 V = [
251                     [0,0,0],
252                     [0,0,0],
253                     [0,0,0]
254                 ],
255                 A = [
256                     [0,0,0],
257                     [0,0,0],
258                     [0,0,0]
259                 ], nloops=0;
260 
261             // Initialization. Set initial Eigenvectors.
262             for (i = 0; i < n; i++) {
263                 for (j = 0; j < n; j++) {
264                     V[i][j] = 0.0;
265                     A[i][j] = Ain[i][j];
266                     sum += Math.abs(A[i][j]);
267                 }
268                 V[i][i] = 1.0;
269             }
270             // Trivial problems
271             if (n == 1) {
272                 return [A,V];
273             }
274             if (sum <= 0.0) {
275                 return [A,V];
276             }
277 
278             sum /= (n * n);
279 
280             // Reduce matrix to diagonal
281             do {
282                 ssum = 0.0;
283                 amax = 0.0;
284                 for (j = 1; j < n; j++) {
285                     for (i = 0; i < j; i++) {
286                         // Check if A[i][j] is to be reduced
287                         aa = Math.abs(A[i][j]);
288                         if (aa > amax) {
289                             amax = aa;
290                         }
291                         ssum += aa;
292                         if (aa >= eps) {
293                             // calculate rotation angle
294                             aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5;
295                             si = Math.sin(aa);
296                             co = Math.cos(aa);
297                             // Modify 'i' and 'j' columns
298                             for (k = 0; k < n; k++) {
299                                 tt = A[k][i];
300                                 A[k][i] = co * tt + si * A[k][j];
301                                 A[k][j] = -si * tt + co * A[k][j];
302                                 tt = V[k][i];
303                                 V[k][i] = co * tt + si * V[k][j];
304                                 V[k][j] = -si * tt + co * V[k][j];
305                             }
306                             // Modify diagonal terms
307                             A[i][i] = co * A[i][i] + si * A[j][i];
308                             A[j][j] = -si * A[i][j] + co * A[j][j];
309                             A[i][j] = 0.0;
310                             // Make 'A' matrix symmetrical
311                             for (k = 0; k < n; k++) {
312                                 A[i][k] = A[k][i];
313                                 A[j][k] = A[k][j];
314                             }
315                             // A[i][j] made zero by rotation
316                         }
317                     }
318                 }
319                 nloops++;
320             } while (Math.abs(ssum) / sum > eps && nloops<2000);
321             return [A,V];
322         },
323 
324         /**
325          * Calculates the integral of function f over interval using Newton-Cotes-algorithm.
326          * @param {Array} interval The integration interval, e.g. [0, 3].
327          * @param {function} f A function which takes one argument of type number and returns a number.
328          * @param {object} [config={number_of_nodes:28,integration_type:'milne'}] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type
329          * with value being either 'trapez', 'simpson', or 'milne'.
330          * @returns {Number} Integral value of f over interval
331          * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use
332          * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4.
333          * @example
334          * function f(x) {
335          *   return x*x;
336          * }
337          *
338          * // calculates integral of <tt>f</tt> from 0 to 2.
339          * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f);
340          *
341          * // the same with an anonymous function
342          * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; });
343          *
344          * // use trapez rule with 16 nodes
345          * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f,
346          *                                   {number_of_nodes: 16, intergration_type: 'trapez'});
347          */
348         NewtonCotes: function(interval, f, config) {
349             var integral_value = 0.0,
350                 number_of_nodes = config && typeof config.number_of_nodes === 'number' ? config.number_of_nodes : 28,
351                 available_types = {trapez: true, simpson: true, milne: true},
352                 integration_type = config && config.integration_type && available_types.hasOwnProperty(config.integration_type) && available_types[config.integration_type] ? config.integration_type : 'milne',
353                 step_size = (interval[1] - interval[0]) / number_of_nodes,
354                 evaluation_point, i, number_of_intervals;
355 
356             switch (integration_type) {
357                 case 'trapez':
358                     integral_value = (f(interval[0]) + f(interval[1])) * 0.5;
359 
360                     evaluation_point = interval[0];
361                     for (i = 0; i < number_of_nodes - 1; i++) {
362                         evaluation_point += step_size;
363                         integral_value += f(evaluation_point);
364                     }
365                     integral_value *= step_size;
366 
367                     break;
368                 case 'simpson':
369                     if (number_of_nodes % 2 > 0) {
370                         throw new Error("JSXGraph:  INT_SIMPSON requires config.number_of_nodes dividable by 2.");
371                     }
372                     number_of_intervals = number_of_nodes / 2.0;
373                     integral_value = f(interval[0]) + f(interval[1]);
374                     evaluation_point = interval[0];
375                     for (i = 0; i < number_of_intervals - 1; i++) {
376                         evaluation_point += 2.0 * step_size;
377                         integral_value += 2.0 * f(evaluation_point);
378                     }
379                     evaluation_point = interval[0] - step_size;
380                     for (i = 0; i < number_of_intervals; i++) {
381                         evaluation_point += 2.0 * step_size;
382                         integral_value += 4.0 * f(evaluation_point);
383                     }
384                     integral_value *= step_size / 3.0;
385                     break;
386                 default:
387                     if (number_of_nodes % 4 > 0) {
388                         throw new Error("JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4");
389                     }
390                     number_of_intervals = number_of_nodes * 0.25;
391                     integral_value = 7.0 * (f(interval[0]) + f(interval[1]));
392                     evaluation_point = interval[0];
393                     for (i = 0; i < number_of_intervals - 1; i++) {
394                         evaluation_point += 4.0 * step_size;
395                         integral_value += 14.0 * f(evaluation_point);
396                     }
397                     evaluation_point = interval[0] - 3.0 * step_size;
398                     for (i = 0; i < number_of_intervals; i++) {
399                         evaluation_point += 4.0 * step_size;
400                         integral_value += 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size));
401                     }
402                     evaluation_point = interval[0] - 2.0 * step_size;
403                     for (i = 0; i < number_of_intervals; i++) {
404                         evaluation_point += 4.0 * step_size;
405                         integral_value += 12.0 * f(evaluation_point);
406                     }
407                     integral_value *= 2.0 * step_size / 45.0;
408             }
409             return integral_value;
410         },
411 
412         /**
413          * Integral of function f over interval.
414          * @param {Array} interval The integration interval, e.g. [0, 3].
415          * @param {function} f A function which takes one argument of type number and returns a number.
416          * @returns {Number} The value of the integral of f over interval
417          * @see JXG.Math.Numerics.NewtonCotes
418          */
419         I: function(interval, f) {
420             return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'});
421         },
422 
423         /**
424          * Newton's method to find roots of a funtion in one variable.
425          * @param {function} f We search for a solution of f(x)=0.
426          * @param {Number} x initial guess for the root, i.e. start value.
427          * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a
428          *                 method of an object and contains a reference to its parent object via "this".
429          * @returns {Number} A root of the function f.
430          */
431         Newton: function(f, x, object) {
432             var i = 0,
433                 h = JXG.Math.eps,
434                 newf = f.apply(object, [x]), // set "this" to "object" in f
435                 nfev = 1, 
436                 df;
437             if (JXG.isArray(x)) {  // For compatibility
438                 x = x[0];
439             }
440             while (i < 50 && Math.abs(newf) > h) {
441                 df = this.D(f, object)(x);  nfev += 2;
442                 if (Math.abs(df) > h) {
443                     x -= newf / df;
444                 } else {
445                     x += (Math.random() * 0.2 - 1.0);
446                 }
447                 newf = f.apply(object, [x]); nfev++;
448                 i++;
449             }
450             //JXG.debug("N nfev="+nfev);
451             return x;
452         },
453 
454         /**
455          * Abstract method to find roots of univariate functions.
456          * @param {function} f We search for a solution of f(x)=0.
457          * @param {Number} x initial guess for the root, i.e. staring value.
458          * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a
459          *                 method of an object and contains a reference to its parent object via "this".
460          * @returns {Number} A root of the function f.
461          */
462         root: function(f, x, object) {
463             //return this.Newton(f, x, object);
464             return this.fzero(f, x, object);
465         },
466 
467         /**
468          * Returns the Lagrange polynomials for curves with equidistant nodes, see
469          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
470          * SIAM Review, Vol 46, No 3, (2004) 501-517.
471          * The graph of the parametric curve [x(t),y(t)] runs through the given points.
472          * @param {Array} p Array of JXG.Points
473          * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve
474          * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero.
475          */
476         Neville: function(p) {
477             var w = [],
478                 makeFct = function(fun) {
479                     return function(t, suspendedUpdate) {
480                         var i, d, s,
481                             bin = JXG.Math.binomial,
482                             len = p.length,
483                             len1 = len - 1,
484                             num = 0.0,
485                             denom = 0.0;
486 
487                         if (!suspendedUpdate) {
488                             s = 1;
489                             for (i = 0; i < len; i++) {
490                                 w[i] = bin(len1, i) * s;
491                                 s *= (-1);
492                             }
493                         }
494 
495                         d = t;
496                         for (i = 0; i < len; i++) {
497                             if (d === 0) {
498                                 return p[i][fun]();
499                             } else {
500                                 s = w[i] / d;
501                                 d--;
502                                 num += p[i][fun]() * s;
503                                 denom += s;
504                             }
505                         }
506                         return num / denom;
507                     };
508                 },
509 
510                 xfct = makeFct('X'),
511                 yfct = makeFct('Y');
512 
513             return [xfct, yfct, 0, function() {
514                 return p.length - 1;
515             }];
516         },
517 
518         /**
519          * Calculates second derivatives at the given knots.
520          * @param {Array} x x values of knots
521          * @param {Array} y y values of knots
522          * @returns {Array} Second derivatives of the interpolated function at the knots.
523          * @see #splineEval
524          */
525         splineDef: function(x, y) {
526             var n = Math.min(x.length, y.length),
527                 pair, i, l,
528                 diag = [],
529                 z = [],
530                 data = [],
531                 dx = [],
532                 delta = [],
533                 F = [];
534 
535             if (n === 2) {
536                 return [0, 0];
537             }
538 
539             for (i = 0; i < n; i++) {
540                 pair = {X: x[i], Y: y[i]};
541                 data.push(pair);
542             }
543             data.sort(function (a, b) {
544                 return a.X - b.X;
545             });
546             for (i = 0; i < n; i++) {
547                 x[i] = data[i].X;
548                 y[i] = data[i].Y;
549             }
550 
551             for (i = 0; i < n - 1; i++) {
552                 dx.push(x[i + 1] - x[i]);
553             }
554             for (i = 0; i < n - 2; i++) {
555                 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i]));
556             }
557 
558             // ForwardSolve
559             diag.push(2 * (dx[0] + dx[1]));
560             z.push(delta[0]);
561 
562             for (i = 0; i < n - 3; i++) {
563                 l = dx[i + 1] / diag[i];
564                 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]);
565                 z.push(delta[i + 1] - l * z[i]);
566             }
567 
568             // BackwardSolve
569             F[n - 3] = z[n - 3] / diag[n - 3];
570             for (i = n - 4; i >= 0; i--) {
571                 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i];
572             }
573 
574             // Generate f''-Vector
575             for (i = n - 3; i >= 0; i--) {
576                 F[i + 1] = F[i];
577             }
578 
579             // natural cubic spline
580             F[0] = 0;
581             F[n - 1] = 0;
582             return F;
583         },
584 
585         /**
586          * Evaluate points on spline.
587          * @param {Number,Array} x0 A single float value or an array of values to evaluate
588          * @param {Array} x x values of knots
589          * @param {Array} y y values of knots
590          * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef}
591          * @see #splineDef
592          * @returns {Number,Array} A single value or an array, depending on what is given as x0.
593          */
594         splineEval: function(x0, x, y, F) {
595             var n = Math.min(x.length, y.length),
596                 l = 1,
597                 asArray = false,
598                 y0 = [],
599                 i, j, a, b, c, d, x_;
600 
601             // number of points to be evaluated
602             if (JXG.isArray(x0)) {
603                 l = x0.length;
604                 asArray = true;
605             } else
606                 x0 = [x0];
607 
608             for (i = 0; i < l; i++) {
609                 // is x0 in defining interval?
610                 if ((x0[i] < x[0]) || (x[i] > x[n - 1]))
611                     return NaN;
612 
613                 // determine part of spline in which x0 lies
614                 for (j = 1; j < n; j++) {
615                     if (x0[i] <= x[j])
616                         break;
617                 }
618                 j--;
619 
620                 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1];
621                 // determine the coefficients of the polynomial in this interval
622                 a = y[j];
623                 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]);
624                 c = F[j] / 2;
625                 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j]));
626                 // evaluate x0[i]
627                 x_ = x0[i] - x[j];
628                 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_);
629                 y0.push(a + (b + (c + d * x_) * x_) * x_);
630             }
631 
632             if (asArray)
633                 return y0;
634             else
635                 return y0[0];
636         },
637 
638         /**
639          * Generate a string containing the function term of a polynomial.
640          * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i.
641          * @param {Number} deg Degree of the polynomial
642          * @param {String} varname Name of the variable (usually 'x')
643          * @param {Number} prec Precision
644          * @returns {String} A string containg the function term of the polynomial.
645          */
646         generatePolynomialTerm: function(coeffs, deg, varname, prec) {
647             var t = [], i;
648             for (i = deg; i >= 0; i--) {
649                 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']);
650                 if (i > 1) {
651                     t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']);
652                 } else if (i === 1) {
653                     t = t.concat(['*', varname, ' + ']);
654                 }
655             }
656             return t.join('');
657         },
658 
659         /**
660          * Computes the polynomial through a given set of coordinates in Lagrange form.
661          * Returns the Lagrange polynomials, see
662          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
663          * SIAM Review, Vol 46, No 3, (2004) 501-517.
664          * @param {Array} p Array of JXG.Points
665          * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points.
666          */
667         lagrangePolynomial: function(p) {
668             var w = [],
669                 fct = function(x, suspendedUpdate) {
670                     var i, k, len, xi, s,
671                         num = 0, denom = 0,
672                         M, j;
673 
674                     len = p.length;
675                     if (!suspendedUpdate) {
676                         for (i = 0; i < len; i++) {
677                             w[i] = 1.0;
678                             xi = p[i].X();
679                             for (k = 0; k < len; k++) if (k != i) {
680                                 w[i] *= (xi - p[k].X());
681                             }
682                             w[i] = 1 / w[i];
683                         }
684 
685                         M = [];
686                         for (j = 0; j < len; j++) {
687                             M.push([1]);
688                         }
689                     }
690 
691                     for (i = 0; i < len; i++) {
692                         xi = p[i].X();
693                         if (x === xi) {
694                             return p[i].Y();
695                         } else {
696                             s = w[i] / (x - xi);
697                             denom += s;
698                             num += s * p[i].Y();
699                         }
700                     }
701                     return num / denom;
702                 };
703 
704             fct.getTerm = function() {
705                 return '';
706             };
707 
708             return fct;
709         },
710         
711         /**
712          * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve
713          * is uniformly parametrized.
714          * Two artificial control points at the beginning and the end are added.
715          * @param {Array} points Array consisting of JXG.Points.
716          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
717          * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply
718          * returning the length of the points array minus three.
719         */
720         CatmullRomSpline: function(points) {
721             var coeffs = [],
722                 p, first = {}, last = {}, // control point at the beginning and at the end
723                 
724                 makeFct = function(which) {
725                     return function(t, suspendedUpdate) {
726                         var len = points.length,
727                             s, c;
728 
729                         if (len < 2 ) { return NaN; }
730                         if (!suspendedUpdate) {
731                             first[which] = function() {return 2*points[0][which]()-points[1][which]();}
732                             last[which] = function() {return 2*points[len-1][which]()-points[len-2][which]();}
733                             p = [first].concat(points, [last]);
734                             coeffs[which] = [];
735                             for (s=0; s < len-1; s++) {
736                                 coeffs[which][s] = [
737                                     2*p[s+1][which](),
738                                      -p[s][which]() +   p[s+2][which](),
739                                     2*p[s][which]() - 5*p[s+1][which]() + 4*p[s+2][which]() - p[s+3][which](),
740                                      -p[s][which]() + 3*p[s+1][which]() - 3*p[s+2][which]() + p[s+3][which]()
741                                     ];
742                             }
743                         }
744                         len += 2;  // add the two control points 
745                         if (isNaN(t)) { return NaN; }
746                         // This is necessay for our advanced plotting algorithm:
747                         if (t < 0) {   
748                             return p[1][which]();
749                         } else if (t >= len - 3) {
750                             return p[len-2][which]();
751                         }
752 
753                         s = Math.floor(t);
754                         if (s==t) {
755                             return p[s][which]();
756                         }
757                         t -= s;
758                         c = coeffs[which][s];
759                         return 0.5*(((c[3]*t + c[2])*t + c[1])*t + c[0]);
760                     };
761                 };
762 
763             return [makeFct('X'),
764                 makeFct('Y'),
765                 0,
766                 function() {
767                     return points.length - 1;
768                 }
769             ];
770         },
771 
772         
773         /**
774          * Computes the regression polynomial of a given degree through a given set of coordinates.
775          * Returns the regression polynomial function.
776          * @param {Number} degree number, function or slider.
777          * Either
778          * @param {Array} dataX array containing the x-coordinates of the data set
779          * @param {Array} dataY array containing the y-coordinates of the data set,
780          * or
781          * @param {Array} dataX Array consisting of JXG.Points.
782          * @param {} dataY Ignored
783          * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree.
784          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
785          */
786         regressionPolynomial: function(degree, dataX, dataY) {
787             var coeffs,
788                 deg, dX, dY,
789                 inputType,
790                 fct,
791                 term = '';
792 
793             if (JXG.isPoint(degree) && typeof degree.Value == 'function') {  // Slider
794                 deg = function() {return degree.Value();};
795             } else if (JXG.isFunction(degree)) {
796                 deg = degree;
797             } else if (JXG.isNumber(degree)) {
798                 deg = function() {return degree;};
799             } else {
800                 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'.");
801             }
802 
803             if (arguments.length == 3 && JXG.isArray(dataX) && JXG.isArray(dataY)) {              // Parameters degree, dataX, dataY
804                 inputType = 0;
805             } else if (arguments.length == 2 && JXG.isArray(dataX) && dataX.length > 0 && JXG.isPoint(dataX[0])) {  // Parameters degree, point array
806                 inputType = 1;
807             } else {
808                 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters.");
809             }
810 
811             fct = function(x, suspendedUpdate) {
812                 var i, j, M, MT, y, B, c, s,
813                     d, len = dataX.length;                        // input data
814 
815                 d = Math.floor(deg());                      // input data
816                 if (!suspendedUpdate) {
817                     if (inputType === 1) {  // point list as input
818                         dX = [];
819                         dY = [];
820                         for (i = 0; i < len; i++) {
821                             dX[i] = dataX[i].X();
822                             dY[i] = dataX[i].Y();
823                         }
824                     }
825 
826                     if (inputType === 0) {  // check for functions
827                         dX = [];
828                         dY = [];
829                         for (i = 0; i < len; i++) {
830                             if (JXG.isFunction(dataX[i]))
831                                 dX.push(dataX[i]());
832                             else
833                                 dX.push(dataX[i]);
834                             if (JXG.isFunction(dataY[i]))
835                                 dY.push(dataY[i]());
836                             else
837                                 dY.push(dataY[i]);
838                         }
839                     }
840 
841                     M = [];
842                     for (j = 0; j < len; j++) {
843                         M.push([1]);
844                     }
845                     for (i = 1; i <= d; i++) {
846                         for (j = 0; j < len; j++) {
847                             M[j][i] = M[j][i - 1] * dX[j];      // input data
848                         }
849                     }
850 
851                     y = dY;                                 // input data
852                     MT = JXG.Math.transpose(M);
853                     B = JXG.Math.matMatMult(MT, M);
854                     c = JXG.Math.matVecMult(MT, y);
855                     coeffs = JXG.Math.Numerics.Gauss(B, c);
856                     term = JXG.Math.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3);
857                 }
858 
859                 // Horner's scheme to evaluate polynomial
860                 s = coeffs[d];
861                 for (i = d - 1; i >= 0; i--) {
862                     s = (s * x + coeffs[i]);
863                 }
864                 return s;
865             };
866 
867             fct.getTerm = function() {
868                 return term;
869             };
870 
871             return fct;
872         },
873 
874         /**
875          * Computes the cubic Bezier curve through a given set of points.
876          * @param {Array} points Array consisting of 3*k+1 JXG.Points.
877          * The points at position k with k mod 3 = 0 are the data points,
878          * points at position k with k mod 3 = 1 or 2 are the control points.
879          * @returns {Array} An array consisting of two functions of one parameter t which return the
880          * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting
881          * no parameters and returning one third of the length of the points. 
882          */
883         bezier: function(points) {
884             var len,
885                 makeFct = function(which) {
886                     return function(t, suspendedUpdate) {
887                         var z = Math.floor(t) * 3,
888                             t0 = t % 1,
889                             t1 = 1 - t0;
890 
891                         if (!suspendedUpdate) {
892                             len = Math.floor(points.length / 3);
893                         }
894 
895                         if (t < 0) {
896                             return points[0][which]();
897                         }
898                         if (t >= len) {
899                             return points[points.length - 1][which]();
900                         }
901                         if (isNaN(t)) {
902                             return NaN;
903                         }
904                         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;
905                     };
906                 };
907 
908             return [
909                 makeFct('X'),
910                 makeFct('Y'),
911                 0,
912                 function() {
913                     return Math.floor(points.length / 3);
914                 }];
915         },
916 
917         /**
918          * Computes the B-spline curve of order k (order = degree+1) through a given set of points.
919          * @param {Array} points Array consisting of JXG.Points.
920          * @param {Number} order Order of the B-spline curve.
921          * @todo closed B-spline curves
922          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
923          * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply
924          * returning the length of the points array minus one.
925          */
926         bspline: function(points, order) {
927             var knots, N = [],
928                 _knotVector = function(n, k) {
929                     var j, kn = [];
930                     for (j = 0; j < n + k + 1; j++) {
931                         if (j < k) {
932                             kn[j] = 0.0;
933                         } else if (j <= n) {
934                             kn[j] = j - k + 1;
935                         } else {
936                             kn[j] = n - k + 2;
937                         }
938                     }
939                     return kn;
940                 },
941 
942                 _evalBasisFuncs = function(t, kn, n, k, s) {
943                     var i,j,a,b,den,
944                         N = [];
945 
946                     if (kn[s] <= t && t < kn[s + 1]) {
947                         N[s] = 1;
948                     } else {
949                         N[s] = 0;
950                     }
951                     for (i = 2; i <= k; i++) {
952                         for (j = s - i + 1; j <= s; j++) {
953                             if (j <= s - i + 1 || j < 0) {
954                                 a = 0.0;
955                             } else {
956                                 a = N[j];
957                             }
958                             if (j >= s) {
959                                 b = 0.0;
960                             } else {
961                                 b = N[j + 1];
962                             }
963                             den = kn[j + i - 1] - kn[j];
964                             if (den == 0) {
965                                 N[j] = 0;
966                             } else {
967                                 N[j] = (t - kn[j]) / den * a;
968                             }
969                             den = kn[j + i] - kn[j + 1];
970                             if (den != 0) {
971                                 N[j] += (kn[j + i] - t) / den * b;
972                             }
973                         }
974                     }
975                     return N;
976                 },
977                 makeFct = function(which) {
978                     return function(t, suspendedUpdate) {
979                         var len = points.length,
980                             y, j, s,
981                             n = len - 1,
982                             k = order;
983 
984                         if (n <= 0) return NaN;
985                         if (n + 2 <= k) k = n + 1;
986                         if (t <= 0) return points[0][which]();
987                         if (t >= n - k + 2) return points[n][which]();
988 
989                         knots = _knotVector(n, k);
990                         s = Math.floor(t) + k - 1;
991                         N = _evalBasisFuncs(t, knots, n, k, s);
992 
993                         y = 0.0;
994                         for (j = s - k + 1; j <= s; j++) {
995                             if (j < len && j >= 0) y += points[j][which]() * N[j];
996                         }
997                         return y;
998                     };
999                 };
1000 
1001             return [makeFct('X'),
1002                 makeFct('Y'),
1003                 0,
1004                 function() {
1005                     return points.length - 1;
1006                 }
1007             ];
1008         },
1009 
1010         /**
1011          * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve}
1012          * and {@link JXG.Curve#hasPoint}.
1013          * @param {function} f Function in one variable to be differentiated.
1014          * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a
1015          * method of an object and contains a reference to its parent object via "this".
1016          * @returns {function} Derivative function of a given function f.
1017          */
1018         D: function(f, obj) {
1019             var h = 0.00001,
1020                 h2 = 1.0 / (h * 2.0);
1021 
1022             if (arguments.length == 1 || (arguments.length > 1 && !JXG.exists(arguments[1]))) {
1023                 return function(x, suspendUpdate) {
1024                     return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2;
1025                 };
1026             } else { // set "this" to "obj" in f
1027                 return function(x, suspendUpdate) {
1028                     return (f.apply(obj, [x + h,suspendUpdate]) - f.apply(obj, [x - h,suspendUpdate])) * h2;
1029                 };
1030             }
1031         },
1032 
1033         /**
1034          * Helper function to create curve which displays Riemann sums.
1035          * Compute coordinates for the rectangles showing the Riemann sum.
1036          * @param {function} f Function, whose integral is approximated by the Riemann sum.
1037          * @param {Number} n number of rectangles.
1038          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', or 'trapezodial'.
1039          * @param {Number} start Left border of the approximation interval
1040          * @param {Number} end Right border of the approximation interval
1041          * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This
1042          * 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
1043          * rectangles.
1044          */
1045         riemann: function(f, n, type, start, end) {
1046             var xarr = [],
1047                 yarr = [],
1048                 i, j = 0,
1049                 delta,
1050                 x = start, y,
1051                 x1, y1, delta1, 
1052                 sum = 0;
1053 
1054             n = Math.round(n);
1055 
1056             xarr[j] = x;
1057             yarr[j] = 0.0;
1058 
1059             if (n > 0) {
1060                 delta = (end - start) / n;
1061                 delta1 = delta * 0.01; // for 'lower' and 'upper'
1062 
1063                 for (i = 0; i < n; i++) {
1064                     if (type === 'right') {
1065                         y = f(x + delta);
1066                     } else if (type === 'middle') {
1067                         y = f(x + delta * 0.5);
1068                     } else if ((type === 'left') || (type === 'trapezodial')) {
1069                         y = f(x);
1070                     } else if (type === 'lower') {
1071                         y = f(x);
1072                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1073                             y1 = f(x1);
1074                             if (y1 < y) {
1075                                 y = y1;
1076                             }
1077                         }
1078                     } else if (type === 'upper') {
1079                         y = f(x);
1080                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1081                             y1 = f(x1);
1082                             if (y1 > y) {
1083                                 y = y1;
1084                             }
1085                         }
1086                     } else {  // if (type === 'random')
1087                         y = f(x + delta * Math.random());
1088                     }
1089 
1090                     j++;
1091                     xarr[j] = x;
1092                     yarr[j] = y;
1093                     j++;
1094                     x += delta;
1095                     if (type === 'trapezodial') {
1096                         y = f(x);
1097                     }
1098                     xarr[j] = x;
1099                     yarr[j] = y;
1100                     j++;
1101                     xarr[j] = x;
1102                     yarr[j] = 0.0;
1103                     sum += y*delta;
1104                 }
1105             }
1106             return [xarr,yarr,sum];
1107         },
1108 
1109         /**
1110          * Approximate the integral by Riemann sums.
1111          * Compute the area described by the riemann sum rectangles.
1112          * @deprecated Replaced by JXG.Curve.Value(), @see JXG.Curve#riemannsum
1113          * @param {function} f Function, whose integral is approximated by the Riemann sum.
1114          * @param {Number} n number of rectangles.
1115          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random'. or 'trapezodial'.
1116          * @param {Number} start Left border of the approximation interval
1117          * @param {Number} end Right border of the approximation interval
1118          * @returns {Number} The sum of the areas of the rectangles.
1119          */
1120         riemannsum: function(f, n, type, start, end) {
1121             var sum = .0,
1122                 i, delta,
1123                 x = start, y,
1124                 x1, y1, delta1;
1125 
1126             n = Math.floor(n);
1127             if (n > 0) {
1128                 delta = (end - start) / n;
1129                 delta1 = delta * 0.01; // for 'lower' and 'upper'
1130                 for (i = 0; i < n; i++) {
1131                     if (type === 'right') {
1132                         y = f(x + delta);
1133                     } else if (type === 'middle') {
1134                         y = f(x + delta * 0.5);
1135                     } else if (type === 'trapezodial') {
1136                         y = 0.5 * (f(x + delta) + f(x));
1137                     } else if (type === 'left') {
1138                         y = f(x);
1139                     } else if (type === 'lower') {
1140                         y = f(x);
1141                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1142                             y1 = f(x1);
1143                             if (y1 < y) {
1144                                 y = y1;
1145                             }
1146                         }
1147                     } else if (type === 'upper') {
1148                         y = f(x);
1149                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1150                             y1 = f(x1);
1151                             if (y1 > y) {
1152                                 y = y1;
1153                             }
1154                         }
1155                     } else {  // if (type === 'random')
1156                         y = f(x + delta * Math.random());
1157                     }
1158                     sum += delta * y;
1159                     x += delta;
1160                 }
1161             }
1162             return sum;
1163         },
1164 
1165         /**
1166          * Solve initial value problems numerically using Runge-Kutta-methods.
1167          * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
1168          * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
1169          * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
1170          * <pre>
1171          * {
1172          *     s: <Number>,
1173          *     A: <matrix>,
1174          *     b: <Array>,
1175          *     c: <Array>
1176          * }
1177          * </pre>
1178          * 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
1179          * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array.
1180          * @param {Array} I Interval on which to integrate.
1181          * @param {Number} N Number of evaluation points.
1182          * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
1183          * 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
1184          * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has.
1185          * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
1186          * @example
1187          * // A very simple autonomous system dx(t)/dt = x(t);
1188          * function f(t, x) {
1189          *     return x;
1190          * }
1191          *
1192          * // Solve it with initial value x(0) = 1 on the interval [0, 2]
1193          * // with 20 evaluation points.
1194          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
1195          *
1196          * // Prepare data for plotting the solution of the ode using a curve.
1197          * var dataX = [];
1198          * var dataY = [];
1199          * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
1200          * for(var i=0; i<data.length; i++) {
1201          *     dataX[i] = i*h;
1202          *     dataY[i] = data[i][0];
1203          * }
1204          * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
1205          * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
1206          * <script type="text/javascript">
1207          * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
1208          * function f(t, x) {
1209          *     // we have to copy the value.
1210          *     // return x; would just return the reference.
1211          *     return [x[0]];
1212          * }
1213          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
1214          * var dataX = [];
1215          * var dataY = [];
1216          * var h = 0.1;
1217          * for(var i=0; i<data.length; i++) {
1218          *     dataX[i] = i*h;
1219          *     dataY[i] = data[i][0];
1220          * }
1221          * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
1222          * </script><pre>
1223          */
1224         rungeKutta: function(butcher, x0, I, N, f) {
1225             var x = [],
1226                 y = [],
1227                 h = (I[1] - I[0]) / N,
1228                 t = I[0],
1229                 e, i, j,
1230                 k, l,
1231                 dim = x0.length,
1232                 s,
1233                 result = [],
1234                 r = 0;
1235 
1236             if (JXG.isString(butcher)) {
1237                 butcher = predefinedButcher[butcher] || predefinedButcher.euler;
1238             }
1239             s = butcher.s;
1240 
1241             // don't change x0, so copy it
1242             for (e = 0; e < dim; e++)
1243                 x[e] = x0[e];
1244 
1245             for (i = 0; i < N; i++) {
1246                 // Optimization doesn't work for ODEs plotted using time
1247                 //        if((i % quotient == 0) || (i == N-1)) {
1248                 result[r] = [];
1249                 for (e = 0; e < dim; e++)
1250                     result[r][e] = x[e];
1251                 r++;
1252                 //        }
1253                 // init k
1254                 k = [];
1255                 for (j = 0; j < s; j++) {
1256                     // init y = 0
1257                     for (e = 0; e < dim; e++)
1258                         y[e] = 0.;
1259 
1260                     // Calculate linear combination of former k's and save it in y
1261                     for (l = 0; l < j; l++) {
1262                         for (e = 0; e < dim; e++) {
1263                             y[e] += (butcher.A[j][l]) * h * k[l][e];
1264                         }
1265                     }
1266 
1267                     // add x(t) to y
1268                     for (e = 0; e < dim; e++) {
1269                         y[e] += x[e];
1270                     }
1271 
1272                     // calculate new k and add it to the k matrix
1273                     k.push(f(t + butcher.c[j] * h, y));
1274                 }
1275 
1276                 // init y = 0
1277                 for (e = 0; e < dim; e++)
1278                     y[e] = 0.;
1279 
1280                 for (l = 0; l < s; l++) {
1281                     for (e = 0; e < dim; e++)
1282                         y[e] += butcher.b[l] * k[l][e];
1283                 }
1284 
1285                 for (e = 0; e < dim; e++) {
1286                     x[e] = x[e] + h * y[e];
1287                 }
1288 
1289                 t += h;
1290             }
1291 
1292             return result;
1293         },
1294 
1295 
1296         /*
1297          * Maximum number of iterations in @see #fzero
1298          */
1299         maxIterationsRoot: 80, 
1300 
1301         /*
1302          * Maximum number of iterations in @see #fminbr
1303          */
1304         maxIterationsMinimize: 500, 
1305         
1306         /**
1307          *
1308          * Find zero of an univariate function f.
1309          * @param {function} f Function, whose root is to be found
1310          * @param {array or number} x0  Start value or start interval enclosing the root
1311          * @param {object} object Parent object in case f is method of it
1312          * @return {number} the approximation of the root
1313          * Algorithm:
1314          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
1315          *  computations. M., Mir, 1980, p.180 of the Russian edition
1316          * 
1317          * if x0 is an array containing lower and upper bound for the zero
1318          * algorithm 748 is applied. Otherwise, if x0 is a number,
1319          * the algorithm tries to bracket a zero of f starting from x0.
1320          * If this fails, we fall back to Newton's method.
1321          *
1322          **/
1323         fzero: function(f, x0, object) {
1324             var tol = JXG.Math.eps,
1325                 maxiter = this.maxIterationsRoot, niter = 0,
1326                 nfev = 0,
1327                 eps = tol,
1328                 a,b,c, 
1329                 fa, fb, fc,
1330                 aa, blist, i, len, u, fu, 
1331                 prev_step,
1332                 tol_act,         // Actual tolerance
1333                 p,               // Interpolation step is calcu-
1334                 q,               // lated in the form p/q; divi- 
1335                                  // sion operations is delayed  
1336                                  // until the last moment   
1337                 new_step,        // Step at this iteration
1338                 t1, cb, t2;
1339 
1340             if (JXG.isArray(x0)) {
1341                 if (x0.length<2) 
1342                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
1343                 a = x0[0]; fa = f.apply(object,[a]); nfev++;
1344                 b = x0[1]; fb = f.apply(object,[b]); nfev++;
1345             } else {
1346                 a = x0; fa = f.apply(object,[a]); nfev++;
1347                 // Try to get b.
1348                 if (a == 0) {
1349                     aa = 1;
1350                 } else {
1351                     aa = a;
1352                 }
1353                 blist = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa, 1.5*aa, -aa, 2*aa, -10*aa, 10*aa];
1354                 len = blist.length;
1355                 for (i=0;i<len;i++) {
1356                     b = blist[i];
1357                     fb = f.apply(object,[b]); nfev++;
1358                     if (fa*fb<=0) {
1359                         break;
1360                     }
1361                 }
1362                 if (b < a) {
1363                     u = a; a = b; b = u;
1364                     fu = fa; fa = fb; fb = fu;
1365                 }
1366             }
1367 
1368             if (fa*fb > 0) {
1369                 // Bracketing not successful, fall back to Newton's method or to fminbr
1370                 if (JXG.isArray(x0)) {
1371                     //JXG.debug("fzero falls back to fminbr");
1372                     return this.fminbr(f, [a,b], object);
1373                 } else {
1374                     //JXG.debug("fzero falls back to Newton");
1375                     return this.Newton(f, a, object);
1376                 }
1377             }
1378 
1379             // OK, we have enclosed a zero of f.
1380             // Now we can start Brent's method
1381 
1382             c = a;   
1383             fc = fa;
1384             while (niter<maxiter) {   // Main iteration loop
1385                 prev_step = b-a;      // Distance from the last but one
1386                                       // to the last approximation 
1387 
1388                 if ( Math.abs(fc) < Math.abs(fb) ) {
1389                                                 // Swap data for b to be the  
1390                     a = b;  b = c;  c = a;      // best approximation   
1391                     fa=fb;  fb=fc;  fc=fa;
1392                 }
1393                 tol_act = 2*eps*Math.abs(b) + tol*0.5;
1394                 new_step = (c-b)*0.5;
1395 
1396                 if ( Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps ) {
1397                     //JXG.debug("nfev="+nfev);
1398                     return b;                           //  Acceptable approx. is found 
1399                 }
1400                     
1401                 // Decide if the interpolation can be tried 
1402                 if ( Math.abs(prev_step) >= tol_act     // If prev_step was large enough
1403                     && Math.abs(fa) > Math.abs(fb) ) {  // and was in true direction,  
1404                                                         // Interpolatiom may be tried  
1405                     cb = c-b;
1406                     if ( a==c ) {                       // If we have only two distinct 
1407                                                         // points linear interpolation 
1408                         t1 = fb/fa;                     // can only be applied    
1409                         p = cb*t1;
1410                         q = 1.0 - t1;
1411                     } else {                            // Quadric inverse interpolation
1412                         q = fa/fc;  t1 = fb/fc;  t2 = fb/fa;
1413                         p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
1414                         q = (q-1.0) * (t1-1.0) * (t2-1.0);
1415                     }
1416                     if ( p>0 ) {                        // p was calculated with the op-
1417                         q = -q;                         // posite sign; make p positive 
1418                     } else {                            // and assign possible minus to
1419                         p = -p;                         // q 
1420                     }
1421                     
1422                     if( p < (0.75*cb*q-Math.abs(tol_act*q)*0.5)     // If b+p/q falls in [b,c]
1423                         && p < Math.abs(prev_step*q*0.5) ) {        // and isn't too large 
1424                         new_step = p/q;                 // it is accepted   
1425                     }                                   // If p/q is too large then the 
1426                                                         // bissection procedure can
1427                                                         // reduce [b,c] range to more
1428                                                         // extent
1429                 }
1430 
1431                 if ( Math.abs(new_step) < tol_act ) {   // Adjust the step to be not less
1432                     if ( new_step > 0 ) {               // than tolerance
1433                         new_step = tol_act;
1434                     } else {
1435                         new_step = -tol_act;
1436                     }
1437                 }
1438 
1439                 a = b;  fa = fb;                        // Save the previous approx.
1440                 b += new_step;  
1441                 fb = f.apply(object,[b]); nfev++;       // Do step to a new approxim.
1442                 if ( (fb>0 && fc>0) || (fb<0 && fc<0) ) {
1443                                                         // Adjust c for it to have a sign
1444                     c = a;  fc = fa;                    // opposite to that of b 
1445                 }
1446                 niter++;
1447             }                                           // End while
1448             
1449             //JXG.debug("fzero: maxiter="+maxiter+" reached.");
1450             return b;
1451         },
1452 
1453      
1454         /**
1455          *
1456          * Find minimum of an univariate function f.
1457          * @param {function} f Function, whose minimum is to be found
1458          * @param {array} x0  Start interval enclosing the minimum
1459          * @param {object} object Parent object in case f is method of it
1460          * @return {number} the approximation of the minimum
1461          * Algorithm:
1462          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
1463          *  computations. M., Mir, 1980, p.180 of the Russian edition
1464          * x0 
1465          **/
1466         fminbr: function(f, x0, object) {              // An estimate to the min location
1467             var a, b, x, v, w,
1468                 fx, fv, fw,
1469                 r = (3.-Math.sqrt(5.0))*0.5,            // Golden section ratio   
1470                 tol = JXG.Math.eps,
1471                 sqrteps = Math.sqrt(JXG.Math.eps),
1472                 maxiter = this.maxIterationsMinimize, 
1473                 niter = 0,
1474                 range, middle_range, tol_act, new_step,
1475                 p, q, t, ft,
1476                 nfev = 0;
1477 
1478             if (!JXG.isArray(x0) || x0.length<2) {
1479                 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two.");
1480             }
1481             a = x0[0];
1482             b = x0[1];
1483             v = a + r*(b-a);  
1484             fv = f.apply(object,[v]); nfev++;           // First step - always gold section
1485             x = v;  w = v;
1486             fx=fv;  fw=fv;
1487 
1488             while (niter<maxiter) {
1489                 range = b-a;                            // Range over which the minimum 
1490                                                         // is seeked for
1491                 middle_range = (a+b)*0.5;
1492                 tol_act = sqrteps*Math.abs(x) + tol/3;  // Actual tolerance  
1493                 if( Math.abs(x-middle_range) + range*0.5 <= 2*tol_act ) {
1494                     //JXG.debug(nfev);
1495                     return x;                           // Acceptable approx. is found
1496                 }
1497                                                         // Obtain the golden section step 
1498                 new_step = r * ( x<middle_range ? b-x : a-x );
1499                                                         // Decide if the interpolation can be tried 
1500                 if ( Math.abs(x-w) >= tol_act  ) {      // If x and w are distinct     
1501                                                         // interpolatiom may be tried 
1502                     // Interpolation step is calculated as p/q; 
1503                     // division operation is delayed until last moment 
1504                     t = (x-w) * (fx-fv);
1505                     q = (x-v) * (fx-fw);
1506                     p = (x-v)*q - (x-w)*t;
1507                     q = 2*(q-t);
1508 
1509                     if ( q>0 ) {                        // q was calculated with the op-
1510                         p = -p;                         // posite sign; make q positive 
1511                     } else {                            // and assign possible minus to
1512                         q = -q;                         // p
1513                     }
1514                     if ( Math.abs(p) < Math.abs(new_step*q) &&      // If x+p/q falls in [a,b]
1515                          p > q*(a-x+2*tol_act) &&                   //  not too close to a and
1516                          p < q*(b-x-2*tol_act)  ) {                 // b, and isn't too large */
1517                          new_step = p/q;                            // it is accepted        
1518                     }
1519                     // If p/q is too large then the 
1520                     // golden section procedure can   
1521                     // reduce [a,b] range to more   
1522                     // extent           
1523                 }
1524 
1525                 if ( Math.abs(new_step) < tol_act ) {    // Adjust the step to be not less
1526                     if( new_step > 0 ) {                 // than tolerance     
1527                         new_step = tol_act;
1528                     } else {
1529                         new_step = -tol_act;
1530                     }
1531                 }
1532                 
1533                 // Obtain the next approximation to min 
1534                 // and reduce the enveloping range
1535                 t = x + new_step;                       // Tentative point for the min
1536                 ft = f.apply(object,[t]); nfev++;
1537                 if ( ft <= fx ) {                       // t is a better approximation 
1538                     if ( t < x ) {                      // Reduce the range so that
1539                         b = x;                          // t would fall within it 
1540                     } else {
1541                         a = x;
1542                     }
1543                     v = w;  w = x;  x = t;              // Assign the best approx to x 
1544                     fv=fw;  fw=fx;  fx=ft;
1545                 } else {                                // x remains the better approx
1546                     if ( t < x ) {                      // Reduce the range enclosing x 
1547                         a = t;                   
1548                     } else {
1549                         b = t;
1550                     }
1551                     if ( ft <= fw || w==x ) {
1552                         v = w;  w = t;
1553                         fv=fw;  fw=ft;
1554                     } else if ( ft<=fv || v==x || v==w ) {
1555                         v = t;
1556                         fv=ft;
1557                     }
1558                 }
1559                 niter++;
1560             } 
1561             //JXG.debug("fminbr: maxiter="+maxiter+" reached.");
1562             return x;
1563         },
1564 
1565         /**
1566          * Helper function to create curve which displays Reuleaux polygons.
1567          * @param {array} points Array of points which should be the vertices of the Reuleaux polygon. Typically,
1568          *                       these point list is the array vrtices of a regular polygon.
1569          * @param {number} nr Number of vertices
1570          * @returns {array} An array containing the two functions defining the Reuleaux polygon and the two values
1571          * for the start and the end of the paramtric curve.
1572          * array may be used as parent array of a {@link JXG.Curve}.
1573          *
1574          * @example
1575          * var A = brd.create('point',[-2,-2]);
1576          * var B = brd.create('point',[0,1]);
1577          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 
1578          * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 
1579          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
1580          *
1581          * </pre><div id="2543a843-46a9-4372-abc1-94d9ad2db7ac" style="width: 300px; height: 300px;"></div>
1582          * <script type="text/javascript">
1583          * var brd = JXG.JSXGraph.initBoard('2543a843-46a9-4372-abc1-94d9ad2db7ac', {boundingbox: [-5, 5, 5, -5], axis: true, showcopyright:false, shownavigation: false});
1584          * var A = brd.create('point',[-2,-2]);
1585          * var B = brd.create('point',[0,1]);
1586          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 
1587          * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 
1588          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
1589          * </script><pre>
1590          */
1591         reuleauxPolygon: function(points, nr) {
1592             var pi2 = Math.PI*2,
1593                 pi2_n = pi2/nr,
1594                 diag = (nr-1)/2,
1595                 beta, d = 0,
1596                 makeFct = function(which, trig) {
1597                     return function(t, suspendUpdate) {
1598                         if (!suspendUpdate) {
1599                             d = points[0].Dist(points[diag]);
1600                             beta = JXG.Math.Geometry.rad([points[0].X()+1,points[0].Y()],points[0],points[(diag)%nr]);
1601                         }
1602                         var t1 = (t%pi2 + pi2) % pi2;
1603                         var j = Math.floor(t1 / pi2_n)%nr;
1604                         if (isNaN(j)) return j;
1605                         //t1 = (t1-j*pi2_n)*0.5 + beta+j*pi2_n;
1606                         t1 = t1*0.5+j*pi2_n*0.5 + beta;
1607                         return points[j][which]()+d*Math[trig](t1);
1608                     };
1609                 };
1610             return [
1611                 makeFct('X','cos'),
1612                 makeFct('Y','sin'),
1613                 0,
1614                 Math.PI*2
1615             ];
1616         },
1617 
1618         /**
1619          * Implements the Ramer-Douglas-Peuker algorithm.
1620          * It discards points which are not necessary from the polygonal line defined by the point array
1621          * pts. The computation is done in screen coordinates.
1622          * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points.
1623          * @param {Array} pts Array of {@link JXG.Coords}
1624          * @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>.
1625          * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points.
1626          */
1627         RamerDouglasPeuker: function(pts, eps) {
1628             var newPts = [], i, k, len,
1629                 /**
1630                  * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}.
1631                  * It runs recursively through the point set and searches the
1632                  * point which has the largest distance from the line between the first point and
1633                  * the last point. If the distance from the line is greater than eps, this point is
1634                  * included in our new point set otherwise it is discarded.
1635                  * If it is taken, we recursively apply the subroutine to the point set before
1636                  * and after the chosen point.
1637                  * @param {Array} pts Array of {@link JXG.Coords}
1638                  * @param {Number} i Index of an element of pts
1639                  * @param {Number} j Index of an element of pts
1640                  * @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>.
1641                  * @param {Array} newPts Array of {@link JXG.Coords}
1642                  * @private
1643                  */
1644                 RDP = function(pts, i, j, eps, newPts) {
1645                     var result = findSplit(pts, i, j);
1646 
1647                     if (result[0] > eps) {
1648                         RDP(pts, i, result[1], eps, newPts);
1649                         RDP(pts, result[1], j, eps, newPts);
1650                     } else {
1651                         newPts.push(pts[j]);
1652                     }
1653                 },
1654                 /**
1655                  * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}.
1656                  * It searches for the point between index i and j which
1657                  * has the largest distance from the line between the points i and j.
1658                  * @param {Array} pts Array of {@link JXG.Coords}
1659                  * @param {Number} i Index of a point in pts
1660                  * @param {Number} j Index of a point in pts
1661                  **/
1662                 findSplit = function(pts, i, j) {
1663                     var dist = 0,
1664                         f = i,
1665                         d, k, ci, cj, ck,
1666                         x0, y0, x1, y1,
1667                         den, lbda;
1668 
1669                     if (j - i < 2) return [-1.0,0];
1670 
1671                     ci = pts[i].scrCoords;
1672                     cj = pts[j].scrCoords;
1673                     if (isNaN(ci[1] + ci[2] + cj[1] + cj[2])) return [NaN,j];
1674 
1675                     for (k = i + 1; k < j; k++) {
1676                         ck = pts[k].scrCoords;
1677                         x0 = ck[1] - ci[1];
1678                         y0 = ck[2] - ci[2];
1679                         x1 = cj[1] - ci[1];
1680                         y1 = cj[2] - ci[2];
1681                         den = x1 * x1 + y1 * y1;
1682                         /*
1683                         if (den >= JXG.Math.eps) {
1684                             lbda = (x0 * x1 + y0 * y1) / den;
1685                             d = x0 * x0 + y0 * y0 - lbda * (x0 * x1 + y0 * y1);
1686                         } else {
1687                             lbda = 0.0;
1688                             d = x0 * x0 + y0 * y0;
1689                         }
1690                         if (lbda < 0.0) {
1691                             d = x0 * x0 + y0 * y0;
1692                         } else if (lbda > 1.0) {
1693                             x0 = ck[1] - cj[1];
1694                             y0 = ck[2] - cj[2];
1695                             d = x0 * x0 + y0 * y0;
1696                         }
1697                         */
1698                         if (den >= JXG.Math.eps) {
1699                             lbda = (x0 * x1 + y0 * y1) / den;
1700                             //d = x0 * x0 + y0 * y0 - lbda * (x0 * x1 + y0 * y1);
1701                             if (lbda<0.0) {
1702                                 lbda = 0.0;
1703                             } else if (lbda>1.0) {
1704                                 lbda = 1.0;
1705                             }
1706                             x0 = x0-lbda*x1;
1707                             y0 = y0-lbda*y1;
1708                             d = x0*x0+y0*y0;
1709                         } else {
1710                             lbda = 0.0;
1711                             d = x0 * x0 + y0 * y0;
1712                         }
1713                         
1714                         if (d > dist) {
1715                             dist = d;
1716                             f = k;
1717                         }
1718                     }
1719                     return [Math.sqrt(dist),f];
1720                 };
1721 
1722             len = pts.length;
1723 
1724             // Search for the left most point woithout NaN coordinates
1725             i = 0;
1726             while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) {
1727                 i++;
1728             }
1729             // Search for the right most point woithout NaN coordinates
1730             k = len - 1;
1731             while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) {
1732                 k--;
1733             }
1734 
1735             // Only proceed if something is left
1736             if (!(i > k || i == len)) {
1737                 newPts[0] = pts[i];
1738                 RDP(pts, i, k, eps, newPts);
1739             }
1740 
1741             return newPts;
1742         }
1743 
1744     }
1745 })(JXG, Math);
1746 
1747