Branch data Line data Source code
1 : : /*
2 : : * fspecial.cpp - special functions implementation
3 : : *
4 : : * Copyright (C) 2006, 2008 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : : #include <stdlib.h>
31 : : #include <string.h>
32 : : #include <cmath>
33 : : #include <float.h>
34 : :
35 : : #include "compat.h"
36 : : #include "constants.h"
37 : : #include "fspecial.h"
38 : :
39 : : #include <limits>
40 : :
41 : : /* The function computes the complete elliptic integral of first kind
42 : : K() and the second kind E() using the arithmetic-geometric mean
43 : : algorithm (AGM) by Abramowitz and Stegun. */
44 : 0 : void fspecial::ellip_ke (nr_double_t arg, nr_double_t &k, nr_double_t &e) {
45 : 0 : int iMax = 16;
46 [ # # ]: 0 : if (arg == 1.0) {
47 : 0 : k = std::numeric_limits<nr_double_t>::infinity();
48 : 0 : e = 0;
49 : : }
50 [ # # ][ # # ]: 0 : else if (std::isinf (arg) && arg < 0) {
[ # # ]
51 : 0 : k = 0;
52 : 0 : e = std::numeric_limits<nr_double_t>::infinity();
53 : : }
54 : : else {
55 : 0 : nr_double_t a, b, c, f, s, fk = 1, fe = 1, t, da = arg;
56 : : int i;
57 [ # # ]: 0 : if (arg < 0) {
58 : 0 : fk = 1 / sqrt (1 - arg);
59 : 0 : fe = sqrt (1 - arg);
60 : 0 : da = -arg / (1 - arg);
61 : : }
62 : 0 : a = 1;
63 : 0 : b = sqrt (1 - da);
64 : 0 : c = sqrt (da);
65 : 0 : f = 0.5;
66 : 0 : s = f * c * c;
67 [ # # ]: 0 : for (i = 0; i < iMax; i++) {
68 : 0 : t = (a + b) / 2;
69 : 0 : c = (a - b) / 2;
70 : 0 : b = sqrt (a * b);
71 : 0 : a = t;
72 : 0 : f *= 2;
73 : 0 : s += f * c * c;
74 [ # # ]: 0 : if (c / a < std::numeric_limits<nr_double_t>::epsilon()) break;
75 : : }
76 [ # # ]: 0 : if (i >= iMax) {
77 : 0 : k = 0; e = 0;
78 : : }
79 : : else {
80 : 0 : k = M_PI_2 / a;
81 : 0 : e = M_PI_2 * (1 - s) / a;
82 [ # # ]: 0 : if (arg < 0) {
83 : 0 : k *= fk;
84 : 0 : e *= fe;
85 : : }
86 : : }
87 : : }
88 : 0 : }
89 : :
90 : : const nr_double_t SN_ACC = 1e-5; // Accuracy of sn(x) is SN_ACC^2
91 : : const nr_double_t K_ERR = 1e-8; // Accuracy of K(k)
92 : :
93 : : // Computes Carlson's elliptic integral of the first kind.
94 : 0 : nr_double_t fspecial::ellip_rf (nr_double_t x, nr_double_t y, nr_double_t z) {
95 : : nr_double_t al, av, dx, dy, dz, e2, e3;
96 : : nr_double_t sx, sy, sz, xt, yt, zt;
97 : :
98 : : // constants
99 : 0 : const nr_double_t c1 = 1.0 / 24.0;
100 : 0 : const nr_double_t c2 = 0.1;
101 : 0 : const nr_double_t c3 = 3.0 / 44.0;
102 : 0 : const nr_double_t c4 = 1.0 / 14.0;
103 : :
104 : 0 : xt = x;
105 : 0 : yt = y;
106 : 0 : zt = z;
107 [ # # ][ # # ]: 0 : do {
[ # # ][ # # ]
108 : 0 : sx = sqrt (xt);
109 : 0 : sy = sqrt (yt);
110 : 0 : sz = sqrt (zt);
111 : 0 : al = sx * (sy + sz) + sy * sz;
112 : 0 : xt = 0.25 * (xt + al);
113 : 0 : yt = 0.25 * (yt + al);
114 : 0 : zt = 0.25 * (zt + al);
115 : 0 : av = (xt + yt + zt) / 3.0;
116 : 0 : dx = (av - xt) / av;
117 : 0 : dy = (av - yt) / av;
118 : 0 : dz = (av - zt) / av;
119 : : }
120 : 0 : while (MAX (MAX (fabs (dx), fabs (dy)), fabs (dz)) > K_ERR);
121 : :
122 : 0 : e2 = dx * dy - dz * dz;
123 : 0 : e3 = dx * dy * dz;
124 : 0 : return (1 + (c1 * e2 - c2 - c3 * e3) * e2 + c4 * e3) / sqrt (av);
125 : : }
126 : :
127 : : // Compute the Jacobian elliptic functions sn (u,k), cn (u,k) and dn (u,k).
128 : 0 : nr_double_t fspecial::ellip_sncndn (nr_double_t u, nr_double_t k,
129 : : nr_double_t& sn, nr_double_t& cn, nr_double_t& dn) {
130 : : nr_double_t a, b, c, d;
131 : : nr_double_t fn[14], en[14];
132 : : int l;
133 : : bool bo;
134 : :
135 : 0 : d = 1 - k;
136 [ # # ]: 0 : if (k) {
137 : 0 : bo = (k < 0);
138 [ # # ]: 0 : if (bo) {
139 : 0 : k /= -1 / d;
140 : 0 : u *= (d = sqrt (d));
141 : : }
142 : 0 : a = 1;
143 : 0 : dn = 1;
144 [ # # ]: 0 : for (int i = 1; i < 14; i++) {
145 : 0 : l = i;
146 : 0 : fn[i] = a;
147 : 0 : en[i] = (k = sqrt (k));
148 : 0 : c = 0.5 * (a + k);
149 [ # # ]: 0 : if (fabs (a - k) <= SN_ACC * a)
150 : 0 : break;
151 : 0 : k *= a;
152 : 0 : a = c;
153 : : }
154 : 0 : u *= c;
155 : 0 : sn = sin (u);
156 : 0 : cn = cos (u);
157 [ # # ]: 0 : if (sn) {
158 : 0 : a = cn / sn;
159 : 0 : c *= a;
160 [ # # ]: 0 : for (int ii = l; ii > 0; ii--) {
161 : 0 : b = fn[ii];
162 : 0 : a *= c;
163 : 0 : c *= dn;
164 : 0 : dn = (en[ii] + a) / (b + a);
165 : 0 : a = c / b;
166 : : }
167 : 0 : a = 1 / sqrt (c * c + 1);
168 [ # # ]: 0 : sn = (sn >= 0 ? a : -a);
169 : 0 : cn = sn * c;
170 : : }
171 [ # # ]: 0 : if (bo) {
172 : 0 : a = dn;
173 : 0 : dn = cn;
174 : 0 : cn = a;
175 : 0 : sn /= d;
176 : : }
177 : : }
178 : : else {
179 : 0 : cn = 1 / cosh (u);
180 : 0 : dn = cn;
181 : 0 : sn = tanh (u);
182 : : }
183 : 0 : return sn;
184 : : }
185 : :
186 : : /* data for a Chebyshev series over a given interval */
187 : : struct cheb_series_t {
188 : : nr_double_t * c; /* coefficients */
189 : : int order; /* order of expansion */
190 : : nr_double_t a; /* lower interval point */
191 : : nr_double_t b; /* upper interval point */
192 : : };
193 : : typedef struct cheb_series_t cheb_series;
194 : :
195 : 0 : static nr_double_t cheb_eval (const cheb_series * cs, const nr_double_t x) {
196 : 0 : nr_double_t d = 0.0;
197 : 0 : nr_double_t dd = 0.0;
198 : 0 : nr_double_t y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a);
199 : 0 : nr_double_t y2 = 2.0 * y;
200 [ # # ]: 0 : for (int i = cs->order; i >= 1; i--) {
201 : 0 : nr_double_t t = d;
202 : 0 : d = y2 * d - dd + cs->c[i];
203 : 0 : dd = t;
204 : : }
205 : 0 : d = y * d - dd + 0.5 * cs->c[0];
206 : 0 : return d;
207 : : }
208 : :
209 : : #if !defined (HAVE_ERF) || !defined (HAVE_ERFC)
210 : :
211 : : /* Chebyshev fit for erfc ((t+1)/2), -1 < t < 1 */
212 : : static nr_double_t erfc_xlt1_data[20] = {
213 : : 1.06073416421769980345174155056,
214 : : -0.42582445804381043569204735291,
215 : : 0.04955262679620434040357683080,
216 : : 0.00449293488768382749558001242,
217 : : -0.00129194104658496953494224761,
218 : : -0.00001836389292149396270416979,
219 : : 0.00002211114704099526291538556,
220 : : -5.23337485234257134673693179020e-7,
221 : : -2.78184788833537885382530989578e-7,
222 : : 1.41158092748813114560316684249e-8,
223 : : 2.72571296330561699984539141865e-9,
224 : : -2.06343904872070629406401492476e-10,
225 : : -2.14273991996785367924201401812e-11,
226 : : 2.22990255539358204580285098119e-12,
227 : : 1.36250074650698280575807934155e-13,
228 : : -1.95144010922293091898995913038e-14,
229 : : -6.85627169231704599442806370690e-16,
230 : : 1.44506492869699938239521607493e-16,
231 : : 2.45935306460536488037576200030e-18,
232 : : -9.29599561220523396007359328540e-19
233 : : };
234 : : static cheb_series erfc_xlt1_cs = {
235 : : erfc_xlt1_data, 19, -1, 1
236 : : };
237 : :
238 : : /* Chebyshev fit for erfc (x) * exp (x^2), 1 < x < 5, x = 2t + 3, -1 < t < 1 */
239 : : static nr_double_t erfc_x15_data[25] = {
240 : : 0.44045832024338111077637466616,
241 : : -0.143958836762168335790826895326,
242 : : 0.044786499817939267247056666937,
243 : : -0.013343124200271211203618353102,
244 : : 0.003824682739750469767692372556,
245 : : -0.001058699227195126547306482530,
246 : : 0.000283859419210073742736310108,
247 : : -0.000073906170662206760483959432,
248 : : 0.000018725312521489179015872934,
249 : : -4.62530981164919445131297264430e-6,
250 : : 1.11558657244432857487884006422e-6,
251 : : -2.63098662650834130067808832725e-7,
252 : : 6.07462122724551777372119408710e-8,
253 : : -1.37460865539865444777251011793e-8,
254 : : 3.05157051905475145520096717210e-9,
255 : : -6.65174789720310713757307724790e-10,
256 : : 1.42483346273207784489792999706e-10,
257 : : -3.00141127395323902092018744545e-11,
258 : : 6.22171792645348091472914001250e-12,
259 : : -1.26994639225668496876152836555e-12,
260 : : 2.55385883033257575402681845385e-13,
261 : : -5.06258237507038698392265499770e-14,
262 : : 9.89705409478327321641264227110e-15,
263 : : -1.90685978789192181051961024995e-15,
264 : : 3.50826648032737849245113757340e-16
265 : : };
266 : : static cheb_series erfc_x15_cs = {
267 : : erfc_x15_data, 24, -1, 1
268 : : };
269 : :
270 : : /* Chebyshev fit for erfc(x) * exp(x^2),
271 : : 5 < x < 10, x = (5t + 15)/2, -1 < t < */
272 : : static nr_double_t erfc_x510_data[20] = {
273 : : 1.11684990123545698684297865808,
274 : : 0.003736240359381998520654927536,
275 : : -0.000916623948045470238763619870,
276 : : 0.000199094325044940833965078819,
277 : : -0.000040276384918650072591781859,
278 : : 7.76515264697061049477127605790e-6,
279 : : -1.44464794206689070402099225301e-6,
280 : : 2.61311930343463958393485241947e-7,
281 : : -4.61833026634844152345304095560e-8,
282 : : 8.00253111512943601598732144340e-9,
283 : : -1.36291114862793031395712122089e-9,
284 : : 2.28570483090160869607683087722e-10,
285 : : -3.78022521563251805044056974560e-11,
286 : : 6.17253683874528285729910462130e-12,
287 : : -9.96019290955316888445830597430e-13,
288 : : 1.58953143706980770269506726000e-13,
289 : : -2.51045971047162509999527428316e-14,
290 : : 3.92607828989125810013581287560e-15,
291 : : -6.07970619384160374392535453420e-16,
292 : : 9.12600607264794717315507477670e-17
293 : : };
294 : : static cheb_series erfc_x510_cs = {
295 : : erfc_x510_data, 19, -1, 1
296 : : };
297 : :
298 : : /* Estimates erfc (x) valid for 8 < x < 100, this is based on index 5725
299 : : in Hart et al. */
300 : : static nr_double_t erfc8 (nr_double_t x) {
301 : : static nr_double_t p[] = {
302 : : 2.97886562639399288862,
303 : : 7.409740605964741794425,
304 : : 6.1602098531096305440906,
305 : : 5.019049726784267463450058,
306 : : 1.275366644729965952479585264,
307 : : 0.5641895835477550741253201704
308 : : };
309 : : static nr_double_t q[] = {
310 : : 3.3690752069827527677,
311 : : 9.608965327192787870698,
312 : : 17.08144074746600431571095,
313 : : 12.0489519278551290360340491,
314 : : 9.396034016235054150430579648,
315 : : 2.260528520767326969591866945,
316 : : 1.0
317 : : };
318 : : nr_double_t n = 0.0, d = 0.0;
319 : : int i;
320 : : n = p[5];
321 : : for (i = 4; i >= 0; --i) n = x * n + p[i];
322 : : d = q[6];
323 : : for (i = 5; i >= 0; --i) d = x * d + q[i];
324 : :
325 : : return exp (-x * x) * (n / d);
326 : : }
327 : :
328 : : #endif /* !HAVE_ERF || !HAVE_ERFC */
329 : :
330 : : #ifndef HAVE_ERFC
331 : :
332 : : nr_double_t fspecial::erfc (nr_double_t x) {
333 : : const nr_double_t ax = fabs (x);
334 : : nr_double_t val;
335 : :
336 : : if (ax <= 1.0) {
337 : : nr_double_t t = 2.0 * ax - 1.0;
338 : : val = cheb_eval (&erfc_xlt1_cs, t);
339 : : }
340 : : else if (ax <= 5.0) {
341 : : nr_double_t ex2 = exp (-x * x);
342 : : nr_double_t t = 0.5 * (ax - 3.0);
343 : : val = ex2 * cheb_eval (&erfc_x15_cs, t);
344 : : }
345 : : else if (ax < 10.0) {
346 : : nr_double_t ex = exp(-x * x) / ax;
347 : : nr_double_t t = (2.0 * ax - 15.0) / 5.0;
348 : : val = ex * cheb_eval (&erfc_x510_cs, t);
349 : : }
350 : : else {
351 : : val = erfc8 (ax);
352 : : }
353 : : return (x < 0.0) ? 2.0 - val : val;
354 : : }
355 : : #else
356 : 0 : nr_double_t fspecial::erfc (nr_double_t d) {
357 : 0 : return ::erfc (d);
358 : : }
359 : : #endif /* HAVE_ERFC */
360 : :
361 : : #ifndef HAVE_ERF
362 : :
363 : : /* Abramowitz + Stegun, 7.1.5 */
364 : : static nr_double_t erfseries (nr_double_t x) {
365 : : nr_double_t c = x;
366 : : nr_double_t e = c;
367 : : nr_double_t d;
368 : : for (int k = 1; k < 30; ++k) {
369 : : c *= -x * x / k;
370 : : d = c / (2.0 * k + 1.0);
371 : : e += d;
372 : : }
373 : : return 2.0 / M_SQRTPI * e;
374 : : }
375 : :
376 : : nr_double_t fspecial::erf (nr_double_t x) {
377 : : if (fabs (x) < 1.0) {
378 : : return erfseries (x);
379 : : }
380 : : return 1.0 - erfc (x);
381 : : }
382 : :
383 : : #else
384 : 0 : nr_double_t fspecial::erf (nr_double_t d) {
385 : 0 : return ::erf (d);
386 : : }
387 : : #endif /* HAVE_ERF */
388 : :
389 : : // Inverse of the error function erf().
390 : 0 : nr_double_t fspecial::erfinv (nr_double_t y) {
391 : 0 : nr_double_t x = 0.0; // The output
392 : 0 : nr_double_t z = 0.0; // Intermediate variable
393 : 0 : nr_double_t y0 = 0.7; // Central range variable
394 : :
395 : : // Coefficients in rational approximations.
396 : 0 : nr_double_t a[4] = { 0.886226899, -1.645349621, 0.914624893, -0.140543331};
397 : 0 : nr_double_t b[4] = {-2.118377725, 1.442710462, -0.329097515, 0.012229801};
398 : 0 : nr_double_t c[4] = {-1.970840454, -1.624906493, 3.429567803, 1.641345311};
399 : 0 : nr_double_t d[2] = { 3.543889200, 1.637067800};
400 : :
401 [ # # ][ # # ]: 0 : if (y < -1.0 || 1.0 < y) {
402 : 0 : x = log (-1.0);
403 : : }
404 [ # # ][ # # ]: 0 : else if (y == -1.0 || 1.0 == y) {
405 : 0 : x = -y * log(0.0);
406 : : }
407 [ # # ][ # # ]: 0 : else if (-1.0 < y && y < -y0) {
408 : 0 : z = sqrt(-log((1.0+y)/2.0));
409 : 0 : x = -(((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0);
410 : : }
411 : : else {
412 [ # # ][ # # ]: 0 : if (-y0 < y && y < y0) {
413 : 0 : z = y * y;
414 : 0 : x = y*(((a[3]*z+a[2])*z+a[1])*z+a[0]) /
415 : 0 : ((((b[3]*z+b[3])*z+b[1])*z+b[0])*z+1.0);
416 : : }
417 [ # # ][ # # ]: 0 : else if (y0 < y && y < 1.0) {
418 : 0 : z = sqrt(-log((1.0-y)/2.0));
419 : 0 : x = (((c[3]*z+c[2])*z+c[1])*z+c[0])/((d[1]*z+d[0])*z+1.0);
420 : : }
421 : :
422 : : // Two steps of Newton-Raphson correction to full accuracy.
423 [ # # ]: 0 : x = x - (erf (x) - y) / (2.0 / M_SQRTPI * exp (-x * x));
424 [ # # ]: 0 : x = x - (erf (x) - y) / (2.0 / M_SQRTPI * exp (-x * x));
425 : : }
426 : 0 : return x;
427 : : }
428 : :
429 : : static nr_double_t bi0_data[12] = {
430 : : -.07660547252839144951,
431 : : 1.92733795399380827000,
432 : : .22826445869203013390,
433 : : .01304891466707290428,
434 : : .00043442709008164874,
435 : : .00000942265768600193,
436 : : .00000014340062895106,
437 : : .00000000161384906966,
438 : : .00000000001396650044,
439 : : .00000000000009579451,
440 : : .00000000000000053339,
441 : : .00000000000000000245
442 : : };
443 : : static cheb_series bi0_cs = {
444 : : bi0_data, 11, -1, 1
445 : : };
446 : :
447 : : static nr_double_t ai0_data[21] = {
448 : : .07575994494023796,
449 : : .00759138081082334,
450 : : .00041531313389237,
451 : : .00001070076463439,
452 : : -.00000790117997921,
453 : : -.00000078261435014,
454 : : .00000027838499429,
455 : : .00000000825247260,
456 : : -.00000001204463945,
457 : : .00000000155964859,
458 : : .00000000022925563,
459 : : -.00000000011916228,
460 : : .00000000001757854,
461 : : .00000000000112822,
462 : : -.00000000000114684,
463 : : .00000000000027155,
464 : : -.00000000000002415,
465 : : -.00000000000000608,
466 : : .00000000000000314,
467 : : -.00000000000000071,
468 : : .00000000000000007
469 : : };
470 : : static cheb_series ai0_cs = {
471 : : ai0_data, 20, -1, 1
472 : : };
473 : :
474 : : static nr_double_t ai02_data[22] = {
475 : : .05449041101410882,
476 : : .00336911647825569,
477 : : .00006889758346918,
478 : : .00000289137052082,
479 : : .00000020489185893,
480 : : .00000002266668991,
481 : : .00000000339623203,
482 : : .00000000049406022,
483 : : .00000000001188914,
484 : : -.00000000003149915,
485 : : -.00000000001321580,
486 : : -.00000000000179419,
487 : : .00000000000071801,
488 : : .00000000000038529,
489 : : .00000000000001539,
490 : : -.00000000000004151,
491 : : -.00000000000000954,
492 : : .00000000000000382,
493 : : .00000000000000176,
494 : : -.00000000000000034,
495 : : -.00000000000000027,
496 : : .00000000000000003
497 : : };
498 : : static cheb_series ai02_cs = {
499 : : ai02_data, 21, -1, 1
500 : : };
501 : :
502 : : // Modified Bessel function of order zero.
503 : 0 : nr_double_t fspecial::i0 (nr_double_t x) {
504 : 0 : nr_double_t y = fabs (x);
505 : : nr_double_t val;
506 : :
507 [ # # ]: 0 : if (y < 2.0 * sqrt (std::numeric_limits<nr_double_t>::epsilon())) {
508 : 0 : val = 1.0;
509 : : }
510 [ # # ]: 0 : else if (y <= 3.0) {
511 : 0 : val = 2.75 + cheb_eval (&bi0_cs, y * y / 4.5 - 1.0);
512 : : }
513 [ # # ]: 0 : else if (y <= 8.0) {
514 : 0 : val = cheb_eval (&ai0_cs, (48.0 / y - 11.0) / 5.0);
515 : 0 : val = exp (y) * (0.375 + val) / sqrt (y);
516 : : }
517 : : else {
518 : 0 : val = cheb_eval (&ai02_cs, 16.0 / y - 1.0);
519 : 0 : val = exp (y) * (0.375 + val) / sqrt (y);
520 : : }
521 : 0 : return val;
522 : : }
523 : :
524 : : // Lower tail quantile for the standard normal distribution function.
525 : 0 : nr_double_t fspecial::ltqnorm (nr_double_t x) {
526 : 0 : nr_double_t q, r, e, u, z = 0.0;
527 : : static nr_double_t a[] = {
528 : : -3.969683028665376e+01, 2.209460984245205e+02,
529 : : -2.759285104469687e+02, 1.383577518672690e+02,
530 : : -3.066479806614716e+01, 2.506628277459239e+00 };
531 : : static nr_double_t b[] = {
532 : : -5.447609879822406e+01, 1.615858368580409e+02,
533 : : -1.556989798598866e+02, 6.680131188771972e+01,
534 : : -1.328068155288572e+01 };
535 : : static nr_double_t c[] = {
536 : : -7.784894002430293e-03, -3.223964580411365e-01,
537 : : -2.400758277161838e+00, -2.549732539343734e+00,
538 : : 4.374664141464968e+00, 2.938163982698783e+00 };
539 : : static nr_double_t d[] = {
540 : : 7.784695709041462e-03, 3.224671290700398e-01,
541 : : 2.445134137142996e+00, 3.754408661907416e+00 };
542 : :
543 : : // Define break-points.
544 : 0 : nr_double_t pl = 0.02425;
545 : 0 : nr_double_t ph = 1.0 - pl;
546 : :
547 : : // Rational approximation for central region:
548 [ # # ][ # # ]: 0 : if (pl <= x && x <= ph) {
549 : 0 : q = x - 0.5;
550 : 0 : r = q * q;
551 : 0 : z = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q/
552 : 0 : (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1);
553 : : }
554 : : // Rational approximation for lower region:
555 [ # # ][ # # ]: 0 : else if (0.0 < x && x < pl) {
556 : 0 : q = sqrt(-2*log(x));
557 : 0 : z = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5])/
558 : 0 : ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
559 : : }
560 : : // Rational approximation for upper region:
561 [ # # ][ # # ]: 0 : else if (ph < x && x < 1.0) {
562 : 0 : q = sqrt(-2*log(1-x));
563 : 0 : z = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5])/
564 : 0 : ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
565 : : }
566 : : // Case when X = 0:
567 [ # # ]: 0 : else if (x == 0.0) {
568 : 0 : z = -std::numeric_limits<nr_double_t>::infinity();
569 : : }
570 : : // Case when X = 1:
571 [ # # ]: 0 : else if (x == 1.0) {
572 : 0 : z = +std::numeric_limits<nr_double_t>::infinity();
573 : : }
574 : : // Cases when output will be NaN:
575 [ # # ][ # # ]: 0 : else if (x < 0.0 || x > 1.0 || std::isnan (x)) {
[ # # ][ # # ]
576 : 0 : z = +std::numeric_limits<nr_double_t>::quiet_NaN();
577 : : }
578 : :
579 : : // The relative error of the approximation has absolute value less
580 : : // than 1.15e-9. One iteration of Halley's rational method (third
581 : : // order) gives full machine precision.
582 [ # # ][ # # ]: 0 : if (0.0 < x && x < 1.0) {
583 : 0 : e = 0.5 * erfc (-z / M_SQRT2) - x; // error
584 : 0 : u = e * M_SQRT2 * M_SQRTPI * exp (z * z / 2); // f(z)/df(z)
585 : 0 : z = z - u / (1 + z * u / 2); // Halley's method
586 : : }
587 : 0 : return z;
588 : : }
589 : :
590 : : // Inverse of the error function erfc().
591 : 0 : nr_double_t fspecial::erfcinv (nr_double_t x) {
592 : 0 : return -ltqnorm (x / 2.0) / M_SQRT2;
593 : : }
|
