#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; using ll = long long; using ull = unsigned long long; using ld = long double; #define vall(A) A.begin(), A.end() #define vin(A) for (ll iiii = 0, sz = A.size(); iiii < sz; iiii++){cin >> A[iiii];} #define vout(A) for(ll iiii = 0, sz = A.size(); iiii < sz; iiii++){cout << A[iiii] << " \n"[iiii == sz-1];} #define vout2d(A) for (ll iiii = 0, HHHH = A.size(), WWWW = (A.empty() ? 0 : A[0].size()); iiii < HHHH; iiii++){for (ll jjjj = 0; jjjj < WWWW; jjjj++){cout << A[iiii][jjjj] << " \n"[jjjj==WWWW-1];}} #define vsum(A) [&](const auto &vveecc){ll ssuumm = 0; for(auto &vvaalluuee : vveecc){ssuumm += vvaalluuee;} return ssuumm;}(A) #define adjustedvin(A) for (ll iiii = 1, sz = A.size(); iiii < sz; iiii++){cin >> A[iiii];} #define adjustedvout(A) for(ll iiii = 1, sz = A.size(); iiii < sz; iiii++){cout << A[iiii] << " \n"[iiii == sz-1];} #define adjustedvout2d(A) for (ll iiii = 1, HHHH = A.size(), WWWW = (A.empty() ? 0 : A[0].size()); iiii < HHHH; iiii++){for (ll jjjj = 1; jjjj < WWWW; jjjj++){cout << A[iiii][jjjj] << " \n"[jjjj==WWWW-1];}} vector 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}; vector pow10ll{1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000,10000000000,100000000000,1000000000000,10000000000000,100000000000000,1000000000000000,10000000000000000,100000000000000000,1000000000000000000, 10000000000000000000ull}; vector di{0,1,0,-1}; vector dj{1,0,-1,0}; namespace a1073741824{ //library namespace math{ //math.hpp /// @brief a^bをmで割った余りを返す。bに関して対数時間で計算できる。 /// @param a /// @param b /// @param m /// @return a^b%m ll modpow(ll a, ll b, ll m){ ll t = a%m; ll ans = 1; while (b > 0){ if (b%2){ ans = (ans*t)%m; } b /= 2; t = (t*t)%m; } return ans; } /// @brief a^nを返す。bに関して線形時間で計算できる。 /// @param a /// @param n /// @param m /// @return a^n ll powll(ll a, ll n){ ll r = 1; for (ll i = 1; i <= n; i++){ r *= a; } return r; } /// @brief 正の整数Nを素因数分解する /// @param N /// @return vector {{素因数1,個数}, {素因数2,個数}, {素因数3,個数}...} vector> p_fact(ll N){ if (N == 1){ return vector> {{1,1}}; } vector> R;//戻り値用リスト const ll M = N; for (ll i = 2; i <= sqrtl(M); i++){ if (N % i == 0){ ll divide_count = 0; while (N % i == 0){ divide_count++; N /= i; } R.push_back(vector {i,divide_count}); } } if (N != 1){ R.push_back(vector {N,1}); } return R; } /// @brief 二次元ベクトル表記の素因数分解リストを受け取って約数の和を求める /// @param vv /// @return 約数の和 ll sum_of_divisor(vector> vv){ if (vv[0][0] == 1){ return 1; } ll R = 1; for (vector x : vv){ R *= ((ll)powl(x[0],x[1]+1) - 1)/(x[0] - 1); } return R; } /// @brief 二次元ベクトル表記の素因数分解リストを受け取って約数の種類を求める /// @param vv /// @return 約数の種類 ll kind_of_divisor(vector> vv){ ll R = 1; for (vector x : vv){ R *= x[1]+1; } return R; } /// @brief 1のみを含むデフォルトリストdと、素因数分解の結果pを受け取って、dを約数リストに変換する。 /// @return ない void search_divisor(vector &d, vector> &p, ll depth = 0){ if (depth == p.size()){ sort(vall(d)); return; } ll t = d.size(); for (ll i = 1; i <= p[depth][1]; i++){ for (ll j = 0; j < t; j++){ d.push_back(d[j]*powll(p[depth][0], i)); } } search_divisor(d, p, depth+1); } /// @brief 有理数のceilを求める /// @param y /// @param x /// @return ceil(y/x) ll cefrac(ll y, ll x); /// @brief 有理数のfloorを求める /// @param y /// @param x /// @return floor(y/x) ll flfrac(ll y, ll x){ if ((x^y) > 0){ x = abs(x); y = abs(y); return (y-(y%x))/x; } else if ((x^y) < 0){ x = abs(x); y = abs(y); return -cefrac(y,x); } else{ return y/x; } } ll cefrac(ll y, ll x){ if ((x^y) > 0){ x = abs(x); y = abs(y); if (y%x == 0){ return y/x; } else return 1 + (y-(y%x))/x; } else if ((x^y) < 0){ x = abs(x); y = abs(y); return -flfrac(y,x); } else{ return y/x; } } ll flsqrt(ll N){ if (N){ ll ok = 1; ll ng = min(N,2000000000LL); while (ng - ok >= 2){ ll mid = (ok+ng)/2; if (mid*mid <= N){ ok = mid; } else{ ng = mid; } } return ok; } else{return 0;} } /// @brief エラトステネスの篩 struct Eratosthenes{ vector P_List; vector is_prime; vector min_factor; Eratosthenes(ll N) : is_prime(N+1,1), min_factor(N+1,-1){ is_prime[1] = 0; min_factor[1] = 1; for (ll i = 2; i <= N; i++){ if (is_prime[i]){ P_List.push_back(i); min_factor[i] = i; for (ll j = 2*i; j <= N; j += i){ is_prime[j] = 0; if (min_factor[j] == -1){ min_factor[j] = i; } } } } } void chase_prime(const ll &reference, ll x, vector>> &r){ if (r[reference].empty() || min_factor[x] != r[reference].back()[0]){ r[reference].push_back({min_factor[x], 1}); } else{ r[reference].back()[1]++; } if (x != min_factor[x]){ chase_prime(reference, x/min_factor[x], r); } } vector>> p_fact(ll N){ vector>> r(N+1); r[1].push_back({1,1}); for (ll i = 2; i <= N; i++){ chase_prime(i, i, r); } return r; } }; /// @brief 一次不定方程式ax+by=1の解を1つ見つける /// @param a /// @param b /// @return {x,y} pair axby1(ll a, ll b){ if (a == 1 or b == 1){ if (a == 1){ return {1-b,1}; } else{ return {1,1-a}; } } else if (a>b){ auto p = axby1(a%b, b); return {p.first, p.second - p.first*(a/b)}; } else{ swap(a,b); auto p = axby1(a%b, b); return {p.second - p.first*(a/b), p.first}; } } /// @brief 1/a mod Mを求める /// @param a /// @param M /// @return 1/a mod M ll inverse_mod(ll a, ll M){ if (__gcd(a,M) != 1){ return -1; } return (M+(axby1(a,M).first)%M)%M; } /// @brief modint template struct mll{ ll val = 0; mll(const ll &x){ val = x%M; } mll(){ val = 0; } void operator=(const ll &x){ val = x%M; } void operator=(const mll &a){ val = a.val; } mll operator+(const mll &a){ return mll(val+a.val); } void operator+=(const mll &a){ val = (val+a.val)%M; } mll operator+(const ll &a){ return mll(val+a); } void operator+=(const ll &a){ val = (val+a)%M; } mll operator-(const mll &a){ return mll(M+val-a.val); } void operator-=(const mll &a){ val = (M+val-a.val)%M; } mll operator-(const ll &a){ return mll(M+val-a); } void operator-=(const ll &a){ val = (M+val-a)%M; } mll operator*(const mll &a){ return mll(val*a.val); } void operator*=(const mll &a){ val = (val*a.val)%M; } mll operator*(const ll &a){ return mll(val*a); } void operator*=(const ll &a){ val = (val*a)%M; } mll operator/(const mll &a){ return mll(val*inverse_mod(a.val,M)); } void operator/=(const mll &a){ val = (val*inverse_mod(a.val,M))%M; } mll operator/(const ll &a){ return mll(val*inverse_mod(a,M)); } void operator/=(const ll &a){ val = (val*inverse_mod(a,M))%M; } }; //階乗前計算による二項係数mod998244353 struct factorialncr{ private: vector factorialmod; private: ll N_MAX_N_MAX; public: factorialncr(const ll &N_MAX){ N_MAX_N_MAX = N_MAX; factorialmod = vector(N_MAX+1); factorialmod[0] = 1; for (ll i = 1; i <= N_MAX; i++){ factorialmod[i] = (i*factorialmod[i-1])%998244353; } } ll nCr(ll n, ll r){ if (r < 0 || n < r || n > N_MAX_N_MAX){ return 0; } return (factorialmod[n]*inverse_mod((factorialmod[n-r]*factorialmod[r])%998244353,998244353))%998244353; } }; //表の前計算による二項係数modM struct tablencr{ private: vector> ncrmodlist; private: ll N_MAX_N_MAX; public: tablencr(const ll &N_MAX, const ll &M){ N_MAX_N_MAX = N_MAX; ncrmodlist = vector> (5001, vector(5001,0)); ncrmodlist[0][0] = 1; for (ll i = 1; i <= 5000; i++){ for (ll j = 0; j <= i; j++){ if (j == 0 || j == i){ ncrmodlist[i][j] = 1; } else{ ncrmodlist[i][j] = (ncrmodlist[i-1][j-1] + ncrmodlist[i-1][j])%M; } } } } ll nCr(ll n, ll r){ if (r < 0 || n < r || n > N_MAX_N_MAX){ return 0; } return ncrmodlist[n][r]; } }; //matrix.hpp /// @brief 64bit整数行列。自動でmod998244353などを取ってくれる。 /// @tparam MMOODD template struct matrixll{ vector> M; ll H,W; /// @brief N次単位行列を生成 /// @param N matrixll(ll N){ H = N;W = N; M = vector>(N,vector(N,0)); for (ll i = 0; i < N; i++){ M[i][i] = 1; } } /// @brief h×w型の、全要素がvの行列を生成 /// @param h /// @param w /// @param v matrixll(ll h, ll w, ll v){ H = h; W = w; M = vector>(H,vector(W,v)); } /// @brief 2次元配列を用いて行列を生成 /// @param A matrixll(vector> &A){ M = A; H = A.size(); W = A[0].size(); } matrixll operator+(const matrixll &T){ if (H != T.H || W != T.W){ cerr << "size error\n"; assert(false); } matrixll ans(M); for (ll i = 0; i < H; i++){ for (ll j = 0; j < W; j++){ ans.M[i][j] += T.M[i][j]; ans.M[i][j] %= MMOODD; } } return ans; } void operator+=(const matrixll &T){ if (H != T.H || W != T.W){ cerr << "size error\n"; assert(false); } for (ll i = 0; i < H; i++){ for (ll j = 0; j < W; j++){ M[i][j] += T.M[i][j]; M[i][j] %= MMOODD; } } } matrixll operator-(const matrixll &T){ if (H != T.H || W != T.W){ cerr << "size error\n"; assert(false); } matrixll ans(M); for (ll i = 0; i < H; i++){ for (ll j = 0; j < W; j++){ ans.M[i][j] -= T.M[i][j]; ans.M[i][j] += MMOODD; ans.M[i][j] %= MMOODD; } } return ans; } void operator-=(const matrixll &T){ if (H != T.H || W != T.W){ cerr << "size error\n"; assert(false); } for (ll i = 0; i < H; i++){ for (ll j = 0; j < W; j++){ M[i][j] -= T.M[i][j]; M[i][j] += MMOODD; M[i][j] %= MMOODD; } } } matrixll operator*(const matrixll &T){ if (W != T.H){ cerr << "size error\n"; assert(false); } matrixll ans(H,T.W,0); for (ll i = 0; i < H; i++){ for (ll j = 0; j < T.W; j++){ for (ll k = 0; k < W; k++){ ans.M[i][j] += M[i][k]*T.M[k][j]; ans.M[i][j] %= MMOODD; } } } return ans; } void operator*=(const matrixll &T){ if (W != T.H){ cerr << "size error\n"; assert(false); } matrixll ans(H,T.W,0); for (ll i = 0; i < H; i++){ for (ll j = 0; j < T.W; j++){ for (ll k = 0; k < W; k++){ ans.M[i][j] += M[i][k]*T.M[k][j]; ans.M[i][j] %= MMOODD; } } } W = T.W; M = ans.M; } matrixll operator*(const ll &c){ matrixll ans(H,W,0); for (ll i = 0; i < H; i++){ for (ll j = 0; j < W; j++){ ans.M[i][j] = M[i][j] * c; ans.M[i][j] %= MMOODD; } } return ans; } void operator*=(const ll &c){ for (ll i = 0; i < H; i++){ for (ll j = 0; j < W; j++){ M[i][j] *= c; M[i][j] %= MMOODD; } } } /// @brief A^Nを返す。 /// @param A /// @param N /// @return A^N matrixll matrixmodpow(ll N){ matrixll R(H); matrixll A(M); for (ll i = 0; i < H; i++){ for (ll j = 0; j < H; j++){ A.M[i][j] %= MMOODD; } } if (N){ while (N){ if (N%2){ R *= A; } A *= A; N /= 2; } return R; } else{ return R; } } }; //quatient_range.hpp /// @brief sum[i=1 ~ i=N] i^a * (floor(M/i))^b を計算する /// @param N /// @param M /// @param a /// @param b /// @return ll inv_floor_sum(ll N, ll M, ll a ,ll b){ auto floorsqrt = [](ll N){ ll ok = 1; ll ng = min(N,2000000000LL); while (ng - ok >= 2){ ll mid = (ok+ng)/2; if (mid*mid <= N){ ok = mid; } else{ ng = mid; } } return ok; }; auto lpowl = [](ll x, ll N){ ll r = 1; for (ll i = 1; i <= N; i++){ r *= x; } return r; }; vector> nCrtable(a+2,vector(a+2,0)); nCrtable[0][0] = 1; for (ll i = 1; i <= a+1; i++){ for (ll j = 0; j <= i; j++){ if (j == 0 || j == i){ nCrtable[i][j] = 1; } else{ nCrtable[i][j] = nCrtable[i-1][j-1] + nCrtable[i-1][j]; } } } function sum = [&](ll n, ll l){ if (l == 0){ return n; } if (l == 1){ return n*(n+1)/2; } if (l == 2){ return n*(n+1)*(2*n+1)/6; } ll T = 0; for (ll i = 0; i < l; i++){ T += nCrtable[l+1][i]*sum(n,i); } return (lpowl(n+1,l+1)-1-T)/(l+1); }; ll ans = 0; ll rootM = floorsqrt(M); if (N <= M / rootM){//Nが小さいときはそのまま計算して値を返す for (ll i = 1; i <= N; i++){ ans += lpowl(i,a)*lpowl(M/i,b); } return ans; } //それ以外の場合は主客転倒を行う for (ll i = 1; i <= M/rootM; i++){//最初の方をそのまま足す。 ans += lpowl(i,a)*lpowl(M/i,b); } for (ll i = rootM-1; i > M/N; i--){ ans += lpowl(i,b)*(sum(M/i,a) - sum(M/(i+1),a)); } ans += lpowl(M/N,b)*(sum(N,a) - sum(M/((M/N)+1), a)); return ans; } //convolution.hpp vector powroot998244353{1LL, 998244352LL, 911660635LL, 372528824LL, 69212480LL, 381091786LL, 515872972LL, 273395035LL, 469477650LL, 503794470LL, 464513504LL, 558899782LL, 504969456LL, 840897051LL, 539927980LL, 417009993LL, 725461291LL, 456548895LL, 712494065LL, 542639453LL, 768214923LL, 303507821LL, 438914733LL, 761881641}; vector powrootinv998244353{1LL, 998244352LL, 86583718LL, 509520358LL, 661054123LL, 863921598LL, 323451984LL, 689146517LL, 379690232LL, 240519260LL, 899368279LL, 920065956LL, 588792270LL, 118574449LL, 847593593LL, 858760106LL, 987418793LL, 177938570LL, 753608159LL, 786906984LL, 331540181LL, 655706071LL, 268754442LL, 191076546}; vector powroot1224736769{1LL, 1224736768LL, 24506215LL, 992888437LL, 853017291LL, 235319402LL, 269744380LL, 157861287LL, 894223137LL, 600648668LL, 1091208103LL, 382541006LL, 12818149LL, 149218560LL, 746299392LL, 405692663LL, 633223136LL, 672327338LL, 992917013LL, 758198491LL, 1079610480LL, 1056667043LL, 1039331205LL, 1026803890LL, 449603200}; vector powrootinv1224736769{1LL, 1224736768LL, 1200230554LL, 961581489LL, 796105727LL, 1023008969LL, 406386483LL, 251411652LL, 655108271LL, 1145368249LL, 780593535LL, 38041180LL, 816166160LL, 659160240LL, 1200430513LL, 392462252LL, 15561184LL, 893027826LL, 928760284LL, 497993173LL, 529117122LL, 637457654LL, 931394937LL, 823596420LL, 55047034}; vector DFT998244353(vector X, ll K, bool inverse = false){ if (K == 1){ return vector{(X[0]+X[1])%998244353, (998244353+X[0]-X[1])%998244353}; } vector even(1LL<<(K-1)); vector odd(1LL<<(K-1)); for (ll i = 0; i < (1LL<<(K-1)); i++){ even[i] = (X[i] + X[(1LL<<(K-1))+i])%998244353; } ll temp = 1; if (inverse){ for (ll i = 0; i < (1LL<<(K-1)); i++){ odd[i] = (temp*(998244353 + X[i] - X[(1LL<<(K-1))+i]))%998244353; temp = (temp*powrootinv998244353[K])%998244353; } } else{ for (ll i = 0; i < (1LL<<(K-1)); i++){ odd[i] = (temp*(998244353 + X[i] - X[(1LL<<(K-1))+i]))%998244353; temp = (temp*powroot998244353[K])%998244353; } } even = DFT998244353(even,K-1,inverse); odd = DFT998244353(odd,K-1,inverse); for (ll i = 0; i < (1LL< DFT1224736769(vector X, ll K, bool inverse = false){ if (K == 1){ return vector{(X[0]+X[1])%1224736769, (1224736769+X[0]-X[1])%1224736769}; } vector even(1LL<<(K-1)); vector odd(1LL<<(K-1)); for (ll i = 0; i < (1LL<<(K-1)); i++){ even[i] = (X[i] + X[(1LL<<(K-1))+i])%1224736769; } ll temp = 1; if (inverse){ for (ll i = 0; i < (1LL<<(K-1)); i++){ odd[i] = (temp*(1224736769 + X[i] - X[(1LL<<(K-1))+i]))%1224736769; temp = (temp*powrootinv1224736769[K])%1224736769; } } else{ for (ll i = 0; i < (1LL<<(K-1)); i++){ odd[i] = (temp*(1224736769 + X[i] - X[(1LL<<(K-1))+i]))%1224736769; temp = (temp*powroot1224736769[K])%1224736769; } } even = DFT1224736769(even,K-1,inverse); odd = DFT1224736769(odd,K-1,inverse); for (ll i = 0; i < (1LL< convolution998244353(vector A, vector B){ if (min(A.size(), B.size()) <= 60) { if (A.size() < B.size()) { swap(A, B); } std::vector ans(A.size() + B.size() - 1); for (int i = 0; i < A.size(); i++) { for (int j = 0; j < B.size(); j++) { ans[i + j] += A[i] * B[j]; } } for (auto &v : ans){ v %= 998244353; } return ans; } ll N = A.size()+B.size()-1; ll log2N = 0; while ((1LL< C((1LL< convolution1224736769(vector A, vector B){ if (min(A.size(), B.size()) <= 60) { if (A.size() < B.size()) { swap(A, B); } std::vector ans(A.size() + B.size() - 1); for (int i = 0; i < A.size(); i++) { for (int j = 0; j < B.size(); j++) { ans[i + j] += A[i] * B[j]; } } for (auto &v : ans){ v %= 1224736769; } return ans; } ll N = A.size()+B.size()-1; ll log2N = 0; while ((1LL< C((1LL< tree; vector B; ll log2N = 0; void fill_tree(vector &ruisekiwaA, ll L, ll R, ll index){ if (index >= pow2ll[log2N]){return;} if (index == 0){ tree[index] = ruisekiwaA.back(); fill_tree(ruisekiwaA,L,(L+R)/2,1); return; } tree[index] = ruisekiwaA[R] - (L == 0 ? 0 : ruisekiwaA[L-1]); fill_tree(ruisekiwaA,L,(L+R)/2,2*index); fill_tree(ruisekiwaA,R+1,(3*R-L+1)/2,2*index+1); } FenwickTree(ll howmany, ll value){ while (pow2ll[log2N] < howmany){ log2N++; } tree = vector(pow2ll[log2N], value); B = vector(pow2ll[log2N], value); for (ll i = pow2ll[log2N-1]-1; i >= 1; i--){ tree[i] = tree[2*i]*2; } tree[0] = 2*tree[1]; } FenwickTree(ll howmany, vector &A){ while (pow2ll[log2N] < howmany){ log2N++; } vector ruisekiwaA(pow2ll[log2N]); ruisekiwaA[0] = A[0]; for (ll i = 1; i < pow2ll[log2N]; i++){ if (i < howmany){ ruisekiwaA[i] = A[i] + ruisekiwaA[i-1]; } else{ ruisekiwaA[i] = ruisekiwaA[i-1]; } } tree = vector(pow2ll[log2N],0); B = A; while (B.size() < pow2ll[log2N]){ B.push_back(0); } fill_tree(ruisekiwaA,0,pow2ll[log2N]-1,0); } ll sum_from_origin(ll index){ index++; if (index == pow2ll[log2N]){ return tree[0]; } vector temp(log2N+1); for (ll i = 0; i <= log2N; i++){ temp[i] = index%2; index /= 2; } reverse(temp.begin(), temp.end()); vector temp2; temp2.push_back(0); temp2.push_back(1); for (ll i = 2; i <= log2N; i++){ temp2.push_back(0); temp2[i] = temp2[i-1]*2 + temp[i-1]; } ll sum = 0; for (ll i = 0; i <= log2N; i++){ sum += temp[i] * tree[temp2[i]]; } return sum; } /// @brief [l,r]の総和を返す。 /// @param l /// @param r /// @return 総和 ll range_sum(ll l, ll r){ if (l == 0){ return sum_from_origin(r); } return sum_from_origin(r)-sum_from_origin(l-1); } void add(ll index, ll value, ll L = 0, ll d = 0, ll index_on_tree = 0){ if (index_on_tree == 0){ tree[index_on_tree] += value; add(index,value,0,1,1); return; } if (d == log2N+1){ return; } ll R = L + pow2ll[log2N-d]-1; if (L <= index && index <= R){ tree[index_on_tree] += value; add(index,value,L,d+1,2*index_on_tree); } else{ add(index,value,R+1,d+1,2*index_on_tree+1); } } /// @brief 一か所変える。 /// @param index void update(ll index, ll value){ add(index,value-B[index]); B[index] = value; } /// @brief 一か所をa倍してbを足す。 /// @param index /// @param a /// @param b void update(ll index, ll a, ll b){ add(index,(a-1)*B[index]+b); B[index] = B[index]*a + b; } }; //SegTree&LazySegTree.hpp /// @brief 抽象化セグメントツリー /// @attention コンストラクタ1 SegTree(A, e, op, mapping) /// @attention コンストラクタ2 SegTree(N, I, e, op, mapping) /// @tparam info セグ木の各ノードに載せる情報をまとめた構造体の型 /// @tparam func 更新に使う変数をまとめた構造体の型(アフィン変換なら、aとbを持つ構造体など) /// @param e 載せたものの単位元(sumなら0, maxなら-infなど) /// @param operation 各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) /// @param mapping infoに対してfuncを作用させた結果を返す関数(アフィン変換ならx -> ax+b) template struct SegTree{ int log2N;//セグ木の高さ-1 info e;///単位元 function operation;//各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) function mapping;//更新を行うとどうなるか?(アフィン変換ならx -> ax+b) /// @brief セグ木を扱うためのノード struct SegNode{ pair range;//自身のカバーする範囲 info I;//自身が持っている情報 SegNode(const int &L, const int &R, const info &eee) : range({L,R}), I(eee){} SegNode(){};//空コンストラクタ }; vector tree;//セグ木本体 /// @brief N個のIで初期化 /// @param I 載せたい構造体 /// @param N 載せた個数 /// @param eee 載せたものの単位元(sumなら0, maxなら-infなど) /// @param op 各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) /// @param m 更新を行うとどうなるか?(アフィン変換ならx -> ax+b) SegTree(int N, info I, info eee, function op, function m){ //基本情報を登録 e = eee; operation = op; mapping = m; //セグ木のサイズを決定 log2N = 0; while ((1<(1<<(log2N+1)); for (int i = 0; i < N; i++){ tree[i+(1<= 1; i--){ tree[i] = SegNode(tree[2*i].range.first,tree[2*i+1].range.second,operation(tree[2*i].I, tree[2*i+1].I)); } } /// @brief vector Aで初期化 /// @param I 載せたい構造体 /// @param N 載せた個数 /// @param eee 載せたものの単位元(sumなら0, maxなら-infなど) /// @param op 各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) /// @param m 更新を行うとどうなるか?(アフィン変換ならx -> ax+b) SegTree(vector &A, info eee, function op, function m){ //基本情報を登録 e = eee; operation = op; mapping = m; //セグ木のサイズを決定 log2N = 0; int sz = A.size(); while ((1ULL<(1<<(log2N+1)); for (int i = 0; i < sz; i++){ tree[i+(1<= 1; i--){ tree[i] = SegNode(tree[2*i].range.first,tree[2*i+1].range.second,operation(tree[2*i].I, tree[2*i+1].I)); } } /// @brief 区間の共通部分が空集合であるか判定 /// @return range1 ⋀ {a,a} == Ø private: bool not_intersect(const pair &range1, const int &a, const int &b){ return max(range1.second,b)-min(range1.first,a) > range1.second+b-range1.first-a; } /// @brief range1がrange2に完全に覆われているかを判定 /// @return range1 ⊂ {a,b} private: bool completely_covered(const pair &range1, const int &a, const int &b){ return a <= range1.first && range1.second <= b; } public: /// @brief index閉区間[L,R]において、集計を行う。 /// @param L 左端(左端を含む) /// @param R 右端(右端を含む) /// @return [L,R]での集計結果 info range_get(const int &L, const int &R){ info ret = e; int left = L; while (left < R+1){ int log2interval = min(__builtin_ctz(left),31-__builtin_clz(R+1-left)); ret = operation(ret, tree[(left+(1<>log2interval].I); left += 1<>= 1; while (start_index){ tree[start_index].I = operation(tree[2*start_index].I,tree[2*start_index+1].I); start_index >>= 1; } } /// @brief 左端をLに固定したとき、条件式Gが成り立つ最小の右端indexを返す。もしなければINF(=2147483647)が返ってくる。 /// @attention 判定関数Gは、区間を広げていったときにfalse,false,false,...false,true,true,true...のように、falseが続いた後にtrueが続くものでなければならない。 /// @param L 左端 /// @param G 判定関数...引数a,bがあり、引数aを動かし、引数bを比較対象tに固定した時に引数aによってtrue,falseが変動する。 /// @param t 比較対象 /// @return Gがtrueになる最小右端indexまたは2147483647 int min_right(int L, const function &G, const info &t){ info current_result = e; checkpoint: int ctz = L == 0 ? log2N : __builtin_ctz(L); if (!G(operation(current_result, tree[((1<>ctz].I), t)){ if (tree[((1<>ctz].range.second+1 == 1<>ctz].I); L = tree[((1<>ctz].range.second+1; goto checkpoint; } for (int i = ctz-1; i >= 0; i--){ if (!G(operation(current_result, tree[((1<>i].I), t)){ current_result = operation(current_result, tree[((1<>i].I); L = tree[((1<>i].range.second+1; goto checkpoint; } } return L; } /// @brief 右端をRに固定したとき、条件式Gが成り立つ最大の左端indexを返す。もしなければ-INF-1(=-2147483648)が返ってくる。 /// @attention 判定関数Gは、区間を広げていったときにfalse,false,false,...false,true,true,true...のように、falseが続いた後にtrueが続くものでなければならない。 /// @param R 右端 /// @param G 判定関数...引数a,bがあり、引数aを動かし、引数bを比較対象tに固定した時に引数aによってtrue,falseが変動するようなものである。 /// @param t 比較対象 /// @return Gがtrueになる最大左端index int max_left(int R, const function &G, const info &t){ info current_result = e; checkpoint: int cto = __builtin_ctz(~R);//cto...count trailing one if (!G(operation(current_result, tree[((1<>cto].I), t)){ if (tree[((1<>cto].range.first == 0){ return -2147483648; } current_result = operation(current_result, tree[((1<>cto].I); R = tree[((1<>cto].range.first-1; goto checkpoint; } for (int i = cto-1; i >= 0; i--){ if (!G(operation(current_result, tree[((1<>i].I), t)){ current_result = operation(current_result, tree[((1<>i].I); R = tree[((1<>i].range.first-1; goto checkpoint; } } return R; } }; /// @brief 抽象化遅延セグメントツリー /// @attention コンストラクタ1 LazySegTree(A, e, op, mapping, composition, id) /// @attention コンストラクタ2 LazySegTree(N, I, e, op, mapping, composition , id) /// @tparam info セグ木の各ノードに載せる情報をまとめた構造体の型 /// @tparam func 更新に使う変数をまとめた構造体の型(アフィン変換なら、aとbを持つ構造体など) /// @param e 載せたものの単位元(sumなら0, maxなら-infなど) /// @param operation 各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) /// @param mapping infoに対してfuncを作用させた結果を返す関数(アフィン変換ならx -> ax+b) /// @param composition func同士の合成結果を1つのfuncにする関数((u,v)でu(v())の結果を返すものとする)(ax+bのあとにcx+dを作用させると実質acx+bc+dになるなど) /// @param id funcの恒等写像(アフィン変換ならx -> 1x+0) template struct LazySegTree{ int log2N;//セグ木の高さ+1 info e;///単位元 function operation;//各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) function mapping;//更新を行うとどうなるか?(アフィン変換ならx -> ax+b) function composition;//mを合成した結果(ax+bのあとにcx+dを作用させると実質acx+bc+dになるなど) func id;//mappingの恒等写像版(アフィン変換ならx -> 1x+0など) /// @brief セグ木を扱うためのノード struct SegNode{ pair range;//自身のカバーする範囲 info I;//自身が持っている情報 func delay;//遅延情報として残っている写像 bool is_delay = 0;//遅延情報があるかどうか SegNode(int L, int R, info eee, func ididid) : range({L,R}), I(eee), delay(ididid){} SegNode(){};//空コンストラクタ }; vector tree;//セグ木本体 /// @brief N個のIで初期化 /// @param I 載せたい構造体 /// @param N 載せた個数 /// @param eee 載せたものの単位元(sumなら0, maxなら-infなど) /// @param op 各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) /// @param m 更新を行うとどうなるか?(アフィン変換ならx -> ax+b) /// @param c mを合成した結果((u,v)でu(v())の結果を返すものとする)(ax+bのあとにcx+dを作用させると実質acx+bc+dになるなど) /// @param ididid mの恒等写像版(アフィン変換ならx -> 1x+0) LazySegTree(int N, info I, info eee, function op, function m, function c, func ididid){ //基本情報を登録 e = eee; operation = op; mapping = m; composition = c; id = ididid; //セグ木のサイズを決定 log2N = 0; while ((1<(1<<(log2N+1)); for (int i = 0; i < N; i++){ tree[i+(1<= 1; i--){ tree[i] = SegNode(tree[2*i].range.first,tree[2*i+1].range.second,operation(tree[2*i].I, tree[2*i+1].I),id); } } /// @brief vector Aで初期化 /// @param I 載せたい構造体 /// @param N 載せた個数 /// @param eee 載せたものの単位元(sumなら0, maxなら-infなど) /// @param op 各ノードに載ってる構造体に対する二項演算をする関数(min,max,sumなど) /// @param m 更新を行うとどうなるか?(アフィン変換ならx -> ax+b) /// @param c mを合成した結果((u,v)でu(v())の結果を返すものとする)(ax+bのあとにcx+dを作用させると実質acx+bc+dになるなど) /// @param ididid mの恒等写像版(アフィン変換ならx -> 1x+0) LazySegTree(const vector &A, info eee, function op, function m, function c, func ididid){ //基本情報を登録 e = eee; operation = op; mapping = m; composition = c; id = ididid; //セグ木のサイズを決定 log2N = 0; int sz = A.size(); while ((1ULL<(1<<(log2N+1)); for (int i = 0; i < sz; i++){ tree[i+(1<= 1; i--){ tree[i] = SegNode(tree[2*i].range.first,tree[2*i+1].range.second,operation(tree[2*i].I, tree[2*i+1].I),id); } } /// @brief 区間の共通部分が空集合であるか判定 /// @return range1 ⋀ {a,b} == Ø private: bool not_intersect(const pair &range1, const int &a, const int &b){ return max(range1.second,b)-min(range1.first,a) > range1.second+b-range1.first-a; } /// @brief range1がrange2に完全に覆われているかを判定 /// @return range1 ⊂ {a,b} private: bool completely_covered(const pair &range1, const int &a, const int &b){ return a <= range1.first && range1.second <= b; } public: /// @brief 遅延情報の伝播を行う void tell_info(int index_on_tree){ if (!tree[index_on_tree].is_delay){return;} if (index_on_tree >= (1< lazy_node;//区間集計、区間更新に使う。 deque lazy_node_flipped;//区間更新で使う。 /// @brief index閉区間[L,R]において、集計を行う。 /// @param L 左端(左端を含む) /// @param R 右端(右端を含む) /// @return [L,R]での集計結果 info range_get(const int &L, const int &R){ int Lstart = L + (1<> 1; int rm = (Rstart / (Rstart & -Rstart)) >> 1; while (Lstart < Rstart){ if (Rstart <= rm){ lazy_node.push_back(Rstart); } if (Lstart <= lm){ lazy_node.push_back(Lstart); } Lstart >>= 1; Rstart >>= 1; } while (Lstart){ lazy_node.push_back(Lstart); Lstart >>=1; } while (!lazy_node.empty()){ tell_info(lazy_node.back()); lazy_node.pop_back(); } info ret = e; int left = L; while (left < R+1){ int log2interval = min(__builtin_ctz(left),31-__builtin_clz(R+1-left)); ret = operation(ret, tree[(left+(1<>log2interval].I); left += 1<> 1; int rm = (Rstart / (Rstart & -Rstart)) >> 1; while (Lstart < Rstart){ if (Rstart <= rm){ lazy_node.push_back(Rstart); } if (Lstart <= lm){ lazy_node.push_back(Lstart); } Lstart >>= 1; Rstart >>= 1; } while (Lstart){ lazy_node.push_back(Lstart); Lstart >>=1; } while (!lazy_node.empty()){ tell_info(lazy_node.back()); lazy_node_flipped.push_back(lazy_node.back()); lazy_node.pop_back(); } int left = L; while (left < R+1){ int log2interval = min(__builtin_ctz(left),31-__builtin_clz(R+1-left)); tree[(left+(1<>log2interval].I = mapping(F,tree[(left+(1<>log2interval].I); tree[(left+(1<>log2interval].delay = composition(F,tree[(left+(1<>log2interval].delay); tree[(left+(1<>log2interval].is_delay = 1; left += 1<= (1< &G, const info &t){ info current_result = e; int ctz_init = L == 0 ? log2N : __builtin_ctz(L); for (int i = log2N; i > ctz_init; i--){ tell_info(((1<>i); } checkpoint: int ctz = L == 0 ? log2N : __builtin_ctz(L); tell_info(((1<>ctz); if (!G(operation(current_result, tree[((1<>ctz].I), t)){ if (tree[((1<>ctz].range.second+1 == 1<>ctz].I); L = tree[((1<>ctz].range.second+1; goto checkpoint; } for (int i = ctz-1; i >= 0; i--){ tell_info(((1<>i); if (!G(operation(current_result, tree[((1<>i].I), t)){ current_result = operation(current_result, tree[((1<>i].I); L = tree[((1<>i].range.second+1; goto checkpoint; } } return L; } /// @brief 右端をRに固定したとき、条件式Gが成り立つ最大の左端indexを返す。もしなければ-INF-1(=-2147483648)が返ってくる。 /// @attention 判定関数Gは、区間を広げていったときにfalse,false,false,...false,true,true,true...のように、falseが続いた後にtrueが続くものでなければならない。 /// @param R 右端 /// @param G 判定関数...引数a,bがあり、引数aを動かし、引数bを比較対象tに固定した時に引数aによってtrue,falseが変動するようなものである。 /// @param t 比較対象のinfo /// @return Gがtrueになる最大左端indexまたは-2147483648 int max_left(int R, const function &G, const info &t){ info current_result = e; int cto_init = __builtin_ctz(~R); for (int i = log2N; i > cto_init; i--){ tell_info(((1<>i); } checkpoint: int cto = __builtin_ctz(~R);//cto...count trailing one tell_info(((1<>cto); if (!G(operation(current_result, tree[((1<>cto].I), t)){ if (tree[((1<>cto].range.first == 0){ return -2147483648; } current_result = operation(current_result, tree[((1<>cto].I); R = tree[((1<>cto].range.first-1; goto checkpoint; } for (int i = cto-1; i >= 0; i--){ tell_info(((1<>i); if (!G(operation(current_result, tree[((1<>i].I), t)){ current_result = operation(current_result, tree[((1<>i].I); R = tree[((1<>i].range.first-1; goto checkpoint; } } return R; } }; ///Trie.hpp /// @brief ただのTrie木 struct Trie{ struct TrieNode{ char c;//このノードの文字 bool forbidden = false;//禁止されたかどうか ll stringnum = 0;//何文字分このノードで終わっているか。 unordered_map child;//子ノードへの参照 TrieNode(char d, ll n){ c = d; stringnum = n; } TrieNode(){} }; TrieNode root;//trie木の根。文字を持たない('.'を持っている) Trie(){ root = TrieNode('.', 0); } /// @brief 文字列を追加する。追加しようとしたけどダメだったらfalseが返ってくる。 /// @param S /// @return 追加できたか。 bool addstring(const string &S){ TrieNode *N = &root; for (ll i = 0; i < S.size(); i++){ if (N->child.count(S[i])){ if (N->child[S[i]]->forbidden){ return false; } } else{ N->child[S[i]] = new TrieNode(S[i],0); } if (i == S.size()-1){ N->child[S[i]]->stringnum++; } N = N->child[S[i]]; } return true; } /// @brief 接頭辞Sを禁止する。 /// @param S /// @return 禁止によっていくつの文字が消えたか。 ll forbid(const string &S){ TrieNode *N = &root; for (ll i = 0; i < S.size(); i++){ if (N->child.count(S[i])){ if (N->child[S[i]]->forbidden){ return 0; } } else{ N->child[S[i]] = new TrieNode(S[i],0); } if (i == S.size()-1){ N->child[S[i]]->forbidden = true; } N = N->child[S[i]]; } ll forbidden_string = 0;//DFSでこれより下を全部禁止する deque Q; Q.push_back(N); while (!Q.empty()){ TrieNode *node = Q.back(); Q.pop_back(); forbidden_string += node->stringnum; for (auto v : node->child){ Q.push_back(v.second); } } N->child.clear(); N->stringnum = 0; return forbidden_string; } }; /// } } using namespace a1073741824; /// @brief n以下の非負整数であって、popcountがちょうどkであるようなものの個数 /// @param n /// @param k /// @param memo_g /// @return ll g(const ll &n, const ll &k, map,ll> &memo_g){ if (memo_g.count({k,n})){ return memo_g[{k,n}]; } if (n <= 0){ if (k == 0){ return 1; } return 0; } if (k < 0){ return 0; } if (n%2){ ll ans = g(n>>1, k, memo_g) + g(n>>1, k-1, memo_g); ans %= 998244353; memo_g[{k,n}] = ans; return ans; } else{ ll ans = g(n>>1, k, memo_g) + g((n>>1)-1, k-1, memo_g);// - (__builtin_popcountll(n>>1) == k-1); ans %= 998244353; memo_g[{k,n}] = ans; return ans; } } /// @brief n以下の数であって、b==n, a<=n,p(b)=l,p(a)=kであるようなa,bに関するa&bの総和を求める。 /// @param n /// @param k /// @param l /// @param memo_g /// @return ll h(ll n, const ll &k, const ll &l, map,ll> &memo_g, map,ll> &memo_h){ if (memo_h.count({k,l,n})){ return memo_h[{k,l,n}]; } if (__builtin_popcountll(n) != l){ return 0; } if (n < 1){ return 0; } if (k < 0 || l < 0){ return 0; } if (n%2){ ll ans = g(n>>1, k-1, memo_g) + 2*h(n>>1, k-1, l-1, memo_g, memo_h) + 2*h(n>>1, k, l-1, memo_g, memo_h); ans %= 998244353; memo_h[{k,l,n}] = ans; return ans; } else{ ll ans = 2*h(n>>1, k-1, l, memo_g, memo_h) - n*(__builtin_popcountll(n>>1) == k-1) + 2*h(n>>1, k, l, memo_g, memo_h); ans %= 998244353; memo_h[{k,l,n}] = ans; return ans; } } ll f(const ll &n, const ll &k, const ll &l, map,ll> &memo_f, map,ll> &memo_g, map,ll> &memo_h){ if (memo_f.count({k,l,n})){ return memo_f[{k,l,n}]; } if (n == 1){ if (k == 1 && l == 1){ return 1; } return 0; } if (k < 0 || l < 0 || n <= 0){ return 0; } if (n%2){ ll ans = 2*f(n>>1, k, l, memo_f, memo_g, memo_h) + g(n>>1, k-1, memo_g)*g(n>>1, l-1, memo_g) + 2*f(n>>1, k-1, l-1, memo_f, memo_g, memo_h) + 2*f(n>>1, k, l-1, memo_f, memo_g, memo_h) + 2*f(n>>1, k-1, l, memo_f, memo_g, memo_h); ans %= 998244353; memo_f[{k,l,n}] = ans; return ans; } else{ ll ans = 2*f(n>>1, k, l, memo_f, memo_g, memo_h) + (g((n>>1)-1, k-1, memo_g)*g((n>>1)-1, l-1, memo_g))%998244353 + 2*f((n>>1)-1, k-1, l-1, memo_f, memo_g, memo_h) + 2*f(n>>1, k, l-1, memo_f, memo_g, memo_h) - 2*h(n>>1, k, l-1, memo_g, memo_h) + 2*f(n>>1, k-1, l, memo_f, memo_g, memo_h) - 2*h(n>>1, l, k-1, memo_g, memo_h); ans %= 998244353; memo_f[{k,l,n}] = ans; return ans; } } int main(){ ios::sync_with_stdio(false); std::cin.tie(nullptr); map,ll> memo_f; map,ll> memo_g; map,ll> memo_h; ll N; cin >> N; ll ans = 0; // for (ll i = 1; i <= 60; i++){ ans += f(N,i,i,memo_f,memo_g,memo_h); } // ans -= ((((N%998244353)*((N+1)%998244353))%998244353)*499122177)%998244353; ans = (998244353+(ans%998244353))%998244353; ans = (ans*499122177)%998244353; ans += ((((N%998244353)*((N+1)%998244353))%998244353)*499122177)%998244353; ans = (998244353+(ans%998244353))%998244353; //ll N = 1024; //cout << f(N,2,2,memo_f,memo_g,memo_h) << "\n"; // //ll ans = 0; //for (ll i = 1; i <= N; i++){ // for (ll j = 1; j <= N; j++){ // if (__builtin_popcountll(i) == __builtin_popcountll(j) && __builtin_popcountll(j) == 2){ // ans += (i&j); // ans %= 998244353; // } // } //} cout << ans << "\n"; }