Submission #19559283


Source Code Expand

Copy
#include <bits/stdc++.h>
using namespace std;
/*
#include <atcoder/all>
using namespace atcoder;
*/
/*
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
using bll = boost::multiprecision::cpp_int;
using bdouble = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<100>>;
using namespace boost::multiprecision;
*/
#if defined(LOCAL_TEST) || defined(LOCAL_DEV)
	#define BOOST_STACKTRACE_USE_ADDR2LINE
	#define BOOST_STACKTRACE_ADDR2LINE_LOCATION /usr/local/opt/binutils/bin/addr2line
	#define _GNU_SOURCE 1
	#include <boost/stacktrace.hpp>
#endif
#ifdef LOCAL_TEST
	namespace std {
		template<typename T> class dvector : public std::vector<T> {
		public:
			dvector() : std::vector<T>() {}
			explicit dvector(size_t n, const T& value = T()) : std::vector<T>(n, value) {}
			dvector(const std::vector<T>& v) : std::vector<T>(v) {}
			dvector(const std::initializer_list<T> il) : std::vector<T>(il) {}
			dvector(const std::string::iterator first, const std::string::iterator last) : std::vector<T>(first, last) {}
			dvector(const typename std::vector<T>::iterator first, const typename std::vector<T>::iterator last) : std::vector<T>(first, last) {}
			dvector(const typename std::vector<T>::reverse_iterator first, const typename std::vector<T>::reverse_iterator last) : std::vector<T>(first, last) {}
			dvector(const typename std::vector<T>::const_iterator first, const typename std::vector<T>::const_iterator last) : std::vector<T>(first, last) {}
			dvector(const typename std::vector<T>::const_reverse_iterator first, const typename std::vector<T>::const_reverse_iterator last) : std::vector<T>(first, last) {}
			T& operator[](size_t n) {
				if (this->size() <= n) { std::cerr << boost::stacktrace::stacktrace() << '\n' << "vector::_M_range_check: __n (which is " << n << ") >= this->size() (which is " << this->size() << ")" << '\n'; } return this->at(n);
			}
			const T& operator[](size_t n) const {
				if (this->size() <= n) { std::cerr << boost::stacktrace::stacktrace() << '\n' << "vector::_M_range_check: __n (which is " << n << ") >= this->size() (which is " << this->size() << ")" << '\n'; } return this->at(n);
			}
		};
	}
	class dbool {
	private:
		bool boolvalue;
	public:
		dbool() : boolvalue(false) {}
		dbool(bool b) : boolvalue(b) {}
		operator bool&() { return boolvalue; }
		operator const bool&() const { return boolvalue; }
	};
	#define vector dvector
	#define bool dbool
	class SIGFPE_exception : std::exception {};
	class SIGSEGV_exception : std::exception {};
	void catch_SIGFPE([[maybe_unused]] int e) { std::cerr << boost::stacktrace::stacktrace() << '\n'; throw SIGFPE_exception(); }
	void catch_SIGSEGV([[maybe_unused]] int e) { std::cerr << boost::stacktrace::stacktrace() << '\n'; throw SIGSEGV_exception(); }
	signed convertedmain();
	signed main() { signal(SIGFPE, catch_SIGFPE); signal(SIGSEGV, catch_SIGSEGV); return convertedmain(); }
	#define main() convertedmain()
#endif
#ifdef LOCAL_DEV
	template<typename T1, typename T2> std::ostream& operator<<(std::ostream& s, const std::pair<T1, T2>& p) {
		return s << "(" << p.first << ", " << p.second << ")"; }
	template <typename T, size_t N> std::ostream& operator<<(std::ostream& s, const std::array<T, N>& a) {
		s << "{ "; for (size_t i = 0; i < N; ++i){ s << a[i] << "\t"; } s << "}"; return s; }
	template<typename T> std::ostream& operator<<(std::ostream& s, const std::set<T>& se) {
		s << "{ "; for (auto itr = se.begin(); itr != se.end(); ++itr){ s << (*itr) << "\t"; } s << "}"; return s; }
	template<typename T> std::ostream& operator<<(std::ostream& s, const std::multiset<T>& se) {
		s << "{ "; for (auto itr = se.begin(); itr != se.end(); ++itr){ s << (*itr) << "\t"; } s << "}"; return s; }
	template<typename T1, typename T2> std::ostream& operator<<(std::ostream& s, const std::map<T1, T2>& m) {
		s << "{\n"; for (auto itr = m.begin(); itr != m.end(); ++itr){ s << "\t" << (*itr).first << " : " << (*itr).second << "\n"; } s << "}"; return s; }
	template<typename T> std::ostream& operator<<(std::ostream& s, const std::deque<T>& v) {
		for (size_t i = 0; i < v.size(); ++i){ s << v[i]; if (i < v.size() - 1) s << "\t"; } return s; }
	template<typename T> std::ostream& operator<<(std::ostream& s, const std::vector<T>& v) {
		for (size_t i = 0; i < v.size(); ++i){ s << v[i]; if (i < v.size() - 1) s << "\t"; } return s; }
	template<typename T> std::ostream& operator<<(std::ostream& s, const std::vector<std::vector<T>>& vv) {
		s << "\\\n"; for (size_t i = 0; i < vv.size(); ++i){ s << vv[i] << "\n"; } return s; }
	void debug_impl() { std::cerr << '\n'; }
	template<typename Head, typename... Tail> void debug_impl(Head head, Tail... tail) { std::cerr << " " << head << (sizeof...(tail) ? "," : ""); debug_impl(tail...); }
	#define debug(...) do { std::cerr << ":" << __LINE__ << " (" << #__VA_ARGS__ << ") ="; debug_impl(__VA_ARGS__); } while (false)
	constexpr inline long long prodlocal([[maybe_unused]] long long prod, [[maybe_unused]] long long local) { return local; }
