把EI科技 【载谈 Binomial Sum】 用人话说出来

Prean

2021-09-13 13:29:03

Algo. & Theory

这个科技是用来 O(k+\log n)\sum_{i=0}^n[x^i]f(x)p^i(x)\mod x^{k+1} 这个多项式的某些项数的线性组合,其中 f(x) 微分有限。

EI 原文讲得过于抽象和简洁,需要有比较硬的基础才能看懂,所以这篇文章只会大概走一遍这个玩意儿的流程。。。毕竟本来就不是什么简单东西。

微分有限

对于 f(x),称其微分有限仅当存在多项式 Q_0(x),Q_1(x),...Q_m(x) 满足:

\sum_{i=0}^{m}f^{(i)}Q_i(x)=0

要求 Q_i(x) 的项数有限。人话就是存在一个关于 f(x) 的 ODE

例题1 CF932E

求:

\sum_{i=0}^n\binom{n}{i}i^k

\binom{n}{i} 看成 [x^i] (1+x)^n , i^k 看成 [\frac{x^k}{k!}] e^{ix},那么原式变成了:

[\frac{x^k}{k!}]\sum_{i=0}^n[x^i](1+x)^ne^{ix}

相当于

[\frac{x^k}{k!}](e^x+1)^n

当然你也可以将 i^k 直接看成 [\frac{x^k}{k!}]e^{ix} 然后使用二项式定理得到上式。

这里已经可以保留 (e^x+1)\mod x^{k+1},通过多项式快速幂 O(k\log k\log n) 计算答案了,但是我们显然不满足于此。

接下来是重要的一步:

F(z)=(z+1)^n,\mathcal F(z+1)=F(z+1)\mod z^{k+1}

为什么这里是 z+1 而不是 z 呢?

实际上是因为这里的 z=e^x-1常数项为 0

为什么常数项一定要为 0

考虑原式,你算的是 \sum_{i=0}^n[x^i]f(x)p^i(x),答案只截取到第 k 项,那么如果 [x^0]p(x)=0,我们就可以只求和到 k

更严谨的解释是把这个东西在 1 处泰勒展开会变成有穷求和。

换而言之,高次项对答案没有贡献。

所以在这个算法中,在这一步中要设的是 \mathcal F(z+G(0))= F(z+G(0))\mod x^{k+1}

因为 \mathcal F(z+1)F(z+1) 的前 k+1 次项,所以答案 [\frac{x^k}{k!}]F(e^x)=[\frac{x^k}{k!}]\mathcal F(e^x)

问题来了,我们如何得到 \mathcal F(z+1)

我们对 F(z) 列出一个微分方程:

(z+1)F'(z)-nF(z)=0

那么对于 \mathcal F(z+1),呢?

(z+2)\mathcal F'(z+1)-n\mathcal F(z+1)=0

好像不对?

是不是多算了什么东西?

考虑到原本不应该有的 [z^k](z+2)F'(z+1) 突然出现,所以应该在右边加上这玩意儿。

(z+2)\mathcal F'(z+1)-n\mathcal F(z+1)=[z^{k+1}](z+1)F'(z)=(k-n)\binom{n}{k}z^k 2^{n-k}

对这个提取系数就可以 O(k) 递推 \mathcal F(z+1) 了。

提取系数后得到:

\mathcal F[i]=\frac{(k-n)\binom{n}{k}\binom{k}{i-1}(-1)^{k-i+1}2^{n-k}+(n-i+1)\mathcal F[i-1]}i

不过要注意 \mathcal F[0] 的值,咱还没有算。。。

注意 \mathcal F[0] 的定义是 \mathcal F(z+1)0 次项,也就是:

[z^0]\sum_{i=0}^k[z^i]\mathcal F(z+1)(z-1)^i

还记得 \mathcal F(z+1)=F(z+1)\mod x^{k+1} 吧?那么 [z^i]\mathcal F(z+1) 就可以使用二项式定理计算了。

\sum_{i=0}^k([z^0](z-1)^i)\times([z^i](z+2)^n) \sum_{i=0}^k(-1)^i\times\binom{n}{i}\times 2^{n-i}

回到原式:

[\frac{x^k}{k!}]\mathcal F(e^x)=\sum_{i=0}^k e^{ix}[x^i]\mathcal F(x)

也就是:

\sum_{i=0}^k i^k[x^i]\mathcal F(x)

线性筛 id^k 就可以做到 O(k+\log n) 啦。

