1 /*
  2     Copyright 2008-2023
  3         Matthias Ehmann,
  4         Michael Gerhaeuser,
  5         Carsten Miller,
  6         Bianca Valentin,
  7         Alfred Wassermann,
  8         Peter Wilfahrt
  9 
 10     This file is part of JSXGraph.
 11 
 12     JSXGraph is free software dual licensed under the GNU LGPL or MIT License.
 13 
 14     You can redistribute it and/or modify it under the terms of the
 15 
 16       * GNU Lesser General Public License as published by
 17         the Free Software Foundation, either version 3 of the License, or
 18         (at your option) any later version
 19       OR
 20       * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT
 21 
 22     JSXGraph is distributed in the hope that it will be useful,
 23     but WITHOUT ANY WARRANTY; without even the implied warranty of
 24     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 25     GNU Lesser General Public License for more details.
 26 
 27     You should have received a copy of the GNU Lesser General Public License and
 28     the MIT License along with JSXGraph. If not, see <https://www.gnu.org/licenses/>
 29     and <https://opensource.org/licenses/MIT/>.
 30  */
 31 
 32 /*global JXG: true, define: true, Float32Array: true */
 33 /*jslint nomen: true, plusplus: true, bitwise: true*/
 34 
 35 /**
 36  * @fileoverview In this file the namespace JXG.Math is defined, which is the base namespace
 37  * for namespaces like JXG.Math.Numerics, JXG.Math.Plot, JXG.Math.Statistics, JXG.Math.Clip etc.
 38  */
 39 import JXG from "../jxg";
 40 import Type from "../utils/type";
 41 
 42 var undef,
 43     /*
 44      * Dynamic programming approach for recursive functions.
 45      * From "Speed up your JavaScript, Part 3" by Nicholas C. Zakas.
 46      * @see JXG.Math.factorial
 47      * @see JXG.Math.binomial
 48      * http://blog.thejit.org/2008/09/05/memoization-in-javascript/
 49      *
 50      * This method is hidden, because it is only used in JXG.Math. If someone wants
 51      * to use it in JSXGraph outside of JXG.Math, it should be moved to jsxgraph.js
 52      */
 53     memoizer = function (f) {
 54         var cache, join;
 55 
 56         if (f.memo) {
 57             return f.memo;
 58         }
 59 
 60         cache = {};
 61         join = Array.prototype.join;
 62 
 63         f.memo = function () {
 64             var key = join.call(arguments);
 65 
 66             // Seems to be a bit faster than "if (a in b)"
 67             return cache[key] !== undef ? cache[key] : (cache[key] = f.apply(this, arguments));
 68         };
 69 
 70         return f.memo;
 71     };
 72 
 73 /**
 74  * Math namespace. Contains mathematics related methods which are
 75  * specific to JSXGraph or which extend the JavaScript Math class.
 76  * @namespace
 77  */
 78 JXG.Math = {
 79     /**
 80      * eps defines the closeness to zero. If the absolute value of a given number is smaller
 81      * than eps, it is considered to be equal to zero.
 82      * @type Number
 83      */
 84     eps: 0.000001,
 85 
 86     /**
 87      * Determine the relative difference between two numbers.
 88      * @param  {Number} a First number
 89      * @param  {Number} b Second number
 90      * @returns {Number}  Relative difference between a and b: |a-b| / max(|a|, |b|)
 91      */
 92     relDif: function (a, b) {
 93         var c = Math.abs(a),
 94             d = Math.abs(b);
 95 
 96         d = Math.max(c, d);
 97 
 98         return d === 0.0 ? 0.0 : Math.abs(a - b) / d;
 99     },
100 
101     /**
102      * The JavaScript implementation of the % operator returns the symmetric modulo.
103      * mod and "%" are both identical if a >= 0 and m >= 0 but the results differ if a or m < 0.
104      * @param {Number} a
105      * @param {Number} m
106      * @returns {Number} Mathematical modulo <tt>a mod m</tt>
107      */
108     mod: function (a, m) {
109         return a - Math.floor(a / m) * m;
110     },
111 
112     /**
113      * Initializes a vector of size <tt>n</tt> wih coefficients set to the init value (default 0)
114      * @param {Number} n Length of the vector
115      * @param {Number} [init=0] Initial value for each coefficient
116      * @returns {Array} An array of length <tt>n</tt>
117      */
118     vector: function (n, init) {
119         var r, i;
120 
121         init = init || 0;
122         r = [];
123 
124         for (i = 0; i < n; i++) {
125             r[i] = init;
126         }
127 
128         return r;
129     },
130 
131     /**
132      * Initializes a matrix as an array of rows with the given value.
133      * @param {Number} n Number of rows
134      * @param {Number} [m=n] Number of columns
135      * @param {Number} [init=0] Initial value for each coefficient
136      * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a
137      * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows.
138      */
139     matrix: function (n, m, init) {
140         var r, i, j;
141 
142         init = init || 0;
143         m = m || n;
144         r = [];
145 
146         for (i = 0; i < n; i++) {
147             r[i] = [];
148 
149             for (j = 0; j < m; j++) {
150                 r[i][j] = init;
151             }
152         }
153 
154         return r;
155     },
156 
157     /**
158      * Generates an identity matrix. If n is a number and m is undefined or not a number, a square matrix is generated,
159      * if n and m are both numbers, an nxm matrix is generated.
160      * @param {Number} n Number of rows
161      * @param {Number} [m=n] Number of columns
162      * @returns {Array} A square matrix of length <tt>n</tt> with all coefficients equal to 0 except a_(i,i), i out of (1, ..., n), if <tt>m</tt> is undefined or not a number
163      * or a <tt>n</tt> times <tt>m</tt>-matrix with a_(i,j) = 0 and a_(i,i) = 1 if m is a number.
164      */
165     identity: function (n, m) {
166         var r, i;
167 
168         if (m === undef && typeof m !== "number") {
169             m = n;
170         }
171 
172         r = this.matrix(n, m);
173 
174         for (i = 0; i < Math.min(n, m); i++) {
175             r[i][i] = 1;
176         }
177 
178         return r;
179     },
180 
181     /**
182      * Generates a 4x4 matrix for 3D to 2D projections.
183      * @param {Number} l Left
184      * @param {Number} r Right
185      * @param {Number} t Top
186      * @param {Number} b Bottom
187      * @param {Number} n Near
188      * @param {Number} f Far
189      * @returns {Array} 4x4 Matrix
190      */
191     frustum: function (l, r, b, t, n, f) {
192         var ret = this.matrix(4, 4);
193 
194         ret[0][0] = (n * 2) / (r - l);
195         ret[0][1] = 0;
196         ret[0][2] = (r + l) / (r - l);
197         ret[0][3] = 0;
198 
199         ret[1][0] = 0;
200         ret[1][1] = (n * 2) / (t - b);
201         ret[1][2] = (t + b) / (t - b);
202         ret[1][3] = 0;
203 
204         ret[2][0] = 0;
205         ret[2][1] = 0;
206         ret[2][2] = -(f + n) / (f - n);
207         ret[2][3] = -(f * n * 2) / (f - n);
208 
209         ret[3][0] = 0;
210         ret[3][1] = 0;
211         ret[3][2] = -1;
212         ret[3][3] = 0;
213 
214         return ret;
215     },
216 
217     /**
218      * Generates a 4x4 matrix for 3D to 2D projections.
219      * @param {Number} fov Field of view in vertical direction, given in rad.
220      * @param {Number} ratio Aspect ratio of the projection plane.
221      * @param {Number} n Near
222      * @param {Number} f Far
223      * @returns {Array} 4x4 Projection Matrix
224      */
225     projection: function (fov, ratio, n, f) {
226         var t = n * Math.tan(fov / 2),
227             r = t * ratio;
228 
229         return this.frustum(-r, r, -t, t, n, f);
230     },
231 
232     /**
233      * Multiplies a vector vec to a matrix mat: mat * vec. The matrix is interpreted by this function as an array of rows. Please note: This
234      * function does not check if the dimensions match.
235      * @param {Array} mat Two dimensional array of numbers. The inner arrays describe the columns, the outer ones the matrix' rows.
236      * @param {Array} vec Array of numbers
237      * @returns {Array} Array of numbers containing the result
238      * @example
239      * var A = [[2, 1],
240      *          [1, 3]],
241      *     b = [4, 5],
242      *     c;
243      * c = JXG.Math.matVecMult(A, b)
244      * // c === [13, 19];
245      */
246     matVecMult: function (mat, vec) {
247         var i,
248             s,
249             k,
250             m = mat.length,
251             n = vec.length,
252             res = [];
253 
254         if (n === 3) {
255             for (i = 0; i < m; i++) {
256                 res[i] = mat[i][0] * vec[0] + mat[i][1] * vec[1] + mat[i][2] * vec[2];
257             }
258         } else {
259             for (i = 0; i < m; i++) {
260                 s = 0;
261                 for (k = 0; k < n; k++) {
262                     s += mat[i][k] * vec[k];
263                 }
264                 res[i] = s;
265             }
266         }
267         return res;
268     },
269 
270     /**
271      * Computes the product of the two matrices mat1*mat2.
272      * @param {Array} mat1 Two dimensional array of numbers
273      * @param {Array} mat2 Two dimensional array of numbers
274      * @returns {Array} Two dimensional Array of numbers containing result
275      */
276     matMatMult: function (mat1, mat2) {
277         var i,
278             j,
279             s,
280             k,
281             m = mat1.length,
282             n = m > 0 ? mat2[0].length : 0,
283             m2 = mat2.length,
284             res = this.matrix(m, n);
285 
286         for (i = 0; i < m; i++) {
287             for (j = 0; j < n; j++) {
288                 s = 0;
289                 for (k = 0; k < m2; k++) {
290                     s += mat1[i][k] * mat2[k][j];
291                 }
292                 res[i][j] = s;
293             }
294         }
295         return res;
296     },
297 
298     /**
299      * Transposes a matrix given as a two dimensional array.
300      * @param {Array} M The matrix to be transposed
301      * @returns {Array} The transpose of M
302      */
303     transpose: function (M) {
304         var MT, i, j, m, n;
305 
306         // number of rows of M
307         m = M.length;
308         // number of columns of M
309         n = M.length > 0 ? M[0].length : 0;
310         MT = this.matrix(n, m);
311 
312         for (i = 0; i < n; i++) {
313             for (j = 0; j < m; j++) {
314                 MT[i][j] = M[j][i];
315             }
316         }
317 
318         return MT;
319     },
320 
321     /**
322      * Compute the inverse of an nxn matrix with Gauss elimination.
323      * @param {Array} Ain
324      * @returns {Array} Inverse matrix of Ain
325      */
326     inverse: function (Ain) {
327         var i,
328             j,
329             k,
330             s,
331             ma,
332             r,
333             swp,
334             n = Ain.length,
335             A = [],
336             p = [],
337             hv = [];
338 
339         for (i = 0; i < n; i++) {
340             A[i] = [];
341             for (j = 0; j < n; j++) {
342                 A[i][j] = Ain[i][j];
343             }
344             p[i] = i;
345         }
346 
347         for (j = 0; j < n; j++) {
348             // pivot search:
349             ma = Math.abs(A[j][j]);
350             r = j;
351 
352             for (i = j + 1; i < n; i++) {
353                 if (Math.abs(A[i][j]) > ma) {
354                     ma = Math.abs(A[i][j]);
355                     r = i;
356                 }
357             }
358 
359             // Singular matrix
360             if (ma <= this.eps) {
361                 return [];
362             }
363 
364             // swap rows:
365             if (r > j) {
366                 for (k = 0; k < n; k++) {
367                     swp = A[j][k];
368                     A[j][k] = A[r][k];
369                     A[r][k] = swp;
370                 }
371 
372                 swp = p[j];
373                 p[j] = p[r];
374                 p[r] = swp;
375             }
376 
377             // transformation:
378             s = 1.0 / A[j][j];
379             for (i = 0; i < n; i++) {
380                 A[i][j] *= s;
381             }
382             A[j][j] = s;
383 
384             for (k = 0; k < n; k++) {
385                 if (k !== j) {
386                     for (i = 0; i < n; i++) {
387                         if (i !== j) {
388                             A[i][k] -= A[i][j] * A[j][k];
389                         }
390                     }
391                     A[j][k] = -s * A[j][k];
392                 }
393             }
394         }
395 
396         // swap columns:
397         for (i = 0; i < n; i++) {
398             for (k = 0; k < n; k++) {
399                 hv[p[k]] = A[i][k];
400             }
401             for (k = 0; k < n; k++) {
402                 A[i][k] = hv[k];
403             }
404         }
405 
406         return A;
407     },
408 
409     /**
410      * Inner product of two vectors a and b. n is the length of the vectors.
411      * @param {Array} a Vector
412      * @param {Array} b Vector
413      * @param {Number} [n] Length of the Vectors. If not given the length of the first vector is taken.
414      * @returns {Number} The inner product of a and b.
415      */
416     innerProduct: function (a, b, n) {
417         var i,
418             s = 0;
419 
420         if (n === undef || !Type.isNumber(n)) {
421             n = a.length;
422         }
423 
424         for (i = 0; i < n; i++) {
425             s += a[i] * b[i];
426         }
427 
428         return s;
429     },
430 
431     /**
432      * Calculates the cross product of two vectors both of length three.
433      * In case of homogeneous coordinates this is either
434      * <ul>
435      * <li>the intersection of two lines</li>
436      * <li>the line through two points</li>
437      * </ul>
438      * @param {Array} c1 Homogeneous coordinates of line or point 1
439      * @param {Array} c2 Homogeneous coordinates of line or point 2
440      * @returns {Array} vector of length 3: homogeneous coordinates of the resulting point / line.
441      */
442     crossProduct: function (c1, c2) {
443         return [
444             c1[1] * c2[2] - c1[2] * c2[1],
445             c1[2] * c2[0] - c1[0] * c2[2],
446             c1[0] * c2[1] - c1[1] * c2[0]
447         ];
448     },
449 
450     /**
451      * Euclidean norm of a vector.
452      *
453      * @param {Array} a Array containing a vector.
454      * @param {Number} n (Optional) length of the array.
455      * @returns {Number} Euclidean norm of the vector.
456      */
457     norm: function (a, n) {
458         var i,
459             sum = 0.0;
460 
461         if (n === undef || !Type.isNumber(n)) {
462             n = a.length;
463         }
464 
465         for (i = 0; i < n; i++) {
466             sum += a[i] * a[i];
467         }
468 
469         return Math.sqrt(sum);
470     },
471 
472     /**
473      * Compute a * x + y for a scalar a and vectors x and y.
474      *
475      * @param {Number} a
476      * @param {Array} x
477      * @param {Array} y
478      * @returns
479      */
480     axpy: function (a, x, y) {
481         var i,
482             le = x.length,
483             p = [];
484         for (i = 0; i < le; i++) {
485             p[i] = a * x[i] + y[i];
486         }
487         return p;
488     },
489 
490     /**
491      * Compute the factorial of a positive integer. If a non-integer value
492      * is given, the fraction will be ignored.
493      * @function
494      * @param {Number} n
495      * @returns {Number} n! = n*(n-1)*...*2*1
496      */
497     factorial: memoizer(function (n) {
498         if (n < 0) {
499             return NaN;
500         }
501 
502         n = Math.floor(n);
503 
504         if (n === 0 || n === 1) {
505             return 1;
506         }
507 
508         return n * this.factorial(n - 1);
509     }),
510 
511     /**
512      * Computes the binomial coefficient n over k.
513      * @function
514      * @param {Number} n Fraction will be ignored
515      * @param {Number} k Fraction will be ignored
516      * @returns {Number} The binomial coefficient n over k
517      */
518     binomial: memoizer(function (n, k) {
519         var b, i;
520 
521         if (k > n || k < 0) {
522             return NaN;
523         }
524 
525         k = Math.round(k);
526         n = Math.round(n);
527 
528         if (k === 0 || k === n) {
529             return 1;
530         }
531 
532         b = 1;
533 
534         for (i = 0; i < k; i++) {
535             b *= n - i;
536             b /= i + 1;
537         }
538 
539         return b;
540     }),
541 
542     /**
543      * Calculates the cosine hyperbolicus of x.
544      * @function
545      * @param {Number} x The number the cosine hyperbolicus will be calculated of.
546      * @returns {Number} Cosine hyperbolicus of the given value.
547      */
548     cosh:
549         Math.cosh ||
550         function (x) {
551             return (Math.exp(x) + Math.exp(-x)) * 0.5;
552         },
553 
554     /**
555      * Sine hyperbolicus of x.
556      * @function
557      * @param {Number} x The number the sine hyperbolicus will be calculated of.
558      * @returns {Number} Sine hyperbolicus of the given value.
559      */
560     sinh:
561         Math.sinh ||
562         function (x) {
563             return (Math.exp(x) - Math.exp(-x)) * 0.5;
564         },
565 
566     /**
567      * Hyperbolic arc-cosine of a number.
568      *
569      * @param {Number} x
570      * @returns {Number}
571      */
572     acosh:
573         Math.acosh ||
574         function (x) {
575             return Math.log(x + Math.sqrt(x * x - 1));
576         },
577 
578     /**
579      * Hyperbolic arcsine of a number
580      * @param {Number} x
581      * @returns {Number}
582      */
583     asinh:
584         Math.asinh ||
585         function (x) {
586             if (x === -Infinity) {
587                 return x;
588             }
589             return Math.log(x + Math.sqrt(x * x + 1));
590         },
591 
592     /**
593      * Computes the cotangent of x.
594      * @function
595      * @param {Number} x The number the cotangent will be calculated of.
596      * @returns {Number} Cotangent of the given value.
597      */
598     cot: function (x) {
599         return 1 / Math.tan(x);
600     },
601 
602     /**
603      * Computes the inverse cotangent of x.
604      * @param {Number} x The number the inverse cotangent will be calculated of.
605      * @returns {Number} Inverse cotangent of the given value.
606      */
607     acot: function (x) {
608         return (x >= 0 ? 0.5 : -0.5) * Math.PI - Math.atan(x);
609     },
610 
611     /**
612      * Compute n-th real root of a real number. n must be strictly positive integer.
613      * If n is odd, the real n-th root exists and is negative.
614      * For n even, for negative valuees of x NaN is returned
615      * @param  {Number} x radicand. Must be non-negative, if n even.
616      * @param  {Number} n index of the root. must be strictly positive integer.
617      * @returns {Number} returns real root or NaN
618      *
619      * @example
620      * nthroot(16, 4): 2
621      * nthroot(-27, 3): -3
622      * nthroot(-4, 2): NaN
623      */
624     nthroot: function (x, n) {
625         var inv = 1 / n;
626 
627         if (n <= 0 || Math.floor(n) !== n) {
628             return NaN;
629         }
630 
631         if (x === 0.0) {
632             return 0.0;
633         }
634 
635         if (x > 0) {
636             return Math.exp(inv * Math.log(x));
637         }
638 
639         // From here on, x is negative
640         if (n % 2 === 1) {
641             return -Math.exp(inv * Math.log(-x));
642         }
643 
644         // x negative, even root
645         return NaN;
646     },
647 
648     /**
649      * Computes cube root of real number
650      * Polyfill for Math.cbrt().
651      *
652      * @function
653      * @param  {Number} x Radicand
654      * @returns {Number} Cube root of x.
655      */
656     cbrt:
657         Math.cbrt ||
658         function (x) {
659             return this.nthroot(x, 3);
660         },
661 
662     /**
663      * Compute base to the power of exponent.
664      * @param {Number} base
665      * @param {Number} exponent
666      * @returns {Number} base to the power of exponent.
667      */
668     pow: function (base, exponent) {
669         if (base === 0) {
670             if (exponent === 0) {
671                 return 1;
672             }
673             return 0;
674         }
675 
676         // exponent is an integer
677         if (Math.floor(exponent) === exponent) {
678             return Math.pow(base, exponent);
679         }
680 
681         // exponent is not an integer
682         if (base > 0) {
683             return Math.exp(exponent * Math.log(base));
684         }
685 
686         return NaN;
687     },
688 
689     /**
690      * Compute base to the power of the rational exponent m / n.
691      * This function first reduces the fraction m/n and then computes
692      * JXG.Math.pow(base, m/n).
693      *
694      * This function is necessary to have the same results for e.g.
695      * (-8)^(1/3) = (-8)^(2/6) = -2
696      * @param {Number} base
697      * @param {Number} m numerator of exponent
698      * @param {Number} n denominator of exponent
699      * @returns {Number} base to the power of exponent.
700      */
701     ratpow: function (base, m, n) {
702         var g;
703         if (m === 0) {
704             return 1;
705         }
706         if (n === 0) {
707             return NaN;
708         }
709 
710         g = this.gcd(m, n);
711         return this.nthroot(this.pow(base, m / g), n / g);
712     },
713 
714     /**
715      * Logarithm to base 10.
716      * @param {Number} x
717      * @returns {Number} log10(x) Logarithm of x to base 10.
718      */
719     log10: function (x) {
720         return Math.log(x) / Math.log(10.0);
721     },
722 
723     /**
724      * Logarithm to base 2.
725      * @param {Number} x
726      * @returns {Number} log2(x) Logarithm of x to base 2.
727      */
728     log2: function (x) {
729         return Math.log(x) / Math.log(2.0);
730     },
731 
732     /**
733      * Logarithm to arbitrary base b. If b is not given, natural log is taken, i.e. b = e.
734      * @param {Number} x
735      * @param {Number} b base
736      * @returns {Number} log(x, b) Logarithm of x to base b, that is log(x)/log(b).
737      */
738     log: function (x, b) {
739         if (b !== undefined && Type.isNumber(b)) {
740             return Math.log(x) / Math.log(b);
741         }
742 
743         return Math.log(x);
744     },
745 
746     /**
747      * The sign() function returns the sign of a number, indicating whether the number is positive, negative or zero.
748      *
749      * @function
750      * @param  {Number} x A Number
751      * @returns {Number}  This function has 5 kinds of return values,
752      *    1, -1, 0, -0, NaN, which represent "positive number", "negative number", "positive zero", "negative zero"
753      *    and NaN respectively.
754      */
755     sign:
756         Math.sign ||
757         function (x) {
758             x = +x; // convert to a number
759             if (x === 0 || isNaN(x)) {
760                 return x;
761             }
762             return x > 0 ? 1 : -1;
763         },
764 
765     /**
766      * A square & multiply algorithm to compute base to the power of exponent.
767      * Implementated by Wolfgang Riedl.
768      *
769      * @param {Number} base
770      * @param {Number} exponent
771      * @returns {Number} Base to the power of exponent
772      */
773     squampow: function (base, exponent) {
774         var result;
775 
776         if (Math.floor(exponent) === exponent) {
777             // exponent is integer (could be zero)
778             result = 1;
779 
780             if (exponent < 0) {
781                 // invert: base
782                 base = 1.0 / base;
783                 exponent *= -1;
784             }
785 
786             while (exponent !== 0) {
787                 if (exponent & 1) {
788                     result *= base;
789                 }
790 
791                 exponent >>= 1;
792                 base *= base;
793             }
794             return result;
795         }
796 
797         return this.pow(base, exponent);
798     },
799 
800     /**
801      * Greatest common divisor (gcd) of two numbers.
802      * @see <a href="https://rosettacode.org/wiki/Greatest_common_divisor#JavaScript">rosettacode.org</a>
803      *
804      * @param  {Number} a First number
805      * @param  {Number} b Second number
806      * @returns {Number}   gcd(a, b) if a and b are numbers, NaN else.
807      */
808     gcd: function (a, b) {
809         var tmp,
810             endless = true;
811 
812         a = Math.abs(a);
813         b = Math.abs(b);
814 
815         if (!(Type.isNumber(a) && Type.isNumber(b))) {
816             return NaN;
817         }
818         if (b > a) {
819             tmp = a;
820             a = b;
821             b = tmp;
822         }
823 
824         while (endless) {
825             a %= b;
826             if (a === 0) {
827                 return b;
828             }
829             b %= a;
830             if (b === 0) {
831                 return a;
832             }
833         }
834     },
835 
836     /**
837      * Least common multiple (lcm) of two numbers.
838      *
839      * @param  {Number} a First number
840      * @param  {Number} b Second number
841      * @returns {Number}   lcm(a, b) if a and b are numbers, NaN else.
842      */
843     lcm: function (a, b) {
844         var ret;
845 
846         if (!(Type.isNumber(a) && Type.isNumber(b))) {
847             return NaN;
848         }
849 
850         ret = a * b;
851         if (ret !== 0) {
852             return ret / this.gcd(a, b);
853         }
854 
855         return 0;
856     },
857 
858     /**
859      * Special use of Math.round function to round not only to integers but also to chosen decimal values.
860      *
861      * @param {Number} value Value to be rounded.
862      * @param {Number} step Distance between the values to be rounded to. (default: 1.0)
863      * @param {Number} [min] If set, it will be returned the maximum of value and min.
864      * @param {Number} [max] If set, it will be returned the minimum of value and max.
865      * @returns {Number} Fitted value.
866      */
867     roundToStep: function (value, step, min, max) {
868         var n = value,
869             tmp, minOr0;
870 
871         // for performance
872         if (!Type.exists(step) && !Type.exists(min) && !Type.exists(max)) {
873             return n;
874         }
875 
876         if (JXG.exists(max)) {
877             n = Math.min(n, max);
878         }
879         if (JXG.exists(min)) {
880             n = Math.max(n, min);
881         }
882 
883         minOr0 = min || 0;
884 
885         if ( JXG.exists(step)) {
886             tmp = (n - minOr0) / step;
887             if (Number.isInteger(tmp)) {
888                 return n;
889             }
890 
891             tmp = Math.round(tmp);
892             n = minOr0 + tmp * step;
893         }
894 
895         if (JXG.exists(max)) {
896             n = Math.min(n, max);
897         }
898         if (JXG.exists(min)) {
899             n = Math.max(n, min);
900         }
901 
902         return n;
903     },
904 
905     /**
906      *  Error function, see {@link https://en.wikipedia.org/wiki/Error_function}.
907      *
908      * @see JXG.Math.PropFunc.erf
909      * @param  {Number} x
910      * @returns {Number}
911      */
912     erf: function (x) {
913         return this.ProbFuncs.erf(x);
914     },
915 
916     /**
917      * Complementary error function, i.e. 1 - erf(x).
918      *
919      * @see JXG.Math.erf
920      * @see JXG.Math.PropFunc.erfc
921      * @param  {Number} x
922      * @returns {Number}
923      */
924     erfc: function (x) {
925         return this.ProbFuncs.erfc(x);
926     },
927 
928     /**
929      * Inverse of error function
930      *
931      * @see JXG.Math.erf
932      * @see JXG.Math.PropFunc.erfi
933      * @param  {Number} x
934      * @returns {Number}
935      */
936     erfi: function (x) {
937         return this.ProbFuncs.erfi(x);
938     },
939 
940     /**
941      * Normal distribution function
942      *
943      * @see JXG.Math.PropFunc.ndtr
944      * @param  {Number} x
945      * @returns {Number}
946      */
947     ndtr: function (x) {
948         return this.ProbFuncs.ndtr(x);
949     },
950 
951     /**
952      * Inverse of normal distribution function
953      *
954      * @see JXG.Math.ndtr
955      * @see JXG.Math.PropFunc.ndtri
956      * @param  {Number} x
957      * @returns {Number}
958      */
959     ndtri: function (x) {
960         return this.ProbFuncs.ndtri(x);
961     },
962 
963     /**
964      * Returns sqrt(a * a + b * b) for a variable number of arguments.
965      * This is a naive implementation which might be faster than Math.hypot.
966      * The latter is numerically more stable.
967      *
968      * @param {Number} a Variable number of arguments.
969      * @returns Number
970      */
971     hypot: function() {
972         var i, le, a, sum;
973 
974         le = arguments.length;
975         for (i = 0, sum = 0.0; i < le; i++) {
976             a = arguments[i];
977             sum += a * a;
978         }
979         return Math.sqrt(sum);
980     },
981 
982     /**
983      * Heaviside unit step function. Returns 0 for x <, 1 for x > 0, and 0.5 for x == 0.
984      *
985      * @param {Number} x
986      * @returns Number
987      */
988     hstep: function(x) {
989         return (x > 0.0) ? 1 :
990             ((x < 0.0) ? 0.0 : 0.5);
991     },
992 
993     /* ********************  Comparisons and logical operators ************** */
994 
995     /**
996      * Logical test: a < b?
997      *
998      * @param {Number} a
999      * @param {Number} b
1000      * @returns {Boolean}
1001      */
1002     lt: function (a, b) {
1003         return a < b;
1004     },
1005 
1006     /**
1007      * Logical test: a <= b?
1008      *
1009      * @param {Number} a
1010      * @param {Number} b
1011      * @returns {Boolean}
1012      */
1013     leq: function (a, b) {
1014         return a <= b;
1015     },
1016 
1017     /**
1018      * Logical test: a > b?
1019      *
1020      * @param {Number} a
1021      * @param {Number} b
1022      * @returns {Boolean}
1023      */
1024     gt: function (a, b) {
1025         return a > b;
1026     },
1027 
1028     /**
1029      * Logical test: a >= b?
1030      *
1031      * @param {Number} a
1032      * @param {Number} b
1033      * @returns {Boolean}
1034      */
1035     geq: function (a, b) {
1036         return a >= b;
1037     },
1038 
1039     /**
1040      * Logical test: a === b?
1041      *
1042      * @param {Number} a
1043      * @param {Number} b
1044      * @returns {Boolean}
1045      */
1046     eq: function (a, b) {
1047         return a === b;
1048     },
1049 
1050     /**
1051      * Logical test: a !== b?
1052      *
1053      * @param {Number} a
1054      * @param {Number} b
1055      * @returns {Boolean}
1056      */
1057     neq: function (a, b) {
1058         return a !== b;
1059     },
1060 
1061     /**
1062      * Logical operator: a && b?
1063      *
1064      * @param {Boolean} a
1065      * @param {Boolean} b
1066      * @returns {Boolean}
1067      */
1068     and: function (a, b) {
1069         return a && b;
1070     },
1071 
1072     /**
1073      * Logical operator: !a?
1074      *
1075      * @param {Boolean} a
1076      * @returns {Boolean}
1077      */
1078     not: function (a) {
1079         return !a;
1080     },
1081 
1082     /**
1083      * Logical operator: a || b?
1084      *
1085      * @param {Boolean} a
1086      * @param {Boolean} b
1087      * @returns {Boolean}
1088      */
1089     or: function (a, b) {
1090         return a || b;
1091     },
1092 
1093     /**
1094      * Logical operator: either a or b?
1095      *
1096      * @param {Boolean} a
1097      * @param {Boolean} b
1098      * @returns {Boolean}
1099      */
1100     xor: function (a, b) {
1101         return (a || b) && !(a && b);
1102     },
1103 
1104     /**
1105      *
1106      * Convert a floating point number to sign + integer + fraction.
1107      * fraction is given as nominator and denominator.
1108      * <p>
1109      * Algorithm: approximate the floating point number
1110      * by a continued fraction and simultaneously keep track
1111      * of its convergents.
1112      * Inspired by {@link https://kevinboone.me/rationalize.html}.
1113      *
1114      * @param {Number} x Number which is to be converted
1115      * @param {Number} [order=0.001] Small number determining the approximation precision.
1116      * @returns {Array} [sign, leading, nominator, denominator] where sign is 1 or -1.
1117      * @see JXG#toFraction
1118      *
1119      * @example
1120      * JXG.Math.decToFraction(0.33333333);
1121      * // Result: [ 1, 0, 1, 3 ]
1122      *
1123      * JXG.Math.decToFraction(0);
1124      * // Result: [ 1, 0, 0, 1 ]
1125      *
1126      * JXG.Math.decToFraction(-10.66666666666667);
1127      * // Result: [-1, 10, 2, 3 ]
1128     */
1129     decToFraction: function(x, order) {
1130         var lead, sign, a,
1131             n, n1, n2,
1132             d, d1, d2,
1133             it = 0,
1134             maxit = 20;
1135 
1136         order = Type.def(order, 0.001);
1137 
1138         // Round the number.
1139         // Otherwise, 0.999999999 would result in [0, 1, 1].
1140         x = Math.round(x * 1.e12) * 1.e-12;
1141 
1142         // Negative numbers:
1143         // The minus sign is handled in sign.
1144         sign = (x < 0) ? -1 : 1;
1145         x = Math.abs(x);
1146 
1147         // From now on we consider x to be nonnegative.
1148         lead = Math.floor(x);
1149         x -= Math.floor(x);
1150         a = 0.0;
1151         n2 = 1.0;
1152         n = n1 = a;
1153         d2 = 0.0;
1154         d = d1 = 1.0;
1155 
1156         while (x - Math.floor(x) > order && it < maxit) {
1157             x = 1 / (x - a);
1158             a = Math.floor(x);
1159             n = n2 + a * n1;
1160             d = d2 + a * d1;
1161             n2 = n1;
1162             d2 = d1;
1163             n1 = n;
1164             d1 = d;
1165             it++;
1166         }
1167         return [sign, lead, n, d];
1168     },
1169 
1170     /* *************************** Normalize *************************** */
1171 
1172     /**
1173      * Normalize the standard form [c, b0, b1, a, k, r, q0, q1].
1174      * @private
1175      * @param {Array} stdform The standard form to be normalized.
1176      * @returns {Array} The normalized standard form.
1177      */
1178     normalize: function (stdform) {
1179         var n,
1180             signr,
1181             a2 = 2 * stdform[3],
1182             r = stdform[4] / a2;
1183 
1184         stdform[5] = r;
1185         stdform[6] = -stdform[1] / a2;
1186         stdform[7] = -stdform[2] / a2;
1187 
1188         if (!isFinite(r)) {
1189             n = this.hypot(stdform[1], stdform[2]);
1190 
1191             stdform[0] /= n;
1192             stdform[1] /= n;
1193             stdform[2] /= n;
1194             stdform[3] = 0;
1195             stdform[4] = 1;
1196         } else if (Math.abs(r) >= 1) {
1197             stdform[0] = (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) / (2 * r);
1198             stdform[1] = -stdform[6] / r;
1199             stdform[2] = -stdform[7] / r;
1200             stdform[3] = 1 / (2 * r);
1201             stdform[4] = 1;
1202         } else {
1203             signr = r <= 0 ? -1 : 1;
1204             stdform[0] =
1205                 signr * (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) * 0.5;
1206             stdform[1] = -signr * stdform[6];
1207             stdform[2] = -signr * stdform[7];
1208             stdform[3] = signr / 2;
1209             stdform[4] = signr * r;
1210         }
1211 
1212         return stdform;
1213     },
1214 
1215     /**
1216      * Converts a two dimensional array to a one dimensional Float32Array that can be processed by WebGL.
1217      * @param {Array} m A matrix in a two dimensional array.
1218      * @returns {Float32Array} A one dimensional array containing the matrix in column wise notation. Provides a fall
1219      * back to the default JavaScript Array if Float32Array is not available.
1220      */
1221     toGL: function (m) {
1222         var v, i, j;
1223 
1224         if (typeof Float32Array === "function") {
1225             v = new Float32Array(16);
1226         } else {
1227             v = new Array(16);
1228         }
1229 
1230         if (m.length !== 4 && m[0].length !== 4) {
1231             return v;
1232         }
1233 
1234         for (i = 0; i < 4; i++) {
1235             for (j = 0; j < 4; j++) {
1236                 v[i + 4 * j] = m[i][j];
1237             }
1238         }
1239 
1240         return v;
1241     },
1242 
1243     /**
1244      * Theorem of Vieta: Given a set of simple zeroes x_0, ..., x_n
1245      * of a polynomial f, compute the coefficients s_k, (k=0,...,n-1)
1246      * of the polynomial of the form. See {@link https://de.wikipedia.org/wiki/Elementarsymmetrisches_Polynom}.
1247      * <p>
1248      *  f(x) = (x-x_0)*...*(x-x_n) =
1249      *  x^n + sum_{k=1}^{n} (-1)^(k) s_{k-1} x^(n-k)
1250      * </p>
1251      * @param {Array} x Simple zeroes of the polynomial.
1252      * @returns {Array} Coefficients of the polynomial.
1253      *
1254      */
1255     Vieta: function (x) {
1256         var n = x.length,
1257             s = [],
1258             m,
1259             k,
1260             y;
1261 
1262         s = x.slice();
1263         for (m = 1; m < n; ++m) {
1264             y = s[m];
1265             s[m] *= s[m - 1];
1266             for (k = m - 1; k >= 1; --k) {
1267                 s[k] += s[k - 1] * y;
1268             }
1269             s[0] += y;
1270         }
1271         return s;
1272     }
1273 };
1274 
1275 export default JXG.Math;
1276