【BZOJ 35164126】国王奇遇记

problem

给定 n , m n,m n,m，求：

∑ i = 1 n i m m i \sum_{i=1}^ni^mm^i i=1∑nimmi

BZOJ 3516： 1 ≤ n ≤ 1 0 9 , 1 ≤ m ≤ 1000 1\le n\le10^9,1\le m\le1000 1≤n≤109,1≤m≤1000。

BZOJ 4126： 1 ≤ n ≤ 1 0 9 , 1 ≤ m ≤ 5 × 1 0 5 1\le n\le10^9,1\le m\le5\times 10^5 1≤n≤109,1≤m≤5×105。

solution

BZOJ 3516：

我们设 f ( i ) = ∑ j = 1 n j i m j f(i)=\sum_{j=1}^nj^im^j f(i)=∑j=1njimj，用微扰法化简可以得到：

( m − 1 ) × f ( i ) = ∑ j = 1 n j i m j + 1 − ∑ j = 1 n j i m j = ∑ j = 1 n + 1 ( j − 1 ) i m j − ∑ j = 1 n j i m j = n i m n + 1 + ∑ j = 1 n m j ( ∑ k = 0 i ( i k ) j k ( − 1 ) i − k − j i ) = n i m n + 1 + ∑ k = 0 i − 1 ( i k ) ( − 1 ) i − k ∑ j = 1 n j k m j = n i m n + 1 + ∑ k = 0 i − 1 ( i k ) ( − 1 ) i − k f ( k ) \begin{aligned} (m-1)\times f(i)&=\sum_{j=1}^nj^im^{j+1}-\sum_{j=1}^nj^im^j\\ &=\sum_{j=1}^{n+1}(j-1)^im^j-\sum_{j=1}^nj^im^j\\ &=n^im^{n+1}+\sum_{j=1}^nm^j\left(\sum_{k=0}^{i}\binom i kj^k(-1)^{i-k}-j^i\right)\\ &=n^im^{n+1}+\sum_{k=0}^{i-1}\binom i k(-1)^{i-k}\sum_{j=1}^nj^km^j\\ &=n^im^{n+1}+\sum_{k=0}^{i-1}\binom i k(-1)^{i-k}f(k) \end{aligned} (m−1)×f(i)=j=1∑njimj+1−j=1∑njimj=j=1∑n+1(j−1)imj−j=1∑njimj=nimn+1+j=1∑nmj(k=0∑i(ki)jk(−1)i−k−ji)=nimn+1+k=0∑i−1(ki)(−1)i−kj=1∑njkmj=nimn+1+k=0∑i−1(ki)(−1)i−kf(k)

于是我们直接 O ( m 2 ) O(m^2) O(m2) d p dp dp 就可以了，答案就是 f ( m ) f(m) f(m)。

注意要特判 m = 1 m=1 m=1 的情况。

BZOJ 4126：

PS：这部分内容和上面没有什么关系。

我们设 S ( n ) = ∑ i = 0 n − 1 i m m i S(n)=\sum_{i=0}^{n-1}i^mm^i S(n)=∑i=0n−1immi，我们要求的就是 S ( n + 1 ) S(n+1) S(n+1)。

有一个结论：存在一个 ≤ m ≤m ≤m 次的多项式 g ( x ) g(x) g(x)，使得 S ( n ) = m n g ( n ) − g ( 0 ) S(n)=m^ng(n)-g(0) S(n)=mng(n)−g(0)。

证明可以参考这篇博客，我就不赘述了。

现在考虑如何求 g ( x ) g(x) g(x)。将 S ( n ) S(n) S(n) 与 S ( n − 1 ) S(n-1) S(n−1) 相减，得：

m n g ( n ) − m n − 1 g ( n − 1 ) = m n − 1 i n − 1 m^ng(n)-m^{n-1}g(n-1)=m^{n-1}i^{n-1} mng(n)−mn−1g(n−1)=mn−1in−1

所以解出来得到 g ( n ) = g ( n − 1 ) + i n − 1 m g(n)=\frac{g(n-1)+i^{n-1}}{m} g(n)=mg(n−1)+in−1。

所以如果我们算出了 g ( 0 ) g(0) g(0)，就可以不断迭代得到 g ( 1 ) , g ( 2 ) , . . . , g ( n ) g(1),g(2),...,g(n) g(1),g(2),...,g(n)。

不妨设 g ( 0 ) = x g(0)=x g(0)=x，那么其实 g ( 1 ) , g ( 2 ) , . . . , g ( n ) g(1),g(2),...,g(n) g(1),g(2),...,g(n) 可以看作是关于 x x x 的一次函数，而且可以把每个的 k k k 和 b b b 算出来。

然看考虑高阶差分的公式（不会的可以参考百度）：