#else
	#define debug(...) do {} while (false)
	constexpr inline long long prodlocal([[maybe_unused]] long long prod, [[maybe_unused]] long long local) { return prod; }
#endif
//#define int long long
using ll = long long;
//INT_MAX = (1<<31)-1 = 2147483647, INT64_MAX = (1LL<<63)-1 = 9223372036854775807
constexpr ll INF = numeric_limits<ll>::max() == INT_MAX ? (ll)1e9 + 7 : (ll)1e18;
constexpr ll MOD = (ll)1e9 + 7; //primitive root = 5
//constexpr ll MOD = 998244353; //primitive root = 3
constexpr double EPS = 1e-9;
constexpr ll dx[4] = {1, 0, -1, 0};
constexpr ll dy[4] = {0, 1, 0, -1};
constexpr ll dx8[8] = {1, 0, -1, 0, 1, 1, -1, -1};
constexpr ll dy8[8] = {0, 1, 0, -1, 1, -1, 1, -1};
#define rep(i, n)   for(ll i=0, i##_length=(n); i< i##_length; ++i)
#define repeq(i, n) for(ll i=1, i##_length=(n); i<=i##_length; ++i)
#define rrep(i, n)   for(ll i=(n)-1; i>=0; --i)
#define rrepeq(i, n) for(ll i=(n)  ; i>=1; --i)
#define all(v) (v).begin(), (v).end()
#define rall(v) (v).rbegin(), (v).rend()
void p() { std::cout << '\n'; }
template<typename Head, typename... Tail> void p(Head head, Tail... tail) { std::cout << head << (sizeof...(tail) ? " " : ""); p(tail...); }
template<typename T> inline void pv(std::vector<T>& v) { for(ll i=0, N=v.size(); i<N; i++) std::cout << v[i] << " \n"[i==N-1]; }
template<typename T> inline bool chmax(T& a, T b) { return a < b && (a = b, true); }
template<typename T> inline bool chmin(T& a, T b) { return a > b && (a = b, true); }
template<typename T> inline void uniq(std::vector<T>& v) { v.erase(std::unique(v.begin(), v.end()), v.end()); }
template<typename T> inline ll sz(T& v) { return v.size(); }

/*-----8<-----template-----8<-----*/

