BinomialCoefficient.cpp

#include <Base/Math/BinomialCoefficient.hh>
#include <Base/Debug/Exception.hh>
#include <Base/Debug/Error.hh>

#include <limits>
#include <limits.h>

using namespace BIAS;
using namespace std;

long long unsigned BinomialCoefficient::
Compute(const unsigned n, const unsigned k)
{
 if ((k==0) || (k==n) || (k==0 && n==0)){
    return 1u; 
  }
  if (2*k>n) {
    return Compute(n, n-k); 
  } 
  long long unsigned prod = 1, old_prod;
  BIASASSERT(n>=k);
  const unsigned n_minus_k = n-k;
  for (unsigned i=1; i<=n_minus_k; i++){
    old_prod = prod;
    prod *= (long long unsigned)(k+i);
    if (prod<=old_prod){
      BEXCEPTION("BinomialCoefficient::Compute("<<n<<", "<<k<<"): Overflow");
    }
    BIASASSERT(prod%i==0);
    prod /= (long long unsigned)(i);
    // changed if statement from  if (prod>=numeric_limits<unsigned>::max())
    if (prod>=UINT_MAX){
      BEXCEPTION("BinomialCoefficient::Compute("<<n<<", "<<k<<"): Overflow");
    }
  }
  return prod;
}