0

How to find sum of evenly spaced Binomial coefficients modulo M?
ie. (nCa + nCa+r + nCa+2r + nCa+3r + ... + nCa+kr) % M = ?
given: 0 <= a < r, a + kr <= n < a + (k+1)r, n < 105, r < 100

My first attempt was: 

int res = 0;
int mod=1000000009;
for (int k = 0; a + r*k <= n; k++) {
    res = (res + mod_nCr(n, a+r*k, mod)) % mod;
}

but this is not efficient. So after reading here and this paper I found out the above sum is equivalent to:
summation[ω-ja * (1 + ωj)n / r], for 0 <= j < r; and ω = ei2π/r is a primitive rth root of unity.
What can be the code to find this sum in Order(r)?

Edit: n can go upto 105 and r can go upto 100.

Original problem source: https://www.codechef.com/APRIL14/problems/ANUCBC
Editorial for the problem from the contest: https://discuss.codechef.com/t/anucbc-editorial/5113
After revisiting this post 6 years later, I'm unable to recall how I transformed the original problem statement into mine version, nonetheless, I shared the link to the original solution incase anyone wants to have a look at the correct solution approach.

vaibhavatul47
  • 2,766
  • 4
  • 29
  • 42

4 Answers4

2

Binomial coefficients are coefficients of the polynomial (1+x)^n. The sum of the coefficients of x^a, x^(a+r), etc. is the coefficient of x^a in (1+x)^n in the ring of polynomials mod x^r-1. Polynomials mod x^r-1 can be specified by an array of coefficients of length r. You can compute (1+x)^n mod (x^r-1, M) by repeated squaring, reducing mod x^r-1 and mod M at each step. This takes about log_2(n)r^2 steps and O(r) space with naive multiplication. It is faster if you use the Fast Fourier Transform to multiply or exponentiate the polynomials.

For example, suppose n=20 and r=5. 

(1+x)    = {1,1,0,0,0}
(1+x)^2  = {1,2,1,0,0}
(1+x)^4  = {1,4,6,4,1}
(1+x)^8  = {1,8,28,56,70,56,28,8,1} 
           {1+56,8+28,28+8,56+1,70}
           {57,36,36,57,70}
(1+x)^16 = {3249,4104,5400,9090,13380,9144,8289,7980,4900}
           {3249+9144,4104+8289,5400+7980,9090+4900,13380}
           {12393,12393,13380,13990,13380}

(1+x)^20 = (1+x)^16 (1+x)^4
         = {12393,12393,13380,13990,13380}*{1,4,6,4,1}
           {12393,61965,137310,191440,211585,203373,149620,67510,13380}
           {215766,211585,204820,204820,211585}

This tells you the sums for the 5 possible values of a. For example, for a=1, 211585 = 20c1+20c6+20c11+20c16 = 20+38760+167960+4845.

Douglas Zare
  • 3,296
  • 1
  • 14
  • 21
0

Something like that, but you have to check a, n and r because I just put anything without regarding about the condition:

#include <complex>
#include <cmath>
#include <iostream>

using namespace std;

int main( void )
{
    const int r = 10;
    const int a = 2;
    const int n = 4;

    complex<double> i(0.,1.), res(0., 0.), w;

    for( int j(0); j<r; ++j )
    {
        w = exp( i * 2. * M_PI / (double)r );

        res += pow( w, -j * a ) * pow( 1. + pow( w, j ), n ) / (double)r;
    }

    return 0;

}
user2076694
  • 806
  • 1
  • 6
  • 10
  • this is working for small values of n. moreover i was surprised by the output because the result should definitely be a proper integer but using this formula we are getting a complex number with real part equal to the desired sum and some imaginary part! – vaibhavatul47 Apr 10 '14 at 23:26
  • I would compute and use `wj = exp( i * 2 * j * M_PI / r )`. – aschepler Apr 10 '14 at 23:34
  • yes, I forgot to say that the result is the Re(`res`) – user2076694 Apr 11 '14 at 11:32
  • try your calculation with `long long` because `int` is capped at `4E9` something like that – user2076694 Apr 11 '14 at 11:33
  • better int64_t or uint64_t than long long. Limits are guaranteed to be 9223372036854775807 and 18446744073709551615 respectively – Inuart Apr 13 '14 at 12:03
0

the mod operation is expensive, try avoiding it as much as possible

uint64_t res = 0;
int mod=1000000009;
for (int k = 0; a + r*k <= n; k++) {
    res += mod_nCr(n, a+r*k, mod);
    if(res > mod)
        res %= mod;
}

I did not test this code

Inuart
  • 1,432
  • 3
  • 17
  • 28
0

I don't know if you reached something or not in this question, but the key to implementing this formula is to actually figure out that w^i are independent and therefore can form a ring. In simpler terms you should think of implement (1+x)^n%(x^r-1) or finding out (1+x)^n in the ring Z[x]/(x^r-1) If confused I will give you an easy implementation right now.

  1. make a vector of size r . O(r) space + O(r) time

  2. initialization this vector with zeros every where O(r) space +O(r) time

  3. make the first two elements of that vector 1 O(1)

  4. calculate (x+1)^n using the fast exponentiation method. each multiplication takes O(r^2) and there are log n multiplications therefore O(r^2 log(n) )

  5. return first element of the vector.O(1) Complexity O(r^2 log(n) ) time and O(r) space. this r^2 can be reduced to r log(r) using fourier transform. How is the multiplication done, this is regular polynomial multiplication with mod in the power

    vector p1(r,0); vector p2(r,0); p1[0]=p1[1]=1; p2[0]=p2[1]=1; now we want to do the multiplication vector res(r,0); for(int i=0;i<r;i++) { for(int j=0;j<r;j++) { res[(i+j)%r]+=(p1[i]*p2[j]); } } return res[0]; I have implemented this part before, if you are still cofused about something let me know. I would prefer that you implement the code yourself, but if you need the code let me know.

Motaz Hammouda
  • 122
  • 1
  • 9