【tyvj1858】xlkxc（拉格朗日插值）

传送门

题意：
求\(\sum_{i=0}^n\sum_{j=1}^{a+id}\sum_{k=1}^{j}k^K,n,a,d\leq 10^9,K\leq 100\)。

思路：
最右边这个和式为一个最高项次数为\(k+1\)的多项式；
中间这个和式加上右边的和式就是一个最高项次数为\(k+2\)的多项式；
然后整个式子为\(k+3\)次的多项式。
然后拉格朗日插一插就行。

/*
 * Author:  heyuhhh
 * Created Time:  2019/11/20 19:00:18
 */
#include <bits/stdc++.h>
#define MP make_pair
#define fi first
#define se second
#define sz(x) (int)(x).size()
#define all(x) (x).begin(), (x).end()
#define INF 0x3f3f3f3f
#define Local
#ifdef Local
  #define dbg(args...) do { cout << #args << " -> "; err(args); } while (0)
  void err() { std::cout << '\n'; }
  template<typename T, typename...Args>
  void err(T a, Args...args) { std::cout << a << ' '; err(args...); }
#else
  #define dbg(...)
#endif
void pt() {std::cout << '\n'; }
template<typename T, typename...Args>
void pt(T a, Args...args) {std::cout << a << ' '; pt(args...); }
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
//head
const int N = 150, MOD = 1234567891;

int k, a, n, d;

ll qpow(ll a, ll b) {
    ll ans = 1;
    while(b) {
        if(b & 1) ans = ans * a % MOD;   
        a = a * a % MOD;
        b >>= 1;
    }
    return ans;   
}

struct Lagrange {
	static const int SIZE = N;
	ll f[SIZE], fac[SIZE], inv[SIZE], pre[SIZE], suf[SIZE];
	int n;
	inline void add(ll &x, int y) {
		x += y;
		if(x >= MOD) x -= MOD;
	}
	void init(int _n) {
		n = _n;
		fac[0] = 1;
		for (int i = 1; i < SIZE; ++i) fac[i] = fac[i - 1] * i % MOD;
	    inv[SIZE - 1] = qpow(fac[SIZE - 1], MOD - 2);
		for (int i = SIZE - 1; i >= 1; --i) inv[i - 1] = inv[i] * i % MOD;
        f[0] = 0;
	}
	ll calc(ll x) {
		if (x <= n) return f[x];
		pre[0] = x % MOD;
		for (int i = 1; i <= n; ++i) pre[i] = pre[i - 1] * ((x - i) % MOD) % MOD;
		suf[n] = (x - n) % MOD;
		for (int i = n - 1; i >= 0; --i) suf[i] = suf[i + 1] * ((x - i) % MOD) % MOD;
		ll res = 0;
		for (int i = 0; i <= n; ++i) {
			ll tmp = f[i] * inv[n - i] % MOD * inv[i] % MOD;
			if (i) tmp = tmp * pre[i - 1] % MOD;
			if (i < n) tmp = tmp * suf[i + 1] % MOD;
			if ((n - i) & 1) tmp = MOD - tmp;
			add(res, tmp);
		}
		return res;
	}
}A, B, C;

void run(){
    cin >> k >> a >> n >> d;
    A.init(k + 1);
    for(int i = 1; i <= k + 1; i++) A.f[i] = (A.f[i - 1] + qpow(i, k)) % MOD;
    B.init(k + 2);
    for(int i = 1; i <= k + 2; i++) B.f[i] = (B.f[i - 1] + A.calc(i)) % MOD;
    C.init(k + 3);
    for(int i = 0; i <= k + 3; i++) {
        if(i == 0) C.f[i] = B.calc(a);
        else C.f[i] = (C.f[i - 1] + B.calc(a + 1ll * i * d)) % MOD;
    }
    ll res = C.calc(n);
    cout << res << '\n';
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0); cout.tie(0);
    cout << fixed << setprecision(20);
    int T; cin >> T;
    while(T--) run();
	return 0;
}