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