C***************************************************************
C
C	THIS SUBPROGRAM CALCULATES THE DOUBLE PRECISION
C	COMPLEMENTARY ERROR FUNCTION OF A DOUBLE PRECISION
C	ARGUMENT.
C
C	W.G.MADISON	75-10-23
C
      DOUBLE PRECISION FUNCTION DERFC(DX)
      IMPLICIT DOUBLE PRECISION (D)
      DATA DSQR2/1.4142 13562 37309 50488 D0/
      DERFC = 2.D0 * DUNCDL(DSQR2 * -DX)
      RETURN
      END
 
C***************************************************************
C
C	THIS SUBPROGRAM CALCULATES THE SINGLE PRECISION
C	COMPLEMENTARY ERROR FUNCTION OF A SINGLE PRECISION
C	ARGUMENT.
C
C	W.G.MADISON	75-10-23
C
      FUNCTION ERFC(X)
      DATA DSQR2/1.4142 13562 37309 50488 E0/
      ERFC = 2.E0 * UNCDL(DSQR2 * -X)
      RETURN
      END
 
C***************************************************************
C
C	THIS SUBPROGRAM CALCULATES THE SINGLE PRECISION
C	COMPLEMENTARY UNIT NORMAL PROBABILITY FUNCTION OF A
C	SINGLE PRECISION ARGUMENT.
C
C	W.G.MADISON	75-10-23
C
      FUNCTION UNCDR(X)
      UNCDR = UNCDL(-X)
      RETURN
      END
 
C***************************************************************
C
C	THIS SUBPROGRAM CALCULATES THE DOUBLE PRECISION
C	COMPLEMENTARY UNIT NORMAL PROBABILITY FUNCTION OF A
C	DOUBLE PRECISION ARGUMENT.
C
C	W.G.MADISON	75-10-23
C
      DOUBLE PRECISION FUNCTION DUNCDR(DX)
      IMPLICIT DOUBLE PRECISION (D)
      DUNCDR = DUNCDL(-DX)
      RETURN
      END
 
C***************************************************************
C
C	THIS SUBPROGRAM CALCULATES THE DOUBLE PRECISION ERROR
C	FUNCTION OF A DOUBLE PRECISION ARGUMENT.
C
C	W.G.MADISON	75-10-22
C
      DOUBLE PRECISION FUNCTION DERF(DX)
      IMPLICIT DOUBLE PRECISION (D)
      DATA DSQR2/1.4142 13562 37309 50488 D0/
      DERF=2.D0 * DUNCDL(DSQR2 * DX) - 1.D0
      RETURN
      END
 
C***************************************************************
C
C	THIS SUBPROGRAM CALCULATES THE SINGLE PRECISION ERROR
C	FUNCTION OF A SINGLE PRECISION ARGUMENT.
C
C	W.G.MADISON	75-10-22
C
      REAL FUNCTION ERF(X)
      DATA DSQR2/1.4142 13562 37309 50488 E0/
      ERF=2.E0 * UNCDL(DSQR2 * DX) - 1.E0
      RETURN
      END
 
C***************************************************************
C
C	THIS SUBPROGRAM CALCULATES THE INTEGRAL OF THE NORMAL
C	CURVE BETWEEN THE LIMITS OF MINUS INF. AND THE ARG.
C	BOTH THE ARGUMENT AND THE RETURN VALUE ARE DOUBLE
C	PRECISION. THIS SUBPROGRAM IS A DIRECT IMPLEMENTATION
C	OF CACM ALGORITHM NUMBER 304, RESTRICTED TO THE LEFT
C	TAIL ONLY.
C
C	W.G.MADISON	75-10-22
C
C	MODIFIED TO INCLUDE THE SUGGESTIONS OF HOLMGREN AND
C	OF ADAMS, AND ALSO TO PROVIDE RANGE CHECKING AT THE
C	INPUT ARGUMENT IN ORDER TO PREVENT UNDER/OVERFLO.
C
C	W.G.MADISON	75-11-06
C
      DOUBLE PRECISION FUNCTION DUNCDL(DXX)
      IMPLICIT DOUBLE PRECISION (D)
      LOGICAL UPPER
      EQUIVALENCE (DXY, RXY), (DX, RX), (DQ2, RQ2)
      DATA DPI/0.39894 22804 01432 67793 99461 D0/
      DATA RUPRX/8.8261 E0/,RLOWRX/-13.1 E0/
C	DPI IS ONE DIVIDED BY (SQUARE ROOT OF TWO*PI)
C
      DX = DXX
      DXY = DXX
C
C     IF (RX .GT. RUPRX)
      IF (RX .GT. RUPRX) GOTO 99999
      GOTO 99998
99999 CONTINUE
        DUNCDL = 1.0 D0
C     ELSE
      GOTO 99997
99998 CONTINUE
C       IF (RX .LT. RLOWRX)
        IF (RX .LT. RLOWRX) GOTO 99996
      GOTO 99995
99996 CONTINUE
          DUNCDL = 0.0 D0
C       ELSE
      GOTO 99994
99995 CONTINUE
C         IF (RX .EQ. 0.0)
          IF (RX .EQ. 0.0) GOTO 99993
      GOTO 99992
99993 CONTINUE
            DUNCDL = 0.5D0
C         ELSE
      GOTO 99991
99992 CONTINUE
            UPPER = (RX .LE. 0.0)
            DX = DABS(DX)
            DX2 = DX * DX
            DY = DPI * DEXP(-0.5D0 * DX2)
            DN = DY / DX
