bin Function

        pure function bin(n, r)
            integer(i8p)        :: bin
            integer, intent(in) :: n
            integer, intent(in) :: r
        
            integer(i8p)        :: num
            integer(i8p)        :: den
            integer             :: i
            integer             :: k
            integer, parameter  :: primes(*) = [2,3,5,7,11,13,17,19]
            num = 1
            den = 1
            do i=0,r-1
                num = num*(n-i)
                den = den*(i+1)
                if (i > 0) then
                    ! Divide out common prime factors
                    do k=1,size(primes)
                        if (mod(i,primes(k)) == 0) then
                            num = num/primes(k)
                            den = den/primes(k)
                        end if
                    end do
                end if
            end do
            bin = num/den
        end function bin