#include "jestatistics.h"

// LICENSE: This code is freeware.  You may use this code for any purpose you
// like so long as this license is attached.  If this code is used in a
// commercial application, you must send me an email so I feel uber-geeky.  No
// other compensation is required.  Email me at jechols@nerdbucket.com

// This is the combination function.  It is useful as an external function
// not part of any classes, so here it is.  See
// http://en.wikipedia.org/wiki/Combinatorics for general combinatorics info.
//
// We rely on the GMP libraries for this, since combin can return insanely
// large values.
//
// Also note this should always return a whole number.  Since combinations are
// a number of ways to pick r objects out of a group of n, it's not possible to
// have a float.
//
// Also note that this might be able to go faster by factoring out some of
// the factorials - I don't know the speed of the mpz_fac_ui function, so it
// may be super-optimized, but 80! / (20! * 60!) is generally going to be
// slower than if you don't bother to compute 60! at all, and just do the
// product of 80 through 61 divided by product of 20 through 2.
mpz_class JEStatistics::combin(int n, int r) {
  mpz_class fac_n, fac_r, fac_n_minus_r, top, bottom;

  // Sanity checking is necessary here.
  if (n < r || n < 0 || r < 0) {return mpz_class(0);}

  mpz_fac_ui(fac_n.get_mpz_t(), n);
  mpz_fac_ui(fac_r.get_mpz_t(), r);
  mpz_fac_ui(fac_n_minus_r.get_mpz_t(), n - r);

  return fac_n / (fac_r * fac_n_minus_r);
}