Description：

给定 $k,a,n,d,p$ ， $f(i)=1^k+2^k+3^k+......+i^k$ ， $g(x)=f(1)+f(2)+f(3)+....+f(x)$ ，求 $(g(a)+g(a+d)+g(a+2d)+......+g(a+nd))\mod p$

$1\leqslant k\leqslant 123,1\leqslant a,n,d\leqslant 123456789,p=1234567891$

Solution：

首先 $f(x)$ 就是一个自然数幂和，是一个 $k+1$ 次多项式，那么 $g(x)$ 就是一个 $k+2$ 次多项式，可以把 $a+nd$ 看成一个一次函数，那嵌套一下 $g(a+wd)$ 还是一个 $k+2$ 次多项式，那么 $\begin{align}\sum_{w=0}^ng(a+wd)\end{align}$ 就是一个 $k+3$ 次多项式，用拉格朗日插值法选 $k+4$ 个点把 $n$ 的值插出来即可。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

typedef unsigned long long ll;

ll k,a,n,d;

#define MAXK 140

ll pows[MAXK];

#define MOD 1234567891

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % MOD;

		a = a * a % MOD;

		b = b >> 1;

	}

	return res;

}

ll p[MAXK],q[MAXK];

ll power_sum(ll n)

{

	if(n < MAXK)return pows[n];

	ll ans = 0;

	p[0] = n;for(ll i = 1;i <= k + 1;++i)p[i] = p[i - 1] * (n - i) % MOD;

	q[0] = 1;for(ll i = 1;i <= k + 1;++i)q[i] = q[i - 1] * i % MOD;

	for(ll i = 0;i <= k + 1;++i)

	{

		ll res = pows[i];

		res = res * p[k + 1] % MOD;

		res = res * power(n - i,MOD - 2) % MOD;

		res = res * power(q[i],MOD - 2) % MOD;

		res = res * power(q[k + 1 - i],MOD - 2) % MOD;

		ans = ((ans + res * (((k + i) & 1) ? 1 : -1)) % MOD + MOD) % MOD;

	}

	return ans;

}

ll fs[MAXK];

ll g[MAXK];

ll calc_g(ll n)

{

	n = n % MOD;

	if(n < MAXK)return g[n];

	ll ans = 0;

	n = n % MOD;

	p[0] = n % MOD;

	for(ll i = 1;i <= k + 2;++i)p[i] = p[i - 1] * ((n - i + MOD) % MOD) % MOD;

	q[0] = 1;for(ll i = 1;i <= k + 2;++i)q[i] = q[i - 1] * i % MOD;

	for(ll i = 0;i <= k + 2;++i)

	{

		ll res = fs[i];

		res = res * p[k + 2] % MOD;

		res = res * power((n - i + MOD) % MOD,MOD - 2) % MOD;

		res = res * power(q[i],MOD - 2) % MOD;

		res = res * power(q[k + 2 - i],MOD - 2) % MOD;

		res = (res * (((k + i) & 1) ? -1 : 1) + MOD) % MOD;

		ans = (ans + res) % MOD;

	}

	return ans;

}

ll sum[MAXK];

ll calc_sum(ll n)

{

	if(n < MAXK)return sum[n];

	ll ans = 0;

	for(ll i = 0;i <= k + 3;++i)

	{

		ll res = sum[i];//cout << i << " " << sum[i] << endl;

		for(ll j = 0;j <= k + 3;++j)

		{

			if(i == j)continue;

			res = res * (n - j) % MOD;

			res = res * power((i - j + MOD) % MOD,MOD - 2) % MOD;

		}

		ans = (ans + res) % MOD;

	}

	return ans;

}

int main()

{

	int testcases;

	scanf("%d",&testcases);

	for(int c = 1;c <= testcases;++c)

	{

		scanf("%lld%lld%lld%lld",&k,&a,&n,&d);

		for(ll i = 1;i < MAXK;++i)pows[i] = (pows[i - 1] + power(i,k)) % MOD;

		for(ll i = 1;i < MAXK;++i)g[i] = (g[i - 1] + pows[i]) % MOD;

		for(ll i = 1;i < MAXK;++i)fs[i] = (power_sum(i) + fs[i - 1]) % MOD;

		sum[0] = calc_g(a);

		for(ll i = 1;i < MAXK;++i)sum[i] = (sum[i - 1] + calc_g(a + i * d)) % MOD;

		//for(ll i = 0;i < MAXK;++i)cout << sum[i] << " ";cout << endl;

		cout << calc_sum(n) << endl;

	}

	return 0;

}