#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);
}