#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; using namespace atcoder; using namespace internal; typedef long long ll; typedef long double LD; typedef double D; typedef pair P; typedef pair PD; typedef map M; #define vp(c,n) vec c(n);rep(i,n){c[i]=i;} #define out(x) cout<= (n); --i) #define drep(i, n) drep2(i, n, 0) #define REP(i, n) FOR(i, 0, n) #define fore(i,a) for(const auto &i:a) #define RREP(var, n) for (auto var = n - 1; var != static_cast(-1); var--)//なんですかこれは #define all(x) (x).begin(),(x).end() #define memset(v, h) memset((v), h, sizeof(v)) #define frac(x) x-floor(x) #define eb emplace_back #define pb push_back #define pbd(x,y) push_back({x,y}); #define pbt(x,y,z) push_back({x,y,z}); #define bs binary_search #define lb lower_bound #define ub upper_bound #define bp(x) __builtin_popcount(x) #define SUM(x) accumulate(all(x),0) #define IP(a,b) inner_product(all(a),begin(b),0.0)//内積の計算 #define PS(a,b) partial_sum(all(a),begin(b))//累積和(具体的にはbは個数以外空にしておいて格納されていく) #define uni(x) sort(all(x));x.erase(unique(all(x)),x.end()) //重複消去&sort #define fr(f,x) fixed<,greater

