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