Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50371/gauss.for
There are 3 other files named gauss.for in the archive. Click here to see a list.
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