1 /*
  2     Copyright 2008-2022
  3         Matthias Ehmann,
  4         Michael Gerhaeuser,
  5         Carsten Miller,
  6         Alfred Wassermann
  7 
  8     This file is part of JSXGraph.
  9 
 10     JSXGraph is free software dual licensed under the GNU LGPL or MIT License.
 11 
 12     You can redistribute it and/or modify it under the terms of the
 13 
 14       * GNU Lesser General Public License as published by
 15         the Free Software Foundation, either version 3 of the License, or
 16         (at your option) any later version
 17       OR
 18       * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT
 19 
 20     JSXGraph is distributed in the hope that it will be useful,
 21     but WITHOUT ANY WARRANTY; without even the implied warranty of
 22     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 23     GNU Lesser General Public License for more details.
 24 
 25     You should have received a copy of the GNU Lesser General Public License and
 26     the MIT License along with JSXGraph. If not, see <http://www.gnu.org/licenses/>
 27     and <http://opensource.org/licenses/MIT/>.
 28  */
 29 
 30 
 31 /*global JXG: true, define: true*/
 32 /*jslint nomen: true, plusplus: true*/
 33 
 34 /* depends:
 35  jxg
 36  math/math
 37  utils/type
 38  */
 39 
 40 define(['jxg', 'base/constants', 'base/coords', 'math/math', 'math/extrapolate', 'math/numerics',
 41         'math/statistics', 'math/geometry', 'math/ia', 'utils/type'],
 42         function (JXG, Const, Coords, Mat, Extrapolate, Numerics, Statistics, Geometry, IntervalArithmetic, Type) {
 43 
 44     "use strict";
 45 
 46     /**
 47      * Functions for plotting of curves.
 48      * @name JXG.Math.Plot
 49      * @exports Mat.Plot as JXG.Math.Plot
 50      * @namespace
 51      */
 52     Mat.Plot = {
 53 
 54         /**
 55          * Check if at least one point on the curve is finite and real.
 56          **/
 57         checkReal: function (points) {
 58             var b = false,
 59                 i, p,
 60                 len = points.length;
 61 
 62             for (i = 0; i < len; i++) {
 63                 p = points[i].usrCoords;
 64                 if (!isNaN(p[1]) && !isNaN(p[2]) && Math.abs(p[0]) > Mat.eps) {
 65                     b = true;
 66                     break;
 67                 }
 68             }
 69             return b;
 70         },
 71 
 72         //----------------------------------------------------------------------
 73         // Plot algorithm v0
 74         //----------------------------------------------------------------------
 75         /**
 76          * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#doadvancedplot} is <tt>false</tt>.
 77          * @param {JXG.Curve} curve JSXGraph curve element
 78          * @param {Number} mi Left bound of curve
 79          * @param {Number} ma Right bound of curve
 80          * @param {Number} len Number of data points
 81          * @returns {JXG.Curve} Reference to the curve object.
 82          */
 83         updateParametricCurveNaive: function (curve, mi, ma, len) {
 84             var i, t,
 85                 suspendUpdate = false,
 86                 stepSize = (ma - mi) / len;
 87 
 88             for (i = 0; i < len; i++) {
 89                 t = mi + i * stepSize;
 90                 // The last parameter prevents rounding in usr2screen().
 91                 curve.points[i].setCoordinates(Const.COORDS_BY_USER, [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], false);
 92                 curve.points[i]._t = t;
 93                 suspendUpdate = true;
 94             }
 95             return curve;
 96         },
 97 
 98         //----------------------------------------------------------------------
 99         // Plot algorithm v1
100         //----------------------------------------------------------------------
101         /**
102          * Crude and cheap test if the segment defined by the two points <tt>(x0, y0)</tt> and <tt>(x1, y1)</tt> is
103          * outside the viewport of the board. All parameters have to be given in screen coordinates.
104          *
105          * @private
106          * @deprecated
107          * @param {Number} x0
108          * @param {Number} y0
109          * @param {Number} x1
110          * @param {Number} y1
111          * @param {JXG.Board} board
112          * @returns {Boolean} <tt>true</tt> if the given segment is outside the visible area.
113          */
114         isSegmentOutside: function (x0, y0, x1, y1, board) {
115             return (y0 < 0 && y1 < 0) || (y0 > board.canvasHeight && y1 > board.canvasHeight) ||
116                 (x0 < 0 && x1 < 0) || (x0 > board.canvasWidth && x1 > board.canvasWidth);
117         },
118 
119         /**
120          * Compares the absolute value of <tt>dx</tt> with <tt>MAXX</tt> and the absolute value of <tt>dy</tt>
121          * with <tt>MAXY</tt>.
122          *
123          * @private
124          * @deprecated
125          * @param {Number} dx
126          * @param {Number} dy
127          * @param {Number} MAXX
128          * @param {Number} MAXY
129          * @returns {Boolean} <tt>true</tt>, if <tt>|dx| < MAXX</tt> and <tt>|dy| < MAXY</tt>.
130          */
131         isDistOK: function (dx, dy, MAXX, MAXY) {
132             return (Math.abs(dx) < MAXX && Math.abs(dy) < MAXY) && !isNaN(dx + dy);
133         },
134 
135          /**
136          * @private
137          * @deprecated
138          */
139         isSegmentDefined: function (x0, y0, x1, y1) {
140             return !(isNaN(x0 + y0) && isNaN(x1 + y1));
141         },
142 
143         /**
144          * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#doadvancedplot} is <tt>true</tt>.
145          * Since 0.99 this algorithm is deprecated. It still can be used if {@link JXG.Curve#doadvancedplotold} is <tt>true</tt>.
146          *
147          * @deprecated
148          * @param {JXG.Curve} curve JSXGraph curve element
149          * @param {Number} mi Left bound of curve
150          * @param {Number} ma Right bound of curve
151          * @returns {JXG.Curve} Reference to the curve object.
152          */
153         updateParametricCurveOld: function (curve, mi, ma) {
154             var i, t, d,
155                 x, y, t0, x0, y0, top, depth,
156                 MAX_DEPTH, MAX_XDIST, MAX_YDIST,
157                 suspendUpdate = false,
158                 po = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false),
159                 dyadicStack = [],
160                 depthStack = [],
161                 pointStack = [],
162                 divisors = [],
163                 distOK = false,
164                 j = 0,
165                 distFromLine = function (p1, p2, p0) {
166                     var lbda, d,
167                         x0 = p0[1] - p1[1],
168                         y0 = p0[2] - p1[2],
169                         x1 = p2[0] - p1[1],
170                         y1 = p2[1] - p1[2],
171                         den = x1 * x1 + y1 * y1;
172 
173                     if (den >= Mat.eps) {
174                         lbda = (x0 * x1 + y0 * y1) / den;
175                         if (lbda > 0) {
176                             if (lbda <= 1) {
177                                 x0 -= lbda * x1;
178                                 y0 -= lbda * y1;
179                             // lbda = 1.0;
180                             } else {
181                                 x0 -= x1;
182                                 y0 -= y1;
183                             }
184                         }
185                     }
186                     d = x0 * x0 + y0 * y0;
187                     return Math.sqrt(d);
188                 };
189 
190             JXG.deprecated('Curve.updateParametricCurveOld()');
191 
192             if (curve.board.updateQuality === curve.board.BOARD_QUALITY_LOW) {
193                 MAX_DEPTH = 15;
194                 MAX_XDIST = 10; // 10
195                 MAX_YDIST = 10; // 10
196             } else {
197                 MAX_DEPTH = 21;
198                 MAX_XDIST = 0.7; // 0.7
199                 MAX_YDIST = 0.7; // 0.7
200             }
201 
202             divisors[0] = ma - mi;
203             for (i = 1; i < MAX_DEPTH; i++) {
204                 divisors[i] = divisors[i - 1] * 0.5;
205             }
206 
207             i = 1;
208             dyadicStack[0] = 1;
209             depthStack[0] = 0;
210 
211             t = mi;
212             po.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], false);
213 
214             // Now, there was a first call to the functions defining the curve.
215             // Defining elements like sliders have been evaluated.
216             // Therefore, we can set suspendUpdate to false, so that these defining elements
217             // need not be evaluated anymore for the rest of the plotting.
218             suspendUpdate = true;
219             x0 = po.scrCoords[1];
220             y0 = po.scrCoords[2];
221             t0 = t;
222 
223             t = ma;
224             po.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], false);
225             x = po.scrCoords[1];
226             y = po.scrCoords[2];
227 
228             pointStack[0] = [x, y];
229 
230             top = 1;
231             depth = 0;
232 
233             curve.points = [];
234             curve.points[j++] = new Coords(Const.COORDS_BY_SCREEN, [x0, y0], curve.board, false);
235 
236             do {
237                 distOK = this.isDistOK(x - x0, y - y0, MAX_XDIST, MAX_YDIST) || this.isSegmentOutside(x0, y0, x, y, curve.board);
238                 while (depth < MAX_DEPTH && (!distOK || depth < 6) && (depth <= 7 || this.isSegmentDefined(x0, y0, x, y))) {
239                     // We jump out of the loop if
240                     // * depth>=MAX_DEPTH or
241                     // * (depth>=6 and distOK) or
242                     // * (depth>7 and segment is not defined)
243 
244                     dyadicStack[top] = i;
245                     depthStack[top] = depth;
246                     pointStack[top] = [x, y];
247                     top += 1;
248 
249                     i = 2 * i - 1;
250                     // Here, depth is increased and may reach MAX_DEPTH
251                     depth++;
252                     // In that case, t is undefined and we will see a jump in the curve.
253                     t = mi + i * divisors[depth];
254 
255                     po.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, suspendUpdate), curve.Y(t, suspendUpdate)], false, true);
256                     x = po.scrCoords[1];
257                     y = po.scrCoords[2];
258                     distOK = this.isDistOK(x - x0, y - y0, MAX_XDIST, MAX_YDIST) || this.isSegmentOutside(x0, y0, x, y, curve.board);
259                 }
260 
261                 if (j > 1) {
262                     d = distFromLine(curve.points[j - 2].scrCoords, [x, y], curve.points[j - 1].scrCoords);
263                     if (d < 0.015) {
264                         j -= 1;
265                     }
266                 }
267 
268                 curve.points[j] = new Coords(Const.COORDS_BY_SCREEN, [x, y], curve.board, false);
269                 curve.points[j]._t = t;
270                 j += 1;
271 
272                 x0 = x;
273                 y0 = y;
274                 t0 = t;
275 
276                 top -= 1;
277                 x = pointStack[top][0];
278                 y = pointStack[top][1];
279                 depth = depthStack[top] + 1;
280                 i = dyadicStack[top] * 2;
281 
282             } while (top > 0 && j < 500000);
283 
284             curve.numberPoints = curve.points.length;
285 
286             return curve;
287         },
288 
289         //----------------------------------------------------------------------
290         // Plot algorithm v2
291         //----------------------------------------------------------------------
292 
293         /**
294          * Add a point to the curve plot. If the new point is too close to the previously inserted point,
295          * it is skipped.
296          * Used in {@link JXG.Curve._plotRecursive}.
297          *
298          * @private
299          * @param {JXG.Coords} pnt Coords to add to the list of points
300          */
301         _insertPoint_v2: function (curve, pnt, t) {
302             var lastReal = !isNaN(this._lastCrds[1] + this._lastCrds[2]),     // The last point was real
303                 newReal = !isNaN(pnt.scrCoords[1] + pnt.scrCoords[2]),        // New point is real point
304                 cw = curve.board.canvasWidth,
305                 ch = curve.board.canvasHeight,
306                 off = 500;
307 
308             newReal = newReal &&
309                         (pnt.scrCoords[1] > -off && pnt.scrCoords[2] > -off &&
310                          pnt.scrCoords[1] < cw + off && pnt.scrCoords[2] < ch + off);
311 
312             /*
313              * Prevents two consecutive NaNs or points wich are too close
314              */
315             if ((!newReal && lastReal) ||
316                     (newReal && (!lastReal ||
317                         Math.abs(pnt.scrCoords[1] - this._lastCrds[1]) > 0.7 ||
318                         Math.abs(pnt.scrCoords[2] - this._lastCrds[2]) > 0.7))) {
319                 pnt._t = t;
320                 curve.points.push(pnt);
321                 this._lastCrds = pnt.copy('scrCoords');
322             }
323         },
324 
325         /**
326          * Check if there is a single NaN function value at t0.
327          * @param {*} curve
328          * @param {*} t0
329          * @returns {Boolean} true if there is a second NaN point close by, false otherwise
330          */
331         neighborhood_isNaN_v2: function(curve, t0) {
332             var is_undef,
333                 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false),
334                 t, p;
335 
336             t = t0 + Mat.eps;
337             pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, true), curve.Y(t, true)], false);
338             p = pnt.usrCoords;
339             is_undef = isNaN(p[1] + p[2]);
340             if (!is_undef) {
341                 t = t0 - Mat.eps;
342                 pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, true), curve.Y(t, true)], false);
343                 p = pnt.usrCoords;
344                 is_undef = isNaN(p[1] + p[2]);
345                 if (!is_undef) {
346                     return false;
347                 }
348             }
349             return true;
350         },
351 
352         /**
353          * Investigate a function term at the bounds of intervals where
354          * the function is not defined, e.g. log(x) at x = 0.
355          *
356          * c is inbetween a and b
357          * @private
358          * @param {JXG.Curve} curve JSXGraph curve element
359          * @param {Array} a Screen coordinates of the left interval bound
360          * @param {Array} b Screen coordinates of the right interval bound
361          * @param {Array} c Screen coordinates of the bisection point at (ta + tb) / 2
362          * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates
363          * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates
364          * @param {Number} tc (ta + tb) / 2 = tc. Parameter which evaluates to b, i.e. [1, X(tc), Y(tc)] = c in screen coordinates
365          * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0.
366          * @returns {JXG.Boolean} true if the point is inserted and the recursion should stop, false otherwise.
367          */
368         _borderCase: function (curve, a, b, c, ta, tb, tc, depth) {
369             var t, pnt, p,
370                 p_good = null,
371                 j,
372                 max_it = 30,
373                 is_undef = false,
374                 t_nan, t_real, t_real2,
375                 vx, vy, vx2, vy2, dx, dy;
376                 // asymptote;
377 
378             if (depth <= 1) {
379                 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false);
380                 // Test if there is a single undefined point.
381                 // If yes, we ignore it.
382                 if (isNaN(a[1] + a[2]) && !isNaN(c[1] + c[2]) && !this.neighborhood_isNaN_v2(curve, ta)) {
383                     return false;
384                 }
385                 if (isNaN(b[1] + b[2]) && !isNaN(c[1] + c[2]) && !this.neighborhood_isNaN_v2(curve, tb)) {
386                     return false;
387                 }
388                 if (isNaN(c[1] + c[2]) && (!isNaN(a[1] + a[2]) || !isNaN(b[1] + b[2])) &&
389                     !this.neighborhood_isNaN_v2(curve, tc)) {
390                     return false;
391                 }
392 
393                 j = 0;
394                 // Bisect a, b and c until the point t_real is inside of the definition interval
395                 // and as close as possible at the boundary.
396                 // t_real2 is the second closest point.
397                 do {
398                     // There are four cases:
399                     //  a  |  c  |  b
400                     // ---------------
401                     // inf | R   | R
402                     // R   | R   | inf
403                     // inf | inf | R
404                     // R   | inf | inf
405                     //
406                     if (isNaN(a[1] + a[2]) && !isNaN(c[1] + c[2])) {
407                         t_nan = ta;
408                         t_real = tc;
409                         t_real2 = tb;
410                     } else if (isNaN(b[1] + b[2]) && !isNaN(c[1] + c[2])) {
411                         t_nan = tb;
412                         t_real = tc;
413                         t_real2 = ta;
414                     } else if (isNaN(c[1] + c[2]) && !isNaN(b[1] + b[2])) {
415                         t_nan = tc;
416                         t_real = tb;
417                         t_real2 = tb + (tb - tc);
418                     } else if (isNaN(c[1] + c[2]) && !isNaN(a[1] + a[2])) {
419                         t_nan = tc;
420                         t_real = ta;
421                         t_real2 = ta - (tc - ta);
422                     } else {
423                         return false;
424                     }
425                     t = 0.5 * (t_nan + t_real);
426                     pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, true), curve.Y(t, true)], false);
427                     p = pnt.usrCoords;
428 
429                     is_undef = isNaN(p[1] + p[2]);
430                     if (is_undef) {
431                         t_nan = t;
432                     } else {
433                         t_real2 = t_real;
434                         t_real = t;
435                     }
436                     ++j;
437                 } while (is_undef && j < max_it);
438 
439                 // If bisection was successful, take this point.
440                 // Useful only for general curves, for function graph
441                 // the code below overwrite p_good from here.
442                 if (j < max_it) {
443                     p_good = p.slice();
444                     c = p.slice();
445                     t_real = t;
446                 }
447 
448                 // OK, bisection has been done now.
449                 // t_real contains the closest inner point to the border of the interval we could find.
450                 // t_real2 is the second nearest point to this boundary.
451                 // Now we approximate the derivative by computing the slope of the line through these two points
452                 // and test if it is "infinite", i.e larger than 400 in absolute values.
453                 //
454                 vx = curve.X(t_real, true) ;
455                 vx2 = curve.X(t_real2, true) ;
456                 dx = (vx - vx2) / (t_real - t_real2);
457                 vy = curve.Y(t_real, true) ;
458                 vy2 = curve.Y(t_real2, true) ;
459                 dy = (vy - vy2) / (t_real - t_real2);
460 
461                 if (p_good !== null) {
462                     this._insertPoint_v2(curve, new Coords(Const.COORDS_BY_USER, p_good, curve.board, false));
463                     return true;
464                 }
465            }
466            return false;
467        },
468 
469         /**
470          * Recursive interval bisection algorithm for curve plotting.
471          * Used in {@link JXG.Curve.updateParametricCurve}.
472          * @private
473          * @deprecated
474          * @param {JXG.Curve} curve JSXGraph curve element
475          * @param {Array} a Screen coordinates of the left interval bound
476          * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates
477          * @param {Array} b Screen coordinates of the right interval bound
478          * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates
479          * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0.
480          * @param {Number} delta If the distance of the bisection point at (ta + tb) / 2 from the point (a + b) / 2 is less then delta,
481          *                 the segment [a,b] is regarded as straight line.
482          * @returns {JXG.Curve} Reference to the curve object.
483          */
484         _plotRecursive_v2: function (curve, a, ta, b, tb, depth, delta) {
485             var tc, c,
486                 ds, mindepth = 0,
487                 isSmooth, isJump, isCusp,
488                 cusp_threshold = 0.5,
489                 jump_threshold = 0.99,
490                 pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false);
491 
492             if (curve.numberPoints > 65536) {
493                 return;
494             }
495 
496             // Test if the function is undefined in an interval
497             if (depth < this.nanLevel && this._isUndefined(curve, a, ta, b, tb)) {
498                 return this;
499             }
500 
501             if (depth < this.nanLevel && this._isOutside(a, ta, b, tb, curve.board)) {
502                 return this;
503             }
504 
505             tc = (ta  + tb) * 0.5;
506             pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(tc, true), curve.Y(tc, true)], false);
507             c = pnt.scrCoords;
508 
509             if (this._borderCase(curve, a, b, c, ta, tb, tc, depth)) {
510                 return this;
511             }
512 
513             ds = this._triangleDists(a, b, c);           // returns [d_ab, d_ac, d_cb, d_cd]
514 
515             isSmooth = (depth < this.smoothLevel) && (ds[3] < delta);
516 
517             isJump = (depth < this.jumpLevel) &&
518                         ((ds[2] > jump_threshold * ds[0]) ||
519                          (ds[1] > jump_threshold * ds[0]) ||
520                         ds[0] === Infinity || ds[1] === Infinity || ds[2] === Infinity);
521 
522             isCusp = (depth < this.smoothLevel + 2) && (ds[0] < cusp_threshold * (ds[1] + ds[2]));
523 
524             if (isCusp) {
525                 mindepth = 0;
526                 isSmooth = false;
527             }
528 
529             --depth;
530 
531             if (isJump) {
532                 this._insertPoint_v2(curve, new Coords(Const.COORDS_BY_SCREEN, [NaN, NaN], curve.board, false), tc);
533             } else if (depth <= mindepth || isSmooth) {
534                 this._insertPoint_v2(curve, pnt, tc);
535                 //if (this._borderCase(a, b, c, ta, tb, tc, depth)) {}
536             } else {
537                 this._plotRecursive_v2(curve, a, ta, c, tc, depth, delta);
538 
539                 if (!isNaN(pnt.scrCoords[1] + pnt.scrCoords[2])) {
540                     this._insertPoint_v2(curve, pnt, tc);
541                 }
542 
543                 this._plotRecursive_v2(curve, c, tc, b, tb, depth, delta);
544             }
545 
546             return this;
547         },
548 
549         /**
550          * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#plotVersion} is <tt>3</tt>.
551          *
552          * @param {JXG.Curve} curve JSXGraph curve element
553          * @param {Number} mi Left bound of curve
554          * @param {Number} ma Right bound of curve
555          * @returns {JXG.Curve} Reference to the curve object.
556          */
557         updateParametricCurve_v2: function (curve, mi, ma) {
558             var ta, tb, a, b,
559                 suspendUpdate = false,
560                 pa = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false),
561                 pb = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false),
562                 depth, delta,
563                 w2, h2, bbox,
564                 ret_arr;
565 
566             //console.time("plot");
567             if (curve.board.updateQuality === curve.board.BOARD_QUALITY_LOW) {
568                 depth = Type.evaluate(curve.visProp.recursiondepthlow) || 13;
569                 delta = 2;
570                 // this.smoothLevel = 5; //depth - 7;
571                 this.smoothLevel = depth - 6;
572                 this.jumpLevel = 3;
573             } else {
574                 depth = Type.evaluate(curve.visProp.recursiondepthhigh) || 17;
575                 delta = 2;
576                 // smoothLevel has to be small for graphs in a huge interval.
577                 // this.smoothLevel = 3; //depth - 7; // 9
578                 this.smoothLevel = depth - 9; // 9
579                 this.jumpLevel = 2;
580             }
581             this.nanLevel = depth - 4;
582 
583             curve.points = [];
584 
585             if (this.xterm === 'x') {
586                 // For function graphs we can restrict the plot interval
587                 // to the visible area +plus margin
588                 bbox = curve.board.getBoundingBox();
589                 w2 = (bbox[2] - bbox[0]) * 0.3;
590                 h2 = (bbox[1] - bbox[3]) * 0.3;
591                 ta = Math.max(mi, bbox[0] - w2);
592                 tb = Math.min(ma, bbox[2] + w2);
593             } else {
594                 ta = mi;
595                 tb = ma;
596             }
597             pa.setCoordinates(Const.COORDS_BY_USER, [curve.X(ta, suspendUpdate), curve.Y(ta, suspendUpdate)], false);
598 
599             // The first function calls of X() and Y() are done. We can now
600             // switch `suspendUpdate` on. If supported by the functions, this
601             // avoids for the rest of the plotting algorithm, evaluation of any
602             // parent elements.
603             suspendUpdate = true;
604 
605             pb.setCoordinates(Const.COORDS_BY_USER, [curve.X(tb, suspendUpdate), curve.Y(tb, suspendUpdate)], false);
606 
607             // Find start and end points of the visible area (plus a certain margin)
608             ret_arr = this._findStartPoint(curve, pa.scrCoords, ta, pb.scrCoords, tb);
609             pa.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false);
610             ta = ret_arr[1];
611             ret_arr = this._findStartPoint(curve, pb.scrCoords, tb, pa.scrCoords, ta);
612             pb.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false);
613             tb = ret_arr[1];
614 
615             // Save the visible area.
616             // This can be used in Curve.hasPoint().
617             this._visibleArea = [ta, tb];
618 
619             // Start recursive plotting algorithm
620             a = pa.copy('scrCoords');
621             b = pb.copy('scrCoords');
622             pa._t = ta;
623             curve.points.push(pa);
624             this._lastCrds = pa.copy('scrCoords');   // Used in _insertPoint
625             this._plotRecursive_v2(curve, a, ta, b, tb, depth, delta);
626             pb._t = tb;
627             curve.points.push(pb);
628 
629             curve.numberPoints = curve.points.length;
630             //console.timeEnd("plot");
631 
632             return curve;
633         },
634 
635         //----------------------------------------------------------------------
636         // Plot algorithm v3
637         //----------------------------------------------------------------------
638         /**
639          *
640          * @param {JXG.Curve} curve JSXGraph curve element
641          * @param {*} pnt
642          * @param {*} t
643          * @param {*} depth
644          * @param {*} limes
645          * @private
646          */
647         _insertLimesPoint: function(curve, pnt, t, depth, limes) {
648             var p0, p1, p2;
649 
650             // Ignore jump point if it follows limes
651             if ((Math.abs(this._lastUsrCrds[1]) === Infinity && Math.abs(limes.left_x) === Infinity) ||
652                 (Math.abs(this._lastUsrCrds[2]) === Infinity && Math.abs(limes.left_y) === Infinity)) {
653                 // console.log("SKIP:", pnt.usrCoords, this._lastUsrCrds, limes);
654                 return;
655             }
656 
657             // // Ignore jump left from limes
658             // if (Math.abs(limes.left_x) > 100 * Math.abs(this._lastUsrCrds[1])) {
659             //     x = Math.sign(limes.left_x) * Infinity;
660             // } else {
661             //     x = limes.left_x;
662             // }
663             // if (Math.abs(limes.left_y) > 100 * Math.abs(this._lastUsrCrds[2])) {
664             //     y = Math.sign(limes.left_y) * Infinity;
665             // } else {
666             //     y = limes.left_y;
667             // }
668             // //pnt.setCoordinates(Const.COORDS_BY_USER, [x, y], false);
669 
670             // Add points at a jump. pnt contains [NaN, NaN]
671             //console.log("Add", t, pnt.usrCoords, limes, depth)
672             p0 = new Coords(Const.COORDS_BY_USER, [limes.left_x, limes.left_y], curve.board);
673             p0._t = t;
674             curve.points.push(p0);
675 
676             if (!isNaN(limes.left_x) && !isNaN(limes.left_y) && !isNaN(limes.right_x) && !isNaN(limes.right_y) &&
677                 (Math.abs(limes.left_x - limes.right_x) > Mat.eps || Math.abs(limes.left_y - limes.right_y) > Mat.eps)) {
678                 p1 = new Coords(Const.COORDS_BY_SCREEN, pnt, curve.board);
679                 p1._t = t;
680                 curve.points.push(p1);
681             }
682 
683             p2 = new Coords(Const.COORDS_BY_USER, [limes.right_x, limes.right_y], curve.board);
684             p2._t = t;
685             curve.points.push(p2);
686             this._lastScrCrds = p2.copy('scrCoords');
687             this._lastUsrCrds = p2.copy('usrCoords');
688 
689         },
690 
691         /**
692          * Add a point to the curve plot. If the new point is too close to the previously inserted point,
693          * it is skipped.
694          * Used in {@link JXG.Curve._plotRecursive}.
695          *
696          * @private
697          * @param {JXG.Curve} curve JSXGraph curve element
698          * @param {JXG.Coords} pnt Coords to add to the list of points
699          */
700         _insertPoint: function (curve, pnt, t, depth, limes) {
701             var last_is_real = !isNaN(this._lastScrCrds[1] + this._lastScrCrds[2]),     // The last point was real
702                 point_is_real  = !isNaN(pnt[1] + pnt[2]),                               // New point is real point
703                 cw = curve.board.canvasWidth,
704                 ch = curve.board.canvasHeight,
705                 p,
706                 near = 0.8,
707                 off = 500;
708 
709             if (Type.exists(limes)) {
710                 this._insertLimesPoint(curve, pnt, t, depth, limes);
711                 return;
712             }
713 
714             // Check if point has real coordinates and
715             // coordinates are not too far away from canvas.
716             point_is_real = point_is_real &&
717                         (pnt[1] > -off     && pnt[2] > -off &&
718                          pnt[1] < cw + off && pnt[2] < ch + off);
719 
720             // Prevent two consecutive NaNs
721             if (!last_is_real && !point_is_real) {
722                 return;
723             }
724 
725             // Prevent two consecutive points which are too close
726             if (point_is_real && last_is_real &&
727                 Math.abs(pnt[1] - this._lastScrCrds[1]) < near &&
728                 Math.abs(pnt[2] - this._lastScrCrds[2]) < near) {
729                 return;
730             }
731 
732             // Prevent two consecutive points at infinity (either direction)
733             if ((Math.abs(pnt[1]) === Infinity &&
734                  Math.abs(this._lastUsrCrds[1]) === Infinity) ||
735                 (Math.abs(pnt[2]) === Infinity &&
736                  Math.abs(this._lastUsrCrds[2]) === Infinity)) {
737                 return;
738             }
739 
740             //console.log("add", t, pnt.usrCoords, depth)
741             // Add regular point
742             p = new Coords(Const.COORDS_BY_SCREEN, pnt, curve.board);
743             p._t = t;
744             curve.points.push(p);
745             this._lastScrCrds = p.copy('scrCoords');
746             this._lastUsrCrds = p.copy('usrCoords');
747         },
748 
749         /**
750          * Compute distances in screen coordinates between the points ab,
751          * ac, cb, and cd, where d = (a + b)/2.
752          * cd is used for the smoothness test, ab, ac, cb are used to detect jumps, cusps and poles.
753          *
754          * @private
755          * @param {Array} a Screen coordinates of the left interval bound
756          * @param {Array} b Screen coordinates of the right interval bound
757          * @param {Array} c Screen coordinates of the bisection point at (ta + tb) / 2
758          * @returns {Array} array of distances in screen coordinates between: ab, ac, cb, and cd.
759          */
760         _triangleDists: function (a, b, c) {
761             var d, d_ab, d_ac, d_cb, d_cd;
762 
763             d = [a[0] * b[0], (a[1] + b[1]) * 0.5, (a[2] + b[2]) * 0.5];
764 
765             d_ab = Geometry.distance(a, b, 3);
766             d_ac = Geometry.distance(a, c, 3);
767             d_cb = Geometry.distance(c, b, 3);
768             d_cd = Geometry.distance(c, d, 3);
769 
770             return [d_ab, d_ac, d_cb, d_cd];
771         },
772 
773         /**
774          * Test if the function is undefined on an interval:
775          * If the interval borders a and b are undefined, 20 random values
776          * are tested if they are undefined, too.
777          * Only if all values are undefined, we declare the function to be undefined in this interval.
778          *
779          * @private
780          * @param {JXG.Curve} curve JSXGraph curve element
781          * @param {Array} a Screen coordinates of the left interval bound
782          * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates
783          * @param {Array} b Screen coordinates of the right interval bound
784          * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates
785          */
786         _isUndefined: function (curve, a, ta, b, tb) {
787             var t, i, pnt;
788 
789             if (!isNaN(a[1] + a[2]) || !isNaN(b[1] + b[2])) {
790                 return false;
791             }
792 
793             pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false);
794 
795             for (i = 0; i < 20; ++i) {
796                 t = ta + Math.random() * (tb - ta);
797                 pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, true), curve.Y(t, true)], false);
798                 if (!isNaN(pnt.scrCoords[0] + pnt.scrCoords[1] + pnt.scrCoords[2])) {
799                     return false;
800                 }
801             }
802 
803             return true;
804         },
805 
806         /**
807          * Decide if a path segment is too far from the canvas that we do not need to draw it.
808          * @private
809          * @param  {Array}  a  Screen coordinates of the start point of the segment
810          * @param  {Array}  ta Curve parameter of a  (unused).
811          * @param  {Array}  b  Screen coordinates of the end point of the segment
812          * @param  {Array}  tb Curve parameter of b (unused).
813          * @param  {JXG.Board} board
814          * @returns {Boolean}   True if the segment is too far away from the canvas, false otherwise.
815          */
816         _isOutside: function (a, ta, b, tb, board) {
817             var off = 500,
818                 cw = board.canvasWidth,
819                 ch = board.canvasHeight;
820 
821             return !!((a[1] < -off && b[1] < -off) ||
822                 (a[2] < -off && b[2] < -off) ||
823                 (a[1] > cw + off && b[1] > cw + off) ||
824                 (a[2] > ch + off && b[2] > ch + off));
825         },
826 
827         /**
828          * Decide if a point of a curve is too far from the canvas that we do not need to draw it.
829          * @private
830          * @param {Array}  a  Screen coordinates of the point
831          * @param {JXG.Board} board
832          * @returns {Boolean}  True if the point is too far away from the canvas, false otherwise.
833          */
834         _isOutsidePoint: function (a, board) {
835             var off = 500,
836                 cw = board.canvasWidth,
837                 ch = board.canvasHeight;
838 
839             return !!(a[1] < -off ||
840                       a[2] < -off ||
841                       a[1] > cw + off ||
842                       a[2] > ch + off);
843         },
844 
845         /**
846          * For a curve c(t) defined on the interval [ta, tb] find the first point
847          * which is in the visible area of the board (plus some outside margin).
848          * <p>
849          * This method is necessary to restrict the recursive plotting algorithm
850          * {@link JXG.Curve._plotRecursive} to the visible area and not waste
851          * recursion to areas far outside of the visible area.
852          * <p>
853          * This method can also be used to find the last visible point
854          * by reversing the input parameters.
855          *
856          * @param {JXG.Curve} curve JSXGraph curve element
857          * @param  {Array}  ta Curve parameter of a.
858          * @param  {Array}  b  Screen coordinates of the end point of the segment (unused)
859          * @param  {Array}  tb Curve parameter of b
860          * @return {Array}  Array of length two containing the screen ccordinates of
861          * the starting point and the curve parameter at this point.
862          * @private
863          */
864         _findStartPoint: function (curve, a, ta, b, tb) {
865             var i, delta, tc,
866                 td, z, isFound,
867                 w2, h2,
868                 pnt =  new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false),
869                 steps = 40,
870                 eps = 0.01,
871                 fnX1, fnX2, fnY1, fnY2,
872                 bbox = curve.board.getBoundingBox();
873 
874             // The code below is too unstable.
875             // E.g. [function(t) { return Math.pow(t, 2) * (t + 5) * Math.pow(t - 5, 2); }, -8, 8]
876             // Therefore, we return here.
877             if (true || !this._isOutsidePoint(a, curve.board)) {
878                 return [a, ta];
879             }
880 
881             w2 = (bbox[2] - bbox[0]) * 0.3;
882             h2 = (bbox[1] - bbox[3]) * 0.3;
883             bbox[0] -= w2;
884             bbox[1] += h2;
885             bbox[2] += w2;
886             bbox[3] -= h2;
887 
888             delta = (tb - ta) / steps;
889             tc = ta + delta;
890             isFound = false;
891 
892             fnX1 = function(t) { return curve.X(t, true) - bbox[0]; };
893             fnY1 = function(t) { return curve.Y(t, true) - bbox[1]; };
894             fnX2 = function(t) { return curve.X(t, true) - bbox[2]; };
895             fnY2 = function(t) { return curve.Y(t, true) - bbox[3]; };
896             for (i = 0; i < steps; ++i) {
897                 // Left border
898                 z = bbox[0];
899                 td = Numerics.root(fnX1, [tc - delta, tc], curve);
900                 // td = Numerics.fzero(fnX1, [tc - delta, tc], this);
901                 // console.log("A", tc - delta, tc, td, Math.abs(this.X(td, true) - z));
902                 if (Math.abs(curve.X(td, true) - z) < eps) { //} * Math.abs(z)) {
903                     isFound = true;
904                     break;
905                 }
906                 // Top border
907                 z = bbox[1];
908                 td = Numerics.root(fnY1, [tc - delta, tc], curve);
909                 // td = Numerics.fzero(fnY1, [tc - delta, tc], this);
910                 // console.log("B", tc - delta, tc, td, Math.abs(this.Y(td, true) - z));
911                 if (Math.abs(curve.Y(td, true) - z) < eps) { // * Math.abs(z)) {
912                     isFound = true;
913                     break;
914                 }
915                 // Right border
916                 z = bbox[2];
917                 td = Numerics.root(fnX2, [tc - delta, tc], curve);
918                 // td = Numerics.fzero(fnX2, [tc - delta, tc], this);
919                 // console.log("C", tc - delta, tc, td, Math.abs(this.X(td, true) - z));
920                 if (Math.abs(curve.X(td, true) - z) < eps) { // * Math.abs(z)) {
921                     isFound = true;
922                     break;
923                 }
924                 // Bottom border
925                 z = bbox[3];
926                 td = Numerics.root(fnY2, [tc - delta, tc], curve);
927                 // td = Numerics.fzero(fnY2, [tc - delta, tc], this);
928                 // console.log("D", tc - delta, tc, td, Math.abs(this.Y(td, true) - z));
929                 if (Math.abs(curve.Y(td, true) - z) < eps) { // * Math.abs(z)) {
930                     isFound = true;
931                     break;
932                 }
933                 tc += delta;
934             }
935             if (isFound) {
936                 pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(td, true), curve.Y(td, true)], false);
937                 return [pnt.scrCoords, td];
938             }
939             console.log("TODO _findStartPoint", curve.Y.toString(), tc);
940             pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(ta, true), curve.Y(ta, true)], false);
941             return [pnt.scrCoords, ta];
942         },
943 
944         /**
945          * Investigate a function term at the bounds of intervals where
946          * the function is not defined, e.g. log(x) at x = 0.
947          *
948          * c is inbetween a and b
949          *
950          * @param {JXG.Curve} curve JSXGraph curve element
951          * @param {Array} a Screen coordinates of the left interval bound
952          * @param {Array} b Screen coordinates of the right interval bound
953          * @param {Array} c Screen coordinates of the bisection point at (ta + tb) / 2
954          * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates
955          * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates
956          * @param {Number} tc (ta + tb) / 2 = tc. Parameter which evaluates to b, i.e. [1, X(tc), Y(tc)] = c in screen coordinates
957          * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0.
958          * @returns {JXG.Boolean} true if the point is inserted and the recursion should stop, false otherwise.
959          *
960          * @private
961          */
962         _getBorderPos: function(curve, ta, a, tc, c, tb, b) {
963             var t, pnt, p,
964                 j,
965                 max_it = 30,
966                 is_undef = false,
967                 t_real2,
968                 t_good, t_bad;
969 
970             pnt = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false);
971             j = 0;
972             // Bisect a, b and c until the point t_real is inside of the definition interval
973             // and as close as possible at the boundary.
974             // t_real2 is the second closest point.
975             // There are four cases:
976             //  a  |  c  |  b
977             // ---------------
978             // inf | R   | R
979             // R   | R   | inf
980             // inf | inf | R
981             // R   | inf | inf
982             //
983             if (isNaN(a[1] + a[2]) && !isNaN(c[1] + c[2])) {
984                 t_bad = ta;
985                 t_good = tc;
986                 t_real2 = tb;
987             } else if (isNaN(b[1] + b[2]) && !isNaN(c[1] + c[2])) {
988                 t_bad = tb;
989                 t_good = tc;
990                 t_real2 = ta;
991             } else if (isNaN(c[1] + c[2]) && !isNaN(b[1] + b[2])) {
992                 t_bad = tc;
993                 t_good = tb;
994                 t_real2 = tb + (tb - tc);
995             } else if (isNaN(c[1] + c[2]) && !isNaN(a[1] + a[2])) {
996                 t_bad = tc;
997                 t_good = ta;
998                 t_real2 = ta - (tc - ta);
999             } else {
1000                 return false;
1001             }
1002             do {
1003                 t = 0.5 * (t_good + t_bad);
1004                 pnt.setCoordinates(Const.COORDS_BY_USER, [curve.X(t, true), curve.Y(t, true)], false);
1005                 p = pnt.usrCoords;
1006                 is_undef = isNaN(p[1] + p[2]);
1007                 if (is_undef) {
1008                     t_bad = t;
1009                 } else {
1010                     t_real2 = t_good;
1011                     t_good = t;
1012                 }
1013                 ++j;
1014             } while (j < max_it && Math.abs(t_good - t_bad) > Mat.eps);
1015             return t;
1016         },
1017 
1018         /**
1019          *
1020          * @param {JXG.Curve} curve JSXGraph curve element
1021          * @param {Number} ta
1022          * @param {Number} tb
1023          */
1024         _getCuspPos: function(curve, ta, tb) {
1025             var a = [curve.X(ta, true), curve.Y(ta, true)],
1026                 b = [curve.X(tb, true), curve.Y(tb, true)],
1027                 max_func = function(t) {
1028                     var c = [curve.X(t, true), curve.Y(t, true)];
1029                     return -(Math.sqrt((a[0] - c[0]) * (a[0] - c[0]) + (a[1] - c[1]) * (a[1] - c[1])) +
1030                             Math.sqrt((b[0] - c[0]) * (b[0] - c[0]) + (b[1] - c[1]) * (b[1] - c[1])));
1031                 };
1032 
1033             return Numerics.fminbr(max_func, [ta, tb], curve);
1034         },
1035 
1036         /**
1037          *
1038          * @param {JXG.Curve} curve JSXGraph curve element
1039          * @param {Number} ta
1040          * @param {Number} tb
1041          */
1042         _getJumpPos: function(curve, ta, tb) {
1043             var max_func = function(t) {
1044                     var e = Mat.eps * Mat.eps,
1045                         c1 = [curve.X(t, true), curve.Y(t, true)],
1046                         c2 = [curve.X(t + e, true), curve.Y(t + e, true)];
1047                     return -Math.abs( (c2[1] - c1[1]) / (c2[0] - c1[0]) );
1048                 };
1049 
1050             return Numerics.fminbr(max_func, [ta, tb], curve);
1051         },
1052 
1053         /**
1054          *
1055          * @param {JXG.Curve} curve JSXGraph curve element
1056          * @param {Number} t
1057          * @private
1058          */
1059         _getLimits: function(curve, t) {
1060             var res,
1061                 step = 2 / (curve.maxX() - curve.minX()),
1062                 x_l, x_r, y_l, y_r;
1063 
1064             // From left
1065             res = Extrapolate.limit(t, -step, curve.X);
1066             x_l = res[0];
1067             if (res[1] === 'infinite') {
1068                 x_l = Math.sign(x_l) * Infinity;
1069             }
1070 
1071             res = Extrapolate.limit(t, -step, curve.Y);
1072             y_l = res[0];
1073             if (res[1] === 'infinite') {
1074                 y_l = Math.sign(y_l) * Infinity;
1075             }
1076 
1077             // From right
1078             res = Extrapolate.limit(t, step, curve.X);
1079             x_r = res[0];
1080             if (res[1] === 'infinite') {
1081                 x_r = Math.sign(x_r) * Infinity;
1082             }
1083 
1084             res = Extrapolate.limit(t, step, curve.Y);
1085             y_r = res[0];
1086             if (res[1] === 'infinite') {
1087                 y_r = Math.sign(y_r) * Infinity;
1088             }
1089 
1090             return {
1091                     left_x: x_l,
1092                     left_y: y_l,
1093                     right_x: x_r,
1094                     right_y: y_r,
1095                     t: t
1096                 };
1097         },
1098 
1099         /**
1100          *
1101          * @param {JXG.Curve} curve JSXGraph curve element
1102          * @param {Array} a
1103          * @param {Number} tc
1104          * @param {Array} c
1105          * @param {Number} tb
1106          * @param {Array} b
1107          * @param {String} may_be_special
1108          * @param {Number} depth
1109          * @private
1110          */
1111         _getLimes: function(curve, ta, a, tc, c, tb, b, may_be_special, depth) {
1112             var t;
1113 
1114             if (may_be_special === 'border') {
1115                 t = this._getBorderPos(curve, ta, a, tc, c, tb, b);
1116             } else if (may_be_special === 'cusp') {
1117                 t = this._getCuspPos(curve, ta, tb);
1118             } else if (may_be_special === 'jump') {
1119                 t = this._getJumpPos(curve, ta, tb);
1120             }
1121             return this._getLimits(curve, t);
1122         },
1123 
1124         /**
1125          * Recursive interval bisection algorithm for curve plotting.
1126          * Used in {@link JXG.Curve.updateParametricCurve}.
1127          * @private
1128          * @param {JXG.Curve} curve JSXGraph curve element
1129          * @param {Array} a Screen coordinates of the left interval bound
1130          * @param {Number} ta Parameter which evaluates to a, i.e. [1, X(ta), Y(ta)] = a in screen coordinates
1131          * @param {Array} b Screen coordinates of the right interval bound
1132          * @param {Number} tb Parameter which evaluates to b, i.e. [1, X(tb), Y(tb)] = b in screen coordinates
1133          * @param {Number} depth Actual recursion depth. The recursion stops if depth is equal to 0.
1134          * @param {Number} delta If the distance of the bisection point at (ta + tb) / 2 from the point (a + b) / 2 is less then delta,
1135          *                 the segment [a,b] is regarded as straight line.
1136          * @returns {JXG.Curve} Reference to the curve object.
1137          */
1138         _plotNonRecursive: function (curve, a, ta, b, tb, d) {
1139             var tc, c, ds,
1140                 mindepth = 0,
1141                 limes = null,
1142                 a_nan, b_nan,
1143                 isSmooth = false,
1144                 may_be_special = '',
1145                 x, y, oc, depth, ds0,
1146                 stack = [],
1147                 stack_length = 0,
1148                 item;
1149 
1150             oc = curve.board.origin.scrCoords;
1151             stack[stack_length++] = [a, ta, b, tb, d, Infinity];
1152             while (stack_length > 0) {
1153                 // item = stack.pop();
1154                 item = stack[--stack_length];
1155                 a = item[0];
1156                 ta = item[1];
1157                 b = item[2];
1158                 tb = item[3];
1159                 depth = item[4];
1160                 ds0 = item[5];
1161 
1162                 isSmooth = false;
1163                 may_be_special = '';
1164                 limes = null;
1165                 //console.log(stack.length, item)
1166 
1167                 if (curve.points.length > 65536) {
1168                     return;
1169                 }
1170 
1171                 if (depth < this.nanLevel) {
1172                     // Test if the function is undefined in the whole interval [ta, tb]
1173                     if (this._isUndefined(curve, a, ta, b, tb)) {
1174                         continue;
1175                     }
1176                     // Test if the graph is far outside the visible are for the interval [ta, tb]
1177                     if (this._isOutside(a, ta, b, tb, curve.board)) {
1178                         continue;
1179                     }
1180                 }
1181 
1182                 tc = (ta  + tb) * 0.5;
1183 
1184                 // Screen coordinates of point at tc
1185                 x = curve.X(tc, true);
1186                 y = curve.Y(tc, true);
1187                 c = [1, oc[1] + x * curve.board.unitX, oc[2] - y * curve.board.unitY];
1188                 ds = this._triangleDists(a, b, c);           // returns [d_ab, d_ac, d_cb, d_cd]
1189 
1190                 a_nan = isNaN(a[1] + a[2]);
1191                 b_nan = isNaN(b[1] + b[2]);
1192                 if ((a_nan && !b_nan) || (!a_nan && b_nan)) {
1193                     may_be_special = 'border';
1194                 } else if (ds[0] > 0.66 * ds0 ||
1195                             ds[0] < this.cusp_threshold * (ds[1] + ds[2]) ||
1196                             ds[1] > 5 * ds[2] ||
1197                             ds[2] > 5 * ds[1]) {
1198                     may_be_special = 'cusp';
1199                 } else if ((ds[2] > this.jump_threshold * ds[0]) ||
1200                            (ds[1] > this.jump_threshold * ds[0]) ||
1201                             ds[0] === Infinity || ds[1] === Infinity || ds[2] === Infinity) {
1202                     may_be_special = 'jump';
1203                 }
1204                 isSmooth = (may_be_special === '' && depth < this.smoothLevel && ds[3] < this.smooth_threshold);
1205 
1206                 if (depth < this.testLevel && !isSmooth) {
1207                     if (may_be_special === '') {
1208                         isSmooth = true;
1209                     } else {
1210                         limes = this._getLimes(curve, ta, a, tc, c, tb, b, may_be_special, depth);
1211                     }
1212                 }
1213 
1214                 if (limes !== null) {
1215                     c = [1, NaN, NaN];
1216                     this._insertPoint(curve, c, tc, depth, limes);
1217                 } else if (depth <= mindepth || isSmooth) {
1218                     this._insertPoint(curve, c, tc, depth, null);
1219                 } else {
1220                     stack[stack_length++] = [c, tc, b, tb, depth - 1, ds[0]];
1221                     stack[stack_length++] = [a, ta, c, tc, depth - 1, ds[0]];
1222                 }
1223             }
1224 
1225             return this;
1226         },
1227 
1228         /**
1229          * Updates the data points of a parametric curve. This version is used if {@link JXG.Curve#plotVersion} is <tt>3</tt>.
1230          * This is an experimental plot version, <b>not recommended</b> to be used.
1231          * @param {JXG.Curve} curve JSXGraph curve element
1232          * @param {Number} mi Left bound of curve
1233          * @param {Number} ma Right bound of curve
1234          * @returns {JXG.Curve} Reference to the curve object.
1235          */
1236         updateParametricCurve_v3: function (curve, mi, ma) {
1237             var ta, tb, a, b,
1238                 suspendUpdate = false,
1239                 pa = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false),
1240                 pb = new Coords(Const.COORDS_BY_USER, [0, 0], curve.board, false),
1241                 depth,
1242                 w2, // h2,
1243                 bbox,
1244                 ret_arr;
1245 
1246             // console.log("-----------------------------------------------------------");
1247             // console.time("plot");
1248             if (curve.board.updateQuality === curve.board.BOARD_QUALITY_LOW) {
1249                 depth = Type.evaluate(curve.visProp.recursiondepthlow) || 14;
1250             } else {
1251                 depth = Type.evaluate(curve.visProp.recursiondepthhigh) || 17;
1252             }
1253 
1254             // smoothLevel has to be small for graphs in a huge interval.
1255             this.smoothLevel = 7; //depth - 10;
1256             this.nanLevel = depth - 4;
1257             this.testLevel = 4;
1258             this.cusp_threshold = 0.5;
1259             this.jump_threshold = 0.99;
1260             this.smooth_threshold = 2;
1261 
1262             curve.points = [];
1263 
1264             if (curve.xterm === 'x') {
1265                 // For function graphs we can restrict the plot interval
1266                 // to the visible area +plus margin
1267                 bbox = curve.board.getBoundingBox();
1268                 w2 = (bbox[2] - bbox[0]) * 0.3;
1269                 //h2 = (bbox[1] - bbox[3]) * 0.3;
1270                 ta = Math.max(mi, bbox[0] - w2);
1271                 tb = Math.min(ma, bbox[2] + w2);
1272             } else {
1273                 ta = mi;
1274                 tb = ma;
1275             }
1276             pa.setCoordinates(Const.COORDS_BY_USER, [curve.X(ta, suspendUpdate), curve.Y(ta, suspendUpdate)], false);
1277 
1278             // The first function calls of X() and Y() are done. We can now
1279             // switch `suspendUpdate` on. If supported by the functions, this
1280             // avoids for the rest of the plotting algorithm, evaluation of any
1281             // parent elements.
1282             suspendUpdate = true;
1283 
1284             pb.setCoordinates(Const.COORDS_BY_USER, [curve.X(tb, suspendUpdate), curve.Y(tb, suspendUpdate)], false);
1285 
1286             // Find start and end points of the visible area (plus a certain margin)
1287             ret_arr = this._findStartPoint(curve, pa.scrCoords, ta, pb.scrCoords, tb);
1288             pa.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false);
1289             ta = ret_arr[1];
1290             ret_arr = this._findStartPoint(curve, pb.scrCoords, tb, pa.scrCoords, ta);
1291             pb.setCoordinates(Const.COORDS_BY_SCREEN, ret_arr[0], false);
1292             tb = ret_arr[1];
1293 
1294             // Save the visible area.
1295             // This can be used in Curve.hasPoint().
1296             this._visibleArea = [ta, tb];
1297 
1298             // Start recursive plotting algorithm
1299             a = pa.copy('scrCoords');
1300             b = pb.copy('scrCoords');
1301             pa._t = ta;
1302             curve.points.push(pa);
1303             this._lastScrCrds = pa.copy('scrCoords');   // Used in _insertPoint
1304             this._lastUsrCrds = pa.copy('usrCoords');   // Used in _insertPoint
1305 
1306             this._plotNonRecursive(curve, a, ta, b, tb, depth);
1307 
1308             pb._t = tb;
1309             curve.points.push(pb);
1310 
1311             curve.numberPoints = curve.points.length;
1312             // console.timeEnd("plot");
1313             // console.log("number of points:", this.numberPoints);
1314 
1315             return curve;
1316         },
1317 
1318         //----------------------------------------------------------------------
1319         // Plot algorithm v4
1320         //----------------------------------------------------------------------
1321 
1322         _criticalInterval: function(vec, le, level) {
1323             var i, j, le1, med,
1324                 sgn, sgnChange,
1325                 isGroup   = false,
1326                 abs_vec,
1327                 last = -Infinity,
1328                 very_small = false,
1329                 smooth    = false,
1330                 group     = 0,
1331                 groups    = [],
1332                 types     = [],
1333                 positions = [];
1334 
1335             abs_vec = Statistics.abs(vec);
1336             med = Statistics.median(abs_vec);
1337 
1338             if (med < 1.0e-7) {
1339                 med = 1.0e-7;
1340                 very_small = true;
1341             } else {
1342                 med *= this.criticalThreshold;
1343             }
1344 
1345             //console.log("Median", med);
1346             for (i = 0; i < le; i++) {
1347                 // Start a group if not yet done and
1348                 // add position to group
1349                 if (abs_vec[i] > med /*&& abs_vec[i] > 0.01*/)  {
1350                     positions.push({i: i, v: vec[i], group: group});
1351                     last = i;
1352                     if (!isGroup) {
1353                         isGroup = true;
1354                     }
1355                 } else {
1356                     if (isGroup && i > last + 4) {
1357                         // End the group
1358                         if (positions.length > 0) {
1359                             groups.push(positions.slice(0));
1360                         }
1361                         positions = [];
1362                         isGroup = false;
1363                         group++;
1364                     }
1365                 }
1366             }
1367             if (isGroup) {
1368                 if (positions.length > 1) {
1369                     groups.push(positions.slice(0));
1370                 }
1371             }
1372 
1373             if (very_small && groups.length === 0) {
1374                 smooth = true;
1375             }
1376 
1377             // Decide if there is a singular critical point
1378             // or if a whole interval is problematic.
1379             // The latter is the case if the differences have many sign changes.
1380             for (j = 0; j < groups.length; j++) {
1381                 types[j] = 'point';
1382                 le1 = groups[j].length;
1383                 if (le1 < 64) {
1384                     continue;
1385                 }
1386                 sgnChange = 0;
1387                 sgn = Math.sign(groups[j][0].v);
1388                 for (i = 1; i < le1; i++) {
1389                     if (Math.sign(groups[j][i].v) !== sgn) {
1390                         sgnChange++;
1391                         sgn = Math.sign(groups[j][i].v);
1392                     }
1393                 }
1394                 if (sgnChange * 6 > le1) {
1395                     types[j] = 'interval';
1396                 }
1397             }
1398 
1399             return {smooth: smooth, groups: groups, types: types};
1400         },
1401 
1402         Component: function() {
1403             this.left_isNaN =  false;
1404             this.right_isNaN = false;
1405             this.left_t = null;
1406             this.right_t = null;
1407             this.t_values = [];
1408             this.x_values = [];
1409             this.y_values = [];
1410             this.len = 0;
1411         },
1412 
1413         findComponents: function(curve, mi, ma, steps) {
1414             var i, t, le, h, x, y,
1415                 components = [],
1416                 comp,
1417                 comp_nr = 0,
1418                 cnt = 0,
1419                 cntNaNs = 0,
1420                 comp_started = false,
1421                 suspended = false;
1422 
1423             h = (ma - mi) / steps;
1424             components[comp_nr] = new this.Component();
1425             comp = components[comp_nr];
1426 
1427             for (i = 0, t = mi; i <= steps; i++, t += h) {
1428                 x = curve.X(t, suspended);
1429                 y = curve.Y(t, suspended);
1430 
1431                 if (isNaN(x) || isNaN(y)) {
1432                     cntNaNs++;
1433                     // Wait for - at least - two consecutive NaNs
1434                     // This avoids starting a new component if
1435                     // the function value has infinity as intermediate value.
1436                     if (cntNaNs > 1 && comp_started) {
1437                         // Finalize a component
1438                         comp.right_isNaN = true;
1439                         comp.right_t = t - h;
1440                         comp.len = cnt;
1441 
1442                         // Prepare a new component
1443                         comp_started = false;
1444                         comp_nr++;
1445                         components[comp_nr] =  new this.Component();
1446                         comp = components[comp_nr];
1447                         cntNaNs = 0;
1448                     }
1449                 } else {
1450                     // Now there is a non-NaN entry.
1451                     if (!comp_started) {
1452                         // Start the component
1453                         comp_started = true;
1454                         cnt = 0;
1455                         if (cntNaNs > 0) {
1456                             comp.left_t = t - h;
1457                             comp.left_isNaN = true;
1458                         }
1459                     }
1460                     cntNaNs = 0;
1461                     // Add the value to the component
1462                     comp.t_values[cnt] = t;
1463                     comp.x_values[cnt] = x;
1464                     comp.y_values[cnt] = y;
1465                     cnt++;
1466                 }
1467                 if (i === 0) {
1468                     suspended = true;
1469                 }
1470             }
1471             if (comp_started) {
1472                 comp.len = cnt;
1473             } else {
1474                 components.pop();
1475             }
1476 
1477             return components;
1478         },
1479 
1480         getPointType: function(curve, pos, t_approx, t_values, x_table, y_table, len) {
1481             var x_values = x_table[0],
1482                 y_values = y_table[0],
1483                 full_len = t_values.length,
1484                 result = {
1485                     idx: pos,
1486                     t: t_approx, //t_values[pos],
1487                     x: x_values[pos],
1488                     y: y_values[pos],
1489                     type: 'other'
1490                 };
1491 
1492             if (pos < 5) {
1493                 result.type = 'borderleft';
1494                 result.idx = 0;
1495                 result.t = t_values[0];
1496                 result.x = x_values[0];
1497                 result.y = y_values[0];
1498 
1499                 // console.log('Border left', result.t);
1500                 return result;
1501             }
1502             if (pos > len - 6) {
1503                 result.type = 'borderright';
1504                 result.idx = full_len - 1;
1505                 result.t = t_values[full_len - 1];
1506                 result.x = x_values[full_len - 1];
1507                 result.y = y_values[full_len - 1];
1508 
1509                 // console.log('Border right', result.t, full_len - 1);
1510                 return result;
1511             }
1512 
1513             return result;
1514         },
1515 
1516         newtonApprox: function(idx, t, h, level, table) {
1517             var i, s = 0.0;
1518             for (i = level; i > 0; i--) {
1519                 s = (s + table[i][idx]) * (t - (i - 1) * h) / i;
1520             }
1521             return s + table[0][idx];
1522         },
1523 
1524         thiele: function(t, recip, t_values, idx, degree) {
1525             var i, v = 0.0;
1526             for (i = degree; i > 1; i--) {
1527                 v = (t - t_values[idx + i]) / (recip[i][idx + 1] - recip[i - 2][idx + 1] + v);
1528             }
1529             return recip[0][idx + 1] + (t - t_values[idx + 1]) / (recip[1][idx + 1] + v);
1530         },
1531 
1532         differenceMethodExperiments: function(component, curve) {
1533             var i, level, le, up,
1534                 t_values = component.t_values,
1535                 x_values = component.x_values,
1536                 y_values = component.y_values,
1537                 x_diffs = [],
1538                 y_diffs = [],
1539                 x_slopes = [],
1540                 y_slopes = [],
1541                 x_table = [],
1542                 y_table = [],
1543                 x_recip = [],
1544                 y_recip = [],
1545                 h, numerator,
1546                 // x_med, y_med,
1547                 foundCriticalPoint = 0,
1548                 pos, ma, j, v,
1549                 groups,
1550                 criticalPoints = [];
1551 
1552             h = t_values[1] - t_values[0];
1553             x_table.push([]);
1554             y_table.push([]);
1555             x_recip.push([]);
1556             y_recip.push([]);
1557             le = y_values.length;
1558             for (i = 0; i < le; i++) {
1559                 x_table[0][i] = x_values[i];
1560                 y_table[0][i] = y_values[i];
1561                 x_recip[0][i] = x_values[i];
1562                 y_recip[0][i] = y_values[i];
1563             }
1564 
1565             x_table.push([]);
1566             y_table.push([]);
1567             x_recip.push([]);
1568             y_recip.push([]);
1569             numerator = h;
1570             le = y_values.length - 1;
1571             for (i = 0; i < le; i++) {
1572                 x_diffs[i] = x_values[i + 1] - x_values[i];
1573                 y_diffs[i] = y_values[i + 1] - y_values[i];
1574                 x_slopes[i] = x_diffs[i];
1575                 y_slopes[i] = y_diffs[i];
1576                 x_table[1][i] = x_diffs[i];
1577                 y_table[1][i] = y_diffs[i];
1578                 x_recip[1][i] = numerator / x_diffs[i];
1579                 y_recip[1][i] = numerator / y_diffs[i];
1580             }
1581             le--;
1582 
1583             up = Math.min(8, y_values.length - 1);
1584             for (level = 1; level < up; level++) {
1585                 x_table.push([]);
1586                 y_table.push([]);
1587                 x_recip.push([]);
1588                 y_recip.push([]);
1589                 numerator *= h;
1590                 for (i = 0; i < le; i++) {
1591                     x_diffs[i] = x_diffs[i + 1] - x_diffs[i];
1592                     y_diffs[i] = y_diffs[i + 1] - y_diffs[i];
1593                     x_table[level + 1][i] = x_diffs[i];
1594                     y_table[level + 1][i] = y_diffs[i];
1595                     x_recip[level + 1][i] = numerator / (x_recip[level][i + 1] - x_recip[level][i]) + x_recip[level - 1][i + 1];
1596                     y_recip[level + 1][i] = numerator / (y_recip[level][i + 1] - y_recip[level][i]) + y_recip[level - 1][i + 1];
1597                 }
1598 
1599                 // if (level == 1) {
1600                 //     console.log("bends level=", level, y_diffs.toString());
1601                 // }
1602 
1603                 // Store point location which may be centered around
1604                 // critical points.
1605                 // If the lebvel is suitable, step out of the loop.
1606                 groups = this._criticalPoints(y_diffs, le, level);
1607                 if (groups === false) {
1608                     // Its seems, the degree of the polynomial is equal to level
1609 console.log("Polynomial of degree", level);
1610                     groups = [];
1611                     break;
1612                 }
1613                 if (groups.length > 0) {
1614                     foundCriticalPoint++;
1615                     if (foundCriticalPoint > 1 && level % 2 === 0) {
1616                         break;
1617                     }
1618                 }
1619                 le--;
1620             }
1621 
1622             // console.log("Last diffs", y_diffs, "level", level);
1623 
1624             // Analyze the groups which have been found.
1625             for (i = 0; i < groups.length; i++) {
1626                 // console.log("Group", i, groups[i])
1627                 // Identify the maximum difference, i.e. the center of the "problem"
1628                 ma = -Infinity;
1629                 for (j = 0; j < groups[i].length; j++) {
1630                     v = Math.abs(groups[i][j].v);
1631                     if (v > ma) {
1632                         ma = v;
1633                         pos = j;
1634                     }
1635                 }
1636                 pos = Math.floor(groups[i][pos].i + level / 2);
1637                 // Analyze the critical point
1638                 criticalPoints.push(this.getPointType(curve, pos, t_values, x_values, y_values, x_slopes, y_slopes, le + 1));
1639             }
1640 
1641             return [criticalPoints, x_table, y_table, x_recip, y_recip];
1642 
1643         },
1644 
1645         getCenterOfCriticalInterval: function(group, degree, t_values) {
1646             var ma, j, pos, v,
1647                 num = 0.0,
1648                 den = 0.0,
1649                 h = t_values[1] - t_values[0],
1650                 pos_mean,
1651                 range = [];
1652 
1653             // Identify the maximum difference, i.e. the center of the "problem"
1654             // If there are several equal maxima, store the positions
1655             // in the array range and determine the center of the array.
1656 
1657             ma = -Infinity;
1658             range = [];
1659             for (j = 0; j < group.length; j++) {
1660                 v = Math.abs(group[j].v);
1661                 if (v > ma) {
1662                     range = [j];
1663                     ma = v;
1664                     pos = j;
1665                 } else if (ma === v) {
1666                     range.push(j);
1667                 }
1668             }
1669             if (range.length > 0) {
1670                 pos_mean = range.reduce(function(total, val) { return total + val; }, 0) / range.length;
1671                 pos = Math.floor(pos_mean);
1672                 pos_mean += group[0].i;
1673             }
1674 
1675             if (ma < Infinity) {
1676                 for (j = 0; j < group.length; j++) {
1677                     num += Math.abs(group[j].v) * group[j].i;
1678                     den += Math.abs(group[j].v);
1679                 }
1680                 pos_mean = num / den;
1681             }
1682             pos_mean += degree / 2;
1683             return [group[pos].i + degree / 2, pos_mean, t_values[Math.floor(pos_mean)] + h * (pos_mean - Math.floor(pos_mean))];
1684         },
1685 
1686         differenceMethod: function(component, curve) {
1687             var i, level, le, up,
1688                 t_values = component.t_values,
1689                 x_values = component.x_values,
1690                 y_values = component.y_values,
1691                 x_table = [],
1692                 y_table = [],
1693                 foundCriticalPoint = 0,
1694                 degree_x = -1,
1695                 degree_y = -1,
1696                 pos, res, res_x, res_y, t_approx,
1697                 groups = [],
1698                 types,
1699                 criticalPoints = [];
1700 
1701             le = y_values.length;
1702             // x_table.push([]);
1703             // y_table.push([]);
1704             // for (i = 0; i < le; i++) {
1705             //     x_table[0][i] = x_values[i];
1706             //     y_table[0][i] = y_values[i];
1707             // }
1708             x_table.push(new Float64Array(x_values));
1709             y_table.push(new Float64Array(y_values));
1710 
1711             le--;
1712             up = Math.min(12, le);
1713             for (level = 0; level < up; level++) {
1714                 // Old style method:
1715                 // x_table.push([]);
1716                 // y_table.push([]);
1717                 // for (i = 0; i < le; i++) {
1718                 //     x_table[level + 1][i] = x_table[level][i + 1] - x_table[level][i];
1719                 //     y_table[level + 1][i] = y_table[level][i + 1] - y_table[level][i];
1720                 // }
1721                 // New method:
1722                 x_table.push(new Float64Array(le));
1723                 y_table.push(new Float64Array(le));
1724                 x_table[level + 1] = x_table[level].map(function(v, idx, arr) { return arr[idx + 1] - v;});
1725                 y_table[level + 1] = y_table[level].map(function(v, idx, arr) { return arr[idx + 1] - v;});
1726 
1727                 // Store point location which may be centered around critical points.
1728                 // If the level is suitable, step out of the loop.
1729                 res_y = this._criticalInterval(y_table[level + 1], le, level);
1730                 if (res_y.smooth === true) {
1731                     // Its seems, the degree of the polynomial is equal to level
1732                     // If the values in level + 1 are zero, it might be a polynomial of degree level.
1733                     // Seems to work numerically stable until degree 6.
1734                     degree_y = level;
1735                     groups = [];
1736                 }
1737                 res_x = this._criticalInterval(x_table[level + 1], le, level);
1738                 if (degree_x === -1 && res_x.smooth === true) {
1739                     // Its seems, the degree of the polynomial is equal to level
1740                     // If the values in level + 1 are zero, it might be a polynomial of degree level.
1741                     // Seems to work numerically stable until degree 6.
1742                     degree_x = level;
1743                 }
1744                 if (degree_y >= 0) {
1745                     break;
1746                 }
1747 
1748                 if (res_y.groups.length > 0) {
1749                     foundCriticalPoint++;
1750                     if (foundCriticalPoint > 2 && (level + 1) % 2 === 0) {
1751                         groups = res_y.groups;
1752                         types = res_y.types;
1753                         break;
1754                     }
1755                 }
1756                 le--;
1757             }
1758 
1759             // console.log("Last diffs", y_table[Math.min(level + 1, up)], "level", level + 1);
1760             // Analyze the groups which have been found.
1761             for (i = 0; i < groups.length; i++) {
1762                 if (types[i] === 'interval') {
1763                     continue;
1764                 }
1765                 // console.log("Group", i, groups[i], types[i], level + 1)
1766                 res = this.getCenterOfCriticalInterval(groups[i], level + 1, t_values);
1767                 pos = res_y[0];
1768                 pos = Math.floor(res[1]);
1769                 t_approx = res[2];
1770                 // console.log("Critical points:", groups, res, pos)
1771 
1772                 // Analyze the type of the critical point
1773                 // Result is of type 'borderleft', borderright', 'other'
1774                 criticalPoints.push(this.getPointType(curve, pos, t_approx, t_values, x_table, y_table, le + 1));
1775             }
1776 
1777             // if (level === up) {
1778             //     console.log("No convergence!");
1779             // } else {
1780             //     console.log("Convergence level", level);
1781             // }
1782             return [criticalPoints, x_table, y_table, degree_x, degree_y];
1783 
1784         },
1785 
1786         _insertPoint_v4: function (curve, crds, t, doLog) {
1787             var p,
1788                 prev = null,
1789                 x, y,
1790                 near = 0.8;
1791 
1792             if (curve.points.length > 0) {
1793                 prev = curve.points[curve.points.length - 1].scrCoords;
1794             }
1795 
1796             // Add regular point
1797             p = new Coords(Const.COORDS_BY_USER, crds, curve.board);
1798 
1799             if (prev !== null) {
1800                 x = p.scrCoords[1] - prev[1];
1801                 y = p.scrCoords[2] - prev[2];
1802                 if (x * x + y * y < near * near) {
1803                 // Math.abs(p.scrCoords[1] - prev[1]) < near &&
1804                 // Math.abs(p.scrCoords[2] - prev[2]) < near) {
1805                     return;
1806                 }
1807             }
1808 
1809             p._t = t;
1810             curve.points.push(p);
1811         },
1812 
1813         getInterval: function(curve, ta, tb) {
1814             var t_int, x_int, y_int;
1815 
1816             //console.log('critical point', ta, tb);
1817             IntervalArithmetic.disable();
1818 
1819             t_int = IntervalArithmetic.Interval(ta, tb);
1820             curve.board.mathLib = IntervalArithmetic;
1821             curve.board.mathLibJXG = IntervalArithmetic;
1822             x_int = curve.X(t_int, true);
1823             y_int = curve.Y(t_int, true);
1824             curve.board.mathLib = Math;
1825             curve.board.mathLibJXG = JXG.Math;
1826 
1827             //console.log(x_int, y_int);
1828             return y_int;
1829         },
1830 
1831         sign: function (v) {
1832             if (v < 0) { return -1; }
1833             if (v > 0) { return 1; }
1834             return 0;
1835         },
1836 
1837         handleBorder: function(curve, comp, group, x_table, y_table) {
1838             var idx = group.idx,
1839                 t, t1, t2,
1840                 size = 32,
1841                 y_int,
1842                 x, y,
1843                 lo, hi, i,
1844                 components2, le, h;
1845 
1846             // console.log("HandleBorder at t =", t_approx);
1847             // console.log("component:", comp)
1848             // console.log("Group:", group);
1849 
1850             h = comp.t_values[1] - comp.t_values[0];
1851             if (group.type === 'borderleft') {
1852                 t = comp.left_isNaN ? comp.left_t : group.t - h;
1853                 t1 = t;
1854                 t2 = t1 + h;
1855             } else if (group.type === 'borderright') {
1856                 t = comp.right_isNaN ? comp.right_t : group.t + h;
1857                 t2 = t;
1858                 t1 = t2 - h;
1859             } else {
1860                 console.log("No bordercase!!!");
1861             }
1862 
1863             components2 = this.findComponents(curve, t1, t2, size);
1864             if (components2.length === 0) {
1865                 return;
1866             }
1867             if (group.type === 'borderleft') {
1868                 t1 = components2[0].left_t;
1869                 t2 = components2[0].t_values[0];
1870                 h = components2[0].t_values[1] - components2[0].t_values[0];
1871                 t1 = (t1 === null) ? t2- h : t1;
1872                 t = t1;
1873                 y_int = this.getInterval(curve, t1, t2);
1874                 if (Type.isObject(y_int)) {
1875                     lo = y_int.lo;
1876                     hi = y_int.hi;
1877 
1878                     x = curve.X(t, true);
1879                     y = (y_table[1][idx] < 0) ? hi : lo;
1880                     this._insertPoint_v4(curve, [1, x, y], t);
1881                 }
1882             }
1883 
1884             le = components2[0].t_values.length;
1885             for (i = 0; i < le; i++) {
1886                 t = components2[0].t_values[i];
1887                 x = components2[0].x_values[i];
1888                 y = components2[0].y_values[i];
1889                 this._insertPoint_v4(curve, [1, x, y], t);
1890             }
1891 
1892             if (group.type === 'borderright') {
1893                 t1 = components2[0].t_values[le - 1];
1894                 t2 = components2[0].right_t;
1895                 h = components2[0].t_values[1] - components2[0].t_values[0];
1896                 t2 = (t2 === null) ? t1 + h : t2;
1897 
1898                 t = t2;
1899                 y_int = this.getInterval(curve, t1, t2);
1900                 if (Type.isObject(y_int)) {
1901                     lo = y_int.lo;
1902                     hi = y_int.hi;
1903                     x = curve.X(t, true);
1904                     y = (y_table[1][idx] > 0) ? hi : lo;
1905                     this._insertPoint_v4(curve, [1, x, y], t);
1906                 }
1907             }
1908 
1909         },
1910 
1911         _seconditeration_v4: function(curve, comp, group, x_table, y_table) {
1912             var i, t1, t2, ret,
1913                 components2, comp2, idx, groups2, g,
1914                 x_table2, y_table2, start, le;
1915 
1916             // Look at two points, hopefully left and right from the critical point
1917             t1 = comp.t_values[group.idx - 2];
1918             t2 = comp.t_values[group.idx + 2];
1919             components2 = this.findComponents(curve, t1, t2, 64);
1920             for (idx = 0; idx < components2.length; idx++) {
1921                 comp2 = components2[idx];
1922                 ret = this.differenceMethod(comp2, curve);
1923                 groups2 = ret[0];
1924                 x_table2 = ret[1];
1925                 y_table2 = ret[2];
1926                 start = 0;
1927                 for (g = 0; g <= groups2.length; g++) {
1928                     if (g === groups2.length) {
1929                         le = comp2.len;
1930                     } else {
1931                         le = groups2[g].idx;
1932                     }
1933 
1934                     // Insert all uncritical points until next critical point
1935                     for (i = start; i < le; i++) {
1936                         if (!isNaN(comp2.x_values[i]) && !isNaN(comp2.y_values[i])) {
1937                             this._insertPoint_v4(curve, [1, comp2.x_values[i], comp2.y_values[i]], comp2.t_values[i]);
1938                         }
1939                     }
1940                     // Handle next critical point
1941                     if (g < groups2.length) {
1942                         this.handleSingularity(curve, comp2, groups2[g], x_table2, y_table2);
1943                         start = groups2[g].idx + 1;
1944                     }
1945                 }
1946                 le = comp2.len;
1947                 if (idx < components2.length - 1) {
1948                     this._insertPoint_v4(curve, [1, NaN, NaN], comp2.right_t);
1949                 }
1950             }
1951             return this;
1952         },
1953 
1954         _recurse_v4: function(curve, t1, t2, x1, y1, x2, y2, level) {
1955             var tol = 2,
1956                 t = (t1 + t2) * 0.5,
1957                 x = curve.X(t, true),
1958                 y = curve.Y(t, true),
1959                 dx, dy;
1960 
1961             //console.log("Level", level)
1962             if (level === 0) {
1963                 this._insertPoint_v4(curve, [1, NaN, NaN], t);
1964                 return;
1965             }
1966             // console.log("R", t1, t2)
1967             dx = (x - x1) * curve.board.unitX;
1968             dy = (y - y1) * curve.board.unitY;
1969             // console.log("D1", Math.sqrt(dx * dx + dy * dy))
1970             if (Math.sqrt(dx * dx + dy * dy) > tol) {
1971                 this._recurse_v4(curve, t1, t, x1, y1, x, y, level - 1);
1972             } else {
1973                 this._insertPoint_v4(curve, [1, x, y], t);
1974             }
1975             dx = (x - x2) * curve.board.unitX;
1976             dy = (y - y2) * curve.board.unitY;
1977             // console.log("D2", Math.sqrt(dx * dx + dy * dy), x-x2, y-y2)
1978             if (Math.sqrt(dx * dx + dy * dy) > tol) {
1979                 this._recurse_v4(curve, t, t2, x, y, x2, y2, level - 1);
1980             } else {
1981                 this._insertPoint_v4(curve, [1, x, y], t);
1982             }
1983         },
1984 
1985         handleSingularity: function(curve, comp, group, x_table, y_table) {
1986             var idx = group.idx,
1987                 t, t1, t2, y_int,
1988                 i1, i2,
1989                 x, y, lo, hi,
1990                 d_lft, d_rgt,
1991                 d_thresh = 100,
1992                 di1 = 5,
1993                 di2 = 3,
1994                 d1, d2;
1995 
1996             t = group.t;
1997             console.log("HandleSingularity at t =", t);
1998             // console.log(comp.t_values[idx - 1], comp.y_values[idx - 1], comp.t_values[idx + 1], comp.y_values[idx + 1]);
1999             // console.log(group);
2000 
2001             // Look at two points, hopefully left and right from the critical point
2002             t1 = comp.t_values[idx - di1];
2003             t2 = comp.t_values[idx + di1];
2004 
2005             y_int = this.getInterval(curve, t1, t2);
2006             if (Type.isObject(y_int)) {
2007                 lo = y_int.lo;
2008                 hi = y_int.hi;
2009             } else {
2010                 if (y_table[0][idx - 1] < y_table[0][idx + 1]) {
2011                     lo = y_table[0][idx - 1];
2012                     hi = y_table[0][idx + 1];
2013                 } else {
2014                     lo = y_table[0][idx + 1];
2015                     hi = y_table[0][idx - 1];
2016                 }
2017             }
2018 
2019             x = curve.X(t, true);
2020 
2021             d_lft = (y_table[0][idx - di2] - y_table[0][idx - di1]) / (comp.t_values[idx - di2] - comp.t_values[idx - di1]);
2022             d_rgt = (y_table[0][idx + di2] - y_table[0][idx + di1]) / (comp.t_values[idx + di2] - comp.t_values[idx + di1]);
2023 
2024             console.log(":::", d_lft, d_rgt);
2025 
2026             //this._insertPoint_v4(curve, [1, NaN, NaN], 0);
2027 
2028             if (d_lft < -d_thresh) {
2029                 // Left branch very steep downwards -> add the minimum
2030                 this._insertPoint_v4(curve, [1, x, lo], t, true);
2031                 if (d_rgt <= d_thresh) {
2032                     // Right branch not very steep upwards -> interrupt the curve
2033                     // I.e. it looks like -infty / (finite or infty) and not like -infty / -infty
2034                     this._insertPoint_v4(curve, [1, NaN, NaN], t);
2035                 }
2036             } else if (d_lft > d_thresh) {
2037                 // Left branch very steep upwards -> add the maximum
2038                 this._insertPoint_v4(curve, [1, x, hi], t);
2039                 if (d_rgt >= -d_thresh) {
2040                     // Right branch not very steep downwards -> interrupt the curve
2041                     // I.e. it looks like infty / (finite or -infty) and not like infty / infty
2042                     this._insertPoint_v4(curve, [1, NaN, NaN], t);
2043                 }
2044             } else {
2045                 if (lo === -Infinity) {
2046                     this._insertPoint_v4(curve, [1, x, lo], t, true);
2047                     this._insertPoint_v4(curve, [1, NaN, NaN], t);
2048                 }
2049                 if (hi === Infinity) {
2050                     this._insertPoint_v4(curve, [1, NaN, NaN], t);
2051                     this._insertPoint_v4(curve, [1, x, hi], t, true);
2052                 }
2053 
2054                 if (group.t < comp.t_values[idx]) {
2055                     i1 = idx - 1;
2056                     i2 = idx;
2057                 } else {
2058                     i1 = idx;
2059                     i2 = idx + 1;
2060                 }
2061                 t1 = comp.t_values[i1];
2062                 t2 = comp.t_values[i2];
2063                 this._recurse_v4(curve, t1, t2,
2064                         x_table[0][i1],
2065                         y_table[0][i1],
2066                         x_table[0][i2],
2067                         y_table[0][i2],
2068                         10
2069                     );
2070 
2071                 // x = (x_table[0][idx] - x_table[0][idx - 1]) * curve.board.unitX;
2072                 // y = (y_table[0][idx] - y_table[0][idx - 1]) * curve.board.unitY;
2073                 // d1 = Math.sqrt(x * x + y * y);
2074                 // x = (x_table[0][idx + 1] - x_table[0][idx]) * curve.board.unitX;
2075                 // y = (y_table[0][idx + 1] - y_table[0][idx]) * curve.board.unitY;
2076                 // d2 = Math.sqrt(x * x + y * y);
2077 
2078                 // console.log("end", t1, t2, t);
2079                 // if (true || (d1 > 2 || d2 > 2)) {
2080 
2081 // console.log(d1, d2, y_table[0][idx])
2082 //                     // Finite jump
2083 //                     this._insertPoint_v4(curve, [1, NaN, NaN], t);
2084 //                 } else {
2085 //                     if (lo !== -Infinity && hi !== Infinity) {
2086 //                         // Critical point which can be ignored
2087 //                         this._insertPoint_v4(curve, [1, x_table[0][idx], y_table[0][idx]], comp.t_values[idx]);
2088 //                     } else {
2089 //                         if (lo === -Infinity) {
2090 //                             this._insertPoint_v4(curve, [1, x, lo], t, true);
2091 //                             this._insertPoint_v4(curve, [1, NaN, NaN], t);
2092 //                         }
2093 //                         if (hi === Infinity) {
2094 //                             this._insertPoint_v4(curve, [1, NaN, NaN], t);
2095 //                             this._insertPoint_v4(curve, [1, x, hi], t, true);
2096 //                         }
2097 //                     }
2098                 // }
2099             }
2100             if (d_rgt < -d_thresh) {
2101                 // Right branch very steep downwards -> add the maximum
2102                 this._insertPoint_v4(curve, [1, x, hi], t);
2103             } else if (d_rgt > d_thresh) {
2104                 // Right branch very steep upwards -> add the minimum
2105                 this._insertPoint_v4(curve, [1, x, lo], t);
2106             }
2107 
2108         },
2109 
2110         /**
2111          * Number of equidistant points where the function is evaluated
2112          */
2113         steps: 1021, //2053, // 1021,
2114 
2115         /**
2116          * If the absolute maximum of the set of differences is larger than
2117          * criticalThreshold * median of these values, it is regarded as critical point.
2118          * @see JXG.Math.Plot#_criticalInterval
2119          */
2120         criticalThreshold: 1000,
2121 
2122         plot_v4: function(curve, ta, tb, steps) {
2123             var i, j, le, components, idx, comp,
2124                 groups, g, start,
2125                 ret, x_table, y_table,
2126                 t, t1, t2,
2127                 good, bad,
2128                 x_int, y_int,
2129                 degree_x, degree_y,
2130                 h  = (tb - ta) / steps,
2131                 Ypl = function(x) { return curve.Y(x, true); },
2132                 Ymi = function(x) { return -curve.Y(x, true); },
2133                 h2 = h * 0.5;
2134 
2135             components = this.findComponents(curve, ta, tb, steps);
2136             for (idx = 0; idx < components.length; idx++) {
2137                 comp = components[idx];
2138                 ret = this.differenceMethod(comp, curve);
2139                 groups = ret[0];
2140                 x_table = ret[1];
2141                 y_table = ret[2];
2142                 degree_x = ret[3];
2143                 degree_y = ret[4];
2144 
2145                 // if (degree_x >= 0) {
2146                 //     console.log("x polynomial of degree", degree_x);
2147                 // }
2148                 // if (degree_y >= 0) {
2149                 //     console.log("y polynomial of degree", degree_y);
2150                 // }
2151                 if (groups.length === 0 || groups[0].type !== 'borderleft') {
2152                     groups.unshift({
2153                         idx: 0,
2154                         t: comp.t_values[0],
2155                         x: comp.x_values[0],
2156                         y: comp.y_values[0],
2157                         type: 'borderleft'
2158                     });
2159                 }
2160                 if (groups[groups.length - 1].type !== 'borderright') {
2161                     le = comp.t_values.length;
2162                     groups.push({
2163                         idx: le - 1,
2164                         t: comp.t_values[le - 1],
2165                         x: comp.x_values[le - 1],
2166                         y: comp.y_values[le - 1],
2167                         type: 'borderright'
2168                     });
2169                 }
2170 
2171 
2172                 start = 0;
2173                 for (g = 0; g <= groups.length; g++) {
2174                     if (g === groups.length) {
2175                         le = comp.len;
2176                     } else {
2177                         le = groups[g].idx - 1;
2178                     }
2179 
2180                     good = 0;
2181                     bad = 0;
2182                     // Insert all uncritical points until next critical point
2183                     for (i = start; i < le - 2; i++) {
2184                         this._insertPoint_v4(curve, [1, comp.x_values[i], comp.y_values[i]], comp.t_values[i]);
2185                         j = Math.max(0, i - 2);
2186                         // Add more points in critical intervals
2187                         if (true &&
2188                             //degree_y === -1 && // No polynomial
2189                             i >= start + 3 &&
2190                             i < le - 3 &&               // Do not do this if too close to a critical point
2191                             y_table.length > 3 &&
2192                             Math.abs(y_table[2][i]) > 0.2 * Math.abs(y_table[0][i])) {
2193                             t = comp.t_values[i];
2194                             h2 = h * 0.25;
2195                             y_int = this.getInterval(curve, t, t + h);
2196                             if (Type.isObject(y_int)) {
2197                                 if (y_table[2][i] > 0) {
2198                                     this._insertPoint_v4(curve, [1, t + h2, y_int.lo], t + h2);
2199                                 } else {
2200                                     this._insertPoint_v4(curve, [1, t + h - h2, y_int.hi], t + h - h2);
2201                                 }
2202                             } else {
2203                                 t1 = Numerics.fminbr(Ypl, [t, t + h]);
2204                                 t2 = Numerics.fminbr(Ymi, [t, t + h]);
2205                                 if (t1 < t2) {
2206                                     this._insertPoint_v4(curve, [1, curve.X(t1, true), curve.Y(t1, true)], t1);
2207                                     this._insertPoint_v4(curve, [1, curve.X(t2, true), curve.Y(t2, true)], t2);
2208                                 } else {
2209                                     this._insertPoint_v4(curve, [1, curve.X(t2, true), curve.Y(t2, true)], t2);
2210                                     this._insertPoint_v4(curve, [1, curve.X(t1, true), curve.Y(t1, true)], t1);
2211                                 }
2212                             }
2213                             bad++;
2214                         } else {
2215                             good++;
2216                         }
2217                     }
2218                     // console.log("GOOD", good, "BAD", bad);
2219 
2220                     // Handle next critical point
2221                     if (g < groups.length) {
2222                         //console.log("critical point / interval", groups[g]);
2223 
2224                         i = groups[g].idx;
2225                         if (groups[g].type === 'borderleft' || groups[g].type === 'borderright') {
2226                             this.handleBorder(curve, comp, groups[g], x_table, y_table);
2227                         } else {
2228                             this._seconditeration_v4(curve, comp, groups[g], x_table, y_table);
2229                         }
2230 
2231                         start = groups[g].idx + 1 + 1;
2232                     }
2233                 }
2234 
2235                 le = comp.len;
2236                 if (idx < components.length - 1) {
2237                     this._insertPoint_v4(curve, [1, NaN, NaN], comp.right_t);
2238                 }
2239             }
2240 
2241 
2242         },
2243 
2244         /**
2245          * Updates the data points of a parametric curve, plotVersion 4. This version is used if {@link JXG.Curve#plotVersion} is <tt>4</tt>.
2246          * @param {JXG.Curve} curve JSXGraph curve element
2247          * @param {Number} mi Left bound of curve
2248          * @param {Number} ma Right bound of curve
2249          * @returns {JXG.Curve} Reference to the curve object.
2250          */
2251         updateParametricCurve_v4: function (curve, mi, ma) {
2252             var ta, tb, w2, bbox;
2253 
2254             if (curve.xterm === 'x') {
2255                 // For function graphs we can restrict the plot interval
2256                 // to the visible area +plus margin
2257                 bbox = curve.board.getBoundingBox();
2258                 w2 = (bbox[2] - bbox[0]) * 0.3;
2259                 // h2 = (bbox[1] - bbox[3]) * 0.3;
2260                 ta = Math.max(mi, bbox[0] - w2);
2261                 tb = Math.min(ma, bbox[2] + w2);
2262             } else {
2263                 ta = mi;
2264                 tb = ma;
2265             }
2266 
2267             curve.points = [];
2268 
2269             //console.log("--------------------");
2270             this.plot_v4(curve, ta, tb, this.steps);
2271 
2272             curve.numberPoints = curve.points.length;
2273             //console.log(curve.numberPoints);
2274         },
2275 
2276         //----------------------------------------------------------------------
2277         // Plot algorithm alias
2278         //----------------------------------------------------------------------
2279 
2280         /**
2281          * Updates the data points of a parametric curve, alias for {@link JXG.Curve#updateParametricCurve_v2}.
2282          * This is needed for backwards compatibility, if this method has been
2283          * used directly in an application.
2284          * @param {JXG.Curve} curve JSXGraph curve element
2285          * @param {Number} mi Left bound of curve
2286          * @param {Number} ma Right bound of curve
2287          * @returns {JXG.Curve} Reference to the curve object.
2288          *
2289          * @see JXG.Curve#updateParametricCurve_v2
2290          */
2291         updateParametricCurve: function (curve, mi, ma) {
2292             return this.updateParametricCurve_v2(curve, mi, ma);
2293         }
2294     };
2295 
2296 
2297     return Mat.Plot;
2298 });
2299