      SUBROUTINE BINCDF(X,P,N,CDF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE CUMULATIVE DISTRIBUTION
C              FUNCTION VALUE AT THE SINGLE PRECISION VALUE X
C              FOR THE BINOMIAL DISTRIBUTION
C              WITH SINGLE PRECISION 'BERNOULLI PROBABILITY'
C              PARAMETER = P, 
C              AND INTEGER 'NUMBER OF BERNOULLI TRIALS'
C              PARAMETER = N. 
C              THE BINOMIAL DISTRIBUTION USED
C              HEREIN HAS MEAN = N*P
C              AND STANDARD DEVIATION = SQRT(N*P*(1-P)).
C              THIS DISTRIBUTION IS DEFINED FOR ALL
C              DISCRETE INTEGER X BETWEEN 0 (INCLUSIVELY)
C              AND N (INCLUSIVELY).
C              THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION
C              F(X) = C(N,X) * P**X * (1-P)**(N-X).
C              WHERE C(N,X) IS THE COMBINATORIAL FUNCTION
C              EQUALING THE NUMBER OF COMBINATIONS OF N ITEMS
C              TAKEN X AT A TIME.
C              THE BINOMIAL DISTRIBUTION IS THE
C              DISTRIBUTION OF THE NUMBER OF
C              SUCCESSES IN N BERNOULLI (0,1)
C              TRIALS WHERE THE PROBABILITY OF SUCCESS
C              IN A SINGLE TRIAL = P.
C     INPUT  ARGUMENTS--X      = THE SINGLE PRECISION VALUE 
C                                AT WHICH THE CUMULATIVE DISTRIBUTION 
C                                FUNCTION IS TO BE EVALUATED.
C                                X SHOULD BE INTEGRAL-VALUED,
C                                AND BETWEEN 0.0 (INCLUSIVELY)
C                                AND N (INCLUSIVELY).
C                     --P      = THE SINGLE PRECISION VALUE 
C                                OF THE 'BERNOULLI PROBABILITY'
C                                PARAMETER FOR THE BINOMIAL 
C                                DISTRIBUTION.
C                                P SHOULD BE BETWEEN
C                                0.0 (EXCLUSIVELY) AND
C                                1.0 (EXCLUSIVELY).
C                     --N      = THE INTEGER VALUE
C                                OF THE 'NUMBER OF BERNOULLI TRIALS'
C                                PARAMETER.
C                                N SHOULD BE A POSITIVE INTEGER.
C     OUTPUT ARGUMENTS--CDF    = THE SINGLE PRECISION CUMULATIVE
C                                DISTRIBUTION FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION CUMULATIVE DISTRIBUTION
C             FUNCTION VALUE CDF
C             FOR THE BINOMIAL DISTRIBUTION
C             WITH 'BERNOULLI PROBABILITY' PARAMETER = P
C             AND 'NUMBER OF BERNOULLI TRIALS' PARAMETER = N.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--X SHOULD BE INTEGRAL-VALUED,
C                   AND BETWEEN 0.0 (INCLUSIVELY) 
C                   AND N (INCLUSIVELY).
C                 --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C                 --N SHOULD BE A POSITIVE INTEGER.
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--DSQRT, DATAN.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--NOTE THAT EVEN THOUGH THE INPUT
C              TO THIS CUMULATIVE
C              DISTRIBUTION FUNCTION SUBROUTINE
C              FOR THIS DISCRETE DISTRIBUTION
C              SHOULD (UNDER NORMAL CIRCUMSTANCES) BE A
C              DISCRETE INTEGER VALUE,
C              THE INPUT VARIABLE X IS SINGLE
C              PRECISION IN MODE.
C              X HAS BEEN SPECIFIED AS SINGLE
C              PRECISION SO AS TO CONFORM WITH THE DATAPAC
C              CONVENTION THAT ALL INPUT ****DATA****
C              (AS OPPOSED TO SAMPLE SIZE, FOR EXAMPLE)
C              VARIABLES TO ALL
C              DATAPAC SUBROUTINES ARE SINGLE PRECISION.
C              THIS CONVENTION IS BASED ON THE BELIEF THAT
C              1) A MIXTURE OF MODES (FLOATING POINT
C              VERSUS INTEGER) IS INCONSISTENT AND
C              AN UNNECESSARY COMPLICATION
C              IN A DATA ANALYSIS; AND
C              2) FLOATING POINT MACHINE ARITHMETIC
C              (AS OPPOSED TO INTEGER ARITHMETIC) 
C              IS THE MORE NATURAL MODE FOR DOING 
C              DATA ANALYSIS. 
C     REFERENCES--HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGE 38.
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 945, FORMULAE 26.5.24 AND
C                 26.5.28, AND PAGE 929.
C               --JOHNSON AND KOTZ, DISCRETE
C                 DISTRIBUTIONS, 1969, PAGES 50-86,
C                 ESPECIALLY PAGES 63-64.
C               --FELLER, AN INTRODUCTION TO PROBABILITY
C                 THEORY AND ITS APPLICATIONS, VOLUME 1,
C                 EDITION 2, 1957, PAGES 135-142. 
C               --KENDALL AND STUART, THE ADVANCED THEORY OF
C                 STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 120-125.
C               --MOOD AND GRABLE, INTRODUCTION TO THE THEORY
C                 OF STATISTICS, EDITION 2, 1963, PAGES 64-69.
C               --OWEN, HANDBOOK OF STATISTICAL
C                 TABLES, 1962, PAGES 264-272.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY (205.03)
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE:  301-921-2315
C     ORIGINAL VERSION--NOVEMBER  1975. 
C     UPDATED         --MAY       1977. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DX,PI,ANU1,ANU2,Z,SUM,TERM,AI,COEF1,COEF2,ARG
      DOUBLE PRECISION COEF
      DOUBLE PRECISION THETA,SINTH,COSTH,A,B
      DOUBLE PRECISION DSQRT,DATAN
      DATA PI/3.14159265358979D0/
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      AN=N
      IF(P.LT.0.0.OR.P.GT.1.0)GOTO50
      IF(N.LT.1)GOTO55
      IF(X.LT.0.0.OR.X.GT.AN)GOTO60
      INTX=X+0.0001 
      FINTX=INTX
      DEL=X-FINTX
      IF(DEL.LT.0.0)DEL=-DEL
      IF(DEL.GT.0.001)GOTO65
      GOTO90
   50 WRITE(IPR,11) 
      WRITE(IPR,46)P
      CDF=0.0
      RETURN
   55 WRITE(IPR,25) 
      WRITE(IPR,47)N
      CDF=0.0
      RETURN
   60 WRITE(IPR,4)N 
      WRITE(IPR,46)X
      IF(X.LT.0.0)CDF=0.0
      IF(X.GT.AN)CDF=1.0
      RETURN
   65 WRITE(IPR,5)
      WRITE(IPR,46)X
   90 CONTINUE
    4 FORMAT(1H ,111H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT TO THE BINCDF SUBROUTINE IS OUTSIDE THE USUAL (0,N) = (0,,I7,
     1 11H,INTERVAL *)
    5 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT TO THE BINCDF SUBROUTINE IS NON-INTEGRAL *****)
   11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 BINCDF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 BINCDF SUBROUTINE IS NON-POSITIVE *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
