【校内训练/BZOJ3453】XLkxc（拉格朗日插值）

题目大意：求
$$\sum _{i=0}^n\sum _{j=1}^{a+id}\sum _{l=1}^j{l}^k$$
模数为 $1234567891$，$1\le k\le 123，0\le a,n,d\le 123456789$
模拟赛的数据 $1\le k\le 3000$

• 先考虑 $k\le 123$ 的数据。
• 我们设 $f(n)={\sum }_{i=1}^n{i}^k，g(n)={\sum }_{i=1}^{n}f(i)，h(n)={\sum }_{i=0}^{n}g(a+id)$。
• 根据拉格朗日插值公式（确定 $n$ 个点，求 $n-1$次多项式）：

$$f(x)=\sum _{i=1}^n{y}_i\prod _{j̸=i}\frac{x-x_j}{x_i-x_j}$$

• 对 $g,h$ 插值即可。
• 显然 $f$ 是 $k+1$ 次多项式，$g$ 是 $k+2$ 次多项式，$h$ 是 $k+3$ 次多项式。
• 然后我们暴力求出 $g$ 的前 $k+4$ 个点值，然后求 $h$ 的前几个点值时，直接代进去对 $g$ 插值即可。最后对 $h$ 插值求 $h(n)$。
• 然后我们考虑一下 $k\le 3000$，这样需要我们实现线性插值。
• 对于这种我们点值的 $x$ 都是 $1$ 到 $n$ 的特殊的插值，我们可以线性插值。
• 具体做法就是我们 $O(1)$ 算 ${\prod }_{j̸=i}\frac{x-x_j}{x_i-x_j}$。
• 分子我们插值之前可以预处理然后 $O(1)$ 计算，分母就是两个阶乘，讨论一下正负就做完了。
• 顺带说一下，模数不是质数，但是和 $1$ 到 $3000$ 的数互质，所以用 $exgcd$ 求逆元。

#include <map>
#include <set>
#include <cmath>
#include <cctype>
#include <vector>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>

typedef long long s64; 

const int MaxN = 4e3 + 5; 
const s64 mod = 1234567891; 

s64 T, K, a, n, d; 
s64 ans = 0; 

s64 pw[MaxN]; 
s64 fac[MaxN], ind[MaxN]; 
s64 f[MaxN], g[MaxN], h[MaxN]; 

inline void add(s64 &a, const s64 &b)
{
	a += b; 
	if (a >= mod) a -= mod; 
}

inline s64 ex_gcd(s64 a, s64 b, s64 &x, s64 &y)
{
	if (!b)
		return x = 1, y = 0, a; 
	s64 res = ex_gcd(b, a % b, y, x); 
	return y -= a / b * x, res; 
}

inline s64 inv(s64 a)
{
	s64 x, y; 
	return ex_gcd(a, mod, x, y), (x % mod + mod) % mod; 
}

inline s64 qpow(s64 a, s64 b)
{
	s64 res = 1; 
	for (; b; b >>= 1, a = 1LL * a * a % mod)
		if (b & 1) res = 1LL * res * a % mod; 
	return res; 
}

inline s64 Lagrange(int n, s64 *y, s64 K)
{
	static s64 pre_mul[MaxN], suf_mul[MaxN]; 
	
	s64 res = 0; 
	pre_mul[0] = suf_mul[n + 1] = 1; 
	
	for (int i = 1; i <= n; ++i)
		pre_mul[i] = (1LL * pre_mul[i - 1] * (K - i) % mod + mod) % mod; 
	for (int i = n; i >= 1; --i)
		suf_mul[i] = (1LL * suf_mul[i + 1] * (K - i) % mod + mod) % mod; 
	
	for (int i = 1; i <= n; ++i)
	{
		s64 p = 1LL * pre_mul[i - 1] * suf_mul[i + 1] % mod; 
		s64 q = 1LL * ind[i - 1] * ((n - i & 1) ? mod - ind[n - i] : ind[n - i]) % mod; 
		
		res = (res + 1LL * y[i] * p % mod * q) % mod; 
	}
	
	return res; 
}

int main()
{
	freopen("study.in", "r", stdin); 
	freopen("study.out", "w", stdout); 
	
	fac[0] = 1; 
	for (int i = 1; i <= 4000; ++i)
		fac[i] = 1LL * fac[i - 1] * i % mod; 
	ind[4000] = inv(fac[4000]); 
	for (int i = 3999; i >= 0; --i)
		ind[i] = 1LL * ind[i + 1] * (i + 1) % mod; 
	
	std::cin >> T; 
	while (T--)
	{
		std::cin >> K >> a >> n >> d; 
		for (int i = 1; i <= K + 4; ++i) add(f[i] = f[i - 1], qpow(i, K)); 
		
		for (int i = 1; i <= K + 4; ++i) add(g[i] = g[i - 1], f[i]); 
		for (int i = 0; i <= K + 4; ++i)
			add(h[i] = i ? h[i - 1] : 0, Lagrange(K + 3, g, (a + i * d) % mod));
		std::cout << Lagrange(K + 4, h, n % mod) << std::endl; 
	}
	
	fclose(stdin); 
	fclose(stdout); 
	return 0; 
}