这个科技是用来 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+1 和 n+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 小练习