Mercurial > octave-nkf

changeset 14786:e70a0c9cada6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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	Tue Jun 19 19:00:55 2012 -0400
+++ b/liboctave/lo-specfun.cc	Thu Jun 21 08:05:17 2012 -0700
@@ -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;