#include using namespace std; typedef long long ll; typedef long double ld; #define inf 1000000007 #define MP make_pair #define MT make_tuple #define PB push_back #define fi first #define se second #define rep(i,n) for(int i = 0; i < (int)(n); ++i) #define rrep(i,n) for(int i = (int)n-1; i >= 0; --i) #define srep(i,a,b) for(int i = (int)a; i < (int)(b); ++i) #define all(x) (x).begin(),(x).end() #define SUM(v) accumulate(all(v), 0LL) #define MIN(v) *min_element(all(v)) #define MAX(v) *max_element(all(v)) #define lb(c, x) distance((c).begin(), lower_bound(all(c), (x))) #define ub(c, x) distance((c).begin(), upper_bound(all(c), (x))) #define UNIQUE(x) sort(all(x)), x.erase(unique(all(x)), x.end()) #define SZ(c) (int)(c).size() template ostream& operator << (ostream& os, vector& vec) { os << "{"; for (int i = 0; i<(int)vec.size(); i++) { os << vec[i] << (i + 1 == (int)vec.size() ? "" : ", "); } os << "}"; return os; } // pair出力 template ostream& operator << (ostream& os, pair pair_var) { os << "(" << pair_var.first << ", " << pair_var.second << ")"; return os; } // map出力 template ostream& operator << (ostream& os, map& map_var) { os << "{"; for(auto itr = map_var.begin(); itr != map_var.end(); itr++){ os << "(" << itr->first << ", " << itr->second << ")"; itr++; if(itr != map_var.end()) os << ", "; itr--; } os << "}"; return os; } // set 出力 template ostream& operator << (ostream& os, set& set_var) { os << "{"; for(auto itr = set_var.begin(); itr != set_var.end(); itr++){ os << (*itr); ++itr; if(itr != set_var.end()) os << ", "; itr--; } os << "}"; return os; } // tuple 出力 template void out(ostream &os,const Tuple &t){} template void out(ostream &os,const Tuple &t){ if(N)os<<", "; os<(t); out(os,t); } template ostream& operator<<(ostream &os, const tuple &t){ os<<"("; out<0,tuple,Ts...>(os,t); os<<")"; return os; } #define overload2(_1, _2, name, ...) name #define vec(type, name, ...) vector name(__VA_ARGS__) #define VEC(type, name, size) \ vector name(size); \ IN(name) #define vv(type, name, h, ...) vector> name(h, vector(__VA_ARGS__)) #define VV(type, name, h, w) \ vector> name(h, vector(w)); \ IN(name) #define vvv(type, name, h, w, ...) vector>> name(h, vector>(w, vector(__VA_ARGS__))) #define vvvv(type, name, a, b, c, ...) \ vector>>> name(a, vector>>(b, vector>(c, vector(__VA_ARGS__)))) #define INT(...) \ int __VA_ARGS__; \ IN(__VA_ARGS__) #define LL(...) \ ll __VA_ARGS__; \ IN(__VA_ARGS__) #define STR(...) \ string __VA_ARGS__; \ IN(__VA_ARGS__) #define CHR(...) \ char __VA_ARGS__; \ IN(__VA_ARGS__) #define DBL(...) \ double __VA_ARGS__; \ IN(__VA_ARGS__) int scan() { return getchar(); } void scan(int &a) { cin >> a; } void scan(long long &a) { cin >> a; } void scan(char &a) { cin >> a; } void scan(double &a) { cin >> a; } void scan(string &a) { cin >> a; } template void scan(pair &p) { scan(p.first), scan(p.second); } template void scan(vector &); template void scan(vector &a) { for(auto &i : a) scan(i); } template void scan(T &a) { cin >> a; } void IN() {} template void IN(Head &head, Tail &...tail) { scan(head); IN(tail...); } const string YESNO[2] = {"NO", "YES"}; const string YesNo[2] = {"No", "Yes"}; const string yesno[2] = {"no", "yes"}; void YES(bool t = 1) { cout << YESNO[t] << endl; } void NO(bool t = 1) { YES(!t); } void Yes(bool t = 1) { cout << YesNo[t] << endl; } void No(bool t = 1) { Yes(!t); } void yes(bool t = 1) { cout << yesno[t] << endl; } void no(bool t = 1) { yes(!t); } #ifdef LOCAL void debug_out() { cerr << endl; } template void debug_out(Head H, Tail... T) { cerr << " " << H; debug_out(T...); } #define dbg(...) \ cerr << __LINE__ << " [" << #__VA_ARGS__ << "]:", debug_out(__VA_ARGS__) #define dump(x) cerr << __LINE__ << " " << #x << " = " << (x) << endl #else #define dbg(...) (void(0)) #define dump(x) (void(0)) #endif template std::enable_if_t::value> fill(A& array, const T& val) { array = val; } template std::enable_if_t::value> fill(A& array, const T& val) { for (auto& a : array) { fill(a, val); } } template T ceil(T x, S y) { assert(y); return (y < 0 ? ceil(-x, -y) : (x > 0 ? (x + y - 1) / y : x / y)); } template T floor(T x, S y) { assert(y); return (y < 0 ? floor(-x, -y) : (x > 0 ? x / y : x / y - (x % y == 0 ? 0 : 1))); } vector iota(int n) {vector a(n);iota(all(a), 0);return a;} template T POW(T x, int n) {T res = 1;for(; n; n >>= 1, x *= x){if(n & 1) res *= x;}return res;} ll pow2(int i) { return 1LL << i; } int topbit(signed t) { return t == 0 ? -1 : 31 - __builtin_clz(t); } int topbit(ll t) { return t == 0 ? -1 : 63 - __builtin_clzll(t); } int lowbit(signed a) { return a == 0 ? 32 : __builtin_ctz(a); } int lowbit(ll a) { return a == 0 ? 64 : __builtin_ctzll(a); } // int allbit(int n) { return (1 << n) - 1; } ll allbit(ll n) { return (1LL << n) - 1; } int popcount(signed t) { return __builtin_popcount(t); } int popcount(ll t) { return __builtin_popcountll(t); } bool ispow2(int i) { return i && (i & -i) == i; } template void fold_in(vector &v) {} template void fold_in(vector &v, Head &&a, Tail &&...tail) { for(auto e : a) v.emplace_back(e); fold_in(v, tail...); } template void renumber(vector &v) {} template void renumber(vector &v, Head &&a, Tail &&...tail) { for(auto &&e : a) e = lb(v, e); renumber(v, tail...); } template void zip(vector &head, Args &&...args) { vector v; fold_in(v, head, args...); sort(all(v)), v.erase(unique(all(v)), v.end()); renumber(v, head, args...); } template inline bool chmax(T &a, T b){ if(a inline bool chmin(T &a, T b){ if(a>b){ a = b; return true; } return false; } #include template struct LazyMontgomeryModInt { using mint = LazyMontgomeryModInt; using i32 = int32_t; using u32 = uint32_t; using u64 = uint64_t; static constexpr u32 get_r() { u32 ret = mod; for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret; return ret; } static constexpr u32 r = get_r(); static constexpr u32 n2 = -u64(mod) % mod; static_assert(r * mod == 1, "invalid, r * mod != 1"); static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30"); static_assert((mod & 1) == 1, "invalid, mod % 2 == 0"); u32 a; constexpr LazyMontgomeryModInt() : a(0) {} constexpr LazyMontgomeryModInt(const int64_t &b) : a(reduce(u64(b % mod + mod) * n2)){}; static constexpr u32 reduce(const u64 &b) { return (b + u64(u32(b) * u32(-r)) * mod) >> 32; } constexpr mint &operator+=(const mint &b) { if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod; return *this; } constexpr mint &operator-=(const mint &b) { if (i32(a -= b.a) < 0) a += 2 * mod; return *this; } constexpr mint &operator*=(const mint &b) { a = reduce(u64(a) * b.a); return *this; } constexpr mint &operator/=(const mint &b) { *this *= b.inverse(); return *this; } constexpr mint operator+(const mint &b) const { return mint(*this) += b; } constexpr mint operator-(const mint &b) const { return mint(*this) -= b; } constexpr mint operator*(const mint &b) const { return mint(*this) *= b; } constexpr mint operator/(const mint &b) const { return mint(*this) /= b; } constexpr bool operator==(const mint &b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); } constexpr bool operator!=(const mint &b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); } constexpr mint operator-() const { return mint() - mint(*this); } constexpr mint pow(u64 n) const { mint ret(1), mul(*this); while (n > 0) { if (n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } constexpr mint inverse() const { return pow(mod - 2); } friend ostream &operator<<(ostream &os, const mint &b) { return os << b.get(); } friend istream &operator>>(istream &is, mint &b) { int64_t t; is >> t; b = LazyMontgomeryModInt(t); return (is); } constexpr u32 get() const { u32 ret = reduce(a); return ret >= mod ? ret - mod : ret; } static constexpr u32 get_mod() { return mod; } }; #line 2 "ntt/ntt-avx2.hpp" #line 2 "modint/simd-montgomery.hpp" #line 4 "modint/simd-montgomery.hpp" __attribute__((target("sse4.2"))) inline __m128i my128_mullo_epu32( const __m128i &a, const __m128i &b) { return _mm_mullo_epi32(a, b); } __attribute__((target("sse4.2"))) inline __m128i my128_mulhi_epu32( const __m128i &a, const __m128i &b) { __m128i a13 = _mm_shuffle_epi32(a, 0xF5); __m128i b13 = _mm_shuffle_epi32(b, 0xF5); __m128i prod02 = _mm_mul_epu32(a, b); __m128i prod13 = _mm_mul_epu32(a13, b13); __m128i prod = _mm_unpackhi_epi64(_mm_unpacklo_epi32(prod02, prod13), _mm_unpackhi_epi32(prod02, prod13)); return prod; } __attribute__((target("sse4.2"))) inline __m128i montgomery_mul_128( const __m128i &a, const __m128i &b, const __m128i &r, const __m128i &m1) { return _mm_sub_epi32( _mm_add_epi32(my128_mulhi_epu32(a, b), m1), my128_mulhi_epu32(my128_mullo_epu32(my128_mullo_epu32(a, b), r), m1)); } __attribute__((target("sse4.2"))) inline __m128i montgomery_add_128( const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) { __m128i ret = _mm_sub_epi32(_mm_add_epi32(a, b), m2); return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("sse4.2"))) inline __m128i montgomery_sub_128( const __m128i &a, const __m128i &b, const __m128i &m2, const __m128i &m0) { __m128i ret = _mm_sub_epi32(a, b); return _mm_add_epi32(_mm_and_si128(_mm_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("avx2"))) inline __m256i my256_mullo_epu32( const __m256i &a, const __m256i &b) { return _mm256_mullo_epi32(a, b); } __attribute__((target("avx2"))) inline __m256i my256_mulhi_epu32( const __m256i &a, const __m256i &b) { __m256i a13 = _mm256_shuffle_epi32(a, 0xF5); __m256i b13 = _mm256_shuffle_epi32(b, 0xF5); __m256i prod02 = _mm256_mul_epu32(a, b); __m256i prod13 = _mm256_mul_epu32(a13, b13); __m256i prod = _mm256_unpackhi_epi64(_mm256_unpacklo_epi32(prod02, prod13), _mm256_unpackhi_epi32(prod02, prod13)); return prod; } __attribute__((target("avx2"))) inline __m256i montgomery_mul_256( const __m256i &a, const __m256i &b, const __m256i &r, const __m256i &m1) { return _mm256_sub_epi32( _mm256_add_epi32(my256_mulhi_epu32(a, b), m1), my256_mulhi_epu32(my256_mullo_epu32(my256_mullo_epu32(a, b), r), m1)); } __attribute__((target("avx2"))) inline __m256i montgomery_add_256( const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) { __m256i ret = _mm256_sub_epi32(_mm256_add_epi32(a, b), m2); return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2), ret); } __attribute__((target("avx2"))) inline __m256i montgomery_sub_256( const __m256i &a, const __m256i &b, const __m256i &m2, const __m256i &m0) { __m256i ret = _mm256_sub_epi32(a, b); return _mm256_add_epi32(_mm256_and_si256(_mm256_cmpgt_epi32(m0, ret), m2), ret); } #line 4 "ntt/ntt-avx2.hpp" namespace ntt_inner { using u64 = uint64_t; constexpr uint32_t get_pr(uint32_t mod) { if (mod == 2) return 1; u64 ds[32] = {}; int idx = 0; u64 m = mod - 1; for (u64 i = 2; i * i <= m; ++i) { if (m % i == 0) { ds[idx++] = i; while (m % i == 0) m /= i; } } if (m != 1) ds[idx++] = m; uint32_t pr = 2; while (1) { int flg = 1; for (int i = 0; i < idx; ++i) { u64 a = pr, b = (mod - 1) / ds[i], r = 1; while (b) { if (b & 1) r = r * a % mod; a = a * a % mod; b >>= 1; } if (r == 1) { flg = 0; break; } } if (flg == 1) break; ++pr; } return pr; } constexpr int SZ_FFT_BUF = 1 << 23; uint32_t _buf1[SZ_FFT_BUF] __attribute__((aligned(64))); uint32_t _buf2[SZ_FFT_BUF] __attribute__((aligned(64))); } // namespace ntt_inner template struct NTT { static constexpr uint32_t mod = mint::get_mod(); static constexpr uint32_t pr = ntt_inner::get_pr(mint::get_mod()); static constexpr int level = __builtin_ctzll(mod - 1); mint dw[level], dy[level]; mint *buf1, *buf2; constexpr NTT() { setwy(level); union raw_cast { mint dat; uint32_t _; }; buf1 = &(((raw_cast *)(ntt_inner::_buf1))->dat); buf2 = &(((raw_cast *)(ntt_inner::_buf2))->dat); } constexpr void setwy(int k) { mint w[level], y[level]; w[k - 1] = mint(pr).pow((mod - 1) / (1 << k)); y[k - 1] = w[k - 1].inverse(); for (int i = k - 2; i > 0; --i) w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1]; dw[0] = dy[0] = w[1] * w[1]; dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2]; for (int i = 3; i < k; ++i) { dw[i] = dw[i - 1] * y[i - 2] * w[i]; dy[i] = dy[i - 1] * w[i - 2] * y[i]; } } __attribute__((target("avx2"))) void ntt(mint *a, int n) { int k = n ? __builtin_ctz(n) : 0; if (k == 0) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; return; } if (k & 1) { int v = 1 << (k - 1); if (v < 8) { for (int j = 0; j < v; ++j) { mint ajv = a[j + v]; a[j + v] = a[j] - ajv; a[j] += ajv; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m2 = _mm256_set1_epi32(mod + mod); int j0 = 0; int j1 = v; for (; j0 < v; j0 += 8, j1 += 8) { __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); __m256i naj = montgomery_add_256(T0, T1, m2, m0); __m256i najv = montgomery_sub_256(T0, T1, m2, m0); _mm256_storeu_si256((__m256i *)(a + j0), naj); _mm256_storeu_si256((__m256i *)(a + j1), najv); } } } int u = 1 << (2 + (k & 1)); int v = 1 << (k - 2 - (k & 1)); mint one = mint(1); mint imag = dw[1]; while (v) { if (v == 1) { mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { ww = xx * xx, wx = ww * xx; mint t0 = a[jh + 0], t1 = a[jh + 1] * xx; mint t2 = a[jh + 2] * ww, t3 = a[jh + 3] * wx; mint t0p2 = t0 + t2, t1p3 = t1 + t3; mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag; a[jh + 0] = t0p2 + t1p3, a[jh + 1] = t0p2 - t1p3; a[jh + 2] = t0m2 + t1m3, a[jh + 3] = t0m2 - t1m3; xx *= dw[__builtin_ctz((jh += 4))]; } } else if (v == 4) { const __m128i m0 = _mm_set1_epi32(0); const __m128i m1 = _mm_set1_epi32(mod); const __m128i m2 = _mm_set1_epi32(mod + mod); const __m128i r = _mm_set1_epi32(mint::r); const __m128i Imag = _mm_set1_epi32(imag.a); mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = j1 + v; int j3 = j2 + v; int je = v; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P2 = montgomery_add_128(T0, T2, m2, m0); const __m128i T1P3 = montgomery_add_128(T1, T3, m2, m0); const __m128i T0M2 = montgomery_sub_128(T0, T2, m2, m0); const __m128i T1M3 = montgomery_mul_128( montgomery_sub_128(T1, T3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_sub_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_add_128(T0M2, T1M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M2, T1M3, m2, m0)); } } else { ww = xx * xx, wx = ww * xx; const __m128i WW = _mm_set1_epi32(ww.a); const __m128i WX = _mm_set1_epi32(wx.a); const __m128i XX = _mm_set1_epi32(xx.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i MT1 = montgomery_mul_128(T1, XX, r, m1); const __m128i MT2 = montgomery_mul_128(T2, WW, r, m1); const __m128i MT3 = montgomery_mul_128(T3, WX, r, m1); const __m128i T0P2 = montgomery_add_128(T0, MT2, m2, m0); const __m128i T1P3 = montgomery_add_128(MT1, MT3, m2, m0); const __m128i T0M2 = montgomery_sub_128(T0, MT2, m2, m0); const __m128i T1M3 = montgomery_mul_128( montgomery_sub_128(MT1, MT3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_sub_128(T0P2, T1P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_add_128(T0M2, T1M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M2, T1M3, m2, m0)); } } xx *= dw[__builtin_ctz((jh += 4))]; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i m2 = _mm256_set1_epi32(mod + mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i Imag = _mm256_set1_epi32(imag.a); mint ww = one, xx = one, wx = one; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = j1 + v; int j3 = j2 + v; int je = v; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P2 = montgomery_add_256(T0, T2, m2, m0); const __m256i T1P3 = montgomery_add_256(T1, T3, m2, m0); const __m256i T0M2 = montgomery_sub_256(T0, T2, m2, m0); const __m256i T1M3 = montgomery_mul_256( montgomery_sub_256(T1, T3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_sub_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_add_256(T0M2, T1M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M2, T1M3, m2, m0)); } } else { ww = xx * xx, wx = ww * xx; const __m256i WW = _mm256_set1_epi32(ww.a); const __m256i WX = _mm256_set1_epi32(wx.a); const __m256i XX = _mm256_set1_epi32(xx.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i MT1 = montgomery_mul_256(T1, XX, r, m1); const __m256i MT2 = montgomery_mul_256(T2, WW, r, m1); const __m256i MT3 = montgomery_mul_256(T3, WX, r, m1); const __m256i T0P2 = montgomery_add_256(T0, MT2, m2, m0); const __m256i T1P3 = montgomery_add_256(MT1, MT3, m2, m0); const __m256i T0M2 = montgomery_sub_256(T0, MT2, m2, m0); const __m256i T1M3 = montgomery_mul_256( montgomery_sub_256(MT1, MT3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_sub_256(T0P2, T1P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_add_256(T0M2, T1M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M2, T1M3, m2, m0)); } } xx *= dw[__builtin_ctz((jh += 4))]; } } u <<= 2; v >>= 2; } } __attribute__((target("avx2"))) void intt(mint *a, int n, int normalize = true) { int k = n ? __builtin_ctz(n) : 0; if (k == 0) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; if (normalize) { a[0] *= mint(2).inverse(); a[1] *= mint(2).inverse(); } return; } int u = 1 << (k - 2); int v = 1; mint one = mint(1); mint imag = dy[1]; while (u) { if (v == 1) { mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { ww = xx * xx, yy = xx * imag; mint t0 = a[jh + 0], t1 = a[jh + 1]; mint t2 = a[jh + 2], t3 = a[jh + 3]; mint t0p1 = t0 + t1, t2p3 = t2 + t3; mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy; a[jh + 0] = t0p1 + t2p3, a[jh + 2] = (t0p1 - t2p3) * ww; a[jh + 1] = t0m1 + t2m3, a[jh + 3] = (t0m1 - t2m3) * ww; xx *= dy[__builtin_ctz(jh += 4)]; } } else if (v == 4) { const __m128i m0 = _mm_set1_epi32(0); const __m128i m1 = _mm_set1_epi32(mod); const __m128i m2 = _mm_set1_epi32(mod + mod); const __m128i r = _mm_set1_epi32(mint::r); const __m128i Imag = _mm_set1_epi32(imag.a); mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = v + v; int j3 = j2 + v; for (; j0 < v; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0); const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0); const __m128i T0M1 = montgomery_sub_128(T0, T1, m2, m0); const __m128i T2M3 = montgomery_mul_128( montgomery_sub_128(T2, T3, m2, m0), Imag, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j2), montgomery_sub_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_add_128(T0M1, T2M3, m2, m0)); _mm_storeu_si128((__m128i *)(a + j3), montgomery_sub_128(T0M1, T2M3, m2, m0)); } } else { ww = xx * xx, yy = xx * imag; const __m128i WW = _mm_set1_epi32(ww.a); const __m128i XX = _mm_set1_epi32(xx.a); const __m128i YY = _mm_set1_epi32(yy.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 4, j1 += 4, j2 += 4, j3 += 4) { const __m128i T0 = _mm_loadu_si128((__m128i *)(a + j0)); const __m128i T1 = _mm_loadu_si128((__m128i *)(a + j1)); const __m128i T2 = _mm_loadu_si128((__m128i *)(a + j2)); const __m128i T3 = _mm_loadu_si128((__m128i *)(a + j3)); const __m128i T0P1 = montgomery_add_128(T0, T1, m2, m0); const __m128i T2P3 = montgomery_add_128(T2, T3, m2, m0); const __m128i T0M1 = montgomery_mul_128( montgomery_sub_128(T0, T1, m2, m0), XX, r, m1); __m128i T2M3 = montgomery_mul_128( montgomery_sub_128(T2, T3, m2, m0), YY, r, m1); _mm_storeu_si128((__m128i *)(a + j0), montgomery_add_128(T0P1, T2P3, m2, m0)); _mm_storeu_si128( (__m128i *)(a + j2), montgomery_mul_128(montgomery_sub_128(T0P1, T2P3, m2, m0), WW, r, m1)); _mm_storeu_si128((__m128i *)(a + j1), montgomery_add_128(T0M1, T2M3, m2, m0)); _mm_storeu_si128( (__m128i *)(a + j3), montgomery_mul_128(montgomery_sub_128(T0M1, T2M3, m2, m0), WW, r, m1)); } } xx *= dy[__builtin_ctz(jh += 4)]; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i m2 = _mm256_set1_epi32(mod + mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i Imag = _mm256_set1_epi32(imag.a); mint ww = one, xx = one, yy = one; u <<= 2; for (int jh = 0; jh < u;) { if (jh == 0) { int j0 = 0; int j1 = v; int j2 = v + v; int j3 = j2 + v; for (; j0 < v; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0); const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0); const __m256i T0M1 = montgomery_sub_256(T0, T1, m2, m0); const __m256i T2M3 = montgomery_mul_256( montgomery_sub_256(T2, T3, m2, m0), Imag, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j2), montgomery_sub_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_add_256(T0M1, T2M3, m2, m0)); _mm256_storeu_si256((__m256i *)(a + j3), montgomery_sub_256(T0M1, T2M3, m2, m0)); } } else { ww = xx * xx, yy = xx * imag; const __m256i WW = _mm256_set1_epi32(ww.a); const __m256i XX = _mm256_set1_epi32(xx.a); const __m256i YY = _mm256_set1_epi32(yy.a); int j0 = jh * v; int j1 = j0 + v; int j2 = j1 + v; int j3 = j2 + v; int je = j1; for (; j0 < je; j0 += 8, j1 += 8, j2 += 8, j3 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); const __m256i T2 = _mm256_loadu_si256((__m256i *)(a + j2)); const __m256i T3 = _mm256_loadu_si256((__m256i *)(a + j3)); const __m256i T0P1 = montgomery_add_256(T0, T1, m2, m0); const __m256i T2P3 = montgomery_add_256(T2, T3, m2, m0); const __m256i T0M1 = montgomery_mul_256( montgomery_sub_256(T0, T1, m2, m0), XX, r, m1); const __m256i T2M3 = montgomery_mul_256( montgomery_sub_256(T2, T3, m2, m0), YY, r, m1); _mm256_storeu_si256((__m256i *)(a + j0), montgomery_add_256(T0P1, T2P3, m2, m0)); _mm256_storeu_si256( (__m256i *)(a + j2), montgomery_mul_256(montgomery_sub_256(T0P1, T2P3, m2, m0), WW, r, m1)); _mm256_storeu_si256((__m256i *)(a + j1), montgomery_add_256(T0M1, T2M3, m2, m0)); _mm256_storeu_si256( (__m256i *)(a + j3), montgomery_mul_256(montgomery_sub_256(T0M1, T2M3, m2, m0), WW, r, m1)); } } xx *= dy[__builtin_ctz(jh += 4)]; } } u >>= 4; v <<= 2; } if (k & 1) { v = 1 << (k - 1); if (v < 8) { for (int j = 0; j < v; ++j) { mint ajv = a[j] - a[j + v]; a[j] += a[j + v]; a[j + v] = ajv; } } else { const __m256i m0 = _mm256_set1_epi32(0); const __m256i m2 = _mm256_set1_epi32(mod + mod); int j0 = 0; int j1 = v; for (; j0 < v; j0 += 8, j1 += 8) { const __m256i T0 = _mm256_loadu_si256((__m256i *)(a + j0)); const __m256i T1 = _mm256_loadu_si256((__m256i *)(a + j1)); __m256i naj = montgomery_add_256(T0, T1, m2, m0); __m256i najv = montgomery_sub_256(T0, T1, m2, m0); _mm256_storeu_si256((__m256i *)(a + j0), naj); _mm256_storeu_si256((__m256i *)(a + j1), najv); } } } if (normalize) { mint invn = mint(n).inverse(); for (int i = 0; i < n; i++) a[i] *= invn; } } __attribute__((target("avx2"))) void inplace_multiply( int l1, int l2, int zero_padding = true) { int l = l1 + l2 - 1; int M = 4; while (M < l) M <<= 1; if (zero_padding) { for (int i = l1; i < M; i++) ntt_inner::_buf1[i] = 0; for (int i = l2; i < M; i++) ntt_inner::_buf2[i] = 0; } const __m256i m0 = _mm256_set1_epi32(0); const __m256i m1 = _mm256_set1_epi32(mod); const __m256i r = _mm256_set1_epi32(mint::r); const __m256i N2 = _mm256_set1_epi32(mint::n2); for (int i = 0; i < l1; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i)); __m256i b = montgomery_mul_256(a, N2, r, m1); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), b); } for (int i = 0; i < l2; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf2 + i)); __m256i b = montgomery_mul_256(a, N2, r, m1); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf2 + i), b); } ntt(buf1, M); ntt(buf2, M); for (int i = 0; i < M; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i)); __m256i b = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf2 + i)); __m256i c = montgomery_mul_256(a, b, r, m1); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), c); } intt(buf1, M, false); const __m256i INVM = _mm256_set1_epi32((mint(M).inverse()).a); for (int i = 0; i < l; i += 8) { __m256i a = _mm256_loadu_si256((__m256i *)(ntt_inner::_buf1 + i)); __m256i b = montgomery_mul_256(a, INVM, r, m1); __m256i c = my256_mulhi_epu32(my256_mullo_epu32(b, r), m1); __m256i d = _mm256_and_si256(_mm256_cmpgt_epi32(c, m0), m1); __m256i e = _mm256_sub_epi32(d, c); _mm256_storeu_si256((__m256i *)(ntt_inner::_buf1 + i), e); } } void ntt(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; ntt(buf1, M); for (int i = 0; i < M; i++) a[i].a = buf1[i].a; } void intt(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; intt(buf1, M, true); for (int i = 0; i < M; i++) a[i].a = buf1[i].a; } vector multiply(const vector &a, const vector &b) { if (a.size() == 0 && b.size() == 0) return vector{}; int l = a.size() + b.size() - 1; if (min(a.size(), b.size()) <= 40) { vector s(l); for (int i = 0; i < (int)a.size(); ++i) for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j]; return s; } assert(l <= ntt_inner::SZ_FFT_BUF); int M = 4; while (M < l) M <<= 1; for (int i = 0; i < (int)a.size(); ++i) buf1[i].a = a[i].a; for (int i = (int)a.size(); i < M; ++i) buf1[i].a = 0; for (int i = 0; i < (int)b.size(); ++i) buf2[i].a = b[i].a; for (int i = (int)b.size(); i < M; ++i) buf2[i].a = 0; ntt(buf1, M); ntt(buf2, M); for (int i = 0; i < M; ++i) buf1[i].a = mint::reduce(uint64_t(buf1[i].a) * buf2[i].a); intt(buf1, M, false); vector s(l); mint invm = mint(M).inverse(); for (int i = 0; i < l; ++i) s[i] = buf1[i] * invm; return s; } void ntt_doubling(vector &a) { int M = (int)a.size(); for (int i = 0; i < M; i++) buf1[i].a = a[i].a; intt(buf1, M); mint r = 1, zeta = mint(pr).pow((mint::get_mod() - 1) / (M << 1)); for (int i = 0; i < M; i++) buf1[i] *= r, r *= zeta; ntt(buf1, M); a.resize(2 * M); for (int i = 0; i < M; i++) a[M + i].a = buf1[i].a; } }; #line 5 "ntt/arbitrary-ntt.hpp" namespace ArbitraryNTT { using i64 = int64_t; using u128 = __uint128_t; constexpr int32_t m0 = 167772161; constexpr int32_t m1 = 469762049; constexpr int32_t m2 = 754974721; using mint0 = LazyMontgomeryModInt; using mint1 = LazyMontgomeryModInt; using mint2 = LazyMontgomeryModInt; constexpr int r01 = mint1(m0).inverse().get(); constexpr int r02 = mint2(m0).inverse().get(); constexpr int r12 = mint2(m1).inverse().get(); constexpr int r02r12 = i64(r02) * r12 % m2; constexpr i64 w1 = m0; constexpr i64 w2 = i64(m0) * m1; template vector mul(const vector &a, const vector &b) { static NTT ntt; vector s(a.size()), t(b.size()); for (int i = 0; i < (int)a.size(); ++i) s[i] = i64(a[i] % submint::get_mod()); for (int i = 0; i < (int)b.size(); ++i) t[i] = i64(b[i] % submint::get_mod()); return ntt.multiply(s, t); } template vector multiply(const vector &s, const vector &t, int mod) { auto d0 = mul(s, t); auto d1 = mul(s, t); auto d2 = mul(s, t); int n = d0.size(); vector ret(n); const int W1 = w1 % mod; const int W2 = w2 % mod; for (int i = 0; i < n; i++) { int n1 = d1[i].get(), n2 = d2[i].get(), a = d0[i].get(); int b = i64(n1 + m1 - a) * r01 % m1; int c = (i64(n2 + m2 - a) * r02r12 + i64(m2 - b) * r12) % m2; ret[i] = (i64(a) + i64(b) * W1 + i64(c) * W2) % mod; } return ret; } template vector multiply(const vector &a, const vector &b) { if (a.size() == 0 && b.size() == 0) return {}; if (min(a.size(), b.size()) < 128) { vector ret(a.size() + b.size() - 1); for (int i = 0; i < (int)a.size(); ++i) for (int j = 0; j < (int)b.size(); ++j) ret[i + j] += a[i] * b[j]; return ret; } vector s(a.size()), t(b.size()); for (int i = 0; i < (int)a.size(); ++i) s[i] = a[i].get(); for (int i = 0; i < (int)b.size(); ++i) t[i] = b[i].get(); vector u = multiply(s, t, mint::get_mod()); vector ret(u.size()); for (int i = 0; i < (int)u.size(); ++i) ret[i] = mint(u[i]); return ret; } template vector multiply_u128(const vector &s, const vector &t) { if (s.size() == 0 && t.size() == 0) return {}; if (min(s.size(), t.size()) < 128) { vector ret(s.size() + t.size() - 1); for (int i = 0; i < (int)s.size(); ++i) for (int j = 0; j < (int)t.size(); ++j) ret[i + j] += i64(s[i]) * t[j]; return ret; } auto d0 = mul(s, t); auto d1 = mul(s, t); auto d2 = mul(s, t); int n = d0.size(); vector ret(n); for (int i = 0; i < n; i++) { i64 n1 = d1[i].get(), n2 = d2[i].get(); i64 a = d0[i].get(); u128 b = (n1 + m1 - a) * r01 % m1; u128 c = ((n2 + m2 - a) * r02r12 + (m2 - b) * r12) % m2; ret[i] = a + b * w1 + c * w2; } return ret; } } // namespace ArbitraryNTT #line 2 "fps/formal-power-series.hpp" template struct FormalPowerSeries : vector { using vector::vector; using FPS = FormalPowerSeries; FPS &operator+=(const FPS &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; } FPS &operator+=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] += r; return *this; } FPS &operator-=(const FPS &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; } FPS &operator-=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] -= r; return *this; } FPS &operator*=(const mint &v) { for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v; return *this; } FPS &operator/=(const FPS &r) { if (this->size() < r.size()) { this->clear(); return *this; } int n = this->size() - r.size() + 1; if ((int)r.size() <= 64) { FPS f(*this), g(r); g.shrink(); mint coeff = g.back().inverse(); for (auto &x : g) x *= coeff; int deg = (int)f.size() - (int)g.size() + 1; int gs = g.size(); FPS quo(deg); for (int i = deg - 1; i >= 0; i--) { quo[i] = f[i + gs - 1]; for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j]; } *this = quo * coeff; this->resize(n, mint(0)); return *this; } return *this = ((*this).rev().pre(n) * r.rev().inv(n)).pre(n).rev(); } FPS &operator%=(const FPS &r) { *this -= *this / r * r; shrink(); return *this; } FPS operator+(const FPS &r) const { return FPS(*this) += r; } FPS operator+(const mint &v) const { return FPS(*this) += v; } FPS operator-(const FPS &r) const { return FPS(*this) -= r; } FPS operator-(const mint &v) const { return FPS(*this) -= v; } FPS operator*(const FPS &r) const { return FPS(*this) *= r; } FPS operator*(const mint &v) const { return FPS(*this) *= v; } FPS operator/(const FPS &r) const { return FPS(*this) /= r; } FPS operator%(const FPS &r) const { return FPS(*this) %= r; } FPS operator-() const { FPS ret(this->size()); for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i]; return ret; } void shrink() { while (this->size() && this->back() == mint(0)) this->pop_back(); } FPS rev() const { FPS ret(*this); reverse(begin(ret), end(ret)); return ret; } FPS dot(FPS r) const { FPS ret(min(this->size(), r.size())); for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i]; return ret; } FPS pre(int sz) const { return FPS(begin(*this), begin(*this) + min((int)this->size(), sz)); } FPS operator>>(int sz) const { if ((int)this->size() <= sz) return {}; FPS ret(*this); ret.erase(ret.begin(), ret.begin() + sz); return ret; } FPS operator<<(int sz) const { FPS ret(*this); ret.insert(ret.begin(), sz, mint(0)); return ret; } FPS diff() const { const int n = (int)this->size(); FPS ret(max(0, n - 1)); mint one(1), coeff(1); for (int i = 1; i < n; i++) { ret[i - 1] = (*this)[i] * coeff; coeff += one; } return ret; } FPS integral() const { const int n = (int)this->size(); FPS ret(n + 1); ret[0] = mint(0); if (n > 0) ret[1] = mint(1); auto mod = mint::get_mod(); for (int i = 2; i <= n; i++) ret[i] = (-ret[mod % i]) * (mod / i); for (int i = 0; i < n; i++) ret[i + 1] *= (*this)[i]; return ret; } mint eval(mint x) const { mint r = 0, w = 1; for (auto &v : *this) r += w * v, w *= x; return r; } FPS log(int deg = -1) const { assert((*this)[0] == mint(1)); if (deg == -1) deg = (int)this->size(); return (this->diff() * this->inv(deg)).pre(deg - 1).integral(); } FPS pow(int64_t k, int deg = -1) const { const int n = (int)this->size(); if (deg == -1) deg = n; for (int i = 0; i < n; i++) { if ((*this)[i] != mint(0)) { if (i * k > deg) return FPS(deg, mint(0)); mint rev = mint(1) / (*this)[i]; FPS ret = (((*this * rev) >> i).log(deg) * k).exp(deg) * ((*this)[i].pow(k)); ret = (ret << (i * k)).pre(deg); if ((int)ret.size() < deg) ret.resize(deg, mint(0)); return ret; } } return FPS(deg, mint(0)); } static void *ntt_ptr; static void set_fft(); FPS &operator*=(const FPS &r); void ntt(); void intt(); void ntt_doubling(); static int ntt_pr(); FPS inv(int deg = -1) const; FPS exp(int deg = -1) const; }; template void *FormalPowerSeries::ntt_ptr = nullptr; /** * @brief 多項式/形式的冪級数ライブラリ * @docs docs/fps/formal-power-series.md */ #line 5 "fps/arbitrary-fps.hpp" template void FormalPowerSeries::set_fft() { ntt_ptr = nullptr; } template void FormalPowerSeries::ntt() { exit(1); } template void FormalPowerSeries::intt() { exit(1); } template void FormalPowerSeries::ntt_doubling() { exit(1); } template int FormalPowerSeries::ntt_pr() { exit(1); } template FormalPowerSeries& FormalPowerSeries::operator*=( const FormalPowerSeries& r) { if (this->empty() || r.empty()) { this->clear(); return *this; } auto ret = ArbitraryNTT::multiply(*this, r); return *this = FormalPowerSeries(ret.begin(), ret.end()); } template FormalPowerSeries FormalPowerSeries::inv(int deg) const { assert((*this)[0] != mint(0)); if (deg == -1) deg = (*this).size(); FormalPowerSeries ret({mint(1) / (*this)[0]}); for (int i = 1; i < deg; i <<= 1) ret = (ret + ret - ret * ret * (*this).pre(i << 1)).pre(i << 1); return ret.pre(deg); } template FormalPowerSeries FormalPowerSeries::exp(int deg) const { assert((*this).size() == 0 || (*this)[0] == mint(0)); if (deg == -1) deg = (int)this->size(); FormalPowerSeries ret({mint(1)}); for (int i = 1; i < deg; i <<= 1) { ret = (ret * (pre(i << 1) + mint(1) - ret.log(i << 1))).pre(i << 1); } return ret.pre(deg); } #line 7 "verify/verify-yuki/yuki-0215-nth-term.test.cpp" using mint = LazyMontgomeryModInt<998244353>; using vm = vector; using vvm = vector; using fps = FormalPowerSeries; #line 2 "fps/nth-term.hpp" #line 2 "fps/berlekamp-massey.hpp" template vector BerlekampMassey(const vector &s) { const int N = (int)s.size(); vector b, c; b.reserve(N + 1); c.reserve(N + 1); b.push_back(mint(1)); c.push_back(mint(1)); mint y = mint(1); for (int ed = 1; ed <= N; ed++) { int l = int(c.size()), m = int(b.size()); mint x = 0; for (int i = 0; i < l; i++) x += c[i] * s[ed - l + i]; b.emplace_back(mint(0)); m++; if (x == mint(0)) continue; mint freq = x / y; if (l < m) { auto tmp = c; c.insert(begin(c), m - l, mint(0)); for (int i = 0; i < m; i++) c[m - 1 - i] -= freq * b[m - 1 - i]; b = tmp; y = x; } else { for (int i = 0; i < m; i++) c[l - 1 - i] -= freq * b[m - 1 - i]; } } reverse(begin(c), end(c)); return c; } #line 2 "fps/kitamasa.hpp" #line 4 "fps/kitamasa.hpp" template mint LinearRecurrence(long long k, FormalPowerSeries Q, FormalPowerSeries P) { Q.shrink(); mint ret = 0; if (P.size() >= Q.size()) { auto R = P / Q; P -= R * Q; P.shrink(); if (k < (int)R.size()) ret += R[k]; } if ((int)P.size() == 0) return ret; FormalPowerSeries::set_fft(); if (FormalPowerSeries::ntt_ptr == nullptr) { P.resize((int)Q.size() - 1); while (k) { auto Q2 = Q; for (int i = 1; i < (int)Q2.size(); i += 2) Q2[i] = -Q2[i]; auto S = P * Q2; auto T = Q * Q2; if (k & 1) { for (int i = 1; i < (int)S.size(); i += 2) P[i >> 1] = S[i]; for (int i = 0; i < (int)T.size(); i += 2) Q[i >> 1] = T[i]; } else { for (int i = 0; i < (int)S.size(); i += 2) P[i >> 1] = S[i]; for (int i = 0; i < (int)T.size(); i += 2) Q[i >> 1] = T[i]; } k >>= 1; } return ret + P[0]; } else { int N = 1; while (N < (int)Q.size()) N <<= 1; P.resize(2 * N); Q.resize(2 * N); P.ntt(); Q.ntt(); vector S(2 * N), T(2 * N); vector btr(N); for (int i = 0, logn = __builtin_ctz(N); i < (1 << logn); i++) { btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (logn - 1)); } mint dw = mint(FormalPowerSeries::ntt_pr()) .inverse() .pow((mint::get_mod() - 1) / (2 * N)); while (k) { mint inv2 = mint(2).inverse(); // even degree of Q(x)Q(-x) T.resize(N); for (int i = 0; i < N; i++) T[i] = Q[(i << 1) | 0] * Q[(i << 1) | 1]; S.resize(N); if (k & 1) { // odd degree of P(x)Q(-x) for (auto &i : btr) { S[i] = (P[(i << 1) | 0] * Q[(i << 1) | 1] - P[(i << 1) | 1] * Q[(i << 1) | 0]) * inv2; inv2 *= dw; } } else { // even degree of P(x)Q(-x) for (int i = 0; i < N; i++) { S[i] = (P[(i << 1) | 0] * Q[(i << 1) | 1] + P[(i << 1) | 1] * Q[(i << 1) | 0]) * inv2; } } swap(P, S); swap(Q, T); k >>= 1; if (k < N) break; P.ntt_doubling(); Q.ntt_doubling(); } P.intt(); Q.intt(); return ret + (P * (Q.inv()))[k]; } } template mint kitamasa(long long N, FormalPowerSeries Q, FormalPowerSeries a) { assert(!Q.empty() && Q[0] != 0); if (N < (int)a.size()) return a[N]; assert((int)a.size() >= int(Q.size()) - 1); auto P = a.pre((int)Q.size() - 1) * Q; P.resize(Q.size() - 1); return LinearRecurrence(N, Q, P); } /** * @brief 線形漸化式の高速計算 * @docs docs/fps/kitamasa.md */ #line 5 "fps/nth-term.hpp" template mint nth_term(long long n, const vector &s) { using fps = FormalPowerSeries; auto bm = BerlekampMassey(s); return kitamasa(n, fps{begin(bm), end(bm)}, fps{begin(s), end(s)}); } const int M = 10200; mint dp[M][101]; mint ng[M][101]; mint s[M]; int main(){ INT(n,m); for(int i=2;i<=m;i++){ dp[0][i] = 1; ng[0][i] = 1; } for(int c=1;c=i)ng[c][i] += ng[c-i][i]; s[c] += dp[c][i] - ng[c][i]; } for(int i=2;i<=m;i++){ mint P = s[c] - dp[c][i] + ng[c][i]; ng[c][i] += P; dp[c][i] += P; } } fps f(M); rep(i,M){ f[i] = s[i]; } f[0] = 1; dbg(f); cout << nth_term(n,f) << endl; return 0; }