TestBinomialCoefficient.cpp

#include <Base/Math/BinomialCoefficient.hh>

#include <Base/Math/Random.hh>
#include <Base/Debug/Exception.hh>

#include <iostream>

using namespace BIAS;
using namespace std;

unsigned binomial_coeff(const unsigned n, const unsigned k);


int main(int argc, char *argv[])
{
  const unsigned num = 1000;
  const unsigned max_n = 25;
  string arg1;
  if (argc==2) arg1 = argv[1];
  const bool verbose = (arg1!="silent");

  Random MyRandom;

  for (unsigned t=0; t<num; t++){
    const unsigned n = MyRandom.GetUniformDistributedInt(0u, max_n);
    const unsigned k = MyRandom.GetUniformDistributedInt(0u, n);
    try {
      #if  __GNUC_PREREQ(4,6)
        #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
      #endif
      long long unsigned bc = 0;
      bc = BinomialCoefficient::Compute(n, k);
      #if  __GNUC_PREREQ(4,6)
        #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
      #endif
      unsigned bc_rec = 0;
      bc_rec = binomial_coeff(n, k);
      BIASASSERT(bc == bc_rec);
    } catch (BaseException &ex) {
      cout << "cannot compute "<<n<<" over "<<k<<": overflow:"<<ex.What()<<endl;
    }
    if (verbose && t%10==0 ){
      cout << "."<<flush;
    }
    
  }
  if (verbose){
    cout << endl;
  }
}

unsigned binomial_coeff(const unsigned n, const unsigned k)
{
  if ((k==0) || (k==n) || (k==0 && n==0)){
    return 1u;
  }
  assert(n>0 && k>0);
  unsigned res1 = binomial_coeff(n-1, k-1);
  unsigned res2 = binomial_coeff(n-1, k);
  unsigned binomial_coefficent = res1 + res2;
  if (binomial_coefficent < res1 || 
      binomial_coefficent < res2){
    BEXCEPTION("Overflow\n");
  }
  return binomial_coefficent;
}