Mercurial > hg > octave-terminal
changeset 3266:3add30d33bd0
[project @ 1999-09-03 17:15:40 by jwe]
author | jwe |
---|---|
date | Fri, 03 Sep 1999 17:15:41 +0000 |
parents | a2b3a1413d28 |
children | 8092e8197ce4 |
files | scripts/ChangeLog scripts/specfun/erfinv.m |
diffstat | 2 files changed, 21 insertions(+), 7 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +Fri Sep 3 00:01:38 1999 John W. Eaton <jwe@bevo.che.wisc.edu> + + * specfun/erfinv.m: Improve stopping criterion. + Add iteration count as second return value. + Mon Aug 30 12:07:00 1999 John W. Eaton <jwe@bevo.che.wisc.edu> * statistics/base/mean.m: Use .', not ' to reorient row vectors.
--- a/scripts/specfun/erfinv.m +++ b/scripts/specfun/erfinv.m @@ -1,4 +1,4 @@ -## Copyright (C) 1995, 1996, 1997 Kurt Hornik +## Copyright (C) 1995, 1996 Kurt Hornik ## ## This program is free software; you can redistribute it and/or modify ## it under the terms of the GNU General Public License as published by @@ -22,16 +22,21 @@ ## Created: 27 September 1994 ## Adapted-By: jwe -function y = erfinv (x) +function [y, iterations] = erfinv (x) if (nargin != 1) usage ("erfinv (x)"); endif - + + maxit = 100; + tol = eps; + + iterations = 0; + [m, n] = size (x); x = reshape (x, m * n, 1); y = zeros (m * n, 1); - + i = find ((x < -1) | (x > 1)); if any (i) y(i) = NaN * ones (length (i), 1); @@ -42,15 +47,19 @@ t = find (x == 1); y (t) = Inf * ones (size (t)); - + i = find ((x > -1) & (x < 1)); if any (i) s = sqrt (pi) / 2; z_old = ones (length (i), 1); - z_new = zeros (length (i), 1); - while (any (any (abs (z_new - z_old) > 2 * eps))) + z_new = sqrt (-log (1 - abs (x(i)))) .* sign (x(i)); + while (any (abs (erf (z_old) - x(i)) > tol * abs (x(i)))) z_old = z_new; z_new = z_old - (erf (z_old) - x(i)) .* exp (z_old.^2) * s; + if (++iterations > maxit) + warning ("erfinv: iteration limit exceeded"); + break; + endif endwhile y(i) = z_new; endif