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 * Initializes a vector as an array with the coefficients set to the given value resp. zero. 77 * @param {Number} n Length of the vector 78 * @param {Number} [init=0] Initial value for each coefficient 79 * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a 80 * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows. 81 */ 82 vector: function(n, init) { 83 var r, i; 84 85 init = init || 0; 86 87 r = new Array(Math.ceil(n)); 88 for(i=0; i<n; i++) { 89 r[i] = init; 90 } 91 92 return r; 93 }, 94 95 /** 96 * Initializes a matrix as an array of rows with the given value. 97 * @param {Number} n Number of rows 98 * @param {Number} [m=n] Number of columns 99 * @param {Number} [init=0] Initial value for each coefficient 100 * @returns {Array} A <tt>n</tt> times <tt>m</tt>-matrix represented by a 101 * two-dimensional array. The inner arrays hold the columns, the outer array holds the rows. 102 */ 103 matrix: function(n, m, init) { 104 var r, i, j; 105 106 init = init || 0; 107 m = m || n; 108 109 r = new Array(Math.ceil(n)); 110 for(i=0; i<n; i++) { 111 r[i] = new Array(Math.ceil(m)); 112 for(j=0; j<m; j++) { 113 r[i][j] = init; 114 } 115 } 116 117 return r; 118 }, 119 120 /** 121 * Generates an identity matrix. If n is a number and m is undefined or not a number, a square matrix is generated, 122 * if n and m are both numbers, an nxm matrix is generated. 123 * @param {Number} n Number of rows 124 * @param {Number} [m=n] Number of columns 125 * @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 126 * 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. 127 */ 128 identity: function(n, m) { 129 var r, i; 130 131 if((m === undefined) && (typeof m !== 'number')) { 132 m = n; 133 } 134 135 r = this.matrix(n, m); 136 for(i=0; i<Math.min(n, m); i++) { 137 r[i][i] = 1; 138 } 139 140 return r; 141 }, 142 143 /** 144 * 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 145 * function does not check if the dimensions match. 146 * @param {Array} mat Two dimensional array of numbers. The inner arrays describe the columns, the outer ones the matrix' rows. 147 * @param {Array} vec Array of numbers 148 * @returns {Array} Array of numbers containing the result 149 * @example 150 * var A = [[2, 1], 151 * [1, 3]], 152 * b = [4, 5], 153 * c; 154 * c = JXG.Math.matVecMult(A, b) 155 * // c === [13, 19]; 156 */ 157 matVecMult: function(mat, vec) { 158 var m = mat.length, 159 n = vec.length, 160 res = [], 161 i, s, k; 162 163 if (n===3) { 164 for (i=0;i<m;i++) { 165 res[i] = mat[i][0]*vec[0] + mat[i][1]*vec[1] + mat[i][2]*vec[2]; 166 } 167 } else { 168 for (i=0;i<m;i++) { 169 s = 0; 170 for (k=0;k<n;k++) { s += mat[i][k]*vec[k]; } 171 res[i] = s; 172 } 173 } 174 return res; 175 }, 176 177 /** 178 * Computes the product of the two matrices mat1*mat2. 179 * @param {Array} mat1 Two dimensional array of numbers 180 * @param {Array} mat2 Two dimensional array of numbers 181 * @returns {Array} Two dimensional Array of numbers containing result 182 */ 183 matMatMult: function(mat1, mat2) { 184 var m = mat1.length, 185 n = m>0 ? mat2[0].length : 0, 186 m2 = mat2.length, 187 res = this.matrix(m,n), 188 i, j, s, k; 189 190 for (i=0;i<m;i++) { 191 for (j=0;j<n;j++) { 192 s = 0; 193 for (k=0;k<m2;k++) { 194 s += mat1[i][k]*mat2[k][j]; 195 } 196 res[i][j] = s; 197 } 198 } 199 return res; 200 }, 201 202 /** 203 * Transposes a matrix given as a two dimensional array. 204 * @param {Array} M The matrix to be transposed 205 * @returns {Array} The transpose of M 206 */ 207 transpose: function(M) { 208 var MT, i, j, 209 m, n; 210 211 m = M.length; // number of rows of M 212 n = M.length>0 ? M[0].length : 0; // number of columns of M 213 MT = this.matrix(n,m); 214 215 for (i=0; i<n; i++) { 216 for (j=0;j<m;j++) { 217 MT[i][j] = M[j][i]; 218 } 219 } 220 return MT; 221 }, 222 223 /** 224 * Compute the inverse of an nxn matrix with Gauss elimination. 225 * @param {Array} Ain 226 * @returns {Array} Inverse matrix of Ain 227 */ 228 inverse: function(Ain) { 229 var i,j,k,s,ma,r,swp, 230 n = Ain.length, 231 A = [], 232 p = [], 233 hv = []; 234 235 for (i=0;i<n;i++) { 236 A[i] = []; 237 for (j=0;j<n;j++) { A[i][j] = Ain[i][j]; } 238 p[i] = i; 239 } 240 241 for (j=0;j<n;j++) { 242 // pivot search: 243 ma = Math.abs(A[j][j]); 244 r = j; 245 for (i=j+1;i<n;i++) { 246 if (Math.abs(A[i][j])>ma) { 247 ma = Math.abs(A[i][j]); 248 r = i; 249 } 250 } 251 if (ma<=JXG.Math.eps) { // Singular matrix 252 return false; 253 } 254 // swap rows: 255 if (r>j) { 256 for (k=0;k<n;k++) { 257 swp = A[j][k]; A[j][k] = A[r][k]; A[r][k] = swp; 258 } 259 swp = p[j]; p[j] = p[r]; p[r] = swp; 260 } 261 // transformation: 262 s = 1.0/A[j][j]; 263 for (i=0;i<n;i++) { 264 A[i][j] *= s; 265 } 266 A[j][j] = s; 267 for (k=0;k<n;k++) if (k!=j) { 268 for (i=0;i<n;i++) if (i!=j) { 269 A[i][k] -= A[i][j]*A[j][k]; 270 } 271 A[j][k] = -s*A[j][k]; 272 } 273 } 274 // swap columns: 275 for (i=0;i<n;i++) { 276 for (k=0;k<n;k++) { hv[p[k]] = A[i][k]; } 277 for (k=0;k<n;k++) { A[i][k] = hv[k]; } 278 } 279 return A; 280 }, 281 282 /** 283 * Inner product of two vectors a and b. n is the length of the vectors. 284 * @param {Array} a Vector 285 * @param {Array} b Vector 286 * @param {Number} [n] Length of the Vectors. If not given the length of the first vector is taken. 287 * @return The inner product of a and b. 288 */ 289 innerProduct: function(a, b, n) { 290 var i, s = 0; 291 292 if((n === undefined) || (typeof n !== 'number')) { 293 n = a.length; 294 } 295 296 for (i=0; i<n; i++) { 297 s += a[i]*b[i]; 298 } 299 300 return s; 301 }, 302 303 /** 304 * Calculates the cross product of two vectors both of length three. 305 * In case of homogeneous coordinates this is either 306 * <ul> 307 * <li>the intersection of two lines</li> 308 * <li>the line through two points</li> 309 * </ul> 310 * @param {Array} c1 Homogeneous coordinates of line or point 1 311 * @param {Array} c2 Homogeneous coordinates of line or point 2 312 * @returns {Array} vector of length 3: homogeneous coordinates of the resulting point / line. 313 */ 314 crossProduct: function(c1,c2) { 315 return [c1[1]*c2[2]-c1[2]*c2[1], 316 c1[2]*c2[0]-c1[0]*c2[2], 317 c1[0]*c2[1]-c1[1]*c2[0]]; 318 }, 319 320 /** 321 * Compute the factorial of a positive integer. If a non-integer value 322 * is given, the fraction will be ignored. 323 * @function 324 * @param {Number} n 325 * @returns {Number} n! = n*(n-1)*...*2*1 326 */ 327 factorial: memoizer(function (n) { 328 if (n<0) return NaN; 329 n = Math.floor(n); 330 if (n===0 || n===1) return 1; 331 return n*arguments.callee(n-1); 332 }), 333 334 /** 335 * Computes the binomial coefficient n over k. 336 * @function 337 * @param {Number} n Fraction will be ignored 338 * @param {Number} k Fraction will be ignored 339 * @returns {Number} The binomial coefficient n over k 340 */ 341 binomial: memoizer(function(n,k) { 342 var b, i; 343 344 if (k>n || k<0) return NaN; 345 346 k = Math.floor(k); 347 n = Math.floor(n); 348 349 if (k===0 || k===n) return 1; 350 351 b = 1; 352 for (i=0; i<k; i++) { 353 b *= (n-i); 354 b /= (i+1); 355 } 356 return b; 357 }), 358 359 /** 360 * Calculates the cosine hyperbolicus of x. 361 * @param {Number} x The number the cosine hyperbolicus will be calculated of. 362 * @returns {Number} Cosine hyperbolicus of the given value. 363 */ 364 cosh: function(x) { 365 return (Math.exp(x)+Math.exp(-x))*0.5; 366 }, 367 368 /** 369 * Sine hyperbolicus of x. 370 * @param {Number} x The number the sine hyperbolicus will be calculated of. 371 * @returns {Number} Sine hyperbolicus of the given value. 372 */ 373 sinh: function(x) { 374 return (Math.exp(x)-Math.exp(-x))*0.5; 375 }, 376 377 /** 378 * Compute base to the power of exponent. 379 * @param {Number} base 380 * @param {Number} exponent 381 * @returns {Number} base to the power of exponent. 382 */ 383 pow: function(base, exponent) { 384 if (base===0) { 385 if (exponent===0) { 386 return 1; 387 } else { 388 return 0; 389 } 390 } 391 392 if (Math.floor(exponent)===exponent) {// a is an integer 393 return Math.pow(base,exponent); 394 } else { // a is not an integer 395 if (base>0) { 396 return Math.exp(exponent*Math.log(Math.abs(base))); 397 } else { 398 return NaN; 399 } 400 } 401 }, 402 403 /** 404 * A square & multiply algorithm to compute base to the power of exponent. 405 * Implementated by Wolfgang Riedl. 406 * @param {Number} base 407 * @param {Number} exponent 408 * @returns {Number} Base to the power of exponent 409 */ 410 squampow: function(base, exponent) { 411 var result; 412 413 if (Math.floor(exponent)===exponent) { // exponent is integer (could be zero) 414 result = 1; 415 if(exponent < 0) { 416 // invert: base 417 base = 1.0 / base; 418 exponent *= -1; 419 } 420 while(exponent != 0) { 421 if(exponent & 1) 422 result *= base; 423 exponent >>= 1; 424 base *= base; 425 } 426 return result; 427 } else { 428 return this.pow(base, exponent); 429 } 430 }, 431 432 /** 433 * Normalize the standard form [c, b0, b1, a, k, r, q0, q1]. 434 * @private 435 * @param {Array} stdform The standard form to be normalized. 436 * @returns {Array} The normalized standard form. 437 */ 438 normalize: function(stdform) { 439 var a2 = 2*stdform[3], 440 r = stdform[4]/(a2), // k/(2a) 441 n, signr; 442 stdform[5] = r; 443 stdform[6] = -stdform[1]/a2; 444 stdform[7] = -stdform[2]/a2; 445 if (r===Infinity || isNaN(r)) { 446 n = Math.sqrt(stdform[1]*stdform[1]+stdform[2]*stdform[2]); 447 stdform[0] /= n; 448 stdform[1] /= n; 449 stdform[2] /= n; 450 stdform[3] = 0; 451 stdform[4] = 1; 452 } else if (Math.abs(r)>=1) { 453 stdform[0] = (stdform[6]*stdform[6]+stdform[7]*stdform[7]-r*r)/(2*r); 454 stdform[1] = -stdform[6]/r; 455 stdform[2] = -stdform[7]/r; 456 stdform[3] = 1/(2*r); 457 stdform[4] = 1; 458 } else { 459 signr = (r<=0)?(-1):(1/*(r==0)?0:1*/); 460 stdform[0] = signr*(stdform[6]*stdform[6]+stdform[7]*stdform[7]-r*r)*0.5; 461 stdform[1] = -signr*stdform[6]; 462 stdform[2] = -signr*stdform[7]; 463 stdform[3] = signr/2; 464 stdform[4] = signr*r; 465 } 466 return stdform; 467 } 468 }; // JXG.Math 469 })(JXG, Math); 470