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

$$\begin{array}{rl}f(n)& =\sum _{i=1}^k{i}^k\\ g(n)& =\sum _{i=1}^{n}f(i)\\ h(n)& =\sum _{i=0}^{n}g(a+id)\end{array}$$

$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;
}