130 |
/* |
/* |
131 |
* Normal cdf for 0.66291 < fabs(x) < sqrt(32). |
* Normal cdf for 0.66291 < fabs(x) < sqrt(32). |
132 |
*/ |
*/ |
133 |
static double gauss_mediumx_Q ( const double x) |
static double gauss_mediumx ( const double x) |
134 |
{ |
{ |
135 |
int n; |
int n; |
136 |
double temp = 0.0; |
double temp = 0.0; |
172 |
|
|
173 |
result = get_del(x, temp); |
result = get_del(x, temp); |
174 |
|
|
|
if ( x < 0.0 ) |
|
|
{ |
|
|
result = 1.0 - result; |
|
|
} |
|
|
return result; |
|
|
} |
|
|
static double gauss_mediumx_P ( const double x) |
|
|
{ |
|
|
int n; |
|
|
double temp = 0.0; |
|
|
double result = 0.0; |
|
|
double xnum; |
|
|
double xden; |
|
|
double absx; |
|
|
|
|
|
const double c[9] = { |
|
|
0.39894151208813466764, |
|
|
8.8831497943883759412, |
|
|
93.506656132177855979, |
|
|
597.27027639480026226, |
|
|
2494.5375852903726711, |
|
|
6848.1904505362823326, |
|
|
11602.651437647350124, |
|
|
9842.7148383839780218, |
|
|
1.0765576773720192317e-8 |
|
|
}; |
|
|
const double d[8] = { |
|
|
22.266688044328115691, |
|
|
235.38790178262499861, |
|
|
1519.377599407554805, |
|
|
6485.558298266760755, |
|
|
18615.571640885098091, |
|
|
34900.952721145977266, |
|
|
38912.003286093271411, |
|
|
19685.429676859990727 |
|
|
}; |
|
|
absx = fabs(x); |
|
|
|
|
|
xnum = c[8] * absx; |
|
|
xden = absx; |
|
|
for (n = 0; n < 7; ++n) { |
|
|
xnum = (xnum + c[n]) * absx; |
|
|
xden = (xden + d[n]) * absx; |
|
|
} |
|
|
temp = (xnum + c[7]) / (xden + d[7]); |
|
|
|
|
|
result = get_del(x, temp); |
|
|
|
|
|
if( x > 0.0 ) |
|
|
{ |
|
|
result = 1.0 - result; |
|
|
} |
|
|
|
|
175 |
return result; |
return result; |
176 |
} |
} |
177 |
|
|
179 |
* Normal cdf for |
* Normal cdf for |
180 |
* {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }. |
* {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }. |
181 |
*/ |
*/ |
182 |
static double gauss_grandex_Q ( const double x ) |
static double gauss_grandex ( const double x ) |
|
{ |
|
|
int n; |
|
|
double result; |
|
|
double xsq; |
|
|
double temp; |
|
|
double xnum; |
|
|
double xden; |
|
|
double absx; |
|
|
|
|
|
const double p[6] = { |
|
|
0.21589853405795699, |
|
|
0.1274011611602473639, |
|
|
0.022235277870649807, |
|
|
0.001421619193227893466, |
|
|
2.9112874951168792e-5, |
|
|
0.02307344176494017303 |
|
|
}; |
|
|
const double q[5] = { |
|
|
1.28426009614491121, |
|
|
0.468238212480865118, |
|
|
0.0659881378689285515, |
|
|
0.00378239633202758244, |
|
|
7.29751555083966205e-5 |
|
|
}; |
|
|
|
|
|
absx = fabs(x); |
|
|
xsq = 1.0 / (x * x); |
|
|
xnum = p[5] * xsq; |
|
|
xden = xsq; |
|
|
for (n = 0; n < 4; ++n) { |
|
|
xnum = (xnum + p[n]) * xsq; |
|
|
xden = (xden + q[n]) * xsq; |
|
|
} |
|
|
temp = xsq * (xnum + p[4]) / (xden + q[4]); |
|
|
temp = (M_1_SQRT2PI - temp) / absx; |
|
|
|
|
|
result = get_del(x, temp); |
|
|
|
|
|
if( x < 0.0 ) |
|
|
{ |
|
|
result = 1.0 - result; |
|
|
} |
|
|
return result; |
|
|
} |
|
|
static double gauss_grandex_P ( const double x ) |
|
183 |
{ |
{ |
184 |
int n; |
int n; |
185 |
double result; |
double result; |
218 |
|
|
219 |
result = get_del(x, temp); |
result = get_del(x, temp); |
220 |
|
|
|
/* If result is the upper tail area, subtract it from 1.0. */ |
|
|
if( x > 0.0 ) |
|
|
{ |
|
|
result = 1.0 - result; |
|
|
} |
|
|
|
|
221 |
return result; |
return result; |
222 |
} |
} |
223 |
|
|
242 |
} |
} |
243 |
else if ( absx < sqrt(32)) |
else if ( absx < sqrt(32)) |
244 |
{ |
{ |
245 |
result = gauss_mediumx_Q(x); |
result = gauss_mediumx(x); |
246 |
|
if( x < 0.0 ) |
247 |
|
{ |
248 |
|
result = 1.0 - result; |
249 |
|
} |
250 |
} |
} |
251 |
else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER ) |
else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER ) |
252 |
{ |
{ |
253 |
result = gauss_grandex_Q (x); |
result = gauss_grandex (x); |
254 |
|
if ( x < 0.0 ) |
255 |
|
{ |
256 |
|
result = 1.0 - result; |
257 |
|
} |
258 |
} |
} |
259 |
else if ( x > GAUSS_XUPPER ) |
else if ( x > GAUSS_XUPPER ) |
260 |
{ |
{ |
285 |
} |
} |
286 |
else if ( absx < sqrt(32)) |
else if ( absx < sqrt(32)) |
287 |
{ |
{ |
288 |
result = gauss_mediumx_P(x); |
result = gauss_mediumx(x); |
289 |
|
if ( x > 0.0 ) |
290 |
|
{ |
291 |
|
result = 1.0 - result; |
292 |
|
} |
293 |
} |
} |
294 |
else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER ) |
else if ( x < GAUSS_XUPPER || x < GAUSS_XLOWER ) |
295 |
{ |
{ |
296 |
result = gauss_grandex_P (x); |
result = gauss_grandex (x); |
297 |
|
if( x > 0.0 ) |
298 |
|
{ |
299 |
|
result = 1.0 - result; |
300 |
|
} |
301 |
} |
} |
302 |
else if ( x > GAUSS_XUPPER ) |
else if ( x > GAUSS_XUPPER ) |
303 |
{ |
{ |