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             return this.gaussBareiss(mat);
232         },
233         
234         /**
235          * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method
236          * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990
237          * @param {Array} Ain A symmetric 3x3 matrix.
238          * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors.
239          */
240         Jacobi: function(Ain) {
241             var i,j,k,aa,si,co,tt,eps = JXG.Math.eps,
242                 sum = 0.0,
243                 ssum, amax,
244                 n = Ain.length,
245                 V = [
246                     [0,0,0],
247                     [0,0,0],
248                     [0,0,0]
249                 ],
250                 A = [
251                     [0,0,0],
252                     [0,0,0],
253                     [0,0,0]
254                 ], nloops=0;
255 
256             // Initialization. Set initial Eigenvectors.
257             for (i = 0; i < n; i++) {
258                 for (j = 0; j < n; j++) {
259                     V[i][j] = 0.0;
260                     A[i][j] = Ain[i][j];
261                     sum += Math.abs(A[i][j]);
262                 }
263                 V[i][i] = 1.0;
264             }
265             // Trivial problems
266             if (n == 1) {
267                 return [A,V];
268             }
269             if (sum <= 0.0) {
270                 return [A,V];
271             }
272 
273             sum /= (n * n);
274 
275             // Reduce matrix to diagonal
276             do {
277                 ssum = 0.0;
278                 amax = 0.0;
279                 for (j = 1; j < n; j++) {
280                     for (i = 0; i < j; i++) {
281                         // Check if A[i][j] is to be reduced
282                         aa = Math.abs(A[i][j]);
283                         if (aa > amax) {
284                             amax = aa;
285                         }
286                         ssum += aa;
287                         if (aa >= eps) {
288                             // calculate rotation angle
289                             aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5;
290                             si = Math.sin(aa);
291                             co = Math.cos(aa);
292                             // Modify 'i' and 'j' columns
293                             for (k = 0; k < n; k++) {
294                                 tt = A[k][i];
295                                 A[k][i] = co * tt + si * A[k][j];
296                                 A[k][j] = -si * tt + co * A[k][j];
297                                 tt = V[k][i];
298                                 V[k][i] = co * tt + si * V[k][j];
299                                 V[k][j] = -si * tt + co * V[k][j];
300                             }
301                             // Modify diagonal terms
302                             A[i][i] = co * A[i][i] + si * A[j][i];
303                             A[j][j] = -si * A[i][j] + co * A[j][j];
304                             A[i][j] = 0.0;
305                             // Make 'A' matrix symmetrical
306                             for (k = 0; k < n; k++) {
307                                 A[i][k] = A[k][i];
308                                 A[j][k] = A[k][j];
309                             }
310                             // A[i][j] made zero by rotation
311                         }
312                     }
313                 }
314                 nloops++;
315             } while (Math.abs(ssum) / sum > eps && nloops<2000);
316             return [A,V];
317         },
318 
319         /**
320          * Calculates the integral of function f over interval using Newton-Cotes-algorithm.
321          * @param {Array} interval The integration interval, e.g. [0, 3].
322          * @param {function} f A function which takes one argument of type number and returns a number.
323          * @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
324          * with value being either 'trapez', 'simpson', or 'milne'.
325          * @returns {Number} Integral value of f over interval
326          * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use
327          * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4.
328          * @example
329          * function f(x) {
330          *   return x*x;
331          * }
332          *
333          * // calculates integral of <tt>f</tt> from 0 to 2.
334          * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f);
335          *
336          * // the same with an anonymous function
337          * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; });
338          *
339          * // use trapez rule with 16 nodes
340          * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f,
341          *                                   {number_of_nodes: 16, intergration_type: 'trapez'});
342          */
343         NewtonCotes: function(interval, f, config) {
344             var integral_value = 0.0,
345                 number_of_nodes = config && typeof config.number_of_nodes === 'number' ? config.number_of_nodes : 28,
346                 available_types = {trapez: true, simpson: true, milne: true},
347                 integration_type = config && config.integration_type && available_types.hasOwnProperty(config.integration_type) && available_types[config.integration_type] ? config.integration_type : 'milne',
348                 step_size = (interval[1] - interval[0]) / number_of_nodes,
349                 evaluation_point, i, number_of_intervals;
350 
351             switch (integration_type) {
352                 case 'trapez':
353                     integral_value = (f(interval[0]) + f(interval[1])) * 0.5;
354 
355                     evaluation_point = interval[0];
356                     for (i = 0; i < number_of_nodes - 1; i++) {
357                         evaluation_point += step_size;
358                         integral_value += f(evaluation_point);
359                     }
360                     integral_value *= step_size;
361 
362                     break;
363                 case 'simpson':
364                     if (number_of_nodes % 2 > 0) {
365                         throw new Error("JSXGraph:  INT_SIMPSON requires config.number_of_nodes dividable by 2.");
366                     }
367                     number_of_intervals = number_of_nodes / 2.0;
368                     integral_value = f(interval[0]) + f(interval[1]);
369                     evaluation_point = interval[0];
370                     for (i = 0; i < number_of_intervals - 1; i++) {
371                         evaluation_point += 2.0 * step_size;
372                         integral_value += 2.0 * f(evaluation_point);
373                     }
374                     evaluation_point = interval[0] - step_size;
375                     for (i = 0; i < number_of_intervals; i++) {
376                         evaluation_point += 2.0 * step_size;
377                         integral_value += 4.0 * f(evaluation_point);
378                     }
379                     integral_value *= step_size / 3.0;
380                     break;
381                 default:
382                     if (number_of_nodes % 4 > 0) {
383                         throw new Error("JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4");
384                     }
385                     number_of_intervals = number_of_nodes * 0.25;
386                     integral_value = 7.0 * (f(interval[0]) + f(interval[1]));
387                     evaluation_point = interval[0];
388                     for (i = 0; i < number_of_intervals - 1; i++) {
389                         evaluation_point += 4.0 * step_size;
390                         integral_value += 14.0 * f(evaluation_point);
391                     }
392                     evaluation_point = interval[0] - 3.0 * step_size;
393                     for (i = 0; i < number_of_intervals; i++) {
394                         evaluation_point += 4.0 * step_size;
395                         integral_value += 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size));
396                     }
397                     evaluation_point = interval[0] - 2.0 * step_size;
398                     for (i = 0; i < number_of_intervals; i++) {
399                         evaluation_point += 4.0 * step_size;
400                         integral_value += 12.0 * f(evaluation_point);
401                     }
402                     integral_value *= 2.0 * step_size / 45.0;
403             }
404             return integral_value;
405         },
406 
407         /**
408          * Integral of function f over interval.
409          * @param {Array} interval The integration interval, e.g. [0, 3].
410          * @param {function} f A function which takes one argument of type number and returns a number.
411          * @returns {Number} The value of the integral of f over interval
412          * @see JXG.Math.Numerics.NewtonCotes
413          */
414         I: function(interval, f) {
415             return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'});
416         },
417 
418         /**
419          * Newton's method to find roots of a funtion in one variable.
420          * @param {function} f We search for a solution of f(x)=0.
421          * @param {Number} x initial guess for the root, i.e. start value.
422          * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a
423          *                 method of an object and contains a reference to its parent object via "this".
424          * @returns {Number} A root of the function f.
425          */
426         Newton: function(f, x, object) {
427             var i = 0,
428                 h = JXG.Math.eps,
429                 newf = f.apply(object, [x]), // set "this" to "object" in f
430                 nfev = 1, 
431                 df;
432             if (JXG.isArray(x)) {  // For compatibility
433                 x = x[0];
434             }
435             while (i < 50 && Math.abs(newf) > h) {
436                 df = this.D(f, object)(x);  nfev += 2;
437                 if (Math.abs(df) > h) {
438                     x -= newf / df;
439                 } else {
440                     x += (Math.random() * 0.2 - 1.0);
441                 }
442                 newf = f.apply(object, [x]); nfev++;
443                 i++;
444             }
445             //JXG.debug("N nfev="+nfev);
446             return x;
447         },
448 
449         /**
450          * Abstract method to find roots of univariate functions.
451          * @param {function} f We search for a solution of f(x)=0.
452          * @param {Number} x initial guess for the root, i.e. staring value.
453          * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a
454          *                 method of an object and contains a reference to its parent object via "this".
455          * @returns {Number} A root of the function f.
456          */
457         root: function(f, x, object) {
458             //return this.Newton(f, x, object);
459             return this.fzero(f, x, object);
460         },
461 
462         /**
463          * Returns the Lagrange polynomials for curves with equidistant nodes, see
464          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
465          * SIAM Review, Vol 46, No 3, (2004) 501-517.
466          * The graph of the parametric curve [x(t),y(t)] runs through the given points.
467          * @param {Array} p Array of JXG.Points
468          * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve
469          * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero.
470          */
471         Neville: function(p) {
472             var w = [],
473                 makeFct = function(fun) {
474                     return function(t, suspendedUpdate) {
475                         var i, d, s,
476                             bin = JXG.Math.binomial,
477                             len = p.length,
478                             len1 = len - 1,
479                             num = 0.0,
480                             denom = 0.0;
481 
482                         if (!suspendedUpdate) {
483                             s = 1;
484                             for (i = 0; i < len; i++) {
485                                 w[i] = bin(len1, i) * s;
486                                 s *= (-1);
487                             }
488                         }
489 
490                         d = t;
491                         for (i = 0; i < len; i++) {
492                             if (d === 0) {
493                                 return p[i][fun]();
494                             } else {
495                                 s = w[i] / d;
496                                 d--;
497                                 num += p[i][fun]() * s;
498                                 denom += s;
499                             }
500                         }
501                         return num / denom;
502                     };
503                 },
504 
505                 xfct = makeFct('X'),
506                 yfct = makeFct('Y');
507 
508             return [xfct, yfct, 0, function() {
509                 return p.length - 1;
510             }];
511         },
512 
513         /**
514          * Calculates second derivatives at the given knots.
515          * @param {Array} x x values of knots
516          * @param {Array} y y values of knots
517          * @returns {Array} Second derivatives of the interpolated function at the knots.
518          * @see #splineEval
519          */
520         splineDef: function(x, y) {
521             var n = Math.min(x.length, y.length),
522                 pair, i, l,
523                 diag = [],
524                 z = [],
525                 data = [],
526                 dx = [],
527                 delta = [],
528                 F = [];
529 
530             if (n === 2) {
531                 return [0, 0];
532             }
533 
534             for (i = 0; i < n; i++) {
535                 pair = {X: x[i], Y: y[i]};
536                 data.push(pair);
537             }
538             data.sort(function (a, b) {
539                 return a.X - b.X;
540             });
541             for (i = 0; i < n; i++) {
542                 x[i] = data[i].X;
543                 y[i] = data[i].Y;
544             }
545 
546             for (i = 0; i < n - 1; i++) {
547                 dx.push(x[i + 1] - x[i]);
548             }
549             for (i = 0; i < n - 2; i++) {
550                 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i]));
551             }
552 
553             // ForwardSolve
554             diag.push(2 * (dx[0] + dx[1]));
555             z.push(delta[0]);
556 
557             for (i = 0; i < n - 3; i++) {
558                 l = dx[i + 1] / diag[i];
559                 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]);
560                 z.push(delta[i + 1] - l * z[i]);
561             }
562 
563             // BackwardSolve
564             F[n - 3] = z[n - 3] / diag[n - 3];
565             for (i = n - 4; i >= 0; i--) {
566                 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i];
567             }
568 
569             // Generate f''-Vector
570             for (i = n - 3; i >= 0; i--) {
571                 F[i + 1] = F[i];
572             }
573 
574             // natural cubic spline
575             F[0] = 0;
576             F[n - 1] = 0;
577             return F;
578         },
579 
580         /**
581          * Evaluate points on spline.
582          * @param {Number,Array} x0 A single float value or an array of values to evaluate
583          * @param {Array} x x values of knots
584          * @param {Array} y y values of knots
585          * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef}
586          * @see #splineDef
587          * @returns {Number,Array} A single value or an array, depending on what is given as x0.
588          */
589         splineEval: function(x0, x, y, F) {
590             var n = Math.min(x.length, y.length),
591                 l = 1,
592                 asArray = false,
593                 y0 = [],
594                 i, j, a, b, c, d, x_;
595 
596             // number of points to be evaluated
597             if (JXG.isArray(x0)) {
598                 l = x0.length;
599                 asArray = true;
600             } else
601                 x0 = [x0];
602 
603             for (i = 0; i < l; i++) {
604                 // is x0 in defining interval?
605                 if ((x0[i] < x[0]) || (x[i] > x[n - 1]))
606                     return NaN;
607 
608                 // determine part of spline in which x0 lies
609                 for (j = 1; j < n; j++) {
610                     if (x0[i] <= x[j])
611                         break;
612                 }
613                 j--;
614 
615                 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1];
616                 // determine the coefficients of the polynomial in this interval
617                 a = y[j];
618                 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]);
619                 c = F[j] / 2;
620                 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j]));
621                 // evaluate x0[i]
622                 x_ = x0[i] - x[j];
623                 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_);
624                 y0.push(a + (b + (c + d * x_) * x_) * x_);
625             }
626 
627             if (asArray)
628                 return y0;
629             else
630                 return y0[0];
631         },
632 
633         /**
634          * Generate a string containing the function term of a polynomial.
635          * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i.
636          * @param {Number} deg Degree of the polynomial
637          * @param {String} varname Name of the variable (usually 'x')
638          * @param {Number} prec Precision
639          * @returns {String} A string containg the function term of the polynomial.
640          */
641         generatePolynomialTerm: function(coeffs, deg, varname, prec) {
642             var t = [], i;
643             for (i = deg; i >= 0; i--) {
644                 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']);
645                 if (i > 1) {
646                     t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']);
647                 } else if (i === 1) {
648                     t = t.concat(['*', varname, ' + ']);
649                 }
650             }
651             return t.join('');
652         },
653 
654         /**
655          * Computes the polynomial through a given set of coordinates in Lagrange form.
656          * Returns the Lagrange polynomials, see
657          * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation,
658          * SIAM Review, Vol 46, No 3, (2004) 501-517.
659          * @param {Array} p Array of JXG.Points
660          * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points.
661          */
662         lagrangePolynomial: function(p) {
663             var w = [],
664                 fct = function(x, suspendedUpdate) {
665                     var i, k, len, xi, s,
666                         num = 0, denom = 0,
667                         M, j;
668 
669                     len = p.length;
670                     if (!suspendedUpdate) {
671                         for (i = 0; i < len; i++) {
672                             w[i] = 1.0;
673                             xi = p[i].X();
674                             for (k = 0; k < len; k++) if (k != i) {
675                                 w[i] *= (xi - p[k].X());
676                             }
677                             w[i] = 1 / w[i];
678                         }
679 
680                         M = [];
681                         for (j = 0; j < len; j++) {
682                             M.push([1]);
683                         }
684                     }
685 
686                     for (i = 0; i < len; i++) {
687                         xi = p[i].X();
688                         if (x === xi) {
689                             return p[i].Y();
690                         } else {
691                             s = w[i] / (x - xi);
692                             denom += s;
693                             num += s * p[i].Y();
694                         }
695                     }
696                     return num / denom;
697                 };
698 
699             fct.getTerm = function() {
700                 return '';
701             };
702 
703             return fct;
704         },
705 
706         /**
707          * Computes the regression polynomial of a given degree through a given set of coordinates.
708          * Returns the regression polynomial function.
709          * @param {Number} degree number, function or slider.
710          * Either
711          * @param {Array} dataX array containing the x-coordinates of the data set
712          * @param {Array} dataY array containing the y-coordinates of the data set,
713          * or
714          * @param {Array} dataX Array consisting of JXG.Points.
715          * @param {} dataY Ignored
716          * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree.
717          * It possesses the method getTerm() which returns the string containing the function term of the polynomial.
718          */
719         regressionPolynomial: function(degree, dataX, dataY) {
720             var coeffs,
721                 deg, dX, dY,
722                 inputType,
723                 fct,
724                 term = '';
725 
726             if (JXG.isPoint(degree) && typeof degree.Value == 'function') {  // Slider
727                 deg = function() {return degree.Value();};
728             } else if (JXG.isFunction(degree)) {
729                 deg = degree;
730             } else if (JXG.isNumber(degree)) {
731                 deg = function() {return degree;};
732             } else {
733                 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'.");
734             }
735 
736             if (arguments.length == 3 && JXG.isArray(dataX) && JXG.isArray(dataY)) {              // Parameters degree, dataX, dataY
737                 inputType = 0;
738             } else if (arguments.length == 2 && JXG.isArray(dataX) && dataX.length > 0 && JXG.isPoint(dataX[0])) {  // Parameters degree, point array
739                 inputType = 1;
740             } else {
741                 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters.");
742             }
743 
744             fct = function(x, suspendedUpdate) {
745                 var i, j, M, MT, y, B, c, s,
746                     d, len = dataX.length;                        // input data
747 
748                 d = Math.floor(deg());                      // input data
749                 if (!suspendedUpdate) {
750                     if (inputType === 1) {  // point list as input
751                         dX = [];
752                         dY = [];
753                         for (i = 0; i < len; i++) {
754                             dX[i] = dataX[i].X();
755                             dY[i] = dataX[i].Y();
756                         }
757                     }
758 
759                     if (inputType === 0) {  // check for functions
760                         dX = [];
761                         dY = [];
762                         for (i = 0; i < len; i++) {
763                             if (JXG.isFunction(dataX[i]))
764                                 dX.push(dataX[i]());
765                             else
766                                 dX.push(dataX[i]);
767                             if (JXG.isFunction(dataY[i]))
768                                 dY.push(dataY[i]());
769                             else
770                                 dY.push(dataY[i]);
771                         }
772                     }
773 
774                     M = [];
775                     for (j = 0; j < len; j++) {
776                         M.push([1]);
777                     }
778                     for (i = 1; i <= d; i++) {
779                         for (j = 0; j < len; j++) {
780                             M[j][i] = M[j][i - 1] * dX[j];      // input data
781                         }
782                     }
783 
784                     y = dY;                                 // input data
785                     MT = JXG.Math.transpose(M);
786                     B = JXG.Math.matMatMult(MT, M);
787                     c = JXG.Math.matVecMult(MT, y);
788                     coeffs = JXG.Math.Numerics.Gauss(B, c);
789                     term = JXG.Math.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3);
790                 }
791 
792                 // Horner's scheme to evaluate polynomial
793                 s = coeffs[d];
794                 for (i = d - 1; i >= 0; i--) {
795                     s = (s * x + coeffs[i]);
796                 }
797                 return s;
798             };
799 
800             fct.getTerm = function() {
801                 return term;
802             };
803 
804             return fct;
805         },
806 
807         /**
808          * Computes the cubic Bezier curve through a given set of points.
809          * @param {Array} points Array consisting of 3*k+1 JXG.Points.
810          * The points at position k with k mod 3 = 0 are the data points,
811          * points at position k with k mod 3 = 1 or 2 are the control points.
812          * @returns {Array} An array consisting of two functions of one parameter t which return the
813          * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting
814          * no parameters and returning one third of the length of the points. 
815          */
816         bezier: function(points) {
817             var len,
818                 makeFct = function(which) {
819                     return function(t, suspendedUpdate) {
820                         var z = Math.floor(t) * 3,
821                             t0 = t % 1,
822                             t1 = 1 - t0;
823 
824                         if (!suspendedUpdate) {
825                             len = Math.floor(points.length / 3);
826                         }
827 
828                         if (t < 0) {
829                             return points[0][which]();
830                         }
831                         if (t >= len) {
832                             return points[points.length - 1][which]();
833                         }
834                         if (isNaN(t)) {
835                             return NaN;
836                         }
837                         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;
838                     };
839                 };
840 
841             return [
842                 makeFct('X'),
843                 makeFct('Y'),
844                 0,
845                 function() {
846                     return Math.floor(points.length / 3);
847                 }];
848         },
849 
850         /**
851          * Computes the B-spline curve of order k (order = degree+1) through a given set of points.
852          * @param {Array} points Array consisting of JXG.Points.
853          * @param {Number} order Order of the B-spline curve.
854          * @todo closed B-spline curves
855          * @returns {Array} An Array consisting of four components: Two functions each of one parameter t
856          * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply
857          * returning the length of the points array minus one.
858          */
859         bspline: function(points, order) {
860             var knots, N = [],
861                 _knotVector = function(n, k) {
862                     var j, kn = [];
863                     for (j = 0; j < n + k + 1; j++) {
864                         if (j < k) {
865                             kn[j] = 0.0;
866                         } else if (j <= n) {
867                             kn[j] = j - k + 1;
868                         } else {
869                             kn[j] = n - k + 2;
870                         }
871                     }
872                     return kn;
873                 },
874 
875                 _evalBasisFuncs = function(t, kn, n, k, s) {
876                     var i,j,a,b,den,
877                         N = [];
878 
879                     if (kn[s] <= t && t < kn[s + 1]) {
880                         N[s] = 1;
881                     } else {
882                         N[s] = 0;
883                     }
884                     for (i = 2; i <= k; i++) {
885                         for (j = s - i + 1; j <= s; j++) {
886                             if (j <= s - i + 1 || j < 0) {
887                                 a = 0.0;
888                             } else {
889                                 a = N[j];
890                             }
891                             if (j >= s) {
892                                 b = 0.0;
893                             } else {
894                                 b = N[j + 1];
895                             }
896                             den = kn[j + i - 1] - kn[j];
897                             if (den == 0) {
898                                 N[j] = 0;
899                             } else {
900                                 N[j] = (t - kn[j]) / den * a;
901                             }
902                             den = kn[j + i] - kn[j + 1];
903                             if (den != 0) {
904                                 N[j] += (kn[j + i] - t) / den * b;
905                             }
906                         }
907                     }
908                     return N;
909                 },
910                 makeFct = function(which) {
911                     return function(t, suspendedUpdate) {
912                         var len = points.length,
913                             y, j, s,
914                             n = len - 1,
915                             k = order;
916 
917                         if (n <= 0) return NaN;
918                         if (n + 2 <= k) k = n + 1;
919                         if (t <= 0) return points[0][which]();
920                         if (t >= n - k + 2) return points[n][which]();
921 
922                         knots = _knotVector(n, k);
923                         s = Math.floor(t) + k - 1;
924                         N = _evalBasisFuncs(t, knots, n, k, s);
925 
926                         y = 0.0;
927                         for (j = s - k + 1; j <= s; j++) {
928                             if (j < len && j >= 0) y += points[j][which]() * N[j];
929                         }
930                         return y;
931                     };
932                 };
933 
934             return [makeFct('X'),
935                 makeFct('Y'),
936                 0,
937                 function() {
938                     return points.length - 1;
939                 }
940             ];
941         },
942 
943         /**
944          * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve}
945          * and {@link JXG.Curve#hasPoint}.
946          * @param {function} f Function in one variable to be differentiated.
947          * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a
948          * method of an object and contains a reference to its parent object via "this".
949          * @returns {function} Derivative function of a given function f.
950          */
951         D: function(f, obj) {
952             var h = 0.00001,
953                 h2 = 1.0 / (h * 2.0);
954 
955             if (arguments.length == 1 || (arguments.length > 1 && !JXG.exists(arguments[1]))) {
956                 return function(x, suspendUpdate) {
957                     return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2;
958                 };
959             } else { // set "this" to "obj" in f
960                 return function(x, suspendUpdate) {
961                     return (f.apply(obj, [x + h,suspendUpdate]) - f.apply(obj, [x - h,suspendUpdate])) * h2;
962                 };
963             }
964         },
965 
966         /**
967          * Helper function to create curve which displays Riemann sums.
968          * Compute coordinates for the rectangles showing the Riemann sum.
969          * @param {function} f Function, whose integral is approximated by the Riemann sum.
970          * @param {Number} n number of rectangles.
971          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', or 'trapezodial'.
972          * @param {Number} start Left border of the approximation interval
973          * @param {Number} end Right border of the approximation interval
974          * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This
975          * array may be used as parent array of a {@link JXG.Curve}.
976          */
977         riemann: function(f, n, type, start, end) {
978             var xarr = [],
979                 yarr = [],
980                 i, j = 0,
981                 delta,
982                 x = start, y,
983                 x1, y1, delta1;
984 
985             n = Math.floor(n);
986             xarr[j] = x;
987             yarr[j] = 0.0;
988 
989             if (n > 0) {
990                 delta = (end - start) / n;
991                 delta1 = delta * 0.01; // for 'lower' and 'upper'
992 
993                 for (i = 0; i < n; i++) {
994                     if (type === 'right') {
995                         y = f(x + delta);
996                     } else if (type === 'middle') {
997                         y = f(x + delta * 0.5);
998                     } else if ((type === 'left') || (type === 'trapezodial')) {
999                         y = f(x);
1000                     } else if (type === 'lower') {
1001                         y = f(x);
1002                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1003                             y1 = f(x1);
1004                             if (y1 < y) {
1005                                 y = y1;
1006                             }
1007                         }
1008                     } else { // (type=='upper')
1009                         y = f(x);
1010                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1011                             y1 = f(x1);
1012                             if (y1 > y) {
1013                                 y = y1;
1014                             }
1015                         }
1016                     }
1017 
1018                     j++;
1019                     xarr[j] = x;
1020                     yarr[j] = y;
1021                     j++;
1022                     x += delta;
1023                     if (type === 'trapezodial') {
1024                         y = f(x);
1025                     }
1026                     xarr[j] = x;
1027                     yarr[j] = y;
1028                     j++;
1029                     xarr[j] = x;
1030                     yarr[j] = 0.0;
1031                 }
1032             }
1033             return [xarr,yarr];
1034         },
1035 
1036         /**
1037          * Approximate the integral by Riemann sums.
1038          * Compute the area described by the riemann sum rectangles.
1039          * @param {function} f Function, whose integral is approximated by the Riemann sum.
1040          * @param {Number} n number of rectangles.
1041          * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', or 'trapezodial'.
1042          * @param {Number} start Left border of the approximation interval
1043          * @param {Number} end Right border of the approximation interval
1044          * @returns {Number} The sum of the areas of the rectangles.
1045          */
1046         riemannsum: function(f, n, type, start, end) {
1047             var sum = .0,
1048                 i, delta,
1049                 x = start, y,
1050                 x1, y1, delta1;
1051 
1052             n = Math.floor(n);
1053             if (n > 0) {
1054                 delta = (end - start) / n;
1055                 delta1 = delta * 0.01; // for 'lower' and 'upper'
1056                 for (i = 0; i < n; i++) {
1057                     if (type === 'right') {
1058                         y = f(x + delta);
1059                     } else if (type === 'middle') {
1060                         y = f(x + delta * 0.5);
1061                     } else if (type === 'trapezodial') {
1062                         y = 0.5 * (f(x + delta) + f(x));
1063                     } else if (type === 'left') {
1064                         y = f(x);
1065                     } else if (type === 'lower') {
1066                         y = f(x);
1067                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1068                             y1 = f(x1);
1069                             if (y1 < y) {
1070                                 y = y1;
1071                             }
1072                         }
1073                     } else { // (type=='upper')
1074                         y = f(x);
1075                         for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) {
1076                             y1 = f(x1);
1077                             if (y1 > y) {
1078                                 y = y1;
1079                             }
1080                         }
1081                     }
1082                     sum += delta * y;
1083                     x += delta;
1084                 }
1085             }
1086             return sum;
1087         },
1088 
1089         /**
1090          * Solve initial value problems numerically using Runge-Kutta-methods.
1091          * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm.
1092          * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing
1093          * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure
1094          * <pre>
1095          * {
1096          *     s: <Number>,
1097          *     A: <matrix>,
1098          *     b: <Array>,
1099          *     c: <Array>
1100          * }
1101          * </pre>
1102          * 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
1103          * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array.
1104          * @param {Array} I Interval on which to integrate.
1105          * @param {Number} N Number of evaluation points.
1106          * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode
1107          * 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
1108          * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has.
1109          * @returns {Array} An array of vectors describing the solution of the ode on the given interval I.
1110          * @example
1111          * // A very simple autonomous system dx(t)/dt = x(t);
1112          * function f(t, x) {
1113          *     return x;
1114          * }
1115          *
1116          * // Solve it with initial value x(0) = 1 on the interval [0, 2]
1117          * // with 20 evaluation points.
1118          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
1119          *
1120          * // Prepare data for plotting the solution of the ode using a curve.
1121          * var dataX = [];
1122          * var dataY = [];
1123          * var h = 0.1;        // (I[1] - I[0])/N  = (2-0)/20
1124          * for(var i=0; i<data.length; i++) {
1125          *     dataX[i] = i*h;
1126          *     dataY[i] = data[i][0];
1127          * }
1128          * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'});
1129          * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div>
1130          * <script type="text/javascript">
1131          * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false});
1132          * function f(t, x) {
1133          *     // we have to copy the value.
1134          *     // return x; would just return the reference.
1135          *     return [x[0]];
1136          * }
1137          * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f);
1138          * var dataX = [];
1139          * var dataY = [];
1140          * var h = 0.1;
1141          * for(var i=0; i<data.length; i++) {
1142          *     dataX[i] = i*h;
1143          *     dataY[i] = data[i][0];
1144          * }
1145          * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'});
1146          * </script><pre>
1147          */
1148         rungeKutta: function(butcher, x0, I, N, f) {
1149             var x = [],
1150                 y = [],
1151                 h = (I[1] - I[0]) / N,
1152                 t = I[0],
1153                 e, i, j,
1154                 k, l,
1155                 dim = x0.length,
1156                 s,
1157                 result = [],
1158                 r = 0;
1159 
1160             if (JXG.isString(butcher)) {
1161                 butcher = predefinedButcher[butcher] || predefinedButcher.euler;
1162             }
1163             s = butcher.s;
1164 
1165             // don't change x0, so copy it
1166             for (e = 0; e < dim; e++)
1167                 x[e] = x0[e];
1168 
1169             for (i = 0; i < N; i++) {
1170                 // Optimization doesn't work for ODEs plotted using time
1171                 //        if((i % quotient == 0) || (i == N-1)) {
1172                 result[r] = [];
1173                 for (e = 0; e < dim; e++)
1174                     result[r][e] = x[e];
1175                 r++;
1176                 //        }
1177                 // init k
1178                 k = [];
1179                 for (j = 0; j < s; j++) {
1180                     // init y = 0
1181                     for (e = 0; e < dim; e++)
1182                         y[e] = 0.;
1183 
1184                     // Calculate linear combination of former k's and save it in y
1185                     for (l = 0; l < j; l++) {
1186                         for (e = 0; e < dim; e++) {
1187                             y[e] += (butcher.A[j][l]) * h * k[l][e];
1188                         }
1189                     }
1190 
1191                     // add x(t) to y
1192                     for (e = 0; e < dim; e++) {
1193                         y[e] += x[e];
1194                     }
1195 
1196                     // calculate new k and add it to the k matrix
1197                     k.push(f(t + butcher.c[j] * h, y));
1198                 }
1199 
1200                 // init y = 0
1201                 for (e = 0; e < dim; e++)
1202                     y[e] = 0.;
1203 
1204                 for (l = 0; l < s; l++) {
1205                     for (e = 0; e < dim; e++)
1206                         y[e] += butcher.b[l] * k[l][e];
1207                 }
1208 
1209                 for (e = 0; e < dim; e++) {
1210                     x[e] = x[e] + h * y[e];
1211                 }
1212 
1213                 t += h;
1214             }
1215 
1216             return result;
1217         },
1218 
1219 
1220         /**
1221          *
1222          * Find zero of an univariate function f.
1223          * @param {function} f Function, whose root is to be found
1224          * @param {array or number} x0  Start value or start interval enclosing the root
1225          * @param {object} object Parent object in case f is method of it
1226          * @return {number} the approximation of the root
1227          * Algorithm:
1228          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
1229          *  computations. M., Mir, 1980, p.180 of the Russian edition
1230          * 
1231          * if x0 is an array containing lower and upper bound for the zero
1232          * algorithm 748 is applied. Otherwise, if x0 is a number,
1233          * the algorithm tries to bracket a zero of f starting from x0.
1234          * If this fails, we fall back to Newton's method.
1235          *
1236          **/
1237         fzero: function(f, x0, object) {
1238             var tol = JXG.Math.eps,
1239                 maxiter = 50, niter = 0,
1240                 nfev = 0,
1241                 eps = tol,
1242                 a,b,c, 
1243                 fa, fb, fc,
1244                 aa, blist, i, len, u, fu, 
1245                 prev_step,
1246                 tol_act,         // Actual tolerance
1247                 p,               // Interpolation step is calcu-
1248                 q,               // lated in the form p/q; divi- 
1249                                  // sion operations is delayed  
1250                                  // until the last moment   
1251                 new_step,        // Step at this iteration
1252                 t1, cb, t2;
1253 
1254             if (JXG.isArray(x0)) {
1255                 if (x0.length<2) 
1256                     throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two.");
1257                 a = x0[0]; fa = f.apply(object,[a]); nfev++;
1258                 b = x0[1]; fb = f.apply(object,[b]); nfev++;
1259             } else {
1260                 a = x0; fa = f.apply(object,[a]); nfev++;
1261                 // Try to get b.
1262                 if (a == 0) {
1263                     aa = 1;
1264                 } else {
1265                     aa = a;
1266                 }
1267                 blist = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa, 1.5*aa, -aa, 2*aa, -10*aa, 10*aa];
1268                 len = blist.length;
1269                 for (i=0;i<len;i++) {
1270                     b = blist[i];
1271                     fb = f.apply(object,[b]); nfev++;
1272                     if (fa*fb<=0) {
1273                         break;
1274                     }
1275                 }
1276                 if (b < a) {
1277                     u = a; a = b; b = u;
1278                     fu = fa; fa = fb; fb = fu;
1279                 }
1280             }
1281 
1282             if (fa*fb > 0) {
1283                 // Bracketing not successful, fall back to Newton's method or to fminbr
1284                 if (JXG.isArray(x0)) {
1285                     //JXG.debug("fzero falls back to fminbr");
1286                     return this.fminbr(f, [a,b], object);
1287                 } else {
1288                     //JXG.debug("fzero falls back to Newton");
1289                     return this.Newton(f, a, object);
1290                 }
1291             }
1292 
1293             // OK, we have enclosed a zero of f.
1294             // Now we can start Brent's method
1295 
1296             c = a;   
1297             fc = fa;
1298             while (niter<maxiter) {   // Main iteration loop
1299                 prev_step = b-a;      // Distance from the last but one
1300                                       // to the last approximation 
1301 
1302                 if ( Math.abs(fc) < Math.abs(fb) ) {
1303                                                 // Swap data for b to be the  
1304                     a = b;  b = c;  c = a;      // best approximation   
1305                     fa=fb;  fb=fc;  fc=fa;
1306                 }
1307                 tol_act = 2*eps*Math.abs(b) + tol*0.5;
1308                 new_step = (c-b)*0.5;
1309 
1310                 if ( Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps ) {
1311                     //JXG.debug("nfev="+nfev);
1312                     return b;                           //  Acceptable approx. is found 
1313                 }
1314                     
1315                 // Decide if the interpolation can be tried 
1316                 if ( Math.abs(prev_step) >= tol_act     // If prev_step was large enough
1317                     && Math.abs(fa) > Math.abs(fb) ) {  // and was in true direction,  
1318                                                         // Interpolatiom may be tried  
1319                     cb = c-b;
1320                     if ( a==c ) {                       // If we have only two distinct 
1321                                                         // points linear interpolation 
1322                         t1 = fb/fa;                     // can only be applied    
1323                         p = cb*t1;
1324                         q = 1.0 - t1;
1325                     } else {                            // Quadric inverse interpolation
1326                         q = fa/fc;  t1 = fb/fc;  t2 = fb/fa;
1327                         p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
1328                         q = (q-1.0) * (t1-1.0) * (t2-1.0);
1329                     }
1330                     if ( p>0 ) {                        // p was calculated with the op-
1331                         q = -q;                         // posite sign; make p positive 
1332                     } else {                            // and assign possible minus to
1333                         p = -p;                         // q 
1334                     }
1335                     
1336                     if( p < (0.75*cb*q-Math.abs(tol_act*q)*0.5)     // If b+p/q falls in [b,c]
1337                         && p < Math.abs(prev_step*q*0.5) ) {        // and isn't too large 
1338                         new_step = p/q;                 // it is accepted   
1339                     }                                   // If p/q is too large then the 
1340                                                         // bissection procedure can
1341                                                         // reduce [b,c] range to more
1342                                                         // extent
1343                 }
1344 
1345                 if ( Math.abs(new_step) < tol_act ) {   // Adjust the step to be not less
1346                     if ( new_step > 0 ) {               // than tolerance
1347                         new_step = tol_act;
1348                     } else {
1349                         new_step = -tol_act;
1350                     }
1351                 }
1352 
1353                 a = b;  fa = fb;                        // Save the previous approx.
1354                 b += new_step;  
1355                 fb = f.apply(object,[b]); nfev++;       // Do step to a new approxim.
1356                 if ( (fb>0 && fc>0) || (fb<0 && fc<0) ) {
1357                                                         // Adjust c for it to have a sign
1358                     c = a;  fc = fa;                    // opposite to that of b 
1359                 }
1360                 niter++;
1361             }                                           // End while
1362             
1363             JXG.debug("fzero: maxiter="+maxiter+" reached.");
1364             return b;
1365         },
1366 
1367         /**
1368          *
1369          * Find minimum of an univariate function f.
1370          * @param {function} f Function, whose minimum is to be found
1371          * @param {array} x0  Start interval enclosing the minimum
1372          * @param {object} object Parent object in case f is method of it
1373          * @return {number} the approximation of the minimum
1374          * Algorithm:
1375          *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
1376          *  computations. M., Mir, 1980, p.180 of the Russian edition
1377          * x0 
1378          **/
1379 
1380         fminbr: function(f, x0, object) {              // An estimate to the min location
1381             var a, b, x, v, w,
1382                 fx, fv, fw,
1383                 r = (3.-Math.sqrt(5.0))*0.5,            // Golden section ratio   
1384                 tol = JXG.Math.eps,
1385                 sqrteps = Math.sqrt(JXG.Math.eps),
1386                 maxiter = 50, niter = 0,
1387                 range, middle_range, tol_act, new_step,
1388                 p, q, t, ft,
1389                 nfev = 0;
1390 
1391             if (!JXG.isArray(x0) || x0.length<2) {
1392                 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two.");
1393             }
1394             
1395             a = x0[0];
1396             b = x0[1];
1397             v = a + r*(b-a);  
1398             fv = f.apply(object,[v]); nfev++;           // First step - always gold section
1399             x = v;  w = v;
1400             fx=fv;  fw=fv;
1401 
1402             while (niter<maxiter) {
1403                 range = b-a;                            // Range over which the minimum 
1404                                                         // is seeked for
1405                 middle_range = (a+b)*0.5;
1406                 tol_act = sqrteps*Math.abs(x) + tol/3;  // Actual tolerance  
1407                 if( Math.abs(x-middle_range) + range*0.5 <= 2*tol_act ) {
1408                     //JXG.debug(nfev);
1409                     return x;                           // Acceptable approx. is found
1410                 }
1411                                                         // Obtain the golden section step 
1412                 new_step = r * ( x<middle_range ? b-x : a-x );
1413                                                         // Decide if the interpolation can be tried 
1414                 if ( Math.abs(x-w) >= tol_act  ) {      // If x and w are distinct     
1415                                                         // interpolatiom may be tried 
1416                     // Interpolation step is calculated as p/q; 
1417                     // division operation is delayed until last moment 
1418                     t = (x-w) * (fx-fv);
1419                     q = (x-v) * (fx-fw);
1420                     p = (x-v)*q - (x-w)*t;
1421                     q = 2*(q-t);
1422 
1423                     if ( q>0 ) {                        // q was calculated with the op-
1424                         p = -p;                         // posite sign; make q positive 
1425                     } else {                            // and assign possible minus to
1426                         q = -q;                         // p
1427                     }
1428                     if ( Math.abs(p) < Math.abs(new_step*q) &&      // If x+p/q falls in [a,b]
1429                          p > q*(a-x+2*tol_act) &&                   //  not too close to a and
1430                          p < q*(b-x-2*tol_act)  ) {                 // b, and isn't too large */
1431                          new_step = p/q;                            // it is accepted        
1432                     }
1433                     // If p/q is too large then the 
1434                     // golden section procedure can   
1435                     // reduce [a,b] range to more   
1436                     // extent           
1437                 }
1438 
1439                 if ( Math.abs(new_step) < tol_act ) {    // Adjust the step to be not less
1440                     if( new_step > 0 ) {                 // than tolerance     
1441                         new_step = tol_act;
1442                     } else {
1443                         new_step = -tol_act;
1444                     }
1445                 }
1446                 
1447                 // Obtain the next approximation to min 
1448                 // and reduce the enveloping range
1449                 t = x + new_step;                       // Tentative point for the min
1450                 ft = f.apply(object,[t]); nfev++;
1451                 if ( ft <= fx ) {                       // t is a better approximation 
1452                     if ( t < x ) {                      // Reduce the range so that
1453                         b = x;                          // t would fall within it 
1454                     } else {
1455                         a = x;
1456                     }
1457                     v = w;  w = x;  x = t;              // Assign the best approx to x 
1458                     fv=fw;  fw=fx;  fx=ft;
1459                 } else {                                // x remains the better approx
1460                     if ( t < x ) {                      // Reduce the range enclosing x 
1461                         a = t;                   
1462                     } else {
1463                         b = t;
1464                     }
1465                     if ( ft <= fw || w==x ) {
1466                         v = w;  w = t;
1467                         fv=fw;  fw=ft;
1468                     } else if ( ft<=fv || v==x || v==w ) {
1469                         v = t;
1470                         fv=ft;
1471                     }
1472                 }
1473                 niter++;
1474             } 
1475             JXG.debug("fminbr: maxiter="+maxiter+" reached.");
1476             return x;
1477         },
1478 
1479         /**
1480          * Helper function to create curve which displays Reuleaux polygons.
1481          * @param {array} points Array of points which should be the vertices of the Reuleaux polygon. Typically,
1482          *                       these point list is the array vrtices of a regular polygon.
1483          * @param {number} nr Number of vertices
1484          * @returns {array} An array containing the two functions defining the Reuleaux polygon and the two values
1485          * for the start and the end of the paramtric curve.
1486          * array may be used as parent array of a {@link JXG.Curve}.
1487          *
1488          * @example
1489          * var A = brd.create('point',[-2,-2]);
1490          * var B = brd.create('point',[0,1]);
1491          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 
1492          * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 
1493          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
1494          *
1495          * </pre><div id="2543a843-46a9-4372-abc1-94d9ad2db7ac" style="width: 300px; height: 300px;"></div>
1496          * <script type="text/javascript">
1497          * var brd = JXG.JSXGraph.initBoard('2543a843-46a9-4372-abc1-94d9ad2db7ac', {boundingbox: [-5, 5, 5, -5], axis: true, showcopyright:false, shownavigation: false});
1498          * var A = brd.create('point',[-2,-2]);
1499          * var B = brd.create('point',[0,1]);
1500          * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 
1501          * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 
1502          *                          {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'});
1503          * </script><pre>
1504          */
1505         reuleauxPolygon: function(points, nr) {
1506     var pi2 = Math.PI*2,
1507         pi2_n = pi2/nr,
1508         diag = (nr-1)/2,
1509         beta, d = 0,
1510         makeFct = function(which, trig) {
1511                 return function(t, suspendUpdate) {
1512                     if (!suspendUpdate) {
1513                         d = points[0].Dist(points[diag]);
1514                         beta = JXG.Math.Geometry.rad([points[0].X()+1,points[0].Y()],points[0],points[(diag)%nr]);
1515                     }
1516                     var t1 = (t%pi2 + pi2) % pi2;
1517                     var j = Math.floor(t1 / pi2_n)%nr;
1518                     if (isNaN(j)) return j;
1519                     //t1 = (t1-j*pi2_n)*0.5 + beta+j*pi2_n;
1520                     t1 = t1*0.5+j*pi2_n*0.5 + beta;
1521                     return points[j][which]()+d*Math[trig](t1);
1522                 };
1523             };
1524     return [
1525             makeFct('X','cos'),
1526             makeFct('Y','sin'),
1527             0,
1528             Math.PI*2
1529         ];
1530         }        
1531 
1532     }
1533 })(JXG, Math);
1534 
1535