m = ∏ (p_i ^ e_i)
F(p, n) = ∏ k = p^(n/p) * (n/p)!
1<=k<=n
p | k
G(p, n) = ∏ k
1<=k<=n
gcd(k,n)=1
n! = F(p, n) * G(p, n)
m
n! = p^K * ∏ G(p, n/(p^j))
j=0
n! / (r! * (n-r)!) = p^M * (∏ G(p, (n/p^j)) * [ ∏ G(p, r/(p^j)) * ∏ G(p, (n-r)/(p^j)) ]^(-1)
H(p, e, n) ≡ (-1)^(n/(p^e)) * H(p, e, n % (p^e)) (mod p^e)
// invert k modulo p, k and p are supposed coprime
unsigned long long invertMod(unsigned long long k, unsigned long long p) {
unsigned long long q, pn = 1, po = 0, r = p, s = k;
unsigned odd = 1;
do {
q = r/s;
q = pn*q + po;
po = pn;
pn = q;
q = r%s;
r = s;
s = q;
odd ^= 1;
}while(pn < p);
return odd ? p-po : po;
}
// Calculate the binomial coefficient (n choose k) modulo (prime^exponent)
// requires prime to be a prime, exponent > 0, and 0 <= k <= n,
// furthermore supposes prime^exponent < 2^32, otherwise intermediate
// computations could have mathematical results out of range.
// If k or (n-k) is small, a direct computation would be more efficient.
unsigned long long binmod(unsigned long long prime, unsigned exponent,
unsigned long long n, unsigned long long k) {
// The modulus, prime^exponent
unsigned long long ppow = 1;
// We suppose exponent is small, so that exponentiation by repeated
// squaring wouldn't gain much.
for(unsigned i = 0; i < exponent; ++i) {
ppow *= prime;
}
// array of remainders of products
unsigned long long *remainders = malloc(ppow * sizeof *remainders);
if (!remainders) {
fprintf(stderr, "Allocation failuren");
exit(EXIT_FAILURE);
}
for(unsigned long long i = 1; i < ppow; ++i) {
remainders[i] = i;
}
for(unsigned long long i = 0; i < ppow; i += prime) {
remainders[i] = 1;
}
for(unsigned long long i = 2; i < ppow; ++i) {
remainders[i] *= remainders[i-1];
remainders[i] %= ppow;
}
// Now to business.
unsigned long long pmult = 0, ntemp = n, ktemp = k, mtemp = n-k,
numer = 1, denom = 1, q, r, f;
if (prime == 2 && exponent > 2) {
f = 0;
} else {
f = 1;
}
while(ntemp) {
r = ntemp % ppow;
q = ntemp / ppow;
numer *= remainders[r];
numer %= ppow;
if (q & f) {
numer = ppow - numer;
}
ntemp /= prime;
pmult += ntemp;
}
while(ktemp) {
r = ktemp % ppow;
q = ktemp / ppow;
denom *= remainders[r];
denom %= ppow;
if (q & f) {
denom = ppow - denom;
}
ktemp /= prime;
pmult -= ktemp;
}
while(mtemp) {
r = mtemp % ppow;
q = mtemp / ppow;
denom *= remainders[r];
denom %= ppow;
if (q & f) {
denom = ppow - denom;
}
mtemp /= prime;
pmult -= mtemp;
}
// free memory before returning, we don't use it anymore
free(remainders);
if (pmult >= exponent) {
return 0;
}
while(pmult > 0) {
numer = (numer * prime) % ppow;
--pmult;
}
return (numer * invertMod(denom, ppow)) % ppow;
}
int factorial_prime_power(int f, int p) {
// When p is prime, returns k such that p^k divides f!
int k = 0;
while (f > p) {
f = f / p;
k = k + f;
}
return k;
}
int binomial_prime_power(int n, int r, int p) {
// when p is prime, returns k such that p^k divides nCr
return factorial_prime_power(n,p) - factorial_prime_power(r,p) - factorial_prime_power(n-r,p);
}
int powmod(int p, int k, int m) {
// quickly calculates p^k mod m
int res = 1;
int q = p;
int j = k;
while (j > 0) {
// invariant: p^k is congruent to res * q^j
if (j is odd) {
res = (res * q) % m;
j = (j-1)/2;
} else {
j = j / 2;
}
q = (q * q) % m;
}
return res;
}
int big_binomial(int n, int r, int m) {
if (n < r or r < 0) {
return 0;
}
int res = 1;
for(p in all primes from 2 to n) {
k = binomial_prime_power(n,r,p);
res = (res * powmod(p,k,m)) % m;
}
return res;
}
int factorial(int x) {
int i;
for(i=1; i<x; i++)
{
x *= i;
}
return x;
}
int main()
{
//use like this
factorial(n)/( factorial(r) - factorial(n-r) );
//do your stuffs
}