C     TREAT IMMEDIATELY THE SPECIAL CASE OF X = N,
C     IN WHICH CASE CDF = 1.0.
C     ALSO TREAT IMMEDIATELY THE SPECIAL CASE OF P = 0.0
C     IN WHICH CASE CDF = 1.0 FOR ALL X.
C     THIRDLY, TREAT THE SPECIAL CASE IN WHICH P = 1.0
C     IN WHICH CASE CDF = 0.0 FOR ALL X SMALLER THAN N
C     AND CDF = 1.0 FOR ALL X EQUAL TO OR LARGER
C     THAN N.
C
      INTX=X+0.0001 
      CDF=1.0
      IF(INTX.EQ.N)RETURN
      IF(P.EQ.0.0)RETURN
      IF(P.EQ.1.0.AND.INTX.GE.N)RETURN
      IF(P.EQ.1.0.AND.INTX.LT.N)CDF=0.0 
      IF(P.EQ.1.0.AND.INTX.LT.N)RETURN
C
C     EXPRESS THE BINOMIAL CUMULATIVE DISTRIBUTION
C     FUNCTION IN TERMS OF THE EQUIVALENT F
C     CUMULATIVE DISTRIBUTION FUNCTION, 
C     AND THEN EVALUATE THE LATTER.
C
      AN=N
      DX=(P/(1.0-P))*((AN-X)/(X+1.0))
      NU1=2.0*(X+1.0)+0.1
      NU2=2.0*(AN-X)+0.1
      ANU1=NU1
      ANU2=NU2
      Z=ANU2/(ANU2+ANU1*DX)