#include<cstdio>
typedef unsigned uint;
const uint M=5005,mod=1e9+7;
uint n,k,top,pri[M],pos[M],idk[M];uint F[M],p2[M],Ck[M],Cn[M],inv[M];
inline uint Add(const uint&a,const uint&b){
    return a+b>=mod?a+b-mod:a+b;
}
inline uint Del(const uint&a,const uint&b){
    return b>a?a-b+mod:a-b;
}
inline uint pow(uint a,uint b){
    uint ans=1;
    for(;b;b>>=1,a=1ull*a*a%mod)if(b&1)ans=1ull*ans*a%mod;
    return ans;
}
inline void sieve(const uint&M){
    register uint i,j,x;idk[1]=1;
    for(i=2;i<=M;++i){
        if(!pos[i])idk[pri[pos[i]=++top]=i]=pow(i,k);
        for(j=1;j<=top&&(x=i*pri[j])<=M;++j){
            idk[x]=1ull*idk[i]*idk[pri[j]]%mod;
            if((pos[x]=j)==pos[i])break;
        }
    }
}
signed main(){
    register uint i,x,ans=0;inv[0]=1;inv[1]=1;
    scanf("%u%u",&n,&k);Cn[0]=Ck[0]=1;Ck[1]=k;Cn[1]=n;sieve(k);
    for(i=2;i<=k;++i){
        inv[i]=1ull*(mod-mod/i)*inv[mod%i]%mod;
        Cn[i]=1ull*Cn[i-1]*(n-i+1)%mod*inv[i]%mod;
        Ck[i]=1ull*Ck[i-1]*(k-i+1)%mod*inv[i]%mod;
    }
    if(n<=k){
        for(i=1;i<=n;++i)ans=Add(ans,1ull*Cn[i]*idk[i]%mod);
        return!printf("%u",ans);
    }
    x=pow(2,n-k);
    for(i=k;i<=k;--i){
        if(i&1)F[0]=Del(F[0],1ull*Cn[i]*x%mod);
        else F[0]=Add(F[0],1ull*Cn[i]*x%mod);x=Add(x,x);
    }
    x=1ull*(n-k)*Cn[k]%mod*pow(2,n-k)%mod;
    for(i=1;i<=k;++i){
        if(k-i&1)F[i]=Del(1ull*(n-i+1)*F[i-1]%mod,1ull*x*Ck[i-1]%mod);
        else F[i]=Add(1ull*(n-i+1)*F[i-1]%mod,1ull*x*Ck[i-1]%mod);
        F[i]=1ull*F[i]*inv[i]%mod;
    }
    for(i=1;i<=k;++i)ans=Add(ans,1ull*idk[i]*F[i]%mod);
    printf("%u",ans);
}

例题2 CF392C

求:

\sum_{i=0}^n fib_i\times i^k

和上面一样:

[\frac{x^k}{k!}]\sum_{i=0}^n fib_i\times e^{ix}

F(x)=\sum_{i=0}^{n}fib_i\times x^i=\frac{1}{1-x-x^2}\mod x^{n+1}

[\frac{x^k}{k!}]F(e^x)

如果不截取的话有:

F(z)=zF(z)+z^2F(z)+1

截取后,右边多了 n+1n+2 次项的贡献,应该减去。也就是:

F(z)=zF(z)+z^2F(z)+1-(fib_{n-1}+fib_n)z^{n+1}-fib_nz^{n+2} F(z)=\frac {1-(fib_{n-1}+fib_n)z^{n+1}-fib_nz^{n+2}}{1-z-z^2}

因为需要求 F(e^x),所以还是设 \mathcal F(z+1)=F(z+1) \mod z^{k+1}

再设一个 G(z)=1-(fib_{n-1}+fib_n)z^{n+1}-fib_nz^{n+2}

F(z+1)=-\frac {G(z+1)} {1+3z+z^2} F(z+1)=-(G(z+1)+3zF(z+1)+z^2F(z+1))

原本在这里可以提取系数直接递推的,但是我们发现我们不会求 G(z+1)。。。

继续截取,再设一个 \mathcal G(z+1)=G(z+1) \mod z^{k+1}

考虑一个:

H_n(z+1)=(z+1)^n \mod x^{k+1}

有:

\mathcal G(z+1)=1-(fib_{n-1}+fib_n)H_{n+1}(z+1)-fib_nH_{n+2}(z+1)

所以只需要将 H_{n+1}(z+1)H_{n+2}(z+1) 递推出来就可以递推 \mathcal G(z+1)

如果你不理解这里的过程,回忆一下如果使用多项式直接计算 \mathcal G(z+1) 时是如何计算的。(多项式在计算的过程中会自动截取)

利用求导有:

(z+1)((z+1)^n)'=n((z+1)^n) (z+1)H_n'(z+1)=nH_n(z+1)+(n-k)\binom n kz^k

于是可以递推 H_{n+1}(z+1),H_{n+2}(z+1)\mathcal G(z+1)

回代:

[z^i]\mathcal F(z+1)=[i=k+1](3([z^k]\mathcal F(z+1))+([z^{k-1}]\mathcal F(z+1)))+[i=k+2]([z^k]\mathcal F(z+1)) -([z^i]\mathcal G(z+1)+3[z^{i-1}]\mathcal F(z+1)+[z^{i-2}]\mathcal F(z+1))

感觉好麻烦啊。。。

为了方便,设 a=3([z^k]\mathcal F(z+1))+([z^{k-1}]\mathcal F(z+1)),b=[z^k]\mathcal F(z+1)

将上式化简:

\mathcal F(z+1)=\frac{az^{k+1}+bz^{k+2}-\mathcal G(z+1)}{1+3z+z^2}

也就是:

\mathcal F(z)=\frac{a(z-1)^{k+1}+b(z-1)^{k+2}-\mathcal G(z)}{1-z-z^2}

这下可以递推 \mathcal F(z) 啦。

回代的结果和例题1一样,都是 \sum_{i=0}^ki^k[x^i]\mathcal F(x)

总结

计算 F(G(x)) 时,设:

\mathcal F(z+G(0))=F(z+G(0)) \mod z^{k+1}

然后通过求导列出关于 F(z+G(0)) 的方程,将 F 替换为 \mathcal F 后加上/减去被少算/多算的项,然后提取系数继续递推。

参考资料:

Elegia 「实验性讲稿」载谭 Binomial Sums 详解

Elegia 载谭 Binomial Sum:多项式复合、插值与泰勒展开

GuidingStar CF932E题解

GuidingStar 载谭 Binomial Sum 小练习