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