C           IF (.NOT. UPPER .AND. ((1.D0-DN) .EQ. 1.D0))
            IF (.NOT. UPPER .AND. ((1.D0-DN) .EQ. 1.D0)) GOTO 99990
      GOTO 99989
99990 CONTINUE
              DUNCDL = 1.D0
C           ELSE
      GOTO 99988
99989 CONTINUE
C             IF ((RXY .LT. -2.32 E0) .OR. (RXY .GT. 3.5 E0))
              IF ((RXY .LT. -2.32 E0) .OR. (RXY .GT. 3.5 E0)) GOTO 99987
      GOTO 99986
99987 CONTINUE
                DA1 = 2.0 D0
                DA2 = 0.0 D0
                DN = DX2 + 3.0 D0
                DP1 = DY
                DQ1 = DX
                DP2 = (DN - 1.0 D0) * DY
                DQ2 = DN * DX
                DM = DP1 / DQ1
                DT = DP2 / DQ2
C               IF (.NOT. UPPER)
                IF (.NOT. UPPER) GOTO 99985
      GOTO 99984
99985 CONTINUE
                  DM = 1.0 D0 - DM
                  DT = 1.0 D0 - DT
C               FI
99984 CONTINUE
C               DO UNTIL, ((DM .EQ. DT) .OR. (DS .EQ. DT))
      GOTO 99983
99982 CONTINUE
      IF((DM .EQ. DT) .OR. (DS .EQ. DT)) GOTO 99981
99983 CONTINUE
                  DN = DN + 4.0 D0
                  DA1 = DA1 - 8.0 D0
                  DA2 = DA1 + DA2
                  DS = DA2 * DP1 + DN * DP2
                  DP1 = DP2
                  DP2 = DS
                  DS = DA2 * DQ1 + DN * DQ2
                  DQ1 = DQ2
                  DQ2 = DS
                  DS = DM
                  DM = DT
C                 IF (RQ2 .GT. 1.0 E30)
                  IF (RQ2 .GT. 1.0 E30) GOTO 99980
      GOTO 99979
99980 CONTINUE
                    DP1 = DP1 * 1.D-30
                    DP2 = DP2 * 1.D-30
                    DQ1 = DQ1 * 1.D-30
                    DQ2 = DQ2 * 1.D-30
C                 FI
99979 CONTINUE
C                 IF (UPPER)
                  IF (UPPER) GOTO 99978
      GOTO 99977
99978 CONTINUE
                    DT = DP2 / DQ2
C                 ELSE
      GOTO 99976
99977 CONTINUE
                    DT = 1.D0 - DP2 / DQ2
C                 FI
99976 CONTINUE
C               OD
      GOTO 99982
99981 CONTINUE
                DUNCDL = DT
C             ELSE
      GOTO 99975
99986 CONTINUE
                DX = DY * DX
                DS = DX
                DN = 1.D0
                DT = 0.D0
C               DO WHILE, (DS .NE. DT)
99974 CONTINUE
      IF(DS .NE. DT) GOTO 99973
      GOTO 99972
99973 CONTINUE
                  DN = DN + 2.D0
                  DT = DS
                  DX = DX * DX2 / DN
                  DS = DS + DX
C               OD
      GOTO 99974
99972 CONTINUE
C               IF (UPPER)
                IF (UPPER) GOTO 99971
      GOTO 99970
99971 CONTINUE
                  DUNCDL = 0.5D0 - DS
C               ELSE
      GOTO 99969
99970 CONTINUE
                  DUNCDL = 0.5D0 + DS
C               FI
99969 CONTINUE
C             FI
99975 CONTINUE
C           FI
99988 CONTINUE
C         FI
99991 CONTINUE
C       FI
99994 CONTINUE
C     FI
99997 CONTINUE
      RETURN
      END
CC     **********
CC
CC                  SINGLE-PRECISION SUBPROGRAM
CC
C      SUBROUTINE UNRCUQ(IPROB,U,DENS,CPRB)
CC	RETURNS UNIT-NORMAL DENSITY, ALSO LEFT-TAIL CUMULATIVE
CC	PROBABILITY UNLESS CALLED WITH IPROB = 0.  NBS HNDBK 26.2.17.
CC	ABSOLUTE ERROR < .0000001 FOR ALL U.
CC	RELATIVE ERROR < .0001	  FOR U > -3.2.
CC	RELATIVE ERROR < .1	  FOR U > -25.7
C
C	THE ABOVE REFLECTS THE MANECON VERSION AS IMPLEMENTED
C	BY SCHLAIFER/SCHLEIFER.
C
C	MODIFIED FOR GENERAL USE BY
C	W.G.MADISON	75-10-28
C
      REAL FUNCTION UNCDL(U)
C
C
C
C
      DATA A/.31938153/, B/-1.1164195437/, C/-4.9962391778/
      DATA D/-1.0223286745/, E/-.7304159575/, P/.2316419/
      DATA ARGBND/13.2/
C       GREATEST X SUCH THAT EXP(-.5*X*X) IS WITHIN PDP-10 CAPACITY.
C
C      DENS = 0.
      CPRB = 0.
      IF(ABS(U).GT.ARGBND) GOTO 30
      DENS = .3989422804*EXP(-.5*U*U)
C
      T = 1./(1.+P*ABS(U))
      S = A*T*(1.+B*T*(1.+C*T*(1.+D*T*(1.+E*T))))
      CPRB = S*DENS
   30 IF(U.GT.0.) CPRB = 1.-CPRB
      UNCDL = CPRB
      RETURN
      END