> #define PQS PQ> #define fi first//大文字になれない #define se second #define cauto const auto& #define bit(n,k) ((n>>k)&1) #define cl clock()/CLOCKS_PER_SEC//時間計測 #define printd(n,x) cout<>a[i];}//試験運転してみる(多分つかわない) #define cinvv(a,n,m); rep(i,n){rep(j,m){cin>>a[i][j];}} #define popcount(n) __builtin_popcount(n) template inline bool chmax(T& a,T b){if(a < b){a=b;return 1;}return 0;} template inline bool chmin(T& a,T b){if(a > b){a=b;return 1;}return 0;} typedef vector vec; typedef vector ves; typedef vector mat; const LD PI=acos(-1); //const ll mod = 1000000009;//手動で切り替える. //const ll mod = 998244353; const ll mod = 1000000007; //const ll mod=1e9; //using mint =modint998244353;//.val()が必要に注意 一旦置いておく //using mint =modint1000000009; using mint =modint1000000007; const ll INF = (1LL<<(60)); //ll dx[4]={0,1,0,-1},dy[4]={-1,0,1,0}; //struct Edge{ll to,co;Edge(ll to=0,ll co=0):to(to),co(co){}};//ABC191で何となくsunukeさんの!! //struct Edge {ll to;ll cost;}; //struct Edge {ll from;ll to;ll cost;}; //struct Edge {ll to;ll w;Edge(ll to,ll w):to(to),w(w){}}; //けんちょん本 G[頂点a].pb(Edge(頂点b,weight));とか //using Graph = vector>; //struct Edge {ll to,t,k;Edge(ll to,ll t,ll k):to(to),t(t),k(k){}};//ABC192で何となく //using Graph =vector>; #define cing(a,b,c) a[b].pb(c);a[c].pb(b); //vector seen; //void gdfs(const Graph &G,ll v){seen[v]=true;fore(next_v,G[v]){if(seen[next_v])continue;gdfs(G,next_v);}}//それぞれの点でif(seen[i])continue;else gdfs(G,i)をする ll modpow(ll a,ll n,ll mod){if(mod==1)return 0;ll ret=1;ll p=a%mod;while(n){if(n&1)ret=ret*p%mod;p=p*p%mod;n>>=1;}return ret; } M factor(ll n) {M ret;for(ll i=2;i*i<=n;i++){while(n%i==0){ret[i]++;n /= i;}}if(n != 1){ret[n]=1;}return ret;}//素因数分解 mapにする vec divisor(ll n){vec K;for(ll i=1;i*i<=n;i++){if(n%i==0){K.pb(i);if(i*i!=n)K.pb(n/i);}}sort(all(K));return K;}//約数列挙 cauto& a=divisor(100);的に書く ll modlog(ll a,ll b,ll p){ll g=1;for(ll i=p;i;i/=2)(g*=a)%=p;g=gcd(g,p);ll t=1,c=0;for(;t%g;c++){if(t==b)return c;(t*=a)%=p;}if(b%g){return -1;}t/=g;b /= g;ll n=p/g,h=0,gs=1;for(;h*hbs;for(ll s=0,e=b;s(NdN2, my_div(N, x2) + 1);ll ix = x;for (ll i = 1; i < imax; ++i) {hl[i] -= (ix < NdN2 ? hl[ix] : hs[my_div(N, ix)]) - pi;ix += x;}for (int n = N2; n >= x2; n--) {hs[n] -= hs[my_div(n, x)] - pi;}++pi;}return hl[1];} // /* constexpr ll MAX = 3000000;ll fac[MAX],finv[MAX],inv[MAX];//int main(){cominit();}と打つのを忘れないように. void cominit(){fac[0]=fac[1]=1;finv[0]=finv[1]=1;inv[1]=1;for(int i=2;iT phi(T n) {T ret=n;for(T i = 2; i * i <= n; i++){if(n % i == 0){ret -= ret / i;while(n % i == 0) n /= i;}}if(n > 1){ret -= ret / n;}return ret;} templateT 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;} mat mat_mul(mat&a,mat&b){mat res(a.size(),vec(b[0].size()));rep(i,a.size()){rep(j,b[0].size()){rep(k,b.size()){(res[i][j]+=a[i][k]*b[k][j])%=mod;}}}return res;} mat matpow(mat a,ll n){mat res(a.size(),vec(a.size()));rep(i,a.size())res[i][i]=1;while(n>0){if(n & 1)res=mat_mul(a,res);a=mat_mul(a,a);n>>=1;}return res;} //#define mant vector> //mant mant_mul(mant&a,mant&b){mant res(a.size(),vector(b[0].size()));rep(i,a.size()){rep(j,b[0].size()){rep(k,b.size()){(res[i][j]+=a[i][k]*b[k][j])}}}return res;} //mant mantpow(mant a,ll n){mant res(a.size(),vector(a.size()));rep(i,a.size())res[i][i]=1;while(n>0){if(n & 1)res=mant_mul(a,res);a=mant_mul(a,a);n>>=1;}return res;} #define cinmat(a,x,y) mat a(x,vec(y));rep(i,x)rep(j,y)cin>>a[i][j];//;をつけなくていいよ!(お気持ちかなり使いやすそう使いにくかったらやめるけども) #define comat(a,x,y) rep(i,x){rep(j,y){cout< //fw(n);としてfw.add(p,x)->fw[p]+x fw.sum(l,r)->Σ[l,r-1]fw[i] (O(logn)!!) #define SA(s) vec suffix_arry(s)//よくわかんない #define LCP(s,a) vec lcp_rray(s,a)//Long Common Prefix #define ZA(s) z_algorithm(s)//上と同じようなやつ #define IV(x,m) inv_mod(x,m)//逆元がなんと対応!! //ついでに pair crt(r,m) 解なし(0,0) n=0 (0,1) CRT!! #define FS(n,m,a,b) floor_sum(n,m,a,b)//Σ[0,n-1]floor((a*i+b)/m) (これがなんとO(log(max(n,m,a,b)))) #define FFT(a,b,c) vec a=convolution(b,c)//FFT,NNT 畳み込み をFFT()の直後に挿入したらmodできる. void ACC(vec a,vec &b){b[0]=0;rrep(i,a.size()){b[i]=b[i-1]+a[i-1];}return;} struct SuccinctIndexableDictionary {size_t length;size_t blocks;vector< unsigned > bit, sum;SuccinctIndexableDictionary() = default;SuccinctIndexableDictionary(size_t length) : length(length), blocks((length + 31) >> 5) {bit.assign(blocks, 0U);sum.assign(blocks, 0U);} void set(int k) {bit[k >> 5] |= 1U << (k & 31);} void build() {sum[0] = 0U;for(int i = 1; i < blocks; i++) {sum[i] = sum[i - 1] + __builtin_popcount(bit[i - 1]);}} bool operator[](int k) {return (bool((bit[k >> 5] >> (k & 31)) & 1));} int rank(int k) {return (sum[k >> 5] + __builtin_popcount(bit[k >> 5] & ((1U << (k & 31)) - 1)));} int rank(bool val, int k) {return (val ? rank(k) : k - rank(k));}}; template< typename T, int MAXLOG > struct WaveletMatrix {size_t length;SuccinctIndexableDictionary matrix[MAXLOG];int mid[MAXLOG];WaveletMatrix() = default;WaveletMatrix(vector< T > v) : length(v.size()) {vector< T > l(length), r(length); for(int level = MAXLOG - 1; level >= 0; level--) {matrix[level] = SuccinctIndexableDictionary(length + 1);int left = 0, right = 0;for(int i = 0; i < length; i++) {if(((v[i] >> level) & 1)) {matrix[level].set(i);r[right++] = v[i];} else {l[left++] = v[i];}} mid[level] = left;matrix[level].build();v.swap(l);for(int i = 0; i < right; i++) {v[left + i] = r[i];}}}pair< int, int > succ(bool f, int l, int r, int level) { return {matrix[level].rank(f, l) + mid[level] * f, matrix[level].rank(f, r) + mid[level] * f};} T access(int k) {T ret = 0;for(int level = MAXLOG - 1; level >= 0; level--) {bool f = matrix[level][k];if(f) ret |= T(1) << level;k = matrix[level].rank(f, k) + mid[level] * f;}return ret;} T operator[](const int &k) {return access(k);}int rank(const T &x, int r) {int l = 0;for(int level = MAXLOG - 1; level >= 0; level--) {tie(l, r) = succ((x >> level) & 1, l, r, level);}return r - l;} T kth_smallest(int l, int r, int k) {assert(0 <= k && k < r - l);T ret = 0; for(int level = MAXLOG - 1; level >= 0; level--) {int cnt = matrix[level].rank(false, r) - matrix[level].rank(false, l);bool f = cnt <= k;if(f) {ret |= T(1) << level;k -= cnt;}tie(l, r) = succ(f, l, r, level);}return ret;} T kth_largest(int l, int r, int k) {return kth_smallest(l, r, r - l - k - 1);} int range_freq(int l, int r, T upper) {int ret = 0; for(int level = MAXLOG - 1; level >= 0; level--) { bool f = ((upper >> level) & 1);if(f) {ret += matrix[level].rank(false, r) - matrix[level].rank(false, l);}tie(l, r) = succ(f, l, r, level);}return ret;} int range_freq(int l, int r, T lower, T upper) { return range_freq(l, r, upper) - range_freq(l, r, lower);} T prev_value(int l, int r, T upper) {int cnt = range_freq(l, r, upper);return cnt == 0 ? T(-1) : kth_smallest(l, r, cnt - 1);} T next_value(int l, int r, T lower) {int cnt = range_freq(l, r, lower);return cnt == r - l ? T(-1) : kth_smallest(l, r, cnt);}}; template< typename T, int MAXLOG > struct CWM {WaveletMatrix< int, MAXLOG > mat;vector< T > ys; CWM(const vector< T > &v) : ys(v) {sort(begin(ys), end(ys));ys.erase(unique(begin(ys), end(ys)), end(ys));vector< int > t(v.size());for(int i = 0; i < v.size(); i++) t[i] = get(v[i]);mat = WaveletMatrix< int, MAXLOG >(t);} inline int get(const T& x) {return lower_bound(begin(ys), end(ys), x) - begin(ys);} T access(int k) {return ys[mat.access(k)];} T operator[](const int &k) {return access(k);} int rank(const T &x, int r) {auto pos = get(x);if(pos == ys.size() || ys[pos] != x) {return 0;}return mat.rank(pos, r);} T kth_smallest(int l, int r, int k) {return ys[mat.kth_smallest(l, r, k)];} T kth_largest(int l, int r, int k) {return ys[mat.kth_largest(l, r, k)];} int range_freq(int l, int r, T upper) {return mat.range_freq(l, r, get(upper));} int range_freq(int l, int r, T lower, T upper) {return mat.range_freq(l, r, get(lower), get(upper));} T prev_value(int l, int r, T upper) {auto ret = mat.prev_value(l, r, get(upper));return ret == -1 ? T(-1) : ys[ret];} T next_value(int l, int r, T lower) {auto ret = mat.next_value(l, r, get(lower));return ret == -1 ? T(-1) : ys[ret];}}; void cum(vec &a,vec b){ ll n=b.size(); a.resize(n); rep(i,n-1){a[i+1]=a[i]+b[i];} } void cum1(vec &a,vec b){ ll n=b.size(); a.resize(n-1); a[0]=b[0]; rrep(i,n-1){a[i]=a[i-1]+b[i];} } namespace inner { using i32 = int32_t; using u32 = uint32_t; using i64 = int64_t; using u64 = uint64_t; template T gcd(T a, T b) { while (b) swap(a %= b, b); return a; } template T inv(T a, T p) { T b = p, x = 1, y = 0; while (a) { T q = b / a; swap(a, b %= a); swap(x, y -= q * x); } assert(b == 1); return y < 0 ? y + p : y; } template T modpow(T a, U n, T p) { T ret = 1 % p; for (; n; n >>= 1, a = U(a) * a % p) if (n & 1) ret = U(ret) * a % p; return ret; } } // namespace inner namespace my_rand { // [0, 2^64 - 1) uint64_t rng() { static uint64_t x_ = uint64_t(chrono::duration_cast( chrono::high_resolution_clock::now().time_since_epoch()) .count()) * 10150724397891781847ULL; x_ ^= x_ << 7; return x_ ^= x_ >> 9; } // [l, r) int64_t randint(int64_t l, int64_t r) { assert(l < r); return l + rng() % (r - l); } // choose n numbers from [l, r) without overlapping vector randset(int64_t l, int64_t r, int64_t n) { assert(l <= r && n <= r - l); unordered_set s; for (int64_t i = n; i; --i) { int64_t m = randint(l, r + 1 - i); if (s.find(m) != s.end()) m = r - i; s.insert(m); } vector ret; for (auto& x : s) ret.push_back(x); return ret; } // [0.0, 1.0) double rnd() { union raw_cast { double t; uint64_t u; }; constexpr uint64_t p = uint64_t(1023 - 64) << 52; return rng() * ((raw_cast*)(&p))->t; } template void randshf(vector& v) { int n = v.size(); for (int loop = 0; loop < 2; loop++) for (int i = 0; i < n; i++) swap(v[i], v[randint(0, n)]); } } // namespace my_rand using my_rand::randint;using my_rand::randset;using my_rand::randshf;using my_rand::rnd;using my_rand::rng; struct ArbitraryLazyMontgomeryModInt { using mint = ArbitraryLazyMontgomeryModInt; using i32 = int32_t; using u32 = uint32_t; using u64 = uint64_t; static u32 mod; static u32 r; static u32 n2; static u32 get_r() { u32 ret = mod; for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret; return ret; } static void set_mod(u32 m) { assert(m < (1 << 30)); assert((m & 1) == 1); mod = m; n2 = -u64(m) % m; r = get_r(); assert(r * mod == 1); } u32 a; ArbitraryLazyMontgomeryModInt() : a(0) {} ArbitraryLazyMontgomeryModInt(const int64_t &b) : a(reduce(u64(b % mod + mod) * n2)){}; static u32 reduce(const u64 &b) { return (b + u64(u32(b) * u32(-r)) * mod) >> 32; } mint &operator+=(const mint &b) { if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod; return *this; } mint &operator-=(const mint &b) { if (i32(a -= b.a) < 0) a += 2 * mod; return *this; } mint &operator*=(const mint &b) { a = reduce(u64(a) * b.a); return *this; } mint &operator/=(const mint &b) { *this *= b.inverse(); return *this; } mint operator+(const mint &b) const { return mint(*this) += b; } mint operator-(const mint &b) const { return mint(*this) -= b; } mint operator*(const mint &b) const { return mint(*this) *= b; } mint operator/(const mint &b) const { return mint(*this) /= b; } bool operator==(const mint &b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); } bool operator!=(const mint &b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); } mint operator-() const { return mint() - mint(*this); } mint pow(u64 n) const { mint ret(1), mul(*this); while (n > 0) { if (n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } friend ostream &operator<<(ostream &os, const mint &b) { return os << b.get(); } friend istream &operator>>(istream &is, mint &b) { int64_t t; is >> t; b = ArbitraryLazyMontgomeryModInt(t); return (is); } mint inverse() const { return pow(mod - 2); } u32 get() const { u32 ret = reduce(a); return ret >= mod ? ret - mod : ret; } static u32 get_mod() { return mod; } }; typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::mod;typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::r;typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::n2; struct montgomery64 { using mint = montgomery64; using i64 = int64_t; using u64 = uint64_t; using u128 = __uint128_t; static u64 mod; static u64 r; static u64 n2; static u64 get_r() { u64 ret = mod; for (i64 i = 0; i < 5; ++i) ret *= 2 - mod * ret; return ret; } static void set_mod(u64 m) { assert(m < (1LL << 62)); assert((m & 1) == 1); mod = m; n2 = -u128(m) % m; r = get_r(); assert(r * mod == 1); } u64 a; montgomery64() : a(0) {} montgomery64(const int64_t &b) : a(reduce((u128(b) + mod) * n2)){}; static u64 reduce(const u128 &b) { return (b + u128(u64(b) * u64(-r)) * mod) >> 64; } mint &operator+=(const mint &b) { if (i64(a += b.a - 2 * mod) < 0) a += 2 * mod; return *this; } mint &operator-=(const mint &b) { if (i64(a -= b.a) < 0) a += 2 * mod; return *this; } mint &operator*=(const mint &b) { a = reduce(u128(a) * b.a); return *this; } mint &operator/=(const mint &b) { *this *= b.inverse(); return *this; } mint operator+(const mint &b) const { return mint(*this) += b; } mint operator-(const mint &b) const { return mint(*this) -= b; } mint operator*(const mint &b) const { return mint(*this) *= b; } mint operator/(const mint &b) const { return mint(*this) /= b; } bool operator==(const mint &b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); } bool operator!=(const mint &b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); } mint operator-() const { return mint() - mint(*this); } mint pow(u128 n) const { mint ret(1), mul(*this); while (n > 0) { if (n & 1) ret *= mul; mul *= mul; n >>= 1; } return ret; } friend ostream &operator<<(ostream &os, const mint &b) { return os << b.get(); } friend istream &operator>>(istream &is, mint &b) { int64_t t; is >> t; b = montgomery64(t); return (is); } mint inverse() const { return pow(mod - 2); } u64 get() const { u64 ret = reduce(a); return ret >= mod ? ret - mod : ret; } static u64 get_mod() { return mod; } }; typename montgomery64::u64 montgomery64::mod, montgomery64::r, montgomery64::n2; namespace fast_factorize { using u64 = uint64_t; template bool miller_rabin(u64 n, vector as) { if (mint::get_mod() != n) mint::set_mod(n); u64 d = n - 1; while (~d & 1) d >>= 1; mint e{1}, rev{int64_t(n - 1)}; for (u64 a : as) { if (n <= a) break; u64 t = d; mint y = mint(a).pow(t); while (t != n - 1 && y != e && y != rev) { y *= y; t *= 2; } if (y != rev && t % 2 == 0) return false; } return true; } bool is_prime(u64 n) { if (~n & 1) return n == 2; if (n <= 1) return false; if (n < (1LL << 30)) return miller_rabin(n, {2, 7, 61}); else return miller_rabin( n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022}); } template T pollard_rho(T n) { if (~n & 1) return 2; if (is_prime(n)) return n; if (mint::get_mod() != n) mint::set_mod(n); mint R, one = 1; auto f = [&](mint x) { return x * x + R; }; auto rnd_ = [&]() { return rng() % (n - 2) + 2; }; while (1) { mint x, y, ys, q = one; R = rnd_(), y = rnd_(); T g = 1; constexpr int m = 128; for (int r = 1; g == 1; r <<= 1) { x = y; for (int i = 0; i < r; ++i) y = f(y); for (int k = 0; g == 1 && k < r; k += m) { ys = y; for (int i = 0; i < m && i < r - k; ++i) q *= x - (y = f(y)); g = inner::gcd(q.get(), n); } } if (g == n) do g = inner::gcd((x - (ys = f(ys))).get(), n); while (g == 1); if (g != n) return g; } exit(1); } using i64 = long long; vector inner_factorize(u64 n) { if (n <= 1) return {}; u64 p; if (n <= (1LL << 30)) p = pollard_rho(n); else p = pollard_rho(n); if (p == n) return {i64(p)}; auto l = inner_factorize(p); auto r = inner_factorize(n / p); copy(begin(r), end(r), back_inserter(l)); return l; } vector factorize(u64 n) { auto ret = inner_factorize(n); sort(begin(ret), end(ret)); return ret; } map factor_count(u64 n) { map mp; for (auto &x : factorize(n)) mp[x]++; return mp; } vector divisors(u64 n) { if (n == 0) return {}; vector> v; for (auto &p : factorize(n)) { if (v.empty() || v.back().first != p) { v.emplace_back(p, 1); } else { v.back().second++; } } vector ret; auto f = [&](auto rc, int i, i64 x) -> void { if (i == (int)v.size()) { ret.push_back(x); return; } for (int j = v[i].second;; --j) { rc(rc, i + 1, x); if (j == 0) break; x *= v[i].first; } }; f(f, 0, 1); sort(begin(ret), end(ret)); return ret; } } // namespace fast_factorize using fast_factorize::divisors;using fast_factorize::factor_count;using fast_factorize::factorize;using fast_factorize::is_prime; using namespace fast_factorize; vec base(ll x,ll b) {vec ret;ll t = 1, k = abs(b); while(x) {ret.emplace_back((x * t) % k);if(ret.back() < 0){ret.back() += k;}x -= ret.back() * t;x /= k;t *= b / k;} if(ret.empty()) {ret.emplace_back(0);}reverse(begin(ret), end(ret));return ret;} struct Barrett { using u32 = unsigned int; using i64 = long long; using u64 = unsigned long long; u32 m; u64 im; Barrett() : m(), im() {} Barrett(int n) : m(n), im(u64(-1) / m + 1) {} constexpr inline i64 quo(u64 n) { u64 x = u64((__uint128_t(n) * im) >> 64); u32 r = n - x * m; return m <= r ? x - 1 : x; } constexpr inline i64 rem(u64 n) { u64 x = u64((__uint128_t(n) * im) >> 64); u32 r = n - x * m; return m <= r ? r + m : r; } constexpr inline pair quorem(u64 n) { u64 x = u64((__uint128_t(n) * im) >> 64); u32 r = n - x * m; if (m <= r) return {x - 1, r + m}; return {x, r}; } constexpr inline i64 pow(u64 n, i64 p) { u32 a = rem(n), r = m == 1 ? 0 : 1; while (p) { if (p & 1) r = rem(u64(r) * a); a = rem(u64(a) * a); p >>= 1; } return r; } }; struct prime_power_binomial { int p, q, M; vector fac, ifac, inv; int delta; Barrett bm, bp; prime_power_binomial(int _p, int _q) : p(_p), q(_q) { assert(1 < p && p <= ((1LL << 30) - 1)); assert(_q > 0); long long m = 1; while (_q--) { m *= p; assert(m <=((1LL << 30) - 1)); } M = m; bm = Barrett(M), bp = Barrett(p); enumerate(); delta = (p == 2 && q >= 3) ? 1 : M - 1; } void enumerate() { int MX = min(M, 20000000 + 10); fac.resize(MX); ifac.resize(MX); inv.resize(MX); fac[0] = ifac[0] = inv[0] = 1; fac[1] = ifac[1] = inv[1] = 1; for (int i = 2; i < MX; i++) { if (i % p == 0) { fac[i] = fac[i - 1]; fac[i + 1] = bm.rem(1LL * fac[i - 1] * (i + 1)); i++; } else { fac[i] = bm.rem(1LL * fac[i - 1] * i); } } ifac[MX - 1] = bm.pow(fac[MX - 1], M / p * (p - 1) - 1); for (int i = MX - 2; i > 1; --i) { if (i % p == 0) { ifac[i] = bm.rem(1LL * ifac[i + 1] * (i + 1)); ifac[i - 1] = ifac[i]; i--; } else { ifac[i] = bm.rem(1LL * ifac[i + 1] * (i + 1)); } } } long long Lucas(long long n, long long m) { int res = 1; while (n) { int n0, m0; tie(n, n0) = bp.quorem(n); tie(m, m0) = bp.quorem(m); if (n0 < m0) return 0; res = bm.rem(1LL * res * fac[n0]); int buf = bm.rem(1LL * ifac[n0 - m0] * ifac[m0]); res = bm.rem(1LL * res * buf); } return res; } long long C(long long n, long long m) { if (n < m || n < 0 || m < 0) return 0; if (q == 1) return Lucas(n, m); long long r = n - m; int e0 = 0, eq = 0, i = 0; int res = 1; while (n) { res = bm.rem(1LL * res * fac[bm.rem(n)]); res = bm.rem(1LL * res * ifac[bm.rem(m)]); res = bm.rem(1LL * res * ifac[bm.rem(r)]); n = bp.quo(n); m = bp.quo(m); r = bp.quo(r); int eps = n - m - r; e0 += eps; if (e0 >= q) return 0; if (++i >= q) eq += eps; } if (eq & 1) res = bm.rem(1LL * res * delta); res = bm.rem(1LL * res * bm.pow(p, e0)); return res; } }; struct binomial_mod { int mod; vector M; vector cs; binomial_mod(long long md) : mod(md) { assert(1 <= md); assert(md <= ((1LL << 30) - 1)); for (int i = 2; i * i <= md; i++) { if (md % i == 0) { int j = 0, k = 1; while (md % i == 0) md /= i, j++, k *= i; M.push_back(k); cs.emplace_back(i, j); assert(M.back() == cs.back().M); } } if (md != 1) { M.push_back(md); cs.emplace_back(md, 1); } assert(M.size() == cs.size()); } long long binom(long long n, long long m) { if (mod == 1) return 0; vector rem, d; for (int i = 0; i < (int)cs.size(); i++) { rem.push_back(cs[i].C(n, m)); d.push_back(M[i]); } return atcoder::crt(rem, d).first; } }; #define ip(x) is_prime_constexpr(x) ll m=10; mat mmat_mul(mat&a,mat&b){mat res(a.size(),vec(b[0].size()));rep(i,a.size()){rep(j,b[0].size()){rep(k,b.size()){(res[i][j]+=a[i][k]*b[k][j])%=m;}}}return res;} mat mmatpow(mat a,ll n){mat res(a.size(),vec(a.size()));rep(i,a.size())res[i][i]=1;while(n>0){if(n & 1)res=mmat_mul(a,res);a=mmat_mul(a,a);n>>=1;}return res;} int main(){ ll p,q,r,k; cin>>p>>q>>r>>k; mat g={{1,1,1},{1,0,0},{0,1,0}}; auto b=mmatpow(g,k-3); cout<<(b[0][0]*r+b[0][1]*q+b[0][2]*p)%10<