BZOJ 3157/3516 国王奇遇记 [数学]

题意

[sum_{i=0}^n i ^ m imes m ^ i ]

其中 (1 leqslant n leqslant 10^9) ， (1 leqslant m leqslant N) ，其中的 (N) 是因为这个题有三个版本

[egin{align*} &bzoj3157 - N leqslant 200 \ &bzoj3516 - N leqslant 1000 \ &b bzoj4126 - N leqslant 500000 \ end{align*} ]

但我只会 (O(m^2)) 的，所以先不提第三题（挖坑待填）

分析

发现 (n) 很大，不能直接枚举 (i) ，考虑从 (m) 入手 递推

记 (f_{n,a} = sum_{i=0}^n i^a m^i) （这里 (i) 从 (0) 开始方便后面转换，(i =0) 对答案没有影响），题目要求 (f_{n,m}) ，考虑转移

[egin{align*} f_{n,a} &= sum_{i=0}^n i^a m^i \ &= m sum_{i=0}^n i^a m^{i-1} //提出一个m \ &= m sum_{i=0}^{n-1} (i+1)^a m^i //用 i 表示 i+1 \ end{align*} ]

当 (i = -1) 时 ((i+1)^a =0) 大力不理

再根据二项式定理 ((a+b)^n = sum_{i=0}^n C_n^i a^i b^{n-i})

[egin{align*} &= m sum_{i=0}^{n-1} m^i sum_{k=0}^a C_a^k i^k \ &= m sum_{k=0}^a C_a^k sum_{i=0}^{n-1} i^k m^i //交换枚举顺序 \ &= m sum_{k=0}^a C_a^k f_{n-1,k} \ &= m sum_{k=0}^a C_a^k (f_{n,k} - n^k m^n) end{align*} ]

得到了(O(m^2)) 的转移方程，完结撒花

好像 (O(m)) 的要多项式插值和快速傅里叶变换求解，不会，不想填坑了

等等！完结不了，注意：

1、上式的 (k) 有取到 (a) 的情况，解个关于 (f_{n,i}) 方程就好了

2、(f[0] = sum_{i=0}^n m^i) 等比数列，先一算，特判 (m=1) 时不能用等比数列求和公式，直接 (ans = sum_{i=1}^n i = frac{n imes (n+1)}{2})

代码

#include<bits/stdc++.h>
using namespace std;

namespace TYC
{
	typedef long long ll;

	const ll p=1e9+7;
	const int N=1010;

	int n,m;
	ll f[N],Pow[N],C[N][N];

	inline ll qpow(ll x,ll tim)
	{
		ll ans=1;
		for(;tim;tim>>=1,x=x*x%p)
			if(tim&1) ans=ans*x%p;
		return ans;
	}

	#define Mod(x) (x)>=p?(x)-p:(x)

	void init()
	{
		C[0][0]=1;
		for(int i=1;i<=m;i++)
		{
			C[i][0]=1;
			for(int j=1;j<=m;j++)
				C[i][j]=Mod(C[i-1][j-1]+C[i-1][j]);
		}
		Pow[0]=1;
		for(int i=1;i<=m;i++) Pow[i]=Pow[i-1]*n%p;
	}

	void work()
	{
		scanf("%d%d",&n,&m);
		if(m==1)
		{
			printf("%lld
",(ll)(n+1)*n/2%p);
			return;
		}
		init();
		ll t1=qpow(m,n),t2=t1*m%p,inv=qpow(m-1,p-2);
		f[0]=(Mod(qpow(m,n+1)-1+p))*qpow(m-1,p-2)%p;
		for(int i=1;i<=m;i++)
		{
			ll sum=0;
			for(int j=0;j<i;j++)
				sum=Mod(sum+(C[i][j]*(f[j]-Pow[j]*t1%p+p)%p));
			if(m!=1) f[i]=(Pow[i]*t2-m*sum+p)%p*inv%p;
		}
		printf("%lld
",f[m]);
	}
}

int main()
{
	TYC::work();
	return 0;
}