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
    }