Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50371/gauss.sf4
There is 1 other file named gauss.sf4 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
IF (RX .GT. RUPRX)
DUNCDL = 1.0 D0
ELSE
IF (RX .LT. RLOWRX)
DUNCDL = 0.0 D0
ELSE
IF (RX .EQ. 0.0)
DUNCDL = 0.5D0
ELSE
UPPER = (RX .LE. 0.0)
DX = DABS(DX)
DX2 = DX * DX
DY = DPI * DEXP(-0.5D0 * DX2)
DN = DY / DX
IF (.NOT. UPPER .AND. ((1.D0-DN) .EQ. 1.D0))
DUNCDL = 1.D0
ELSE
IF ((RXY .LT. -2.32 E0) .OR. (RXY .GT. 3.5 E0))
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
IF (.NOT. UPPER)
DM = 1.0 D0 - DM
DT = 1.0 D0 - DT
FI
DO UNTIL, ((DM .EQ. DT) .OR. (DS .EQ. DT))
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
IF (RQ2 .GT. 1.0 E30)
DP1 = DP1 * 1.D-30
DP2 = DP2 * 1.D-30
DQ1 = DQ1 * 1.D-30
DQ2 = DQ2 * 1.D-30
FI
IF (UPPER)
DT = DP2 / DQ2
ELSE
DT = 1.D0 - DP2 / DQ2
FI
OD
DUNCDL = DT
ELSE
DX = DY * DX
DS = DX
DN = 1.D0
DT = 0.D0
DO WHILE, (DS .NE. DT)
DN = DN + 2.D0
DT = DS
DX = DX * DX2 / DN
DS = DS + DX
OD
IF (UPPER)
DUNCDL = 0.5D0 - DS
ELSE
DUNCDL = 0.5D0 + DS
FI
FI
FI
FI
FI
FI
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