1 /* 2 Copyright 2008-2015 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 <http://www.gnu.org/licenses/> 29 and <http://opensource.org/licenses/MIT/>. 30 */ 31 32 33 /*global JXG: true, define: true, Float32Array: true */ 34 /*jslint nomen: true, plusplus: true, bitwise: true*/ 35 36 /* depends: 37 jxg 38 */ 39 40 /** 41 * @fileoverview In this file the namespace JXG.Math is defined, which is the base namespace 42 * for namespaces like Math.Numerics, Math.Algebra, Math.Statistics etc. 43 */ 44 45 define(['jxg', 'utils/type'], function (JXG, Type) { 46 47 "use strict"; 48 49 var undef, 50 51 /* 52 * Dynamic programming approach for recursive functions. 53 * From "Speed up your JavaScript, Part 3" by Nicholas C. Zakas. 54 * @see JXG.Math.factorial 55 * @see JXG.Math.binomial 56 * http://blog.thejit.org/2008/09/05/memoization-in-javascript/ 57 * 58 * This method is hidden, because it is only used in JXG.Math. If someone wants 59 * to use it in JSXGraph outside of JXG.Math, it should be moved to jsxgraph.js 60 */ 61 memoizer = function (f) { 62 var cache, join; 63 64 if (f.memo) { 65 return f.memo; 66 } 67 68 cache = {}; 69 join = Array.prototype.join; 70 71 f.memo = function () { 72 var key = join.call(arguments); 73 74 // Seems to be a bit faster than "if (a in b)" 75 return (cache[key] !== undef) ? 76 cache[key] : 77 cache[key] = f.apply(this, arguments); 78 }; 79 80 return f.memo; 81 }; 82 83 /** 84 * Math namespace. 85 * @namespace 86 */ 87 JXG.Math = { 88 /** 89 * eps defines the closeness to zero. If the absolute value of a given number is smaller 90 * than eps, it is considered to be equal to zero. 91 * @type number 92 */ 93 eps: 0.000001, 94 95 /** 96 * The JavaScript implementation of the % operator returns the symmetric modulo. 97 * They are both identical if a >= 0 and m >= 0 but the results differ if a or m < 0. 98 * @param {Number} a 99 * @param {Number} m 100 * @returns {Number} Mathematical modulo <tt>a mod m</tt> 101 */ 102 mod: function (a, m) { 103 return a - Math.floor(a / m) * m; 104 }, 105 106 /** 107 * Initializes a vector as an array with the coefficients set to the given value resp. zero. 108 * @param {Number} n Length of the vector 109 * @param {Number} [init=0] Initial value for each coefficient 110 * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a 111 * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows. 112 */ 113 vector: function (n, init) { 114 var r, i; 115 116 init = init || 0; 117 r = []; 118 119 for (i = 0; i < n; i++) { 120 r[i] = init; 121 } 122 123 return r; 124 }, 125 126 /** 127 * Initializes a matrix as an array of rows with the given value. 128 * @param {Number} n Number of rows 129 * @param {Number} [m=n] Number of columns 130 * @param {Number} [init=0] Initial value for each coefficient 131 * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a 132 * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows. 133 */ 134 matrix: function (n, m, init) { 135 var r, i, j; 136 137 init = init || 0; 138 m = m || n; 139 r = []; 140 141 for (i = 0; i < n; i++) { 142 r[i] = []; 143 144 for (j = 0; j < m; j++) { 145 r[i][j] = init; 146 } 147 } 148 149 return r; 150 }, 151 152 /** 153 * Generates an identity matrix. If n is a number and m is undefined or not a number, a square matrix is generated, 154 * if n and m are both numbers, an nxm matrix is generated. 155 * @param {Number} n Number of rows 156 * @param {Number} [m=n] Number of columns 157 * @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 158 * 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. 159 */ 160 identity: function (n, m) { 161 var r, i; 162 163 if ((m === undef) && (typeof m !== 'number')) { 164 m = n; 165 } 166 167 r = this.matrix(n, m); 168 169 for (i = 0; i < Math.min(n, m); i++) { 170 r[i][i] = 1; 171 } 172 173 return r; 174 }, 175 176 /** 177 * Generates a 4x4 matrix for 3D to 2D projections. 178 * @param {Number} l Left 179 * @param {Number} r Right 180 * @param {Number} t Top 181 * @param {Number} b Bottom 182 * @param {Number} n Near 183 * @param {Number} f Far 184 * @returns {Array} 4x4 Matrix 185 */ 186 frustum: function (l, r, b, t, n, f) { 187 var ret = this.matrix(4, 4); 188 189 ret[0][0] = (n * 2) / (r - l); 190 ret[0][1] = 0; 191 ret[0][2] = (r + l) / (r - l); 192 ret[0][3] = 0; 193 194 ret[1][0] = 0; 195 ret[1][1] = (n * 2) / (t - b); 196 ret[1][2] = (t + b) / (t - b); 197 ret[1][3] = 0; 198 199 ret[2][0] = 0; 200 ret[2][1] = 0; 201 ret[2][2] = -(f + n) / (f - n); 202 ret[2][3] = -(f * n * 2) / (f - n); 203 204 ret[3][0] = 0; 205 ret[3][1] = 0; 206 ret[3][2] = -1; 207 ret[3][3] = 0; 208 209 return ret; 210 }, 211 212 /** 213 * Generates a 4x4 matrix for 3D to 2D projections. 214 * @param {Number} fov Field of view in vertical direction, given in rad. 215 * @param {Number} ratio Aspect ratio of the projection plane. 216 * @param {Number} n Near 217 * @param {Number} f Far 218 * @returns {Array} 4x4 Projection Matrix 219 */ 220 projection: function (fov, ratio, n, f) { 221 var t = n * Math.tan(fov / 2), 222 r = t * ratio; 223 224 return this.frustum(-r, r, -t, t, n, f); 225 }, 226 227 /** 228 * 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 229 * function does not check if the dimensions match. 230 * @param {Array} mat Two dimensional array of numbers. The inner arrays describe the columns, the outer ones the matrix' rows. 231 * @param {Array} vec Array of numbers 232 * @returns {Array} Array of numbers containing the result 233 * @example 234 * var A = [[2, 1], 235 * [1, 3]], 236 * b = [4, 5], 237 * c; 238 * c = JXG.Math.matVecMult(A, b) 239 * // c === [13, 19]; 240 */ 241 matVecMult: function (mat, vec) { 242 var i, s, k, 243 m = mat.length, 244 n = vec.length, 245 res = []; 246 247 if (n === 3) { 248 for (i = 0; i < m; i++) { 249 res[i] = mat[i][0] * vec[0] + mat[i][1] * vec[1] + mat[i][2] * vec[2]; 250 } 251 } else { 252 for (i = 0; i < m; i++) { 253 s = 0; 254 for (k = 0; k < n; k++) { 255 s += mat[i][k] * vec[k]; 256 } 257 res[i] = s; 258 } 259 } 260 return res; 261 }, 262 263 /** 264 * Computes the product of the two matrices mat1*mat2. 265 * @param {Array} mat1 Two dimensional array of numbers 266 * @param {Array} mat2 Two dimensional array of numbers 267 * @returns {Array} Two dimensional Array of numbers containing result 268 */ 269 matMatMult: function (mat1, mat2) { 270 var i, j, s, k, 271 m = mat1.length, 272 n = m > 0 ? mat2[0].length : 0, 273 m2 = mat2.length, 274 res = this.matrix(m, n); 275 276 for (i = 0; i < m; i++) { 277 for (j = 0; j < n; j++) { 278 s = 0; 279 for (k = 0; k < m2; k++) { 280 s += mat1[i][k] * mat2[k][j]; 281 } 282 res[i][j] = s; 283 } 284 } 285 return res; 286 }, 287 288 /** 289 * Transposes a matrix given as a two dimensional array. 290 * @param {Array} M The matrix to be transposed 291 * @returns {Array} The transpose of M 292 */ 293 transpose: function (M) { 294 var MT, i, j, 295 m, n; 296 297 // number of rows of M 298 m = M.length; 299 // number of columns of M 300 n = M.length > 0 ? M[0].length : 0; 301 MT = this.matrix(n, m); 302 303 for (i = 0; i < n; i++) { 304 for (j = 0; j < m; j++) { 305 MT[i][j] = M[j][i]; 306 } 307 } 308 309 return MT; 310 }, 311 312 /** 313 * Compute the inverse of an nxn matrix with Gauss elimination. 314 * @param {Array} Ain 315 * @returns {Array} Inverse matrix of Ain 316 */ 317 inverse: function (Ain) { 318 var i, j, k, s, ma, r, swp, 319 n = Ain.length, 320 A = [], 321 p = [], 322 hv = []; 323 324 for (i = 0; i < n; i++) { 325 A[i] = []; 326 for (j = 0; j < n; j++) { 327 A[i][j] = Ain[i][j]; 328 } 329 p[i] = i; 330 } 331 332 for (j = 0; j < n; j++) { 333 // pivot search: 334 ma = Math.abs(A[j][j]); 335 r = j; 336 337 for (i = j + 1; i < n; i++) { 338 if (Math.abs(A[i][j]) > ma) { 339 ma = Math.abs(A[i][j]); 340 r = i; 341 } 342 } 343 344 // Singular matrix 345 if (ma <= this.eps) { 346 return []; 347 } 348 349 // swap rows: 350 if (r > j) { 351 for (k = 0; k < n; k++) { 352 swp = A[j][k]; 353 A[j][k] = A[r][k]; 354 A[r][k] = swp; 355 } 356 357 swp = p[j]; 358 p[j] = p[r]; 359 p[r] = swp; 360 } 361 362 // transformation: 363 s = 1.0 / A[j][j]; 364 for (i = 0; i < n; i++) { 365 A[i][j] *= s; 366 } 367 A[j][j] = s; 368 369 for (k = 0; k < n; k++) { 370 if (k !== j) { 371 for (i = 0; i < n; i++) { 372 if (i !== j) { 373 A[i][k] -= A[i][j] * A[j][k]; 374 } 375 } 376 A[j][k] = -s * A[j][k]; 377 } 378 } 379 } 380 381 // swap columns: 382 for (i = 0; i < n; i++) { 383 for (k = 0; k < n; k++) { 384 hv[p[k]] = A[i][k]; 385 } 386 for (k = 0; k < n; k++) { 387 A[i][k] = hv[k]; 388 } 389 } 390 391 return A; 392 }, 393 394 /** 395 * Inner product of two vectors a and b. n is the length of the vectors. 396 * @param {Array} a Vector 397 * @param {Array} b Vector 398 * @param {Number} [n] Length of the Vectors. If not given the length of the first vector is taken. 399 * @returns {Number} The inner product of a and b. 400 */ 401 innerProduct: function (a, b, n) { 402 var i, 403 s = 0; 404 405 if ((n === undef) || (typeof n !== 'number')) { 406 n = a.length; 407 } 408 409 for (i = 0; i < n; i++) { 410 s += a[i] * b[i]; 411 } 412 413 return s; 414 }, 415 416 /** 417 * Calculates the cross product of two vectors both of length three. 418 * In case of homogeneous coordinates this is either 419 * <ul> 420 * <li>the intersection of two lines</li> 421 * <li>the line through two points</li> 422 * </ul> 423 * @param {Array} c1 Homogeneous coordinates of line or point 1 424 * @param {Array} c2 Homogeneous coordinates of line or point 2 425 * @returns {Array} vector of length 3: homogeneous coordinates of the resulting point / line. 426 */ 427 crossProduct: function (c1, c2) { 428 return [c1[1] * c2[2] - c1[2] * c2[1], 429 c1[2] * c2[0] - c1[0] * c2[2], 430 c1[0] * c2[1] - c1[1] * c2[0]]; 431 }, 432 433 /** 434 * Compute the factorial of a positive integer. If a non-integer value 435 * is given, the fraction will be ignored. 436 * @function 437 * @param {Number} n 438 * @returns {Number} n! = n*(n-1)*...*2*1 439 */ 440 factorial: memoizer(function (n) { 441 if (n < 0) { 442 return NaN; 443 } 444 445 n = Math.floor(n); 446 447 if (n === 0 || n === 1) { 448 return 1; 449 } 450 451 return n * this.factorial(n - 1); 452 }), 453 454 /** 455 * Computes the binomial coefficient n over k. 456 * @function 457 * @param {Number} n Fraction will be ignored 458 * @param {Number} k Fraction will be ignored 459 * @returns {Number} The binomial coefficient n over k 460 */ 461 binomial: memoizer(function (n, k) { 462 var b, i; 463 464 if (k > n || k < 0) { 465 return NaN; 466 } 467 468 k = Math.round(k); 469 n = Math.round(n); 470 471 if (k === 0 || k === n) { 472 return 1; 473 } 474 475 b = 1; 476 477 for (i = 0; i < k; i++) { 478 b *= (n - i); 479 b /= (i + 1); 480 } 481 482 return b; 483 }), 484 485 /** 486 * Calculates the cosine hyperbolicus of x. 487 * @param {Number} x The number the cosine hyperbolicus will be calculated of. 488 * @returns {Number} Cosine hyperbolicus of the given value. 489 */ 490 cosh: function (x) { 491 return (Math.exp(x) + Math.exp(-x)) * 0.5; 492 }, 493 494 /** 495 * Sine hyperbolicus of x. 496 * @param {Number} x The number the sine hyperbolicus will be calculated of. 497 * @returns {Number} Sine hyperbolicus of the given value. 498 */ 499 sinh: function (x) { 500 return (Math.exp(x) - Math.exp(-x)) * 0.5; 501 }, 502 503 /** 504 * Compute base to the power of exponent. 505 * @param {Number} base 506 * @param {Number} exponent 507 * @returns {Number} base to the power of exponent. 508 */ 509 pow: function (base, exponent) { 510 if (base === 0) { 511 if (exponent === 0) { 512 return 1; 513 } 514 515 return 0; 516 } 517 518 if (Math.floor(exponent) === exponent) { 519 // a is an integer 520 return Math.pow(base, exponent); 521 } 522 523 // a is not an integer 524 if (base > 0) { 525 return Math.exp(exponent * Math.log(Math.abs(base))); 526 } 527 528 return NaN; 529 }, 530 531 /** 532 * Logarithm to base 10. 533 * @param {Number} x 534 * @returns {Number} log10(x) Logarithm of x to base 10. 535 */ 536 log10: function (x) { 537 return Math.log(x) / Math.log(10.0); 538 }, 539 540 /** 541 * Logarithm to base 2. 542 * @param {Number} x 543 * @returns {Number} log2(x) Logarithm of x to base 2. 544 */ 545 log2: function (x) { 546 return Math.log(x) / Math.log(2.0); 547 }, 548 549 /** 550 * Logarithm to arbitrary base b. If b is not given, natural log is taken, i.e. b = e. 551 * @param {Number} x 552 * @param {Number} b base 553 * @returns {Number} log(x, b) Logarithm of x to base b, that is log(x)/log(b). 554 */ 555 log: function (x, b) { 556 if (typeof b !== 'undefined' && Type.isNumber(b)) { 557 return Math.log(x) / Math.log(b); 558 } else { 559 return Math.log(x); 560 } 561 }, 562 563 /** 564 * A square & multiply algorithm to compute base to the power of exponent. 565 * Implementated by Wolfgang Riedl. 566 * @param {Number} base 567 * @param {Number} exponent 568 * @returns {Number} Base to the power of exponent 569 */ 570 squampow: function (base, exponent) { 571 var result; 572 573 if (Math.floor(exponent) === exponent) { 574 // exponent is integer (could be zero) 575 result = 1; 576 577 if (exponent < 0) { 578 // invert: base 579 base = 1.0 / base; 580 exponent *= -1; 581 } 582 583 while (exponent !== 0) { 584 if (exponent & 1) { 585 result *= base; 586 } 587 588 exponent >>= 1; 589 base *= base; 590 } 591 return result; 592 } 593 594 return this.pow(base, exponent); 595 }, 596 597 /** 598 * Normalize the standard form [c, b0, b1, a, k, r, q0, q1]. 599 * @private 600 * @param {Array} stdform The standard form to be normalized. 601 * @returns {Array} The normalized standard form. 602 */ 603 normalize: function (stdform) { 604 var n, signr, 605 a2 = 2 * stdform[3], 606 r = stdform[4] / a2; 607 608 stdform[5] = r; 609 stdform[6] = -stdform[1] / a2; 610 stdform[7] = -stdform[2] / a2; 611 612 if (!isFinite(r)) { 613 n = Math.sqrt(stdform[1] * stdform[1] + stdform[2] * stdform[2]); 614 615 stdform[0] /= n; 616 stdform[1] /= n; 617 stdform[2] /= n; 618 stdform[3] = 0; 619 stdform[4] = 1; 620 } else if (Math.abs(r) >= 1) { 621 stdform[0] = (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) / (2 * r); 622 stdform[1] = -stdform[6] / r; 623 stdform[2] = -stdform[7] / r; 624 stdform[3] = 1 / (2 * r); 625 stdform[4] = 1; 626 } else { 627 signr = (r <= 0 ? -1 : 1); 628 stdform[0] = signr * (stdform[6] * stdform[6] + stdform[7] * stdform[7] - r * r) * 0.5; 629 stdform[1] = -signr * stdform[6]; 630 stdform[2] = -signr * stdform[7]; 631 stdform[3] = signr / 2; 632 stdform[4] = signr * r; 633 } 634 635 return stdform; 636 }, 637 638 /** 639 * Converts a two dimensional array to a one dimensional Float32Array that can be processed by WebGL. 640 * @param {Array} m A matrix in a two dimensional array. 641 * @returns {Float32Array} A one dimensional array containing the matrix in column wise notation. Provides a fall 642 * back to the default JavaScript Array if Float32Array is not available. 643 */ 644 toGL: function (m) { 645 var v, i, j; 646 647 if (typeof Float32Array === 'function') { 648 v = new Float32Array(16); 649 } else { 650 v = new Array(16); 651 } 652 653 if (m.length !== 4 && m[0].length !== 4) { 654 return v; 655 } 656 657 for (i = 0; i < 4; i++) { 658 for (j = 0; j < 4; j++) { 659 v[i + 4 * j] = m[i][j]; 660 } 661 } 662 663 return v; 664 } 665 }; 666 667 return JXG.Math; 668 }); 669