3

I would like to optimize a part of my program where I'm calculating the sum of Binomial Coefficients up to K. i.e. 

C(N,0) + C(N,1) + ... + C(N,K)

Since the values go beyond the data type (long long) can support, I'm to calculate values mod M and was looking for procedures to do that.

Currently, I've done it with Pascal's Triangle but it seems to be taking a bit of load. so, I was wondering if there's any other efficient way to do this. I've considered Lucas' Theorem, although M I have is already large enough so that C(N,k) goes out of hand!

Any pointers as how can I do this differently, maybe calculate the whole sum altogether with some other neat expression of teh sum. If not I'll leave it with the Pascal's Triangle method itself.

Thank you,

Here is what I have so far O(N^2) :

#define MAX 1000000007
long long NChooseK_Sum(int N, int K){
    vector<long long> prevV, V;
    prevV.push_back(1);     prevV.push_back(1);
    for(int i=2;i<=N;++i){
            V.clear();
            V.push_back(1);
            for(int j=0;j<(i-1);++j){
                    long long val = prevV[j] + prevV[j+1];
                    if(val >= MAX)
                            val %= MAX;
                    V.push_back(val);
            }
            V.push_back(1);
            prevV = V;
    }
    long long res=0;
    for(int i=0;i<=K;++i){
            res+=V[i];
            if(res >= MAX)
                    res %= MAX;
    }
    return res;
}
srbhkmr
  • 2,074
  • 1
  • 14
  • 19
  • I've always used Pascal's triangle, I think you get rounding errors quite soon especially with 32 bit ints but it worked fine upto about the 20th row so was fine for my needs. – John Feb 21 '12 at 00:04
  • What are the approximate value ranges for *N*, *K* and *M*? Also, can you read SML-caml-F#? I have code in F# if that works for you at all. – kkm inactive - support strike Feb 21 '12 at 00:44
  • If `K` is much smaller than `N`, you can gain quite a bit by stopping the inner loop at `K`, also if `K` is close to `N`, by stopping at `N-K` and using the fact that the sum of all binomial coefficients is `2^N`. But if you really need it fast, part deux' suggestion (with the modular inverses) gets you the sum (modulo `MAX`) in O(K*log(min(K,MAX))) steps. (Some care is needed if `K >= MAX`.) – Daniel Fischer Feb 21 '12 at 01:21
  • You can definitely do better by changing from vector to a static array.. might spend too much time allocating. – Larry Feb 21 '12 at 01:52
  • @kkm - I'm afraid, I've no knowledge of SML-caml-F# but you are welcome to wite your code with the approach here. 1<=K <= N<=100000, M=1000000007 – srbhkmr Feb 21 '12 at 03:26
  • @DanielFischer - Those are good observations. I'll keep them in mind. K in my case is well below MAX. I'll go with the modular inverse approach. – srbhkmr Feb 21 '12 at 03:32
  • @Larry - Thanks, I'll shift to static array. – srbhkmr Feb 21 '12 at 03:33
  • @srbh.kmr: Sorry, I do not think my approach would help you. My problem was different; I have a solution in O(*F* √ *N*), where *F* is the largest prime factor in *M*. but in your case O(*N*) is better. The idea is that for each prime *p* ≤ *F*, you calculate the power of that prime factor in C[*N*, *K*] mod q_i (Goetgheluck 1987), where q_i is the prime factor in *M*. Regular powermod multiplication then calculates C[*N*, *K*] mod q_i. After that, the Chinese reminder calculated the answer given partial moduli. But this fails to improve the case of M>>N. – kkm inactive - support strike Feb 21 '12 at 08:20
  • @kkm Oh I see. nevertheless its' interesting to know. – srbhkmr Feb 21 '12 at 22:08
  • You can actually compute `C(n,k) % m` in `O(n)` time for any values of `m`. And it's not too complicated either, see my answer at: http://stackoverflow.com/a/24500377/205521 – Thomas Ahle Jun 30 '14 at 23:34

1 Answers1

5

An algorithm that performs a linear number of arithmetic bignum operations is

def binom(n):
    nck = 1
    for k in range(n + 1):  # 0..n
        yield nck
        nck = (nck * (n - k)) / (k + 1)

This uses division, but modulo a prime p, you can accomplish much the same thing by multiplying by the solution i to the equation i * (k + 1) = 1 mod p. The value i can be found in a logarithmic number of arithmetic ops via the extended Euclidean algorithm.

part deux
  • 66
  • 1
  • Yes, Thank you, that should work out neatly. Additions as in Pascal's triangle work out okay but with division in this approach, I was kind of hoping there has to be something to straighten it out in the mod ring. Thanks for pointing out modular inverse and the ext euclidean. – srbhkmr Feb 21 '12 at 03:20