题目大意:求
∑ i = 0 n ∑ j = 1 a + i d ∑ l = 1 j l k \sum_{i=0}^n\sum_{j=1}^{a+id}\sum_{l=1}^jl^k i=0∑nj=1∑a+idl=1∑jlk
模数为 1234567891 1234567891 1234567891, 1 ≤ k ≤ 123 , 0 ≤ a , n , d ≤ 123456789 1\le k\le 123,0\le a,n,d\le 123456789 1≤k≤123,0≤a,n,d≤123456789
模拟赛的数据 1 ≤ k ≤ 3000 1\le k\le 3000 1≤k≤3000
- 先考虑 k ≤ 123 k\le 123 k≤123 的数据。
- 我们设 f ( n ) = ∑ i = 1 n i k , g ( n ) = ∑ i = 1 n f ( i ) , h ( n ) = ∑ i = 0 n g ( a + i d ) f(n)=\sum_{i=1}^ni^k,g(n)=\sum_{i=1}^nf(i),h(n)=\sum_{i=0}^ng(a+id) f(n)=∑i=1nik,g(n)=∑i=1nf(i),h(n)=∑i=0ng(a+id)。
- 根据拉格朗日插值公式(确定 n n n 个点,求 n − 1 n-1 n−1次多项式):
f ( x ) = ∑ i = 1 n y i ∏ j ≠ i x − x j x i − x j f(x)=\sum_{i=1}^ny_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} f(x)=i=1∑nyij̸=i∏xi−xjx−xj
- 对 g , h g,h g,h 插值即可。
- 显然 f f f 是 k + 1 k+1 k+1 次多项式, g g g 是 k + 2 k+2 k+2 次多项式, h h h 是 k + 3 k+3 k+3 次多项式。
- 然后我们暴力求出 g g g 的前 k + 4 k+4 k+4 个点值,然后求 h h h 的前几个点值时,直接代进去对 g g g 插值即可。最后对 h h h 插值求 h ( n ) h(n) h(n)。
- 然后我们考虑一下 k ≤ 3000 k\le 3000 k≤3000,这样需要我们实现线性插值。
- 对于这种我们点值的 x x x 都是 1 1 1 到 n n n 的特殊的插值,我们可以线性插值。
- 具体做法就是我们 O ( 1 ) O(1) O(1) 算 ∏ j ≠ i x − x j x i − x j \prod_{j\neq i}\frac{x-x_j}{x_i-x_j} ∏j̸=ixi−xjx−xj。
- 分子我们插值之前可以预处理然后 O ( 1 ) O(1) O(1) 计算,分母就是两个阶乘,讨论一下正负就做完了。
- 顺带说一下,模数不是质数,但是和 1 1 1 到 3000 3000 3000 的数互质,所以用 e x g c d exgcd 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;
}