#line 1 ".lib/template.hpp"


#include <iostream>
#include <string>
#include <vector>
#include <array>
#include <tuple>
#include <stack>
#include <queue>
#include <deque>
#include <algorithm>
#include <set>
#include <map>
#include <unordered_set>
#include <unordered_map>
#include <bitset>
#include <cmath>
#include <functional>
#include <cassert>
#include <climits>
#include <iomanip>
#include <numeric>
#include <memory>
#include <random>
#include <thread>
#include <chrono>
#define allof(obj) (obj).begin(), (obj).end()
#define range(i, l, r) for(int i=l;i<r;i++)
#define unique_elem(obj) obj.erase(std::unique(allof(obj)), obj.end())
#define bit_subset(i, S) for(int i=S, zero_cnt=0;(zero_cnt+=i==S)<2;i=(i-1)&S)
#define bit_kpop(i, n, k) for(int i=(1<<k)-1,x_bit,y_bit;i<(1<<n);x_bit=(i&-i),y_bit=i+x_bit,i=(!i?(1<<n):((i&~y_bit)/x_bit>>1)|y_bit))
#define bit_kth(i, k) ((i >> k)&1)
#define bit_highest(i) (i?63-__builtin_clzll(i):-1)
#define bit_lowest(i) (i?__builtin_ctzll(i):-1)
#define sleepms(t) std::this_thread::sleep_for(std::chrono::milliseconds(t))
using ll = long long;
using ld = long double;
using ul = uint64_t;
using pi = std::pair<int, int>;
using pl = std::pair<ll, ll>;
using namespace std;

template<typename F, typename S>
std::ostream &operator<<(std::ostream &dest, const std::pair<F, S> &p){
  dest << p.first << ' ' << p.second;
  return dest;
}
template<typename T>
std::ostream &operator<<(std::ostream &dest, const std::vector<std::vector<T>> &v){
  int sz = v.size();
  if(sz==0) return dest;
  for(int i=0;i<sz;i++){
    int m = v[i].size();
    for(int j=0;j<m;j++) dest << v[i][j] << (i!=sz-1&&j==m-1?'\n':' ');
  }
  return dest;
}
template<typename T>
std::ostream &operator<<(std::ostream &dest, const std::vector<T> &v){
  int sz = v.size();
  if(sz==0) return dest;
  for(int i=0;i<sz-1;i++) dest << v[i] << ' ';
  dest << v[sz-1];
  return dest;
}
template<typename T, size_t sz>
std::ostream &operator<<(std::ostream &dest, const std::array<T, sz> &v){
  if(sz==0) return dest;
  for(int i=0;i<sz-1;i++) dest << v[i] << ' ';
  dest << v[sz-1];
  return dest;
}
template<typename T>
std::ostream &operator<<(std::ostream &dest, const std::set<T> &v){
  for(auto itr=v.begin();itr!=v.end();){
    dest << *itr;
    itr++;
    if(itr!=v.end()) dest << ' ';
  }
  return dest;
}
template<typename T, typename E>
std::ostream &operator<<(std::ostream &dest, const std::map<T, E> &v){
  for(auto itr=v.begin();itr!=v.end();){
    dest << '(' << itr->first << ", " << itr->second << ')';
    itr++;
    if(itr!=v.end()) dest << '\n';
  }
  return dest;
}
std::ostream &operator<<(std::ostream &dest, __int128_t value) {
  std::ostream::sentry s(dest);
  if (s) {
    __uint128_t tmp = value < 0 ? -value : value;
    char buffer[128];
    char *d = std::end(buffer);
    do {
      --d;
      *d = "0123456789"[tmp % 10];
      tmp /= 10;
    } while (tmp != 0);
    if (value < 0) {
      --d;
      *d = '-';
    }
    int len = std::end(buffer) - d;
    if (dest.rdbuf()->sputn(d, len) != len) {
      dest.setstate(std::ios_base::badbit);
    }
  }
  return dest;
}
template<typename T>
vector<T> make_vec(size_t sz, T val){return std::vector<T>(sz, val);}
template<typename T, typename... Tail>
auto make_vec(size_t sz, Tail ...tail){
  return std::vector<decltype(make_vec<T>(tail...))>(sz, make_vec<T>(tail...));
}
template<typename T>
vector<T> read_vec(size_t sz){
  std::vector<T> v(sz);
  for(int i=0;i<(int)sz;i++) std::cin >> v[i];
  return v;
}
template<typename T, typename... Tail>
auto read_vec(size_t sz, Tail ...tail){
  auto v = std::vector<decltype(read_vec<T>(tail...))>(sz);
  for(int i=0;i<(int)sz;i++) v[i] = read_vec<T>(tail...);
  return v;
}
void io_init(){
  std::cin.tie(nullptr);
  std::ios::sync_with_stdio(false);
}

