ndtr.c
Go to the documentation of this file.
1 /* ndtr.c
2  *
3  * Normal distribution function
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, ndtr();
10  *
11  * y = ndtr( x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Returns the area under the Gaussian probability density
18  * function, integrated from minus infinity to x:
19  *
20  * x
21  * -
22  * 1 | | 2
23  * ndtr(x) = --------- | exp( - t /2 ) dt
24  * sqrt(2pi) | |
25  * -
26  * -inf.
27  *
28  * = ( 1 + erf(z) ) / 2
29  * = erfc(z) / 2
30  *
31  * where z = x/sqrt(2). Computation is via the functions
32  * erf and erfc.
33  *
34  *
35  * ACCURACY:
36  *
37  * Relative error:
38  * arithmetic domain # trials peak rms
39  * IEEE -13,0 30000 3.4e-14 6.7e-15
40  *
41  *
42  * ERROR MESSAGES:
43  *
44  * message condition value returned
45  * erfc underflow x > 37.519379347 0.0
46  *
47  */
48 /* erf.c
49  *
50  * Error function
51  *
52  *
53  *
54  * SYNOPSIS:
55  *
56  * double x, y, erf();
57  *
58  * y = erf( x );
59  *
60  *
61  *
62  * DESCRIPTION:
63  *
64  * The integral is
65  *
66  * x
67  * -
68  * 2 | | 2
69  * erf(x) = -------- | exp( - t ) dt.
70  * sqrt(pi) | |
71  * -
72  * 0
73  *
74  * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
75  * erf(x) = 1 - erfc(x).
76  *
77  *
78  *
79  * ACCURACY:
80  *
81  * Relative error:
82  * arithmetic domain # trials peak rms
83  * IEEE 0,1 30000 3.7e-16 1.0e-16
84  *
85  */
86 /* erfc.c
87  *
88  * Complementary error function
89  *
90  *
91  *
92  * SYNOPSIS:
93  *
94  * double x, y, erfc();
95  *
96  * y = erfc( x );
97  *
98  *
99  *
100  * DESCRIPTION:
101  *
102  *
103  * 1 - erf(x) =
104  *
105  * inf.
106  * -
107  * 2 | | 2
108  * erfc(x) = -------- | exp( - t ) dt
109  * sqrt(pi) | |
110  * -
111  * x
112  *
113  *
114  * For small x, erfc(x) = 1 - erf(x); otherwise rational
115  * approximations are computed.
116  *
117  *
118  *
119  * ACCURACY:
120  *
121  * Relative error:
122  * arithmetic domain # trials peak rms
123  * IEEE 0,26.6417 30000 5.7e-14 1.5e-14
124  */
125 
126 
127 /*
128  * Cephes Math Library Release 2.2: June, 1992
129  * Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
130  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
131  */
132 
133 #include <float.h> /* DBL_EPSILON */
134 #include "mconf.h"
135 
136 extern double MAXLOG;
137 
138 static double P[] = {
139  2.46196981473530512524E-10,
140  5.64189564831068821977E-1,
141  7.46321056442269912687E0,
142  4.86371970985681366614E1,
143  1.96520832956077098242E2,
144  5.26445194995477358631E2,
145  9.34528527171957607540E2,
146  1.02755188689515710272E3,
147  5.57535335369399327526E2
148 };
149 
150 static double Q[] = {
151  /* 1.00000000000000000000E0, */
152  1.32281951154744992508E1,
153  8.67072140885989742329E1,
154  3.54937778887819891062E2,
155  9.75708501743205489753E2,
156  1.82390916687909736289E3,
157  2.24633760818710981792E3,
158  1.65666309194161350182E3,
159  5.57535340817727675546E2
160 };
161 
162 static double R[] = {
163  5.64189583547755073984E-1,
164  1.27536670759978104416E0,
165  5.01905042251180477414E0,
166  6.16021097993053585195E0,
167  7.40974269950448939160E0,
168  2.97886665372100240670E0
169 };
170 
171 static double S[] = {
172  /* 1.00000000000000000000E0, */
173  2.26052863220117276590E0,
174  9.39603524938001434673E0,
175  1.20489539808096656605E1,
176  1.70814450747565897222E1,
177  9.60896809063285878198E0,
178  3.36907645100081516050E0
179 };
180 
181 static double T[] = {
182  9.60497373987051638749E0,
183  9.00260197203842689217E1,
184  2.23200534594684319226E3,
185  7.00332514112805075473E3,
186  5.55923013010394962768E4
187 };
188 
189 static double U[] = {
190  /* 1.00000000000000000000E0, */
191  3.35617141647503099647E1,
192  5.21357949780152679795E2,
193  4.59432382970980127987E3,
194  2.26290000613890934246E4,
195  4.92673942608635921086E4
196 };
197 
198 #define UTHRESH 37.519379347
199 
200 
201 double ndtr(double a)
202 {
203  double x, y, z;
204 
205  if (cephes_isnan(a)) {
206  sf_error("ndtr", SF_ERROR_DOMAIN, NULL);
207  return NAN;
208  }
209 
210  x = a * M_SQRT1_2;
211  z = fabs(x);
212 
213  if (z < M_SQRT1_2) {
214  y = 0.5 + 0.5 * erf(x);
215  }
216  else {
217  y = 0.5 * erfc(z);
218  if (x > 0) {
219  y = 1.0 - y;
220  }
221  }
222 
223  return y;
224 }
225 
226 
227 double erfc(double a)
228 {
229  double p, q, x, y, z;
230 
231  if (cephes_isnan(a)) {
232  sf_error("erfc", SF_ERROR_DOMAIN, NULL);
233  return NAN;
234  }
235 
236  if (a < 0.0) {
237  x = -a;
238  }
239  else {
240  x = a;
241  }
242 
243  if (x < 1.0) {
244  return 1.0 - erf(a);
245  }
246 
247  z = -a * a;
248 
249  if (z < -MAXLOG) {
250  goto under;
251  }
252 
253  z = exp(z);
254 
255  if (x < 8.0) {
256  p = polevl(x, P, 8);
257  q = p1evl(x, Q, 8);
258  }
259  else {
260  p = polevl(x, R, 5);
261  q = p1evl(x, S, 6);
262  }
263  y = (z * p) / q;
264 
265  if (a < 0) {
266  y = 2.0 - y;
267  }
268 
269  if (y != 0.0) {
270  return y;
271  }
272 
273 under:
274  sf_error("erfc", SF_ERROR_UNDERFLOW, NULL);
275  if (a < 0) {
276  return 2.0;
277  }
278  else {
279  return 0.0;
280  }
281 }
282 
283 
284 
285 double erf(double x)
286 {
287  double y, z;
288 
289  if (cephes_isnan(x)) {
290  sf_error("erf", SF_ERROR_DOMAIN, NULL);
291  return NAN;
292  }
293 
294  if (x < 0.0) {
295  return -erf(-x);
296  }
297 
298  if (fabs(x) > 1.0) {
299  return (1.0 - erfc(x));
300  }
301  z = x * x;
302 
303  y = x * polevl(z, T, 4) / p1evl(z, U, 5);
304  return y;
305 }
cephes_isnan
#define cephes_isnan(x)
Definition: mconf.h:99
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
R
static double R[]
Definition: ndtr.c:162
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
P
static double P[]
Definition: ndtr.c:138
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
SF_ERROR_UNDERFLOW
@ SF_ERROR_UNDERFLOW
Definition: sf_error.h:11
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
Eigen::Triplet
A small structure to hold a non zero as a triplet (i,j,value).
Definition: SparseUtil.h:162
M_SQRT1_2
#define M_SQRT1_2
Definition: mconf.h:124
U
static double U[]
Definition: ndtr.c:189
MAXLOG
double MAXLOG
Definition: const.c:57
ndtr
double ndtr(double a)
Definition: ndtr.c:201
y
Scalar * y
Definition: level1_cplx_impl.h:124
Eigen::Quaternion
The quaternion class used to represent 3D orientations and rotations.
Definition: ForwardDeclarations.h:293
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
p
float * p
Definition: Tutorial_Map_using.cpp:9
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
erf
double erf(double x)
Definition: ndtr.c:285
NULL
#define NULL
Definition: ccolamd.c:609
erfc
double erfc(double a)
Definition: ndtr.c:227
S
static double S[]
Definition: ndtr.c:171
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:03:11