/***** A BETTER COMPLEMENTARY ERROR FUNCTION

Author:
	Dr. Allen Wilkinson           (phone) (216) 433-2075
	NASA Lewis Research Center    (fax)   (216) 433-8660
	M.S. 500-102
	21000 Brookpark Road
	Cleveland, OH 44135  USA      (INTERNET) AW@LERC.NASA.GOV
Date: 1 Sep 95
Provided without guarantee for public non-commercial use.

I have found that the native GAUSS erfc function does not do well for x>5
when multiplying exp(x).*erfc(x) in the limit of large x, where the product
should vanish.

I have found a better numerical method for all (x):
*****/

/* Procedure using Numerical Recipes in C ... book 2nd Ed. page 221
   Returns the complementary error function erfc(x) with fractional error
   everywhere less than 1.2e-7 "based on Chebyshev fitting to an inspired
   guess as to the functional form"
*/

proc (1) = nrcerfc(x);

  local z,t,y;

  z = abs(x);
  t = 1/(1 + 0.5*z);

  y = t.*exp(-z.*z - 1.26551223 + t.*(1.00002368 + t.*(0.37409196 +
      t.*(0.09678418 + t.*(-0.18628806 + t.*(0.27886807 + t.*(-1.13520398 +
      t.*(1.48851587 + t.*(-0.82215223 + t.*0.17087277)))))))));

  retp(y);

  @ ? y~erfc(x)~(exp(-x^2)./sqrt(2)) @
  @ The last printed expression is the assymptotic form of erfc() for x > 10 @
endp;