1 /*
  2     Copyright 2008-2024
  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         console.log("erfc", "UNDERFLOW");
313         if (a < 0) {
314             return 2.0;
315         }
316         return 0.0;
317     },
318 
319     /**
320      *
321      *	Complementary error function
322      *
323      * SYNOPSIS:
324      *
325      * double x, y, erfc();
326      *
327      * y = erfc( x );
328      *
329      *
330      *
331      * DESCRIPTION:
332      *
333      *
334      *  1 - erf(x) =
335      *
336      *                           inf.
337      *                             -
338      *                  2         | |          2
339      *   erfc(x)  =  --------     |    exp( - t  ) dt
340      *               sqrt(pi)   | |
341      *                           -
342      *                            x
343      *
344      *
345      * For small x, erfc(x) = 1 - erf(x); otherwise rational
346      * approximations are computed.
347      *
348      * A special function expx2.c is used to suppress error amplification
349      * in computing exp(-x^2).
350      *
351      *
352      * ACCURACY:
353      *
354      *                      Relative error:
355      * arithmetic   domain     # trials      peak         rms
356      *    IEEE      0,26.6417   30000       1.3e-15     2.2e-16
357      *
358      *
359      * ERROR MESSAGES:
360      *
361      *   message         condition              value returned
362      * erfc underflow    x > 9.231948545 (DEC)       0.0
363      *
364      * @param  {Number} a
365      * @returns {Number}
366      */
367     erfc: function (a) {
368         var p, q, x, y, z;
369 
370         if (a < 0.0) {
371             x = -a;
372         } else {
373             x = a;
374         }
375         if (x < 1.0) {
376             return 1.0 - this.erf(a);
377         }
378 
379         z = -a * a;
380         if (z < -this.MAXLOG) {
381             return this._underflow(a);
382         }
383 
384         z = this.expx2(a, -1); // Compute z = exp(z).
385 
386         if (x < 8.0) {
387             p = this.polevl(x, this.P, 8);
388             q = this.p1evl(x, this.Q, 8);
389         } else {
390             p = this.polevl(x, this.R, 5);
391             q = this.p1evl(x, this.S, 6);
392         }
393 
394         y = (z * p) / q;
395 
396         if (a < 0) {
397             y = 2.0 - y;
398         }
399 
400         if (y === 0.0) {
401             return this._underflow(a);
402         }
403 
404         return y;
405     },
406 
407     /**
408      * Exponentially scaled erfc function
409      *   exp(x^2) erfc(x)
410      *   valid for x > 1.
411      *   Use with ndtr and expx2.
412      *
413      * @private
414      * @param {Number} x
415      * @returns {Number}
416      */
417     erfce: function (x) {
418         var p, q;
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         return p / q;
428     },
429 
430     /**
431      *	Error function
432      *
433      * SYNOPSIS:
434      *
435      * double x, y, erf();
436      *
437      * y = erf( x );
438      *
439      *
440      *
441      * DESCRIPTION:
442      *
443      * The integral is
444      *
445      *                           x
446      *                            -
447      *                 2         | |          2
448      *   erf(x)  =  --------     |    exp( - t  ) dt.
449      *              sqrt(pi)   | |
450      *                          -
451      *                           0
452      *
453      * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
454      * erf(x) = 1 - erfc(x).
455      *
456      *
457      * ACCURACY:
458      *
459      *                      Relative error:
460      * arithmetic   domain     # trials      peak         rms
461      *    DEC       0,1         14000       4.7e-17     1.5e-17
462      *    IEEE      0,1         30000       3.7e-16     1.0e-16
463      *
464      * @param  {Number} x
465      * @returns {Number}
466      */
467     erf: function (x) {
468         var y, z;
469 
470         if (Math.abs(x) > 1.0) {
471             return 1.0 - this.erfc(x);
472         }
473         z = x * x;
474         y = (x * this.polevl(z, this.T, 4)) / this.p1evl(z, this.U, 5);
475         return y;
476     },
477 
478     s2pi: 2.50662827463100050242, // sqrt(2pi)
479 
480     // approximation for 0 <= |y - 0.5| <= 3/8 */
481     P0: [
482         -5.99633501014107895267e1, 9.80010754185999661536e1, -5.66762857469070293439e1,
483         1.39312609387279679503e1, -1.23916583867381258016
484     ],
485 
486     Q0: [
487         1.95448858338141759834, 4.67627912898881538453, 8.63602421390890590575e1,
488         -2.25462687854119370527e2, 2.00260212380060660359e2, -8.20372256168333339912e1,
489         1.59056225126211695515e1, -1.18331621121330003142
490     ],
491 
492     //  Approximation for interval z = sqrt(-2 log y ) between 2 and 8
493     //  i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
494     P1: [
495         4.05544892305962419923, 3.15251094599893866154e1, 5.71628192246421288162e1,
496         4.408050738932008347e1, 1.46849561928858024014e1, 2.18663306850790267539,
497         -1.40256079171354495875e-1, -3.50424626827848203418e-2, -8.57456785154685413611e-4
498     ],
499 
500     Q1: [
501         1.57799883256466749731e1, 4.53907635128879210584e1, 4.1317203825467203044e1,
502         1.50425385692907503408e1, 2.50464946208309415979, -1.42182922854787788574e-1,
503         -3.80806407691578277194e-2, -9.33259480895457427372e-4
504     ],
505 
506     // Approximation for interval z = sqrt(-2 log y ) between 8 and 64
507     // i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
508     P2: [
509         3.2377489177694603597, 6.91522889068984211695, 3.93881025292474443415,
510         1.33303460815807542389, 2.01485389549179081538e-1, 1.23716634817820021358e-2,
511         3.01581553508235416007e-4, 2.65806974686737550832e-6, 6.2397453918498329373e-9
512     ],
513 
514     Q2: [
515         6.02427039364742014255, 3.67983563856160859403, 1.37702099489081330271,
516         2.1623699359449663589e-1, 1.34204006088543189037e-2, 3.28014464682127739104e-4,
517         2.89247864745380683936e-6, 6.79019408009981274425e-9
518     ],
519 
520     /**
521      *
522      *	Inverse of Normal distribution function
523      *
524      * SYNOPSIS:
525      *
526      * double x, y, ndtri();
527      *
528      * x = ndtri( y );
529      *
530      * DESCRIPTION:
531      *
532      * Returns the argument, x, for which the area under the
533      * Gaussian probability density function (integrated from
534      * minus infinity to x) is equal to y.
535      *
536      *
537      * For small arguments 0 < y < exp(-2), the program computes
538      * z = sqrt( -2.0 * log(y) );  then the approximation is
539      * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
540      * There are two rational functions P/Q, one for 0 < y < exp(-32)
541      * and the other for y up to exp(-2).  For larger arguments,
542      * w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
543      *
544      *
545      * ACCURACY:
546      *
547      *                      Relative error:
548      * arithmetic   domain        # trials      peak         rms
549      *    DEC      0.125, 1         5500       9.5e-17     2.1e-17
550      *    DEC      6e-39, 0.135     3500       5.7e-17     1.3e-17
551      *    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
552      *    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
553      *
554      *
555      * ERROR MESSAGES:
556      *
557      *   message         condition    value returned
558      * ndtri domain       x <= 0        -MAXNUM
559      * ndtri domain       x >= 1         MAXNUM
560      *
561      * @param  {Number} y0
562      * @returns {Number}
563      */
564     ndtri: function (y0) {
565         var x, y, z, y2, x0, x1, code;
566 
567         if (y0 <= 0.0) {
568             //console.log("ndtri", "DOMAIN ");
569             return -Infinity; // -this.MAXNUM;
570         }
571         if (y0 >= 1.0) {
572             // console.log("ndtri", "DOMAIN");
573             return Infinity; // this.MAXNUM;
574         }
575 
576         code = 1;
577         y = y0;
578         if (y > 1.0 - 0.13533528323661269189) {
579             // 0.135... = exp(-2)
580             y = 1.0 - y;
581             code = 0;
582         }
583 
584         if (y > 0.13533528323661269189) {
585             y = y - 0.5;
586             y2 = y * y;
587             x = y + y * ((y2 * this.polevl(y2, this.P0, 4)) / this.p1evl(y2, this.Q0, 8));
588             x = x * this.s2pi;
589             return x;
590         }
591 
592         x = Math.sqrt(-2.0 * Math.log(y));
593         x0 = x - Math.log(x) / x;
594 
595         z = 1.0 / x;
596         if (x < 8.0) {
597             // y > exp(-32) = 1.2664165549e-14
598             x1 = (z * this.polevl(z, this.P1, 8)) / this.p1evl(z, this.Q1, 8);
599         } else {
600             x1 = (z * this.polevl(z, this.P2, 8)) / this.p1evl(z, this.Q2, 8);
601         }
602         x = x0 - x1;
603         if (code !== 0) {
604             x = -x;
605         }
606         return x;
607     },
608 
609     /**
610      * Inverse of error function erf.
611      *
612      * @param  {Number} x
613      * @returns {Number}
614      */
615     erfi: function (x) {
616         return this.ndtri((x + 1) * 0.5) * this.SQRTH;
617     }
618 };
619 
620 export default Mat.ProbFuncs;
621