#line 1 ".lib/math/mod.hpp"


#line 6 ".lib/math/mod.hpp"
#include <type_traits>
#line 8 ".lib/math/mod.hpp"
#include <ostream>
#line 1 ".lib/math/minior/mod_base.hpp"


#line 4 ".lib/math/minior/mod_base.hpp"
// @param m `1 <= m`
constexpr long long safe_mod(long long x, long long m){
  x %= m;
  if (x < 0) x += m;
  return x;
}
struct barrett{
  unsigned int _m;
  unsigned long long im;
  explicit barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1){}
  unsigned int umod()const{return _m;}
  unsigned int mul(unsigned int a, unsigned int b)const{
    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`
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;
}
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<int n>
constexpr bool is_prime = is_prime_constexpr(n);

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 <int m>
constexpr int primitive_root = primitive_root_constexpr(m);

int ceil_pow2(int n){
  int x = 0;
  while ((1U << x) < (unsigned int)(n)) x++;
  return x;
}
int bsf(unsigned int n){
  return __builtin_ctz(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<long long, long long> inv_gcd(long long a, long long b){
  a = safe_mod(a, b);
  if(a == 0) return {b, 0};
  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;
    auto tmp = s;
    s = t;
    t = tmp;
    tmp = m0;
    m0 = m1;
    m1 = tmp;
  }
  if(m0 < 0) m0 += b / s;
  return {s, m0};
}


#line 13 ".lib/math/mod.hpp"

template<int m>
long long modpow(long long a, long long b){
  assert(0 <= b);
  assert(0 < m);
  a = safe_mod(a, m);
  long long ret = 1;
  while(b){
    if(b & 1) ret = (ret * a) % m;
    a = (a * a) % m;
    b >>= 1;
  }
  return ret;
}
// @param 0 <= b, 0 < m
long long modpow(long long a, long long b, int m){
  assert(0 <= b);
  assert(0 < m);
  a = safe_mod(a, m);
  long long ret = 1;
  while(b){
    if(b & 1) ret = (ret * a) % m;
    a = (a * a) % m;
    b >>= 1;
  }
  return ret;
}

struct modint_base {};
struct static_modint_base : modint_base {};

template <int m, std::enable_if_t<(1 <= m)>* = nullptr>
struct static_modint : static_modint_base{
  using mint = static_modint;
public:
  static constexpr int mod(){return m;}
  static mint raw(int v) {
    mint x;
    x._v = v;
    return x;
  }
  static_modint(): _v(0){}
  template <class T>
  static_modint(T v){
    long long x = v % (long long)umod();
    if (x < 0) x += umod();
    _v = x;
  }
  unsigned int val()const{return _v;}
  mint& operator++(){
    _v++;
    if (_v == umod()) _v = 0;
    return *this;
  }
  mint& operator--(){
    if (_v == 0) _v = umod();
    _v--;
    return *this;
  }
  mint operator++(int){
    mint result = *this;
    ++*this;
    return result;
  }
  mint operator--(int){
    mint result = *this;
    --*this;
    return result;
  }
  mint& operator+=(const mint& rhs){
    _v += rhs._v;
    if (_v >= umod()) _v -= umod();
    return *this;
  }
  mint& operator-=(const mint& rhs){
    _v -= rhs._v;
    if (_v >= umod()) _v += umod();
    return *this;
  }
  mint& operator*=(const mint& rhs){
    unsigned long long z = _v;
    z *= rhs._v;
    _v = (unsigned int)(z % umod());
    return *this;
  }
  mint& operator/=(const mint& rhs){return *this = *this * rhs.inv();}
  mint operator+()const{return *this;}
  mint operator-()const{return mint() - *this;}
  mint pow(long long n)const{
    assert(0 <= n);
    mint x = *this, r = 1;
    while(n){
      if (n & 1) r *= x;
      x *= x;
      n >>= 1;
    }
    return r;
  }
  mint inv()const{
    if(prime){
      assert(_v);
      return pow(umod() - 2);
    }else{
      auto eg = inv_gcd(_v, m);
      assert(eg.first == 1);
      return eg.second;
    }
  }
  friend mint operator+(const mint& lhs, const mint& rhs){return mint(lhs) += rhs;}
  friend mint operator-(const mint& lhs, const mint& rhs){return mint(lhs) -= rhs;}
  friend mint operator*(const mint& lhs, const mint& rhs){return mint(lhs) *= rhs;}
  friend mint operator/(const mint& lhs, const mint& rhs){return mint(lhs) /= rhs;}
  friend bool operator==(const mint& lhs, const mint& rhs){return lhs._v == rhs._v;}
  friend bool operator!=(const mint& lhs, const mint& rhs){return lhs._v != rhs._v;}
private:
  unsigned int _v;
  static constexpr unsigned int umod(){return m;}
  static constexpr bool prime = is_prime<m>;
};

template<int id> 
struct dynamic_modint : modint_base{
  using mint = dynamic_modint;
public:
  static int mod(){return (int)(bt.umod());}
  static void set_mod(int m){
    assert(1 <= m);
    bt = barrett(m);
  }
  static mint raw(int v){
    mint x;
    x._v = v;
    return x;
  }
  dynamic_modint(): _v(0){}
  template <class T>
  dynamic_modint(T v){
    long long x = v % (long long)(mod());
    if (x < 0) x += mod();
    _v = x;
  }
  unsigned int val()const{return _v;}
  mint& operator++(){
    _v++;
    if(_v == umod()) _v = 0;
    return *this;
  }
  mint& operator--(){
    if (_v == 0) _v = umod();
    _v--;
    return *this;
  }
  mint operator++(int){
    mint result = *this;
    ++*this;
    return result;
  }
  mint operator--(int){
    mint result = *this;
    --*this;
    return result;
  }
  mint& operator+=(const mint& rhs){
    _v += rhs._v;
    if(_v >= umod()) _v -= umod();
    return *this;
  }
  mint& operator-=(const mint& rhs){
    _v += mod() - rhs._v;
    if(_v >= umod()) _v -= umod();
    return *this;
  }
  mint& operator*=(const mint& rhs){
    _v = bt.mul(_v, rhs._v);
    return *this;
  }
  mint& operator/=(const mint& rhs){return *this = *this * rhs.inv();}
  mint operator+()const{return *this;}
  mint operator-()const{return mint() - *this;}
  mint pow(long long n)const{
    assert(0 <= n);
    mint x = *this, r = 1;
    while(n){
      if (n & 1) r *= x;
      x *= x;
      n >>= 1;
    }
    return r;
  }
  mint inv()const{
    auto eg = inv_gcd(_v, mod());
    assert(eg.first == 1);
    return eg.second;
  }
  friend mint operator+(const mint& lhs, const mint& rhs){return mint(lhs) += rhs;}
  friend mint operator-(const mint& lhs, const mint& rhs){return mint(lhs) -= rhs;}
  friend mint operator*(const mint& lhs, const mint& rhs){return mint(lhs) *= rhs;}
  friend mint operator/(const mint& lhs, const mint& rhs){return mint(lhs) /= rhs;}
  friend bool operator==(const mint& lhs, const mint& rhs){return lhs._v == rhs._v;}
  friend bool operator!=(const mint& lhs, const mint& rhs){return lhs._v != rhs._v;}
private:
  unsigned int _v;
  static barrett bt;
  static unsigned int umod(){return bt.umod();}
};
template <int id>
barrett dynamic_modint<id>::bt(998244353);
using modint = dynamic_modint<-1>;

using modint998244353 = static_modint<998244353>;
using modint1000000007 = static_modint<1000000007>;
template <class T>
using is_modint = std::is_base_of<modint_base, T>;
template <class T>
using is_modint_t = std::enable_if_t<is_modint<T>::value>;
template <class T>
using is_static_modint = std::is_base_of<static_modint_base, T>;
template <class T>
using is_static_modint_t = std::enable_if_t<is_static_modint<T>::value>;
template <class> struct is_dynamic_modint : public std::false_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint<id>> : public std::true_type {};
template <class T>
using is_dynamic_modint_t = std::enable_if_t<is_dynamic_modint<T>::value>;
template<int m>
std::ostream &operator<<(std::ostream &dest, const static_modint<m> &a){
  dest << a.val();
  return dest;
}
template<int id>
std::ostream &operator<<(std::ostream &dest, const dynamic_modint<id> &a){
  dest << a.val();
  return dest;
}

// 0 <= n < m <= int_max
// 前処理 O(n + log(m))
// 各種計算 O(1)
// 変数 <= n
template<typename mint, is_modint<mint>* = nullptr>
struct modcomb{
private:
  int n;
  std::vector<mint> f, i, fi;
  void init(int _n){
    assert(0 <= _n && _n < mint::mod());
    if(_n < f.size()) return;
    n = _n;
    f.resize(n + 1), i.resize(n + 1), fi.resize(n + 1);
    f[0] = fi[0] = mint(1);
    if(n) f[1] = fi[1] = i[1] = mint(1);
    for(int j = 2; j <= n; j++) f[j] = f[j - 1] * j;
    fi[n] = f[n].inv();
    for(int j = n; j >= 2; j--){
      fi[j - 1] = fi[j] * j;
      i[j] = f[j - 1] * fi[j];
    }
  }
public:
  modcomb(): n(-1){}
  modcomb(int _n){
    init(_n);
  }
  void recalc(int _n){
    init(std::min(mint::mod() - 1, 1 << ceil_pow2(_n)));
  }
  mint comb(int a, int b){
    if((a < 0) || (b < 0) || (a < b)) return 0;
    return f[a] * fi[a - b] * fi[b];
  }
  mint perm(int a, int b){
    if((a < 0) || (b < 0) || (a < b)) return 0;
    return f[a] * fi[a - b];
  }
  mint fac(int x){
    assert(0 <= x && x <= n);
    return f[x];
  }
  mint inv(int x){
    assert(0 < x && x <= n);
    return i[x];
  }
  mint finv(int x){
    assert(0 <= x && x <= n);
    return fi[x];
  }
};
template<typename mint, is_modint<mint>* = nullptr>
struct modpow_table{
  std::vector<mint> v;
  // x^maxkまで計算できる
  modpow_table(){}
  void init(int x, int maxk){
    v.resize(maxk + 1);
    v[0] = 1;
    for(int i = 1; i <= maxk; i++) v[i] = v[i - 1] * x;
  }
  mint pow(int k){
    assert(0 <= k && k < v.size());
    return v[k];
  }
};

#line 1 ".lib/data_structure/range_query/pseudo_tree.hpp"


#line 6 ".lib/data_structure/range_query/pseudo_tree.hpp"

// 0-indexedのセグメントツリーを模した木
template<typename Idx = int>
struct pseudo_segment_tree{
  static constexpr int bitlen = sizeof(Idx) * 8;
  Idx N, M;
  pseudo_segment_tree(){}
  pseudo_segment_tree(Idx n): N(n){
    M = 1;
    while(M < N) M <<= 1;
  }
  // aの深さ
  int depth(Idx a){
    if(bitlen <= 32) return 31 - __builtin_clz(a + 1);
    return 63 - __builtin_clzll(a + 1);
  }
  // aが表す区間の幅
  Idx width(Idx a){
    return M >> depth(a);
  }
  // aが葉か
  bool is_leaf(Idx a){
    return M - 1 <= a;
  }
  // a, bの最短距離
  Idx dist(Idx a, Idx b){
    return depth(a) + depth(b) - 2 * depth(lca(a, b));
  }
  // aのk個親, 深さを超える場合は-1
  Idx la(Idx a, int k){
    if(depth(a) < k) return -1;
    return ((a + 1) >> k) - 1;
  }
  // lca
  Idx lca(Idx a, Idx b){
    a++, b++;
    int da = depth(a), db = depth(b);
    if(da > db) std::swap(a, b), std::swap(da, db);
    b >>= (db - da);
    if(a == b) return a - 1;
    int msb_diff = (bitlen <= 32 ? 31 - __builtin_clz(a ^ b) : 63 - __builtin_clzll(a ^ b)) + 1;
    return (a >> msb_diff) - 1;
  }
  // aが対応する区間
  std::pair<Idx, Idx> index_to_range(Idx a){
    assert(0 <= a && a < 2 * M - 1);
    int dep = depth(a);
    Idx offset = (a + 1) - ((Idx)1 << dep), wid = M >> dep;
    return std::make_pair(offset * wid, (offset + 1) * wid);
  }
  // 区間[l, r)に対応するノード番号(左が先)
  std::vector<Idx> range_to_index(Idx l, Idx r){
    l = std::max(l, 0), r = std::min(r, N);
    assert(l <= r);
    l += M, r += M;
    std::vector<Idx> left, right;
    while(l < r){
      if(l & 1) left.push_back((l++) - 1);
      if(r & 1) right.push_back((--r) - 1);
      l >>= 1;
      r >>= 1;
    }
    std::reverse(right.begin(), right.end());
    left.insert(left.end(), right.begin(), right.end());
    return left;
  }
  // 葉a( < N) から根まで辿るときのノード番号(底が先)
  std::vector<Idx> leaf_to_root(Idx a){
    assert(0 <= a && a < N);
    a += M - 1;
    std::vector<Idx> ret{a};
    while(a){
      a = (a - 1) >> 1;
      ret.push_back(a);
    }
    return ret;
  }
};


// 0-indexedのk分木を模した木
// 頂点iから ki + 1, ki + 2....ki + kに辺が伸びている(nを超える場合はなし)
// (= 頂点iから (i - 1) / kに辺が伸びている(0からはなし))
template<typename Idx, int k>
struct pseudo_k_ary_tree{
  static constexpr int bitlen = sizeof(Idx) * 8;
  Idx N, M;
  std::vector<Idx> Lelem; // 各深さの最左ノード
  std::vector<Idx> kpow;
  pseudo_k_ary_tree(){}
  pseudo_k_ary_tree(Idx n): N(n){
    assert(n);
    M = 1;
    Lelem.push_back(0);
    while(M < N){
      Lelem.push_back(M);
      // Mは最大でNK程度になり, N, kが大きいとMがオーバーフローする可能性がある
      assert((std::numeric_limits<Idx>::max() - 1) / k >= M);
      M = (M * k + 1);
    }
    Idx p = 1;
    for(int i = 0; i < Lelem.size(); i++){
      kpow.push_back(p);
      p *= k;
    }
  }
  int height(){
    return Lelem.size();
  }
  // aの深さ
  int depth(Idx a){
    int ret = 0;
    while(a){
      a = (a - 1) / k;
      ret++;
    }
    return ret;
  }
  // {aの深さ, aと同じ深さのノードでaより小さいものの数}
  std::pair<int, Idx> index_sibling(Idx a){
    int d = depth(a);
    return {d, a - Lelem[d]};
  }
  // 深さが最も深いノードの数
  Idx num_deepest(){
    return N - Lelem.back();
  }
  // 葉の数
  Idx num_leaf(){
    Idx nd = num_deepest();
    Idx ALLLEAF = M - Lelem.back();
    return nd + (ALLLEAF - nd) / k;
  }
  // aの部分木に含まれる最も深いノードの数
  Idx num_subdeepest(Idx a){
    auto [d, si] = index_sibling(a);
    int hdiff = (int)Lelem.size() - d;
    // 完全k分木ならk ^ (h - 1 - d)個の葉がある
    return std::max(Idx(0), num_deepest() - si * kpow[hdiff - 1]);
  }
  // aの部分木に含まれる葉の数
  Idx num_subleaf(Idx a){
    auto [d, si] = index_sibling(a);
    int hdiff = (int)Lelem.size() - d;
    // 完全k分木ならk ^ (h - 1 - d)個の葉がある
    Idx subdeep = std::max(Idx(0), num_deepest() - si * kpow[hdiff - 1]);
    return subdeep + (kpow[hdiff - 1] - subdeep) / k;
  }
  // aの部分木のサイズ
  Idx num_subtree(Idx a){
    auto [d, si] = index_sibling(a);
    int hdiff = (int)Lelem.size() - d;
    // 完全k分木ならk ^ (h - 1 - d)個の葉がある
    Idx subdeep = std::max(Idx(0), num_deepest() - si * kpow[hdiff - 1]);
    return Lelem[hdiff - 1] + subdeep;
  }
  // aが葉か
  bool is_leaf(Idx a){
    return Lelem.back() <= a;
  }
  Idx dist(Idx a, Idx b){
    return dist2(a, b).second;
  }
  // a, bの{lca, 最短距離}
  std::pair<Idx, Idx> dist2(Idx a, Idx b){
    Idx d = 0;
    while(a != b){
      if(a < b) std::swap(a, b);
      a = (a - 1) / k;
      d++;
    }
    return {a, d};
  }
  // 親, ない場合は-1
  Idx parent(Idx a){
    return a ? (a - 1) / k : -1;
  }
  // aのt個親, 深さを超える場合は-1
  Idx la(Idx a, int t){
    for(int i = 0; i < t; i++){
      if(!a) return -1;
      a = (a - 1) / k;
    }
    return a;
  }
  // lca
  Idx lca(Idx a, Idx b){
    return dist2(a, b).first;
  }
  // {ノードの深さ, その部分木の深さh-1のノードがいくつ欠けているか}でノードを分類すると, その種類数は高々3h
  // {個数, ノードの深さ, 深さh-1のノードがいくつ欠けているか}を返す
  std::vector<std::tuple<Idx, Idx, Idx>> depth_frequency_decompose(){
    std::vector<std::tuple<Idx, Idx, Idx>> ret;
    Idx x = N - 1, nd = 1;
    int h = (int)Lelem.size();
    for(int d = h - 1; d >= 0; d--){
      Idx L = x - Lelem[d], R = kpow[d] - 1 - L;
      if(L) ret.push_back({L, d, 0});
      if(R) ret.push_back({R, d, kpow[h - 1 - d]});
      ret.push_back({1, d, kpow[h - 1 - d] - nd});
      if(d){
        x--;
        nd += kpow[h - 1 - d] * (x % k);
        x /= k;
      }
    }
    return ret;
  }
  // aの部分木の頻度テーブル(ans[i] := aの部分木に含まれaとの距離がiのノード数)
  std::vector<Idx> depth_frequency(Idx a){
    if(a >= N) return {};
    std::vector<Idx> ret;
    int t = 0;
    while(a < M){
      if(N <= a){
        ret.push_back(std::max(Idx(0), kpow[t++] - (a - N + 1)));
        return ret;
      }else{
        ret.push_back(kpow[t++]);
      }
      a = a * k + k;
    }
    return ret;
  }
  // aの部分木の頂点でaとの距離がdの頂点の数
  Idx count_dist_subtree(Idx a, int d){
    if(d < 0 || a >= N) return 0;
    int da = depth(a);
    return __count_dist_subtree(a, da, d);
  }
  Idx __count_dist_subtree(Idx a, int da, int d){
    if(d < 0 || a >= N) return 0;
    int h = Lelem.size();
    if(da + d >= h) return 0;
    if(da + d < h - 1){
      return kpow[d];
    }else{
      Idx ldeep = (a - Lelem[da]) * kpow[d];
      return std::min(kpow[d], std::max(Idx(0), num_deepest() - ldeep));
    }
  }
  // aとの距離がdの頂点の数
  Idx count_dist(Idx a, int d){
    if(d < 0 || a >= N) return 0;
    Idx ans = 0;
    int da = depth(a);
    while(d >= 0){
      ans += __count_dist_subtree(a, da, d);
      if(!a) return ans;
      ans -= __count_dist_subtree(a, da, d - 2);
      d--, da--;
      a = (a - 1) / k;
    }
    return ans;
  }
};


#line 4 "a.cpp"
using mint = modint1000000007;

int main(){
  io_init();
  int d, l, r, k;
  std::cin >> d >> l >> r >> k;
  l--, r--;
  pseudo_k_ary_tree<int, 2> t((1 << d) - 1);
  int d1 = t.depth(l), d2 = t.depth(r);
  if(d1 > d2){
    std::swap(l, r);
    std::swap(d1, d2);
  }
  int diff_d = k - (d2 - d1);
  int d3 = d1 - diff_d / 2;
  if(diff_d % 2 == 1 || d1 + k < d2 || d3 < 0){
    std::cout << 0 << '\n';
    return 0;
  }
  mint base = 1;
  modcomb<mint> mcb(1 << d);
  range(i, 0, t.height()){
    int dnum = t.count_dist_subtree(0, i);
    if(i == d1) dnum--;
    if(i == d2) dnum--;
    assert(dnum >= 0);
    base *= mcb.fac(dnum);
  }
  mint ans = 0;
  d1 -= d3, d2 -= d3;
  assert(d2 > 0);
  int ri = (d3 == t.height() - 1 ? (1 << d) : t.Lelem[d3 + 1]);
  for(int i = t.Lelem[d3]; i < ri; i++){
    if(d1 == 0){
      ans += t.__count_dist_subtree(i, d3, d2);
    }else{
      ans += (mint)t.__count_dist_subtree(i * 2 + 1, d3 + 1, d1 - 1) * t.__count_dist_subtree(i * 2 + 2, d3 + 1, d2 - 1);
      ans += (mint)t.__count_dist_subtree(i * 2 + 1, d3 + 1, d2 - 1) * t.__count_dist_subtree(i * 2 + 2, d3 + 1, d1 - 1);
    }
  }
  std::cout << ans * base << '\n';
}