Mercurial > octave-nkf

changeset 14770:cb85e836d035

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New function: erfcinv (bug #36607)
author Axel Mathéi <axel.mathei@gmail.com>
date Thu, 14 Jun 2012 23:31:59 +0200
parents 56b7aa9cd021
children 10ed11922f19
files NEWS doc/interpreter/arith.txi liboctave/lo-specfun.cc liboctave/lo-specfun.h scripts/help/unimplemented.m scripts/statistics/distributions/stdnormal_inv.m src/mappers.cc src/ov-base.cc src/ov-base.h src/ov-float.cc src/ov-flt-re-mat.cc src/ov-re-mat.cc src/ov-re-sparse.cc src/ov-scalar.cc src/ov.h
diffstat 15 files changed, 126 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS	Thu Jun 14 20:24:08 2012 -0400
+++ b/NEWS	Thu Jun 14 23:31:59 2012 +0200
@@ -77,7 +77,7 @@
       colorcube   splinefit
       lines       tetramesh
       rgbplot     shrinkfaces
-      findfigs
+      findfigs    erfcinv
 
  ** Deprecated functions.
 
--- a/doc/interpreter/arith.txi	Thu Jun 14 20:24:08 2012 -0400
+++ b/doc/interpreter/arith.txi	Thu Jun 14 23:31:59 2012 +0200
@@ -271,6 +271,8 @@
 
 @DOCSTRING(erfinv)
 
+@DOCSTRING(erfcinv)
+
 @DOCSTRING(gamma)
 
 @DOCSTRING(gammainc)
--- a/liboctave/lo-specfun.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/liboctave/lo-specfun.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -3221,6 +3221,76 @@
   return do_erfinv (x, false);
 }
 
