module Erf(erf, erfc, erfcx, normcdf) where -- Transcribed from FORTRAN by Lennart Augustsson, lennart@augustsson.net -------------------------------------------------------------------- -- -- This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) -- for a real argument x. It contains three FUNCTION type -- subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), -- and one SUBROUTINE type subprogram, CALERF. The calling -- statements for the primary entries are: -- -- Y=ERF(X) (or Y=DERF(X)), -- -- Y=ERFC(X) (or Y=DERFC(X)), -- and -- Y=ERFCX(X) (or Y=DERFCX(X)). -- -- The routine CALERF is intended for internal packet use only, -- all computations within the packet being concentrated in this -- routine. The function subprograms invoke CALERF with the -- statement -- -- CALL CALERF(ARG,RESULT,JINT) -- -- where the parameter usage is as follows -- -- Function Parameters for CALERF -- call ARG Result JINT -- -- ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 -- ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 -- ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 -- -- The main computation evaluates near-minimax approximations -- from "Rational Chebyshev approximations for the error function" -- by W. J. Cody, Math. Comp., 1969, PP. 631-638. This -- transportable program uses rational functions that theoretically -- approximate erf(x) and erfc(x) to at least 18 significant -- decimal digits. The accuracy achieved depends on the arithmetic -- system, the compiler, the intrinsic functions, and proper -- selection of the machine-dependent constants. -- --------------------------------------------------------------------- --------------------------------------------------------------------- -- -- Explanation of machine-dependent constants -- -- XMIN = the smallest positive floating-point number. -- XINF = the largest positive finite floating-point number. -- XNEG = the largest negative argument acceptable to ERFCX; -- the negative of the solution to the equation -- 2*exp(x*x) = XINF. -- XSMALL = argument below which erf(x) may be represented by -- 2*x/sqrt(pi) and above which x*x will not underflow. -- A conservative value is the largest machine number X -- such that 1.0 + X = 1.0 to machine precision. -- XBIG = largest argument acceptable to ERFC; solution to -- the equation: W(x) * (1-0.5/x**2) = XMIN, where -- W(x) = exp(-x*x)/[x*sqrt(pi)]. -- XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to -- machine precision. A conservative value is -- 1/[2*sqrt(XSMALL)] -- XMAX = largest acceptable argument to ERFCX; the minimum -- of XINF and 1/[sqrt(pi)*XMIN]. -- -- Approximate values for some important machines are: -- -- XMIN XINF XNEG XSMALL -- -- CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 -- CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 -- IEEE (IBM/XT, -- SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 -- IEEE (IBM/XT, -- SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 -- IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 -- UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 -- VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 -- VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 -- -- -- XBIG XHUGE XMAX -- -- CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 -- CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 -- IEEE (IBM/XT, -- SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 -- IEEE (IBM/XT, -- SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 -- IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 -- UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 -- VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 -- VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 -- --------------------------------------------------------------------- --------------------------------------------------------------------- -- -- Error returns -- -- The program returns ERFC = 0 for ARG >= XBIG; -- -- ERFCX = XINF for ARG < XNEG; -- and -- ERFCX = 0 for ARG >= XMAX. -- -- -- Intrinsic functions required are: -- -- abs, truncate, exp -- -- -- Author: W. J. Cody -- Mathematics and Computer Science Division -- Argonne National Laboratory -- Argonne, IL 60439 -- -- Latest modification: March 19, 1990 -- -------------------------------------------------------------------- -- Horner's scheme for polynomial evaluation horner :: (Num a) => [a] -> a -> a horner ks x = foldl (\ r k -> r*x + k) 0 ks -- Like FORTRAN's AINT, truncate an FP number aint :: (RealFrac a) => a -> a aint = fromInteger . truncate data Mode = Erf | Erfc | Erfcx deriving (Eq) {-# SPECIALIZE calerf :: Mode -> Double -> Double #-} calerf :: (RealFrac a, Floating a) => Mode -> a -> a calerf mode x = {- private static final double X_MIN = Double.MIN_VALUE; private static final double X_INF = Double.MAX_VALUE; private static final double X_NEG = -9.38241396824444; private static final double X_SMALL = 1.110223024625156663E-16; private static final double X_BIG = 9.194E0; private static final double X_HUGE = 1.0D / (2.0D * Math.sqrt(X_SMALL)); private static final double X_MAX = Math.min(X_INF, (1 / (Math.sqrt(Math.PI) * X_MIN))); -} let -- FP format dependent constants xinf = 1.79e308 xneg = -26.628 xsmall = 1.11e-16 xbig = 26.543 xhuge = 6.71e7 xmax = 2.53e307 in let y = abs x in if y <= 0.46875 then -------------------------------------------------------------------- -- Evaluate erf for |X| <= 0.46875 -------------------------------------------------------------------- let a = [1.85777706184603153e-1, 3.16112374387056560e00, 1.13864154151050156e02, 3.77485237685302021e02, 3.20937758913846947e03 ] b = [1.0, 2.36012909523441209e01, 2.44024637934444173e02, 1.28261652607737228e03, 2.84423683343917062e03] ysq = if y > xsmall then y * y else 0.0 xnum = horner a ysq xden = horner b ysq result = x * xnum / xden result' = if mode /= Erf then 1.0 - result else result in if mode == Erfcx then exp(ysq) * result' else result' else let reslt = -------------------------------------------------------------------- -- Evaluate erfcx for 0.46875 <= |X| <= 4.0 -------------------------------------------------------------------- if y <= 4.0 then let c = [2.15311535474403846e-8, 5.64188496988670089e-1, 8.88314979438837594e0, 6.61191906371416295e01,2.98635138197400131e02, 8.81952221241769090e02,1.71204761263407058e03, 2.05107837782607147e03,1.23033935479799725e03 ] d = [1.0, 1.57449261107098347e01,1.17693950891312499e02, 5.37181101862009858e02,1.62138957456669019e03, 3.29079923573345963e03,4.36261909014324716e03, 3.43936767414372164e03,1.23033935480374942e03] xnum = horner c y xden = horner d y result = xnum / xden in if mode /= Erfcx then let ysq = aint(y*16.0)/16.0 del = (y-ysq)*(y+ysq) in exp(-ysq*ysq) * exp(-del) * result else result -------------------------------------------------------------------- -- Evaluate erfcx for |X| > 4.0 -------------------------------------------------------------------- else let sqrpi = 0.56418958354775628695 in if y >= xbig && (mode /= Erfcx || y >= xmax) then 0.0 else if y >= xbig && y >= xhuge then sqrpi / y else let p = [1.63153871373020978e-2, 3.05326634961232344e-1,3.60344899949804439e-1, 1.25781726111229246e-1,1.60837851487422766e-2, 6.58749161529837803e-4] q = [1.0, 2.56852019228982242e00,1.87295284992346047e00, 5.27905102951428412e-1,6.05183413124413191e-2, 2.33520497626869185e-3] ysq = 1.0 / (y * y) xnum = horner p ysq xden = horner q ysq result = ysq * xnum / xden result' = (sqrpi - result) / y in if mode /= Erfcx then let x' = aint(y*16.0)/16.0 del = (y-x')*(y+x') in exp(-x'*x') * exp(-del) * result' else result' -------------------------------------------------------------------- -- Fix up for negative argument, erf, etc. -------------------------------------------------------------------- in if mode == Erf then let result = (0.5 - reslt) + 0.5 in if x < 0.0 then -result else result else if mode == Erfc then if x < 0.0 then 2.0 - reslt else reslt else if x < 0.0 then if x < xneg then xinf else let x' = aint(x*16.0)/16.0 del = (x-x')*(x+x') y' = exp(x'*x') * exp(del) in (y'+y') - reslt else reslt -- Computes approximate values for erf(x) = 2/sqrt(pi) * integral(0,x,\t.exp(-t*t)) {-# SPECIALIZE erf :: Double -> Double #-} erf :: (RealFrac a, Floating a) => a -> a erf x = calerf Erf x -- Computes approximate values for erfc(x) = 1 - erf(x). {-# SPECIALIZE erfc :: Double -> Double #-} erfc :: (RealFrac a, Floating a) => a -> a erfc x = calerf Erfc x -- Computes approximate values for erfcx(x) = exp(x*x) * erfc(x). {-# SPECIALIZE erfcx :: Double -> Double #-} erfcx :: (RealFrac a, Floating a) => a -> a erfcx x = calerf Erfcx x -- Computes approximate values for the cumulative normal distribution function. {-# SPECIALIZE normcdf :: Double -> Double #-} normcdf :: (RealFrac a, Floating a) => a -> a normcdf x = 0.5 * erfc(-x / sqrt 2)