Mercurial > hg > octave-nkf
changeset 14786:e70a0c9cada6
Pre-compute bounds (constant folding) for erfcinv function.
* lo-specfun.cc (erfcinv): Do constant folding for ranges of function
validity. Add explanation about algorithm source.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Thu, 21 Jun 2012 08:05:17 -0700 |
parents | 566cf544d020 |
children | acb09716fc94 |
files | liboctave/lo-specfun.cc |
diffstat | 1 files changed, 10 insertions(+), 6 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/lo-specfun.cc +++ b/liboctave/lo-specfun.cc @@ -3221,6 +3221,10 @@ return do_erfinv (x, false); } +// The algorthim for erfcinv is an adaptation of the erfinv algorithm above +// from P. J. Acklam. It has been modified to run over the different input +// domain of erfcinv. See the notes for erfinv for an explanation. + static double do_erfcinv (double x, bool refine) { // Coefficients of rational approximation. @@ -3241,12 +3245,12 @@ 2.445134137142996e+00, 3.754408661907416e+00 }; static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. - static const double pi = 3.14159265358979323846; - static const double pbreak = 0.95150; + static const double pbreak_lo = 0.04850; // 1-pbreak + static const double pbreak_hi = 1.95150; // 1+pbreak double y; // Select case. - if (x <= 1+pbreak && x >= 1-pbreak) + if (x >= pbreak_lo && x <= pbreak_hi) { // Middle region. const double q = 0.5*(1-x), r = q*q; @@ -3254,15 +3258,15 @@ const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0; y = yn / yd; } - else if (x < 2.0 && x > 0.0) + else if (x > 0.0 && x < 2.0) { // Tail region. const double q = x < 1 ? sqrt (-2*log (0.5*x)) : sqrt (-2*log (0.5*(2-x))); const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]; const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0; y = yn / yd; - if (x < 1-pbreak) - y *= -1; + if (x < pbreak_lo) + y = -y; } else if (x == 0.0) return octave_Inf;