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