1      DOUBLE PRECISION FUNCTION DERF(X)
2C***BEGIN PROLOGUE  DERF
3C***DATE WRITTEN   770701   (YYMMDD)
4C***REVISION DATE  840425   (YYMMDD)
5C***REVISION HISTORY (YYMMDD)
6C   000601   Added EXTERNAL statement for DERFC
7C***CATEGORY NO.  C8A,L5A1E
8C***KEYWORDS  DOUBLE PRECISION,ERROR FUNCTION,SPECIAL FUNCTION
9C***AUTHOR  FULLERTON, W., (LANL)
10C***PURPOSE  Computes the d.p. error function, ERF, of X.
11C***DESCRIPTION
12C
13C DERF(X) calculates the double precision error function for double
14C precision argument X.
15C
16C Series for ERF        on the interval  0.          to  1.00000E+00
17C                                        with weighted error   1.28E-32
18C                                         log weighted error  31.89
19C                               significant figures required  31.05
20C                                    decimal places required  32.55
21C***REFERENCES  (NONE)
22C***ROUTINES CALLED  D1MACH,DCSEVL,DERFC,INITDS
23C***END PROLOGUE  DERF
24      DOUBLE PRECISION X, ERFCS(21), SQEPS, SQRTPI, XBIG, Y, D1MACH,
25     1  DCSEVL, DERFC
26      EXTERNAL DERFC
27      DATA ERF CS(  1) / -.4904612123 4691808039 9845440333 76 D-1     /
28      DATA ERF CS(  2) / -.1422612051 0371364237 8247418996 31 D+0     /
29      DATA ERF CS(  3) / +.1003558218 7599795575 7546767129 33 D-1     /
30      DATA ERF CS(  4) / -.5768764699 7674847650 8270255091 67 D-3     /
31      DATA ERF CS(  5) / +.2741993125 2196061034 4221607914 71 D-4     /
32      DATA ERF CS(  6) / -.1104317550 7344507604 1353812959 05 D-5     /
33      DATA ERF CS(  7) / +.3848875542 0345036949 9613114981 74 D-7     /
34      DATA ERF CS(  8) / -.1180858253 3875466969 6317518015 81 D-8     /
35      DATA ERF CS(  9) / +.3233421582 6050909646 4029309533 54 D-10    /
36      DATA ERF CS( 10) / -.7991015947 0045487581 6073747085 95 D-12    /
37      DATA ERF CS( 11) / +.1799072511 3961455611 9672454866 34 D-13    /
38      DATA ERF CS( 12) / -.3718635487 8186926382 3168282094 93 D-15    /
39      DATA ERF CS( 13) / +.7103599003 7142529711 6899083946 66 D-17    /
40      DATA ERF CS( 14) / -.1261245511 9155225832 4954248533 33 D-18    /
41      DATA ERF CS( 15) / +.2091640694 1769294369 1705002666 66 D-20    /
42      DATA ERF CS( 16) / -.3253973102 9314072982 3641600000 00 D-22    /
43      DATA ERF CS( 17) / +.4766867209 7976748332 3733333333 33 D-24    /
44      DATA ERF CS( 18) / -.6598012078 2851343155 1999999999 99 D-26    /
45      DATA ERF CS( 19) / +.8655011469 9637626197 3333333333 33 D-28    /
46      DATA ERF CS( 20) / -.1078892517 7498064213 3333333333 33 D-29    /
47      DATA ERF CS( 21) / +.1281188399 3017002666 6666666666 66 D-31    /
48      DATA SQRTPI / 1.772453850 9055160272 9816748334 115D0 /
49      DATA NTERF, XBIG, SQEPS / 0, 2*0.D0 /
50C***FIRST EXECUTABLE STATEMENT  DERF
51      IF (NTERF.NE.0) GO TO 10
52      NTERF = INITDS (ERFCS, 21, 0.1*SNGL(D1MACH(3)))
53      XBIG = DSQRT (-DLOG(SQRTPI*D1MACH(3)))
54      SQEPS = DSQRT (2.0D0*D1MACH(3))
55C
56 10   Y = DABS(X)
57      IF (Y.GT.1.D0) GO TO 20
58C
59C ERF(X) = 1.0 - ERFC(X)  FOR  -1.0 .LE. X .LE. 1.0
60C
61      IF (Y.LE.SQEPS) DERF = 2.0D0*X/SQRTPI
62      IF (Y.GT.SQEPS) DERF = X*(1.0D0 + DCSEVL (2.D0*X*X-1.D0,
63     1  ERFCS, NTERF))
64      RETURN
65C
66C ERF(X) = 1.0 - ERFC(X) FOR ABS(X) .GT. 1.0
67C
68 20   IF (Y.LE.XBIG) DERF = DSIGN (1.0D0-DERFC(Y), X)
69      IF (Y.GT.XBIG) DERF = DSIGN (1.0D0, X)
70C
71      RETURN
72      END
73