Logo ROOT  
Reference Guide
SpecFuncCephesInv.cxx
Go to the documentation of this file.
1// inverse of gamma and beta from Cephes library
2// see: http://www.netlib.org/cephes
3//
4// Copyright 1985, 1987, 2000 by Stephen L. Moshier
5
6
7
8#include "Math/Error.h"
9
10#include "SpecFuncCephes.h"
11
12#include <cmath>
13
14#include <limits>
15
16namespace ROOT {
17namespace Math {
18
19namespace Cephes {
20
21
22
23/*
24 *
25 * Inverse of Normal distribution function
26 *
27 *
28 *
29 * SYNOPSIS:
30 *
31 * double x, y, ndtri();
32 *
33 * x = ndtri( y );
34 *
35 *
36 *
37 * DESCRIPTION:
38 *
39 * Returns the argument, x, for which the area under the
40 * Gaussian probability density function (integrated from
41 * minus infinity to x) is equal to y.
42 *
43 *
44 * For small arguments 0 < y < exp(-2), the program computes
45 * z = sqrt( -2.0 * log(y) ); then the approximation is
46 * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
47 * There are two rational functions P/Q, one for 0 < y < exp(-32)
48 * and the other for y up to exp(-2). For larger arguments,
49 * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
50 *
51 *
52 * ACCURACY:
53 *
54 * Relative error:
55 * arithmetic domain # trials peak rms
56 * DEC 0.125, 1 5500 9.5e-17 2.1e-17
57 * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17
58 * IEEE 0.125, 1 20000 7.2e-16 1.3e-16
59 * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
60 *
61 *
62 * ERROR MESSAGES:
63 *
64 * message condition value returned
65 * ndtri domain x <= 0 -MAXNUM
66 * ndtri domain x >= 1 MAXNUM
67 *
68 */
69
70/*
71Cephes Math Library Release 2.8: June, 2000
72Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
73*/
74
75
76static double s2pi = 2.50662827463100050242E0;
77
78static double P0[5] = {
79-5.99633501014107895267E1,
80 9.80010754185999661536E1,
81-5.66762857469070293439E1,
82 1.39312609387279679503E1,
83-1.23916583867381258016E0,
84};
85static double Q0[8] = {
86 1.95448858338141759834E0,
87 4.67627912898881538453E0,
88 8.63602421390890590575E1,
89-2.25462687854119370527E2,
90 2.00260212380060660359E2,
91-8.20372256168333339912E1,
92 1.59056225126211695515E1,
93-1.18331621121330003142E0,
94};
95static double P1[9] = {
96 4.05544892305962419923E0,
97 3.15251094599893866154E1,
98 5.71628192246421288162E1,
99 4.40805073893200834700E1,
100 1.46849561928858024014E1,
101 2.18663306850790267539E0,
102-1.40256079171354495875E-1,
103-3.50424626827848203418E-2,
104-8.57456785154685413611E-4,
105};
106static double Q1[8] = {
107 1.57799883256466749731E1,
108 4.53907635128879210584E1,
109 4.13172038254672030440E1,
110 1.50425385692907503408E1,
111 2.50464946208309415979E0,
112-1.42182922854787788574E-1,
113-3.80806407691578277194E-2,
114-9.33259480895457427372E-4,
115};
116static double P2[9] = {
117 3.23774891776946035970E0,
118 6.91522889068984211695E0,
119 3.93881025292474443415E0,
120 1.33303460815807542389E0,
121 2.01485389549179081538E-1,
122 1.23716634817820021358E-2,
123 3.01581553508235416007E-4,
124 2.65806974686737550832E-6,
125 6.23974539184983293730E-9,
126};
127static double Q2[8] = {
128 6.02427039364742014255E0,
129 3.67983563856160859403E0,
130 1.37702099489081330271E0,
131 2.16236993594496635890E-1,
132 1.34204006088543189037E-2,
133 3.28014464682127739104E-4,
134 2.89247864745380683936E-6,
135 6.79019408009981274425E-9,
136};
137double ndtri( double y0 )
138{
139 double x, y, z, y2, x0, x1;
140 int code;
141 if( y0 <= 0.0 )
142 return( - std::numeric_limits<double>::infinity() );
143 if( y0 >= 1.0 )
144 return( + std::numeric_limits<double>::infinity() );
145 code = 1;
146 y = y0;
147 if( y > (1.0 - 0.13533528323661269189) )
148 {
149 y = 1.0 - y;
150 code = 0;
151 }
152 if( y > 0.13533528323661269189 )
153 {
154 y = y - 0.5;
155 y2 = y * y;
156 x = y + y * (y2 * Polynomialeval( y2, P0, 4)/ Polynomial1eval( y2, Q0, 8 ));
157 x = x * s2pi;
158 return(x);
159 }
160 x = std::sqrt( -2.0 * std::log(y) );
161 x0 = x - std::log(x)/x;
162 z = 1.0/x;
163 if( x < 8.0 )
164 x1 = z * Polynomialeval( z, P1, 8 )/ Polynomial1eval ( z, Q1, 8 );
165 else
166 x1 = z * Polynomialeval( z, P2, 8 )/ Polynomial1eval( z, Q2, 8 );
167 x = x0 - x1;
168 if( code != 0 )
169 x = -x;
170 return( x );
171}
172
173
174
175
176/*
177 *
178 * Inverse of complemented imcomplete gamma integral
179 *
180 *
181 *
182 * SYNOPSIS:
183 *
184 * double a, x, p, igami();
185 *
186 * x = igami( a, p );
187 *
188 * DESCRIPTION:
189 *
190 * Given p, the function finds x such that
191 *
192 * igamc( a, x ) = p.
193 *
194 * Starting with the approximate value
195 *
196 * 3
197 * x = a t
198 *
199 * where
200 *
201 * t = 1 - d - ndtri(p) sqrt(d)
202 *
203 * and
204 *
205 * d = 1/9a,
206 *
207 * the routine performs up to 10 Newton iterations to find the
208 * root of igamc(a,x) - p = 0.
209 *
210 * ACCURACY:
211 *
212 * Tested at random a, p in the intervals indicated.
213 *
214 * a p Relative error:
215 * arithmetic domain domain # trials peak rms
216 * IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
217 * IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
218 * IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
219 */
220/*
221Cephes Math Library Release 2.8: June, 2000
222Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
223*/
224
225double igami( double a, double y0 )
226{
227 double x0, x1, x, yl, yh, y, d, lgm, dithresh;
228 int i, dir;
229
230 // check the domain
231 if (a<= 0) {
232 MATH_ERROR_MSG("Cephes::igami","Wrong domain for parameter a (must be > 0)");
233 return 0;
234 }
235 if (y0 <= 0) {
236 //if (y0<0) MATH_ERROR_MSG("Cephes::igami","Wrong domain for y (must be in [0,1])");
237 return std::numeric_limits<double>::infinity();
238 }
239 if (y0 >= 1) {
240 //if (y0>1) MATH_ERROR_MSG("Cephes::igami","Wrong domain for y (must be in [0,1])");
241 return 0;
242 }
243
244
245/* bound the solution */
246 static double kMAXNUM = std::numeric_limits<double>::max();
247 x0 = kMAXNUM;
248 yl = 0;
249 x1 = 0;
250 yh = 1.0;
251 dithresh = 5.0 * kMACHEP;
252
253/* approximation to inverse function */
254 d = 1.0/(9.0*a);
255 y = ( 1.0 - d - ndtri(y0) * std::sqrt(d) );
256 x = a * y * y * y;
257
258 lgm = lgam(a);
259
260 for( i=0; i<10; i++ )
261 {
262 if( x > x0 || x < x1 )
263 goto ihalve;
264 y = igamc(a,x);
265 if( y < yl || y > yh )
266 goto ihalve;
267 if( y < y0 )
268 {
269 x0 = x;
270 yl = y;
271 }
272 else
273 {
274 x1 = x;
275 yh = y;
276 }
277/* compute the derivative of the function at this point */
278 d = (a - 1.0) * std::log(x) - x - lgm;
279 if( d < -kMAXLOG )
280 goto ihalve;
281 d = -std::exp(d);
282/* compute the step to the next approximation of x */
283 d = (y - y0)/d;
284 if( std::fabs(d/x) < kMACHEP )
285 goto done;
286 x = x - d;
287 }
288
289/* Resort to interval halving if Newton iteration did not converge. */
290ihalve:
291
292 d = 0.0625;
293 if( x0 == kMAXNUM )
294 {
295 if( x <= 0.0 )
296 x = 1.0;
297 while( x0 == kMAXNUM )
298 {
299 x = (1.0 + d) * x;
300 y = igamc( a, x );
301 if( y < y0 )
302 {
303 x0 = x;
304 yl = y;
305 break;
306 }
307 d = d + d;
308 }
309 }
310 d = 0.5;
311 dir = 0;
312
313 for( i=0; i<400; i++ )
314 {
315 x = x1 + d * (x0 - x1);
316 y = igamc( a, x );
317 lgm = (x0 - x1)/(x1 + x0);
318 if( std::fabs(lgm) < dithresh )
319 break;
320 lgm = (y - y0)/y0;
321 if( std::fabs(lgm) < dithresh )
322 break;
323 if( x <= 0.0 )
324 break;
325 if( y >= y0 )
326 {
327 x1 = x;
328 yh = y;
329 if( dir < 0 )
330 {
331 dir = 0;
332 d = 0.5;
333 }
334 else if( dir > 1 )
335 d = 0.5 * d + 0.5;
336 else
337 d = (y0 - yl)/(yh - yl);
338 dir += 1;
339 }
340 else
341 {
342 x0 = x;
343 yl = y;
344 if( dir > 0 )
345 {
346 dir = 0;
347 d = 0.5;
348 }
349 else if( dir < -1 )
350 d = 0.5 * d;
351 else
352 d = (y0 - yl)/(yh - yl);
353 dir -= 1;
354 }
355 }
356
357// if( x == 0.0 )
358// mtherr( "igami", UNDERFLOW );
359
360done:
361 return( x );
362}
363
364
365/*
366 *
367 * Inverse of imcomplete beta integral
368 *
369 *
370 *
371 * SYNOPSIS:
372 *
373 * double a, b, x, y, incbi();
374 *
375 * x = incbi( a, b, y );
376 *
377 *
378 *
379 * DESCRIPTION:
380 *
381 * Given y, the function finds x such that
382 *
383 * incbet( a, b, x ) = y .
384 *
385 * The routine performs interval halving or Newton iterations to find the
386 * root of incbet(a,b,x) - y = 0.
387 *
388 *
389 * ACCURACY:
390 *
391 * Relative error:
392 * x a,b
393 * arithmetic domain domain # trials peak rms
394 * IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
395 * IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
396 * IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
397 * VAX 0,1 .5,100 25000 3.5e-14 1.1e-15
398 * With a and b constrained to half-integer or integer values:
399 * IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
400 * IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
401 * With a = .5, b constrained to half-integer or integer values:
402 * IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
403 */
404
405/*
406Cephes Math Library Release 2.8: June, 2000
407Copyright 1984, 1996, 2000 by Stephen L. Moshier
408*/
409
410
411double incbi( double aa, double bb, double yy0 )
412{
413 double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt;
414 int i, rflg, dir, nflg;
415
416 // check the domain
417 if (aa<= 0) {
418 MATH_ERROR_MSG("Cephes::incbi","Wrong domain for parameter a (must be > 0)");
419 return 0;
420 }
421 if (bb<= 0) {
422 MATH_ERROR_MSG("Cephes::incbi","Wrong domain for parameter b (must be > 0)");
423 return 0;
424 }
425
426
427 i = 0;
428 if( yy0 <= 0 )
429 return(0.0);
430 if( yy0 >= 1.0 )
431 return(1.0);
432 x0 = 0.0;
433 yl = 0.0;
434 x1 = 1.0;
435 yh = 1.0;
436 nflg = 0;
437
438 if( aa <= 1.0 || bb <= 1.0 )
439 {
440 dithresh = 1.0e-6;
441 rflg = 0;
442 a = aa;
443 b = bb;
444 y0 = yy0;
445 x = a/(a+b);
446 y = incbet( a, b, x );
447 goto ihalve;
448 }
449 else
450 {
451 dithresh = 1.0e-4;
452 }
453/* approximation to inverse function */
454
455 yp = -ndtri(yy0);
456
457 if( yy0 > 0.5 )
458 {
459 rflg = 1;
460 a = bb;
461 b = aa;
462 y0 = 1.0 - yy0;
463 yp = -yp;
464 }
465 else
466 {
467 rflg = 0;
468 a = aa;
469 b = bb;
470 y0 = yy0;
471 }
472
473 lgm = (yp * yp - 3.0)/6.0;
474 x = 2.0/( 1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0) );
475 d = yp * std::sqrt( x + lgm ) / x
476 - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
477 * (lgm + 5.0/6.0 - 2.0/(3.0*x));
478 d = 2.0 * d;
479 if( d < kMINLOG )
480 {
481 x = 1.0;
482 goto under;
483 }
484 x = a/( a + b * std::exp(d) );
485 y = incbet( a, b, x );
486 yp = (y - y0)/y0;
487 if( std::fabs(yp) < 0.2 )
488 goto newt;
489
490/* Resort to interval halving if not close enough. */
491ihalve:
492
493 dir = 0;
494 di = 0.5;
495 for( i=0; i<100; i++ )
496 {
497 if( i != 0 )
498 {
499 x = x0 + di * (x1 - x0);
500 if( x == 1.0 )
501 x = 1.0 - kMACHEP;
502 if( x == 0.0 )
503 {
504 di = 0.5;
505 x = x0 + di * (x1 - x0);
506 if( x == 0.0 )
507 goto under;
508 }
509 y = incbet( a, b, x );
510 yp = (x1 - x0)/(x1 + x0);
511 if( std::fabs(yp) < dithresh )
512 goto newt;
513 yp = (y-y0)/y0;
514 if( std::fabs(yp) < dithresh )
515 goto newt;
516 }
517 if( y < y0 )
518 {
519 x0 = x;
520 yl = y;
521 if( dir < 0 )
522 {
523 dir = 0;
524 di = 0.5;
525 }
526 else if( dir > 3 )
527 di = 1.0 - (1.0 - di) * (1.0 - di);
528 else if( dir > 1 )
529 di = 0.5 * di + 0.5;
530 else
531 di = (y0 - y)/(yh - yl);
532 dir += 1;
533 if( x0 > 0.75 )
534 {
535 if( rflg == 1 )
536 {
537 rflg = 0;
538 a = aa;
539 b = bb;
540 y0 = yy0;
541 }
542 else
543 {
544 rflg = 1;
545 a = bb;
546 b = aa;
547 y0 = 1.0 - yy0;
548 }
549 x = 1.0 - x;
550 y = incbet( a, b, x );
551 x0 = 0.0;
552 yl = 0.0;
553 x1 = 1.0;
554 yh = 1.0;
555 goto ihalve;
556 }
557 }
558 else
559 {
560 x1 = x;
561 if( rflg == 1 && x1 < kMACHEP )
562 {
563 x = 0.0;
564 goto done;
565 }
566 yh = y;
567 if( dir > 0 )
568 {
569 dir = 0;
570 di = 0.5;
571 }
572 else if( dir < -3 )
573 di = di * di;
574 else if( dir < -1 )
575 di = 0.5 * di;
576 else
577 di = (y - y0)/(yh - yl);
578 dir -= 1;
579 }
580 }
581 //mtherr( "incbi", PLOSS );
582 if( x0 >= 1.0 )
583 {
584 x = 1.0 - kMACHEP;
585 goto done;
586 }
587 if( x <= 0.0 )
588 {
589 under:
590 //mtherr( "incbi", UNDERFLOW );
591 x = 0.0;
592 goto done;
593 }
594
595newt:
596
597 if( nflg )
598 goto done;
599 nflg = 1;
600 lgm = lgam(a+b) - lgam(a) - lgam(b);
601
602 for( i=0; i<8; i++ )
603 {
604 /* Compute the function at this point. */
605 if( i != 0 )
606 y = incbet(a,b,x);
607 if( y < yl )
608 {
609 x = x0;
610 y = yl;
611 }
612 else if( y > yh )
613 {
614 x = x1;
615 y = yh;
616 }
617 else if( y < y0 )
618 {
619 x0 = x;
620 yl = y;
621 }
622 else
623 {
624 x1 = x;
625 yh = y;
626 }
627 if( x == 1.0 || x == 0.0 )
628 break;
629 /* Compute the derivative of the function at this point. */
630 d = (a - 1.0) * std::log(x) + (b - 1.0) * std::log(1.0-x) + lgm;
631 if( d < kMINLOG )
632 goto done;
633 if( d > kMAXLOG )
634 break;
635 d = std::exp(d);
636 /* Compute the step to the next approximation of x. */
637 d = (y - y0)/d;
638 xt = x - d;
639 if( xt <= x0 )
640 {
641 y = (x - x0) / (x1 - x0);
642 xt = x0 + 0.5 * y * (x - x0);
643 if( xt <= 0.0 )
644 break;
645 }
646 if( xt >= x1 )
647 {
648 y = (x1 - x) / (x1 - x0);
649 xt = x1 - 0.5 * y * (x1 - x);
650 if( xt >= 1.0 )
651 break;
652 }
653 x = xt;
654 if( std::fabs(d/x) < 128.0 * kMACHEP )
655 goto done;
656 }
657/* Did not converge. */
658 dithresh = 256.0 * kMACHEP;
659 goto ihalve;
660
661done:
662
663 if( rflg )
664 {
665 if( x <= kMACHEP )
666 x = 1.0 - kMACHEP;
667 else
668 x = 1.0 - x;
669 }
670 return( x );
671}
672
673} // end namespace Cephes
674
675} // end namespace Math
676} // end namespace ROOT
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:83
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
static const double x1[5]
#define kMACHEP
#define kMAXLOG
#define kMINLOG
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
static double Q0[8]
static double P0[5]
double incbet(double aa, double bb, double xx)
DESCRIPTION:
static double P1[9]
double lgam(double x)
double igamc(double a, double x)
incomplete complementary gamma function igamc(a, x) = 1 - igam(a, x)
double igami(double a, double y)
double incbi(double a, double b, double y)
double ndtri(double y)
double Polynomial1eval(double x, double *a, unsigned int N)
double Polynomialeval(double x, double *a, unsigned int N)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
auto * a
Definition: textangle.C:12