結果
| 問題 | No.3619 Compositional Power with Schröder Coordinate |
| コンテスト | |
| ユーザー |
 2251799813685248
|
| 提出日時 | 2026-08-12 00:07:02 |
| 言語 | C++23(gcc16) (gcc 16.1.0 + boost 1.90.0) |
| 結果 |
AC
|
| 実行時間 | 1,738 ms / 10,000 ms |
| + 961µs | |
| コード長 | 38,699 bytes |
| 記録 | |
| コンパイル時間 | 7,160 ms |
| コンパイル使用メモリ | 488,080 KB |
| 実行使用メモリ | 104,000 KB |
| 最終ジャッジ日時 | 2026-08-12 00:07:23 |
| 合計ジャッジ時間 | 17,931 ms |
|
ジャッジサーバーID （参考情報） |
judge2_1 / judge3_1 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 2 |
| other | AC * 6 |
ソースコード
#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include <unordered_set>
#include <unordered_map>
#include <queue>
#include <algorithm>
#include <iomanip>
#include <cassert>
#include <functional>
#include <random>
#include <bitset>
#include <unistd.h>
using namespace std;
using ll = long long;
using lll = __int128_t;
using ull = unsigned long long;
using ld = long double;
using pii = array<int,2>;
using pll = array<ll,2>;
using plll = array<lll,2>;
#define vall(A) A.begin(), A.end()
template<typename T> inline void vin(T& A){for (int i = 0, sz = A.size(); i < sz; i++){cin >> A[i];}}
template<typename T> inline void vout(const T& A){for (int i = 0, sz = A.size(); i < sz; i++){cout << A[i] << " \n"[i == sz-1];}}
template<typename T> inline void vout2d(const T& A){for (int i = 0, H = A.size(); i < H; i++){vout(A[i]);}}
template<typename T> inline void adjvin(T& A){for (int i = 1, sz = A.size(); i < sz; i++){cin >> A[i];}}
template<typename T> inline void adjvout(const T& A){for (int i = 1, sz = A.size(); i < sz; i++){cout << A[i] << " \n"[i == sz-1];}}
template<typename T> inline void adjvout2d(const T& A){for (int i = 1, H = A.size(); i < H; i++){adjvout(A[i]);}}
template<typename T> inline bool btest(T K, int i){return K&(1ull<<i);}
template<typename T> void print(T object){cout << (object) << "\n";}
template<typename T, typename U> void print(T object1, U object2){cout << (object1) << " " << (object2) << "\n";}
template<typename T, typename U, typename V> void print(T object1, U object2, V object3){cout << (object1) << " " << (object2) << " " << (object3) << "\n";}
template<typename T, typename U, typename V, typename W> void print(T object1, U object2, V object3, W object4){cout << (object1) << " " << (object2) << " " << (object3) << " " << (object4) << "\n";}
const vector<ull> pow2ll{1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,32768,65536,131072,262144,524288,1048576,2097152,4194304,8388608,16777216,33554432,67108864,134217728,268435456,536870912,1073741824,2147483648,4294967296,8589934592,17179869184,34359738368,68719476736,137438953472,274877906944,549755813888,1099511627776,2199023255552,4398046511104,8796093022208,17592186044416,35184372088832,70368744177664,140737488355328,281474976710656,562949953421312,1125899906842624,2251799813685248,4503599627370496,9007199254740992,18014398509481984,36028797018963968,72057594037927936,144115188075855872,288230376151711744,576460752303423488,1152921504606846976,2305843009213693952,4611686018427387904, 9223372036854775808ull};
const vector<ull> pow10ll{1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,10000000000,100000000000,1000000000000,10000000000000,100000000000000,1000000000000000,10000000000000000,100000000000000000,1000000000000000000, 10000000000000000000ull};
const vector<ll> di{0,1,0,-1};
const vector<ll> di8{0,1,1,1,0,-1,-1,-1};
const vector<ll> dj{1,0,-1,0};
const vector<ll> dj8{1,1,0,-1,-1,-1,0,1};
#ifndef INV_TABLE_HPP_
#define INV_TABLE_HPP_
#include <vector>
#include <cassert>
using uint = unsigned;
using ull = unsigned long long;
using ll = long long;
using namespace std;
/// @brief mod M上での階乗,逆元テーブルを保持する構造体
/// @tparam M
template<uint M>
struct mod_table{
vector<uint> invmodlist;
vector<uint> factorialmodlist;
vector<uint> factorialmodinvlist;
uint N_MAX;
constexpr mod_table(const uint N_MAX__){
N_MAX = max(1u, N_MAX__);
invmodlist = vector<uint>(N_MAX+1);
factorialmodlist = vector<uint>(N_MAX+1);
factorialmodinvlist = vector<uint>(N_MAX+1);
invmodlist[1] = 1;
for (uint i = 2; i <= N_MAX; i++){
invmodlist[i] = (M-M/i)*(ull)invmodlist[M%i]%M;
}
factorialmodinvlist[0] = 1;
factorialmodlist[0] = 1;
for (ull i = 1; i <= N_MAX; i++){
factorialmodinvlist[i] = (invmodlist[i]*(ull)factorialmodinvlist[i-1])%M;
factorialmodlist[i] = (factorialmodlist[i-1]*i)%M;
}
}
ll inverse_mod(ll x){
assert(0 <= x and x < N_MAX);
return invmodlist[x];
}
ll factorialmod(ll x){
assert(0 <= x and x < N_MAX);
return factorialmodlist[x];
}
ll factorialmodinv(ll x){
assert(0 <= x and x < N_MAX);
return factorialmodinvlist[x];
}
};
#endif /* INV_TABLE_HPP_ */
#ifndef FPS__HPP_
#define FPS__HPP_
#include <iostream>
#include <iterator>
#include <vector>
#include <cassert>
#ifndef MODINT__HPP_
#define MODINT__HPP_
#include <iostream>
#ifndef MATH_FUNCTION_HPP_
#define MATH_FUNCTION_HPP_
#include <array>
#include <cmath>
#include <cassert>
using namespace std;
using ll = long long;
using ld = long double;
using ull = unsigned long long;
using ulll = __uint128_t;
#ifdef BOOST_VERSION
using bll = boost::multiprecision::cpp_int;
#endif
/// @brief a^bをmで割った余りを返す。bに関して対数時間で計算できる
constexpr ll modpow(ll a, ull b, const ll m){
ll t = a%m;
ll ans = (m == 1 ? 0 : 1);
while (b > 0){
if (b&1){
ans = (ans*t)%m;
}
b >>= 1;
t = (t*t)%m;
}
return ans;
}
/// @brief a^bをmで割った余りを返す。bに関して対数時間で計算できる。mはコンパイル時に決定している必要がある
template<ll m> constexpr ll modpow(ll a, ull b){
ll t = a%m;
ll ans = (m == 1 ? 0 : 1);
while (b > 0){
if (b&1){
ans = (ans*t)%m;
}
b >>= 1;
t = (t*t)%m;
}
return ans;
}
/// @brief a^nを返す。bに関して対数時間で計算できる。
template<typename T> constexpr T powll(T a, T n){
T t = a;
T ans = 1;
while (n > 0){
if (n%2){
ans *= t;
}
n >>= 1;
t *= t;
}
return ans;
}
/// @brief a^nを返す。bに関して対数時間で計算できる
template<> constexpr ll powll(ll a, ll n){
ll t = a;
ll ans = 1;
while (n > 0){
if (n&1){
ans *= t;
}
n >>= 1;
t *= t;
}
return ans;
}
/// @brief floor(sqrt(N))を返す。1.5×10^19まで対応
constexpr ll isqrt(ull N){
assert(N <= 15000000000000000000ull);
ull ret = sqrt(N);
while (ret*ret > N){
ret--;
}
while ((ret+1)*(ret+1) <= N){
ret++;
}
return ret;
}
#ifdef BOOST_VERSION
/// @brief floor(sqrt(N))を返す。多倍長整数に対応
constexpr bll isqrt_large(bll N){
bll ret = sqrt(N);
while (ret*ret > N){
ret--;
}
while ((ret+1)*(ret+1) <= N){
ret++;
}
return ret;
}
#endif
/// @brief floor(log_a(L))を返す
constexpr ll ilog(ll a, ll L){
__int128_t t = 1;
ll ans = 0;
while (t <= L){
ans++;
t *= a;
}
return ans-1;
}
/// @brief 有理数のfloorを求める。 floor(y/x)
template<typename T> constexpr inline T floor2(T y, T x){
if ((x^y) > 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return y/x;
}
else if ((x^y) < 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return -((y+x-1)/x);
}
else{
return y/x;
}
}
/// @brief 有理数のceilを求める。 ceil(y/x)
template<typename T> constexpr inline T ceil2(T y, T x){
if ((x^y) > 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return (y+x-1)/x;
}
else if ((x^y) < 0){
x = x > 0 ? x : -x;
y = y > 0 ? y : -y;
return -(y/x);
}
else{
return y/x;
}
}
/// @brief 一次不定方程式ax+by=gcd(a,b)の解を1つ見つける
/// @param a `a>=0`である必要がある
/// @param b `b>=0`である必要がある
/// @return {x,y,gcd(a,b)}
template<typename T>
constexpr array<T,3> axby1(T a, T b){
T x = 1, y = 0;
T z = 0, w = 1;
T tmp = 0;
while (b){
T p = a/b, q = a%b;
tmp = x - y * p; x = y; y = tmp;
tmp = z - w * p; z = w; w = tmp;
a = b; b = q;
}
return {x, z, a};
}
/// @brief 1/a mod Mを求める
template<typename T, typename U>
constexpr T inverse_mod(T a, U M){
auto temp = axby1(a,(T)M);
assert(temp[2] == 1);
return (M+temp[0])%M;
}
/// @brief sqrt(a) mod Mを求める。ないなら-1が返される。
template<ll M>
inline constexpr ll cipolla(ll a){
a %= M;
if (M == 2) return a;
if (a == 0) return 0;
ll z = (M-1)/2;
if (modpow<M>(a, z) != 1){return -1;}
int b = 0;
while (modpow<M>((b*b+M-a)%M, z) == 1){
b++;
}
array<ll,2> x{1,0};
array<ll,2> y{b, 1};
ll w = (b*b+M-a)%M;
z++;
while (z){
if (z&1){
ll temp = x[0];
x[0] = x[0]*y[0]%M+x[1]*y[1]%M*w%M;
if (x[0] >= M){x[0] -= M;}
x[1] = temp*y[1]%M+x[1]*y[0]%M;
if (x[1] >= M){x[1] -= M;}
}
ll temp = y[0];
y[0] = y[0]*y[0]%M+y[1]*y[1]%M*w%M;
if (y[0] >= M){y[0] -= M;}
y[1] = 2*temp*y[1]%M;
z >>= 1;
}
return x[0];
}
inline constexpr ll cipolla(ll a, const ll M){
a %= M;
if (M == 2) return a;
if (a == 0) return 0;
ll z = (M-1)/2;
if (modpow(a, z, M) != 1){return -1;}
int b = 0;
while (modpow((b*b+M-a)%M, z, M) == 1){
b++;
}
array<ll,2> x{1,0};
array<ll,2> y{b, 1};
ll w = (b*b+M-a)%M;
z++;
while (z){
if (z&1){
ll temp = x[0];
x[0] = x[0]*y[0]%M+x[1]*y[1]%M*w%M;
if (x[0] >= M){x[0] -= M;}
x[1] = temp*y[1]%M+x[1]*y[0]%M;
if (x[1] >= M){x[1] -= M;}
}
ll temp = y[0];
y[0] = y[0]*y[0]%M+y[1]*y[1]%M*w%M;
if (y[0] >= M){y[0] -= M;}
y[1] = 2*temp*y[1]%M;
z >>= 1;
}
return x[0];
}
/// @brief x以下の最大の2冪を返す。0は0が返る。
constexpr ull lowerpow2(ull x){
if (x == 0){return 0;}
return 1ull<<(63-__builtin_clzll(x));
}
/// @brief x以上の最小の2冪を返す。0は0が返る。
constexpr ull upperpow2(ull x){
if (x == 0){return 0;}
if (x == 1){return 1;}
return 1ull<<(64-__builtin_clzll(x-1));
}
/// @brief xのpopcountを求める
constexpr int popcount(ull x){
return __builtin_popcountll(x);
}
/// @brief xのbit lengthを求める。
constexpr int bit_length(ull x){
if (x == 0){return 0;}
return 64-__builtin_clzll(x);
}
#endif /* MATH_FUNCTION_HPP_ */
using ll = long long;
using uint = unsigned int;
using ull = unsigned long long;
using namespace std;
template <uint M>
struct constant_modint {
uint val;
constant_modint() : val(0) {}
template <class T> constant_modint(T x) {
ll y = (ll)(x % (ll)M);
if (y < 0) y += M;
val = y;
}
constant_modint& operator+=(const constant_modint& rhs) {
val += rhs.val;
if (val >= M) val -= M;
return *this;
}
constant_modint& operator-=(const constant_modint& rhs) {
if (val < rhs.val) val += M;
val -= rhs.val;
return *this;
}
constant_modint& operator*=(const constant_modint& rhs) {
val = (ull)val * rhs.val % M;
return *this;
}
constant_modint& operator/=(const constant_modint& rhs) {
return *this *= rhs.inv();
}
constant_modint operator+() const { return *this; }
constant_modint operator-() const { return constant_modint(0) - *this; }
friend constant_modint operator+(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) += rhs; }
friend constant_modint operator-(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) -= rhs; }
friend constant_modint operator*(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) *= rhs; }
friend constant_modint operator/(const constant_modint& lhs, const constant_modint& rhs) { return constant_modint(lhs) /= rhs; }
friend bool operator==(const constant_modint& lhs, const constant_modint& rhs) { return lhs.val == rhs.val; }
friend bool operator!=(const constant_modint& lhs, const constant_modint& rhs) { return lhs.val != rhs.val; }
constant_modint pow(ull n) const {
constant_modint res = 1, a = *this;
while (n) {
if (n&1) res *= a;
a *= a;
n >>= 1;
}
return res;
}
constant_modint inv() const {
return inverse_mod((ll)val, M);
}
friend std::ostream& operator<<(std::ostream& os, const constant_modint& m) {
return os << m.val;
}
friend std::istream& operator>>(std::istream& is, constant_modint& m) {
ll x;
is >> x;
m = constant_modint(x);
return is;
}
};
template <int id = -1>
struct dynamic_modint {
uint val;
static uint& mod() {
static uint M = 998244353;
return M;
}
static void set_mod(uint m) {
mod() = m;
}
dynamic_modint() : val(0) {}
template <class T> dynamic_modint(T x) {
ll y = x % (ll)mod();
if (y < 0) y += mod();
val = (uint)y;
}
dynamic_modint& operator+=(const dynamic_modint& rhs) {
val += rhs.val;
if (val >= mod()) val -= mod();
return *this;
}
dynamic_modint& operator-=(const dynamic_modint& rhs) {
if (val < rhs.val) val += mod();
val -= rhs.val;
return *this;
}
dynamic_modint& operator*=(const dynamic_modint& rhs) {
val = (ull)val * rhs.val % mod();
return *this;
}
dynamic_modint& operator/=(const dynamic_modint& rhs) {
return *this *= rhs.inv();
}
dynamic_modint operator+() const { return *this; }
dynamic_modint operator-() const { return dynamic_modint(0) - *this; }
friend dynamic_modint operator+(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) += rhs; }
friend dynamic_modint operator-(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) -= rhs; }
friend dynamic_modint operator*(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) *= rhs; }
friend dynamic_modint operator/(const dynamic_modint& lhs, const dynamic_modint& rhs) { return dynamic_modint(lhs) /= rhs; }
friend bool operator==(const dynamic_modint& lhs, const dynamic_modint& rhs) { return lhs.val == rhs.val; }
friend bool operator!=(const dynamic_modint& lhs, const dynamic_modint& rhs) { return lhs.val != rhs.val; }
dynamic_modint pow(ull n) const {
dynamic_modint res = 1, a = *this;
while (n) {
if (n & 1) res *= a;
a *= a;
n >>= 1;
}
return res;
}
dynamic_modint inv() const {
ll a = val, b = mod(), u = 1, v = 0;
while (b) {
ll t = a / b;
a -= t * b; swap(a, b);
u -= t * v; swap(u, v);
}
return dynamic_modint(u);
}
friend std::ostream& operator<<(std::ostream& os, const dynamic_modint& m) {
return os << m.val;
}
friend std::istream& operator>>(std::istream& is, dynamic_modint& m) {
ll x;
is >> x;
m = dynamic_modint(x);
return is;
}
};
#endif /* MODINT__HPP_ */
using namespace std;
using ll = long long;
using uint = unsigned;
template<uint T>
struct FormalPowerSeries{
using mint = constant_modint<T>;
vector<mint> f;
uint sz;
FormalPowerSeries(){
sz = 0;
}
/// @brief 指定したサイズのfpsを作成する。
FormalPowerSeries(uint _init_sz){
f.resize(_init_sz, 0);
sz = _init_sz;
}
FormalPowerSeries(uint _init_sz, mint _init_val){
f.resize(_init_sz, _init_val);
sz = _init_sz;
}
template<typename U> FormalPowerSeries(const vector<U>& _init){
f.assign(_init.begin(), _init.end());
sz = f.size();
}
template<typename U> FormalPowerSeries(initializer_list<U> _init) : FormalPowerSeries(vector<U>(_init)){}
uint size() const {return sz;}
void resize(uint _new_size){
if (sz == _new_size){return;}
f.resize(_new_size, 0);
sz = _new_size;
}
inline mint& operator[](uint deg) {return f[deg];}
inline const mint& operator[](uint deg) const {return f[deg];}
template<typename U> operator vector<U>() const {
vector<U> res(sz);
for (uint i = 0; i < sz; i++){
res[i] = f[i].val;
}
return res;
}
void differential(){
if (sz == 0){return;}
for (uint i = 0; i < sz-1; i++){
f[i].val = (i+1)*(ull)f[i+1].val%T;
}
f.pop_back();
sz--;
}
void integral(const mod_table<T>& mtable){
if (sz == 0){return;}
sz++;
f.push_back(0);
for (uint i = sz; i > 0; i--){
f[i].val = mtable.invmodlist[i]*(ull)f[i-1].val%T;
}
}
};
#endif /* FPS__HPP_ */
#ifndef FPS_OPERATION_HPP_
#define FPS_OPERATION_HPP_
#include <vector>
#include <algorithm>
#include <immintrin.h>
using namespace std;
using ll = long long;
using uint = unsigned;
using ull = unsigned long long;
template<typename T> inline bool btest_for_fps(T K, int i){return K&(1ull<<i);}
/// @brief mod M上での形式的冪級数の計算を行う構造体
template<uint M>
struct fps_operator{
uint sum_e[30];
uint sum_ie[30];
uint log_max_length;
uint last_powroot;
// --- モンゴメリ乗算用の定数と関数 ---
static constexpr uint get_r() {
uint res = M;
for (int i = 0; i < 4; ++i) res *= 2 - M * res;
return res;
}
static constexpr uint R = get_r();
static constexpr uint R2 = -ull(M) % M;
static constexpr uint NEG_INV = 0 - R;
static inline constexpr uint reduce(ull x) {
uint res = (x + ull(uint(x) * NEG_INV) * M) >> 32;
return res >= M ? res - M : res;
}
static inline constexpr uint to_montgomery(uint x) {
return reduce(ull(x) * R2);
}
static inline constexpr uint from_montgomery(uint x) {
return reduce(x);
}
static inline constexpr uint montgomery_mul(uint x, uint y) {
return reduce(ull(x) * y);
}
static inline constexpr uint mod_add(uint x, uint y) {
return x + y >= M ? x + y - M : x + y;
}
static inline constexpr uint mod_sub(uint x, uint y) {
return x < y ? x + M - y : x - y;
}
// --- AVX-512 / AVX2 SIMD Intrinsics ---
#if defined(__AVX512F__) && defined(__AVX512DQ__)
static inline __m512i montgomery_mul_simd(__m512i a, __m512i b, __m512i mod, __m512i neg_inv) {
__m512i mul0 = _mm512_mul_epu32(a, b);
__m512i mul1 = _mm512_mul_epu32(_mm512_srli_epi64(a, 32), _mm512_srli_epi64(b, 32));
__m512i q0 = _mm512_mul_epu32(_mm512_mullo_epi32(mul0, neg_inv), mod);
__m512i q1 = _mm512_mul_epu32(_mm512_mullo_epi32(mul1, neg_inv), mod);
__m512i res0 = _mm512_srli_epi64(_mm512_add_epi64(mul0, q0), 32);
__m512i res1 = _mm512_srli_epi64(_mm512_add_epi64(mul1, q1), 32);
__m512i res = _mm512_or_si512(res0, _mm512_slli_epi64(res1, 32));
__mmask16 cmp = _mm512_cmpge_epu32_mask(res, mod);
return _mm512_mask_sub_epi32(res, cmp, res, mod);
}
static inline __m512i mod_add_simd(__m512i a, __m512i b, __m512i mod) {
__m512i res = _mm512_add_epi32(a, b);
__mmask16 cmp = _mm512_cmpge_epu32_mask(res, mod);
return _mm512_mask_sub_epi32(res, cmp, res, mod);
}
static inline __m512i mod_sub_simd(__m512i a, __m512i b, __m512i mod) {
__mmask16 cmp = _mm512_cmplt_epu32_mask(a, b);
__m512i res = _mm512_sub_epi32(a, b);
return _mm512_mask_add_epi32(res, cmp, res, mod);
}
#elif defined(__AVX2__)
// AVX2用 モンゴメリ乗算・剰余加減算
static inline __m256i montgomery_mul_simd(__m256i a, __m256i b, __m256i mod, __m256i neg_inv) {
__m256i mul0 = _mm256_mul_epu32(a, b);
__m256i mul1 = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
__m256i q0 = _mm256_mul_epu32(_mm256_mullo_epi32(mul0, neg_inv), mod);
__m256i q1 = _mm256_mul_epu32(_mm256_mullo_epi32(mul1, neg_inv), mod);
__m256i res0 = _mm256_srli_epi64(_mm256_add_epi64(mul0, q0), 32);
__m256i res1 = _mm256_srli_epi64(_mm256_add_epi64(mul1, q1), 32);
__m256i res = _mm256_blend_epi32(res0, _mm256_slli_epi64(res1, 32), 0xAA);
__m256i diff = _mm256_sub_epi32(res, mod);
// diff < 0 なら -1 (全ビット1)、そうでないなら 0 となるマスク
__m256i mask = _mm256_srai_epi32(diff, 31);
return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
}
static inline __m256i mod_add_simd(__m256i a, __m256i b, __m256i mod) {
__m256i diff = _mm256_sub_epi32(_mm256_add_epi32(a, b), mod);
__m256i mask = _mm256_srai_epi32(diff, 31);
return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
}
static inline __m256i mod_sub_simd(__m256i a, __m256i b, __m256i mod) {
__m256i diff = _mm256_sub_epi32(a, b);
__m256i mask = _mm256_srai_epi32(diff, 31);
return _mm256_add_epi32(diff, _mm256_and_si256(mask, mod));
}
#endif
constexpr fps_operator(){
vector<ll> powroot{1};
vector<ll> powrootinv;
while (powroot.back() >= 0){
powroot.push_back(cipolla<M>(powroot.back()));
if (powroot.back() == powroot[powroot.size()-2]){
powroot.back() = M-powroot.back();
}
}
powroot.pop_back();
log_max_length = powroot.size()-1;
last_powroot = powroot.back();
for (auto v : powroot){
powrootinv.push_back(inverse_mod<ll,ll>(v,M));
}
int cnt2 = powroot.size()-1;
uint now = reduce(R2);
for (int i = 0; i <= cnt2-2; i++){
sum_e[i] = montgomery_mul(to_montgomery(powroot[i+2]), now);
now = montgomery_mul(now, to_montgomery(powrootinv[i+2]));
}
uint inow = reduce(R2);
for (int i = 0; i <= cnt2-2; i++){
sum_ie[i] = montgomery_mul(to_montgomery(powrootinv[i+2]), inow);
inow = montgomery_mul(inow, to_montgomery(powroot[i+2]));
}
for (int i = cnt2-1; i < 30; i++){
sum_e[i] = 0;
sum_ie[i] = 0;
}
}
void inplaceDFT(FormalPowerSeries<M>& F) const {
F.resize(upperpow2(F.sz));
int n = F.sz;
if (n == 0) return;
int h = __builtin_ctz(n);
for (int ph = 1; ph <= h; ph++) {
int w = 1 << (ph - 1), p = 1 << (h - ph);
uint now = to_montgomery(1);
for (int s = 0; s < w; s++) {
int offset = s << (h - ph + 1);
#if defined(__AVX512F__) && defined(__AVX512DQ__)
if (p >= 16) {
__m512i vmod = _mm512_set1_epi32(M);
__m512i vneg_inv = _mm512_set1_epi32(NEG_INV);
__m512i vnow = _mm512_set1_epi32(now);
for (int i = 0; i < p; i += 16) {
__m512i l = _mm512_loadu_si512((__m512i*)&F[i + offset].val);
__m512i r = _mm512_loadu_si512((__m512i*)&F[i + offset + p].val);
__m512i r_now = montgomery_mul_simd(r, vnow, vmod, vneg_inv);
_mm512_storeu_si512((__m512i*)&F[i + offset].val, mod_add_simd(l, r_now, vmod));
_mm512_storeu_si512((__m512i*)&F[i + offset + p].val, mod_sub_simd(l, r_now, vmod));
}
} else
#elif defined(__AVX2__)
if (p >= 8) {
__m256i vmod = _mm256_set1_epi32(M);
__m256i vneg_inv = _mm256_set1_epi32(NEG_INV);
__m256i vnow = _mm256_set1_epi32(now);
for (int i = 0; i < p; i += 8) {
__m256i l = _mm256_loadu_si256((__m256i*)&F[i + offset].val);
__m256i r = _mm256_loadu_si256((__m256i*)&F[i + offset + p].val);
__m256i r_now = montgomery_mul_simd(r, vnow, vmod, vneg_inv);
_mm256_storeu_si256((__m256i*)&F[i + offset].val, mod_add_simd(l, r_now, vmod));
_mm256_storeu_si256((__m256i*)&F[i + offset + p].val, mod_sub_simd(l, r_now, vmod));
}
} else
#endif
{
for (int i = 0; i < p; i++) {
uint l = F[i + offset].val;
uint r = montgomery_mul(F[i + offset + p].val, now);
F[i + offset].val = mod_add(l, r);
F[i + offset + p].val = mod_sub(l, r);
}
}
now = montgomery_mul(now, sum_e[__builtin_ctz(~s)]);
}
}
}
void inplaceIDFT(FormalPowerSeries<M>& F) const {
F.resize(upperpow2(F.sz));
int n = F.sz;
if (n == 0) return;
int h = __builtin_ctz(n);
for (int ph = h; ph >= 1; ph--) {
int w = 1 << (ph - 1), p = 1 << (h - ph);
uint inow = to_montgomery(1);
for (int s = 0; s < w; s++) {
int offset = s << (h - ph + 1);
#if defined(__AVX512F__) && defined(__AVX512DQ__)
if (p >= 16) {
__m512i vmod = _mm512_set1_epi32(M);
__m512i vneg_inv = _mm512_set1_epi32(NEG_INV);
__m512i vinow = _mm512_set1_epi32(inow);
for (int i = 0; i < p; i += 16) {
__m512i l = _mm512_loadu_si512((__m512i*)&F[i + offset].val);
__m512i r = _mm512_loadu_si512((__m512i*)&F[i + offset + p].val);
_mm512_storeu_si512((__m512i*)&F[i + offset].val, mod_add_simd(l, r, vmod));
_mm512_storeu_si512((__m512i*)&F[i + offset + p].val, montgomery_mul_simd(mod_sub_simd(l, r, vmod), vinow, vmod, vneg_inv));
}
} else
#elif defined(__AVX2__)
if (p >= 8) {
__m256i vmod = _mm256_set1_epi32(M);
__m256i vneg_inv = _mm256_set1_epi32(NEG_INV);
__m256i vinow = _mm256_set1_epi32(inow);
for (int i = 0; i < p; i += 8) {
__m256i l = _mm256_loadu_si256((__m256i*)&F[i + offset].val);
__m256i r = _mm256_loadu_si256((__m256i*)&F[i + offset + p].val);
_mm256_storeu_si256((__m256i*)&F[i + offset].val, mod_add_simd(l, r, vmod));
_mm256_storeu_si256((__m256i*)&F[i + offset + p].val, montgomery_mul_simd(mod_sub_simd(l, r, vmod), vinow, vmod, vneg_inv));
}
} else
#endif
{
for (int i = 0; i < p; i++) {
uint l = F[i + offset].val;
uint r = F[i + offset + p].val;
F[i + offset].val = mod_add(l, r);
F[i + offset + p].val = montgomery_mul(mod_sub(l, r), inow);
}
}
inow = montgomery_mul(inow, sum_ie[__builtin_ctz(~s)]);
}
}
}
void inplaceDFT_T(FormalPowerSeries<M>& F) const {
F.resize(upperpow2(F.sz));
int n = F.sz;
if (n == 0) return;
int h = __builtin_ctz(n);
for (int ph = h; ph >= 1; ph--) {
int w = 1 << (ph - 1), p = 1 << (h - ph);
uint32_t now = to_montgomery(1);
for (int s = 0; s < w; s++) {
int offset = s << (h - ph + 1);
#if defined(__AVX512F__) && defined(__AVX512DQ__)
if (p >= 16) {
__m512i vmod = _mm512_set1_epi32(M);
__m512i vneg_inv = _mm512_set1_epi32(NEG_INV);
__m512i vnow = _mm512_set1_epi32(now);
for (int i = 0; i < p; i += 16) {
__m512i l = _mm512_loadu_si512((__m512i*)&F[i + offset].val);
__m512i r = _mm512_loadu_si512((__m512i*)&F[i + offset + p].val);
_mm512_storeu_si512((__m512i*)&F[i + offset].val, mod_add_simd(l, r, vmod));
_mm512_storeu_si512((__m512i*)&F[i + offset + p].val, montgomery_mul_simd(mod_sub_simd(l, r, vmod), vnow, vmod, vneg_inv));
}
} else
#elif defined(__AVX2__)
if (p >= 8) {
__m256i vmod = _mm256_set1_epi32(M);
__m256i vneg_inv = _mm256_set1_epi32(NEG_INV);
__m256i vnow = _mm256_set1_epi32(now);
for (int i = 0; i < p; i += 8) {
__m256i l = _mm256_loadu_si256((__m256i*)&F[i + offset].val);
__m256i r = _mm256_loadu_si256((__m256i*)&F[i + offset + p].val);
_mm256_storeu_si256((__m256i*)&F[i + offset].val, mod_add_simd(l, r, vmod));
_mm256_storeu_si256((__m256i*)&F[i + offset + p].val, montgomery_mul_simd(mod_sub_simd(l, r, vmod), vnow, vmod, vneg_inv));
}
} else
#endif
{
for (int i = 0; i < p; i++) {
uint32_t l = F[i + offset].val;
uint32_t r = F[i + offset + p].val;
F[i + offset].val = mod_add(l, r);
F[i + offset + p].val = montgomery_mul(mod_sub(l, r), now);
}
}
now = montgomery_mul(now, sum_e[__builtin_ctz(~s)]);
}
}
}
void inplaceIDFT_T(FormalPowerSeries<M>& F) const {
F.resize(upperpow2(F.sz));
int n = F.sz;
if (n == 0) return;
int h = __builtin_ctz(n);
for (int ph = 1; ph <= h; ph++) {
int w = 1 << (ph - 1), p = 1 << (h - ph);
uint32_t inow = to_montgomery(1);
for (int s = 0; s < w; s++) {
int offset = s << (h - ph + 1);
#if defined(__AVX512F__) && defined(__AVX512DQ__)
if (p >= 16) {
__m512i vmod = _mm512_set1_epi32(M);
__m512i vneg_inv = _mm512_set1_epi32(NEG_INV);
__m512i vinow = _mm512_set1_epi32(inow);
for (int i = 0; i < p; i += 16) {
__m512i l = _mm512_loadu_si512((__m512i*)&F[i + offset].val);
__m512i r = _mm512_loadu_si512((__m512i*)&F[i + offset + p].val);
__m512i r_inow = montgomery_mul_simd(r, vinow, vmod, vneg_inv);
_mm512_storeu_si512((__m512i*)&F[i + offset].val, mod_add_simd(l, r_inow, vmod));
_mm512_storeu_si512((__m512i*)&F[i + offset + p].val, mod_sub_simd(l, r_inow, vmod));
}
} else
#elif defined(__AVX2__)
if (p >= 8) {
__m256i vmod = _mm256_set1_epi32(M);
__m256i vneg_inv = _mm256_set1_epi32(NEG_INV);
__m256i vinow = _mm256_set1_epi32(inow);
for (int i = 0; i < p; i += 8) {
__m256i l = _mm256_loadu_si256((__m256i*)&F[i + offset].val);
__m256i r = _mm256_loadu_si256((__m256i*)&F[i + offset + p].val);
__m256i r_inow = montgomery_mul_simd(r, vinow, vmod, vneg_inv);
_mm256_storeu_si256((__m256i*)&F[i + offset].val, mod_add_simd(l, r_inow, vmod));
_mm256_storeu_si256((__m256i*)&F[i + offset + p].val, mod_sub_simd(l, r_inow, vmod));
}
} else
#endif
{
for (int i = 0; i < p; i++) {
uint32_t l = F[i + offset].val;
uint32_t r = montgomery_mul(F[i + offset + p].val, inow);
F[i + offset].val = mod_add(l, r);
F[i + offset + p].val = mod_sub(l, r);
}
}
inow = montgomery_mul(inow, sum_ie[__builtin_ctz(~s)]);
}
}
}
FormalPowerSeries<M> convolution(FormalPowerSeries<M> F1, FormalPowerSeries<M> F2) const {
int n = F1.size();
int m = F2.size();
if (n == 0 || m == 0) return FormalPowerSeries<M>(0);
if (std::min(n, m) <= 60) {
if (n < m) {
std::swap(n, m);
std::swap(F1.sz, F2.sz);
std::swap(F1.f, F2.f);
}
FormalPowerSeries<M> ans(n + m - 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < m; j++) {
ans[i + j].val = (ans[i + j].val + ull(F1[i].val) * F2[j].val) % M;
}
}
return ans;
}
int reference_size = upperpow2(n + m - 1);
F1.resize(reference_size);
F2.resize(reference_size);
for (int i = 0; i < reference_size; i++) {
F1[i].val = to_montgomery(F1[i].val);
F2[i].val = to_montgomery(F2[i].val);
}
inplaceDFT(F1);
inplaceDFT(F2);
int i = 0;
#if defined(__AVX512F__) && defined(__AVX512DQ__)
__m512i vmod512 = _mm512_set1_epi32(M);
__m512i vneg_inv512 = _mm512_set1_epi32(NEG_INV);
for (; i + 15 < reference_size; i += 16) {
__m512i a = _mm512_loadu_si512((__m512i*)&F1[i].val);
__m512i b = _mm512_loadu_si512((__m512i*)&F2[i].val);
_mm512_storeu_si512((__m512i*)&F1[i].val, montgomery_mul_simd(a, b, vmod512, vneg_inv512));
}
#elif defined(__AVX2__)
__m256i vmod256 = _mm256_set1_epi32(M);
__m256i vneg_inv256 = _mm256_set1_epi32(NEG_INV);
for (; i + 7 < reference_size; i += 8) {
__m256i a = _mm256_loadu_si256((__m256i*)&F1[i].val);
__m256i b = _mm256_loadu_si256((__m256i*)&F2[i].val);
_mm256_storeu_si256((__m256i*)&F1[i].val, montgomery_mul_simd(a, b, vmod256, vneg_inv256));
}
#endif
for (; i < reference_size; i++) {
F1[i].val = montgomery_mul(F1[i].val, F2[i].val);
}
inplaceIDFT(F1);
uint iz_normal = modpow<M>(reference_size, M-2);
i = 0;
#if defined(__AVX512F__) && defined(__AVX512DQ__)
__m512i viz_normal512 = _mm512_set1_epi32(iz_normal);
for (; i + 15 < reference_size; i += 16) {
__m512i a = _mm512_loadu_si512((__m512i*)&F1[i].val);
_mm512_storeu_si512((__m512i*)&F1[i].val, montgomery_mul_simd(a, viz_normal512, vmod512, vneg_inv512));
}
#elif defined(__AVX2__)
__m256i viz_normal256 = _mm256_set1_epi32(iz_normal);
for (; i + 7 < reference_size; i += 8) {
__m256i a = _mm256_loadu_si256((__m256i*)&F1[i].val);
_mm256_storeu_si256((__m256i*)&F1[i].val, montgomery_mul_simd(a, viz_normal256, vmod256, vneg_inv256));
}
#endif
for (; i < reference_size; i++) {
F1[i].val = reduce(ull(F1[i].val) * iz_normal);
}
F1.resize(n + m - 1);
return F1;
}
// f(g(x)) mod x^len(f) を求める
FormalPowerSeries<M> composition(FormalPowerSeries<M> f, FormalPowerSeries<M> g, const mod_table<M>& mtable) const {
f.resize(max(f.size(),g.size()));
int N = f.size();
if (N == 0) return FormalPowerSeries<M>(0);
int n = 1;
while (n < (int)f.size()) n *= 2;
f.resize(n), g.resize(n);
FormalPowerSeries<M> W(2*n);
{
// bit reverse order
vector<int> btr(2*n);
int log = 31-__builtin_clz(2*n);
for (int i = 0; i < 2*n; i++){btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (log - 1));}
int t = log_max_length;
uint r = last_powroot;
ull dw = modpow<M>(inverse_mod<ll,ll>(r, M), (1<<t)/(4*n));
uint w = 1;
for (auto i: btr) { W[i].val = w, w = w*dw%M; }
}
auto rec = [&](auto &rec, int n, int k, FormalPowerSeries<M> &Q) -> FormalPowerSeries<M> {
if (n == 1) {
reverse(f.f.begin(), f.f.end());
FormalPowerSeries<M> p(2*k);
for (int i = 0; i < k; i++) p[2 * i].val = f[i].val;
return p;
}
Q.resize(4*n*k);
Q[2*n*k].val = 1;
inplaceDFT(Q);
FormalPowerSeries<M> nxt_Q(2*n*k);
for (int i = 0; i < 2*n*k; i++) nxt_Q[i].val = Q[2*i].val*(ull)Q[2*i+1].val%M;
inplaceIDFT(nxt_Q);
ull iz = inverse_mod<ll,ll>(nxt_Q.size(), M);
for (uint i = 0; i < nxt_Q.size(); i++) nxt_Q[i].val = nxt_Q[i].val*iz%M;
for (int j = 0; j < 2*k; j++) for (int i = n/2; i < n; i++) nxt_Q[n*j+i].val = 0;
nxt_Q[0].val = 0;
FormalPowerSeries p = rec(rec, n/2, 2*k, nxt_Q);
for (int j = 0; j < 2*k; j++) for (int i = n/2; i < n; i++) p[n * j + i].val = 0;
iz = mtable.invmodlist[p.size()];
for (uint i = 0; i < p.size(); i++) p[i].val = p[i].val*iz%M;
inplaceIDFT_T(p);
p.resize(4*n*k);
for (int i = 2*n*k-1; i >= 0; i--){
p[2*i+1].val = (ull)(mtable.invmodlist[2]-1+M)%M*W[i].val%M*Q[2*i].val%M*p[i].val%M;
p[2*i].val = (ull)mtable.invmodlist[2]*W[i].val%M*Q[2*i+1].val%M*p[i].val%M;
}
inplaceDFT_T(p);
p.resize(2*n*k);
return p;
};
FormalPowerSeries<M> Q(2*n);
for(int i = 0; i < n; i++) Q[i].val = g[i].val == 0 ? 0 : M-g[i].val;
FormalPowerSeries p = rec(rec, n, 1, Q);
p.resize(n);
reverse(p.f.begin(), p.f.end());
p.resize(N);
return p;
}
};
#endif /* FPS_OPERATION_HPP_ */
using mint = constant_modint<998244353>;
fps_operator<998244353> op;
mod_table<998244353> mt(1048576);
int main(){
ios::sync_with_stdio(false);
std::cin.tie(nullptr);
ll N,M;
cin >> N >> M;
FormalPowerSeries<998244353> F(N),G(N), H(N);
vin(F);
vin(G);
vin(H);
mint temp = F[1].pow(M);
for (int i = 0; i < N; i++){
G[i] *= temp;
}
auto ans = op.composition(H,G,mt);
vout(ans);
}