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