【BZOJ3453】XLkxc（拉格朗日插值）

\[\begin{aligned} f(n)&=\sum_{i=1}^k i^k\\ g(n)&=\sum_{i=1}^n f(i)\\ h(n)&=\sum_{i=0}^ng(a+id) \end{aligned} \]

\(f(n)\) 是自然数幂求和，为 \(k+1\) 次多项式。

\(g(n)\) 的差分 \(g(n)-g(n-1)=f(n)\) 为 \(k+1\) 次多项式，故 \(g(n)\) 为 \(k+2\) 次多项式。

\(h(n)\) 的差分：\(h(n)-h(n-1)=g(a+nd)\) 为 \(k+2\) 次多项式，故 \(h(n)\) 为 \(k+3\) 次多项式。

所以找 \(k+4\) 个值代入 \(h(n)\) 再拉格朗日插值即可。

取 \(n=1,2,\cdots,k+4\) 代入 \(h(n)\) 的话，就需要 \(g(a),g(a+d),\cdots,g(a+(k+4)d)\)。

于是我们也要对 \(g\) 进行 \(k+4\) 次插值（

时间复杂度 \(O(k^2)\)。

#include<bits/stdc++.h> 

#define K 150

using namespace std;

namespace modular
{
	const int mod=1234567891;
	inline int add(const int x,const int y){return 0ll+x+y>=mod?0ll+x+y-mod:x+y;}
	inline int dec(const int x,const int y){return 0ll+x-y<0?0ll+x-y+mod:x-y;}
	inline int mul(const int x,const int y){return 1ll*x*y%mod;}
}using namespace modular;

inline int poww(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

int T,k,a,n,d;
int fac[K],ifac[K];
int f[K],g[K],h[K];

void init()
{
	const int maxn=130;
	fac[0]=1;
	for(int i=1;i<=maxn;i++) fac[i]=mul(fac[i-1],i);
	ifac[maxn]=poww(fac[maxn],mod-2);
	for(int i=maxn;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
}

int calc(int k,int n,int *y)
{
	static int pre[K],suf[K];
	pre[0]=1; 
	for(int i=1;i<=k;i++) pre[i]=mul(pre[i-1],dec(n,i));
	suf[k+1]=1;
	for(int i=k;i>=1;i--) suf[i]=mul(suf[i+1],dec(n,i));
	int ans=0;
	for(int i=1;i<=k;i++)
		ans=add(ans,mul(y[i],mul(mul(pre[i-1],suf[i+1]),mul(((k-i)&1)?(mod-1):1,mul(ifac[i-1],ifac[k-i])))));
	return ans;
}

int main()
{
	T=read();
	init();
	while(T--)
	{
		k=read(),a=read(),n=read(),d=read();
		for(int i=1;i<=k+3;i++) f[i]=poww(i,k);
		for(int i=1;i<=k+3;i++) f[i]=add(f[i-1],f[i]);
		for(int i=1;i<=k+3;i++) g[i]=add(g[i-1],f[i]);
		for(int i=0;i<=k+4;i++) h[i]=calc(k+3,add(a,mul(i,d)),g);
		for(int i=1;i<=k+4;i++) h[i]=add(h[i-1],h[i]);
		printf("%d\n",calc(k+4,n,h));
	}
	return 0;
}