【ybt金牌导航8-3-5】多重前缀和（拉格朗日插值）

多重前缀和

题目链接：ybt金牌导航8-3-5

题目大意

求这个式子的值。
sum i=0~n sum j=1~a+id sum l=1~j l^k
其中 n,a,d,k 都给出。

思路

重新写式子：
\(\sum\limits_{i=0}^n\sum\limits_{j=1}^{a+i*d}\sum\limits_{l=1}^jl^k\)

其实看到 \(l^k\) 和 \(k\) 的小范围我们不难你想到用拉格朗日插值来搞这个 \(k\) 项的多项式。
然后你发现它前面还有两个前缀和。

那其实 \(\sum\limits_{l=1}^jl^k\) 这个是 \(k+1\) 项的，或者你通过证明 \(l^k\) 是 \(k\) 项的多项式也可以看出，每次前缀和，它的项数就 \(+1\)，那这道题的就是 \(k+3\) 项。

那你就一层一层的求下去即可。
表示出 \(l^k\)，然后拉格朗日插值可以求出 \(\sum\limits_{l=1}^jl^k\) 的，然后表示出之后又可以拉格朗日插值求 \(\sum\limits_{j=1}^{a+i*d}\sum\limits_{l=1}^j l^k\) 的，然后再表示出来再求就可以得到 \(\sum\limits_{i=0}^n\sum\limits_{j=1}^{a+i*d}\sum\limits_{l=1}^jl^k\) 了。

代码

#include<cstdio>
#define ll long long
#define mo 1234567891

using namespace std;

int T;
ll k, a, n, d, y[131], f[131], g[131];
ll pre[131], suf[131], jc[131];

ll ksm(ll x, ll y) {
	ll re = 1;
	while (y) {
		if (y & 1) re = re * x % mo;
		x = x * x % mo;
		y >>= 1;
	}
	return re;
}

ll query(ll *f, ll k, ll n) {//拉格朗日插值
	ll re = 0;
	pre[0] = 1;
	for (int i = 1; i <= k; i++) pre[i] = pre[i - 1] * ((n - i + mo) % mo) % mo;//这里 n-i 要单独取模了再运算，因为你的 n 在求 y 数组的时候是 long long 级别的
	suf[k + 1] = 1;
	for (int i = k; i >= 1; i--) suf[i] = suf[i + 1] * ((n - i + mo) % mo) % mo;
	
	for (int i = 1; i <= k; i++) {
		re = (re + f[i] * pre[i - 1] % mo * suf[i + 1] % mo * ksm(jc[i - 1] * jc[k - i] % mo * (((k - i) & 1) ? mo - 1 : 1) % mo, mo - 2) % mo) % mo;
	}
	return re;
}

int main() {
	jc[0] = 1;
	for (int i = 1; i <= 130; i++) jc[i] = jc[i - 1] * i % mo;
	
	scanf("%d", &T);
	while (T--) {
		scanf("%lld %lld %lld %lld", &k, &a, &n, &d);
		
		for (int i = 1; i <= k + 2; i++)
			f[i] = (f[i - 1] + ksm(i, k)) % mo;
		for (int i = 1; i <= k + 3; i++)
			g[i] = (g[i - 1] + query(f, k + 2, i)) % mo;
		y[0] = query(g, k + 3, a);//记得这个前缀是从 0 开始
		for (int i = 1; i <= k + 4; i++)
			y[i] = (y[i - 1] + query(g, k + 3, a + d * i)) % mo;
		
		printf("%lld\n", query(y, k + 4, n));
	}
	
	return 0;
}