Mercurial > octave-nkf
changeset 9835:1bb1ed717d2f
implement built-in erfinv
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Thu, 19 Nov 2009 10:30:57 +0100 |
| parents | 92d8f35ff217 |
| children | 804c21f3659b |
| files | liboctave/ChangeLog liboctave/lo-specfun.cc liboctave/lo-specfun.h scripts/ChangeLog scripts/specfun/erfinv.m scripts/specfun/module.mk src/ChangeLog 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 | 16 files changed, 140 insertions(+), 74 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog Wed Nov 18 23:14:09 2009 -0500 +++ b/liboctave/ChangeLog Thu Nov 19 10:30:57 2009 +0100 @@ -1,3 +1,9 @@ +2009-11-19 Jaroslav Hajek <highegg@gmail.com> + + * lo-mappers.cc (do_erfinv): New static function. + (erfinv (double), erfinv (float)): New mappers. + * lo-mappers.h: Declare them. + 2009-11-18 David Grundberg <davidg@cs.umu.se> * str-vec.cc (string_vector::list_in_columns): Avoid crash on
--- a/liboctave/lo-specfun.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/liboctave/lo-specfun.cc Thu Nov 19 10:30:57 2009 +0100 @@ -3101,6 +3101,78 @@ return x < -1.0f ? FloatComplex (logf (-(1.0f + x)), pi) : FloatComplex (log1pf (x)); } +// This algorithm is due to P. Jacklam +// See http:// +// The rational approximation has relative accuracy 1e-9 in the whole region. +// For doubles, it is refined by a single step of Higham's 3rd order method. +// For single precision, the accuracy is already OK, so we skip it to get +// faster evaluation. + +static double do_erfinv (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 pbreak = 0.95150; + double ax = fabs (x), y; + + // Select case. + if (ax <= pbreak) + { + // Middle region. + const double q = 0.5 * 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 (ax < 1.0) + { + // Tail region. + const double q = sqrt (-2*log (0.5*(1-ax))); + 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 * signum (-x); + } + else if (ax == 1.0) + return octave_Inf * signum (x); + else + return octave_NaN; + + if (refine) + { + // One iteration of Halley's method gives full precision. + double u = (erf(y) - x) * spi2 * exp (y*y); + y -= u / (1 + y*u); + } + + return y; +} + +double erfinv (double x) +{ + return do_erfinv (x, true); +} + +float erfinv (float x) +{ + return do_erfinv (x, false); +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/lo-specfun.h Wed Nov 18 23:14:09 2009 -0500 +++ b/liboctave/lo-specfun.h Thu Nov 19 10:30:57 2009 +0100 @@ -583,6 +583,9 @@ extern OCTAVE_API Complex rc_log1p (double); extern OCTAVE_API FloatComplex rc_log1p (float); +extern OCTAVE_API double erfinv (double x); +extern OCTAVE_API float erfinv (float x); + #endif /*
--- a/scripts/ChangeLog Wed Nov 18 23:14:09 2009 -0500 +++ b/scripts/ChangeLog Thu Nov 19 10:30:57 2009 +0100 @@ -1,3 +1,8 @@ +2009-11-19 Jaroslav Hajek <highegg@gmail.com> + + * specfun/erfinv.m: Remove. + * specfun/module.mk: Update. + 2009-11-18 Ben Abbott <bpabbott@mac.com> * plot/orient.m: Flip papersize and paperposition when orientation
--- a/scripts/specfun/erfinv.m Wed Nov 18 23:14:09 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,73 +0,0 @@ -## Copyright (C) 1995, 1996, 1997, 1999, 2000, 2002, 2004, 2005, 2006, -## 2007 Kurt Hornik -## Copyright (C) 2008 Alois Schloegl -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Mapping Function} {} erfinv (@var{z}) -## Computes the inverse of the error function. -## @seealso{erf, erfc} -## @end deftypefn - -## Author: KH <Kurt.Hornik@wu-wien.ac.at> -## Created: 27 September 1994 -## Adapted-By: jwe - -function [y, iterations] = erfinv (x) - - if (nargin != 1) - print_usage (); - endif - - maxit = 100; - if (isa (x, "single")) - tol = eps ("single"); - else - tol = eps; - endif - - iterations = 0; - - sz = size (x); - nel = numel (x); - - x = reshape (x, nel, 1); - y = zeros (nel, 1); - - ## x < -1 or x > 1 ==> NaN - y(abs (x) >= 1) = NaN; - y(x == -1) = -Inf; - y(x == +1) = +Inf; - - i = find ((x > -1) & (x < 1)); - if (any (i)) - s = sqrt (pi) / 2; - z = sqrt (-log (1 - abs (x(i)))) .* sign (x(i)); - while (any (abs (erf (z) - x(i)) > tol * abs (x(i)))) - z = z - (erf (z) - x(i)) .* exp (z.^2) * s; - if (++iterations > maxit) - warning ("erfinv: iteration limit exceeded"); - break; - endif - endwhile - y(i) = z; - endif - - y = reshape (y, sz); - -endfunction
--- a/scripts/specfun/module.mk Wed Nov 18 23:14:09 2009 -0500 +++ b/scripts/specfun/module.mk Thu Nov 19 10:30:57 2009 +0100 @@ -5,7 +5,6 @@ specfun/beta.m \ specfun/betai.m \ specfun/betaln.m \ - specfun/erfinv.m \ specfun/factor.m \ specfun/factorial.m \ specfun/gammai.m \
--- a/src/ChangeLog Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ChangeLog Thu Nov 19 10:30:57 2009 +0100 @@ -1,3 +1,15 @@ +2009-11-19 Jaroslav Hajek <highegg@gmail.com> + + * ov-base.h (unary_mapper_t): New member: umap_erfinv. + * ov-base.c (octave_base_value::get_umap_name): Add "erfinv" here. + * ov.h (octave_value::erfinv): New method. + * ov-scalar.cc (octave_scalar::map): Handle umap_erfinv. + * ov-float.cc (octave_float::map): Ditto. + * ov-re-mat.cc (octave_matrix::map): Ditto. + * ov-flt-re-mat.cc (octave_float_matrix::map): Ditto. + * ov-re-sparse.cc (octave_sparse_matrix::map): Ditto. + * mappers.cc (Ferfinv): New DEFUN. + 2009-11-14 Shai Ayal <shaiay@users.sourceforge.net> * gl-render.cc (opengl_renderer::text_to_pixels):
--- a/src/mappers.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/src/mappers.cc Thu Nov 19 10:30:57 2009 +0100 @@ -568,6 +568,40 @@ */ +DEFUN (erfinv, args, , + "-*- texinfo -*-\n\ +@deftypefn {Mapping Function} {} erfinv (@var{x})\n\ +Computes the inverse error function, i.e. @var{y} such that\n\ +@example\n\ + erf(@var{y}) == @var{x}\n\ +@end example\n\ +@seealso{erfc, erf}\n\ +@end deftypefn") +{ + octave_value retval; + if (args.length () == 1) + retval = args(0).erfinv (); + else + print_usage (); + + return retval; +} + +/* +%% middle region +%!assert (erf (erfinv ([-0.9 -0.3 0 0.4 0.8])), [-0.9 -0.3 0 0.4 0.8], 1e-16) +%!assert (erf (erfinv (single ([-0.9 -0.3 0 0.4 0.8]))), single ([-0.9 -0.3 0 0.4 0.8]), 1e-8) +%% tail region +%!assert (erf (erfinv ([-0.999 -0.99 0.9999 0.99999])), [-0.999 -0.99 0.9999 0.99999], 1e-16) +%!assert (erf (erfinv (single ([-0.999 -0.99 0.9999 0.99999]))), single ([-0.999 -0.99 0.9999 0.99999]), 1e-8) +%% backward - loss of accuracy +%!assert (erfinv (erf ([-3 -1 -0.4 0.7 1.3 2.8])), [-3 -1 -0.4 0.7 1.3 2.8], -1e-12) +%!assert (erfinv (erf (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 (erfinv ([-1, 1, 1.1, -2.1]), [-Inf, Inf, NaN, NaN]) +%!error erfinv (1+2i) +*/ + DEFUN (erfc, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} erfc (@var{z})\n\
--- a/src/ov-base.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ov-base.cc Thu Nov 19 10:30:57 2009 +0100 @@ -1137,6 +1137,7 @@ "cos", "cosh", "erf", + "erfinv", "erfc", "exp", "expm1",
--- a/src/ov-base.h Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ov-base.h Thu Nov 19 10:30:57 2009 +0100 @@ -649,6 +649,7 @@ umap_cos, umap_cosh, umap_erf, + umap_erfinv, umap_erfc, umap_exp, umap_expm1,
--- a/src/ov-float.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ov-float.cc Thu Nov 19 10:30:57 2009 +0100 @@ -278,6 +278,7 @@ SCALAR_MAPPER (atan, ::atanf); SCALAR_MAPPER (atanh, rc_atanh); SCALAR_MAPPER (erf, ::erff); + SCALAR_MAPPER (erfinv, ::erfinv); SCALAR_MAPPER (erfc, ::erfcf); SCALAR_MAPPER (gamma, xgamma); SCALAR_MAPPER (lgamma, rc_lgamma);
--- a/src/ov-flt-re-mat.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ov-flt-re-mat.cc Thu Nov 19 10:30:57 2009 +0100 @@ -758,6 +758,7 @@ ARRAY_MAPPER (atan, float, ::atanf); RC_ARRAY_MAPPER (atanh, FloatComplex, rc_atanh); ARRAY_MAPPER (erf, float, ::erff); + ARRAY_MAPPER (erfinv, float, ::erfinv); ARRAY_MAPPER (erfc, float, ::erfcf); ARRAY_MAPPER (gamma, float, xgamma); RC_ARRAY_MAPPER (lgamma, FloatComplex, rc_lgamma);
--- a/src/ov-re-mat.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ov-re-mat.cc Thu Nov 19 10:30:57 2009 +0100 @@ -786,6 +786,7 @@ ARRAY_MAPPER (atan, double, ::atan); RC_ARRAY_MAPPER (atanh, Complex, rc_atanh); ARRAY_MAPPER (erf, double, ::erf); + ARRAY_MAPPER (erfinv, double, ::erfinv); ARRAY_MAPPER (erfc, double, ::erfc); ARRAY_MAPPER (gamma, double, xgamma); RC_ARRAY_MAPPER (lgamma, Complex, rc_lgamma);
--- a/src/ov-re-sparse.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ov-re-sparse.cc Thu Nov 19 10:30:57 2009 +0100 @@ -840,6 +840,7 @@ ARRAY_MAPPER (atan, double, ::atan); ARRAY_MAPPER (atanh, Complex, rc_atanh); ARRAY_MAPPER (erf, double, ::erf); + ARRAY_MAPPER (erfinv, double, ::erfinv); ARRAY_MAPPER (erfc, double, ::erfc); ARRAY_MAPPER (gamma, double, xgamma); ARRAY_MAPPER (lgamma, Complex, rc_lgamma);
--- a/src/ov-scalar.cc Wed Nov 18 23:14:09 2009 -0500 +++ b/src/ov-scalar.cc Thu Nov 19 10:30:57 2009 +0100 @@ -293,6 +293,7 @@ SCALAR_MAPPER (atan, ::atan); SCALAR_MAPPER (atanh, rc_atanh); SCALAR_MAPPER (erf, ::erf); + SCALAR_MAPPER (erfinv, ::erfinv); SCALAR_MAPPER (erfc, ::erfc); SCALAR_MAPPER (gamma, xgamma); SCALAR_MAPPER (lgamma, rc_lgamma);