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