!     Last change:  DOS  29 Jul 2000   11:57 am
!                                       *** copyright 2000 ***
! *** filename ibicof.f95               *** John F. Monahan **
!                                       **********************
      PROGRAM PIBICOF
!  TEST PROGRAM FOR IBICOF - COMPUTES BINOMIAL COEFFICIENTS -
!     USUALLY FROM TABLE
      IMPLICIT NONE
      INTEGER S,P,I,J,KRED,KBLUE
      INTEGER, PARAMETER :: NSMALL = 5  ! NUMBER OF BALLS DRAWN
      INTEGER IBICOF
!                                                                       
  21  FORMAT(' FIRST TEST -- IS SUM OF BINOMIAL COEFFICIENTS 2**N ?' )
  22  FORMAT(' SECOND TEST -- HYPERGEOMETRIC SUM ' )
  23  FORMAT(I6,2I20)
  24  FORMAT(2I4,2I20)
!
      OPEN( UNIT=6, FILE='ibicof.out' )
!                       FIRST TEST - DO BINOM COEFFS SUM TO 2**?        
!
      P = 1
      DO I = 1,25
      P = P*2                                                           
      S = IBICOF(I,0)                                                   
      DO J = 1,I
      S = S + IBICOF(I,J)                                               
      END DO  ! LOOP ON J
      WRITE(6,23) I,S,P
      END DO  ! LOOP ON I
!                       NEGATIVE BINOMIAL TEST - PROBS SUM TO 1 ?
      DO KRED = 3,15,3
      DO KBLUE = 4,16,4
      S = 0
!                   DON'T CHECK BOUNDS ON J -- IBICOF SHOULD DO IT OK
      DO J = 0,KRED              ! HOW MANY WAYS TO GET J RED BALLS
      S = S + IBICOF(KRED,J)*IBICOF(KBLUE,NSMALL-J)
      END DO  ! LOOP ON J
      P = IBICOF(KRED+KBLUE,NSMALL)
      WRITE(6,24) KRED,KBLUE,S,P
      END DO  ! LOOP ON KRED
      END DO  ! LOOP ON KBLUE
!
      STOP
      END PROGRAM PIBICOF
      INTEGER FUNCTION IBICOF(I,J)                                      
!  COMPUTES BINOMIAL COEFFICIENTS, USUALLY BY TABLE                     
!
      INTEGER, INTENT(IN) :: I,J
      INTEGER K,L,N
      INTEGER, DIMENSION(21) :: PUSH    ! TABLE OF POINTERS
      INTEGER, DIMENSION(122) :: BICFTB
!                                                                       
!     TEST WHETHER IT IS IN TABLE -
!         REMEMBER THE RESULT WILL GO HAYWIRE IN INTEGERS IF TOO BIG
!
      DATA PUSH/1,2,4,6,8,11,14,18,22,27,32,38,44,51,58,66,74,83,92,  &
     &          102,112/
      DATA BICFTB/1,1,1,1,2,1,3,1,4,6,1,5,10,1,6,15,20,1,7,21,35,     &
     &      1,8,28,56,70,1,9,36,84,126,1,10,45,120,210,252,           &
     &      1,11,55,165,330,462,1,12,66,220,495,792,924,              &
     &      1,13,78,286,715,1287,1716,1,14,91,364,1001,2002,3003,3432,&
     &  1,15,105,455,1365,3003,5005,6435,1,16,120,560,1820,4368,8008, &
     & 11440,12870, 1,17,136,680,2380,6188,12376,19448,24310,1,18,153,&
     & 816,3060,8568,18564,31824,43758,48620, 1,19,171,969,3876,11628,&
     & 27132,50388,75582,92378, 1,20,190,1140,4845,15504,38760,77520, &
     & 125970, 167960, 184756 /
!
!                   TAKE CARE OF SPECIAL CASES
      IBICOF = 0
      IF( ( J .LT. 0 ) .OR. ( J .GT. I ) ) RETURN
      IBICOF = 1
      IF( ( J .EQ. 0 ) .OR. ( J .EQ. I ) ) RETURN
!                   GENERAL CASE
      L = MIN(J,I-J)
      IF( I .LE. 20 ) THEN    ! TABLE LOOKUP
      K = PUSH(I+1)           !  FIND START POINT FOR I
      IBICOF = BICFTB(K+L)    !  PULL NUMBER FROM TABLE
                      ELSE    ! MULTIPLY
!                         FOR BIG VALUES COMPUTE BY MULTIPLICATION
      IBICOF = 1
      DO N = 1,L
      IBICOF = (IBICOF*(I-N+1)) / N                                     
      END DO  ! LOOP ON N
                      END IF ! ( I .LE. 20 )
      RETURN
      END FUNCTION IBICOF
! *** end of filename ibicof.f95        *********************