【题意】给定k<=123，a,n,d<=10^9，求：
$$f(n)=\sum_{i=0}^{n}\sum_{j=1}^{a+id}\sum_{x=1}^{j}x^k$$
【算法】拉格朗日插值
【题解】参考：拉格朗日插值法及应用 by DZYO
虽然式子很复杂，但一点一点化简有条理的化简后就可以做了。
首先最后是一个自然数幂和：
$$\sum_{x=1}^{j}x^k$$
这是一个k+1次多项式，可以理解为k+一个Σ（一般一个Σ增加一次项）。
然后会发现最后部分和第二部分之间不需要插值，因为第二部分的前若干小项的计算只需要第一部分的前若干小项，那么：
$$g(m)=\sum_{j=1}^{m}\sum_{x=1}^{j}x^k$$
这是一个k+2次多项式，对g(m)可以O(k)计算前k+2项，因为后面的自然数幂和与j无关，所以可以处理成前缀和（不过对最终复杂度没有意义）。
剩余的部分：
$$f(n)=\sum_{i=0}^{n}g(a+id)$$
这是一个k+3次多项式，强制插值。对于前k+3项，每一项需要O(k)的插值运算，复杂度O(k^2)。
从这道题可以看出相邻插值嵌套的计算是O(k^2)，而且适用于多层嵌套复杂度不变。
这道题有一点比较坑，模数*2会爆int，要使用unsigned int。
#include<cstdio>
#include<cstring>
#include<algorithm>
#define int unsigned int
using namespace std;
const int maxn=,MOD=;
int f[maxn],g[maxn],a,b,k,n,mx,fac[maxn],fav[maxn];
int M(int x){return x>=MOD?x-MOD:x;}
int power(int x,int k){int ans=;while(k){if(k&)ans=1ll*ans*x%MOD;x=1ll*x*x%MOD;k>>=;}return ans;}
int calc(int *g,int u){
if(u<=mx)return g[u];//?
int v=;
for(int i=;i<=mx;i++)v=1ll*v*(u-i+MOD)%MOD;
int ans=;
for(int i=;i<=mx;i++){
int t=1ll*fav[i-]*fav[mx-i]%MOD;
if((mx-i)&)t=M(MOD-t);
ans=M(ans+1ll*g[i]*v%MOD*power(M(u-i+MOD),MOD-)%MOD*t%MOD);
}
return ans;
}
#undef int
int main(){
#define int unsigned int
int T;scanf("%d",&T);
fac[]=;for(int i=;i<=;i++)fac[i]=1ll*fac[i-]*i%MOD;
fav[]=power(fac[],MOD-);for(int i=;i>=;i--)fav[i-]=1ll*fav[i]*i%MOD;
while(T--){
scanf("%d%d%d%d",&k,&a,&n,&b);mx=k+;
g[]=;
for(int i=;i<=mx;i++){
g[i]=g[i-];
for(int j=;j<=i;j++){
g[i]=M(g[i]+power(j,k));
}
}
f[]=calc(g,a);//!!!
for(int i=;i<=mx;i++){
f[i]=M(f[i-]+calc(g,M(a+1ll*i*b%MOD)));
}
printf("%d\n",calc(f,n));
}
return ;
}