// https://github.com/yuki-4869/Algorithm/blob/main/template.cpp #include #include using namespace std; using namespace atcoder; using ll=long long; using LL=__int128; using mint=modint998244353; using ld = long double; using pll=pair; using tu=tuple; using tu4 = tuple;using vl=vector; using vb=vector; using vvb=vector; using vvl=vector; using vs=vector; using vld=vector; using vc=vector; using vmi=vector; using vp = vector; using vt = vector;using graph = vvl; using Graph = vector; template auto vec2(ll n,ll m,T x=T()){return vector(n,vector(m,x));} template auto vec3(ll n,ll m,ll l,T x=T()){return vector(n,vector(m,vector(l,x)));} template auto vec4(ll a,ll b,ll c,ll d,T x=T()){return vector(a,vector(b,vector(c,vector(d,x))));} template using rpq=priority_queue,greater>; // 小さい順priority_queue using inverse_priority_queue=rpq; 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)<std::ostream&{auto flags=std::cout.flags();auto precision=std::cout.precision();std::cout<(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)< inline ll sz(const T& x){return (ll)x.size();} // サイズ取得 ostream& operator<<(ostream& os,const mint& x){return os<>(istream& is,mint& x){ll v; is>>v; x=mint(v); return is;} // mint入力 template inline bool chmin(T& a,const T& b){if(a>b){a=b; return true;} return false;} // 最小値更新 template inline bool chmax(T& a,const T& b){if(a ll LB(const vector& v,const U& a){return lower_bound(all(v),a)-v.begin();} // lower_boundのindex template ll UB(const vector& v,const U& a){return upper_bound(all(v),a)-v.begin();} // upper_boundのindex template T vec_min(const vector& v){do { if (!(!v.empty())) std::abort(); } while (false); return *min_element(all(v));} // vectorの最小値 template T vec_max(const vector& v){do { if (!(!v.empty())) std::abort(); } while (false); return *max_element(all(v));} // vectorの最大値 template T vec_sum(const vector& v){T res=T(0); for(const auto& x:v) res+=x; return res;} // vectorの総和 template vector cin_vl(ll n){vector ret(n); rep(i,n) cin>>ret[i]; return ret;} // vector入力、型省略時ll template vector> cin_vvl(ll h,ll w){vector> ret(h,vector(w)); rep(i,h) rep(j,w) cin>>ret[i][j]; return ret;} // h*w入力、型省略時ll template pair,vector> cin_xy(ll n){vector x(n),y(n); rep(i,n) cin>>x[i]>>y[i]; return {x,y};} // x,y入力、型省略時ll template tuple,vector,vector> cin_xyz(ll n){vector 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 vector> Cin_vl(ll n,ll k){vector> ret(k,vector(n)); rep(i,n) rep(j,k) cin>>ret[j][i]; return ret;} // n行k列入力を列ごと保持、型省略時ll template vector> cin_pairs(ll n){vector> ret(n); rep(i,n) cin>>ret[i].first>>ret[i].second; return ret;} // pair列入力、型省略時ll template vector> cin_tuples(ll n){vector> ret(n); rep(i,n){T a,b,c; cin>>a>>b>>c; ret[i]={a,b,c};} return ret;} // 3要素tuple列入力、型省略時ll template void cout_vl(const vector& a,char sep=' ',char end='\n'){rep(i,sz(a)){if(i) cout< void cout_ver(const vector& a){for(const auto& value:a) cout< void Cout_vl(const vector& a,const vector& b){rep(i,sz(a)) cout< void cout_vvl(const vector>& a){for(const auto& e:a) cout_vl(e);} // 2次元vector出力 template void Swap(vector>& v,ll& h,ll& w){vector> ret(w,vector(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 vector sort_unique_vec(vector v){sort(all(v)); v.erase(unique(all(v)),v.end()); return v;} // sortして重複削除した新vector template vector compress(const vector& v){auto xs=sort_unique_vec(v); vector ret; for(auto& x:v) ret.pb(LB(xs,x)); return ret;} // 座標圧縮後のindex列 template vector prefix_sum(const vector& v){vector s(sz(v)+1,0); rep(i,sz(v)) s[i+1]=s[i]+v[i]; return s;} // 1次元累積和 template vector> prefix_sum2(const vector>& a){int h=sz(a),w=sz(a[0]); vector> s(h+1,vector(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 T range_sum(const vector& s,ll l,ll r){return s[r]-s[l];} // 累積和から[l,r)の和 template vector> run_length(const vector& v){vector> ret; for(auto& x:v){if(ret.empty()||ret.back().first!=x) ret.pb({x,1}); else ret.back().second++;} return ret;} // ランレングス圧縮 vector> run_length(const string& s){vector> ret; for(char c:s){if(ret.empty()||ret.back().first!=c) ret.pb({c,1}); else ret.back().second++;} return ret;} // 文字列ランレングス圧縮 template ll argmin(const vector& v){do { if (!(!v.empty())) std::abort(); } while (false); return min_element(all(v))-v.begin();} // 最小値のindex template ll argmax(const vector& v){do { if (!(!v.empty())) std::abort(); } while (false); return max_element(all(v))-v.begin();} // 最大値のindex template bool in_vec(const vector& 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 T square(T x){return x*x;} // 二乗 template void Uni_erase(vector& v){sort(all(v)); v.erase(unique(all(v)),v.end());} // sortして重複削除 inline bool OutIn(ll x,ll n){return 0<=x&&x 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 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 starts){return dijkstra(g,vl(starts));} vl dijkstra(const Graph& g,const vl& starts){ vl dist(g.size(),inf); rpq 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 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 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; namespace fps_detail { using coefficient=mint; using coefficient_vector=vector; inline void require(bool ok){if(!ok)std::abort();} inline const coefficient_vector& inverses(ll n){ require(0<=n&&n=0); return FPS(begin(),begin()+min(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()=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<>=(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 first(length),second(length); copy(f.begin(),f.begin()+min(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=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 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(-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 divmod(FPS divisor)const{ FPS dividend=*this; dividend.shrink();divisor.shrink(); fps_detail::require(!divisor.empty()); if(dividend.size()(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。 基本テンプレートのll、vector等を使用する。code/codetではACLが利用可能。 codecfでは の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は入力項数 FPS::mod_sqrt(a) 定数の平方根の有無と値 / O(log mod+v2(mod-1)^2) pair 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 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; struct FPSPointSet{ private: ll n,sz; FPSVector x; vector 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;ivoid{ if(left>=n||polynomial.empty())return; ll last=min(right,n); if(last-left<=32||polynomial.size()<=32){ for(ll i=left;i partial(sz<<1); for(ll i=0;i> 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); }