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