Mercurial > octave-nkf

changeset 10391:59e34bcdff13

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement built-in erfcx
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 04 Mar 2010 09:35:38 +0100
parents ad0b54ae206a
children b4e5dcf023c9
files liboctave/ChangeLog liboctave/lo-specfun.cc liboctave/lo-specfun.h src/ChangeLog src/mappers.cc src/ov-base.h src/ov-float.cc src/ov-flt-re-mat.cc src/ov-re-mat.cc src/ov-scalar.cc src/ov.h
diffstat 11 files changed, 150 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog	Wed Mar 03 13:01:44 2010 -0500
+++ b/liboctave/ChangeLog	Thu Mar 04 09:35:38 2010 +0100
@@ -1,3 +1,8 @@
+2010-03-04  Jaroslav Hajek  <highegg@gmail.com>
+
+	* lo-specfun.cc (erfcx, erfcx_impl): New functions.
+	* lo-specfun.h: Declare erfcx.
+
 2010-03-03  John W. Eaton  <jwe@octave.org>
 
 	* oct-convn.cc (convolve): Cast int constant to octave_idx_type in
--- a/liboctave/lo-specfun.cc	Wed Mar 03 13:01:44 2010 -0500
+++ b/liboctave/lo-specfun.cc	Thu Mar 04 09:35:38 2010 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1996, 1998, 2002, 2003, 2004, 2005, 2006, 2007, 2008
               John W. Eaton
+Copyright (C) 2010 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -3170,3 +3171,103 @@
 {
   return do_erfinv (x, false); 
 }
+
+// Implementation based on the Fortran code by W.J.Cody
+// see http://www.netlib.org/specfun/erf.
+// Templatized and simplified workflow.
+
+// FIXME: Maybe this should be globally visible.
+static inline float erfc (float x) { return erfcf (x); }
+
+template <class T>
+static T 
+erfcx_impl (T x)
+{
+  static const T c[] = 
+    { 
+      5.64188496988670089e-1,8.88314979438837594,
+      6.61191906371416295e+1,2.98635138197400131e+2,
+      8.81952221241769090e+2,1.71204761263407058e+3,
+      2.05107837782607147e+3,1.23033935479799725e+3,
+      2.15311535474403846e-8 
+    };
+
+  static const T d[] = 
+    { 
+      1.57449261107098347e+1,1.17693950891312499e+2,
+      5.37181101862009858e+2,1.62138957456669019e+3,
+      3.29079923573345963e+3,4.36261909014324716e+3,
+      3.43936767414372164e+3,1.23033935480374942e+3
+    };
+
+  static const T p[] = 
+    {
+      3.05326634961232344e-1,3.60344899949804439e-1,
+      1.25781726111229246e-1,1.60837851487422766e-2,
+      6.58749161529837803e-4,1.63153871373020978e-2
+    };
+
+  static const T q[] = 
+    {
+      2.56852019228982242,1.87295284992346047,
+      5.27905102951428412e-1,6.05183413124413191e-2,
+      2.33520497626869185e-3
+    };
+
+  static const T sqrpi = 5.6418958354775628695e-1;
+  static const T xhuge = sqrt (1.0 / std::numeric_limits<T>::epsilon ());
+  static const T xneg = -sqrt (log (std::numeric_limits<T>::max ()/2.0));
+
+  double y = fabs (x), result;
+  if (x < xneg)
+    result = octave_Inf;
+  else if (y <= 0.46875)
+    result = std::exp (x*x) * erfc (x);
+  else
+    {
+      if (y <= 4.0)
+        {
+          double xnum = c[8]*y, xden = y;
+          for (int i = 0; i < 7; i++)
+            {
+              xnum = (xnum + c[i]) * y;
+              xden = (xden + d[i]) * y;
+            }
+
+          result = (xnum + c[7]) / (xden + d[7]);
+        }
+      else if (y <= xhuge)
+        {
+          double y2 = 1/(y*y), xnum = p[5]*y2, xden = y2;
+          for (int i = 0; i < 4; i++)
+            {
+              xnum = (xnum + p[i]) * y2;
+              xden = (xden + q[i]) * y2;
+            }
+
+          result = y2 * (xnum + p[4]) / (xden + q[4]);
+          result = (sqrpi - result) / y;
+        }
+      else
+        result = sqrpi / y;
+
+      // Fix up negative argument.
+      if (x < 0)
+        {
+          double y2 = ceil (x / 16.0) * 16.0, del = (x-y2)*(x+y2);
+          result = 2*(std::exp(y2*y2) * std::exp(del)) - result;
+        }
+    }
+
+  return result;
+}
+
+double erfcx (double x)
+{
+  return erfcx_impl (x);
+}
+
+float erfcx (float x)
+{
+  return erfcx_impl (x);
+}
--- a/liboctave/lo-specfun.h	Wed Mar 03 13:01:44 2010 -0500
+++ b/liboctave/lo-specfun.h	Thu Mar 04 09:35:38 2010 +0100
@@ -585,4 +585,7 @@
 extern OCTAVE_API double erfinv (double x);
 extern OCTAVE_API float erfinv (float x);
 
+extern OCTAVE_API double erfcx (double x);
+extern OCTAVE_API float erfcx (float x);
+
 #endif
