# HG changeset patch
# User Rik <octave@nomad.inbox5.com>
# Date 1340291117 25200
# Node ID e70a0c9cada6065e95d04b054e7e82c96c0fae3d
# Parent  566cf544d0206d0dd76da9a287cd74d3f07eb0dc
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.

diff --git a/liboctave/lo-specfun.cc b/liboctave/lo-specfun.cc
--- 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;