home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / ERFINV.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.5 KB  |  59 lines
  1. function x = erfinv(y);
  2. %ERFINV    Inverse of the error function.
  3. %    x = erfinv(y) satisfies y = erf(x),
  4. %        -1 <= y < 1, -inf <= x <= inf.
  5. %
  6. %    See also ERF.
  7.  
  8. %    Cleve Moler, 12-31-90, 5-7-91, 7-25-91.
  9. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  10.  
  11. x = NaN*ones(size(y));
  12.  
  13. % Coefficients in rational approximations.
  14.  
  15. a = [ 0.886226899 -1.645349621  0.914624893 -0.140543331];
  16. b = [-2.118377725  1.442710462 -0.329097515  0.012229801];
  17. c = [-1.970840454 -1.624906493  3.429567803  1.641345311];
  18. d = [ 3.543889200  1.637067800];
  19.  
  20. % Central range.
  21.  
  22. y0 = .7;
  23. k = find(abs(y) <= y0);
  24. if any(k)
  25.     z = y(k).*y(k);
  26.     x(k) = y(k) .* (((a(4)*z+a(3)).*z+a(2)).*z+a(1)) ./ ...
  27.          ((((b(4)*z+b(3)).*z+b(2)).*z+b(1)).*z+1);
  28. end;
  29.  
  30. % Near end points of range.
  31.  
  32. k = find(( y0 < y ) & (y <  1));
  33. if any(k)
  34.     z = sqrt(-log((1-y(k))/2));
  35.     x(k) = (((c(4)*z+c(3)).*z+c(2)).*z+c(1)) ./ ((d(2)*z+d(1)).*z+1);
  36. end
  37.  
  38. k = find((-y0 > y ) & (y > -1));
  39. if any(k)
  40.     z = sqrt(-log((1+y(k))/2));
  41.     x(k) = -(((c(4)*z+c(3)).*z+c(2)).*z+c(1)) ./ ((d(2)*z+d(1)).*z+1);
  42. end
  43.  
  44. % Two steps of Newton-Raphson correction to full accuracy.
  45. % Without these steps, erfinv(y) would be about 3 times
  46. % faster to compute, but accurate to only about 6 digits.
  47.  
  48. x = x - (erf(x) - y) ./ (2/sqrt(pi) * exp(-x.^2));
  49. x = x - (erf(x) - y) ./ (2/sqrt(pi) * exp(-x.^2));
  50.  
  51. % Exceptional cases.
  52.  
  53. k = find(y == -1);
  54. if any(k), x(k) = -inf*k; end
  55. k = find(y == 1);
  56. if any(k), x(k) = inf*k; end
  57. k = find(abs(y) > 1);
  58. if any(k), x(k) = NaN*k; end
  59.