1 /* 2 Copyright 2008-2021 3 Matthias Ehmann, 4 Carsten Miller, 5 Andreas Walter, 6 Alfred Wassermann 7 8 This file is part of JSXGraph. 9 10 JSXGraph is free software dual licensed under the GNU LGPL or MIT License. 11 12 You can redistribute it and/or modify it under the terms of the 13 14 * GNU Lesser General Public License as published by 15 the Free Software Foundation, either version 3 of the License, or 16 (at your option) any later version 17 OR 18 * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT 19 20 JSXGraph is distributed in the hope that it will be useful, 21 but WITHOUT ANY WARRANTY; without even the implied warranty of 22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 23 GNU Lesser General Public License for more details. 24 25 You should have received a copy of the GNU Lesser General Public License and 26 the MIT License along with JSXGraph. If not, see <http://www.gnu.org/licenses/> 27 and <http://opensource.org/licenses/MIT/>. 28 */ 29 30 31 /*global JXG: true, define: true*/ 32 /*jslint nomen: true, plusplus: true*/ 33 34 /* depends: 35 jxg 36 math/math 37 utils/type 38 */ 39 40 define(['math/math', 'utils/type'], function (Mat, Type) { 41 42 "use strict"; 43 44 /** 45 * Probability functions, e.g. error function, 46 * see: https://en.wikipedia.org/wiki/Error_function 47 * Ported from 48 * by https://github.com/jeremybarnes/cephes/blob/master/cprob/ndtr.c, 49 * 50 * Cephes Math Library Release 2.9: November, 2000 51 * Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier 52 * 53 * @name JXG.Math.ProbFuncs 54 * @exports Mat.ProbFuncs as JXG.Math.ProbFuncs 55 * @namespace 56 */ 57 Mat.ProbFuncs = { 58 MAXNUM: 1.701411834604692317316873e38, // 2**127 59 SQRTH: 7.07106781186547524401E-1, // sqrt(2)/2 60 SQRT2: 1.41421356237309504880, // sqrt(2) 61 MAXLOG: 7.08396418532264106224E2, // log 2**1022 62 63 P: [ 64 2.46196981473530512524E-10, 65 5.64189564831068821977E-1, 66 7.46321056442269912687E0, 67 4.86371970985681366614E1, 68 1.96520832956077098242E2, 69 5.26445194995477358631E2, 70 9.34528527171957607540E2, 71 1.02755188689515710272E3, 72 5.57535335369399327526E2 73 ], 74 75 Q: [ 76 1.32281951154744992508E1, 77 8.67072140885989742329E1, 78 3.54937778887819891062E2, 79 9.75708501743205489753E2, 80 1.82390916687909736289E3, 81 2.24633760818710981792E3, 82 1.65666309194161350182E3, 83 5.57535340817727675546E2 84 ], 85 86 R: [ 87 5.64189583547755073984E-1, 88 1.27536670759978104416E0, 89 5.01905042251180477414E0, 90 6.16021097993053585195E0, 91 7.40974269950448939160E0, 92 2.97886665372100240670E0 93 ], 94 95 S: [ 96 2.26052863220117276590E0, 97 9.39603524938001434673E0, 98 1.20489539808096656605E1, 99 1.70814450747565897222E1, 100 9.60896809063285878198E0, 101 3.36907645100081516050E0 102 ], 103 104 T: [ 105 9.60497373987051638749E0, 106 9.00260197203842689217E1, 107 2.23200534594684319226E3, 108 7.00332514112805075473E3, 109 5.55923013010394962768E4 110 ], 111 112 U: [ 113 3.35617141647503099647E1, 114 5.21357949780152679795E2, 115 4.59432382970980127987E3, 116 2.26290000613890934246E4, 117 4.92673942608635921086E4 118 ], 119 120 // UTHRESH: 37.519379347, 121 M: 128.0, 122 MINV: 0.0078125, 123 124 /** 125 * 126 * Exponential of squared argument 127 * 128 * SYNOPSIS: 129 * 130 * double x, y, expx2(); 131 * int sign; 132 * 133 * y = expx2( x, sign ); 134 * 135 * 136 * 137 * DESCRIPTION: 138 * 139 * Computes y = exp(x*x) while suppressing error amplification 140 * that would ordinarily arise from the inexactness of the 141 * exponential argument x*x. 142 * 143 * If sign < 0, the result is inverted; i.e., y = exp(-x*x) . 144 * 145 * 146 * ACCURACY: 147 * 148 * Relative error: 149 * arithmetic domain # trials peak rms 150 * IEEE -26.6, 26.6 10^7 3.9e-16 8.9e-17 151 * 152 * @private 153 * @param {Number} x 154 * @param {Number} sign (int) 155 * @returns {Number} 156 */ 157 expx2: function(x, sign) { 158 // double x; 159 // int sign; 160 var u, u1, m, f; 161 162 x = Math.abs(x); 163 if (sign < 0) { 164 x = -x; 165 } 166 167 // Represent x as an exact multiple of M plus a residual. 168 // M is a power of 2 chosen so that exp(m * m) does not overflow 169 // or underflow and so that |x - m| is small. 170 m = this.MINV * Math.floor(this.M * x + 0.5); 171 f = x - m; 172 173 // x^2 = m^2 + 2mf + f^2 174 u = m * m; 175 u1 = 2 * m * f + f * f; 176 177 if (sign < 0) { 178 u = -u; 179 u1 = -u1; 180 } 181 182 if ( u + u1 > this.MAXLOG) { 183 return Infinity; 184 } 185 186 // u is exact, u1 is small. 187 u = Math.exp(u) * Math.exp(u1); 188 return u; 189 }, 190 191 /** 192 * 193 * Evaluate polynomial 194 * 195 * SYNOPSIS: 196 * 197 * int N; 198 * double x, y, coef[N+1], polevl[]; 199 * 200 * y = polevl( x, coef, N ); 201 * 202 * DESCRIPTION: 203 * 204 * Evaluates polynomial of degree N: 205 * 206 * 2 N 207 * y = C + C x + C x +...+ C x 208 * 0 1 2 N 209 * 210 * Coefficients are stored in reverse order: 211 * 212 * coef[0] = C , ..., coef[N] = C . 213 * N 0 214 * 215 * The function p1evl() assumes that coef[N] = 1.0 and is 216 * omitted from the array. Its calling arguments are 217 * otherwise the same as polevl(). 218 * 219 * 220 * SPEED: 221 * 222 * In the interest of speed, there are no checks for out 223 * of bounds arithmetic. This routine is used by most of 224 * the functions in the library. Depending on available 225 * equipment features, the user may wish to rewrite the 226 * program in microcode or assembly language. 227 * 228 * @private 229 * @param {Number} x 230 * @param {Number} coef 231 * @param {Number} N 232 * @returns {Number} 233 */ 234 polevl: function(x, coef, N) { 235 var ans, i; 236 237 if (Type.exists(coef.reduce)) { 238 return coef.reduce(function(acc, c) { 239 return acc * x + c; 240 }, 0); 241 } 242 // Polyfill 243 for (i = 0, ans = 0; i <= N; i++) { 244 ans = ans * x + coef[i]; 245 } 246 return ans; 247 248 }, 249 250 /** 251 * Evaluate polynomial when coefficient of x is 1.0. 252 * Otherwise same as polevl. 253 * 254 * @private 255 * @param {Number} x 256 * @param {Number} coef 257 * @param {Number} N 258 * @returns {Number} 259 */ 260 p1evl: function(x, coef, N) { 261 var ans, i; 262 263 if (Type.exists(coef.reduce)) { 264 return coef.reduce(function(acc, c) { 265 return acc * x + c; 266 }, 1); 267 } 268 // Polyfill 269 for (i = 0, ans = 1; i < N; i++) { 270 ans = ans * x + coef[i]; 271 } 272 return ans; 273 }, 274 275 /** 276 * 277 * Normal distribution function 278 * 279 * SYNOPSIS: 280 * 281 * y = ndtr( x ); 282 * 283 * DESCRIPTION: 284 * 285 * Returns the area under the Gaussian probability density 286 * function, integrated from minus infinity to x: 287 * 288 * x 289 * - 290 * 1 | | 2 291 * ndtr(x) = --------- | exp( - t /2 ) dt 292 * sqrt(2pi) | | 293 * - 294 * -inf. 295 * 296 * = ( 1 + erf(z) ) / 2 297 * = erfc(z) / 2 298 * 299 * where z = x/sqrt(2). Computation is via the functions 300 * erf and erfc with care to avoid error amplification in computing exp(-x^2). 301 * 302 * 303 * ACCURACY: 304 * 305 * Relative error: 306 * arithmetic domain # trials peak rms 307 * IEEE -13,0 30000 1.3e-15 2.2e-16 308 * 309 * 310 * ERROR MESSAGES: 311 * 312 * message condition value returned 313 * erfc underflow x > 37.519379347 0.0 314 * 315 * @param {Number} a 316 * @returns {Number} 317 */ 318 ndtr: function(a) { 319 // a: double, return double 320 var x, y, z; 321 322 x = a * this.SQRTH; 323 z = Math.abs(x); 324 325 if (z < 1.0) { 326 y = 0.5 + 0.5 * this.erf(x); 327 } else { 328 y = 0.5 * this.erfce(z); 329 /* Multiply by exp(-x^2 / 2) */ 330 z = this.expx2(a, -1); 331 y = y * Math.sqrt(z); 332 if (x > 0) { 333 y = 1.0 - y; 334 } 335 } 336 return y; 337 }, 338 339 /** 340 * @private 341 * @param {Number} a 342 * @returns {Number} 343 */ 344 _underflow: function(a) { 345 console.log('erfc', 'UNDERFLOW'); 346 if (a < 0) { 347 return 2.0; 348 } 349 return 0.0; 350 }, 351 352 /** 353 * 354 * Complementary error function 355 * 356 * SYNOPSIS: 357 * 358 * double x, y, erfc(); 359 * 360 * y = erfc( x ); 361 * 362 * 363 * 364 * DESCRIPTION: 365 * 366 * 367 * 1 - erf(x) = 368 * 369 * inf. 370 * - 371 * 2 | | 2 372 * erfc(x) = -------- | exp( - t ) dt 373 * sqrt(pi) | | 374 * - 375 * x 376 * 377 * 378 * For small x, erfc(x) = 1 - erf(x); otherwise rational 379 * approximations are computed. 380 * 381 * A special function expx2.c is used to suppress error amplification 382 * in computing exp(-x^2). 383 * 384 * 385 * ACCURACY: 386 * 387 * Relative error: 388 * arithmetic domain # trials peak rms 389 * IEEE 0,26.6417 30000 1.3e-15 2.2e-16 390 * 391 * 392 * ERROR MESSAGES: 393 * 394 * message condition value returned 395 * erfc underflow x > 9.231948545 (DEC) 0.0 396 * 397 * @param {Number} a 398 * @returns {Number} 399 */ 400 erfc: function(a) { 401 var p, q, x, y, z; 402 403 if (a < 0.0) { 404 x = -a; 405 } else { 406 x = a; 407 } 408 if (x < 1.0) { 409 return 1.0 - this.erf(a); 410 } 411 412 z = -a * a; 413 if (z < -this.MAXLOG) { 414 return this._underflow(a); 415 } 416 417 z = this.expx2(a, -1); // Compute z = exp(z). 418 419 if (x < 8.0) { 420 p = this.polevl(x, this.P, 8); 421 q = this.p1evl(x, this.Q, 8); 422 } else { 423 p = this.polevl(x, this.R, 5); 424 q = this.p1evl(x, this.S, 6); 425 } 426 427 y = (z * p) / q; 428 429 if (a < 0) { 430 y = 2.0 - y; 431 } 432 433 if (y === 0.0) { 434 return this._underflow(a); 435 } 436 437 return y; 438 }, 439 440 /** 441 * Exponentially scaled erfc function 442 * exp(x^2) erfc(x) 443 * valid for x > 1. 444 * Use with ndtr and expx2. 445 * 446 * @private 447 * @param {Number} x 448 * @returns {Number} 449 */ 450 erfce: function(x) { 451 var p, q; 452 453 if (x < 8.0) { 454 p = this.polevl(x, this.P, 8); 455 q = this.p1evl(x, this.Q, 8); 456 } else { 457 p = this.polevl( x, this.R, 5 ); 458 q = this.p1evl( x, this.S, 6 ); 459 } 460 return p / q; 461 }, 462 463 /** 464 * Error function 465 * 466 * SYNOPSIS: 467 * 468 * double x, y, erf(); 469 * 470 * y = erf( x ); 471 * 472 * 473 * 474 * DESCRIPTION: 475 * 476 * The integral is 477 * 478 * x 479 * - 480 * 2 | | 2 481 * erf(x) = -------- | exp( - t ) dt. 482 * sqrt(pi) | | 483 * - 484 * 0 485 * 486 * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise 487 * erf(x) = 1 - erfc(x). 488 * 489 * 490 * ACCURACY: 491 * 492 * Relative error: 493 * arithmetic domain # trials peak rms 494 * DEC 0,1 14000 4.7e-17 1.5e-17 495 * IEEE 0,1 30000 3.7e-16 1.0e-16 496 * 497 * @param {Number} x 498 * @returns {Number} 499 */ 500 erf: function(x) { 501 var y, z; 502 503 if (Math.abs(x) > 1.0) { 504 return 1.0 - this.erfc(x); 505 } 506 z = x * x; 507 y = x * this.polevl(z, this.T, 4) / this.p1evl(z, this.U, 5); 508 return y; 509 }, 510 511 s2pi: 2.50662827463100050242E0, // sqrt(2pi) 512 513 // approximation for 0 <= |y - 0.5| <= 3/8 */ 514 P0: [ 515 -5.99633501014107895267E1, 516 9.80010754185999661536E1, 517 -5.66762857469070293439E1, 518 1.39312609387279679503E1, 519 -1.23916583867381258016E0 520 ], 521 522 Q0: [ 523 1.95448858338141759834E0, 524 4.67627912898881538453E0, 525 8.63602421390890590575E1, 526 -2.25462687854119370527E2, 527 2.00260212380060660359E2, 528 -8.20372256168333339912E1, 529 1.59056225126211695515E1, 530 -1.18331621121330003142E0, 531 ], 532 533 // Approximation for interval z = sqrt(-2 log y ) between 2 and 8 534 // i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14. 535 P1: [ 536 4.05544892305962419923E0, 537 3.15251094599893866154E1, 538 5.71628192246421288162E1, 539 4.40805073893200834700E1, 540 1.46849561928858024014E1, 541 2.18663306850790267539E0, 542 -1.40256079171354495875E-1, 543 -3.50424626827848203418E-2, 544 -8.57456785154685413611E-4 545 ], 546 547 Q1: [ 548 1.57799883256466749731E1, 549 4.53907635128879210584E1, 550 4.13172038254672030440E1, 551 1.50425385692907503408E1, 552 2.50464946208309415979E0, 553 -1.42182922854787788574E-1, 554 -3.80806407691578277194E-2, 555 -9.33259480895457427372E-4 556 ], 557 558 // Approximation for interval z = sqrt(-2 log y ) between 8 and 64 559 // i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890. 560 P2: [ 561 3.23774891776946035970E0, 562 6.91522889068984211695E0, 563 3.93881025292474443415E0, 564 1.33303460815807542389E0, 565 2.01485389549179081538E-1, 566 1.23716634817820021358E-2, 567 3.01581553508235416007E-4, 568 2.65806974686737550832E-6, 569 6.23974539184983293730E-9 570 ], 571 572 Q2: [ 573 6.02427039364742014255E0, 574 3.67983563856160859403E0, 575 1.37702099489081330271E0, 576 2.16236993594496635890E-1, 577 1.34204006088543189037E-2, 578 3.28014464682127739104E-4, 579 2.89247864745380683936E-6, 580 6.79019408009981274425E-9 581 ], 582 583 /** 584 * 585 * Inverse of Normal distribution function 586 * 587 * SYNOPSIS: 588 * 589 * double x, y, ndtri(); 590 * 591 * x = ndtri( y ); 592 * 593 * DESCRIPTION: 594 * 595 * Returns the argument, x, for which the area under the 596 * Gaussian probability density function (integrated from 597 * minus infinity to x) is equal to y. 598 * 599 * 600 * For small arguments 0 < y < exp(-2), the program computes 601 * z = sqrt( -2.0 * log(y) ); then the approximation is 602 * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). 603 * There are two rational functions P/Q, one for 0 < y < exp(-32) 604 * and the other for y up to exp(-2). For larger arguments, 605 * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). 606 * 607 * 608 * ACCURACY: 609 * 610 * Relative error: 611 * arithmetic domain # trials peak rms 612 * DEC 0.125, 1 5500 9.5e-17 2.1e-17 613 * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17 614 * IEEE 0.125, 1 20000 7.2e-16 1.3e-16 615 * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 616 * 617 * 618 * ERROR MESSAGES: 619 * 620 * message condition value returned 621 * ndtri domain x <= 0 -MAXNUM 622 * ndtri domain x >= 1 MAXNUM 623 * 624 * @param {Number} y0 625 * @returns {Number} 626 */ 627 ndtri: function(y0) { 628 var x, y, z, y2, x0, x1, code; 629 630 if (y0 <= 0.0) { 631 //console.log("ndtri", "DOMAIN "); 632 return -Infinity; // -this.MAXNUM; 633 } 634 if (y0 >= 1.0) { 635 // console.log("ndtri", "DOMAIN"); 636 return Infinity; // this.MAXNUM; 637 } 638 639 code = 1; 640 y = y0; 641 if (y > (1.0 - 0.13533528323661269189)) { // 0.135... = exp(-2) 642 y = 1.0 - y; 643 code = 0; 644 } 645 646 if (y > 0.13533528323661269189) { 647 y = y - 0.5; 648 y2 = y * y; 649 x = y + y * (y2 * this.polevl(y2, this.P0, 4) / this.p1evl(y2, this.Q0, 8)); 650 x = x * this.s2pi; 651 return x; 652 } 653 654 x = Math.sqrt( -2.0 * Math.log(y) ); 655 x0 = x - Math.log(x) / x; 656 657 z = 1.0 / x; 658 if (x < 8.0) { // y > exp(-32) = 1.2664165549e-14 659 x1 = z * this.polevl(z, this.P1, 8 ) / this.p1evl(z, this.Q1, 8); 660 } else { 661 x1 = z * this.polevl(z, this.P2, 8) / this.p1evl(z, this.Q2, 8); 662 } 663 x = x0 - x1; 664 if (code !== 0) { 665 x = -x; 666 } 667 return x; 668 }, 669 670 /** 671 * Inverse of error function erf. 672 * 673 * @param {Number} x 674 * @returns {Number} 675 */ 676 erfi: function(x) { 677 return this.ndtri((x + 1) * 0.5) * this.SQRTH; 678 } 679 }; 680 681 return Mat.ProbFuncs; 682 }); 683