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