Mercurial > octave-nkf

comparison scripts/specfun/erfinv.m @ 7410:8a3b2ccc4e11

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2008-01-22 21:34:24 by jwe]
author jwe
date Tue, 22 Jan 2008 21:34:24 +0000
parents a1dbe9d80eee
children ce6adb34ecf8
comparison
equal deleted inserted replaced
7409:73036cdd855d 7410:8a3b2ccc4e11
1 ## Copyright (C) 1995, 1996, 1997, 1999, 2000, 2002, 2004, 2005, 2006, 1 ## Copyright (C) 1995, 1996, 1997, 1999, 2000, 2002, 2004, 2005, 2006,
2 ## 2007 Kurt Hornik 2 ## 2007 Kurt Hornik
3 ## Copyright (C) 2008 Alois Schloegl
3 ## 4 ##
4 ## This file is part of Octave. 5 ## This file is part of Octave.
5 ## 6 ##
6 ## Octave is free software; you can redistribute it and/or modify it 7 ## Octave is free software; you can redistribute it and/or modify it
7 ## under the terms of the GNU General Public License as published by 8 ## under the terms of the GNU General Public License as published by
42 nel = numel (x); 43 nel = numel (x);
43 44
44 x = reshape (x, nel, 1); 45 x = reshape (x, nel, 1);
45 y = zeros (nel, 1); 46 y = zeros (nel, 1);
46 47
47 i = find ((x < -1) | (x > 1) | isnan(x)); 48 ## x < 1 or x > 1 ==> NaN
48 if any (i) 49 y(abs (x) >= 1) = NaN;
49 y(i) = NaN * ones (length (i), 1); 50 y(x == -1) = -Inf;
50 endif 51 y(x == +1) = +Inf;
51
52 t = find (x == -1);
53 y (t) = (-Inf) * ones (size (t));
54
55 t = find (x == 1);
56 y (t) = Inf * ones (size (t));
57 52
58 i = find ((x > -1) & (x < 1)); 53 i = find ((x > -1) & (x < 1));
59 if any (i) 54 if (any (i))
60 s = sqrt (pi) / 2; 55 s = sqrt (pi) / 2;
61 z_old = ones (length (i), 1); 56 z = sqrt (-log (1 - abs (x(i)))) .* sign (x(i));
62 z_new = sqrt (-log (1 - abs (x(i)))) .* sign (x(i)); 57 while (any (abs (erf (z) - x(i)) > tol * abs (x(i))))
63 while (any (abs (erf (z_new) - x(i)) > tol * abs (x(i)))) 58 z = z - (erf (z) - x(i)) .* exp (z.^2) * s;
64 z_old = z_new;
65 z_new = z_old - (erf (z_old) - x(i)) .* exp (z_old.^2) * s;
66 if (++iterations > maxit) 59 if (++iterations > maxit)
67 warning ("erfinv: iteration limit exceeded"); 60 warning ("erfinv: iteration limit exceeded");
68 break; 61 break;
69 endif 62 endif
70 endwhile 63 endwhile
71 y(i) = z_new; 64 y(i) = z;
72 endif 65 endif
73 66
74 y = reshape (y, sz); 67 y = reshape (y, sz);
75 68
76 endfunction 69 endfunction