Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-03 - decus/20-0095/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