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