C
C     DETERMINE IF NU1 AND NU2 ARE EVEN OR ODD
C
      IFLAG1=NU1-2*(NU1/2)
      IFLAG2=NU2-2*(NU2/2)
      IF(IFLAG1.EQ.0)GOTO120
      IF(IFLAG2.EQ.0)GOTO150
      GOTO250
C
C     DO THE NU1 EVEN AND NU2 EVEN OR ODD CASE
C
  120 SUM=0.0D0
      TERM=1.0D0
      IMAX=(NU1-2)/2
      IF(IMAX.LE.0)GOTO110
      DO100I=1,IMAX 
      AI=I
      COEF1=2.0D0*(AI-1.0D0)
      COEF2=2.0D0*AI
      TERM=TERM*((ANU2+COEF1)/COEF2)*(1.0D0-Z)
      SUM=SUM+TERM
  100 CONTINUE
C
  110 SUM=SUM+1.0D0 
      SUM=(Z**(ANU2/2.0D0))*SUM
      CDF=SUM
      RETURN
C
C     DO THE NU1 ODD AND NU2 EVEN CASE
C
  150 SUM=0.0D0
      TERM=1.0D0
      IMAX=(NU2-2)/2
      IF(IMAX.LE.0)GOTO210
      DO200I=1,IMAX 
      AI=I
      COEF1=2.0D0*(AI-1.0D0)
      COEF2=2.0D0*AI
      TERM=TERM*((ANU1+COEF1)/COEF2)*Z
      SUM=SUM+TERM
  200 CONTINUE
C
  210 SUM=SUM+1.0D0 
      CDF=1.0D0-((1.0D0-Z)**(ANU1/2.0D0))*SUM
      RETURN
C
C     DO THE NU1 ODD AND NU2 ODD CASE
C
  250 SUM=0.0D0
      TERM=1.0D0
      ARG=DSQRT((ANU1/ANU2)*DX)
      THETA=DATAN(ARG)
      SINTH=ARG/DSQRT(1.0D0+ARG*ARG)
      COSTH=1.0D0/DSQRT(1.0D0+ARG*ARG)
      IF(NU2.EQ.1)GOTO320
      IF(NU2.EQ.3)GOTO310
      IMAX=NU2-2
      DO300I=3,IMAX,2
      AI=I
      COEF1=AI-1.0D0
      COEF2=AI
      TERM=TERM*(COEF1/COEF2)*(COSTH*COSTH)
      SUM=SUM+TERM
  300 CONTINUE
C
  310 SUM=SUM+1.0D0 
      SUM=SUM*SINTH*COSTH
C
  320 A=(2.0D0/PI)*(THETA+SUM)
      SUM=0.0D0
      TERM=1.0D0
      IF(NU1.EQ.1)B=0.0D0
      IF(NU1.EQ.1)GOTO450
      IF(NU1.EQ.3)GOTO410
      IMAX=NU1-3
      DO400I=1,IMAX,2
      AI=I
      COEF1=AI
      COEF2=AI+2.0D0
      TERM=TERM*((ANU2+COEF1)/COEF2)*(SINTH*SINTH)
      SUM=SUM+TERM
  400 CONTINUE
C
  410 SUM=SUM+1.0D0 
      SUM=SUM*SINTH*(COSTH**N)
      COEF=1.0D0
      IEVODD=NU2-2*(NU2/2)
      IMIN=3
      IF(IEVODD.EQ.0)IMIN=2
      IF(IMIN.GT.NU2)GOTO420
      DO430I=IMIN,NU2,2
      AI=I
      COEF=((AI-1.0D0)/AI)*COEF
  430 CONTINUE
C
  420 COEF=COEF*ANU2
      IF(IEVODD.EQ.0)GOTO440
      COEF=COEF*(2.0D0/PI)
C
  440 B=COEF*SUM
C
  450 CDF=1.0D0-(A-B)
      RETURN
C
      END 
