[BZOJ3453][tyvj1858]XLkxc-拉格朗日插值

XLkxc

Description

神犇LYD虐完HEOI之后给傻×XLk出了一题：

SHY是某国的公主，平时的一大爱好是读诗…(中间略)…结果mod p就可以了

简明题意

给定 k,a,n,d,p

f(i)=1^k+2^k+3^k+……+i^k

g(x)=f(1)+f(2)+f(3)+….+f(x)

求(g(a)+g(a+d)+g(a+2d)+……+g(a+nd))mod p

对于所有数据

1<=k<=123

0<=a,n,d<=123456789

p==1234567891

Input

第一行数据组数，(保证小于6)

以下每行四个整数 k,a,n,d

Output

每行一个结果。

Sample Input

5

1 1 1 1

1 1 1 1

1 1 1 1

1 1 1 1

1 1 1 1

Sample Output

5

5

5

5

5

这个样例啊，亦可赛艇。

思路:

可以发现,f(x)f(x)就是一个自然数幂和。

那么通过差分可知，f(x)f(x)为一个k+1k+1次多项式。

同理根据差分可得g(x)g(x)为一个k+2k+2次多项式。

令h(x)=g(a)+g(a+d)+g(a+2∗d)+......+g(a+n∗d)h(x)=g(a)+g(a+d)+g(a+2∗d)+......+g(a+n∗d)

那么同样根据差分可以得到h(x)h(x)为一个k+3k+3次多项式。

于是考虑拉格朗日插值。

首先预处理处f(x)f(x)，从而算出g(x)g(x)。

然后，插出g(x)g(x)以计算h(x)h(x)的前k+4k+4项。

最后，插出h(x)h(x)并计算答案。

然后就没有了……

记录一下重心拉格朗日插值的公式:

l(x)=∏i=1n(x−xi)l(x)=∏i=1n(x−xi)

wi=1∏j≠i(xi−xj)wi=1∏j≠i(xi−xj)

f(x)=l(x)∑i=1nwix−xif(x)=l(x)∑i=1nwix−xi

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

typedef long long ll;
const ll md=1234567891;
const ll N=139;

ll k,a,n,d;
ll g
,f
,fac
,inv
,tl
,w
;

inline ll qpow(ll a,ll b)
{
ll ret=1;
while(b)
{
if(b&1)ret=ret*a%md;
a=a*a%md;b>>=1;
}
return ret;
}

inline void init()
{
fac[0]=1;
for(ll i=1;i<N;i++)
fac[i]=fac[i-1]*i%md;
inv[N-1]=qpow(fac[N-1],md-2);
for(ll i=N-1;i>=1;i--)
inv[i-1]=inv[i]*i%md;
}

inline ll cg(ll x,ll m,ll *y)
{
if(x<=m)return y[x];
ll l=1,sum=0;
for(ll i=1;i<=m;i++)
l=l*(x-i)%md;
for(ll i=1;i<=m;i++)
(sum+=w[i]*qpow(x-i,md-2)%md*y[i]%md)%=md;
return sum*l%md;
}

inline ll initw(ll m)
{
for(ll i=1;i<=m;i++)
w[i]=((((m-i)&1)?-1ll:1ll)*inv[m-i]*inv[i-1]%md+md)%md;
}

inline void mina()
{
scanf("%lld%lld%lld%lld",&k,&a,&n,&d);
for(ll i=0;i<=k+3;i++)
g[i]=qpow(i,k);
for(ll i=1;i<=k+3;i++)
(g[i]+=g[i-1])%=md;
for(ll i=1;i<=k+3;i++)
(g[i]+=g[i-1])%=md;

initw(k+3);
for(ll i=0;i<=k+4;i++)
f[i]=cg((a+i*d)%md,k+3,g);
for(ll i=1;i<=k+4;i++)
f[i]=(f[i]+f[i-1])%md;
initw(k+4);
printf("%lld\n",cg(n,k+4,f));
}

int main()
{
init();
int T;scanf("%d",&T);
while(T--)
mina();

return 0;
}