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.

Is this for some code challenge? The constant `10000007` looks somewhat familiar: http://stackoverflow.com/questions/14189713/how-to-calculate-abc10000007-where-a-b-c-can-be-at-maximum-1018 — Codor, Jun 27 '14 at 06:02
Very inefficient way to calculate binom.coeff. You say you need to do calculations for *large numbers* – have you tried to estimate how many times your `choose()` must be called to calculate `choose(n,n/2)` for large n? — CiaPan, Jun 27 '14 at 06:10
Are you sure you have the modulus right? `100000007` is a more common modulus is programming contests since you are often expected to exploit the fact that it is prime. Using a memoized version of your recursion, your run time will be `O(a*b)`. If the modulus is a prime, p, you can solve it in `O((a+b)*log p)`. — Pradhan, Jun 27 '14 at 06:14

score 5 · Answer 1 · edited May 23 '17 at 12:26

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:

x/y mod p == x*(y inverse) mod p; and
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.

score 4 · Answer 2 · edited Sep 16 '19 at 00:01

4

nCr = n! / (r! * (n-r)!) {! = factorial}

now choose r or n - r in such a way that any of them is minimum

#include <cstdio>
#include <cmath>

#define MOD 10000007

int main()
{
    int n, r, i, x = 1;
    long long int res = 1;
    scanf("%d%d", &n, &r);
    int mini = fmin(r, (n - r));//minimum of r,n-r

    for (i = n;i > mini;i--) {
        res = (res * i) / x;
        x++;
    }
    printf("%lld\n", res % MOD);
    return 0;
}

it will work for most cases as required by programming competitions if the value of n and r are not too high

Time complexity :- O(min(r, n - r))

Limitation :- for languages like C/C++ etc. there will be overflow if

n > 60 (approximately)

as no datatype can store the final value..

edited Sep 16 '19 at 00:01

Vineel

436
6
12

answered Jun 27 '14 at 06:20

rock321987

10,942
1
30
43

The central binomial coefficient grows [exponentially](http://en.wikipedia.org/wiki/Binomial_coefficient#Bounds_and_asymptotic_formulas). So, `res` will quickly overflow. You can't simply do `res%MOD` in the end and expect it to work. – Pradhan Jun 27 '14 at 06:43
that's the reason why I had mentioned if the constraints on `n` and `k` are not too high – rock321987 Jun 27 '14 at 06:45

score 0 · Answer 3 · answered Oct 25 '19 at 12:19

The expansion of nCr can always be reduced to product of integers. This is done by canceling out terms in denominator. This approach is applied in the function given below.

This function has time complexity of O(n^2 * log(n)). This will calculate nCr % m for n<=10000 under 1 sec.

#include <numeric>
#include <algorithm>
int M=1e7+7;
int ncr(int n, int r)
{
    r=min(r,n-r);
    int A[r],i,j,B[r];
    iota(A,A+r,n-r+1);  //initializing A starting from n-r+1 to n
    iota(B,B+r,1);      //initializing B starting from 1 to r

    int g;
    for(i=0;i<r;i++)
    for(j=0;j<r;j++)
    {
        if(B[i]==1)
            break;
        g=__gcd(B[i], A[j] );
        A[j]/=g;
        B[i]/=g;
    }
    long long ans=1;
    for(i=0;i<r;i++)
        ans=(ans*A[i])%M;
    return ans;
}

To find combination value of large numbers

3 Answers3

Linked