結果

問題 No.3182 recurrence relation’s intersection sum
ユーザー kwm_t
提出日時 2025-06-14 14:44:14
言語 C++23
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 5 ms / 2,000 ms
コード長 4,882 bytes
コンパイル時間 5,685 ms
コンパイル使用メモリ 334,268 KB
実行使用メモリ 7,844 KB
最終ジャッジ日時 2025-06-14 14:44:22
合計ジャッジ時間 7,147 ms
ジャッジサーバーID
(参考情報)
judge3 / judge5
このコードへのチャレンジ
(要ログイン)
ファイルパターン 結果
sample AC * 3
other AC * 40
権限があれば一括ダウンロードができます

ソースコード

diff #

#include <bits/stdc++.h>
#include <atcoder/all>
using namespace std;
using namespace atcoder;
//using mint = modint1000000007;
//const int mod = 1000000007;
using mint = modint998244353;
const int mod = 998244353;
//const int INF = 1e9;
//const long long LINF = 1e18;
#define rep(i, n) for (int i = 0; i < (n); ++i)
#define rep2(i,l,r)for(int i=(l);i<(r);++i)
#define rrep(i, n) for (int i = (n) - 1; i >= 0; --i)
#define rrep2(i,l,r)for(int i=(r) - 1;i>=(l);--i)
#define all(x) (x).begin(),(x).end()
#define allR(x) (x).rbegin(),(x).rend()
#define P pair<int,int>
template<typename A, typename B> inline bool chmax(A & a, const B & b) { if (a < b) { a = b; return true; } return false; }
template<typename A, typename B> inline bool chmin(A & a, const B & b) { if (a > b) { a = b; return true; } return false; }
#include <vector>
template<class T>
struct Formal_power_series : std::vector<T> {
	using std::vector<T>::vector;
	using std::vector<T>::operator=;
	using F = Formal_power_series;
	F operator-() const {
		F res(*this);
		for (auto &e : res) e = -e;
		return res;
	}
	F inv(int d = -1) const {
		int n = this->size();
		assert(n != 0 && (*this)[0] != 0);
		if (d == -1) d = n;
		assert(d >= 0);
		F res{ (*this)[0].inv() };
		for (int m = 1; m < d; m *= 2) {
			F f(this->begin(), this->begin() + min(n, 2 * m));
			F g(res);
			f.resize(2 * m), internal::butterfly(f);
			g.resize(2 * m), internal::butterfly(g);
			for (int i = 0; i < 2 * m; ++i) f[i] *= g[i];
			internal::butterfly_inv(f);
			f.erase(f.begin(), f.begin() + m);
			f.resize(2 * m), internal::butterfly(f);
			for (int i = 0; i < 2 * m; ++i) f[i] *= g[i];
			internal::butterfly_inv(f);
			T iz = T(2 * m).inv(); iz *= -iz;
			for (int i = 0; i < m; ++i) f[i] *= iz;
			res.insert(res.end(), f.begin(), f.begin() + m);
		}
		res.resize(d);
		return res;
	}
	F &divide_inplace(const F &g, int d = -1) {
		int n = this->size();
		if (d == -1) d = n;
		assert(d >= 0);
		*this = convolution(move(*this), g.inv(d));
		this->resize(d);
		return *this;
	}
	F divide(const F &g, const int d = -1) const { return F(*this).divide_inplace(g, d); }
};

template<class T>
struct Bostan_mori {
	std::vector<T> p, q;
	Bostan_mori(std::vector<T> &_p, std::vector<T> &_q) : p(_p), q(_q) {}
	void rev(std::vector<T> &f) const {
		int d = f.size();
		for (int i = 0; i < d; ++i) if (i & 1) f[i] = -f[i];
	}
	void even(std::vector<T> &f) const {
		int d = (f.size() + 1) >> 1;
		for (int i = 0; i < d; ++i) f[i] = f[i << 1];
		f.resize(d);
	}
	void odd(std::vector<T> &f) const {
		int d = f.size() >> 1;
		for (int i = 0; i < d; ++i) f[i] = f[i << 1 | 1];
		f.resize(d);
	}
	T operator[] (long long n) const {
		std::vector<T> _p(p), _q(q), _q_rev(q);
		rev(_q_rev);
		for (; n; n >>= 1) {
			_p = convolution(move(_p), _q_rev);
			if (n & 1) odd(_p);
			else     even(_p);
			_q = convolution(move(_q), move(_q_rev));
			even(_q);
			_q_rev = _q; rev(_q_rev);
		}
		return _p[0] / _q[0];
	}
};
// https://yukicoder.me/submissions/1098946
// 線形漸化式の初めのn項の特性方程式を返す k次ならO(k^2)
// a_{i} = Σ(c_{j} * a_{i-j})で表されるような数列の先頭前k項(2*k程度で十分)を与えると
// {1,-c[1],-c[2],-c[3],-c[4]....}を返す
// ex	a_{i} = a[i-1] * 3 + a[i-2] * 1-> {1,-3,-1}
std::vector<mint> BarlekampMassey(const std::vector<mint> &A) {
	int M = A.size();
	mint y = 1; // 1つ前のxの値.
	std::vector<mint> B, C;// Cはret Bは1つ前のC.
	B.reserve(M + 1), C.reserve(M + 1);
	B.emplace_back(mint(1)), C.emplace_back(mint(1));

	for (int i = 1; i <= M; i++) {
		int lc = C.size(), lb = B.size();
		mint x = 0;
		for (int k = 0; k < lc; k++) x += C.at(k)*A.at(i - lc + k);
		B.emplace_back(mint(0)); lb++;
		if (x == 0) continue;
		mint Multi = x / y;
		if (lc < lb) {
			auto memo = C;
			if (lb > lc)C.insert(C.begin(), (lb - lc), 0);
			for (int k = 0; k < lb; k++) C.at(k) -= Multi * B.at(k);
			B = memo; y = x;
		}
		else for (int k = 0; k < lb; k++) C.at(lc - 1 - k) -= Multi * B.at(lb - 1 - k);
	}
	std::reverse(C.begin(), C.end());
	return C;
}
// 線形漸化式の初めのd項からn項目を求めるための前処理 
// d次漸化式ならO(d^2+dlogdlogN).
Bostan_mori<mint> BM_BM(const std::vector<mint> &a) {

	auto b = BarlekampMassey(a);// O(d^2)
	auto c = a;
	c.resize(b.size() - 1);
	auto d = atcoder::convolution(c, b);// O(dlogd)
	d.resize(b.size() - 1);

	Bostan_mori<mint> bm(d, b);
	return bm;// O(dlogdlogN)で取得
}
int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	//test();
	//return 0;
	long long K, L, R; cin >> K >> L >> R;
	int n = 5000;
	std::vector<mint> a(n);
	a[0] = 1;
	for (int i = 1; i < n; i++) a[i] = a[i - 1] * K + mint(i - 1).pow(K) + mint(K).pow(i - 1);
	for (int i = 1; i < n; ++i) a[i] += a[i - 1];
	auto bm = BM_BM(a);
	auto ans = bm[R] - (L ? bm[L - 1] : 0);
	cout << ans.val() << endl;
	return 0;
}
0