Commit 819d77bd authored by ulrich_y's avatar ulrich_y

Added function to calculate nCr

This is taken from RosettaCode
parent ee4a63a3
......@@ -117,6 +117,32 @@ CONTAINS
enddo
END FUNCTION factorial
! Adapted from https://rosettacode.org/wiki/Evaluate_binomial_coefficients#Fortran
! published under GNU Free Documentation License 1.2
FUNCTION binom(n, r)
integer, intent(in) :: n, r
integer :: binom
integer(16) :: num, den
integer i, k
integer, parameter :: primes(20) = (/2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71/)
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
binom = int(num/den)
END FUNCTION binom
! subroutine print(s1,s2,s3,s4,s5)
! character(len = *), intent(in), optional :: s1, s2, s3, s4, s5
! if(print_enabled) then
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment