Re: Patch for erfinv.m

[Alois Schloegl]

David Bateman wrote:
> Alois Schloegl wrote:
>   
> > Consider the following patch for erfinv.
> >
> > It (1) replaces z_old and z_new by a single variable z, and
> > (2) makes the initial checks simpler.
> >
> > The code becomes leaner, and a bit faster.
> >
> >
> > Alois
> >     
> Alois,
>
> Although this will probably apply, could you send diffs with context in
> the future. Either use the "-u" or "-c" flag to diff. This allows easier
> application of your patch to code that might have changed relative to
> the version you are working on. Apart from that giving the patch a quick
> look, it looks fine and a good idea to apply it..
>
> Regards
> David
>
>   
David,

sure, here is the patch using the -u flag against the  CVS repository.

Regards,
  Alois


Index: erfinv.m
===================================================================
RCS file: /cvs/octave/scripts/specfun/erfinv.m,v
retrieving revision 1.22
diff -u -r1.22 erfinv.m
--- erfinv.m    12 Oct 2007 21:27:26 -0000      1.22
+++ erfinv.m    21 Jan 2008 14:00:51 -0000
@@ -1,5 +1,6 @@
## Copyright (C) 1995, 1996, 1997, 1999, 2000, 2002, 2004, 2005, 2006,
##               2007 Kurt Hornik
+##               2008 Alois Schloegl
##
## This file is part of Octave.
##
@@ -44,31 +45,22 @@
  x = reshape (x, nel, 1);
  y = zeros (nel, 1);

-  i = find ((x < -1) | (x > 1) | isnan(x));
-  if any (i)
-    y(i) = NaN * ones (length (i), 1);
-  endif
-
-  t = find (x == -1);
-  y (t) = (-Inf) * ones (size (t));
-
-  t = find (x == 1);
-  y (t) = Inf * ones (size (t));
+  y(~(abs(x) < 1)) = NaN;  %% x<1, x>1, x=NaN
+  y(x == -1) = -Inf;
+  y(x == +1) = +Inf;

  i = find ((x > -1) & (x < 1));
  if any (i)
    s = sqrt (pi) / 2;
-    z_old = ones (length (i), 1);
-    z_new = sqrt (-log (1 - abs (x(i)))) .* sign (x(i));
-    while (any (abs (erf (z_new) - x(i)) > tol * abs (x(i))))
-      z_old = z_new;
-      z_new = z_old - (erf (z_old) - x(i)) .* exp (z_old.^2) * s;
+    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_new;
+    y(i) = z;
  endif

  y = reshape (y, sz);

===================================================================