#line 1 "lib/template.hpp" #ifdef TEMPLATE #else #define TEMPLATE # pragma GCC optimize("O3") using namespace std; #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using uint=unsigned; using ll=long long; using ull=unsigned long long; using ld=long double; using pii=pair; using pll=pair; using i128=__int128; using u128=unsigned __int128; templateusing vc=vector; templateusing vvc=vc>; templateusing vvvc=vvc>; templateusing smpq=priority_queue,greater>; templateusing bipq=priority_queue; #define rep(i,n) for(ll i=0;i<(ll)(n);i++) #define REP(i,j,n) for(ll i=(j);i<(ll)(n);i++) #define DREP(i,n,m) for(ll i=(n);i>=(m);i--) #define drep(i,n) for(ll i=((n)-1);i>=0;i--) #define rall(x) x.rbegin(),x.rend() #define mp make_pair #define pb push_back #define fi first #define se second #define is insert #define bg begin() #define ed end() #define all(x) x.begin(),x.end() void scan(int&a) { cin >> a; } void scan(ll&a) { cin >> a; } void scan(string&a) { cin >> a; } void scan(char&a) { cin >> a; } void scan(uint&a) { cin >> a; } void scan(ull&a) { cin >> a; } void scan(bool&a) { cin >> a; } void scan(ld&a){ cin>> a;} template void scan(vector&a) { for(auto&x:a) scan(x); } void read() {} template void read(Head&head, Tail&... tail) { scan(head); read(tail...); } #define INT(...) int __VA_ARGS__; read(__VA_ARGS__); #define LL(...) ll __VA_ARGS__; read(__VA_ARGS__); #define ULL(...) ull __VA_ARGS__; read(__VA_ARGS__); #define STR(...) string __VA_ARGS__; read(__VA_ARGS__); #define VC(type, name, ...) vector name(__VA_ARGS__); read(name); #define VVC(type, name, size, ...) vector> name(size, vector(__VA_ARGS__)); read(name); templatevoid print(T a) { cout << a; } template void print(vectora) { for(int i=0;i<(int)a.size();i++){if(i)cout<<" ";print(a[i]);}cout< void PRT(T a) { print(a); cout < void PRT(Head head, Tail ... tail) { print(head); cout << " "; PRT(tail...); return; } template bool chmin(T &x, F y){ if(x>y){ x=y; return true; } return false; } template bool chmax(T &x, F y){ if(x T floor(T a, T b) { return a / b - (a % b && (a ^ b) < 0); } template T ceil(T x, T y) { return floor(x + y - 1, y); } template T bmod(T x, T y) { return x - y * floor(x, y); } template pair divmod(T x, T y) { T q = floor(x, y); return {q, x - q * y}; } void YesNo(bool b){ cout<<(b?"Yes":"No")<stovi(const string&s,const string&S){ vcv(s.size()); rep(i,s.size()){ auto t=S.find(s[i]); assert(t!=string::npos); v[i]=t; } return v; } template T isqrt(T x){ T F=sqrtl(x); while((F+1)*(F+1)<=x)F++; while(F*F>x)F--; return F; } template vvctrans(const vvc&a){ assert(a.size()&&a[0].size()); vvcb(a[0].size(),vc(a.size())); rep(i,a.size())rep(j,a[0].size()){ b[j][i]=a[i][j]; } return b; } template vctrans(const vc&a){ assert(a.size()&&a[0].size()); vcb(a[0].size(),string(a.size(),0)); rep(i,a.size())rep(j,a[0].size()){ b[j][i]=a[i][j]; } return b; } template int popcount(T n){ return __builtin_popcountll(n); } template L sum(vc&a){ return accumulate(all(a),L(0)); } template vcsubset(T S){ vcans; for(T x=S;x>0;x=(x-1)&S)ans.pb(x); ans.pb(0); return ans; } template T max(vc&a){ return *max_element(all(a)); } template T min(vc&a){ return *min_element(all(a)); } #ifndef COMPRESSER_STRUCT #define COMPRESSER_STRUCT template struct Compresser{ vcx; Compresser(int n=0){x.reserve(n);} Compresser(const vc&xs){ x=xs; } void push(T p){built=false;x.pb(p);} bool built=false; void build(){ if(!chmax(built,1))return; sort(all(x)); x.erase(unique(all(x)),x.end()); } int find(T v){ build(); auto itr=lower_bound(all(x),v)-x.begin(); if(itr==x.size()||x[itr]!=v)return -1; return itr; } int find_next(T v){ build(); return lower_bound(all(x),v)-x.begin(); } int size(){ build(); return x.size(); } T operator[](int i)const{ assert(0<=i&&i vc presum(vc &a){ vc ret(a.size()+1); rep(i,a.size())ret[i+1]=ret[i]+a[i]; return ret; } template vc &operator+=(vc &a,F b){ for (auto&v:a)v += b; return a; } template vc &operator-=(vc&a,F b){ for (auto&v:a)v-=b; return a; } template vc &operator*=(vc&a,F b){ for (auto&v:a)v*=b; return a; } template constexpr T POW(T a,T b){ T res=1; while(b){ if(b&1)res*=a; a*=a; b/=2; } return res; } constexpr ll ten(ll a){ return POW(10,a); } templateconstexpr T inf=numeric_limits::max()/2-1; template int tbit(T x){ using U=make_unsigned_t; U y=(U)x; return y?(int)bit_width(y)-1:-1; } template int lbit(T x){ using U=make_unsigned_t; U y=(U)x; return y?(int)countr_zero(y):-1; } template int tbit(T x,int p){ using U=make_unsigned_t; constexpr int W=numeric_limits::digits; U y=(U)x; if(p<0)return -1; if(p>=W-1)return tbit(y); return tbit(y&((U(1)<<(p+1))-1)); } template int lbit(T x,int p){ using U=make_unsigned_t; constexpr int W=numeric_limits::digits; U y=(U)x; if(p<0)return lbit(y); if(p>=W)return -1; return lbit(y&(~U(0)<>(istream&is,i128&x){ string s;is>>s; x=0; int i=0,neg=0; if(s[0]=='-')neg=1,i=1; for(;i<(int)s.size();i++)x=x*10+s[i]-'0'; if(neg)x=-x; return is; } ostream& operator<<(ostream&os,i128 x){ if(x==0)return os<<0; if(x<0)os<<"-"; using u128=__uint128_t; u128 y=x<0?-(u128)x:(u128)x; string s; while(y)s.pb('0'+y%10),y/=10; reverse(all(s)); return os< ostream& operator<<(ostream&os,const pair&p){ return os<<"("< ostream& operator<<(ostream&os,const array&a){ os<<"["; rep(i,N){ if(i)os<<", "; os< ostream& operator<<(ostream&os,const vc&a){ os<<"["; rep(i,a.size()){ if(i)os<<", "; os< ostream& operator<<(ostream&os,const deque&a){ os<<"["; rep(i,a.size()){ if(i)os<<", "; os< ostream& operator<<(ostream&os,const set&s){ os<<"{"; bool f=0; for(auto&x:s){ if(f)os<<", "; f=1; os< ostream& operator<<(ostream&os,const multiset&s){ os<<"{"; bool f=0; for(auto&x:s){ if(f)os<<", "; f=1; os< ostream& operator<<(ostream&os,const unordered_set&s){ os<<"{"; bool f=0; for(auto&x:s){ if(f)os<<", "; f=1; os< ostream& operator<<(ostream&os,const map&m){ os<<"{"; bool f=0; for(auto&x:m){ if(f)os<<", "; f=1; os< ostream& operator<<(ostream&os,const unordered_map&m){ os<<"{"; bool f=0; for(auto&x:m){ if(f)os<<", "; f=1; os< ostream& operator<<(ostream&os,queueq){ vca; while(q.size())a.pb(q.front()),q.pop(); return os< ostream& operator<<(ostream&os,stacks){ vca; while(s.size())a.pb(s.top()),s.pop(); return os< ostream& operator<<(ostream&os,priority_queueq){ vca; while(q.size())a.pb(q.top()),q.pop(); return os< void debug_out(const T&x,const Ts&...xs){ cout<sync_with_stdio(0); #ifdef LOCAL cout<0); mod=mod_; m=(i128(1)<<64)/mod; } unsigned reduce(uint64_t x){ assert(mod>0); x-=(((i128)x*m)>>64)*mod; return x struct DynamicModInt{ using u32=uint32_t; using u64=uint64_t; u32 val; DynamicModInt():val(0){} DynamicModInt(ll x){ ll v=x%get_mod(); if(v<0)v+=get_mod(); val=v; } static DynamicModInt raw(int v){ assert(v>=0); DynamicModInt mi; mi.val=v; return mi; } DynamicModInt &operator+=(const DynamicModInt&m){ if((val+=m.val)>=get_mod())val-=get_mod(); return *this; } DynamicModInt &operator-=(const DynamicModInt&m){ if((val+=(get_mod()-m.val))>=get_mod())val-=get_mod(); return *this; } DynamicModInt &operator*=(const DynamicModInt&m){ val=rem(u64(val)*m.val); return *this; } DynamicModInt &operator/=(const DynamicModInt&m){ val=rem(u64(val)*m.inv().val); return *this; } DynamicModInt operator-() const{ return DynamicModInt(val?get_mod()-val:0); } DynamicModInt operator+() const { return *this; } friend DynamicModInt operator+(DynamicModInt lhs, const DynamicModInt& rhs){ return lhs+=rhs; } friend DynamicModInt operator-(DynamicModInt lhs, const DynamicModInt& rhs){ return lhs-=rhs; } friend DynamicModInt operator*(DynamicModInt lhs, const DynamicModInt& rhs){ return lhs*=rhs; } friend DynamicModInt operator/(DynamicModInt lhs,const DynamicModInt&rhs){ return lhs/=rhs; } bool operator==(const DynamicModInt&p) const{ return p.val==val; } bool operator!=(const DynamicModInt&p) const{ return p.val!=val; } DynamicModInt pow(int64_t n) const{ DynamicModInt res(1),mul(val); while(n){ if(n%2)res*=mul; mul*=mul; n/=2; } return res; } friend ostream&operator<<(ostream&os,const DynamicModInt&p){ os<>(istream&is,DynamicModInt&p){ int64_t x; is>>x; p=DynamicModInt(x); return is; } DynamicModInt inv()const{ int64_t a=val,b=get_mod(),u=1,v=0,t; #ifdef LOCAL assert(gcd(a,b)==1); #endif while(b>0){ t=a/b; swap(a-=t*b,b); swap(u-=t*v,v); } return DynamicModInt(u); } inline static u32 rem(u64 x){return BarrettReduction().reduce(x);} static inline int &get_mod(){ static int mod=0; return mod; } static void set_mod(int md){ assert(0 struct BinomHasGetMod:false_type{}; template struct BinomHasGetMod>:true_type{}; template struct Binom{ private: static vector&fact_table(){static vectorv={1};return v;} static vector&invfact_table(){static vectorv={1};return v;} static vector&invs_table(){static vectorv={0};return v;} static int&built_mod(){static int mod=-1;return mod;} public: static void build(int n){ auto&_fact=fact_table(); auto&_invfact=invfact_table(); auto&_invs=invs_table(); if constexpr(BinomHasGetMod::value){ auto mod=mint::get_mod(); if(built_mod()!=mod){ _fact={1}; _invfact={1}; _invs={0}; built_mod()=mod; } } if(n<(int)_fact.size())return; int old=_fact.size(); _fact.resize(n+1); _invfact.resize(n+1); _invs.resize(n+1); if constexpr(BinomHasGetMod::value){ auto mod=mint::get_mod(); for(int i=old;i<=n;i++){ _fact[i]=_fact[i-1]*i; if(i==1)_invs[i]=1; else _invs[i]=-_invs[mod%i]*(mod/i); _invfact[i]=_invfact[i-1]*_invs[i]; } }else{ for(int i=old;i<=n;i++){ _fact[i]=_fact[i-1]*i; _invs[i]=mint(1)/i; _invfact[i]=_invfact[i-1]*_invs[i]; } } } static mint fact(int i){ assert(i>=0); build(i); return fact_table()[i]; } static mint invfact(int i){ assert(i>=0); build(i); return invfact_table()[i]; } static mint inv(int i){ assert(i>0); build(i); return invs_table()[i]; } static mint C(int a,int b){//aCb if(a<0||b<0||a-b<0)return mint(0); build(a); auto&_fact=fact_table(); auto&_invfact=invfact_table(); return _fact[a]*_invfact[b]*_invfact[a-b]; } static mint iC(int a,int b){//1/aCb if(a<0||b<0||a-b<0)return mint(0); build(a); auto&_fact=fact_table(); auto&_invfact=invfact_table(); return _fact[b]*_fact[a-b]*_invfact[a]; } static mint P(int a,int b){ if(a T extgcd(T a, T b, T &x, T &y) { T d = a; if(b != 0) { d = extgcd(b, a % b, y, x); y -= (a / b) * x; } else { x = 1; y = 0; } return d; } template pair inv(T x,T m){ T a1,a2; T res=extgcd(x,m,a1,a2); T md=m/res; a1=(a1%md+md)%md; return {a1,md}; } template pair mod_solve(T a,T b,T m){//return x s.t. ax=b mod m a%=m,b%=m;if(a<0)a+=m;if(b<0)b+=m; T g=gcd(gcd(a,b),m); a/=g,b/=g,m/=g; if(gcd(a,m)>1)return {-1,-1}; return {(inv(a,m).first*b)%m,inv(a,m).second}; } //https://nyaannyaan.github.io/library/modulo/mod-sqrt.hpp.html int64_t mod_sqrt(const int64_t &a, const int64_t &p) { assert(0 <= a && a < p); if (a < 2) return a; using Mint = DynamicModInt<409075245>; Mint::set_mod(p); if (Mint(a).pow((p - 1) >> 1) != 1) return -1; Mint b = 1, one = 1; while (b.pow((p - 1) >> 1) == 1) b += one; int64_t m = p - 1, e = 0; while (m % 2 == 0) m >>= 1, e += 1; Mint x = Mint(a).pow((m - 1) >> 1); Mint y = Mint(a) * x * x; x *= a; Mint z = Mint(b).pow(m); while (y != 1) { int64_t j = 0; Mint t = y; while (t != one) { j += 1; t *= t; } z = z.pow(int64_t(1) << (e - j - 1)); x *= z; z *= z; y *= z; e = j; } return x.val; } #line 3 "lib/math/static-mod-int.hpp" template struct StaticModInt{ static_assert(0=0); StaticModInt mi; mi.val=v; return mi; } StaticModInt &operator+=(const StaticModInt&m){ if((val+=m.val)>=mod)val-=mod; return *this; } StaticModInt &operator-=(const StaticModInt&m){ if((val+=(mod-m.val))>=mod)val-=mod; return *this; } StaticModInt &operator*=(const StaticModInt&m){ val=u64(val)*m.val%mod; return *this; } StaticModInt &operator/=(const StaticModInt&m){ val=u64(val)*m.inv().val%mod; return *this; } StaticModInt operator-() const{ return StaticModInt(mod-val); } StaticModInt operator+() const { return *this; } friend StaticModInt operator+(StaticModInt lhs, const StaticModInt& rhs){ return lhs+=rhs; } friend StaticModInt operator-(StaticModInt lhs, const StaticModInt& rhs){ return lhs-=rhs; } friend StaticModInt operator*(StaticModInt lhs, const StaticModInt& rhs){ return lhs*=rhs; } friend StaticModInt operator/(StaticModInt lhs,const StaticModInt&rhs){ return lhs/=rhs; } bool operator==(const StaticModInt&p) const{ return p.val==val; } bool operator!=(const StaticModInt&p) const{ return p.val!=val; } StaticModInt pow(int64_t n) const{ StaticModInt res(1),mul(val); while(n){ if(n%2)res*=mul; mul*=mul; n/=2; } return res; } friend ostream&operator<<(ostream&os,const StaticModInt&p){ os<>(istream&is,StaticModInt&p){ int64_t x; is>>x; p=StaticModInt(x); return is; } StaticModInt inv()const{ int64_t a=val,b=mod,u=1,v=0,t; #ifdef LOCAL assert(gcd(a,b)==1); #endif while(b>0){ t=a/b; swap(a-=t*b,b); swap(u-=t*v,v); } return StaticModInt(u); } }; #line 3 "lib/math/conv.hpp" #ifndef ATCODER_CONVOLUTION_HPP #define ATCODER_CONVOLUTION_HPP 1 #ifndef ATCODER_MATH_HPP #define ATCODER_MATH_HPP 1 #include #include #include #include #ifndef ATCODER_INTERNAL_MATH_HPP #define ATCODER_INTERNAL_MATH_HPP 1 #include #ifdef _MSC_VER #include #endif namespace atcoder { namespace internal { // @param m `1 <= m` // @return x mod m constexpr long long safe_mod(long long x, long long m) { x %= m; if (x < 0) x += m; return x; } // Fast modular multiplication by BarrettReduction reduction // Reference: https://en.wikipedia.org/wiki/Barrett_reduction // NOTE: reconsider after Ice Lake struct BarrettReduction { unsigned int _m; unsigned long long im; // @param m `1 <= m` explicit BarrettReduction(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {} // @return m unsigned int umod() const { return _m; } // @param a `0 <= a < m` // @param b `0 <= b < m` // @return `a * b % m` unsigned int mul(unsigned int a, unsigned int b) const { // [1] m = 1 // a = b = im = 0, so okay // [2] m >= 2 // im = ceil(2^64 / m) // -> im * m = 2^64 + r (0 <= r < m) // let z = a*b = c*m + d (0 <= c, d < m) // a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im // c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2 // ((ab * im) >> 64) == c or c + 1 unsigned long long z = a; z *= b; #ifdef _MSC_VER unsigned long long x; _umul128(z, im, &x); #else unsigned long long x = (unsigned long long)(((unsigned __int128)(z)*im) >> 64); #endif unsigned long long y = x * _m; return (unsigned int)(z - y + (z < y ? _m : 0)); } }; // @param n `0 <= n` // @param m `1 <= m` // @return `(x ** n) % m` constexpr long long pow_mod_constexpr(long long x, long long n, int m) { if (m == 1) return 0; unsigned int _m = (unsigned int)(m); unsigned long long r = 1; unsigned long long y = safe_mod(x, m); while (n) { if (n & 1) r = (r * y) % _m; y = (y * y) % _m; n >>= 1; } return r; } // Reference: // M. Forisek and J. Jancina, // Fast Primality Testing for Integers That Fit into a Machine Word // @param n `0 <= n` constexpr bool is_prime_constexpr(int n) { if (n <= 1) return false; if (n == 2 || n == 7 || n == 61) return true; if (n % 2 == 0) return false; long long d = n - 1; while (d % 2 == 0) d /= 2; constexpr long long bases[3] = {2, 7, 61}; for (long long a : bases) { long long t = d; long long y = pow_mod_constexpr(a, t, n); while (t != n - 1 && y != 1 && y != n - 1) { y = y * y % n; t <<= 1; } if (y != n - 1 && t % 2 == 0) { return false; } } return true; } template constexpr bool is_prime = is_prime_constexpr(n); // @param b `1 <= b` // @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g constexpr std::pair inv_gcd(long long a, long long b) { a = safe_mod(a, b); if (a == 0) return {b, 0}; // Contracts: // [1] s - m0 * a = 0 (mod b) // [2] t - m1 * a = 0 (mod b) // [3] s * |m1| + t * |m0| <= b long long s = b, t = a; long long m0 = 0, m1 = 1; while (t) { long long u = s / t; s -= t * u; m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b // [3]: // (s - t * u) * |m1| + t * |m0 - m1 * u| // <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u) // = s * |m1| + t * |m0| <= b auto tmp = s; s = t; t = tmp; tmp = m0; m0 = m1; m1 = tmp; } // by [3]: |m0| <= b/g // by g != b: |m0| < b/g if (m0 < 0) m0 += b / s; return {s, m0}; } // Compile time primitive root // @param m must be prime // @return primitive root (and minimum in now) constexpr int primitive_root_constexpr(int m) { if (m == 2) return 1; if (m == 167772161) return 3; if (m == 469762049) return 3; if (m == 754974721) return 11; if (m == 998244353) return 3; int divs[20] = {}; divs[0] = 2; int cnt = 1; int x = (m - 1) / 2; while (x % 2 == 0) x /= 2; for (int i = 3; (long long)(i)*i <= x; i += 2) { if (x % i == 0) { divs[cnt++] = i; while (x % i == 0) { x /= i; } } } if (x > 1) { divs[cnt++] = x; } for (int g = 2;; g++) { bool ok = true; for (int i = 0; i < cnt; i++) { if (pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1) { ok = false; break; } } if (ok) return g; } } template constexpr int primitive_root = primitive_root_constexpr(m); // @param n `n < 2^32` // @param m `1 <= m < 2^32` // @return sum_{i=0}^{n-1} floor((ai + b) / m) (mod 2^64) unsigned long long floor_sum_unsigned(unsigned long long n, unsigned long long m, unsigned long long a, unsigned long long b) { unsigned long long ans = 0; while (true) { if (a >= m) { ans += n * (n - 1) / 2 * (a / m); a %= m; } if (b >= m) { ans += n * (b / m); b %= m; } unsigned long long y_max = a * n + b; if (y_max < m) break; // y_max < m * (n + 1) // floor(y_max / m) <= n n = (unsigned long long)(y_max / m); b = (unsigned long long)(y_max % m); std::swap(m, a); } return ans; } } // namespace internal } // namespace atcoder #endif // ATCODER_INTERNAL_MATH_HPP namespace atcoder { long long pow_mod(long long x, long long n, int m) { assert(0 <= n && 1 <= m); if (m == 1) return 0; internal::BarrettReduction bt((unsigned int)(m)); unsigned int r = 1, y = (unsigned int)(internal::safe_mod(x, m)); while (n) { if (n & 1) r = bt.mul(r, y); y = bt.mul(y, y); n >>= 1; } return r; } long long inv_mod(long long x, long long m) { assert(1 <= m); auto z = internal::inv_gcd(x, m); assert(z.first == 1); return z.second; } // (rem, mod) std::pair crt(const std::vector& r, const std::vector& m) { assert(r.size() == m.size()); int n = int(r.size()); // Contracts: 0 <= r0 < m0 long long r0 = 0, m0 = 1; for (int i = 0; i < n; i++) { assert(1 <= m[i]); long long r1 = internal::safe_mod(r[i], m[i]), m1 = m[i]; if (m0 < m1) { std::swap(r0, r1); std::swap(m0, m1); } if (m0 % m1 == 0) { if (r0 % m1 != r1) return {0, 0}; continue; } // assume: m0 > m1, lcm(m0, m1) >= 2 * max(m0, m1) // (r0, m0), (r1, m1) -> (r2, m2 = lcm(m0, m1)); // r2 % m0 = r0 // r2 % m1 = r1 // -> (r0 + x*m0) % m1 = r1 // -> x*u0*g = r1-r0 (mod u1*g) (u0*g = m0, u1*g = m1) // -> x = (r1 - r0) / g * inv(u0) (mod u1) // im = inv(u0) (mod u1) (0 <= im < u1) long long g, im; std::tie(g, im) = internal::inv_gcd(m0, m1); long long u1 = (m1 / g); // |r1 - r0| < (m0 + m1) <= lcm(m0, m1) if ((r1 - r0) % g) return {0, 0}; // u1 * u1 <= m1 * m1 / g / g <= m0 * m1 / g = lcm(m0, m1) long long x = (r1 - r0) / g % u1 * im % u1; // |r0| + |m0 * x| // < m0 + m0 * (u1 - 1) // = m0 + m0 * m1 / g - m0 // = lcm(m0, m1) r0 += x * m0; m0 *= u1; // -> lcm(m0, m1) if (r0 < 0) r0 += m0; } return {r0, m0}; } long long floor_sum(long long n, long long m, long long a, long long b) { assert(0 <= n && n < (1LL << 32)); assert(1 <= m && m < (1LL << 32)); unsigned long long ans = 0; if (a < 0) { unsigned long long a2 = internal::safe_mod(a, m); ans -= 1ULL * n * (n - 1) / 2 * ((a2 - a) / m); a = a2; } if (b < 0) { unsigned long long b2 = internal::safe_mod(b, m); ans -= 1ULL * n * ((b2 - b) / m); b = b2; } return ans + internal::floor_sum_unsigned(n, m, a, b); } } // namespace atcoder #endif // ATCODER_MATH_HPP #include #include #include #include #include #ifndef ATCODER_INTERNAL_BITOP_HPP #define ATCODER_INTERNAL_BITOP_HPP 1 #ifdef _MSC_VER #include #endif #if __cplusplus >= 202002L #include #endif namespace atcoder { namespace internal { #if __cplusplus >= 202002L using std::bit_ceil; #else // @return same with std::bit::bit_ceil unsigned int bit_ceil(unsigned int n) { unsigned int x = 1; while (x < (unsigned int)(n)) x *= 2; return x; } #endif // @param n `1 <= n` // @return same with std::bit::countr_zero int countr_zero(unsigned int n) { #ifdef _MSC_VER unsigned long index; _BitScanForward(&index, n); return index; #else return __builtin_ctz(n); #endif } // @param n `1 <= n` // @return same with std::bit::countr_zero constexpr int countr_zero_constexpr(unsigned int n) { int x = 0; while (!(n & (1 << x))) x++; return x; } } // namespace internal } // namespace atcoder #endif // ATCODER_INTERNAL_BITOP_HPP namespace atcoder { namespace internal { template > struct FastFourierTransformInfo { static constexpr int rank2 = countr_zero_constexpr(mint::get_mod() - 1); std::array root; // root[i]^(2^i) == 1 std::array iroot; // root[i] * iroot[i] == 1 std::array rate2; std::array irate2; std::array rate3; std::array irate3; FastFourierTransformInfo() { root[rank2] = mint(g).pow((mint::get_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 - 2; i++) { rate2[i] = root[i + 2] * prod; irate2[i] = iroot[i + 2] * iprod; prod *= iroot[i + 2]; iprod *= root[i + 2]; } } { 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]; } } } }; template void butterfly(std::vector& a) { int n = int(a.size()); int h = internal::countr_zero((unsigned int)n); static const FastFourierTransformInfo info; int len = 0; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len < h) { if (h - len == 1) { int p = 1 << (h - len - 1); mint rot = 1; for (int s = 0; s < (1 << len); s++) { int offset = s << (h - len); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p] * rot; a[i + offset] = l + r; a[i + offset + p] = l - r; } if (s + 1 != (1 << len)) rot *= info.rate2[countr_zero(~(unsigned int)(s))]; } len++; } else { // 4-base int p = 1 << (h - len - 2); mint rot = 1, imag = info.root[2]; for (int s = 0; s < (1 << len); s++) { mint rot2 = rot * rot; mint rot3 = rot2 * rot; int offset = s << (h - len); for (int i = 0; i < p; i++) { auto mod2 = 1ULL * mint::get_mod() * mint::get_mod(); auto a0 = 1ULL * a[i + offset].val; auto a1 = 1ULL * a[i + offset + p].val * rot.val; auto a2 = 1ULL * a[i + offset + 2 * p].val * rot2.val; auto a3 = 1ULL * a[i + offset + 3 * p].val * rot3.val; auto a1na3imag = 1ULL * mint(a1 + mod2 - a3).val * imag.val; auto na2 = mod2 - a2; a[i + offset] = a0 + a2 + a1 + a3; a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3)); a[i + offset + 2 * p] = a0 + na2 + a1na3imag; a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag); } if (s + 1 != (1 << len)) rot *= info.rate3[countr_zero(~(unsigned int)(s))]; } len += 2; } } } template void butterfly_inv(std::vector& a) { int n = int(a.size()); int h = internal::countr_zero((unsigned int)n); static const FastFourierTransformInfo info; int len = h; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len) { if (len == 1) { int p = 1 << (h - len); mint irot = 1; for (int s = 0; s < (1 << (len - 1)); s++) { int offset = s << (h - len + 1); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p]; a[i + offset] = l + r; a[i + offset + p] = (unsigned long long)((unsigned int)(l.val - r.val) + mint::get_mod()) * irot.val; ; } if (s + 1 != (1 << (len - 1))) irot *= info.irate2[countr_zero(~(unsigned int)(s))]; } len--; } else { // 4-base int p = 1 << (h - len); mint irot = 1, iimag = info.iroot[2]; for (int s = 0; s < (1 << (len - 2)); s++) { mint irot2 = irot * irot; mint irot3 = irot2 * irot; int offset = s << (h - len + 2); for (int i = 0; i < p; i++) { auto a0 = 1ULL * a[i + offset + 0 * p].val; auto a1 = 1ULL * a[i + offset + 1 * p].val; auto a2 = 1ULL * a[i + offset + 2 * p].val; auto a3 = 1ULL * a[i + offset + 3 * p].val; auto a2na3iimag = 1ULL * mint((mint::get_mod() + a2 - a3) * iimag.val).val; a[i + offset] = a0 + a1 + a2 + a3; a[i + offset + 1 * p] = (a0 + (mint::get_mod() - a1) + a2na3iimag) * irot.val; a[i + offset + 2 * p] = (a0 + a1 + (mint::get_mod() - a2) + (mint::get_mod() - a3)) * irot2.val; a[i + offset + 3 * p] = (a0 + (mint::get_mod() - a1) + (mint::get_mod() - a2na3iimag)) * irot3.val; } if (s + 1 != (1 << (len - 2))) irot *= info.irate3[countr_zero(~(unsigned int)(s))]; } len -= 2; } } } template std::vector convolution_naive(const std::vector& a, const std::vector& b) { int n = int(a.size()), m = int(b.size()); std::vector ans(n + m - 1); if (n < m) { for (int j = 0; j < m; j++) { for (int i = 0; i < n; i++) { ans[i + j] += a[i] * b[j]; } } } else { for (int i = 0; i < n; i++) { for (int j = 0; j < m; j++) { ans[i + j] += a[i] * b[j]; } } } return ans; } template std::vector convolution_fft(std::vector a, std::vector b) { int n = int(a.size()), m = int(b.size()); int z = (int)internal::bit_ceil((unsigned int)(n + m - 1)); a.resize(z); internal::butterfly(a); b.resize(z); internal::butterfly(b); for (int i = 0; i < z; i++) { a[i] *= b[i]; } internal::butterfly_inv(a); a.resize(n + m - 1); mint iz = mint(z).inv(); for (int i = 0; i < n + m - 1; i++) a[i] *= iz; return a; } } // namespace internal template std::vector convolution(std::vector&& a, std::vector&& b) { int n = int(a.size()), m = int(b.size()); if (!n || !m) return {}; int z = (int)internal::bit_ceil((unsigned int)(n + m - 1)); assert((mint::get_mod() - 1) % z == 0); if (std::min(n, m) <= 60) return internal::convolution_naive(std::move(a), std::move(b)); return internal::convolution_fft(std::move(a), std::move(b)); } template std::vector convolution(const std::vector& a, const std::vector& b) { int n = int(a.size()), m = int(b.size()); if (!n || !m) return {}; int z = (int)internal::bit_ceil((unsigned int)(n + m - 1)); assert((mint::get_mod() - 1) % z == 0); if (std::min(n, m) <= 60) return internal::convolution_naive(a, b); return internal::convolution_fft(a, b); } template std::vector convolution(const std::vector& a, const std::vector& b) { int n = int(a.size()), m = int(b.size()); if (!n || !m) return {}; using mint = StaticModInt; int z = (int)internal::bit_ceil((unsigned int)(n + m - 1)); assert((mint::get_mod() - 1) % z == 0); std::vector a2(n), b2(m); for (int i = 0; i < n; i++) { a2[i] = mint(a[i]); } for (int i = 0; i < m; i++) { b2[i] = mint(b[i]); } auto c2 = convolution(std::move(a2), std::move(b2)); std::vector c(n + m - 1); for (int i = 0; i < n + m - 1; i++) { c[i] = c2[i].val; } return c; } std::vector convolution_ll(const std::vector& a, const std::vector& b) { int n = int(a.size()), m = int(b.size()); if (!n || !m) return {}; static constexpr unsigned long long MOD1 = 754974721; // 2^24 static constexpr unsigned long long MOD2 = 167772161; // 2^25 static constexpr unsigned long long MOD3 = 469762049; // 2^26 static constexpr unsigned long long M2M3 = MOD2 * MOD3; static constexpr unsigned long long M1M3 = MOD1 * MOD3; static constexpr unsigned long long M1M2 = MOD1 * MOD2; static constexpr unsigned long long M1M2M3 = MOD1 * MOD2 * MOD3; static constexpr unsigned long long i1 = internal::inv_gcd(MOD2 * MOD3, MOD1).second; static constexpr unsigned long long i2 = internal::inv_gcd(MOD1 * MOD3, MOD2).second; static constexpr unsigned long long i3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second; static constexpr int MAX_AB_BIT = 24; static_assert(MOD1 % (1ull << MAX_AB_BIT) == 1, "MOD1 isn't enough to support an array length of 2^24."); static_assert(MOD2 % (1ull << MAX_AB_BIT) == 1, "MOD2 isn't enough to support an array length of 2^24."); static_assert(MOD3 % (1ull << MAX_AB_BIT) == 1, "MOD3 isn't enough to support an array length of 2^24."); assert(n + m - 1 <= (1 << MAX_AB_BIT)); auto c1 = convolution(a, b); auto c2 = convolution(a, b); auto c3 = convolution(a, b); std::vector c(n + m - 1); for (int i = 0; i < n + m - 1; i++) { unsigned long long x = 0; x += (c1[i] * i1) % MOD1 * M2M3; x += (c2[i] * i2) % MOD2 * M1M3; x += (c3[i] * i3) % MOD3 * M1M2; // B = 2^63, -B <= x, r(real value) < B // (x, x - M, x - 2M, or x - 3M) = r (mod 2B) // r = c1[i] (mod MOD1) // focus on MOD1 // r = x, x - M', x - 2M', x - 3M' (M' = M % 2^64) (mod 2B) // r = x, // x - M' + (0 or 2B), // x - 2M' + (0, 2B or 4B), // x - 3M' + (0, 2B, 4B or 6B) (without mod!) // (r - x) = 0, (0) // - M' + (0 or 2B), (1) // -2M' + (0 or 2B or 4B), (2) // -3M' + (0 or 2B or 4B or 6B) (3) (mod MOD1) // we checked that // ((1) mod MOD1) mod 5 = 2 // ((2) mod MOD1) mod 5 = 3 // ((3) mod MOD1) mod 5 = 4 long long diff = c1[i] - internal::safe_mod((long long)(x), (long long)(MOD1)); if (diff < 0) diff += MOD1; static constexpr unsigned long long offset[5] = { 0, 0, M1M2M3, 2 * M1M2M3, 3 * M1M2M3}; x -= offset[diff % 5]; c[i] = x; } return c; } } // namespace atcoder #endif // ATCODER_CONVOLUTION_HPP #line 3 "lib/math/fps.hpp" template struct FormalPowerSeries:vc{ #define BF atcoder::internal::butterfly #define IBF atcoder::internal::butterfly_inv using vc::vc; using vc::operator[]; using poly=FormalPowerSeries; FormalPowerSeries(initializer_list init):vc(init){} FormalPowerSeries(const vc& v):vc(v){} FormalPowerSeries(vc&& v):vc(std::move(v)){} poly& operator=(initializer_list init){vc::operator=(init);return *this;} poly& operator=(const vc& v){vc::operator=(v);return *this;} poly& operator=(vc&& v){vc::operator=(std::move(v));return *this;} poly& operator+=(const poly& g){ int n=g.size(); if(n>(int)this->size())this->resize(n); for(int i=0;i(int)this->size())this->resize(n); for(int i=0;iempty()||g.empty()){this->clear();return *this;} auto res=atcoder::convolution(*this,g); *this=poly(res.begin(),res.end()); return *this; } poly& operator+=(const mint& v){if(this->empty())this->resize(1);(*this)[0]+=v;return *this;} poly& operator-=(const mint& v){if(this->empty())this->resize(1);(*this)[0]-=v;return *this;} poly& operator*=(const mint& v){for(auto& x:*this)x*=v;return *this;} poly& operator/=(const mint& v){mint inv=v.inv();for(auto& x:*this)x*=inv;return *this;} poly& operator>>=(const ll&v){*this=poly(this->begin()+min(this->size(),v),this->end());return *this;} poly& operator<<=(const ll&v){poly np(v+this->size());rep(i,this->size())np[i+v]=(*this)[i];*this=np;return *this;} poly operator+(const poly& g)const{return poly(*this)+=g;} poly operator-(const poly& g)const{return poly(*this)-=g;} poly operator*(const poly& g)const{return poly(*this)*=g;} poly operator+(const mint& v)const{return poly(*this)+=v;} friend poly operator+(const mint& v,const poly& g){return poly(g)+=v;} poly operator-(const mint& v)const{return poly(*this)-=v;} friend poly operator-(const mint& v,const poly& g){return -g+=v;} poly operator*(const mint& v)const{return poly(*this)*=v;} friend poly operator*(const mint& v,const poly& g){return poly(g)*=v;} poly operator/(const mint& v)const{return poly(*this)/=v;} friend poly operator/(const mint& v,const poly& g){return poly(g)/=v;} poly operator>>(const ll&v)const{return poly(*this)>>=v;} poly operator<<(const ll&v)const{return poly(*this)<<=v;} poly operator-()const{ poly res=*this; for(auto& x:res)x=-x; return res; } poly pre(int sz)const{ return poly(this->begin(),this->begin()+min((int)this->size(),sz)); } //[0,n) poly inv(int deg=-1)const{ assert((*this).size()&&(*this)[0].val!=0); if(deg==-1)deg=this->size(); poly res{(*this)[0].inv()}; for(int len=1;leninv_table{0,1}; while((int)inv_table.size()<=this->size()){ int i=inv_table.size(); inv_table.push_back(-inv_table[mint::get_mod()%i]*mint(mint::get_mod()/i)); } poly res(this->size()+1); rep(i,this->size())res[i+1]=(*this)[i]*inv_table[i+1]; return res; } poly log(int deg)const{ if(deg==-1)deg=this->size(); assert(this->size()&&(*this)[0]==1); return (this->diff().pre(deg)*this->inv(deg)).integral().pre(deg); } poly exp(int deg=-1)const{ assert((*this).size()&&(*this)[0]==0); if(deg==-1)deg=this->size(); static vcinv_table{0,1}; while((int)inv_table.size()<=max((int)this->size(),deg)*2){ int i=inv_table.size(); inv_table.push_back(-inv_table[mint::get_mod()%i]*mint(mint::get_mod()/i)); } poly res{1},h{1}; for(int len=1;lensize())f_prime[i]=(*this)[i+1]*(i+1); BF(f_prime); poly w_prime(nxt); rep(i,nxt)w_prime[i]=g_prime[i]*H[i]+f_prime[i]*V[i]; IBF(w_prime);rep(i,nxt)w_prime[i]*=inv; for(int i=nxt-1;i>=len;--i){ mint fi=(i<(int)this->size())?(*this)[i]:0; w_prime[i]=fi-w_prime[i-1]*inv_table[i]; } for(int i=0;isize(); if(m==0){poly res(deg);res[0]=1;return res;} rep(i,this->size()){ if(m*i>deg)return poly(deg); if((*this)[i].val){ auto target=(*this)>>i; mint coef=mint((*this)[i]); mint inv=1/coef; rep(i,target.size())target[i]*=inv; return ((m*target.log(deg)).exp(deg)*mint(coef).pow(m))<<(i*m); } } return poly(deg); } void extend(int n){ if(this->size()resize(n); } poly shrink(){ auto res=(*this); while(res.size()&&res.back()==0)res.pop_back(); return res; } pairdiv(poly f)const{ if(this->size()size()-f.size()+1; auto rthis=*this;reverse(all(rthis)); auto rf=f;reverse(all(rf)); auto rq=(rthis.pre(qsize)*rf.inv(qsize)).pre(qsize); reverse(all(rq)); return {rq,(*this-rq*f).shrink()}; } poly div_only(poly f)const{ if(this->size()size()-f.size()+1; auto rthis=*this;reverse(all(rthis)); auto rf=f;reverse(all(rf)); auto rq=(rthis.pre(qsize)*rf.inv(qsize)).pre(qsize); reverse(all(rq)); return rq; } optional sqrt(int deg=-1)const{ if(deg==-1)deg=this->size(); if((this)->size()==0)return poly(deg); if((*this)[0]==0){ REP(i,1,this->size()){ if((*this)[i].val){ if(i%2)return nullopt; if(deg-i/2<=0)break; auto res=((*this)>>i).sqrt(deg-i/2); if(!res.has_value())return nullopt; auto ret=res.value(); ret<<=(i/2); ret.resize(deg); return ret; } } return poly(deg); } ll s=mod_sqrt((*this)[0].val,mint::get_mod()); if(s==-1)return nullopt; poly res{s}; mint inv=1/mint(2); for(int i=1;i f(vx) poly substitute(mint v){ mint coef=1; poly nxt(*this);rep(i,nxt.size()){ nxt[i]*=coef; coef*=v; } return nxt; } }; template FormalPowerSeries(vc)->FormalPowerSeries; template using fps=FormalPowerSeries; template FormalPowerSeriesto_fps(vc v){ FormalPowerSeriesres(v.size());rep(i,v.size())res[i]=v[i]; return res; } template FormalPowerSeriesto_fps(FormalPowerSeries v){ return v; } template vcto_vec(FormalPowerSeries v){ vcres(v.size());rep(i,v.size())res[i]=v[i]; return res; } template vcto_vec(vc v){ return v; } template FormalPowerSeries all_prod(vc> v){ assert(mint::get_mod()); using poly=FormalPowerSeries; if(v.empty())return {1}; auto cmp=[&](auto&a,auto&b){return a.size()>b.size();}; priority_queue,decltype(cmp)>pq(cmp); for(auto&p:v)pq.push(std::move(p)); while(pq.size()>=2){ auto p1=std::move(pq.top());pq.pop(); auto p2=std::move(pq.top());pq.pop(); pq.push(std::move(p1*p2)); } return pq.top(); } template FormalPowerSeries all_prod(vc> v){ vc> ps; ps.reserve(v.size()); for(auto&p:v)ps.push_back(std::move(p)); return all_prod(ps); } template struct SubproductTree{ using poly=FormalPowerSeries; int M; vc queries; vc T; bool built; SubproductTree(int n=0):built(false){ assert(mint::get_mod()); queries.reserve(n); } void add_query(mint x){ queries.push_back(x); } void build_dfs(int node,int l,int r){ if(r-l==1){ T[node]=poly{mint(1),-queries[l]}; return; } int m=(l+r)/2; build_dfs(2*node,l,m); build_dfs(2*node+1,m,r); T[node]=T[2*node]*T[2*node+1]; } void build_tree(){ M=queries.size(); if(M==0)return; T.assign(4*M,poly{}); build_dfs(1,0,M); built=true; } vc middle_product(const vc& A,const vc& B){ int n=A.size(); int m=B.size(); if(m==0)return vc(n+1,mint(0)); int S=1;while(S a=A;a.resize(S,mint(0)); vc b(S,mint(0)); for(int i=0;i res(n-m+1); for(int i=0;i& U,vc& res){ if(r-l==1){ res[l]=U[0]; return; } int m=(l+r)/2; const poly& T_left=T[2*node]; const poly& T_right=T[2*node+1]; int n=U.size(); int m_l=T_left.size(); int m_r=T_right.size(); int S=1;while(S u_ntt=U;u_ntt.resize(S,mint(0)); atcoder::internal::butterfly(u_ntt); vc t_r_ntt(S,mint(0)); for(int i=0;i u_l(S); for(int i=0;i t_l_ntt(S,mint(0)); for(int i=0;i u_r(S); for(int i=0;i Ul(n-m_r+1); for(int i=0;i Ur(n-m_l+1); for(int i=0;i calc(const poly& target){ if(!built)build_tree(); if(M==0)return {}; int N=target.size(); if(N==0)return vc(M,mint(0)); poly V=T[1].inv(N); vc F(target.begin(),target.end()); F.resize(N+M-1,mint(0)); vc U=middle_product(F,V); vc res(M); evaluate(1,0,M,U,res); return res; } }; template vc multipoint_evaluation(std::type_identity_t> f,vc query){ assert(mint::get_mod()); SubproductTree st(query.size()); for(auto& x:query)st.add_query(x); auto ans=st.calc(f); return ans; } //f(a),f(ar),...,f(ar^{m-1}) template vcchirp_z(std::type_identity_t>f,mint a,mint r,int m){ assert(mint::get_mod()); if(m==0)return {}; if(r==0){ vcans(m); mint coef=1; rep(i,f.size()){ ans[0]+=f[i]*coef; coef*=a; } REP(i,1,m)ans[i]=f[0]; return ans; } FormalPowerSeries A(f.size()),B(f.size()+m); vcir(f.size()+m+1),rp(f.size()+m+1); mint tmp=1;mint IR=1; mint inv_r=1/r; rep(i,ir.size()){ ir[i]=tmp; tmp*=IR; IR*=inv_r; } tmp=1;IR=1; rep(i,rp.size()){ rp[i]=tmp; tmp*=IR; IR*=r; } mint AA=1; rep(i,f.size()){ A[f.size()-i-1]=f[i]*ir[i]*AA; AA*=a; } rep(i,f.size()+m)B[i]=rp[i]; auto res=A*B; vcans(m); rep(i,m){ ans[i]=res[i+f.size()-1]*ir[i]; } return ans; } template FormalPowerSeriesinterporate(vcx,vcy){ assert(mint::get_mod()); assert(x.size()==y.size()); int n=x.size(); int sz=1; while(sz>tree(sz*2,{1}); rep(i,n)tree[sz+i]={-x[i],1}; DREP(i,sz-1,1)tree[i]=tree[i*2]*tree[i*2+1]; auto gdeval=multipoint_evaluation(tree[1].diff(),x); vca(n);rep(i,n)a[i]=y[i]/gdeval[i]; auto dfs=[&](auto&dfs,int l,int r,int idx)->FormalPowerSeries{ if(r-l==1){ if(l>=n)return {0}; return{a[l]}; } int mid=(l+r)>>1; auto L=dfs(dfs,l,mid,idx*2); auto R=dfs(dfs,mid,r,idx*2+1); return L*tree[idx*2+1]+R*tree[idx*2]; }; auto ans=dfs(dfs,0,sz,1); return ans; } //f(a)=y_0 f(ar)=y_1 .... template vcgeo_interporate(mint a,mint r,vcy){ assert(mint::get_mod()); int n=y.size(); if(n==0)return {}; if(n==1)return {y[0]}; if(r==0){ mint c0=y[1]; mint c1=(y[0]-y[1])/a; return {c0,c1}; } //(x-a)(x-ar)...(x-ar^{n-1}) auto get=[&](auto&get,int n)->FormalPowerSeries{ if(n==0)return {1}; auto down=get(get,n/2); auto res=down*down.substitute(1/r.pow(n/2))*r.pow(1ll*(n/2)*(n/2)); if(n%2)res*={-a*r.pow(n-1),1}; return res; }; vcx(n);x[0]=a;REP(i,1,n)x[i]=x[i-1]*r; auto all_mul=get(get,n); auto gdeval=chirp_z(all_mul.diff(),a,r,n); FormalPowerSeriesA(n);rep(i,n)A[i]=y[i]/gdeval[i]; auto coef=chirp_z(-A,r.inv(),r.inv(),n); FormalPowerSeriesC(n);rep(i,n)C[i]=coef[i]; C=C.substitute(1/a);C/=a; return (C*all_mul).pre(n); } template FormalPowerSeriestaylor_shift(FormalPowerSeriesf,int c){ assert(mint::get_mod()); Binombin;bin.build(f.size()+2); vccp(f.size()+1);cp[0]=1;REP(i,1,f.size()+1)cp[i]=cp[i-1]*c; rep(i,f.size())f[i]*=bin.fact(i); FormalPowerSeriesrc(f.size());rep(i,f.size())rc[i]=cp[f.size()-i-1]*bin.invfact(f.size()-i-1); auto ans=(f*rc); ans=FormalPowerSeries(ans.begin()+f.size()-1,ans.end()); rep(i,f.size())ans[i]*=bin.invfact(i); return ans; } template auto taylor_shift(F f,int c)->FormalPowerSeries{ using mint=typename F::value_type; return taylor_shift(FormalPowerSeries(std::move(f)),c); } #line 4 "A/main.cpp" using mint=StaticModInt<998244353>; using f=FormalPowerSeries; void solve(){ INT(n); STR(S); vcs=stovi(S,"in"s); vccnt(n); f ans(n+1); auto bb=[&](int x){ f re(x+1);rep(i,x+1)re[i]=Binom::C(x,i); return re; }; auto dfs=[&](auto&dfs,int l,int r)->pair{ if(l+1==r){ if(s[l]==0)return {{0,1},{0}}; else return {{0},{0,1}}; } int mid=l+r>>1; auto[a,b]=dfs(dfs,l,mid); auto[c,d]=dfs(dfs,mid,r); ans+=a*d; return {a*bb(r-mid)+c,d*bb(mid-l)+b}; }; dfs(dfs,0,n); dbg(ans); ll res=0;rep(i,ans.size())res^=ans[i].val; PRT(res); } signed main(){ int t=1; // cin >> t; while(t--)solve(); }