home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / ERFCORE.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  4.2 KB  |  121 lines
  1. function result = erfcore(x,jint)
  2. %ERFCORE This function is used by erf, erfc, and erfcx.
  3. %    erf(x) = erfcore(x,0)
  4. %    erfc(x) = erfcore(x,1)
  5. %    erfcx(x) = exp(x^2)*erfc(x) = erfcore(x,2)
  6.  
  7. %    C. Moler, 2-1-91.
  8. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  9.  
  10. %   This is a translation of a FORTRAN program by W. J. Cody,
  11. %   Argonne National Laboratory, NETLIB/SPECFUN, March 19, 1990.
  12. %   The main computation evaluates near-minimax approximations
  13. %   from "Rational Chebyshev approximations for the error function"
  14. %   by W. J. Cody, Math. Comp., 1969, PP. 631-638.
  15.  
  16.     if any(imag(x) ~= 0)
  17.        error('Input argument must be real.')
  18.     end
  19.     result = NaN*ones(size(x));
  20. %
  21. %   evaluate  erf  for  |x| <= 0.46875
  22. %
  23.     xbreak = 0.46875;
  24.     k = find(abs(x) <= xbreak);
  25.     if any(k)
  26.         a = [3.16112374387056560e00; 1.13864154151050156e02;
  27.              3.77485237685302021e02; 3.20937758913846947e03;
  28.              1.85777706184603153e-1];
  29.         b = [2.36012909523441209e01; 2.44024637934444173e02;
  30.              1.28261652607737228e03; 2.84423683343917062e03];
  31.  
  32.             y = abs(x(k));
  33.             z = y .* y;
  34.             xnum = a(5)*z;
  35.             xden = z;
  36.             for i = 1:3
  37.                xnum = (xnum + a(i)) .* z;
  38.                xden = (xden + b(i)) .* z;
  39.             end
  40.             result(k) = x(k) .* (xnum + a(4)) ./ (xden + b(4));
  41.             if jint ~= 0, result(k) = 1 - result(k); end
  42.             if jint == 2, result(k) = exp(z) .* result(k); end
  43.     end
  44. %
  45. %   evaluate  erfc  for 0.46875 <= |x| <= 4.0
  46. %
  47.     k = find((abs(x) > xbreak) &  (abs(x) <= 4.));
  48.     if any(k)
  49.         c = [5.64188496988670089e-1; 8.88314979438837594e00;
  50.              6.61191906371416295e01; 2.98635138197400131e02;
  51.              8.81952221241769090e02; 1.71204761263407058e03;
  52.              2.05107837782607147e03; 1.23033935479799725e03;
  53.              2.15311535474403846e-8];
  54.         d = [1.57449261107098347e01; 1.17693950891312499e02;
  55.              5.37181101862009858e02; 1.62138957456669019e03;
  56.              3.29079923573345963e03; 4.36261909014324716e03;
  57.              3.43936767414372164e03; 1.23033935480374942e03];
  58.  
  59.             y = abs(x(k));
  60.             xnum = c(9)*y;
  61.             xden = y;
  62.             for i = 1:7
  63.                xnum = (xnum + c(i)) .* y;
  64.                xden = (xden + d(i)) .* y;
  65.             end
  66.             result(k) = (xnum + c(8)) ./ (xden + d(8));
  67.             if jint ~= 2
  68.                z = fix(y*16)/16;
  69.                del = (y-z).*(y+z);
  70.                result(k) = exp(-z.*z) .* exp(-del) .* result(k);
  71.             end
  72.     end
  73. %
  74. %   evaluate  erfc  for |x| > 4.0
  75. %
  76.     k = find(abs(x) > 4.0);
  77.     if any(k)
  78.         p = [3.05326634961232344e-1; 3.60344899949804439e-1;
  79.              1.25781726111229246e-1; 1.60837851487422766e-2;
  80.              6.58749161529837803e-4; 1.63153871373020978e-2];
  81.         q = [2.56852019228982242e00; 1.87295284992346047e00;
  82.              5.27905102951428412e-1; 6.05183413124413191e-2;
  83.              2.33520497626869185e-3];
  84.  
  85.             y = abs(x(k));
  86.             z = 1 ./ (y .* y);
  87.             xnum = p(6).*z;
  88.             xden = z;
  89.             for i = 1:4
  90.                xnum = (xnum + p(i)) .* z;
  91.                xden = (xden + q(i)) .* z;
  92.             end
  93.             result(k) = z .* (xnum + p(5)) ./ (xden + q(5));
  94.             result(k) = (1/sqrt(pi) -  result(k)) ./ y;
  95.             if jint ~= 2
  96.                z = fix(y*16)/16;
  97.                del = (y-z).*(y+z);
  98.                result(k) = exp(-z.*z) .* exp(-del) .* result(k);
  99.                k = find(~finite(result));
  100.                result(k) = 0*k;
  101.             end
  102.     end
  103. %
  104. %   fix up for negative argument, erf, etc.
  105. %
  106.     if jint == 0
  107.             k = find(x > xbreak);
  108.             result(k) = (0.5 - result(k)) + 0.5;
  109.             k = find(x < -xbreak);
  110.             result(k) = (-0.5 + result(k)) - 0.5;
  111.     elseif jint == 1
  112.             k = find(x < -xbreak);
  113.             result(k) = 2. - result(k);
  114.     else  % jint must = 2 
  115.             k = find(x < -xbreak);
  116.             z = fix(x(k)*16)/16;
  117.             del = (x(k)-z).*(x(k)+z);
  118.             y = exp(z.*z) .* exp(del);
  119.             result(k) = (y+y) - result(k);
  120.     end
  121.