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