ndtri.c
Go to the documentation of this file.
1 /* ndtri.c
2  *
3  * Inverse of Normal distribution function
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, ndtri();
10  *
11  * x = ndtri( y );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Returns the argument, x, for which the area under the
18  * Gaussian probability density function (integrated from
19  * minus infinity to x) is equal to y.
20  *
21  *
22  * For small arguments 0 < y < exp(-2), the program computes
23  * z = sqrt( -2.0 * log(y) ); then the approximation is
24  * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
25  * There are two rational functions P/Q, one for 0 < y < exp(-32)
26  * and the other for y up to exp(-2). For larger arguments,
27  * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
28  *
29  *
30  * ACCURACY:
31  *
32  * Relative error:
33  * arithmetic domain # trials peak rms
34  * IEEE 0.125, 1 20000 7.2e-16 1.3e-16
35  * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
36  *
37  *
38  * ERROR MESSAGES:
39  *
40  * message condition value returned
41  * ndtri domain x < 0 NAN
42  * ndtri domain x > 1 NAN
43  *
44  */
45 
46 
47 /*
48  * Cephes Math Library Release 2.1: January, 1989
49  * Copyright 1984, 1987, 1989 by Stephen L. Moshier
50  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
51  */
52 
53 #include "mconf.h"
54 
55 /* sqrt(2pi) */
56 static double s2pi = 2.50662827463100050242E0;
57 
58 /* approximation for 0 <= |y - 0.5| <= 3/8 */
59 static double P0[5] = {
60  -5.99633501014107895267E1,
61  9.80010754185999661536E1,
62  -5.66762857469070293439E1,
63  1.39312609387279679503E1,
64  -1.23916583867381258016E0,
65 };
66 
67 static double Q0[8] = {
68  /* 1.00000000000000000000E0, */
69  1.95448858338141759834E0,
70  4.67627912898881538453E0,
71  8.63602421390890590575E1,
72  -2.25462687854119370527E2,
73  2.00260212380060660359E2,
74  -8.20372256168333339912E1,
75  1.59056225126211695515E1,
76  -1.18331621121330003142E0,
77 };
78 
79 /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
80  * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
81  */
82 static double P1[9] = {
83  4.05544892305962419923E0,
84  3.15251094599893866154E1,
85  5.71628192246421288162E1,
86  4.40805073893200834700E1,
87  1.46849561928858024014E1,
88  2.18663306850790267539E0,
89  -1.40256079171354495875E-1,
90  -3.50424626827848203418E-2,
91  -8.57456785154685413611E-4,
92 };
93 
94 static double Q1[8] = {
95  /* 1.00000000000000000000E0, */
96  1.57799883256466749731E1,
97  4.53907635128879210584E1,
98  4.13172038254672030440E1,
99  1.50425385692907503408E1,
100  2.50464946208309415979E0,
101  -1.42182922854787788574E-1,
102  -3.80806407691578277194E-2,
103  -9.33259480895457427372E-4,
104 };
105 
106 /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
107  * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
108  */
109 
110 static double P2[9] = {
111  3.23774891776946035970E0,
112  6.91522889068984211695E0,
113  3.93881025292474443415E0,
114  1.33303460815807542389E0,
115  2.01485389549179081538E-1,
116  1.23716634817820021358E-2,
117  3.01581553508235416007E-4,
118  2.65806974686737550832E-6,
119  6.23974539184983293730E-9,
120 };
121 
122 static double Q2[8] = {
123  /* 1.00000000000000000000E0, */
124  6.02427039364742014255E0,
125  3.67983563856160859403E0,
126  1.37702099489081330271E0,
127  2.16236993594496635890E-1,
128  1.34204006088543189037E-2,
129  3.28014464682127739104E-4,
130  2.89247864745380683936E-6,
131  6.79019408009981274425E-9,
132 };
133 
134 double ndtri(double y0)
135 {
136  double x, y, z, y2, x0, x1;
137  int code;
138 
139  if (y0 == 0.0) {
140  return -INFINITY;
141  }
142  if (y0 == 1.0) {
143  return INFINITY;
144  }
145  if (y0 < 0.0 || y0 > 1.0) {
146  sf_error("ndtri", SF_ERROR_DOMAIN, NULL);
147  return NAN;
148  }
149  code = 1;
150  y = y0;
151  if (y > (1.0 - 0.13533528323661269189)) { /* 0.135... = exp(-2) */
152  y = 1.0 - y;
153  code = 0;
154  }
155 
156  if (y > 0.13533528323661269189) {
157  y = y - 0.5;
158  y2 = y * y;
159  x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
160  x = x * s2pi;
161  return (x);
162  }
163 
164  x = sqrt(-2.0 * log(y));
165  x0 = x - log(x) / x;
166 
167  z = 1.0 / x;
168  if (x < 8.0) /* y > exp(-32) = 1.2664165549e-14 */
169  x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
170  else
171  x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
172  x = x0 - x1;
173  if (code != 0)
174  x = -x;
175  return (x);
176 }
P1
static double P1[9]
Definition: ndtri.c:82
Q0
static double Q0[8]
Definition: ndtri.c:67
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
s2pi
static double s2pi
Definition: ndtri.c:56
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
y0
double y0(double x)
Definition: j0.c:220
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
x1
Pose3 x1
Definition: testPose3.cpp:663
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
x0
static Symbol x0('x', 0)
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
P2
static double P2[9]
Definition: ndtri.c:110
y
Scalar * y
Definition: level1_cplx_impl.h:124
Q2
static double Q2[8]
Definition: ndtri.c:122
mconf.h
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
Q1
static double Q1[8]
Definition: ndtri.c:94
ndtri
double ndtri(double y0)
Definition: ndtri.c:134
NULL
#define NULL
Definition: ccolamd.c:609
P0
static double P0[5]
Definition: ndtri.c:59
conf.code
code
Definition: wrap/pybind11/docs/conf.py:73
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:33:49