--- a/src/ChangeLog	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/ChangeLog	Thu Mar 04 09:35:38 2010 +0100
@@ -1,3 +1,12 @@
+2010-03-04  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-base.h (unary_mapper_t::umap_erfcx): New enum member.
+	* ov.h (octave_value::erfcx): New forwarder method.
+	* ov-scalar.cc (octave_scalar::map): Handle erfcx.
+	* ov-float.cc (octave_float_scalar::map): Ditto.
+	* ov-re-mat.cc (octave_matrix::map): Ditto.
+	* ov-flt-re-mat.cc (octave_float_matrix::map): Ditto.
+
 2010-03-03  John W. Eaton  <jwe@octave.org>
 
 	* DLD-FUNCTIONS/__magick_read__.cc (F__magick_read__):
--- a/src/mappers.cc	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/mappers.cc	Thu Mar 04 09:35:38 2010 +0100
@@ -632,6 +632,32 @@
 
 */
 
+DEFUN (erfcx, args, ,
+    "-*- texinfo -*-\n\
+@deftypefn {Mapping Function} {} erfcx (@var{z})\n\
+Computes the scaled complementary error function,\n\
+@tex\n\
+$z^2 (1 - {\\rm erf} (z))$.\n\
+@end tex\n\
+@ifnottex\n\
+@code{z^2*(1 - erf (@var{z}))}.\n\
+@end ifnottex\n\
+@seealso{erfc, erf, erfinv}\n\
+@end deftypefn")
+{
+  octave_value retval;
+  if (args.length () == 1)
+    retval = args(0).erfcx ();
+  else
+    print_usage ();
+
+  return retval;
+}
+
+/*
+
+*/
+
 DEFUN (exp, args, ,
     "-*- texinfo -*-\n\
 @deftypefn {Mapping Function} {} exp (@var{x})\n\
--- a/src/ov-base.h	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/ov-base.h	Thu Mar 04 09:35:38 2010 +0100
@@ -665,6 +665,7 @@
       umap_erf,
       umap_erfinv,
       umap_erfc,
+      umap_erfcx,
       umap_exp,
       umap_expm1,
       umap_finite,
--- a/src/ov-float.cc	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/ov-float.cc	Thu Mar 04 09:35:38 2010 +0100
@@ -287,6 +287,7 @@
       SCALAR_MAPPER (erf, ::erff);
       SCALAR_MAPPER (erfinv, ::erfinv);
       SCALAR_MAPPER (erfc, ::erfcf);
+      SCALAR_MAPPER (erfcx, ::erfcx);
       SCALAR_MAPPER (gamma, xgamma);
       SCALAR_MAPPER (lgamma, rc_lgamma);
       SCALAR_MAPPER (ceil, ::ceilf);
--- a/src/ov-flt-re-mat.cc	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/ov-flt-re-mat.cc	Thu Mar 04 09:35:38 2010 +0100
@@ -767,6 +767,7 @@
       ARRAY_MAPPER (erf, float, ::erff);
       ARRAY_MAPPER (erfinv, float, ::erfinv);
       ARRAY_MAPPER (erfc, float, ::erfcf);
+      ARRAY_MAPPER (erfcx, float, ::erfcx);
       ARRAY_MAPPER (gamma, float, xgamma);
       RC_ARRAY_MAPPER (lgamma, FloatComplex, rc_lgamma);
       ARRAY_MAPPER (ceil, float, ::ceilf);
--- a/src/ov-re-mat.cc	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/ov-re-mat.cc	Thu Mar 04 09:35:38 2010 +0100
@@ -890,6 +890,7 @@
       ARRAY_MAPPER (erf, double, ::erf);
       ARRAY_MAPPER (erfinv, double, ::erfinv);
       ARRAY_MAPPER (erfc, double, ::erfc);
+      ARRAY_MAPPER (erfcx, double, ::erfcx);
       ARRAY_MAPPER (gamma, double, xgamma);
       RC_ARRAY_MAPPER (lgamma, Complex, rc_lgamma);
       ARRAY_MAPPER (ceil, double, ::ceil);
--- a/src/ov-scalar.cc	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/ov-scalar.cc	Thu Mar 04 09:35:38 2010 +0100
@@ -303,6 +303,7 @@
       SCALAR_MAPPER (erf, ::erf);
       SCALAR_MAPPER (erfinv, ::erfinv);
       SCALAR_MAPPER (erfc, ::erfc);
+      SCALAR_MAPPER (erfcx, ::erfcx);
       SCALAR_MAPPER (gamma, xgamma);
       SCALAR_MAPPER (lgamma, rc_lgamma);
       SCALAR_MAPPER (ceil, ::ceil);
--- a/src/ov.h	Wed Mar 03 13:01:44 2010 -0500
+++ b/src/ov.h	Thu Mar 04 09:35:38 2010 +0100
@@ -1089,6 +1089,7 @@
   MAPPER_FORWARD (erf)
   MAPPER_FORWARD (erfinv)
   MAPPER_FORWARD (erfc)
+  MAPPER_FORWARD (erfcx)
   MAPPER_FORWARD (exp)
   MAPPER_FORWARD (expm1)
   MAPPER_FORWARD (finite)