[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

SUBROUTINE X,NMAX,ACC,FZERO,

Found at: ftp.icm.edu.pl:70/packages/netlib/toms-2014-06-10/521

      SUBROUTINE INERFC(X, NMAX, ACC, FZERO, F, IFLAG)                  INE   10
C THIS PROGRAM GENERATES FN(X) FOR N=1,2,...,NMAX, WHERE
C FN(X)=EXP(X**2)*INERFC(X), IF X.GE.0., AND FN(X)=INERFC(X), IF
C X.LT.0., INERFC(X) BEING THE N-TH REPEATED INTEGRAL OF THE COERROR
C FUNCTION ERFC(X), EXTENDED FROM X TO INFINITY. THE PROGRAM ALSO
C SUPPLIES EXP(X**2)*ERFC(X),IF X.GE.0., AND ERFC(X), IF X.LT.0..
C THE USER HAS TO SPECIFY THE DESIRED ACCURACY, AS WELL AS THE PRECISION
C OF THE HIGHEST PRECISION MODE HE IS PREPARED TO USE. THE FORMER IS
C INPUT VIA THE PARAMETER ACC IN THE CALL LIST, WHICH DESIGNATES THE
C NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS DESIRED IN THE RESULTS.
C THE LATTER IS INPUT VIA THE PARAMETER PREC IN THE DATA STATEMENT. THE
C VALUE OF PREC SHOULD BE SET APPROXIMATELY EQUAL TO BETA*ALOG(2.)/
C ALOG(10.), WHERE BETA IS THE NUMBER OF BINARY DIGITS AVAILABLE IN THE
C MANTISSA OF THE HIGHEST PRECISION FLOATING-POINT WORD. ORDINARILY, THE
C HIGHEST PRECISION IS DOUBLE PRECISION, EITHER HARDWARE GENERATED OR
C SOFTWARE GENERATED, BUT IT MAY ALSO BE SINGLE PRECISION. IN ANY CASE,
C ACC HAS TO BE LESS THAN PREC, THE DIFFERENCE BEING AT LEAST EQUAL TO
C ONE, PREFERABLY SEVERAL UNITS LARGER. EVIDENTLY, ACC SHOULD NOT BE
C SPECIFIED TO EXCEED THE PRECISION OF THE LOWEST PRECISION MODE USED.
C   THE DATA STATEMENT CONTAINS ANOTHER MACHINE-DEPENDABLE PARAMETER,
C BOTEXP, WHICH IS, APPROXIMATELY, THE SMALLEST NEGATIVE NUMBER SUCH
C THAT 10.**BOTEXP IS STILL REPRESENTABLE ON THE COMPUTER IN SINGLE PRE-
C CISION FLOATING-POINT ARITHMETIC. THE TYPE OF BOTEXP IS REAL. IN THE
C DATA STATEMENT BELOW THE PARAMETERS PREC AND BOTEXP ARE SET TO CORRE-
C SPOND TO THE MACHINE CHARACTERISTICS OF THE CDC 6500 COMPUTER.
C   PARAMETER LIST
C       X - - THE ARGUMENT OF INERFC. TYPE REAL.
C    NMAX - - THE UPPER LIMIT ON N. TYPE INTEGER.
C     ACC - - THE NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS DESIRED
C             IN THE RESULTS. TYPE REAL.
C   FZERO - - AN OUTPUT VARIABLE RETURNING EXP(X**2)*ERFC(X), IF
C             X.GE.0., AND ERFC(X), IF X.LT.0. TYPE REAL.
C       F - - THE NAME OF A ONE-DIMENSIONAL ARRAY HOLDING THE RESULTS
C             FN(X), N=1,2,...,NMAX. THE ARRAY HAS DIMENSION NMAX.
C   IFLAG - - AN ERROR FLAG INDICATING A NUMBER OF FATAL ERROR CON-
C             DITIONS UPON EXECUTION. TYPE INTEGER. THE VALUES OF THIS
C             VARIABLE HAVE THE FOLLOWING MEANINGS.
C             IFLAG=0 - NO ERROR CONDITION.
C             IFLAG=1 - ILLEGAL NEGATIVE OR ZERO VALUE OF NMAX.
C             IFLAG=2 - EXCESSIVELY LARGE VALUE OF NMAX. REDIMENSION THE
C                       ARRAYS R0 AND R1 TO HAVE DIMENSION .GE. NMAX, IF
C                       INDEED NMAX MUST BE THAT LARGE.
C             IFLAG=3 - INAPPROPRIATE ACCURACY SPECIFICATION.
C             IFLAG=4 - UNUSUAL UNDERFLOW OF EXP(-X**2). TRY REDUCING
C                       THE ERROR TOLERANCE PREC-ACC, EITHER BY DE-
C                       CREASING PREC OR BY INCREASING ACC.
C             IFLAG=5 - TAYLOR SERIES FOR ERFC(X) DOES NOT CONVERGE
C                       WITHIN 200 TERMS FOR SOME UNKNOWN REASON.
C             IFLAG=6 - NU IN THE BACKWARD RECURRENCE ALGORITHM EXCEEDS
C                       5000, PROBABLY BECAUSE X IS A SMALL POSITIVE
C                       NUMBER, NMAX RELATIVELY LARGE, AND THE ERROR
C                       TOLERANCE SMALL. TRY INCREASING THE ERROR TOLER-
C                       ANCE PREC-ACC, EITHER BY INCREASING PREC OR BY
C                       DECREASING ACC.
C    THE SUBROUTINE BELOW IS WRITTEN WITH THE ASSUMPTION THAT PART OF
C THE COMPUTATION IS DONE IN DOUBLE PRECISION AND THE REST IN SINGLE
C PRECISION. TO OBTAIN A VERSION OF THIS SUBROUTINE, THAT OPERATES
C ENTIRELY IN SINGLE PRECISION, THE FOLLOWING CHANGES MUST BE MADE.
C   (1) DELETE THE DOUBLE PRECISION TYPE DECLARATION.
C   (2) REPLACE THE SECTION OF PROGRAM FROM STATEMENT LABELED 30
C (INCLUSIVE) TO STATEMENT LABELED 60 (INCLUSIVE) BY THE FOLLOWING
C PIECE OF PROGRAM.
C          30 DEPS=.5*10.**(-PREC)
C             TE=C*X
C             SUM=1.-TE
C             N=0
C             N1=1
C          40 N=N+1
C             N1=N1+2
C             IF(N.GT.200) GOTO 80
C             TE=-TE*XSQ/FLOAT(N)
C             TERM=TE/FLOAT(N1)
C             SUM=SUM-TERM
C             IF(ABS(TERM).GT.DEPS*ABS(SUM)) GOTO 40
C             F1=T*SUM
C    THE ARRAYS R0 AND R1 IN THE DIMENSION STATEMENT ARE LOCAL TO THE
C SUBROUTINE. THEY BOTH NEED TO HAVE DIMENSION LARGER OR EQUAL TO NMAX.
C WE DECLARED THEIR DIMENSION TO BE 500, ASSUMING THAT NMAX WILL NOT
C EXCEED 500. IN THE EVENT IT DOES, THE SUBROUTINE EXITS WITH THE ERROR
C FLAG AT 2. IF THE USER EXPECTS HIS VALUES OF NMAX TO BE CONSISTENTLY
C MUCH LESS THAN 500, HE CAN SAVE STORAGE BY LOWERING THE DIMENSION OF
C R0 AND R1 AND AT THE SAME TIME ADJUSTING THE ERROR EXIT STATEMENT (THE
C THIRD EXECUTABLE STATEMENT IN THE PROGRAM).
C    WITH THE EXCEPTION NOTED UNDER IFLAG=4, ANY VARIABLE THAT UNDER-
C FLOWS MAY BE SET EQUAL TO ZERO. NO PRECAUTION IS TAKEN AGAINST OVER-
C FLOW, WHICH MAY OCCUR IN THE DO-LOOP AFTER STATEMENT 15 IF X.LT.0.
C AND ABS(X) AND/OR NMAX IS LARGE.
C    REFERENCE - W. GAUTSCHI, EVALUATION OF THE REPEATED INTEGRALS OF
C THE COERROR FUNCTION, TRANS. ON MATH. SOFTWARE.
      DIMENSION F(NMAX), R0(500), R1(500)
      DOUBLE PRECISION DEPS, DC, DX, DXSQ, DFM1, DF0, DF1, DN, DN1,
     *  DTE, DTERM, DSUM
      LOGICAL FRSTTM
      DATA FRSTTM, PREC, BOTEXP /.TRUE.,28.8989,-293./
      IFLAG = 0
      IF (NMAX.LT.1) GO TO 230
      IF (NMAX.GT.500) GO TO 240
      TOL = PREC - ACC
      IF (TOL.LT.1.) GO TO 250
      I = 1
      EPS = .5*10.**(-ACC)
      XSQ = X*X
      IF (.NOT.FRSTTM) GO TO 10
      P = ATAN(1.)
      AL10 = ALOG(10.)
      C = SQRT(1./P)
      DC = DSQRT(1.D0/DATAN(1.D0))
      FRSTTM = .FALSE.
   10 S = 0.
      IF (-XSQ.GT.AL10*BOTEXP) S = EXP(-XSQ)
      B = .5*AL10*TOL + .25*ALOG(2.*P)
      B1 = .5*AL10*(TOL-1.)
      IF (XSQ.GT.B) GO TO 90
      IF (X.LT.0.) GO TO 20
      F0 = C
      IF (S.EQ.0.) GO TO 260
      T = 1./S
      SQ = SQRT(2.*FLOAT(NMAX))
      IF ((X.GE.1.12 .AND. XSQ+SQ*X.GT.B) .OR. (X.LT.1.12 .AND.
     *  SQ*X.GT.B1)) GO TO 100
      GO TO 30
   20 F0 = C*S
      T = 1.
C
C EVALUATION OF ERFC(X) BY TAYLOR SERIES IN DOUBLE PRECISION
C
   30 DEPS = .5D0*10.D0**(-PREC)
      DX = DBLE(X)
      DXSQ = DX*DX
      DTE = DC*DX
      DSUM = 1.D0 - DTE
      DN = 0.D0
      DN1 = 1.D0
   40 DN = DN + 1.D0
      DN1 = DN1 + 2.D0
      IF (DN.GT.2.D2) GO TO 270
      DTE = -DTE*DXSQ/DN
      DTERM = DTE/DN1
      DSUM = DSUM - DTERM
      IF (DABS(DTERM).GT.DEPS*DABS(DSUM)) GO TO 40
      IF (X.LT.0.) GO TO 60
C
C EVALUATION OF EXP(X**2)*INERFC(X),X.GE.0.,BY FORWARD RECURSION IN
C DOUBLE PRECISION
C
      DF0 = DC
      DF1 = DEXP(DXSQ)*DSUM
      FZERO = SNGL(DF1)
      DO 50 N=1,NMAX
        DFM1 = DF0
        DF0 = DF1
        DF1 = (-DX*DF0+.5D0*DFM1)/DBLE(FLOAT(N))
        F(N) = SNGL(DF1)
   50 CONTINUE
      RETURN
C
C EVALUATION OF INERFC(X),X.LT.0.,BY FORWARD RECURSION
C
   60 F1 = SNGL(DSUM)
   70 FZERO = F1
      DO 80 N=1,NMAX
        FM1 = F0
        F0 = F1
        F1 = (-X*F0+.5*FM1)/FLOAT(N)
        F(N) = F1
   80 CONTINUE
      RETURN
   90 IF (X.LT.0.) GO TO 110
  100 N0 = NMAX
      X0 = X
      GO TO 130
  110 IF (XSQ.LT.(ACC+1.)*AL10-.572) GO TO 120
      F0 = C*S
      F1 = 2.
      GO TO 70
  120 I = 2
      N0 = 1
      X0 = -X
C
C BACKWARD RECURRENCE ALGORITHM FOR EVALUATING FN(ABS(X)) FOR N=1,2,
C ...,N0, WHERE N0=NMAX IF X.GT.0. AND N0=1 IF X.LT.0.
C
  130 DO 140 N=1,N0
        R1(N) = 0.
  140 CONTINUE
      FN0 = FLOAT(N0)
      NU = IFIX((SQRT(FN0)+(2.3026*ACC+1.3863)/(2.8284*X0))**2-5.)
  150 NU = NU + 10
      IF (NU.GT.5000) GO TO 280
      N0P1 = N0 + 1
      NUM = NU - N0P1
      DO 160 N=1,N0
        R0(N) = R1(N)
  160 CONTINUE
      R = 0.
      DO 170 K=1,NUM
        N = NU - K
        R = .5/(X0+FLOAT(N+1)*R)
  170 CONTINUE
      DO 180 K=1,N0
        N = N0P1 - K
        R = .5/(X0+FLOAT(N+1)*R)
        R1(N) = R
  180 CONTINUE
      DO 190 N=1,N0
        IF (ABS(R1(N)-R0(N)).GT.EPS*ABS(R1(N))) GO TO 150
  190 CONTINUE
      F0 = .5*C/(X0+R1(1))
      FZERO = F0
      FF = F0
      DO 200 N=1,N0
        FF = R1(N)*FF
        F(N) = FF
  200 CONTINUE
      GO TO (210, 220), I
  210 RETURN
C
C STARTING VALUE FOR INERFC(X),X.LT.0.,PRIOR TO FORWARD RECURSION
C
  220 F1 = 2. - S*F0
      F0 = C*S
      GO TO 70
  230 IFLAG = 1
      RETURN
  240 IFLAG = 2
      RETURN
  250 IFLAG = 3
      RETURN
  260 IFLAG = 4
      RETURN
  270 IFLAG = 5
      RETURN
  280 IFLAG = 6
      RETURN
      END

		
.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]