c++ - To find combination value of large numbers

Question

I want to find (n choose r) for large integers, and I also have to find out the mod of that number.

long long int choose(int a,int b)
{
    if (b > a)
        return (-1);
    if(b==0 || a==1 || b==a)
        return(1);
    else
    {
        long long int r = ((choose(a-1,b))%10000007+(choose(a-1,b-  1))%10000007)%10000007;
        return r;
    }
}

I am using this piece of code, but I am getting TLE. If there is some other method to do that please tell me.

Answer 1

I don't have the reputation to comment yet, but I wanted to point out that the answer by rock321987 works pretty well:

It is fast and correct up to and including C(62, 31)

but cannot handle all inputs that have an output that fits in a uint64_t. As proof, try:

C(67, 33) = 14,226,520,737,620,288,370 (verify correctness and size)

Unfortunately, the other implementation spits out 8,829,174,638,479,413 which is incorrect. There are other ways to calculate nCr which won't break like this, however the real problem here is that there is no attempt to take advantage of the modulus.

Notice that p = 10000007 is prime, which allows us to leverage the fact that all integers have an inverse mod p, and that inverse is unique. Furthermore, we can find that inverse quite quickly. Another question has an answer on how to do that here, which I've replicated below.

This is handy since:

1. x/y mod p == x*(y inverse) mod p; and
2. xy mod p == (x mod p)(y mod p)

Modifying the other code a bit, and generalizing the problem we have the following:

#include <iostream>
#include <assert.h>

// p MUST be prime and less than 2^63
uint64_t inverseModp(uint64_t a, uint64_t p) {
    assert(p < (1ull << 63));
    assert(a < p);
    assert(a != 0);
    uint64_t ex = p-2, result = 1;
    while (ex > 0) {
        if (ex % 2 == 1) {
            result = (result*a) % p;
        }
        a = (a*a) % p;
        ex /= 2;
    }
    return result;
}

// p MUST be prime
uint32_t nCrModp(uint32_t n, uint32_t r, uint32_t p)
{
    assert(r <= n);
    if (r > n-r) r = n-r;
    if (r == 0) return 1;
    if(n/p - (n-r)/p > r/p) return 0;

    uint64_t result = 1; //intermediary results may overflow 32 bits

    for (uint32_t i = n, x = 1; i > r; --i, ++x) {
        if( i % p != 0) {
            result *= i % p;
            result %= p;
        }
        if( x % p != 0) {
            result *= inverseModp(x % p, p);
            result %= p;
        }
    }
    return result;
}

int main() {
    uint32_t smallPrime = 17;
    uint32_t medNum = 3001;
    uint32_t halfMedNum = medNum >> 1;
    std::cout << nCrModp(medNum, halfMedNum, smallPrime) << std::endl;

    uint32_t bigPrime = 4294967291ul; // 2^32-5 is largest prime < 2^32
    uint32_t bigNum = 1ul << 24;
    uint32_t halfBigNum = bigNum >> 1;
    std::cout << nCrModp(bigNum, halfBigNum, bigPrime) << std::endl;
}

Which should produce results for any set of 32-bit inputs if you are willing to wait. To prove a point, I've included the calculation for a 24-bit n, and the maximum 32-bit prime. My modest PC took ~13 seconds to calculate this. Check the answer against wolfram alpha, but beware that it may exceed the 'standard computation time' there.

There is still room for improvement if p is much smaller than (n-r) where r <= n-r. For example, we could precalculate all the inverses mod p instead of doing it on demand several times over.