map<ll,ll> inv_cache;
struct Modint{
	unsigned long long num = 0;
	constexpr Modint() noexcept {}
	//constexpr Modint(const Modint &x) noexcept : num(x.num){}
	inline constexpr operator ll() const noexcept { return num; }
	inline constexpr Modint& operator+=(Modint x) noexcept { num += x.num; if(num >= MOD) num -= MOD; return *this; }
	inline constexpr Modint& operator++() noexcept { if(num == MOD - 1) num = 0; else num++; return *this; }
	inline constexpr Modint operator++(int) noexcept { Modint ans(*this); operator++(); return ans; }
	inline constexpr Modint operator-() const noexcept { return Modint(0) -= *this; }
	inline constexpr Modint& operator-=(Modint x) noexcept { if(num < x.num) num += MOD; num -= x.num; return *this; }
	inline constexpr Modint& operator--() noexcept { if(num == 0) num = MOD - 1; else num--; return *this; }
	inline constexpr Modint operator--(int) noexcept { Modint ans(*this); operator--(); return ans; }
	inline constexpr Modint& operator*=(Modint x) noexcept { num = (unsigned long long)(num) * x.num % MOD; return *this; }
	inline Modint& operator/=(Modint x) noexcept { return operator*=(x.inv()); }
	template<class T> constexpr Modint(T x) noexcept {
		using U = typename conditional<sizeof(T) >= 4, T, int>::type;
		U y = x; y %= U(MOD); if(y < 0) y += MOD; num = (unsigned long long)(y);
	}
	template<class T> inline constexpr Modint operator+(T x) const noexcept { return Modint(*this) += x; }
	template<class T> inline constexpr Modint& operator+=(T x) noexcept { return operator+=(Modint(x)); }
	template<class T> inline constexpr Modint operator-(T x) const noexcept { return Modint(*this) -= x; }
	template<class T> inline constexpr Modint& operator-=(T x) noexcept { return operator-=(Modint(x)); }
	template<class T> inline constexpr Modint operator*(T x) const noexcept { return Modint(*this) *= x; }
	template<class T> inline constexpr Modint& operator*=(T x) noexcept { return operator*=(Modint(x)); }
	template<class T> inline constexpr Modint operator/(T x) const noexcept { return Modint(*this) /= x; }
	template<class T> inline constexpr Modint& operator/=(T x) noexcept { return operator/=(Modint(x)); }
	inline Modint inv() const noexcept { return inv_cache.count(num) ? inv_cache[num] : inv_cache[num] = inv_calc(); }
	inline constexpr ll inv_calc() const noexcept { ll x = 0, y = 0; extgcd(num, MOD, x, y); return x; }
	static inline constexpr ll extgcd(ll a, ll b, ll &x, ll &y) noexcept { ll g = a; x = 1; y = 0; if(b){ g = extgcd(b, a % b, y, x); y -= a / b * x; } return g; }
	inline constexpr Modint pow(ll x) const noexcept { Modint ans = 1, cnt = x>=0 ? *this : inv(); if(x<0) x = -x; while(x){ if(x & 1) ans *= cnt; cnt *= cnt; x /= 2; } return ans; }
	static inline constexpr ll get_mod() { return MOD; }
};
std::istream& operator>>(std::istream& is, Modint& x){ ll a; is>>a; x = a; return is; }
inline constexpr Modint operator""_M(unsigned long long x) noexcept { return Modint(x); }
std::vector<Modint> fac(1, 1), inv(1, 1);
inline void reserve(size_t a){
	if(fac.size() >= a) return;
	if(a < fac.size() * 2) a = fac.size() * 2;
	if(a >= MOD) a = MOD;
	fac.reserve(a);
	while(fac.size() < a) fac.push_back(fac.back() * Modint(fac.size()));
	inv.resize(fac.size());
	inv.back() = fac.back().inv();
	for(ll i = inv.size() - 1; !inv[i - 1]; i--) inv[i - 1] = inv[i] * i;
}
inline Modint factorial(ll n){ if(n < 0) return 0; reserve(n + 1); return fac[n]; }
inline Modint nPk(ll n, ll r){
    if(r < 0 || n < r) return 0;
    if(n >> 24){ Modint ans = 1; for(ll i = 0; i < r; i++) ans *= n--; return ans; }
    reserve(n + 1); return fac[n] * inv[n - r];
}
inline Modint nCk(ll n, ll r){ if(r < 0 || n < r) return 0; r = min(r, n - r); reserve(r + 1); return inv[r] * nPk(n, r); }
inline Modint nHk(ll n, ll r){ return nCk(n + r - 1, n - 1); } //n種類のものから重複を許してr個選ぶ=玉r個と仕切りn-1個
inline Modint catalan(ll n){ reserve(n * 2 + 1); return fac[n * 2] * inv[n] * inv[n + 1]; }

////

//[lib]形式的冪級数FPS.cpp
//[depends on]modint.cpp
template< typename T >
struct FormalPowerSeries : vector< T > {
	using vector< T >::vector;
	using P = FormalPowerSeries;

	using MULT = function< vector< T >(P, P) >;
	using FFT = function< void(P &) >;
	using SQRT = function< T(T) >;

	static MULT &get_mult() {
		static MULT mult = nullptr;
		return mult;
	}

	static void set_mult(MULT f) {
		get_mult() = f;
	}

	static FFT &get_fft() {
		static FFT fft = nullptr;
		return fft;
	}

	static FFT &get_ifft() {
		static FFT ifft = nullptr;
		return ifft;
	}

	static void set_fft(FFT f, FFT g) {
		get_fft() = f;
		get_ifft() = g;
		if(get_mult() == nullptr) {
			auto mult = [&](P a, P b) {
				int need = a.size() + b.size() - 1;
				int nbase = 1;
				while((1 << nbase) < need) nbase++;
				int sz = 1 << nbase;
				a.resize(sz, T(0));
				b.resize(sz, T(0));
				get_fft()(a);
				get_fft()(b);
				for(int i = 0; i < sz; i++) a[i] *= b[i];
				get_ifft()(a);
				a.resize(need);
				return a;
			};
			set_mult(mult);
		}
	}

	static SQRT &get_sqrt() {
		static SQRT sqr = nullptr;
		return sqr;
	}

	static void set_sqrt(SQRT sqr) {
		get_sqrt() = sqr;
	}

	void shrink() {
		while(this->size() && this->back() == T(0)) this->pop_back();
	}

	P operator+(const P &r) const { return P(*this) += r; }

	P operator+(const T &v) const { return P(*this) += v; }

	P operator-(const P &r) const { return P(*this) -= r; }

	P operator-(const T &v) const { return P(*this) -= v; }

	P operator*(const P &r) const { return P(*this) *= r; }

	P operator*(const T &v) const { return P(*this) *= v; }

	P operator/(const P &r) const { return P(*this) /= r; }

	P operator%(const P &r) const { return P(*this) %= r; }