∑ i = 0 k + 1 ( − 1 ) i ( k + 1 i ) g ( k + 1 − i ) = 0 \sum_{i=0}^{k+1}(-1)^i\binom{k+1}{i}g(k+1-i)=0 i=0∑k+1(−1)i(ik+1)g(k+1−i)=0

我们就可以解出 g ( 0 ) g(0) g(0) 了，之后再用拉格朗日插值求 g ( n + 1 ) g(n+1) g(n+1) 就可以求出 S ( n + 1 ) S(n+1) S(n+1) 了。

时间复杂度 O ( m log ⁡ m ) O(m\log m) O(mlogm)。带个 log ⁡ \log log 是因为有快速幂。

code

BZOJ 3516：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 1005
#define P 1000000007
using namespace std;
int n,m,Pow[N],fac[N],ifac[N],f[N];
int add(int x,int y)  {return x+y>=P?x+y-P:x+y;}
int dec(int x,int y)  {return x-y< 0?x-y+P:x-y;}
int mul(int x,int y)  {return 1ll*x*y>=P?1ll*x*y%P:x*y;}
int power(int a,int b,int ans=1){for(;b;b>>=1,a=mul(a,a))if(b&1)  ans=mul(ans,a);return ans;
}
void init(){Pow[0]=1,fac[0]=fac[1]=1;for(int i=1;i<N;++i)  Pow[i]=mul(Pow[i-1],n);for(int i=2;i<N;++i)  fac[i]=mul(fac[i-1],i);ifac[N-1]=power(fac[N-1],P-2);for(int i=N-2;~i;--i)  ifac[i]=mul(ifac[i+1],i+1);
}
int C(int n,int m)  {return mul(fac[n],mul(ifac[m],ifac[n-m]));}
int main(){scanf("%d%d",&n,&m),init();if(m==1)  return printf("%d",(1ll*n*(n+1)/2)%P),0;int val=power(m,n+1),inv=power(m-1,P-2);f[0]=mul(dec(val,m),inv);for(int i=1;i<=m;++i){for(int k=0;k<i;++k)f[i]=(i-k)&1?dec(f[i],mul(C(i,k),f[k])):add(f[i],mul(C(i,k),f[k]));f[i]=add(f[i],mul(Pow[i],val)),f[i]=mul(f[i],inv);}printf("%d\n",f[m]);return 0;
}

BZOJ 4126：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 500005
#define P 1000000007
using namespace std;
int n,m,up,fac[N],ifac[N],f[N],g[N],k[N],b[N];
int add(int x,int y)  {return x+y>=P?x+y-P:x+y;}
int dec(int x,int y)  {return x-y< 0?x-y+P:x-y;}
int mul(int x,int y)  {return 1ll*x*y%P;}
int power(int a,int b,int ans=1){for(;b;b>>=1,a=mul(a,a))if(b&1)  ans=mul(ans,a);return ans;
}
void init(){up=m+1,fac[0]=fac[1]=1;for(int i=2;i<=up;++i)  fac[i]=mul(fac[i-1],i);ifac[up]=power(fac[up],P-2);for(int i=up-1;~i;--i)  ifac[i]=mul(ifac[i+1],i+1);
}
int C(int n,int m)  {return mul(fac[n],mul(ifac[m],ifac[n-m]));}
int Pre[N],Suf[N];
int lagrange(int x){if(x<=up)  return g[x];int ans=0;Pre[0]=1,Suf[up+1]=1;for(int i=1;i<=up;++i)  Pre[i]=mul(Pre[i-1],x-i);for(int i=up;i>=1;--i)  Suf[i]=mul(Suf[i+1],x-i);for(int i=1;i<=up;++i){int temp=mul(g[i],mul(Pre[i-1],Suf[i+1]));temp=mul(temp,mul(ifac[i-1],(up-i)&1?(P-ifac[up-i]):ifac[up-i]));ans=add(ans,temp);}return ans;
}
int main(){scanf("%d%d",&n,&m),init();if(m==1)  return printf("%d\n",(1ll*n*(n+1)/2)%P),0;for(int i=1;i<=m;++i)  f[i]=power(i,m);k[0]=1,b[0]=0;int inv=power(m,P-2),K=0,B=0;for(int i=1;i<=up;++i)  k[i]=mul(k[i-1],inv),b[i]=mul(add(b[i-1],f[i-1]),inv);for(int i=0;i<=up;++i){K=add(K,mul((i&1)?1:P-1,mul(C(m+1,i),k[m+1-i])));B=add(B,mul((i&1)?1:P-1,mul(C(m+1,i),b[m+1-i])));}g[0]=mul(P-B,power(K,P-2));for(int i=1;i<=up;++i)  g[i]=mul(add(g[i-1],f[i-1]),inv);printf("%d\n",dec(mul(power(m,n+1),lagrange(n+1)),g[0]));return 0;
}