結果
| 問題 | No.3619 Compositional Power with Schröder Coordinate |
| コンテスト | |
| ユーザー |
 edon8618
|
| 提出日時 | 2026-08-10 23:51:40 |
| 言語 | C++23 (gcc 15.2.0 + boost 1.90.0) |
| 結果 |
AC
|
| 実行時間 | 776 ms / 10,000 ms |
| + 440µs | |
| コード長 | 37,866 bytes |
| 記録 | |
| コンパイル時間 | 4,909 ms |
| コンパイル使用メモリ | 398,028 KB |
| 実行使用メモリ | 99,720 KB |
| 最終ジャッジ日時 | 2026-08-10 23:58:07 |
| 合計ジャッジ時間 | 10,849 ms |
|
ジャッジサーバーID （参考情報） |
judge1_0 / judge3_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 2 |
| other | AC * 6 |
ソースコード
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2")
#include <bits/stdc++.h>
using namespace std;
#include <atcoder/all>
using namespace atcoder;
// #include <boost/multiprecision/cpp_int.hpp>
// using namespace boost::multiprecision;
#define ll long long
// #define ld long double
#define rep(i, n) for (ll i = 0; i < (ll)(n); ++i)
#define vi vector<int>
#define vl vector<ll>
#define vd vector<double>
#define vb vector<bool>
#define vs vector<string>
#define vc vector<char>
#define ull unsigned long long
#define all(a) (a).begin(), (a).end()
#define rall(a) (a).rbegin(), (a).rend()
template<class T, class U>
inline bool chmax(T &a, const U &b) {
if (a < b) {
a = b;
return true;
}
return false;
}
template<class T, class U>
inline bool chmin(T &a, const U &b) {
if (a > b) {
a = b;
return true;
}
return false;
}
// #define ll int
// #define ll int128_t
// #define ll int256_t
// #define ll cpp_int
constexpr ll inf = (1ll << 61);
// constexpr ll inf = (1 << 30);
// const double PI=3.1415926535897932384626433832795028841971;
// uint32_t xor_x = 123456789, xor_y = 362436069, xor_z = 521288629, xor_w = 88675123;
// inline uint32_t xor_next() {
// uint32_t t = xor_x ^ (xor_x << 11);
// xor_x = xor_y; xor_y = xor_z; xor_z = xor_w;
// return xor_w = (xor_w ^ (xor_w >> 19)) ^ (t ^ (t >> 8));
// }
// inline int rnd(int max_val) { return xor_next() % max_val; }
// struct Timer {
// std::chrono::steady_clock::time_point start_time;
// Timer() {
// reset();
// }
// // 測定の起点リセット用
// void reset() {
// start_time = std::chrono::steady_clock::now();
// }
// // スタートからの経過時間をミリ秒(msec)で返す
// long long get_ms() const {
// auto now = std::chrono::steady_clock::now();
// return std::chrono::duration_cast<std::chrono::milliseconds>(now - start_time).count();
// }
// };
// ll rui(ll a,ll b){
// if(b==0)return 1;
// if(b%2==1) return a*rui(a*a,b/2);
// return rui(a*a,b/2);
// }
// vl fact;
// ll kai(ll n){
// fact.resize(n,1);
// rep(i,n-1)fact[i+1]=fact[i]*(i+1);
// }
// using mint = ld;
using mint = modint998244353;//static_modint<998244353>
// using mint = modint1000000007;//static_modint<1000000007>
// using mint = static_modint<922267487>; // 多分落とされにくい NOT ntt-friendly
// using mint = static_modint<469762049>; // ntt-friendly
// using mint = static_modint<167772161>; // ntt-friendly
// using mint = modint;//mint::set_mod(mod);
// ll const mod=1000000007ll;
// ll const mod=998244353ll;
// ll modrui(ll a,ll b,ll mod){
// a%=mod;
// if(b==0)return 1;
// if(b%2==1) return a*modrui(a*a%mod,b/2,mod)%mod;
// return modrui(a*a%mod,b/2,mod)%mod;
// }
// void incr(vl &v,ll n){// n進法
// ll k=v.size();
// v[k-1]++;
// ll now=k-1;
// while (v[now]>=n)
// {
// v[now]=0;
// if(now==0)break;
// v[now-1]++;
// now--;
// }
// return;
// }
// vector<mint> fact,invf;
// void init_modfact(ll sz){
// fact.resize(sz);
// invf.resize(sz);
// fact[0]=1;
// rep(i,sz-1){
// fact[i+1]=fact[i]*(i+1);
// }
// invf[sz-1]=1/fact[sz-1];
// for(ll i=sz-2; i>=0; i--){
// invf[i]=invf[i+1]*(i+1);
// }
// }
// mint choose(ll n,ll r){
// if(n<r || r<0 || n<0)return 0;
// return fact[n]*invf[r]*invf[n-r];
// }
// vector<mint> modpow,invpow;
// void init_modpow(ll x,ll sz){
// mint inv=1/mint(x);
// modpow.assign(sz,1);
// invpow.assign(sz,1);
// rep(i,sz-1){
// modpow[i+1]=modpow[i]*x;
// invpow[i+1]=invpow[i]*inv;
// }
// }
// long long phi(long long n) {// O(sqrt(n))
// long long res = n;
// for (long long i = 2; i * i <= n; i++) {
// if (n % i == 0) {
// res -= res / i;
// while (n % i == 0) n /= i;
// }
// }
// if (n > 1) res -= res / n;
// return res;
// }
namespace NachiaNtt {
using u32 = unsigned int;
using u64 = unsigned long long;
using mint = atcoder::modint998244353;
inline int LsbIndex(u64 x) {
#ifdef __GNUC__
return __builtin_ctzll(x);
#else
int res = 0;
while ((x & 1) == 0) { x >>= 1; res++; }
return res;
#endif
}
inline int ceil_pow2(int n) {
int x = 0;
while ((1U << x) < (u32)(n)) x++;
return x;
}
struct fft_info {
static constexpr u32 g = 3;
static constexpr int rank2 = 23;
std::array<mint, rank2 + 1> root, iroot, rate3, irate3;
fft_info() {
root[rank2] = mint(g).pow((mint::mod() - 1) >> rank2);
iroot[rank2] = root[rank2].inv();
for (int i = rank2 - 1; i >= 0; i--) {
root[i] = root[i + 1] * root[i + 1];
iroot[i] = iroot[i + 1] * iroot[i + 1];
}
mint prod = 1, iprod = 1;
for (int i = 0; i <= rank2 - 3; i++) {
rate3[i] = root[i + 3] * prod;
irate3[i] = iroot[i + 3] * iprod;
prod *= iroot[i + 3];
iprod *= root[i + 3];
}
}
};
void ButterflyLayered(mint* a, int n, int stride, int repeat) {
static const fft_info info;
int h = n * stride;
while (repeat--) {
int len = 1;
int p = h;
if (ceil_pow2(n) % 2 == 1) {
p >>= 1;
for (int i = 0; i < p; i++) {
mint l = a[i], r = a[i + p];
a[i] = l + r; a[i + p] = l - r;
}
len <<= 1;
}
for (; p > stride;) {
p >>= 2;
mint rot = 1, imag = info.root[2];
u64 mod2 = u64(mint::mod()) * mint::mod();
int offset = p;
for (int s = 0; s < len; s++) {
if (s) rot *= info.rate3[LsbIndex(~(u32)(s - 1))];
mint rot2 = rot * rot;
mint rot3 = rot2 * rot;
for (int i = offset - p; i < offset; i++) {
u64 a0 = u64(a[i].val());
u64 a1 = u64(a[i + p].val()) * rot.val();
u64 a2 = u64(a[i + 2 * p].val()) * rot2.val();
u64 a3 = u64(a[i + 3 * p].val()) * rot3.val();
u64 a1na3imag = u64(mint(a1 + mod2 - a3).val()) * imag.val();
u64 na2 = mod2 - a2;
a[i] = a0 + a2 + a1 + a3;
a[i + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3));
a[i + 2 * p] = a0 + na2 + a1na3imag;
a[i + 3 * p] = a0 + na2 + (mod2 - a1na3imag);
}
offset += p << 2;
}
len <<= 2;
}
a += h;
}
}
void IButterflyLayered(mint* a, int n, int stride, int repeat) {
static const fft_info info;
constexpr int MOD = 998244353;
while (repeat--) {
int len = n;
int p = stride;
for (; 2 < len;) {
len >>= 2;
mint irot = 1, iimag = info.iroot[2];
int offset = p;
for (int s = 0; s < len; s++) {
if (s) irot *= info.irate3[LsbIndex(~(u32)(s - 1))];
mint irot2 = irot * irot;
mint irot3 = irot2 * irot;
for (int i = offset - p; i < offset; i++) {
u64 a0 = a[i].val();
u64 a1 = a[i + p].val();
u64 a2 = a[i + 2 * p].val();
u64 a3 = a[i + 3 * p].val();
u64 a2na3iimag = mint((a2 + MOD - a3) * iimag.val()).val();
a[i] = a0 + a1 + a2 + a3;
a[i + p] = (a0 + (MOD - a1) + a2na3iimag) * irot.val();
a[i + 2 * p] = (a0 + a1 + (MOD - a2) + (MOD - a3)) * irot2.val();
a[i + 3 * p] = (a0 + (MOD - a1) + (MOD - a2na3iimag)) * irot3.val();
}
offset += p << 2;
}
p <<= 2;
}
if (len == 2) {
for (int i = 0; i < p; i++) {
mint l = a[i], r = a[i + p];
a[i] = l + r; a[i + p] = l - r;
}
p <<= 1;
}
a += p;
}
}
} // namespace NachiaNtt
/**
* dense
* inv https://judge.yosupo.jp/submission/362241
* log https://judge.yosupo.jp/submission/362244
* exp https://judge.yosupo.jp/submission/362243
* pow https://judge.yosupo.jp/submission/362252
* sqrt https://judge.yosupo.jp/submission/362246
*
* sparse
* inv https://judge.yosupo.jp/submission/362247
* log https://judge.yosupo.jp/submission/362249
* exp https://judge.yosupo.jp/submission/362248
* pow https://judge.yosupo.jp/submission/362250
* sqrt https://judge.yosupo.jp/submission/362251
*
* div,mod https://judge.yosupo.jp/submission/362266
*
*
* composition https://judge.yosupo.jp/submission/370782, based on nachia https://judge.yosupo.jp/submission/199636
* compositional_inverse https://judge.yosupo.jp/submission/370576, based on nachia https://judge.yosupo.jp/submission/199637
* Taylor Shift https://judge.yosupo.jp/submission/362264
*/
// FPS (Formal Power Series) Library
// 依存: atcoder::convolution, modint998244353
struct FPS : std::vector<mint> {
using std::vector<mint>::vector; // コンストラクタを継承
// 指定したサイズに切り詰める/伸ばす (次数制限用)
FPS pre(int sz) const {
FPS res(begin(), begin() + std::min((int)size(), sz));
res.resize(sz, 0);
return res;
}
// --- 次数削減 (末尾の0を削除して正しい次数にする) ---
void trim() {
while (!empty() && back() == 0) pop_back();
}
// --- デバッグ出力 (dump) ---
void dump(const std::string& name = "", bool as_poly = false) const {
std::cerr << (name.empty() ? "FPS" : name) << " = ";
if (empty()) {
std::cerr << "[]\n";
return;
}
if (!as_poly) {
// 配列表記: [c0, c1, c2, ...]
std::cerr << "[";
for (int i = 0; i < (int)size(); i++) {
std::cerr << (*this)[i].val() << (i == (int)size() - 1 ? "" : ", ");
}
std::cerr << "]\n";
} else {
// 多項式表記: c0 + c1*x + c2*x^2 + ...
bool first = true;
for (int i = 0; i < (int)size(); i++) {
if ((*this)[i].val() == 0) continue;
if (!first) std::cerr << " + ";
std::cerr << (*this)[i].val();
if (i == 1) std::cerr << "x";
else if (i > 1) std::cerr << "x^" << i;
first = false;
}
if (first) std::cerr << "0";
std::cerr << "\n";
}
}
// --- 四則演算 ---
FPS& operator+=(const FPS& r) {
if (r.size() > size()) resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
return *this;
}
FPS& operator-=(const FPS& r) {
if (r.size() > size()) resize(r.size());
for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
return *this;
}
FPS& operator*=(const FPS& r) {
if (empty() || r.empty()) { clear(); return *this; }
auto res = atcoder::convolution(*this, r);
return *this = FPS(res.begin(), res.end());
}
FPS& operator*=(mint c) {
for (auto& x : *this) x *= c;
return *this;
}
FPS operator/(const FPS& r) const {
FPS a = *this; a.trim();
FPS b = r; b.trim();
if (a.size() < b.size()) return FPS();
assert(!b.empty());
int k = a.size() - b.size() + 1;
std::reverse(a.begin(), a.end());
std::reverse(b.begin(), b.end());
// ここで下の方で定義している自動切り替え inv() が呼ばれる
FPS q = (a.pre(k) * b.inv(k)).pre(k);
std::reverse(q.begin(), q.end());
q.trim();
return q;
}
FPS operator%(const FPS& r) const {
FPS a = *this; a.trim();
FPS b = r; b.trim();
if (a.size() < b.size()) return a;
assert(!b.empty());
FPS q = a / b;
FPS rem = a - (q * b);
rem.resize(b.size() - 1);
rem.trim();
return rem;
}
FPS operator+(const FPS& r) const { return FPS(*this) += r; }
FPS operator-(const FPS& r) const { return FPS(*this) -= r; }
FPS operator*(const FPS& r) const { return FPS(*this) *= r; }
FPS operator*(mint c) const { return FPS(*this) *= c; }
FPS operator-() const {
FPS res = *this;
for (auto& x : res) x = -x;
return res;
}
FPS& operator/=(const FPS& r) { return *this = *this / r; }
FPS& operator%=(const FPS& r) { return *this = *this % r; }
// --- 微分・積分 ---
FPS derivative() const {
if (size() <= 1) return {0};
FPS res(size() - 1);
for (int i = 0; i < (int)res.size(); i++) res[i] = (*this)[i + 1] * (i + 1);
return res;
}
FPS integral() const {
if (empty()) return {0};
FPS res(size() + 1, 0);
for (int i = 0; i < (int)size(); i++) res[i + 1] = (*this)[i] / (i + 1);
return res;
}
// --- テイラーシフト (Taylor Shift) f(x + c) ---
// O(N \log N)
FPS shift(mint c) const {
int n = size();
if (n <= 1) return *this;
// 階乗と逆元の事前計算 (O(N))
std::vector<mint> fact(n), inv_fact(n);
fact[0] = 1;
for (int i = 1; i < n; i++) fact[i] = fact[i - 1] * i;
inv_fact[n - 1] = fact[n - 1].inv();
for (int i = n - 1; i > 0; i--) inv_fact[i - 1] = inv_fact[i] * i;
FPS A(n), C(n);
mint c_pow = 1;
for (int i = 0; i < n; i++) {
A[n - 1 - i] = (*this)[i] * fact[i]; // A を左右反転して配置
C[i] = c_pow * inv_fact[i];
c_pow *= c;
}
// 畳み込み
FPS R = A * C;
// 結果を抽出して反転を戻し、1/j! を掛ける
FPS res(n);
for (int i = 0; i < n; i++) {
res[i] = R[n - 1 - i] * inv_fact[i];
}
return res;
}
// --- 逆元 (inv) ---
FPS dense_inv(int deg = -1) const {
if (deg == -1) deg = size();
FPS res{mint(1) / (*this)[0]};
for (int d = 1; d < deg; d *= 2) {
FPS f = pre(2 * d);
FPS g = res.pre(2 * d);
FPS tmp = (f * g).pre(2 * d);
for (int i = 0; i < (int)tmp.size(); i++) tmp[i] = -tmp[i];
tmp[0] += 2;
res = (g * tmp).pre(2 * d);
}
return res.pre(deg);
}
// --- 対数関数 (log) ---
FPS dense_log(int deg = -1) const {
if (deg == -1) deg = size();
return (derivative() * dense_inv(deg)).pre(deg - 1).integral();
}
// --- 指数関数 (exp) ---
FPS dense_exp(int deg = -1) const {
if (deg == -1) deg = size();
FPS res{mint(1)};
for (int d = 1; d < deg; d *= 2) {
FPS f = pre(2 * d);
FPS g = res.pre(2 * d);
FPS diff = g.dense_log(2 * d);
for (int i = 0; i < (int)diff.size(); i++) diff[i] = -diff[i];
for (int i = 0; i < std::min((int)f.size(), (int)diff.size()); i++) diff[i] += f[i];
diff[0] += 1;
res = (g * diff).pre(2 * d);
}
return res.pre(deg);
}
// --- 累乗 (pow) ---
FPS dense_pow(long long k, int deg = -1) const {
if (deg == -1) deg = size();
if (k == 0) {
FPS res(deg, 0);
if (deg > 0) res[0] = 1;
return res;
}
int n = size();
int shift = 0;
while (shift < n && (*this)[shift] == 0) shift++;
if (shift == n || shift > (deg - 1) / k) return FPS(deg, 0);
mint c0 = (*this)[shift];
FPS b(n - shift);
for (int i = 0; i < (int)b.size(); i++) b[i] = (*this)[i + shift] / c0;
int rem = deg - shift * k;
FPS lg = b.dense_log(rem);
for (auto& x : lg) x *= mint(k);
FPS b_pow = lg.dense_exp(rem);
FPS res(deg, 0);
mint c0_pow = c0.pow(k);
for (int i = 0; i < rem; i++) res[i + shift * k] = b_pow[i] * c0_pow;
return res;
}
// --- 平方根 (sqrt) ---
// 条件を満たす g(x) が存在しない場合は、空の FPS (size = 0) を返す
FPS dense_sqrt(int deg = -1) const {
if (deg == -1) deg = size();
if (deg == 0) return FPS();
int n = size();
int shift = 0;
// 1. 最小の非ゼロ項を探す (くくり出し)
while (shift < n && (*this)[shift] == 0) shift++;
// 全て 0 の場合は 0 を返す
if (shift == n) return FPS(deg, 0);
// 最小項の次数が奇数の場合、多項式としての平方根は存在しない
if (shift % 2 != 0) return FPS();
long long c = (*this)[shift].val();
long long p = mint::mod();
long long r = -1;
// 2. Tonelli-Shanks のアルゴリズムで定数項の平方根 r を求める
auto pow_mod = [&](long long x, long long k) {
long long res = 1; x %= p;
while (k > 0) { if (k & 1) res = res * x % p; x = x * x % p; k >>= 1; }
return res;
};
if (pow_mod(c, (p - 1) / 2) == 1) { // 平方剰余であるか判定
long long q = p - 1, s = 0;
while ((q & 1) == 0) { q >>= 1; s++; }
long long z = 2;
while (pow_mod(z, (p - 1) / 2) == 1) z++;
long long cq = pow_mod(z, q);
r = pow_mod(c, (q + 1) / 2);
long long t = pow_mod(c, q);
long long m = s;
while (t != 1) {
long long t2 = t, i = 0;
for (; i < m; i++) { if (t2 == 1) break; t2 = t2 * t2 % p; }
long long b = cq;
for (int j = 0; j < m - i - 1; j++) b = b * b % p;
r = r * b % p;
cq = b * b % p;
t = t * cq % p;
m = i;
}
}
// 平方剰余が存在しない場合 (r == -1)
if (r == -1) return FPS();
// 3. ニュートン法 (Newton's Method) による FPS 平方根の復元
// x_{k+1} = (x_k + a / x_k) / 2
mint inv2 = mint(1) / 2;
FPS f_shift(begin() + shift, end());
FPS res{mint(r)};
int rem_deg = deg - shift / 2;
for (int d = 1; d < rem_deg; d *= 2) {
FPS f_cut = f_shift.pre(2 * d);
FPS g_inv = res.dense_inv(2 * d);
FPS tmp = (f_cut * g_inv).pre(2 * d);
res += tmp;
res *= inv2;
res = res.pre(2 * d);
}
// 4. シフトの復元
FPS final_res(deg, 0);
for (int i = 0; i < rem_deg && i + shift / 2 < deg; i++) {
final_res[i + shift / 2] = res[i];
}
return final_res;
}
// --- 疎な逆元 (sparse inv) ---
FPS sparse_inv(int deg = -1) const {
if (deg == -1) deg = size();
assert(!empty());
int shift = 0;
while (shift < (int)size() && (*this)[shift] == 0) shift++;
assert(shift == 0); // 定数項が0の場合は逆元が存在しない
mint inv_c0 = (*this)[0].inv();
std::vector<std::pair<int, mint>> sparse_terms;
for (int i = 1; i < std::min((int)size(), deg); i++) {
if ((*this)[i] != 0) sparse_terms.emplace_back(i, (*this)[i]);
}
FPS res(deg, 0);
res[0] = inv_c0;
for (int i = 1; i < deg; i++) {
mint sum = 0;
for (auto& p : sparse_terms) {
int j = p.first;
if (i - j < 0) continue;
sum += p.second * res[i - j];
}
res[i] = -sum * inv_c0;
}
return res;
}
// --- 疎な対数関数 (sparse log) ---
FPS sparse_log(int deg = -1) const {
if (deg == -1) deg = size();
assert(!empty() && (*this)[0] == 1);
std::vector<std::pair<int, mint>> sparse_terms;
for (int i = 1; i < std::min((int)size(), deg); i++) {
if ((*this)[i] != 0) sparse_terms.emplace_back(i, (*this)[i]);
}
std::vector<mint> inv(deg);
if (deg > 1) inv[1] = 1;
int MOD = mint::mod();
for (int i = 2; i < deg; i++) inv[i] = -inv[MOD % i] * (MOD / i);
FPS res(deg, 0);
for (int i = 1; i < deg; i++) {
mint sum = 0;
for (auto& p : sparse_terms) {
int j = p.first;
if (i - j < 0) continue;
if (i == j) sum += p.second * mint(i);
else sum -= p.second * res[i - j] * mint(i - j);
}
res[i] = sum * inv[i];
}
return res;
}
// --- 疎な指数関数 (sparse exp) ---
FPS sparse_exp(int deg = -1) const {
if (deg == -1) deg = size();
assert(empty() || (*this)[0] == 0);
std::vector<std::pair<int, mint>> sparse_terms;
for (int i = 1; i < std::min((int)size(), deg); i++) {
if ((*this)[i] != 0) sparse_terms.emplace_back(i, (*this)[i]);
}
std::vector<mint> inv(deg);
if (deg > 1) inv[1] = 1;
int MOD = mint::mod();
for (int i = 2; i < deg; i++) inv[i] = -inv[MOD % i] * (MOD / i);
FPS res(deg, 0);
if (deg > 0) res[0] = 1;
for (int i = 1; i < deg; i++) {
mint sum = 0;
for (auto& p : sparse_terms) {
int j = p.first;
if (i - j < 0) continue;
sum += p.second * res[i - j] * mint(j);
}
res[i] = sum * inv[i];
}
return res;
}
// --- 疎な累乗 (sparse pow) ---
FPS sparse_pow(long long k, int deg = -1) const {
if (deg == -1) deg = size();
if (k == 0) {
FPS res(deg, 0);
if (deg > 0) res[0] = 1;
return res;
}
int n = size();
int shift = 0;
while (shift < n && (*this)[shift] == 0) shift++;
if (shift == n) return FPS(deg, 0);
assert(!(k < 0 && shift > 0));
if (k > 0 && shift > (deg - 1) / k) return FPS(deg, 0);
mint c0 = (*this)[shift];
mint inv_c0 = mint(1) / c0;
std::vector<std::pair<int, mint>> sparse_terms;
for (int i = shift + 1; i < n; i++) {
if ((*this)[i] != 0) sparse_terms.emplace_back(i - shift, (*this)[i] * inv_c0);
}
int rem_deg = deg;
if (k > 0) rem_deg -= shift * k;
std::vector<mint> inv(rem_deg);
if (rem_deg > 1) inv[1] = 1;
int MOD = mint::mod();
for (int i = 2; i < rem_deg; i++) inv[i] = -inv[MOD % i] * (MOD / i);
FPS res(rem_deg, 0);
res[0] = 1;
mint mk = mint(k);
for (int i = 1; i < rem_deg; i++) {
mint sum = 0;
for (auto& p : sparse_terms) {
int j = p.first;
if (i - j < 0) continue;
sum += p.second * res[i - j] * (mk * mint(j) - mint(i - j));
}
res[i] = sum * inv[i];
}
FPS final_res(deg, 0);
mint c0_pow = (k < 0) ? inv_c0.pow(-k) : c0.pow(k);
int final_shift = (k > 0) ? shift * k : 0;
for (int i = 0; i < rem_deg && i + final_shift < deg; i++) {
final_res[i + final_shift] = res[i] * c0_pow;
}
return final_res;
}
// --- 疎な平方根 (sparse sqrt) ---
FPS sparse_sqrt(int deg = -1) const {
if (deg == -1) deg = size();
if (deg == 0) return FPS();
int n = size();
int shift = 0;
while (shift < n && (*this)[shift] == 0) shift++;
if (shift == n) return FPS(deg, 0);
if (shift % 2 != 0) return FPS();
long long c_val = (*this)[shift].val();
long long p = mint::mod();
long long r = -1;
// Tonelli-Shanks
auto pow_mod = [&](long long x, long long k) {
long long res = 1; x %= p;
while (k > 0) { if (k & 1) res = res * x % p; x = x * x % p; k >>= 1; }
return res;
};
if (pow_mod(c_val, (p - 1) / 2) == 1) {
long long q = p - 1, s = 0;
while ((q & 1) == 0) { q >>= 1; s++; }
long long z = 2;
while (pow_mod(z, (p - 1) / 2) == 1) z++;
long long cq = pow_mod(z, q);
r = pow_mod(c_val, (q + 1) / 2);
long long t = pow_mod(c_val, q);
long long m = s;
while (t != 1) {
long long t2 = t, i = 0;
for (; i < m; i++) { if (t2 == 1) break; t2 = t2 * t2 % p; }
long long b = cq;
for (int j = 0; j < m - i - 1; j++) b = b * b % p;
r = r * b % p;
cq = b * b % p;
t = t * cq % p;
m = i;
}
}
if (r == -1) return FPS();
mint inv_2c0 = (mint(2) * mint(c_val)).inv();
std::vector<std::pair<int, mint>> sparse_terms;
for (int i = shift + 1; i < n; i++) {
if ((*this)[i] != 0) sparse_terms.emplace_back(i - shift, (*this)[i]);
}
int rem_deg = deg - shift / 2;
std::vector<mint> inv(rem_deg);
if (rem_deg > 1) inv[1] = 1;
for (int i = 2; i < rem_deg; i++) inv[i] = -inv[p % i] * (p / i);
FPS res(rem_deg, 0);
res[0] = mint(r);
for (int i = 1; i < rem_deg; i++) {
mint sum = 0;
for (auto& [j, coef] : sparse_terms) {
if (i - j < 0) continue;
if (i == j) sum += mint(i) * coef * res[0];
else sum += coef * res[i - j] * mint(3 * j - 2 * i);
}
res[i] = sum * inv[i] * inv_2c0;
}
FPS final_res(deg, 0);
for (int i = 0; i < rem_deg && i + shift / 2 < deg; i++) {
final_res[i + shift / 2] = res[i];
}
return final_res;
}
// --- 自動切り替えフロントエンド (Auto-Switching API) ---
// 非ゼロ項の数 (K) をカウントし、自動で最適なアルゴリズムを選択する
// 切り替えの閾値 (K <= SPARSE_THRESHOLD なら sparse を使用)
static constexpr int SPARSE_THRESHOLD = 60;
int count_nonzero() const {
int cnt = 0;
for (const auto& x : *this) if (x != 0) cnt++;
return cnt;
}
FPS inv(int deg = -1) const {
if (deg == -1) deg = size();
if (count_nonzero() <= SPARSE_THRESHOLD) return sparse_inv(deg);
return dense_inv(deg);
}
FPS log(int deg = -1) const {
if (deg == -1) deg = size();
if (count_nonzero() <= SPARSE_THRESHOLD) return sparse_log(deg);
return dense_log(deg);
}
FPS exp(int deg = -1) const {
if (deg == -1) deg = size();
if (count_nonzero() <= SPARSE_THRESHOLD) return sparse_exp(deg);
return dense_exp(deg);
}
FPS pow(long long k, int deg = -1) const {
if (deg == -1) deg = size();
if (count_nonzero() <= SPARSE_THRESHOLD) return sparse_pow(k, deg);
return dense_pow(k, deg);
}
FPS sqrt(int deg = -1) const {
if (deg == -1) deg = size();
if (count_nonzero() <= SPARSE_THRESHOLD) return sparse_sqrt(deg);
return dense_sqrt(deg);
}
// --- 逆関数 (compositional inverse) F^{-1}(x) ---
// 制約: f[0] == 0, f[1] != 0
// 計算量: O(N log N) (2D Bostan-Mori, nachia's original algorithm)
FPS compositional_inverse(int deg = -1) const {
if (deg == -1) deg = size();
assert(deg > 1 && !empty() && (*this)[0] == 0 && (*this)[1] != 0);
mint t = (*this)[1].inv();
FPS f = pre(deg);
for (auto& x : f) x *= t;
int coeffAt = deg - 1;
int maxPower = deg - 1;
int n = 1;
while (n < std::max(coeffAt + 1, maxPower + 1)) n *= 2;
n *= 2;
std::vector<mint> p(n / 2, 0);
p[n / 2 - 1 - coeffAt] = 1;
std::vector<mint> q(n / 2, 0);
for (int i = 0; i < std::min((int)f.size(), n / 2); i++) q[i] = -f[i];
q[0] += 1;
std::vector<mint> p1(n, 0), p2(n * 2, 0), q1(n, 0), q2(n * 2, 0);
int d = n / 2, e = 1;
mint pcoeff = 1, invn = mint(n).inv();
auto bat = [&](std::vector<mint>& a, int s, int z, int r) {
NachiaNtt::ButterflyLayered(a.data(), z, s, r);
};
auto ibat = [&](std::vector<mint>& a, int s, int z, int r) {
NachiaNtt::IButterflyLayered(a.data(), z, s, r);
};
auto rotbuf = [&]() {
std::swap(p, p1); std::swap(p1, p2);
std::swap(q, q1); std::swap(q1, q2);
};
while (true) {
int fg = (d <= e ? 1 : 0);
int m = n / 2;
for (int iter = 0; iter < 2; iter++) {
if (fg ^ iter) {
ibat(q, d, e, 1);
ibat(p, d, e, 1);
if (d == 1) break;
for (int i = 0; i < m; i++) q1[i] = q[i];
for (int i = 0; i < m; i++) q1[m + i] = 0;
bat(q1, d, e * 2, 1);
for (int i = 0; i < m; i++) p1[i] = p[i];
for (int i = 0; i < m; i++) p1[m + i] = 0;
bat(p1, d, e * 2, 1);
e *= 2; m *= 2;
} else {
for (int i = 0; i < m * 2; i += d * 2) {
for (int j = 0; j < d; j++) q1[i + j] = q[i / 2 + j];
for (int j = 0; j < d; j++) q1[i + d + j] = 0;
}
for (int i = 0; i < m * 2; i += d * 2) {
for (int j = 0; j < d * 2; j++) p1[i + j] = 0;
for (int j = 0; j < d; j++) p1[i + j + 1] = p[i / 2 + j];
}
bat(q1, 1, d * 2, e); bat(p1, 1, d * 2, e);
d *= 2; m *= 2;
}
rotbuf();
}
if (d == 1) break;
pcoeff *= mint(n * 2);
rotbuf();
for (int i = 0; i < n; i++) q2[n + i] -= 2;
d /= 2;
for (int i = 0; i < n; i++) p1[i] = p2[i * 2] * q2[i * 2 + 1] + p2[i * 2 + 1] * q2[i * 2];
ibat(p1, 1, d, e);
for (int i = 0; i < n; i++) q1[i] = q2[i * 2] * q2[i * 2 + 1];
ibat(q1, 1, d, e);
d /= 2;
for (int i = 0; i < n; i += d * 2) for (int j = 0; j < d; j++) p[i / 2 + j] = p1[i + j + 1];
for (int i = 0; i < n; i += d * 2) for (int j = 0; j < d; j++) q[i / 2 + j] = q1[i + j];
for (auto& x : q) x *= invn;
}
FPS g(deg, 0);
mint inv_pcoeff = pcoeff.inv();
for (int i = 0; i <= maxPower && i < (int)p.size(); i++) {
if (p.size() - 1 - i < p.size()) {
g[i] = p[p.size() - 1 - i] * inv_pcoeff;
}
}
std::vector<mint> inv(deg);
if (deg > 1) inv[1] = 1;
for (int i = 2; i < deg; i++) inv[i] = -inv[mint::mod() % i] * (mint::mod() / i);
mint K = deg - 1;
for (int i = 1; i < deg; i++) g[i] *= K * inv[i];
std::reverse(g.begin(), g.end());
g = g.pow((-K).inv().val(), deg);
for (int i = deg - 1; i >= 1; i--) g[i] = g[i - 1];
mint tt = t;
for (int i = 1; i < deg; i++) { g[i] *= tt; tt *= t; }
g[0] = 0;
return g;
}
// --- 合成 (Composition) f(g(x)) ---
// 制約: deg(f) == deg(g) == N (内部でよしなに切り捨て/延長されます)
// 計算量: O(N log N) (2D Bostan-Mori / nachia's original algorithm)
FPS composition(FPS g, int deg = -1) const {
if (deg == -1) deg = size();
FPS f = this->pre(deg);
g = g.pre(deg);
if (deg == 0) return FPS();
// g[0] != 0 の場合は f(x) を f(x + g[0]) にテイラーシフトして g[0] = 0 に帰着
if (g[0] != 0) {
f = f.shift(g[0]);
g[0] = 0;
}
int n = 1; while (n < deg) n *= 2;
n *= 2;
std::vector<mint> q(n, 0);
for (int i = 0; i < std::min((int)g.size(), n / 2); i++) q[i] = -g[i];
q[0] += 1;
std::vector<mint> q2(n * 2, 0), q1(n, 0);
int d = n / 2;
int e = 1;
std::vector<std::vector<mint>> G;
mint invn = mint(n).inv();
auto bat = [&](std::vector<mint>& a, int s, int z, int r) {
NachiaNtt::ButterflyLayered(a.data(), z, s, r);
};
auto ibat = [&](std::vector<mint>& a, int s, int z, int r) {
NachiaNtt::IButterflyLayered(a.data(), z, s, r);
};
auto rotbuf = [&]() {
std::swap(q, q1); std::swap(q1, q2);
};
while (true) {
int fg = (d <= e ? 1 : 0);
int m = n / 2;
for (int t = 0; t < 2; t++) {
if (fg ^ t) {
ibat(q, d, e, 1);
if (d == 1) break;
for (int i = 0; i < d; i++) q1[i] = q[i];
for (int i = 0; i < m; i++) q1[d + i] = 0;
for (int i = d; i < m; i++) q1[m + i] = q[i];
bat(q1, d, e * 2, 1);
e *= 2; m *= 2;
} else {
for (int i = 0; i < m * 2; i += d * 2) {
for (int j = 0; j < d; j++) q1[i + j] = q[i / 2 + j];
for (int j = 0; j < d; j++) q1[i + d + j] = 0;
}
bat(q1, 1, d * 2, e);
d *= 2; m *= 2;
}
rotbuf();
}
if (d == 1) break;
for (int i = 0; i < n; i++) q[n + i] -= 2;
for (int i = 0; i < n; i++) std::swap(q[i * 2], q[i * 2 + 1]);
std::vector<mint> q_copy = q;
q_copy.resize(n * 2, 0);
G.push_back(std::move(q_copy));
for (int x = 0; x < n * 2; x += d) {
for (int y = 1; y < d; y *= 2) {
std::reverse(G.back().begin() + x + y, G.back().begin() + x + y * 2);
}
}
rotbuf();
d /= 2;
for (int i = 0; i < n; i++) q1[i] = q2[i * 2] * q2[i * 2 + 1];
ibat(q1, 1, d, e);
d /= 2;
for (int i = 0; i < n; i += d * 2) {
for (int j = 0; j < d; j++) q[i / 2 + j] = q1[i + j];
}
for (auto& x : q) x *= invn;
}
// 後退代入用に配列サイズを整える
q.assign(n * 2, 0);
q1.assign(n, 0);
q2.assign(n * 2, 0);
for (int i = 0; i < n / 2; i++) {
q[n / 2 - 1 - i] = (i < (int)f.size() ? f[i] : mint(0));
}
while (G.size()) {
bat(q, d, e, 1);
for (int i = 0; i < n; i += d * 2) {
q1[i] = 0;
for (int j = 0; j < d; j++) q1[i + j + 1] = q[i / 2 + j];
for (int j = 1; j < d; j++) q1[i + d + j] = 0;
}
d *= 2;
bat(q1, 1, d, e);
auto qnntt = std::move(G.back()); G.pop_back();
for (int i = 0; i < n; i++) {
q2[i * 2] = q1[i] * qnntt[i * 2];
q2[i * 2 + 1] = q1[i] * qnntt[i * 2 + 1];
}
mint inv2 = mint(2).inv() * invn;
for (int i = 0; i < n * 2; i++) q2[i] *= inv2;
ibat(q2, 1, d * 2, e);
for (int i = 0; i < n * 2; i += d * 2) {
for (int j = 0; j < d; j++) q1[i / 2 + j] = q2[i + j + 1];
}
ibat(q1, d, e, 1);
e /= 2;
for (int i = 0; i < n / 2; i++) q[i] = q1[i];
}
FPS res(deg, 0);
for (int i = 0; i < deg; i++) {
res[i] = q[n / 2 - 1 - i];
}
return res;
}
};
void solve(){
ll n,m;
cin >> n >> m;
vl a(n),b(n),c(n);
rep(i,n)cin >> a[i];
rep(i,n)cin >> b[i];
rep(i,n)cin >> c[i];
FPS f(all(a)),g(all(b)),h(all(c));
// FPS ans=f;
// rep(i,m-1)ans=ans.composition(f);
FPS ag=g*mint(a[1]).pow(m);
// FPS gf=g.composition(ans);
// FPS gh=g.composition(h);
// rep(i,n)cout << ans[i].val() << " \n"[i==n-1];
// ag.dump();
// gf.dump();
// gh.dump();
// FPS hgf=h.composition(gf);
FPS hag=h.composition(ag);
// hgf.dump();
// hag.dump();
rep(i,n)cout << hag[i].val() << " \n"[i==n-1];
}
int main(){
ios::sync_with_stdio(false);
std::cin.tie(nullptr);
// ll mx=450;
// vc fl(mx+1,0);
// for(ll d=2;d<=mx;d++){
// if(fl[d])continue;
// ll x=d;
// ps.push_back(x);
// while(x<=mx){
// fl[x]=1;
// x+=d;
// }
// }
ll t = 1;
// cin >> t;
while (t--){
solve();
}
}