数论 差分+拉格朗日插值 P5430题解
Prean
·
·
题解
新的 O(k+\log n) 做法。
考虑计算每个猴子对答案的贡献。
打个表：
1 1 2 4 8 16 32 ...
可以看出第 i 个猴子对答案的贡献是 i^k \times 2^{n-i-1} ，特别地，最后一只猴子对答案的贡献是 n^k 。
写成柿子：
n^k+\sum_{i=1}^{n-1}i^k \times 2^{n-i-1}
n^k+2^{n-1} \times (\sum_{i=1}^{n-1}i^k \times (2^{-1})^i)
我们只需要计算出 \sum_{i=1}^{n-1}i^k \times (2^{-1})^i 即可。
然后我们发现这个柿子是 CODECHEF qpolysum，然后就做完了。
还是把做法写一遍吧
qpolysum 和本题有不同之处，即 i 从 0 开始，不过没啥区别，因为你要减掉的是一个 0 （
S(n)=\sum_{i=0}^{n-1}i^k \times m^i
在本题中相当于 m=2^{-1} 。
不过这个做法是猜了一个很奇怪的结论，并且做法来自校 OJ 讨论区（
我们猜 S(n)=m^n G(n) - G(0) ，其中 G(x) 是一个不超过 k 次的多项式。
证明可以看这个 blog 我才不告诉你是我看不懂
然后差分一下：
S(n)-S(n-1)=m^n G(n) - m^{n-1} G(n-1) = (n-1)^k \times m^{n-1}
G(n) = \frac {(n-1)^k + G(n-1)} m
设 G(0)=x ，那么我们在 n 为任何值的时候用 x 表示 G(n) 。
因为这个多项式的次数最高为 k ，而对一个 k 次多项式差分 k+1 次后为 0 ，所以我们把 G(x) 差分 k+1 次后得到：
\sum_{i=0}^{k+1} (-1)^{i+1} \binom {k+1} i G(k+1-i) = 0
我们可以用 x 表示 G(0) \sim G(k+1) ，然后解一个一元一次方程即可得到 x ，带入可得到 G(1) \sim G(k+1) 的值。
现在我们可以使用拉格朗日插值计算 G(n) 了，答案就是 m^nG(n)-G(0) 。
什么？你卡我空间？$ 10 $ 个数组太多了？？？
实际上可以缩成 $ 7 $ 个 $ \rm int $ 数组和 $ 1 $ 个 $ \rm bool $ 数组。
我的线性筛使用的是直接纪录最小质因数而非纪录是否为质数，可以改成后者。
然后就是 $ q $ 和 $ p $ 数组在使用的时候，$ x $ 和 $ y $ 数组已经不会使用了，所以直接使用 $ x $ 和 $ y $ 代替 $ q $ 和 $ p $ 即可。
code:
```cpp
#include<cstdio>
#include<cctype>
const int M=2e7+5,mod=1e9+7;
int k,n1,n2,top,x[M],y[M],idk[M],pri[M],fac[M],ifac[M];bool zhi[M];
int G[M];
inline void read(){
char s;long long t1,t2;
while(isdigit(s=getchar())){
t1=n1*10ll+(s^48);n1=t1>=mod?t1%mod:t1;
t2=n2*10ll+(s^48);n2=t2>=mod-1?t2%(mod-1):t2;
}
}
inline int Add(const int&a,const int&b){
return a+b>=mod?a+b-mod:a+b;
}
inline int Del(const int&a,const int&b){
return b>a?a-b+mod:a-b;
}
inline int C(const int&n,const int&m){
return 1ll*fac[n]*ifac[m]%mod*ifac[n-m]%mod;
}
inline int pow(int a,int b){
int ans=1;
for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;
return ans;
}
inline void sieve(const int&M){
register int i,j,x;idk[1]=1;
for(i=2;i<=M;++i){
if(!zhi[i])pri[++top]=i,idk[i]=pow(i,k);
for(j=1;j<=top&&(x=i*pri[j])<=M;++j){
idk[x]=1ll*idk[i]*idk[pri[j]]%mod;zhi[x]=true;
if(!(i%pri[j]))break;
}
}
}
inline int Inter(const int&n){
register int i,tmp,ans=0;
x[0]=y[k+2]=1;
for(i=1;i<=k+1;++i)x[i]=1ll*x[i-1]*Del(n,i)%mod;
for(i=k+1;i>=1;--i)y[i]=1ll*y[i+1]*Del(n,i)%mod;
for(i=1;i<=k+1;++i){
if(k+1-i&1)ans=Del(ans,1ll*1ll*x[i-1]*y[i+1]%mod*G[i]%mod*ifac[i-1]%mod*ifac[k+1-i]%mod);
else ans=Add(ans,1ll*1ll*x[i-1]*y[i+1]%mod*G[i]%mod*ifac[i-1]%mod*ifac[k+1-i]%mod);
}
return ans;
}
signed main(){
register int i,X=0,Y=0;
fac[0]=fac[1]=ifac[0]=ifac[1]=1;read();scanf("%d",&k);sieve(k+1);x[0]=1;y[0]=0;
for(i=1;i<=k+1;++i)x[i]=Add(x[i-1],x[i-1]),y[i]=Add(y[i-1],idk[i-1]),y[i]=Add(y[i],y[i]);
for(i=2;i<=k+1;++i)fac[i]=1ll*fac[i-1]*i%mod,ifac[i]=1ll*(mod-mod/i)*ifac[mod%i]%mod;
for(i=2;i<=k+1;++i)ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;
for(i=0;i<=k+1;++i){
if(i&1){
X=Add(X,1ll*C(k+1,i)*x[k+1-i]%mod);
Y=Add(Y,1ll*C(k+1,i)*y[k+1-i]%mod);
}
else{
X=Del(X,1ll*C(k+1,i)*x[k+1-i]%mod);
Y=Del(Y,1ll*C(k+1,i)*y[k+1-i]%mod);
}
}
G[0]=mod-1ll*Y*pow(X,mod-2)%mod;
for(i=1;i<=k+1;++i)G[i]=Add(1ll*x[i]*G[0]%mod,y[i]);
printf("%d",Add(1ll*Del(1ll*pow(500000004,n2)*Inter(n1)%mod,G[0])*pow(2,n2-1)%mod,pow(n1,k)));
}
```