# HG changeset patch
# User jwe
# Date 1201037664 0
# Node ID 8a3b2ccc4e11a0cfeea480b17c8975cc3f9e02c8
# Parent  73036cdd855d8769dfecee1447d33b47d00f6a4a
[project @ 2008-01-22 21:34:24 by jwe]

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,8 @@
+2008-01-22  Schloegl Alois  <alois.schloegl@tugraz.at>
+
+	* specfun/erfinv.m: Replace z_old and z_new by a single variable z.
+	Simplify initial checks on argument values.
+
 2008-01-22  Michael Goffioul  <michael.goffioul@gmail.com>
 
 	* plot/gnuplot_drawnow.m: New function corresponding to the
diff --git a/scripts/specfun/erfinv.m b/scripts/specfun/erfinv.m
--- a/scripts/specfun/erfinv.m
+++ b/scripts/specfun/erfinv.m
@@ -1,5 +1,6 @@
 ## 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.
 ##
@@ -44,31 +45,23 @@
   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));
+  ## 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)
+  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);