+static double do_erfcinv (double x, bool refine)
+{
+  // Coefficients of rational approximation.
+  static const double a[] =
+    { -2.806989788730439e+01,  1.562324844726888e+02,
+      -1.951109208597547e+02,  9.783370457507161e+01,
+      -2.168328665628878e+01,  1.772453852905383e+00 };
+  static const double b[] =
+    { -5.447609879822406e+01,  1.615858368580409e+02,
+      -1.556989798598866e+02,  6.680131188771972e+01,
+      -1.328068155288572e+01 };
+  static const double c[] =
+    { -5.504751339936943e-03, -2.279687217114118e-01,
+      -1.697592457770869e+00, -1.802933168781950e+00,
+       3.093354679843505e+00,  2.077595676404383e+00 };
+  static const double d[] =
+    {  7.784695709041462e-03,  3.224671290700398e-01,
+       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;
+  double y;
+
+  // Select case.
+  if (x <= 1+pbreak && x >= 1-pbreak)
+    {
+      // Middle region.
+      const double q = 0.5*(1-x), r = q*q;
+      const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
+      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)
+    {
+      // 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;
+    }
+  else if (x == 0.0)
+    return octave_Inf;
+  else if (x == 2.0)
+    return -octave_Inf;
+  else
+    return octave_NaN;
+
+  if (refine)
+    {
+      // One iteration of Halley's method gives full precision.
+      double u = (erf(y) - (1-x)) * spi2 * exp (y*y);
+      y -= u / (1 + y*u);
+    }
+
+  return y;
+}
+
+double erfcinv (double x)
+{
+  return do_erfcinv (x, true);
+}
+
+float erfcinv (float x)
+{
+  return do_erfcinv (x, false);
+}
+
 // Implementation based on the Fortran code by W.J.Cody
 // see http://www.netlib.org/specfun/erf.
 // Templatized and simplified workflow.
--- a/liboctave/lo-specfun.h	Thu Jun 14 20:24:08 2012 -0400
+++ b/liboctave/lo-specfun.h	Thu Jun 14 23:31:59 2012 +0200
@@ -593,6 +593,9 @@
 extern OCTAVE_API double erfinv (double x);
 extern OCTAVE_API float erfinv (float x);
 
+extern OCTAVE_API double erfcinv (double x);
+extern OCTAVE_API float erfcinv (float x);
+
 extern OCTAVE_API double erfcx (double x);
 extern OCTAVE_API float erfcx (float x);
 
--- a/scripts/help/unimplemented.m	Thu Jun 14 20:24:08 2012 -0400
+++ b/scripts/help/unimplemented.m	Thu Jun 14 23:31:59 2012 +0200
@@ -169,7 +169,6 @@
   "echodemo",
   "ellipj",
   "ellipke",
-  "erfcinv",
   "errordlg",
   "evalc",
   "exifread",
--- a/scripts/statistics/distributions/stdnormal_inv.m	Thu Jun 14 20:24:08 2012 -0400
+++ b/scripts/statistics/distributions/stdnormal_inv.m	Thu Jun 14 23:31:59 2012 +0200
@@ -37,7 +37,7 @@
     error ("stdnormal_inv: X must not be complex");
   endif
 
-  inv = sqrt (2) * erfinv (2 * x - 1);
+  inv = - sqrt (2) * erfcinv (2 * x);
 
 endfunction
 
--- a/src/mappers.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/mappers.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -547,7 +547,7 @@
 @end example\n\
 \n\
 @end ifnottex\n\
-@seealso{erfc, erfcx, erfinv}\n\
+@seealso{erfc, erfcx, erfinv, erfcinv}\n\
 @end deftypefn")
 {
   octave_value retval;
@@ -596,7 +596,7 @@
 @example\n\
 erf (@var{y}) == @var{x}\n\
 @end example\n\
-@seealso{erf, erfc, erfcx}\n\
+@seealso{erf, erfc, erfcinv, erfcx}\n\
 @end deftypefn")
 {
   octave_value retval;
@@ -626,6 +626,44 @@
 %!error erfinv (1, 2)
 */
 
+DEFUN (erfcinv, args, ,
+    "-*- texinfo -*-\n\
+@deftypefn {Mapping Function} {} erfcinv (@var{x})\n\
+Compute the inverse complementary error function, i.e., @var{y} such that\n\
+\n\
+@example\n\
+erfc (@var{y}) == @var{x}\n\
+@end example\n\
+@seealso{erfinv, erf, erfc, erfcx}\n\
+@end deftypefn")
+{
+  octave_value retval;
+  if (args.length () == 1)
+    retval = args(0).erfcinv ();
+  else
+    print_usage ();
+
+  return retval;
+}
+
+/*
+## middle region
+%!assert (erfc (erfcinv ([1.9 1.3 1 0.6 0.2])), [1.9 1.3 1 0.6 0.2], eps)
+%!assert (erfc (erfcinv (single ([1.9 1.3 1 0.6 0.2]))), single ([1.9 1.3 1 0.6 0.2]), 1e-8)
+## tail region
+%!assert (erfc (erfcinv ([0.001 0.01 1.9999 1.99999])), [0.001 0.01 1.9999 1.99999], eps)
+%!assert (erfc (erfcinv (single ([0.001 0.01 1.9999 1.99999]))), single ([0.001 0.01 1.9999 1.99999]), 1e-8)
+## backward - loss of accuracy
+%!assert (erfcinv (erfc ([-3 -1 -0.4 0.7 1.3 2.8])), [-3 -1 -0.4 0.7 1.3 2.8], -1e-12)
+%!assert (erfcinv (erfc (single ([-3 -1 -0.4 0.7 1.3 2.8]))), single ([-3 -1 -0.4 0.7 1.3 2.8]), -1e-4)
+## exceptional
+%!assert (erfcinv ([2, 0, -0.1, 2.1]), [-Inf, Inf, NaN, NaN])
+%!error erfcinv (1+2i)
+
+%!error erfcinv ()
+%!error erfcinv (1, 2)
+*/
+
 DEFUN (erfc, args, ,
     "-*- texinfo -*-\n\
 @deftypefn {Mapping Function} {} erfc (@var{z})\n\
@@ -636,7 +674,7 @@
 @ifnottex\n\
 @w{@code{1 - erf (@var{z})}}.\n\
 @end ifnottex\n\
-@seealso{erfcx, erf, erfinv}\n\
+@seealso{erfcx, erf, erfinv, erfcinv}\n\
 @end deftypefn")
 {
   octave_value retval;
--- a/src/ov-base.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov-base.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -1191,6 +1191,7 @@
       "cosh",
       "erf",
       "erfinv",
+      "erfcinv",
       "erfc",
       "exp",
       "expm1",
--- a/src/ov-base.h	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov-base.h	Thu Jun 14 23:31:59 2012 +0200
@@ -688,6 +688,7 @@
       umap_cosh,
       umap_erf,
       umap_erfinv,
+      umap_erfcinv,
       umap_erfc,
       umap_erfcx,
       umap_exp,
--- a/src/ov-float.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov-float.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -291,6 +291,7 @@
       SCALAR_MAPPER (atanh, rc_atanh);
       SCALAR_MAPPER (erf, ::erff);
       SCALAR_MAPPER (erfinv, ::erfinv);
+      SCALAR_MAPPER (erfcinv, ::erfcinv);
       SCALAR_MAPPER (erfc, ::erfcf);
       SCALAR_MAPPER (erfcx, ::erfcx);
       SCALAR_MAPPER (gamma, xgamma);
--- a/src/ov-flt-re-mat.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov-flt-re-mat.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -784,6 +784,7 @@
       RC_ARRAY_MAPPER (atanh, FloatComplex, rc_atanh);
       ARRAY_MAPPER (erf, float, ::erff);
       ARRAY_MAPPER (erfinv, float, ::erfinv);
+      ARRAY_MAPPER (erfcinv, float, ::erfcinv);
       ARRAY_MAPPER (erfc, float, ::erfcf);
       ARRAY_MAPPER (erfcx, float, ::erfcx);
       ARRAY_MAPPER (gamma, float, xgamma);
--- a/src/ov-re-mat.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov-re-mat.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -907,6 +907,7 @@
       RC_ARRAY_MAPPER (atanh, Complex, rc_atanh);
       ARRAY_MAPPER (erf, double, ::erf);
       ARRAY_MAPPER (erfinv, double, ::erfinv);
+      ARRAY_MAPPER (erfcinv, double, ::erfcinv);
       ARRAY_MAPPER (erfc, double, ::erfc);
       ARRAY_MAPPER (erfcx, double, ::erfcx);
       ARRAY_MAPPER (gamma, double, xgamma);
--- a/src/ov-re-sparse.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov-re-sparse.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -910,6 +910,7 @@
       ARRAY_MAPPER (atanh, Complex, rc_atanh);
       ARRAY_MAPPER (erf, double, ::erf);
       ARRAY_MAPPER (erfinv, double, ::erfinv);
+      ARRAY_MAPPER (erfcinv, double, ::erfcinv);
       ARRAY_MAPPER (erfc, double, ::erfc);
       ARRAY_MAPPER (gamma, double, xgamma);
       ARRAY_MAPPER (lgamma, Complex, rc_lgamma);
--- a/src/ov-scalar.cc	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov-scalar.cc	Thu Jun 14 23:31:59 2012 +0200
@@ -307,6 +307,7 @@
       SCALAR_MAPPER (atanh, rc_atanh);
       SCALAR_MAPPER (erf, ::erf);
       SCALAR_MAPPER (erfinv, ::erfinv);
+      SCALAR_MAPPER (erfcinv, ::erfcinv);
       SCALAR_MAPPER (erfc, ::erfc);
       SCALAR_MAPPER (erfcx, ::erfcx);
       SCALAR_MAPPER (gamma, xgamma);
--- a/src/ov.h	Thu Jun 14 20:24:08 2012 -0400
+++ b/src/ov.h	Thu Jun 14 23:31:59 2012 +0200
@@ -1120,6 +1120,7 @@
   MAPPER_FORWARD (cosh)
   MAPPER_FORWARD (erf)
   MAPPER_FORWARD (erfinv)
+  MAPPER_FORWARD (erfcinv)
   MAPPER_FORWARD (erfc)
   MAPPER_FORWARD (erfcx)
   MAPPER_FORWARD (exp)