SUBROUTINE ERRINT (X,ERF,ERFC)
C                (ASA FORTRAN)
C
C AUTOMATIC COMPUTING METHODS FOR SPECIAL FUNCTIONS
C    PART 1. ERROR, PROBABILITY, AND RELATED FUNCTIONS
C BY IRENE A. STEGUN AND RUTH ZUCKER
C J. RES. NBS, VOL 74B, NO. 3, 211-224, (1970)
C
C USAGE.     CALL ERRINT (X,ERF,ERFC)
C            ENTERING WITH X(A REAL VARIABLE), ONE OBTAINS BOTH ERF X 
C            AND ERFC X TO THE MAXIMUM MACHINE ACCURACY POSSIBLE.
C
C            THE PROGRAM, AS SET UP BELOW
C            (1) WITH THE DOUBLE PRECISION TYPE STATEMENT,
C            (2) NBC=11, NBM=60 ON THE FIRST DATA CARD AND
C            (3) THE REMAINING CONSTANTS ON THE SECOND AND THIRD DATA 
C            CARDS GIVEN AS DOUBLE PRECISION CONSTANTS, YIELDS DOUBLE 
C            PRECISION RESULTS ON THE UNIVAC 1108 SYSTEM.
C
C            SINGLE PRECISION RESULTS ARE OBTAINED BY
C            (1) DELETING THE TYPE STATEMENT
C            (2) SETTING NBC=8, NBM=27 ON THE FIRST DATA CARD
C            (3) CHANGING THE D'S TO E'S ON THE SECOND AND THIRD DATA 
C            CARDS. 
C            (4) CHANGING THE FUNCTION TYPE - DABS TO ABS, DEXP TO EXP
C
C            FOR OTHER VALUES OF NBM, THE CONSTANT ON THE THIRD DATA
C            CARD SHOULD BE GIVEN TO THE CORRESPONDING PRECISION.
C
C      CAUTION.    THE SUBROUTINE CANNOT READILY BE ADAPTED TO COMPUTE
C                  THE ERROR FUNCTION FOR A COMPLEX ARGUMENT.
C
C   ERF(X)  =(2/SQRT(PI))*INTEGRAL FROM 0 TO X OF EXP(-T**2) DT
C   ERFC(X) =(2/SQRT(PI))*INTEGRAL FROM X TO INFINITY OF EXP(-T**2) DT
C           =1-ERF(X)
C
C   ERF(-X)  = -ERF(X)
C   ERFC(-X) = 2-ERFC(X)
C
C   NBC=NO. OF BINARY DIGITS IN THE CHARACTERISTIC OF A FLOATING POINT
C       NUMBER.
C   NBM=ACCURACY DESIRED OR MAX. NO. OF BINARY DIGITS IN THE MANTISSA OF
C       A FLOATING POINT NUMBER.
C
C METHOD.    POWER SERIES FOR ABS(X) .LE. 1,ULPS,UPPER LIMIT FOR SERIES
C            CONTINUED FRACTION FOR ABS(X) .GT. 1 AND .LE. ULCF (UPPER
C                              LIMIT OF CONTINUED FRACTION) 
C
C RANGE.     ABS(X) .LE. ULCF           ULCF=.83*(2.**((NBC-1)/2))
C            ABS(ULCF) APPROX.  9.3, NBC=8
C                              26.5, NBC=11
C            BEYOND THIS RANGE THE LIMITING VALUES OF THE FUNCTIONS ARE
C            SUPPLIED - ERF(INF)=1.0,   ERFC(INF)=0.
C
C ACCURACY.  ROUTINE WILL RETURN (NBM-I-3) SIGNIFICANT BINARY DIGITS
C            WHERE I IS THE NUMBER OF BINARY DIGITS REPRESENTING THE
C            INTEGER PART OF X**2. (THIS IS ESSENTIALLY THE ACCURACY OF
C            THE EXPONENTIAL ROUTINE.)
C
C PRECISION. VARIABLE - BY SETTING DESIRED NBM.
C
C MAXIMUM TIME.            UNIVAC 1108-EXEC 11
C   (SECONDS)                NBM=27    NBM=60
C                             .003      .015
C
C
C
C
      DOUBLE PRECISION AN,BN,CONS,C1,DN,ERF,ERFC,F,FN,FNM1,FNM2,FOUR, 
     1   GN,GNM1,GNM2,ONE,PREV,PWR,RNBC,SCF,SUM,TN,TOLER,TRRTPI,TWO,
     2   ULCF,ULPS,WN,X,Y,YSQ 
C
      DATA NBC,NBM/ 11,60/
      DATA ONE,TWO,FOUR,ULPS,CONS / 1.D0,2.D0,4.D0,1.D0,.83D0 /
      DATA TRRTPI/1.128379167095512574D0 /
C
      RNBC=NBC
      TOLER = TWO ** (-NBM)
C
C             TEST ON ZERO
C
      IF (X) 2,1,2
C
  1   ERF=0
      ERFC=ONE
C
      RETURN
C
  2   Y=DABS(X)
      YSQ=Y ** 2
C
      IF (Y-ULPS) 3,3,4
C
C             MAXIMUM ARGUMENT
C
  4   C1=TWO ** ((RNBC-ONE)/TWO)
      ULCF=CONS * C1
C
C             SCALE FACTOR
C
      SCF=TWO ** (C1 ** 2 - RNBC)
C
C             LIMITING VALUE
C
      IF (Y-ULCF) 10,10,11
C
  11  ERF=ONE
      ERFC=0
      GO TO 7
C
C         METHOD --- POWER SERIES
C
  3   SUM=0
      DN=ONE
      TN=ONE
      PWR=TWO*YSQ
  6   DN=DN + TWO
      TN=PWR*TN/DN
      SUM=TN+SUM
C
C             TOLERANCE CHECK 
C
      IF (TN-TOLER) 5,6,6
C
  5   ERF=(SUM+ONE)*TRRTPI*Y*DEXP(-YSQ) 
      ERFC=ONE - ERF
C
C             NEGATIVE ARGUMENT
C
  7   IF (X) 8,9,9
C
  8   ERF=-ERF
      ERFC=TWO - ERFC
C
  9   RETURN
C
C         METHOD --- CONTINUED FRACTION 
C
  10  FNM2=0
      GNM2=ONE
      FNM1=TWO * Y
      GNM1=TWO * YSQ + ONE
C
      PREV=FNM1/GNM1
C
      WN=ONE
      BN=GNM1 + FOUR
C
  14  AN=-WN * (WN + ONE)
C
      FN=BN * FNM1 + AN * FNM2
      GN=BN * GNM1 + AN * GNM2
C
      F=FN/GN
C
C             TOLERANCE CHECK 
C
      IF (DABS( ONE - (F/PREV) ) - TOLER) 12,12,13
  13  IF (PREV - F) 17,17,18
C
C       BOTH FN AND GN MUST BE TESTED IF ABS (X) .LT. .61
  17  IF (GN .LT. SCF) GO TO 16
C
C             SCALING
C
      FN=FN/SCF
      GN=GN/SCF
      FNM1=FNM1/SCF 
      GNM1=GNM1/SCF 
C
  16  FNM2=FNM1
      GNM2=GNM1
      FNM1=FN
      GNM1=GN
C
      WN=WN + TWO
      BN=BN + FOUR
      PREV=F
      GO TO 14
C
  18  F=PREV
  12  ERFC=F*DEXP(-YSQ)*TRRTPI/TWO
      ERF=ONE - ERFC
C
      GO TO 7
C
      END