题意即求$$sum_{i=0}^nsum_{j=1}^{a+id}sum_{x=1}^jx^k$$
我们知道最后一个(sum)是自然数幂和，设(f(n)=sum_{x=1}^nx^k)，这是一个(k+1)次多项式，可以插值求出（当然本题只需要求出任意(k+3)个值即可不需要插值）。
令(g(n)=sum_{i=1}^nf(i))，(打表)差分可知这是一个(k+2)次多项式。
同样令(h(n)=sum_{i=0}^ng(a+id))，同样差分可知这是一个(k+3)次多项式。
所以用拉格朗日插值我们代入(k+4)个值就可以求出(h(n))了。
也就是先求(k+3)个(f(x))的值，再求出(k+3)个(g(x))的值，然后对(g)插值求(k+4)个(g(a+xd))，前缀和一下就有了(h)的(k+4)个值，然后再插一次就得到(h(n))了。（嵌套好鬼畜...）
注意(f,g,h)都是个前缀和...
注意这题两倍的模数会爆int。可以用unsigned int。
从这题可以看出：
- 多一个(sum)一般会使多项式次数+1。
- 插值可以嵌套，且复杂度不变，仍是(O(k^2))。
这是伯努利数的做法：https://blog.csdn.net/qq_20669971/article/details/65938763。
//824kb 28ms
#include <cstdio>
#include <algorithm>
#define mod 1234567891
#define Mod(x) x>=mod&&(x-=mod)
#define Add(x,v) (x+=v)>=mod&&(x-=mod)
typedef long long LL;
typedef unsigned int uint;
const int N=130;
uint g[N],h[N],ifac[N];
inline uint FP(uint x,uint k)
{
uint t=1;
for(; k; k>>=1,x=1ll*x*x%mod)
if(k&1) t=1ll*t*x%mod;
return t;
}
uint Lagrange(uint *y,const int m,uint x)
{
static uint pre[N],suf[N];
pre[0]=x, suf[m+1]=1;
for(int i=1; i<m; ++i) pre[i]=1ll*pre[i-1]*(x+mod-i)%mod;
for(int i=m; i; --i) suf[i]=1ll*suf[i+1]*(x+mod-i)%mod;
LL ans=0;
for(int i=0,up,down; i<=m; ++i)
{
if(i) up=1ll*pre[i-1]*suf[i+1]%mod*y[i]%mod;
else up=1ll*suf[i+1]*y[i]%mod;
down=(m-i)&1?mod-1ll*ifac[i]*ifac[m-i]%mod:1ll*ifac[i]*ifac[m-i]%mod;
ans+=1ll*up*down%mod;
}
return ans%mod;
}
int main()
{
ifac[N-1]=1119688141;//是129!的逆元不是129的!!!今天错了两次这个真是醉了=-=
for(int i=N-1; i; --i) ifac[i-1]=1ll*ifac[i]*i%mod;
int T,K; uint a,n,d; g[0]=0;
for(scanf("%d",&T); T--; )
{
scanf("%d%u%u%u",&K,&a,&n,&d);
for(int i=1; i<=K+2; ++i) g[i]=g[i-1]+FP(i,K), Mod(g[i]); //f = sum i^k
for(int i=1; i<=K+2; ++i) Add(g[i],g[i-1]); //g = sum f(i)
h[0]=Lagrange(g,K+2,a);
for(int i=1; i<=K+3; ++i) Add(a,d), h[i]=h[i-1]+Lagrange(g,K+2,a), Mod(h[i]);
printf("%d
",(int)Lagrange(h,K+3,n));
}
return 0;
}