HDU 6333 莫队分块 + 逆元打表求组合数

HDU 6333 莫队分块 + 逆元打表求组合数

题解

在n个苹果中取最多m个苹果，问有多少中取法，即求\(\sum_{i=0}^m C_n^i\)，但是样例个数和n的范围最高到1e5, 一个一个求铁定会超时。

假设\(S(n, m) = \sum_{i=0}^m C_n^i = C_n^0 + C_n^1+C_n^2+...+C_n^m\)，则可以得到

\((1)\ S(n, m+1) = C_n^0 + C_n1+C_n2+...+C_n^m + C_n^{m+1}=S(n, m) + C_n^{m+1}\)
\((2)\ 2 * S(n, m) = C_n^0 + C_n^0 +C_n^1+C_n^1+C_n^2+C_n^2+...+C_n^m+C_n^m =C_n^0+(C_n^0+C_n^1)+(C_n^1+C_n^2)+(C_n^2+C_n^3) + ... +(C_n^{m-1}+C_n^{m})+C_n^m=C_{n+1}^0+C_{n+1}^1 +C_{n+1}^2 +C_{n+1}^3+...+C_{n+1}^m+C_n^m(杨辉三角)=S(n+1, m)\)
即\(S(n+1,m)=2*S(n,m)-C_n^m\)

由公式(1)(2)可得
\((3)\ S(n,m-1)=S(n,m)-C_n^m\)
\((4)\ S(n-1,m)=(S(n,m)+C_{n-1}^m)\ /\ 2\)

也就是说，我们可以通过\(S(n,m)\)和组合数\(C_n^m\)很容易的求出\(S(n, m+1), S(n+1, m), S(n,m-1),S(n-1,m)\)，于是我们就可以离线后，用莫队分块算法求解，将\(n\)分为\(sqrt(maxn)\)块，并且每一块内\(m\)以升序排列，莫队的具体细节参考百度。

还有一个问题是如何求\(C_n^m\)，考虑组合数的求解公式\(C_n^m=\frac{n!}{m!(n-m)!}\),最后结果对\(MOD=1e9+7\)取余，那么\(C_n^m\ \%\ MOD= fact(n)*inv(m)*inv(n-m)\%MOD\)，其中\(fact(x)\)表示x的阶乘，\(inv(x)\)表示\(x\)阶乘的逆元。我们可以用\(O(n)\)的算法打表求出\(frac[]\)和\(inv[]\)那么就可以在\(O(1)\)的复杂度内求出\(C_n^m\)

代码

#include <bits/stdc++.h>

using namespace std;

#define rep(i, a, b) for (int i = (a); i <= (b); ++i)
#define red(i, a, b) for (int i = (a); i >= (b); --i)
#define clr( x , y ) memset(x,y,sizeof(x))

typedef long long ll;

const int maxn = 1e5 + 5;
const int MOD  = 1e9 + 7;

ll l, r, m;
ll sum, ans[maxn], f[maxn + 1], inv[maxn + 1];

struct nodd
{
    ll l , r , n , k;
} b[maxn];

// 拓展欧几里得
ll exgcd(ll m, ll n, ll &x, ll &y){
    ll x1, y1, x0, y0;
    x0 = 1; y0 = 0;
    x1 = 0; y1 = 1;
    ll r = (m % n + n) % n;
    ll q = (m - r) / n;
    x = 0; y = 1;
    while(r){
        x = x0 - q * x1; y = y0 - q * y1; x0 = x1 ;y0 = y1;
        x1 = x; y1 = y;
        m = n; n = r; r = m % n;
        q = (m - r) / n;
    }
    return n;
}

void cal(){//阶乘及其逆元打表
    f[0]=1;
    for(int i=1;i<=maxn;i++){
        f[i]=f[i-1]*i%MOD;
    }

    ll x,y;
    exgcd(f[maxn],MOD,x,y);//先求出f[N]的逆元，再循环求出f[1~N-1]的逆元
    inv[maxn]=(x%MOD+MOD)%MOD;
    for(int i=maxn-1;i>=0;i--){
        inv[i]=inv[i+1]*(i+1)%MOD;
    }
}


// C(x, y);
ll c(ll x, ll y)
{
    return (f[x] * inv[y] % MOD * inv[x - y]) % MOD;
}

//同一块内按r的升序进行排序
bool cmp(nodd a,nodd b) { return a.k == b.k ? a.r < b.r : a.k < b.k; }

void init()
{
    cal();
    /*
    for (int i = 0; i < 100; ++i)
    {
        printf("%lld, %lld\n", f[i], inv[i]);
    }
    */
    scanf("%lld", &m);
    rep(i,1,m) scanf("%lld%lld",&b[i].l,&b[i].r);
    l = sqrt(maxn);
    //nodd的n属性记录其输入顺序, k属性记录其所在块数
    rep(i,1,m) b[i].n = i , b[i].k = b[i].l / l;
    //同一块内按r的升序进行排序
    sort( b + 1 , b + m + 1 , cmp );
}

// S(n,m)=S(n,m-1)+C(n,m)
void incm(ll n, ll m)
{
    //printf("c(%lld,%lld)=%lld\n", c(n,m));
    sum = (sum + c(n, m)) % MOD;
}
// S(n,m)=S(n,m+1)-C(n,m+1)
void decm(ll n, ll m)
{
    sum = (sum - c(n, m + 1) + MOD) % MOD;
}
// S(n,m)=(S(n+1,m)+C(n,m) / 2)
void decn(ll n, ll m)
{
    sum = ((sum + c(n, m)) * inv[2]) % MOD;
}
// S(n,m)=(S(n-1,m) * 2 - C(n-1, m))
void incn(ll n, ll m)
{
    sum = (sum * 2 - c(n - 1, m) + MOD) % MOD;
}

void work()
{
    rep(i,1,m) {
        sum = 0;
        l = b[i].l; r = b[i].r;
        //printf("%lld, %lld\n", l, r);
        rep(j, 0, r) incm(l, j);// printf("%lld\n", sum);
        ans[ b[i].n ] = sum;
        while ( i < m && b[i].k == b[i+1].k ) {
            i++;
            while ( l < b[i].l ) incn(++l, r);
            while ( l > b[i].l ) decn(--l, r);
            while ( r < b[i].r ) incm(l, ++r);
            ans[ b[i].n ] = sum;
        }
    }
    rep(i,1,m) {
        printf("%lld\n", ans[i]);
    }
}

int main()
{
    init();
    work();
    return 0;
}