【BZOJ 3516&4126】国王奇遇记

problem

给定 $n,m$，求：

$$\sum _{i=1}^n{i}^m{m}^i$$

BZOJ 3516：$1\le n\le 10^9,1\le m\le 1000$。

BZOJ 4126：$1\le n\le 10^9,1\le m\le 5\times 10^5$。

---

solution

BZOJ 3516：

我们设 $f(i)={\sum }_{j=1}^n{j}^i{m}^j$，用微扰法化简可以得到：

$$\begin{array}{rl}(m-1)\times f(i)& =\sum _{j=1}^n{j}^i{m}^{j+1}-\sum _{j=1}^n{j}^i{m}^j\\ & =\sum _{j=1}^{n+1}(j-1{)}^i{m}^j-\sum _{j=1}^n{j}^i{m}^j\\ & =n^i{m}^{n+1}+\sum _{j=1}^n{m}^j\left(\sum _{k=0}^i\left(\genfrac{}{}{0ex}{}{i}{k}\right)j^k(-1{)}^{i-k}-j^i\right)\\ & =n^i{m}^{n+1}+\sum _{k=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i}{k}\right)(-1{)}^{i-k}\sum _{j=1}^n{j}^k{m}^j\\ & =n^i{m}^{n+1}+\sum _{k=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i}{k}\right)(-1{)}^{i-k}f(k)\end{array}$$

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

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

BZOJ 4126：

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

我们设 $S(n)={\sum }_{i=0}^{n-1}{i}^m{m}^i$，我们要求的就是 $S(n+1)$。

有一个结论：存在一个 $\le m$ 次的多项式 $g(x)$，使得 $S(n)=m^{n}g(n)-g(0)$。

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

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

$$m^{n}g(n)-m^{n-1}g(n-1)=m^{n-1}{i}^{n-1}$$

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

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

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

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

$$\sum _{i=0}^{k+1}(-1{)}^i(\genfrac{}{}{0ex}{}{k+1}{i})g(k+1-i)=0$$

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

时间复杂度 $O(m\log m)$。带个 $\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;
}