Logo

Programming-Idioms

  • Fortran
  • Rust
  • D

Idiom #67 Binomial coefficient "n choose k"

Calculate binom(n, k) = n! / (k! * (n-k)!). Use an integer type able to handle huge numbers.

import std.bigint;
BigInt binom(uint n, uint k)
{
   assert(n >= k);
   BigInt r = 1;
   for(uint i = 0; i <= k; ++i)
   {
      r *= n-i;
      r /= i+1;
   }
   return r;
}
integer, parameter :: i8 = selected_int_kind(18)
integer, parameter :: dp = selected_real_kind(15)
n = 100
k = 5
print *,nint(exp(log_gamma(n+1.0_dp)-log_gamma(n-k+1.0_dp)-log_gamma(k+1.0_dp)),kind=i8)

Sorry, Fortran does not do huge integers. The next best thing is probably to calculate the binomial coefficient using the Gamma function, or rather its logarithm, to get an approximation which it is then possible to round. Works for a rather large range.
extern crate num;

use num::bigint::BigInt;
use num::bigint::ToBigInt;
use num::traits::One;
fn binom(n: u64, k: u64) -> BigInt {
    let mut res = BigInt::one();
    for i in 0..k {
        res = (res * (n - i).to_bigint().unwrap()) /
              (i + 1).to_bigint().unwrap();
    }
    res
}
long long binom(long long n, long long k){
  long long outp = 1;
  for(long long i = n - k + 1; i <= n; ++i) outp *= i;
  for(long long i = 2; i <= k; ++i) outp /= i;
  return outp;
}

c does not have a BigInt type, as a substitute long long is used.
To increase the range as much as possible, the formula is slightly different than the one presented in the idiom, but it gives the same results.
New implementation...
< > ⌨
programming-idioms.org 
BigInt binom(uint n, uint k)
{
   assert(n >= k);
   BigInt r = 1;
   for(uint i = 0; i <= k; ++i)
   {
      r *= n-i;
      r /= i+1;
   }
   return r;
}
integer, parameter :: i8 = selected_int_kind(18)
integer, parameter :: dp = selected_real_kind(15)
n = 100
k = 5
print *,nint(exp(log_gamma(n+1.0_dp)-log_gamma(n-k+1.0_dp)-log_gamma(k+1.0_dp)),kind=i8)
fn binom(n: u64, k: u64) -> BigInt {
    let mut res = BigInt::one();
    for i in 0..k {
        res = (res * (n - i).to_bigint().unwrap()) /
              (i + 1).to_bigint().unwrap();
    }
    res
}
long long binom(long long n, long long k){
  long long outp = 1;
  for(long long i = n - k + 1; i <= n; ++i) outp *= i;
  for(long long i = 2; i <= k; ++i) outp /= i;
  return outp;
}
(defn binom [n k]
  (let [fact #(apply * (range 1 (inc %)))]
    (/ (fact n)
       (* (fact k) (fact (- n k))))))
public BigInteger binom(int n, int k)
{
	return factorial(n)/(factorial(k) * factorial(n-k));
}

public BigInteger factorial(int x)
{
	BigInteger result = 1;
	for(int i=1;i<=x;i++)
	{
	result = result * i;
	}
	return result;
}

int binom(int n, int k) {
  int result = 1;
  for (int i = 0; i < k; i++) {
    result = result * (n - i) ~/ (i + 1);
  }
  return result;
}
z := new(big.Int)
z.Binomial(n, k)
binom n k = product [1+n-k..n] `div` product [1..k]

const fac = x => x ? x * fac (x - 1) : x + 1
const binom = (n, k) => fac (n) / fac (k) / fac (n - k >= 0 ? n - k : NaN)
static BigInteger binom(int N, int K) {
    BigInteger ret = BigInteger.ONE;
    for (int k = 0; k < K; k++) {
        ret = ret.multiply(BigInteger.valueOf(N-k))
                 .divide(BigInteger.valueOf(k+1));
    }
    return ret;
}
gmp_binomial($n, $k);
var
  N, K, Res: TValue;
begin
  Res := BinomN(N, K);
end.
sub binom {
  my ($n, $k) = @_;
  if ($k > $n - $k) { $k = $n - $k }
  my $r = 1;
  for ( my $i = $n/$n ; $i <= $k;) {
    $r *= $n-- / $i++
  }
  return $r
}
sub binom {
   my ($n, $k) = @_;
   my $fact = sub {
      my $n = shift;
      return $n<2 ? 1 : $n * $fact->($n-1);
   };
   return $fact->($n) / ($fact->($k) * ($fact->($n-$k)));
}
def binom(n, k):
    return math.factorial(n) // math.factorial(k) // math.factorial(n - k)
def binom(n, k):
    return math.comb(n, k)
def binom(n,k)
  (1+n-k..n).inject(:*)/(1..k).inject(:*)
end