组合数前缀和 & 题解：P5388 [Cnoi2019] 最终幻想

· · 题解

题解怎么都在分段打表啊，这里给出一种更加友好的求解方式（指对电脑友好，不需要打表）。

首先，我们要明确这道题要求的是 \sum\limits^n_{i=0}\binom ki，但是这篇文章的重点不在推导出这个答案，详细的推导过程可以看其他题解。

这里给出一种时间复杂度为 O(\sqrt n\log n) 的组合数前缀和算法。

思路与快速阶乘算法一致，同样先分块，我们发现 \binom n{m+1}=\frac{n-m}{m+1}\binom nm，再把它写成矩阵形式，带上前缀和，可以得到

\binom n{m+1}\\ \sum\limits^m_{i=0}\binom ni \end{bmatrix}=\left(\prod\limits^m_{i=0} \begin{bmatrix} \frac{n-i}{i+1}&0\\ 1&1 \end{bmatrix}\right) \begin{bmatrix} 1\\ 0 \end{bmatrix}

注意上面的乘积顺序运算是从右向左的（即每次把矩阵放到左边和之前的结果相乘）。

但是这个矩阵中带有分数，并不好计算，所以先通分，再把分母提出来，即 (m+1)!，化简后得到

\begin{bmatrix} n-i&0\\ i+1&i+1 \end{bmatrix}\right) \begin{bmatrix} 1\\ 0 \end{bmatrix}

现在只需要维护 3 个点值 f,g,h，倍增计算

f_{2m}(x)&0\\ g_{2m}(x)&h_{2m}(x) \end{bmatrix}=\begin{bmatrix} f_{m}(x+m)&0\\ g_{m}(x+m)&h_{m}(x+m) \end{bmatrix} \begin{bmatrix} f_{m}(x)&0\\ g_{m}(x)&h_{m}(x) \end{bmatrix}

手动计算一下矩阵乘法可以得到

f_m(x)f_m(x+m)&0\\ f_m(x)g_m(x+m)+g_m(x)h_m(x+m)&h_m(x)h_m(x+m) \end{bmatrix}

直接每轮跑 6 次点值平移，倍增即可，实现时取块长为 2^{\lceil\log_2\sqrt m\rceil}，这样可以避免出现从 m 转移到 m+1 的情况。

最后答案要除以 (m+1)!，但此时我们不需要单独计算一次阶乘，上面的 h 实际上维护的就是阶乘的点值，同时算出来就好了。

参考实现：

int sumcomb(int n,int m,const int mod){
    m=min(n,m)+1;
    // 注意这里，上面式子中乘积下标是从 0 开始的，但是这里计算从 1 开始，所以加上 1 再计算
    MPoly f(2),g(2),h(2);
    f.setmod(mod),g.setmod(f.fm),h.setmod(f.fm);
    int B=1<<int(ceil(0.5*log2(m))),r=qpow(B,mod-2,f.fm),siz=m/B;
    f.a[0]=n,f.a[1]=n-B;
    g.a[0]=1,g.a[1]=B+1;
    h.a[0]=1,h.a[1]=B+1; // 初始值
    for (int l=2;l<=B;l<<=1){
        f.extend(pointValueTranslation(f,(l>>1)+1));
        g.extend(pointValueTranslation(g,(l>>1)+1));
        h.extend(pointValueTranslation(h,(l>>1)+1));
        MPoly u=pointValueTranslation(f,r*ll(l>>1)%f.fm);
        MPoly v=pointValueTranslation(g,r*ll(l>>1)%f.fm);
        MPoly w=pointValueTranslation(h,r*ll(l>>1)%f.fm);
        g.mulcoefficient(w);
        g+=v.mulcoefficient(f);
        f.mulcoefficient(u);
        h.mulcoefficient(w);
        f.resize(l+1),g.resize(l+1),h.resize(l+1);
    }
    ll res1=1,res2=0,res3=1;
    for (int i=0;i<siz;i++){
        res2=(res1*g.a[i]+res2*h.a[i])%f.fm;
        res1=res1*f.a[i]%f.fm;
        res3=res3*h.a[i]%f.fm;
    }
    for (int i=siz*B+1;i<=m;i++){
        res2=(res1+res2)*i%f.fm;
        res1=res1*(n-i+1)%f.fm;
        res3=res3*i%f.fm;
    }
    // 计算零散部分
    return qpow(res3,f.mod-2,f.fm,res2);
}

时间复杂度为 O(\sqrt n\log n)，可以轻松通过。

前面的多项式部分可以看这里。

更多题目：loj6386，U562050（我自己的题，需要任意模数）。