	P &operator+=(const P &r) {
		if(r.size() > this->size()) this->resize(r.size());
		for(int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
		return *this;
	}

	P &operator+=(const T &r) {
		if(this->empty()) this->resize(1);
		(*this)[0] += r;
		return *this;
	}

	P &operator-=(const P &r) {
		if(r.size() > this->size()) this->resize(r.size());
		for(int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
		shrink();
		return *this;
	}

	P &operator-=(const T &r) {
		if(this->empty()) this->resize(1);
		(*this)[0] -= r;
		shrink();
		return *this;
	}

	P &operator*=(const T &v) {
		const int n = (int) this->size();
		for(int k = 0; k < n; k++) (*this)[k] *= v;
		return *this;
	}

	P &operator*=(const P &r) {
		if(this->empty() || r.empty()) {
			this->clear();
			return *this;
		}
		assert(get_mult() != nullptr);
		auto ret = get_mult()(*this, r);
		return *this = P(begin(ret), end(ret));
	}

	P &operator%=(const P &r) {
		return *this -= *this / r * r;
	}

	P operator-() const {
		P ret(this->size());
		for(int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i];
		return ret;
	}

	P &operator/=(const P &r) {
		if(this->size() < r.size()) {
			this->clear();
			return *this;
		}
		int n = this->size() - r.size() + 1;
		return *this = (rev().pre(n) * r.rev().inv(n)).pre(n).rev(n);
	}

	P dot(P r) const {
		P ret(min(this->size(), r.size()));
		for(int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i];
		return ret;
	}

	P pre(int sz) const {
		return P(begin(*this), begin(*this) + min((int) this->size(), sz));
	}

	P operator>>(int sz) const {
		if(this->size() <= sz) return {};
		P ret(*this);
		ret.erase(ret.begin(), ret.begin() + sz);
		return ret;
	}

	P operator<<(int sz) const {
		P ret(*this);
		ret.insert(ret.begin(), sz, T(0));
		return ret;
	}

	P rev(int deg = -1) const {
		P ret(*this);
		if(deg != -1) ret.resize(deg, T(0));
		reverse(begin(ret), end(ret));
		return ret;
	}

	T operator()(T x) const {
		T r = 0, w = 1;
		for(auto &v : *this) {
			r += w * v;
			w *= x;
		}
		return r;
	}

	// https://opt-cp.com/fps-implementation/
	// multiply and divide (1 + cz^d)
	P mul(const ll d, const T c) {
		P ret(*this);
		int n = ret.size();
		if (c == T(1)) for(int i=n-d-1; i>=0; --i) ret[i+d] += ret[i];
		else if (c == T(-1)) for(int i=n-d-1; i>=0; --i) ret[i+d] -= ret[i];
		else for(int i=n-d-1; i>=0; --i) ret[i+d] += ret[i] * c;
		return ret;
	}
	P div(const ll d, const T c) {
		P ret(*this);
		int n = ret.size();
		if (c == T(1)) for(int i=0; i<n-d; ++i) ret[i+d] -= ret[i];
		else if (c == T(-1)) for(int i=0; i<n-d; ++i) ret[i+d] += ret[i];
		else for(int i=0; i<n-d; ++i) ret[i+d] -= ret[i] * c;
		return ret;
	}
	// sparse
	P mul(vector<pair<ll, T>> g) {
		if ((int)g.size() == 2 && g[0] == pair<ll, T>(0, 1))
			return mul(g[1].first, g[1].second);
		P ret(*this);
		int n = ret.size();
		auto [d, c] = g.front();
		if (d == 0) g.erase(g.begin());
		else c = 0;
		for(int i=n-1; i>=0; i--){
			ret[i] *= c;
			for (auto&& [j, b] : g) {
				if (j > i) break;
				ret[i] += ret[i-j] * b;
			}
		}
		return ret;
	}
	// sparse, required: "g[0] == (0, c)" and "c != 0"
	P div(vector<pair<ll, T>> g) {
		if ((int)g.size() == 2 && g[0] == pair<ll, T>(0, 1))
			return div(g[1].first, g[1].second);
		P ret(*this);
		int n = ret.size();
		auto [d, c] = g.front();
		assert(d == 0 && c != T(0));
		g.erase(g.begin());
		for(int i=0; i<n; i++) {
			for (auto&& [j, b] : g) {
				if (j > i) break;
				ret[i] -= ret[i-j] * b;
			}
			ret[i] /= c;
		}
		return ret;
	}

	P diff() const;

	P integral() const;

	// F(0) must not be 0
	P inv_fast() const;

	P inv(int deg = -1) const;

	// F(0) must be 1
	P log(int deg = -1) const;

	P sqrt(int deg = -1) const;

	// F(0) must be 0
	P exp_fast(int deg = -1) const;

	P exp(int deg = -1) const;

	P pow(int64_t k, int deg = -1) const;

	P mod_pow(int64_t k, P g) const;

	P taylor_shift(T c) const;
};
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::diff() const {
	const int n = (int) this->size();
	P ret(max(0, n - 1));
	for(int i = 1; i < n; i++) ret[i - 1] = (*this)[i] * T(i);
	return ret;
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::integral() const {
	const int n = (int) this->size();
	P ret(n + 1);
	ret[0] = T(0);
	for(int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / T(i + 1);
	return ret;
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::inv_fast() const {
	assert(((*this)[0]) != T(0));

	const int n = (int) this->size();
	P res{T(1) / (*this)[0]};

	for(int d = 1; d < n; d <<= 1) {
		P f(2 * d), g(2 * d);
		for(int j = 0; j < min(n, 2 * d); j++) f[j] = (*this)[j];
		for(int j = 0; j < d; j++) g[j] = res[j];
		get_fft()(f);
		get_fft()(g);
		for(int j = 0; j < 2 * d; j++) f[j] *= g[j];
		get_ifft()(f);
		for(int j = 0; j < d; j++) {
			f[j] = 0;
			f[j + d] = -f[j + d];
		}
		get_fft()(f);
		for(int j = 0; j < 2 * d; j++) f[j] *= g[j];
		get_ifft()(f);
		for(int j = 0; j < d; j++) f[j] = res[j];
		res = f;
	}
	return res.pre(n);
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::inv(int deg) const {
	assert(((*this)[0]) != T(0));
	const int n = (int) this->size();
	if(deg == -1) deg = n;
	if(get_fft() != nullptr) {
		P ret(*this);
		ret.resize(deg, T(0));
		return ret.inv_fast();
	}
	P ret({T(1) / (*this)[0]});
	for(int i = 1; i < deg; i <<= 1) {
		ret = (ret + ret - ret * ret * pre(i << 1)).pre(i << 1);
	}
	return ret.pre(deg);
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::log(int deg) const {
	assert((*this)[0] == 1);
	const int n = (int) this->size();
	if(deg == -1) deg = n;
	return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::sqrt(int deg) const {
	const int n = (int) this->size();
	if(deg == -1) deg = n;
	if((*this)[0] == T(0)) {
		for(int i = 1; i < n; i++) {
			if((*this)[i] != T(0)) {
				if(i & 1) return {};
				if(deg - i / 2 <= 0) break;
				auto ret = (*this >> i).sqrt(deg - i / 2);
				if(ret.empty()) return {};
				ret = ret << (i / 2);
				if(ret.size() < deg) ret.resize(deg, T(0));
				return ret;
			}
		}
		return P(deg, 0);
	}

	P ret;
	if(get_sqrt() == nullptr) {
		assert((*this)[0] == T(1));
		ret = {T(1)};
	} else {
		auto sqr = get_sqrt()((*this)[0]);
		if(sqr * sqr != (*this)[0]) return {};
		ret = {T(sqr)};
	}

	T inv2 = T(1) / T(2);
	for(int i = 1; i < deg; i <<= 1) {
		ret = (ret + pre(i << 1) * ret.inv(i << 1)) * inv2;
	}
	return ret.pre(deg);
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::exp_fast(int deg) const {
	if(deg == -1) deg = this->size();
	assert((*this)[0] == T(0));

	P inv;
	inv.reserve(deg + 1);
	inv.push_back(T(0));
	inv.push_back(T(1));

	auto inplace_integral = [&](P &F) -> void {
		const int n = (int) F.size();
		auto mod = T::get_mod();
		while((int) inv.size() <= n) {
			int i = inv.size();
			inv.push_back((-inv[mod % i]) * (mod / i));
		}
		F.insert(begin(F), T(0));
		for(int i = 1; i <= n; i++) F[i] *= inv[i];
	};

	auto inplace_diff = [](P &F) -> void {
		if(F.empty()) return;
		F.erase(begin(F));
		T coeff = 1, one = 1;
		for(int i = 0; i < (int) F.size(); i++) {
			F[i] *= coeff;
			coeff += one;
		}
	};

	P b{1, 1 < (int) this->size() ? (*this)[1] : T(0)}, c{1}, z1, z2{1, 1};
	for(int m = 2; m < deg; m *= 2) {
		auto y = b;
		y.resize(2 * m);
		get_fft()(y);
		z1 = z2;
		P z(m);
		for(int i = 0; i < m; ++i) z[i] = y[i] * z1[i];
		get_ifft()(z);
		fill(begin(z), begin(z) + m / 2, T(0));
		get_fft()(z);
		for(int i = 0; i < m; ++i) z[i] *= -z1[i];
		get_ifft()(z);
		c.insert(end(c), begin(z) + m / 2, end(z));
		z2 = c;
		z2.resize(2 * m);
		get_fft()(z2);
		P x(begin(*this), begin(*this) + min< int >(this->size(), m));
		inplace_diff(x);
		x.push_back(T(0));
		get_fft()(x);
		for(int i = 0; i < m; ++i) x[i] *= y[i];
		get_ifft()(x);
		x -= b.diff();
		x.resize(2 * m);
		for(int i = 0; i < m - 1; ++i) x[m + i] = x[i], x[i] = T(0);
		get_fft()(x);
		for(int i = 0; i < 2 * m; ++i) x[i] *= z2[i];
		get_ifft()(x);
		x.pop_back();
		inplace_integral(x);
		for(int i = m; i < min< int >(this->size(), 2 * m); ++i) x[i] += (*this)[i];
		fill(begin(x), begin(x) + m, T(0));
		get_fft()(x);
		for(int i = 0; i < 2 * m; ++i) x[i] *= y[i];
		get_ifft()(x);
		b.insert(end(b), begin(x) + m, end(x));
	}
	return P(begin(b), begin(b) + deg);
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::exp(int deg) const {
	assert((*this)[0] == T(0));
	const int n = (int) this->size();
	if(deg == -1) deg = n;
	if(get_fft() != nullptr) {
		P ret(*this);
		ret.resize(deg, T(0));
		return ret.exp_fast(deg);
	}
	P ret({T(1)});
	for(int i = 1; i < deg; i <<= 1) {
		ret = (ret * (pre(i << 1) + T(1) - ret.log(i << 1))).pre(i << 1);
	}
	return ret.pre(deg);
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::pow(int64_t k, int deg) const {
	const int n = (int) this->size();
	if(deg == -1) deg = n;
	for(int i = 0; i < n; i++) {
		if((*this)[i] != T(0)) {
			T rev = T(1) / (*this)[i];
			P ret = (((*this * rev) >> i).log() * k).exp() * ((*this)[i].pow(k));
			if(i * k > deg) return P(deg, T(0));
			ret = (ret << (i * k)).pre(deg);
			if(ret.size() < deg) ret.resize(deg, T(0));
			return ret;
		}
	}
	return *this;
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::mod_pow(int64_t k, P g) const {
	P modinv = g.rev().inv();
	auto get_div = [&](P base) {
		if(base.size() < g.size()) {
			base.clear();
			return base;
		}
		int n = base.size() - g.size() + 1;
		return (base.rev().pre(n) * modinv.pre(n)).pre(n).rev(n);
	};
	P x(*this), ret{1};
	while(k > 0) {
		if(k & 1) {
			ret *= x;
			ret -= get_div(ret) * g;
		}
		x *= x;
		x -= get_div(x) * g;
		k >>= 1;
	}
	return ret;
}
template< typename T >
typename FormalPowerSeries< T >::P FormalPowerSeries< T >::taylor_shift(T c) const {
	int n = (int) this->size();
	vector< T > fact(n), rfact(n);
	fact[0] = rfact[0] = T(1);
	for(int i = 1; i < n; i++) fact[i] = fact[i - 1] * T(i);
	rfact[n - 1] = T(1) / fact[n - 1];
	for(int i = n - 1; i > 1; i--) rfact[i - 1] = rfact[i] * T(i);
	P p(*this);
	for(int i = 0; i < n; i++) p[i] *= fact[i];
	p = p.rev();
	P bs(n, T(1));
	for(int i = 1; i < n; i++) bs[i] = bs[i - 1] * c * rfact[i] * fact[i - 1];
	p = (p * bs).pre(n);
	p = p.rev();
	for(int i = 0; i < n; i++) p[i] *= rfact[i];
	return p;
}

////

template< typename Mint >
struct NumberTheoreticTransformFriendlyModInt {

	vector< Mint > dw, idw;
	int max_base;
	Mint root;

	NumberTheoreticTransformFriendlyModInt() {
		const unsigned mod = Mint::get_mod();
		assert(mod >= 3 && mod % 2 == 1);
		auto tmp = mod - 1;
		max_base = 0;
		while(tmp % 2 == 0) tmp >>= 1, max_base++;
		root = 2;
		while(root.pow((mod - 1) >> 1) == 1) root += 1;
		assert(root.pow(mod - 1) == 1);
		dw.resize(max_base);
		idw.resize(max_base);
		for(int i = 0; i < max_base; i++) {
			dw[i] = -root.pow((mod - 1) >> (i + 2));
			idw[i] = Mint(1) / dw[i];
		}
	}

	void ntt(vector< Mint > &a) {
		const int n = (int) a.size();
		assert((n & (n - 1)) == 0);
		assert(__builtin_ctz(n) <= max_base);
		for(int m = n; m >>= 1;) {
			Mint w = 1;
			for(int s = 0, k = 0; s < n; s += 2 * m) {
				for(int i = s, j = s + m; i < s + m; ++i, ++j) {
					auto x = a[i], y = a[j] * w;
					a[i] = x + y, a[j] = x - y;
				}
				w *= dw[__builtin_ctz(++k)];
			}
		}
	}

	void intt(vector< Mint > &a, bool f = true) {
		const int n = (int) a.size();
		assert((n & (n - 1)) == 0);
		assert(__builtin_ctz(n) <= max_base);
		for(int m = 1; m < n; m *= 2) {
			Mint w = 1;
			for(int s = 0, k = 0; s < n; s += 2 * m) {
				for(int i = s, j = s + m; i < s + m; ++i, ++j) {
					auto x = a[i], y = a[j];
					a[i] = x + y, a[j] = (x - y) * w;
				}
				w *= idw[__builtin_ctz(++k)];
			}
		}
		if(f) {
			Mint inv_sz = Mint(1) / n;
			for(int i = 0; i < n; i++) a[i] *= inv_sz;
		}
	}

	vector< Mint > multiply(vector< Mint > a, vector< Mint > b) {
		int need = a.size() + b.size() - 1;
		int nbase = 1;
		while((1 << nbase) < need) nbase++;
		int sz = 1 << nbase;
		a.resize(sz, 0);
		b.resize(sz, 0);
		ntt(a);
		ntt(b);
		Mint inv_sz = Mint(1) / sz;
		for(int i = 0; i < sz; i++) a[i] *= b[i] * inv_sz;
		intt(a, false);
		a.resize(need);
		return a;
	}
};

//高速フーリエ変換
//計算量 O((n+m)log(n+m))
namespace FastFourierTransform {
	using real = double;

	struct C {
		real x, y;
		C() : x(0), y(0) {}
		C(real x, real y) : x(x), y(y) {}
		inline C operator+(const C &c) const { return C(x + c.x, y + c.y); }
		inline C operator-(const C &c) const { return C(x - c.x, y - c.y); }
		inline C operator*(const C &c) const { return C(x * c.x - y * c.y, x * c.y + y * c.x); }
		inline C conj() const { return C(x, -y); }
	};

	int base = 1;
	vector< C > rts = {{0, 0},{1, 0}};
	vector< int > rev = {0, 1};

	void ensure_base(int nbase) {
		if(nbase <= base) return;
		rev.resize(1 << nbase);
		rts.resize(1 << nbase);
		for(int i = 0; i < (1 << nbase); i++) {
			rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (nbase - 1));
		}
		while(base < nbase) {
			real angle = M_PI * 2.0 / (1 << (base + 1));
			for(int i = 1 << (base - 1); i < (1 << base); i++) {
				rts[i << 1] = rts[i];
				real angle_i = angle * (2 * i + 1 - (1 << base));
				rts[(i << 1) + 1] = C(cos(angle_i), sin(angle_i));
			}
			++base;
		}
	}

	void fft(vector< C > &a, int n) {
		assert((n & (n - 1)) == 0);
		int zeros = __builtin_ctz(n);
		ensure_base(zeros);
		int shift = base - zeros;
		for(int i = 0; i < n; i++) {
			if(i < (rev[i] >> shift)) {
				swap(a[i], a[rev[i] >> shift]);
			}
		}
		for(int k = 1; k < n; k <<= 1) {
			for(int i = 0; i < n; i += 2 * k) {
				for(int j = 0; j < k; j++) {
					C z = a[i + j + k] * rts[j + k];
					a[i + j + k] = a[i + j] - z;
					a[i + j] = a[i + j] + z;
				}
			}
		}
	}

	template<typename T>
	vector< ll > multiply(vector< T > &a, vector< T > &b) {
		int need = (int) a.size() + (int) b.size() - 1;
		int nbase = 1;
		while((1 << nbase) < need) nbase++;
		ensure_base(nbase);
		int sz = 1 << nbase;
		vector< C > fa(sz);
		for(int i = 0; i < sz; i++) {
			real x = (i < (int) a.size() ? a[i] : 0);
			real y = (i < (int) b.size() ? b[i] : 0);
			fa[i] = C(x, y);
		}
		fft(fa, sz);
		C r(0, -0.25 / (sz >> 1)), s(0, 1), t(0.5, 0);
		for(int i = 0; i <= (sz >> 1); i++) {
			int j = (sz - i) & (sz - 1);
			C z = (fa[j] * fa[j] - (fa[i] * fa[i]).conj()) * r;
			fa[j] = (fa[i] * fa[i] - (fa[j] * fa[j]).conj()) * r;
			fa[i] = z;
		}
		for(int i = 0; i < (sz >> 1); i++) {
			C A0 = (fa[i] + fa[i + (sz >> 1)]) * t;
			C A1 = (fa[i] - fa[i + (sz >> 1)]) * t * rts[(sz >> 1) + i];
			fa[i] = A0 + A1 * s;
		}
		fft(fa, sz >> 1);
		vector< ll > ret(need);
		for(int i = 0; i < need; i++) {
			ret[i] = llround(i & 1 ? fa[i >> 1].y : fa[i >> 1].x);
		}
		return ret;
	}
};
//任意mod畳み込み(Arbitrary-Mod-Convolution)
template< typename T >
struct ArbitraryModConvolution {
	using real = FastFourierTransform::real;
	using C = FastFourierTransform::C;

	ArbitraryModConvolution() = default;

	vector< T > multiply(const vector< T > &a, const vector< T > &b, int need = -1) {
		if(need == -1) need = a.size() + b.size() - 1;
		int nbase = 0;
		while((1 << nbase) < need) nbase++;
		FastFourierTransform::ensure_base(nbase);
		int sz = 1 << nbase;
		vector< C > fa(sz);
		for(int i = 0; i < (int)a.size(); i++) {
			fa[i] = C(a[i].num & ((1 << 15) - 1), a[i].num >> 15);
		}
		fft(fa, sz);
		vector< C > fb(sz);
		if(a == b) {
			fb = fa;
		} else {
			for(int i = 0; i < (int)b.size(); i++) {
				fb[i] = C(b[i].num & ((1 << 15) - 1), b[i].num >> 15);
			}
			fft(fb, sz);
		}
		real ratio = 0.25 / sz;
		C r2(0, -1), r3(ratio, 0), r4(0, -ratio), r5(0, 1);
		for(int i = 0; i <= (sz >> 1); i++) {
			int j = (sz - i) & (sz - 1);
			C a1 = (fa[i] + fa[j].conj());
			C a2 = (fa[i] - fa[j].conj()) * r2;
			C b1 = (fb[i] + fb[j].conj()) * r3;
			C b2 = (fb[i] - fb[j].conj()) * r4;
			if(i != j) {
				C c1 = (fa[j] + fa[i].conj());
				C c2 = (fa[j] - fa[i].conj()) * r2;
				C d1 = (fb[j] + fb[i].conj()) * r3;
				C d2 = (fb[j] - fb[i].conj()) * r4;
				fa[i] = c1 * d1 + c2 * d2 * r5;
				fb[i] = c1 * d2 + c2 * d1;
			}
			fa[j] = a1 * b1 + a2 * b2 * r5;
			fb[j] = a1 * b2 + a2 * b1;
		}
		fft(fa, sz);
		fft(fb, sz);
		vector< T > ret(need);
		for(int i = 0; i < need; i++) {
			ll aa = llround(fa[i].x);
			ll bb = llround(fb[i].x);
			ll cc = llround(fa[i].y);
			aa = ll(T(aa)), bb = ll(T(bb)), cc = ll(T(cc));
			ret[i] = T(aa + (bb << 15) + (cc << 30));
		}
		return ret;
	}
};

////

//[lib]高速きたまさ法.cpp
//[depends on]形式的冪級数FPS.cpp
// f に a * x^n + b を掛ける
void mul_simple(FormalPowerSeries<Modint> &f, Modint a, ll n, Modint b){
	for(ll i = f.size() - 1 ; i >= 0 ; i--){
		f[i] *= b;
		if(i >= n) f[i] += f[i - n] * a; 
	}
}
// f から a * x^n + b を割る
void div_simple(FormalPowerSeries<Modint> &f, Modint a, ll n, Modint b){
	for(ll i = 0 ; i < (ll)f.size() ; i++){
		f[i] /= b;
		if(i + n < (ll)f.size() ) f[n + i] -= f[i] * a; 
	}
}
// f/g を deg(f)次まで求める
FormalPowerSeries<Modint> div_(FormalPowerSeries<Modint> &f, FormalPowerSeries<Modint> g){
	ll n = f.size();
	return (f * g.inv(n)).pre(n);
}

// 隣接K+1項漸化式 A[n] = c[1]*A[n-1]+c[2]*A[n-2]+...+c[K]*A[n-K] のa[N]を求める
// N ... 求めたい項a[N] (0-indexed)
// c ... 漸化式の係数c
// a ... 初期解 (a[0] , a[1] , ... , a[K-1])
// x^N を fでわった剰余を求め、aと内積を取る
//計算量:O(KlogKlogN)
Modint kitamasa(ll N, vector<Modint> c, vector<Modint> a){
	ArbitraryModConvolution< Modint > fft;
	using FPS = FormalPowerSeries< Modint >;
	auto mult = [&](const FPS::P& a, const FPS::P& b) { return fft.multiply(a, b); };
	FPS::set_mult(mult);

	FormalPowerSeries<Modint> Q(c.size()+1);
	Q[0]=1;
	for(ll i=0; i<(ll)c.size(); i++) Q[i+1]=-c[i];

	FormalPowerSeries<Modint> af(a.size());
	for(ll i=0; i<(ll)a.size(); i++) af[i]=a[i];

	ll k = Q.size() - 1;
	assert( (ll)af.size() == k );
	FormalPowerSeries<Modint> P = af * Q; P.resize(k);
	while(N){
		// 初期化
		// N & 1のときはQoを、そうでないときはQxを1シフトする
		FormalPowerSeries<Modint> Pe( k - k / 2 ) , Po( k / 2 );
		FormalPowerSeries<Modint> Qe( (k + 1) - (k + 1) / 2) ;
		FormalPowerSeries<Modint> Qo( (k + 1) / 2 +  (N & 1) );
		FormalPowerSeries<Modint> Qx( (k + 1) / 2 + !(N & 1) );
		
		for(int i = 0 ; i <= k ; i++){
			if(i & 1){
				if(i != k) Po[i >> 1] = P[i];
				Qo[(i >> 1) + (N&1)] = Q[i];
				Qx[(i >> 1) +!(N&1)] = (Q[i] ? MOD - ll(Q[i]) : 0);
			}else {
				if(i != k) Pe[i >> 1] = P[i];
				Qe[i >> 1] = Q[i];
			}
		}
		if(N & 1){
			P = (Pe * Qx) + (Po * Qe) , Q = (Qe * Qe) +(Qo * Qx);
		}else{
			P = (Pe * Qe) + (Po * Qx) , Q = (Qe * Qe) + (Qo * Qx);
		}
		N >>= 1;
	}
	return P[0];
}

/*-----8<-----library-----8<-----*/

//https://atcoder.jp/contests/tdpc/tasks/tdpc_fibonacci
void solve() {
	ll K,N;
	cin>>K>>N;
	N--;

	//漸化式 A[n] = c[1]*A[n-1]+c[2]*A[n-2]+...+c[K]*A[n-K] の係数c
	vector<Modint> c(K);
	rep(i, K) c[i] = 1;

	//初期解 a[0]-a[K-1]
	vector<Modint> a(K);
	rep(i, K) a[i] = 1;

	Modint ans=kitamasa(N, c, a);
	cout << ans << endl;
}

signed main() {
	solve();
	return 0;
}

Submission Info

Submission Time
Task T - フィボナッチ
User kyon2326
Language C++ (GCC 9.2.1)
Score 8
Code Size 32828 Byte
Status AC
Exec Time 22 ms
Memory 4076 KB

Judge Result

Set Name All
Score / Max Score 8 / 8
Status
AC × 7
Set Name Test Cases
All 00, 01, 02, 03, 04, 90, 91
Case Name Status Exec Time Memory
00 AC 22 ms 4044 KB
01 AC 20 ms 4052 KB
02 AC 20 ms 4004 KB
03 AC 13 ms 3948 KB
04 AC 9 ms 4076 KB
90 AC 2 ms 3964 KB
91 AC 3 ms 3808 KB