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