結果
| 問題 | No.3793 Multishift |
| コンテスト | |
| ユーザー |
 yuki4869.
|
| 提出日時 | 2026-11-07 02:09:30 |
| 言語 | C++23 (gcc 15.3.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 1,474 ms / 7,000 ms |
| + 335µs | |
| コード長 | 33,015 bytes |
| 記録 | |
| コンパイル時間 | 5,831 ms |
| コンパイル使用メモリ | 416,528 KB |
| 実行使用メモリ | 62,064 KB |
| 最終ジャッジ日時 | 2026-11-07 02:09:50 |
| 合計ジャッジ時間 | 19,075 ms |
|
ジャッジサーバーID （参考情報） |
judge4_0 / judge1_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 1 |
| other | AC * 13 |
ソースコード
// https://github.com/yuki-4869/Algorithm/blob/main/template.cpp
#include <bits/stdc++.h>
#include <atcoder/all>
using namespace std; using namespace atcoder;
using ll=long long; using LL=__int128; using mint=modint998244353;
using ld = long double;
using pll=pair<ll,ll>; using tu=tuple<ll,ll,ll>; using tu4 = tuple<ll,ll,ll,ll>;using vl=vector<ll>; using vb=vector<bool>; using vvb=vector<vb>; using vvl=vector<vl>; using vs=vector<string>; using vld=vector<long double>; using vc=vector<char>; using vmi=vector<mint>; using vp = vector<pll>; using vt = vector<tu>;using graph = vvl; using Graph = vector<vp>;
template<class T=ll> auto vec2(ll n,ll m,T x=T()){return vector(n,vector<T>(m,x));} template<class T=ll> auto vec3(ll n,ll m,ll l,T x=T()){return vector(n,vector(m,vector<T>(l,x)));}
template<class T=ll> auto vec4(ll a,ll b,ll c,ll d,T x=T()){return vector(a,vector(b,vector(c,vector<T>(d,x))));}
template<class T> using rpq=priority_queue<T,vector<T>,greater<T>>; // 小さい順priority_queue
using inverse_priority_queue=rpq<ll>;
const ll Inf=2147483647LL, Inf9=1000000000LL, inf=9223372036854775807LL, inf18=1000000000000000000LL, mod=998244353LL;
#define all(v) (v).begin(),(v).end()
#define rall(v) (v).rbegin(),(v).rend()
#define Re(n) for(ll i = 0; i<(ll)(n);i++)
#define rep(i,n) for(ll i=0;i<(ll)(n);i++)
#define rrep(i,n) for(ll i=(ll)(n)-1;i>=0;i--)
#define rep1(i,n) for(ll i=1;i<=(ll)(n);i++)
#define rep3(i,s,t) for(ll i=(ll)(s);i<(ll)(t);i++)
#define pb push_back
#define p(s) cout<<(s)<<'\n'
#define p2(s,t) cout<<(s)<<" "<<(t)<<'\n'
#define p3(s,t,u) cout<<(s)<<" "<<(t)<<" "<<(u)<<'\n'
#define p4(s,t,u,v) cout<<(s)<<" "<<(t)<<" "<<(u)<<" "<<(v)<<'\n'
#define pe(s) cout<<(s)<<endl
#define pld(...) ([](const auto&... values)->std::ostream&{auto flags=std::cout.flags();auto precision=std::cout.precision();std::cout<<std::fixed<<std::setprecision(12);size_t i=0;((std::cout<<(i++?" ":"")<<static_cast<long double>(values)),...);std::cout<<'\n';std::cout.flags(flags);std::cout.precision(precision);return std::cout;}(__VA_ARGS__))
#define pld2(s,t) pld((s),(t))
#define pld3(s,t,u) pld((s),(t),(u))
#define pld4(s,t,u,v) pld((s),(t),(u),(v))
#define pe2(s,t) cout<<(s)<<" "<<(t)<<endl
#define pe3(s,t,u) cout<<(s)<<" "<<(t)<<" "<<(u)<<endl
#define pe4(s,t,u,v) cout<<(s)<<" "<<(t)<<" "<<(u)<<" "<<(v)<<endl
#define pval(s) cout<<(s).val()<<'\n'
#define Sort(A) sort(all(A))
#define RSort(A) sort(rall(A))
#define p_yes() p("Yes")
#define p_no() p("No")
template<class T> inline ll sz(const T& x){return (ll)x.size();} // サイズ取得
ostream& operator<<(ostream& os,const mint& x){return os<<x.val();} // mint出力
istream& operator>>(istream& is,mint& x){ll v; is>>v; x=mint(v); return is;} // mint入力
template<class T> inline bool chmin(T& a,const T& b){if(a>b){a=b; return true;} return false;} // 最小値更新
template<class T> inline bool chmax(T& a,const T& b){if(a<b){a=b; return true;} return false;} // 最大値更新
template<class T,class U> ll LB(const vector<T>& v,const U& a){return lower_bound(all(v),a)-v.begin();} // lower_boundのindex
template<class T,class U> ll UB(const vector<T>& v,const U& a){return upper_bound(all(v),a)-v.begin();} // upper_boundのindex
template<class T> T vec_min(const vector<T>& v){do { if (!(!v.empty())) std::abort(); } while (false); return *min_element(all(v));} // vectorの最小値
template<class T> T vec_max(const vector<T>& v){do { if (!(!v.empty())) std::abort(); } while (false); return *max_element(all(v));} // vectorの最大値
template<class T> T vec_sum(const vector<T>& v){T res=T(0); for(const auto& x:v) res+=x; return res;} // vectorの総和
template<class T=ll> vector<T> cin_vl(ll n){vector<T> ret(n); rep(i,n) cin>>ret[i]; return ret;} // vector入力、型省略時ll
template<class T=ll> vector<vector<T>> cin_vvl(ll h,ll w){vector<vector<T>> ret(h,vector<T>(w)); rep(i,h) rep(j,w) cin>>ret[i][j]; return ret;} // h*w入力、型省略時ll
template<class T=ll> pair<vector<T>,vector<T>> cin_xy(ll n){vector<T> x(n),y(n); rep(i,n) cin>>x[i]>>y[i]; return {x,y};} // x,y入力、型省略時ll
template<class T=ll> tuple<vector<T>,vector<T>,vector<T>> cin_xyz(ll n){vector<T> x(n),y(n),z(n); rep(i,n) cin>>x[i]>>y[i]>>z[i]; return {x,y,z};} // x,y,z入力、型省略時ll
template<class T=ll> vector<vector<T>> Cin_vl(ll n,ll k){vector<vector<T>> ret(k,vector<T>(n)); rep(i,n) rep(j,k) cin>>ret[j][i]; return ret;} // n行k列入力を列ごと保持、型省略時ll
template<class T=ll> vector<pair<T,T>> cin_pairs(ll n){vector<pair<T,T>> ret(n); rep(i,n) cin>>ret[i].first>>ret[i].second; return ret;} // pair列入力、型省略時ll
template<class T=ll> vector<tuple<T,T,T>> cin_tuples(ll n){vector<tuple<T,T,T>> ret(n); rep(i,n){T a,b,c; cin>>a>>b>>c; ret[i]={a,b,c};} return ret;} // 3要素tuple列入力、型省略時ll
template<class T> void cout_vl(const vector<T>& a,char sep=' ',char end='\n'){rep(i,sz(a)){if(i) cout<<sep; cout<<a[i];} cout<<end;} // vector出力
template<class T> void cout_ver(const vector<T>& a){for(const auto& value:a) cout<<value<<'\n';}
template<class T> void Cout_vl(const vector<T>& a,const vector<T>& b){rep(i,sz(a)) cout<<a[i]<<" "<<b[i]<<'\n';} // 2vector縦出力
template<class T> void cout_vvl(const vector<vector<T>>& a){for(const auto& e:a) cout_vl(e);} // 2次元vector出力
template<class T=ll> void Swap(vector<vector<T>>& v,ll& h,ll& w){vector<vector<T>> ret(w,vector<T>(h));rep(i,h) rep(j,w) ret[j][i]=v[i][j];v.swap(ret);swap(h,w);} // h*w行列をw*h行列に変換 hとwを自動でswap
ll ceil_div(ll a,ll b){do { if (!(b!=0)) std::abort(); } while (false); if(b<0) a=-a,b=-b; if(a>=0) return (a+b-1)/b; return -((-a)/b);} // 切り上げ除算
ll floor_div(ll a,ll b){do { if (!(b!=0)) std::abort(); } while (false); if(b<0) a=-a,b=-b; if(a>=0) return a/b; return -((-a+b-1)/b);} // 切り下げ除算
ll isqrt(ll n){ll x=sqrtl(n); while((LL)(x+1)*(x+1)<=n) x++; while((LL)x*x>n) x--; return x;} // floor(sqrt(n))
template<class T> vector<T> sort_unique_vec(vector<T> v){sort(all(v)); v.erase(unique(all(v)),v.end()); return v;} // sortして重複削除した新vector
template<class T> vector<ll> compress(const vector<T>& v){auto xs=sort_unique_vec(v); vector<ll> ret; for(auto& x:v) ret.pb(LB(xs,x)); return ret;} // 座標圧縮後のindex列
template<class T> vector<T> prefix_sum(const vector<T>& v){vector<T> s(sz(v)+1,0); rep(i,sz(v)) s[i+1]=s[i]+v[i]; return s;} // 1次元累積和
template<class T> vector<vector<T>> prefix_sum2(const vector<vector<T>>& a){int h=sz(a),w=sz(a[0]); vector<vector<T>> s(h+1,vector<T>(w+1,0)); rep(i,h)rep(j,w)s[i+1][j+1]=s[i+1][j]+s[i][j+1]-s[i][j]+a[i][j]; return s;} // 2次元累積和
template<class T> T range_sum(const vector<T>& s,ll l,ll r){return s[r]-s[l];} // 累積和から[l,r)の和
template<class T> vector<pair<T,ll>> run_length(const vector<T>& v){vector<pair<T,ll>> ret; for(auto& x:v){if(ret.empty()||ret.back().first!=x) ret.pb({x,1}); else ret.back().second++;} return ret;} // ランレングス圧縮
vector<pair<char,ll>> run_length(const string& s){vector<pair<char,ll>> ret; for(char c:s){if(ret.empty()||ret.back().first!=c) ret.pb({c,1}); else ret.back().second++;} return ret;} // 文字列ランレングス圧縮
template<class T> ll argmin(const vector<T>& v){do { if (!(!v.empty())) std::abort(); } while (false); return min_element(all(v))-v.begin();} // 最小値のindex
template<class T> ll argmax(const vector<T>& v){do { if (!(!v.empty())) std::abort(); } while (false); return max_element(all(v))-v.begin();} // 最大値のindex
template<class T> bool in_vec(const vector<T>& v,const T& x){return binary_search(all(v),x);} // sort済vectorに存在するか
bool even(ll x){return x%2==0;}
bool odd(ll x){return x%2!=0;}
ll grid_id(ll x,ll y,ll w){return x*w+y;} // 座標からGridのID
pll grid_xy(ll id,ll w){return {id/w,id%w};} // GridのIDから座標
ll gcd(ll a,ll b){if(b==0) return abs(a); return gcd(b,a%b);} // 最大公約数
ll lcm(ll a,ll b){return a/gcd(a,b)*b;} // 最小公倍数
ll Exp(ll x,ll y){ll cnt=-1; while(x>0) x/=y,cnt++; return cnt;} // y進桁数-1
ll powmod(ll x,ll k,ll m){if(k==0) return 1%m; ll r=powmod(x,k/2,m); r=(LL)r*r%m; if(k&1) r=(LL)r*x%m; return r;} // mod累乗
template<class T> T square(T x){return x*x;} // 二乗
template<class T> void Uni_erase(vector<T>& v){sort(all(v)); v.erase(unique(all(v)),v.end());} // sortして重複削除
inline bool OutIn(ll x,ll n){return 0<=x&&x<n;}
inline bool OutIn(ll x,ll y,ll h,ll w){return 0<=x&&x<h&&0<=y&&y<w;}
inline bool OutIn(ll x,ll y,ll z,ll a,ll b,ll c){return 0<=x&&x<a&&0<=y&&y<b&&0<=z&&z<c;} // グリッド範囲内判定
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
ll Random(ll l,ll r){uniform_int_distribution<ll> dist(l,r-1); return dist(rng);} // [l,r)乱数
void YesNo(bool x){x?p_yes():p_no();} // Yes/No出力
vl vl_n(ll n){vl ret(n); rep(i,n) ret[i]=i; return ret;} // 0..n-1のvector生成
vl vl_N(ll n){vl ret(n); rep(i,n) ret[i]=i+1; return ret;} // 1..nのvector生成
vl dijkstra(const graph& g,const vl& starts){ vl dist(g.size(),inf); queue<ll> q; for(ll s:starts) if(dist[s]==inf){dist[s]=0; q.push(s);} while(!q.empty()){ ll v=q.front(); q.pop(); for(ll to:g[v]) if(dist[to]==inf){dist[to]=dist[v]+1; q.push(to);} } return dist; }
vl dijkstra(const graph& g,ll s){return dijkstra(g,vl{s});}
vl dijkstra(const graph& g,initializer_list<ll> starts){return dijkstra(g,vl(starts));}
vl dijkstra(const Graph& g,const vl& starts){ vl dist(g.size(),inf); rpq<pll> q; for(ll s:starts) if(dist[s]==inf){dist[s]=0; q.push({0,s});} while(!q.empty()){ auto [d,v]=q.top(); q.pop(); if(d!=dist[v]) continue; for(auto [to,w]:g[v]){ if(w>=inf-d) continue; ll nd=d+w; if(nd<dist[to]){dist[to]=nd; q.push({nd,to});} } } return dist; }
vl dijkstra(const Graph& g,ll s){return dijkstra(g,vl{s});}
vl dijkstra(const Graph& g,initializer_list<ll> starts){return dijkstra(g,vl(starts));}
const vl Pow2={1LL,2LL,4LL,8LL,16LL,32LL,64LL,128LL,256LL,512LL,1024LL,2048LL,4096LL,8192LL,16384LL,32768LL,65536LL,131072LL,262144LL,524288LL,1048576LL,2097152LL,4194304LL,8388608LL,16777216LL,33554432LL,67108864LL,134217728LL,268435456LL,536870912LL,1073741824LL,2147483648LL,4294967296LL,8589934592LL,17179869184LL,34359738368LL,68719476736LL,137438953472LL,274877906944LL,549755813888LL,1099511627776LL,2199023255552LL,4398046511104LL,8796093022208LL,17592186044416LL,35184372088832LL,70368744177664LL,140737488355328LL,281474976710656LL,562949953421312LL,1125899906842624LL,2251799813685248LL,4503599627370496LL,9007199254740992LL,18014398509481984LL,36028797018963968LL,72057594037927936LL,144115188075855872LL,288230376151711744LL,576460752303423488LL,1152921504606846976LL};
vl dx={0,1,0,-1}, dy={1,0,-1,0};
vl dx8={1,1,0,-1,-1,-1,0,1}, dy8={0,1,1,1,0,-1,-1,-1}; // 8近傍
vector<mint> fac,finv,pow2;
void factorial(ll N=1000000LL){ll M=N+100; fac.assign(M+1,1); finv.assign(M+1,1); pow2.assign(M+1,1); rep1(i,M) fac[i]=fac[i-1]*i; finv[M]=fac[M].inv(); for(ll i=M;i>=1;i--) finv[i-1]=finv[i]*i; rep1(i,M) pow2[i]=pow2[i-1]*2;} // N+100まで階乗・逆階乗・2の累乗前計算
void ensure_factorial(ll N){if((ll)fac.size()>N) return; factorial(N);} // 階乗表の不足分確保
mint nCr(ll N,ll r){if(N<0||r<0||N<r) return 0; ensure_factorial(N); return fac[N]*finv[r]*finv[N-r];} // 二項係数
mint nPr(ll N,ll r){if(N<0||r<0||N<r) return 0; ensure_factorial(N); return fac[N]*finv[N-r];} // 順列数
using mint=atcoder::modint998244353;
using vmi=vector<mint>;
namespace fps_detail {
using coefficient=mint;
using coefficient_vector=vector<coefficient>;
inline void require(bool ok){if(!ok)std::abort();}
inline const coefficient_vector& inverses(ll n){
require(0<=n&&n<coefficient::mod());
static coefficient_vector table={0,1};
for(ll i=table.size();i<=n;i++)
table.push_back(coefficient(coefficient::mod()-(ll)(coefficient::mod()/i)*table[coefficient::mod()%i].val()%coefficient::mod()));
return table;
}
struct Factorials{
coefficient_vector fac={1},ifac={1};
void ensure(ll n){
const auto& inverse=inverses(n);
for(ll i=fac.size();i<=n;i++){
fac.push_back(fac.back()*i);
ifac.push_back(ifac.back()*inverse[i]);
}
}
};
inline const Factorials& factorials(ll n){
static Factorials table;
table.ensure(n);
return table;
}
inline coefficient_vector naive(const coefficient_vector& a,const coefficient_vector& b){
coefficient_vector result(a.size()+b.size()-1);
if(a.size()<b.size()){
for(ll j=0;j<(ll)b.size();j++)if(b[j]!=coefficient(0))
for(ll i=0;i<(ll)a.size();i++)result[i+j]+=a[i]*b[j];
}else{
for(ll i=0;i<(ll)a.size();i++)if(a[i]!=coefficient(0))
for(ll j=0;j<(ll)b.size();j++)result[i+j]+=a[i]*b[j];
}
return result;
}
inline coefficient_vector convolve(const coefficient_vector& a,const coefficient_vector& b){
if(a.empty()||b.empty())return {};
require(a.size()<=INT_MAX&&b.size()<=INT_MAX&&a.size()+b.size()-1<=INT_MAX);
if(min(a.size(),b.size())<=60)return naive(a,b);
require(a.size()+b.size()-1<=(1U<<23));
return atcoder::convolution(a,b);
}
}
struct FPS:vmi{
using vmi::vmi;
FPS(const vmi& a):vmi(a){}
FPS(vmi&& a):vmi(move(a)){}
private:
using C=mint;
static void precision(ll n){fps_detail::require(0<=n&&n<=C::mod());}
static FPS multiply_prefix(const FPS& a,const FPS& b,ll n){
if(!n||a.empty()||b.empty())return {};
return (a.pre(n)*b.pre(n)).pre(n);
}
public:
FPS pre(ll n)const{
fps_detail::require(n>=0);
return FPS(begin(),begin()+min<ll>(n,size()));
}
FPS rev()const{
FPS result=*this;
reverse(result.begin(),result.end());
return result;
}
FPS& shrink(){
while(!empty()&&back()==C(0))pop_back();
return *this;
}
FPS diff()const{
if(size()<2)return {};
FPS result(size()-1);
for(ll i=1;i<(ll)size();i++)result[i-1]=(*this)[i]*i;
return result;
}
FPS integral()const{
fps_detail::require(size()<(size_t)C::mod());
const auto& inverse=fps_detail::inverses(size());
FPS result(size()+1);
for(ll i=0;i<(ll)size();i++)result[i+1]=(*this)[i]*inverse[i+1];
return result;
}
FPS& operator+=(const FPS& other){
if(size()<other.size())resize(other.size());
for(ll i=0;i<(ll)other.size();i++)(*this)[i]+=other[i];
return *this;
}
FPS& operator-=(const FPS& other){
if(size()<other.size())resize(other.size());
for(ll i=0;i<(ll)other.size();i++)(*this)[i]-=other[i];
return *this;
}
FPS& operator*=(const FPS& other){return *this=*this*other;}
FPS& operator+=(C value){
if(empty())resize(1);
(*this)[0]+=value;
return *this;
}
FPS& operator-=(C value){
if(empty())resize(1);
(*this)[0]-=value;
return *this;
}
FPS& operator*=(C value){
for(auto& x:*this)x*=value;
return *this;
}
FPS& operator/=(C value){
fps_detail::require(value!=C(0));
return *this*=value.inv();
}
FPS operator-()const{
FPS result=*this;
for(auto& x:result)x=-x;
return result;
}
friend FPS operator+(FPS a,const FPS& b){return a+=b;}
friend FPS operator-(FPS a,const FPS& b){return a-=b;}
friend FPS operator*(const FPS& a,const FPS& b){
return FPS(fps_detail::convolve(a,b));
}
friend FPS operator+(FPS a,C b){return a+=b;}
friend FPS operator+(C a,FPS b){return b+=a;}
friend FPS operator-(FPS a,C b){return a-=b;}
friend FPS operator-(C a,const FPS& b){return -b+a;}
friend FPS operator*(FPS a,C b){return a*=b;}
friend FPS operator*(C a,FPS b){return b*=a;}
friend FPS operator/(FPS a,C b){return a/=b;}
FPS operator<<(ll k)const{
fps_detail::require(k>=0&&(size_t)k<=max_size()-size());
FPS result((size_t)k+size());
copy(begin(),end(),result.begin()+k);
return result;
}
FPS operator>>(ll k)const{
fps_detail::require(k>=0);
if((size_t)k>=size())return {};
return FPS(begin()+k,end());
}
FPS& operator<<=(ll k){return *this=*this<<k;}
FPS& operator>>=(ll k){return *this=*this>>k;}
FPS inv(ll n)const{
precision(n);
if(!n)return {};
fps_detail::require(!empty()&&(*this)[0]!=C(0));
FPS f=pre(n);
f.shrink();
C inverse=f[0].inv();
if(f.size()==1){FPS result(n);result[0]=inverse;return result;}
ll initial=min(n,32LL);
FPS result(initial);
result[0]=inverse;
for(ll i=1;i<initial;i++){
C sum=0;
for(ll j=1;j<=i&&j<(ll)f.size();j++)sum+=f[j]*result[i-j];
result[i]=-sum*inverse;
}
for(ll m=initial;m<n;m<<=1){
ll k=min(n,m<<1);
ll length=m<<1;
fps_detail::require(length<=(1LL<<23));
vector<C> first(length),second(length);
copy(f.begin(),f.begin()+min<ll>(k,f.size()),first.begin());
copy(result.begin(),result.end(),second.begin());
atcoder::internal::butterfly(first);
atcoder::internal::butterfly(second);
for(ll i=0;i<length;i++)first[i]*=second[i];
atcoder::internal::butterfly_inv(first);
fill(first.begin(),first.begin()+m,C(0));
fill(first.begin()+k,first.end(),C(0));
atcoder::internal::butterfly(first);
for(ll i=0;i<length;i++)first[i]*=second[i];
atcoder::internal::butterfly_inv(first);
C inverse_length=C(length).inv();
C scale=-inverse_length*inverse_length;
result.resize(k);
for(ll i=m;i<k;i++)result[i]=first[i]*scale;
}
return result;
}
FPS log(ll n)const{
precision(n);
if(!n)return {};
fps_detail::require(!empty()&&(*this)[0]==C(1));
FPS f=pre(n);
f.shrink();
if(n==1||f.size()==1)return FPS(n);
FPS result=multiply_prefix(f.diff(),f.inv(n-1),n-1).integral();
result.resize(n);
return result;
}
FPS exp(ll n)const{
precision(n);
if(!n)return {};
fps_detail::require(empty()||(*this)[0]==C(0));
FPS f=pre(n);
f.shrink();
if(f.empty()){FPS result(n);result[0]=1;return result;}
FPS result={1};
for(ll m=1;m<n;m<<=1){
ll k=min(n,m<<1);
FPS correction=f.pre(k)-result.log(k);
correction.resize(k);
correction[0]+=1;
result=multiply_prefix(result,correction,k);
result.resize(k);
}
return result;
}
FPS pow(ll exponent,ll n)const{
precision(n);
fps_detail::require(exponent>=0);
if(!n)return {};
if(!exponent){FPS result(n);result[0]=1;return result;}
FPS f=pre(n);
f.shrink();
if(f.empty())return FPS(n);
if(exponent==1){f.resize(n);return f;}
ll order=0;
while(f[order]==C(0))order++;
if(order&&exponent>(n-1)/order)return FPS(n);
ll shift=order*exponent;
C leading=f[order];
FPS unit=FPS(f.begin()+order,f.end()).pre(n-shift)/leading;
FPS result=(unit.log(n-shift)*C(exponent)).exp(n-shift)*leading.pow(exponent);
result<<=shift;
result.resize(n);
return result;
}
static pair<bool,C> mod_sqrt(C value){
const ll modulus=C::mod();
if(value==C(0))return {true,0};
if(value.pow((modulus-1)/2)!=C(1))return {false,0};
if(modulus%4==3){
C root=value.pow((modulus+1)/4);
if(root.val()>(-root).val())root=-root;
return {true,root};
}
ll odd=modulus-1,power=0;
while(!(odd&1)){odd>>=1;power++;}
static const C initial=[odd,modulus]{
C nonresidue=2;
while(nonresidue.pow((modulus-1)/2)==C(1))nonresidue++;
return nonresidue.pow(odd);
}();
C base=initial,root=value.pow((odd+1)/2),rest=value.pow(odd);
ll rank=power;
while(rest!=C(1)){
ll index=1;
C squared=rest*rest;
while(index<rank&&squared!=C(1)){squared*=squared;index++;}
fps_detail::require(index<rank);
C correction=base.pow(1LL<<(rank-index-1));
root*=correction;
base=correction*correction;
rest*=base;
rank=index;
}
if(root.val()>(-root).val())root=-root;
return {true,root};
}
FPS sqrt(ll n)const{
precision(n);
if(!n)return {};
FPS f=pre(n);
f.shrink();
if(f.empty())return FPS(n);
ll order=0;
while(f[order]==C(0))order++;
if(order&1)return {};
auto [ok,leading]=mod_sqrt(f[order]);
if(!ok)return {};
ll shift=order/2,terms=n-order;
FPS unit=FPS(f.begin()+order,f.end())/f[order];
FPS result={1};
C inverse_two=C(2).inv();
for(ll m=1;m<terms;m<<=1){
ll k=min(terms,m<<1);
result=(result+multiply_prefix(unit,result.inv(k),k))*inverse_two;
result.resize(k);
}
result*=leading;
result<<=shift;
result.resize(n);
return result;
}
pair<FPS,FPS> divmod(FPS divisor)const{
FPS dividend=*this;
dividend.shrink();divisor.shrink();
fps_detail::require(!divisor.empty());
if(dividend.size()<divisor.size())return {{},move(dividend)};
if(divisor.size()==1){dividend/=divisor[0];return {move(dividend),{}};}
ll count=dividend.size()-divisor.size()+1;
FPS quotient;
if(min<ll>(count,divisor.size())<=32){
quotient.resize(count);
C inverse=divisor.back().inv();
for(ll i=count;i--;){
C value=dividend[i+divisor.size()-1]*inverse;
quotient[i]=value;
if(value!=C(0))for(ll j=0;j<(ll)divisor.size();j++)dividend[i+j]-=value*divisor[j];
}
dividend.resize(divisor.size()-1);
}else{
quotient=multiply_prefix(dividend.rev().pre(count),divisor.rev().inv(count),count);
quotient.resize(count);
reverse(quotient.begin(),quotient.end());
dividend=(dividend-quotient*divisor).pre(divisor.size()-1);
}
quotient.shrink();dividend.shrink();
return {move(quotient),move(dividend)};
}
FPS& operator/=(const FPS& divisor){return *this=divmod(divisor).first;}
FPS& operator%=(const FPS& divisor){return *this=divmod(divisor).second;}
friend FPS operator/(FPS a,const FPS& b){return a/=b;}
friend FPS operator%(FPS a,const FPS& b){return a%=b;}
FPS taylor_shift(C value)const{
if(empty()||value==C(0))return *this;
ll length=size();
FPS f=*this;
f.shrink();
if(f.size()<=1)return *this;
ll count=f.size();
const auto& table=fps_detail::factorials(count-1);
FPS a(count),b(count);
C power=1;
for(ll i=0;i<count;i++){
a[count-1-i]=f[i]*table.fac[i];
b[i]=power*table.ifac[i];power*=value;
}
FPS product=a*b,result(length);
for(ll i=0;i<count;i++)result[i]=product[count-1-i]*table.ifac[i];
return result;
}
C eval(C value)const{
C result=0;
for(ll i=size();i--;)result=result*value+(*this)[i];
return result;
}
};
/*
細かい注意点
法998244353専用。mint=atcoder::modint998244353、vmi=vector<mint>。
基本テンプレートのll、vector等を使用する。code/codetではACLが利用可能。
codecfでは <atcoder/convolution> のincludeを有効にする。ライブラリ本体にはincludeを置かない。
998244353版と1000000007版はともにFPSを定義するので、どちらか一方を選ぶ。
係数は昇順で、添字iはX^iの係数。空配列も零多項式を表す。
pre(n)は先頭min(n,size())項の切り出しで、0埋めしない。nとシフト量は非負。
inv/log/exp/pow/sqrtのnは返す項数で、mod X^nを計算する。0<=n<=mod。
n=0なら空を返す。n>0ではinvはf[0]!=0、logはf[0]=1、expはf[0]=0（空も可）。
powの指数kは非負ll。sqrtは解がなければ空、零級数ならn項の零列を返す。
sqrtは一つの解を返し、先頭非零係数はmod_sqrtで得る小さい方の平方根を使う。
divmodと多項式に対する /, % は多項式除算。除数は非零。商・余りの末尾0は除去する。
integralは入力項数<mod。Taylor shiftは実際の次数<mod。整数の逆元が必要な範囲に注意。
乗算はACL convolution、逆元の倍長化はACL内部NTTを使用する。
NTTを使う各中間畳み込みの必要長（両入力項数の和-1）は2^23以下。
この上限は返り値の項数の上限ではない。次数打切り前の中間積にも適用する。
短い方の入力が60項以下の乗算は直接計算するので、NTT長の制限を受けない。
整数逆元・階乗・逆階乗は共有キャッシュ。最大使用サイズまで合計O(N)で拡張し、O(N)メモリを保持。
無効な引数はNDEBUG時も標準出力・標準エラーを出さず停止する。
vectorを継承し、size/resize、[]、begin/end、push_back、clear等は標準の意味で使える。
使える関数一覧
返り値 呼び出し方 内容 / 計算量
構築 FPS(), FPS(n), FPS(n,c) 空・n項の0・n項のc / O(1), O(n), O(n)
構築 FPS{a,b,...}, FPS(l,r) 初期化リスト・反復子範囲 / O(項数)
構築 FPS(v), FPS(move(v)) vmiからコピー・移動 / O(項数), O(1)
FPS f.pre(n) 先頭最大n項 / O(min(n,N))
FPS f.rev() 係数反転 / O(N)
FPS& f.shrink() 末尾0を除去 / O(N)
FPS f.diff() 微分 / O(N)
FPS f.integral() 積分定数0の積分 / O(N)
FPS f.inv(n) 逆数 mod X^n / O(M(n))
FPS f.log(n) 対数 mod X^n / O(M(n))
FPS f.exp(n) 指数 mod X^n / O(M(n))
FPS f.pow(k,n) k乗 mod X^n / O(M(n)+log k)
FPS f.sqrt(n) 平方根 mod X^n / O(M(n)+log mod+v2(mod-1)^2)
pair<bool,mint> FPS::mod_sqrt(a) 定数の平方根の有無と値 / O(log mod+v2(mod-1)^2)
pair<FPS,FPS> f.divmod(g) 多項式の商・余り / O(M(A+B))
FPS f.taylor_shift(c) f(X+c)、項数を保持 / O(M(N))
mint f.eval(x) Horner法でf(x) / O(N)
FPS f+g, f-g, -f 加減算・符号反転 / O(max(A,B)), O(N)
FPS f*g 畳み込み / O(M(A+B))
FPS f/g, f%g 多項式の商・余り / O(M(A+B))
FPS f+c, c+f, f-c, c-f 定数項との加減算 / O(N)
FPS f*c, c*f, f/c スカラー乗除算 / O(N)
FPS f<<k, f>>k X^k倍・低次k項削除 / O(N+k), O(N)
FPS& f+=g, f-=g 加減算して自己更新 / O(max(A,B))
FPS& f*=g, f/=g, f%=g 積・多項式の商・余りで自己更新 / O(M(A+B))
FPS& f+=c, f-=c 定数項を自己更新 / O(1)
FPS& f*=c, f/=c スカラーで自己更新 / O(N)
FPS& f<<=k, f>>=k シフトして自己更新 / O(N+k), O(N)
N=f.size()、AとBは二つの多項式の項数。M(n)=O(n log n)は乗算の上界。
法を固定した体演算での計算量。スカラーの逆元計算にはO(log mod)が加わる。
v2(t)はtを割る2の指数。空入力・定数・小さい積や除算は直接処理する。
基本の関数
FPS f={1,2,3},g={1,1};
FPS product=f*g;
FPS derivative=f.diff(),integrated=derivative.integral();
FPS inverse=f.inv(10);
FPS logarithm=f.log(10);
FPS exponential=FPS{0,1,2}.exp(10);
FPS power=f.pow(1000000000000LL,10);
FPS root=FPS{1,2,1}.sqrt(10);
auto [quotient,remainder]=f.divmod(g);
FPS shifted=f.taylor_shift(5);
mint value=f.eval(3);
*/
using FPSMint=FPS::value_type;
using FPSVector=vector<FPSMint>;
struct FPSPointSet{
private:
ll n,sz;
FPSVector x;
vector<FPS> p;
mutable bool weights_ready=false;
mutable FPSVector weights;
void prepare_weights()const{
if(weights_ready)return;
FPSVector derivative=eval(p[1].diff()),prefix(n+1,1);
for(ll i=0;i<n;i++){
fps_detail::require(derivative[i]!=FPSMint(0));
prefix[i+1]=prefix[i]*derivative[i];
}
weights.resize(n);
FPSMint inverse=prefix[n].inv();
for(ll i=n;i--;){weights[i]=inverse*prefix[i];inverse*=derivative[i];}
weights_ready=true;
}
public:
explicit FPSPointSet(const FPSVector& points):n(points.size()),sz(1),x(points){
fps_detail::require(n<=(ll)INT_MAX/2);
while(sz<n)sz<<=1;
p.assign(sz<<1,FPS{1});
for(ll i=0;i<n;i++)p[sz+i]={-x[i],1};
for(ll i=sz-1;i;i--)p[i]=p[i<<1]*p[i<<1|1];
}
FPS product()const{return p[1];}
FPSVector eval(const FPS& f)const{
FPSVector answer(n);
if(!n)return answer;
auto visit=[&](auto&& self,ll node,ll left,ll right,const FPS& polynomial)->void{
if(left>=n||polynomial.empty())return;
ll last=min(right,n);
if(last-left<=32||polynomial.size()<=32){
for(ll i=left;i<last;i++)answer[i]=polynomial.eval(x[i]);
return;
}
ll middle=left+(right-left)/2;
self(self,node<<1,left,middle,polynomial%p[node<<1]);
self(self,node<<1|1,middle,right,polynomial%p[node<<1|1]);
};
if(f.size()<p[1].size())visit(visit,1,0,sz,f);
else visit(visit,1,0,sz,f%p[1]);
return answer;
}
FPS interpolate(const FPSVector& values)const{
fps_detail::require((ll)values.size()==n);
if(!n)return {};
prepare_weights();
vector<FPS> partial(sz<<1);
for(ll i=0;i<n;i++)if(values[i]!=FPSMint(0))partial[sz+i]={values[i]*weights[i]};
for(ll i=sz-1;i;i--)
partial[i]=partial[i<<1]*p[i<<1|1]+partial[i<<1|1]*p[i<<1];
partial[1].shrink();
return move(partial[1]);
}
};
FPSVector multipoint_eval(const FPS& f,const FPSVector& x){return FPSPointSet(x).eval(f);}
FPS interpolate(const FPSVector& x,const FPSVector& y){return FPSPointSet(x).interpolate(y);}
/*
細かい注意点
fps.cppまたはfps_mod1000000007.cppを先に貼る。FPSMint/FPSVectorは選択したFPSの型。
多点評価は重複点にも対応。補間は法上で全評価点が相異なり、x.size()==y.size()であること。
空点集合のproductは{1}、evalとinterpolateは空。補間結果の末尾0は除去する。
評価点は構築時に固定。積木は内部情報として保持し、外部から変更できない。
同じ点集合での補間は微分値の逆元を初回にまとめて求めてキャッシュする。
少数の点・低次数の部分問題はHorner法へ切り替える。各NTT長はFPSの制限内。
無効な引数はNDEBUG時も標準エラーを出さず停止する。
使える関数一覧
構造体 FPSPointSet T(x) 固定評価点の積木を構築 / O(M(N)log N)、メモリO(N log N)
FPS T.product() Π(X-x[i])のコピー / O(N)
FPSVector T.eval(f) 点ごとのf(x[i]) / O(M(D+N)log(D+N))
FPS T.interpolate(y) 次数<Nの補間 / O(M(N)log N)、初回は逆元表も構築
FPSVector multipoint_eval(f,x) 一度きりの多点評価 / 構築+eval
FPS interpolate(x,y) 一度きりの補間 / 構築+interpolate
Nは評価点数、Dは入力多項式の項数。Mは選択したFPSの乗算計算量。
基本の関数
FPSVector x={1,2,3,4};
FPS f={5,3,2};
FPSPointSet tree(x);
auto y=tree.eval(f);
FPS restored=tree.interpolate(y);
FPS roots=tree.product();
auto y2=multipoint_eval(f,x);
FPS restored2=interpolate(x,y);
*/
signed main(){
ios::sync_with_stdio(false); cin.tie(nullptr);
ll n,m;
cin >> n >> m;
auto a = cin_vl(n+1);
auto b = cin_vl(m+1);
FPS A(n+1), B(m+1);
rep(i,n+1) A[i] = a[i];
rep(i,m+1) B[i] = b[i];
vmi v(n+m+1);
rep(i,n+m+1) v[i] = i-m;
FPS u = multipoint_eval(A,v);
FPS C = u*B;
v.resize(n+1);
rep(i,n+1) v[i] = i;
FPS U(n+1);
rep(i,n+1) U[i] = C[m+i];
FPS ans = interpolate(v,U);
ans.resize(n+1);
cout_vl(ans);
}