library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub maspypy/library

:heavy_check_mark: convex/lattice_point_sum_polynomial_pq.hpp

Depends on

Verified with

Code

#include "seq/geometric_sequence_sum_formula.hpp"
#include "convex/line_min_function.hpp"
#include "convex/lattice_point_sum_polynomial.hpp"
#include "mod/floor_sum_of_linear_polynomial.hpp"
#include "mod/floor_sum_of_linear_polynomial_pq.hpp"

// 格子点 (x,y) に対して p^xq^yx^iy^j の sum. i<=K, j<=L
template <typename mint, int K1, int K2>
array<array<mint, K2 + 1>, K1 + 1> lattice_point_sum_polynomial_pq(mint p, mint q, vc<tuple<ll, ll, ll>> LINE) {
  assert(p != 0 && q != 0);
  if (p == 1 && q == 1) return lattice_point_sum_polynomial<mint, K1, K2>(LINE);
  if (q == 1) {
    for (auto& [a, b, c]: LINE) { swap(a, b); }
    auto tmp = lattice_point_sum_polynomial_pq<mint, K2, K1>(q, p, LINE);
    array<array<mint, K2 + 1>, K1 + 1> ANS{};
    FOR(i, K1 + 1) FOR(j, K2 + 1) ANS[i][j] = tmp[j][i];
    return ANS;
  }
  assert(q != 1);
  ll L = -infty<ll>, R = infty<ll>;
  vc<tuple<ll, ll, ll>> LINE1, LINE2;
  for (auto& [a, b, c]: LINE) {
    if (b == 0) {
      assert(a != 0);
      if (a > 0) { chmin(R, floor<ll>(c, a) + 1); }
      elif (a < 0) { chmax(L, ceil<ll>(-c, -a)); }
    } else {
      if (b > 0) { LINE2.eb(-a, c, b); }
      if (b < 0) { LINE1.eb(a, -c, -b); }
    }
  }
  if (L >= R) { return {}; }
  if (LINE1.empty() || LINE2.empty()) return {};

  auto LOWER = line_max_function_rational(LINE1, L, R);
  auto UPPER = line_min_function_rational(LINE2, L, R);

  array<array<mint, K2 + 1>, K1 + 1> S{};
  array<mint, K1 + 1> T{};

  bool bad = 0;
  auto wk = [&](ll L, ll R, ll a1, ll b1, ll c1, ll a2, ll b2, ll c2) -> void {
    if (bad) return;
    // 交点 t/s
    i128 s = i128(a2) * c1 - i128(a1) * c2;
    i128 t = i128(b1) * c2 - i128(b2) * c1;
    if (s == 0) {
      if (t > 0) return;
    }
    if (s > 0) { chmax(L, ceil<i128>(t, s)); }
    if (s < 0) { chmin(R, floor<i128>(-t, -s) + 1); }
    if (L >= R) return;
    if (L == -infty<ll> || R == infty<ll>) {
      bad = 1;
      return;
    }
    b2 = b2 + c2, b1 = b1 - 1 + c1;
    auto ADD_S = floor_sum_of_linear_polynomial_pq<mint, K1, K2, ll>(p, q, L, R, a2, b2, c2);
    auto SUB_S = floor_sum_of_linear_polynomial_pq<mint, K1, K2, ll>(p, q, L, R, a1, b1, c1);
    auto ADD_T = floor_sum_of_linear_polynomial_pq<mint, K1, 0, ll>(p, 1, L, R, a2, b2, c2);
    auto SUB_T = floor_sum_of_linear_polynomial_pq<mint, K1, 0, ll>(p, 1, L, R, a1, b1, c1);
    FOR(i, K1 + 1) FOR(j, K2 + 1) S[i][j] += ADD_S[i][j] - SUB_S[i][j];
    FOR(i, K1 + 1) T[i] += ADD_T[i][0] - SUB_T[i][0];
  };

  merge_58(LOWER, UPPER, wk);
  array<array<mint, K2 + 1>, K1 + 1> ANS{};
  if (bad) return ANS;

  FOR(k, K2 + 1) {
    auto [c, f] = geometric_sequence_sum_formula(q, k);
    FOR(i, K1 + 1) {
      ANS[i][k] += c * T[i];
      FOR(j, k + 1) { ANS[i][k] += f[j] * S[i][j]; }
    }
  }

  return ANS;
}
#line 1 "convex/lattice_point_sum_polynomial_pq.hpp"

#line 2 "poly/count_terms.hpp"
template<typename mint>
int count_terms(const vc<mint>& f){
  int t = 0;
  FOR(i, len(f)) if(f[i] != mint(0)) ++t;
  return t;
}
#line 2 "mod/modint_common.hpp"

struct has_mod_impl {
  template <class T>
  static auto check(T &&x) -> decltype(x.get_mod(), std::true_type{});
  template <class T>
  static auto check(...) -> std::false_type;
};

template <class T>
class has_mod : public decltype(has_mod_impl::check<T>(std::declval<T>())) {};

template <typename mint>
mint inv(int n) {
  static const int mod = mint::get_mod();
  static vector<mint> dat = {0, 1};
  assert(0 <= n);
  if (n >= mod) n %= mod;
  while (len(dat) <= n) {
    int k = len(dat);
    int q = (mod + k - 1) / k;
    dat.eb(dat[k * q - mod] * mint::raw(q));
  }
  return dat[n];
}

template <typename mint>
mint fact(int n) {
  static const int mod = mint::get_mod();
  assert(0 <= n && n < mod);
  static vector<mint> dat = {1, 1};
  while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * mint::raw(len(dat)));
  return dat[n];
}

template <typename mint>
mint fact_inv(int n) {
  static vector<mint> dat = {1, 1};
  if (n < 0) return mint(0);
  while (len(dat) <= n) dat.eb(dat[len(dat) - 1] * inv<mint>(len(dat)));
  return dat[n];
}

template <class mint, class... Ts>
mint fact_invs(Ts... xs) {
  return (mint(1) * ... * fact_inv<mint>(xs));
}

template <typename mint, class Head, class... Tail>
mint multinomial(Head &&head, Tail &&... tail) {
  return fact<mint>(head) * fact_invs<mint>(std::forward<Tail>(tail)...);
}

template <typename mint>
mint C_dense(int n, int k) {
  assert(n >= 0);
  if (k < 0 || n < k) return 0;
  static vvc<mint> C;
  static int H = 0, W = 0;
  auto calc = [&](int i, int j) -> mint {
    if (i == 0) return (j == 0 ? mint(1) : mint(0));
    return C[i - 1][j] + (j ? C[i - 1][j - 1] : 0);
  };
  if (W <= k) {
    FOR(i, H) {
      C[i].resize(k + 1);
      FOR(j, W, k + 1) { C[i][j] = calc(i, j); }
    }
    W = k + 1;
  }
  if (H <= n) {
    C.resize(n + 1);
    FOR(i, H, n + 1) {
      C[i].resize(W);
      FOR(j, W) { C[i][j] = calc(i, j); }
    }
    H = n + 1;
  }
  return C[n][k];
}

template <typename mint, bool large = false, bool dense = false>
mint C(ll n, ll k) {
  assert(n >= 0);
  if (k < 0 || n < k) return 0;
  if constexpr (dense) return C_dense<mint>(n, k);
  if constexpr (!large) return multinomial<mint>(n, k, n - k);
  k = min(k, n - k);
  mint x(1);
  FOR(i, k) x *= mint(n - i);
  return x * fact_inv<mint>(k);
}

template <typename mint, bool large = false>
mint C_inv(ll n, ll k) {
  assert(n >= 0);
  assert(0 <= k && k <= n);
  if (!large) return fact_inv<mint>(n) * fact<mint>(k) * fact<mint>(n - k);
  return mint(1) / C<mint, 1>(n, k);
}

// [x^d](1-x)^{-n}
template <typename mint, bool large = false, bool dense = false>
mint C_negative(ll n, ll d) {
  assert(n >= 0);
  if (d < 0) return mint(0);
  if (n == 0) { return (d == 0 ? mint(1) : mint(0)); }
  return C<mint, large, dense>(n + d - 1, d);
}
#line 3 "mod/modint.hpp"

template <int mod>
struct modint {
  static constexpr u32 umod = u32(mod);
  static_assert(umod < u32(1) << 31);
  u32 val;

  static modint raw(u32 v) {
    modint x;
    x.val = v;
    return x;
  }
  constexpr modint() : val(0) {}
  constexpr modint(u32 x) : val(x % umod) {}
  constexpr modint(u64 x) : val(x % umod) {}
  constexpr modint(u128 x) : val(x % umod) {}
  constexpr modint(int x) : val((x %= mod) < 0 ? x + mod : x){};
  constexpr modint(ll x) : val((x %= mod) < 0 ? x + mod : x){};
  constexpr modint(i128 x) : val((x %= mod) < 0 ? x + mod : x){};
  bool operator<(const modint &other) const { return val < other.val; }
  modint &operator+=(const modint &p) {
    if ((val += p.val) >= umod) val -= umod;
    return *this;
  }
  modint &operator-=(const modint &p) {
    if ((val += umod - p.val) >= umod) val -= umod;
    return *this;
  }
  modint &operator*=(const modint &p) {
    val = u64(val) * p.val % umod;
    return *this;
  }
  modint &operator/=(const modint &p) {
    *this *= p.inverse();
    return *this;
  }
  modint operator-() const { return modint::raw(val ? mod - val : u32(0)); }
  modint operator+(const modint &p) const { return modint(*this) += p; }
  modint operator-(const modint &p) const { return modint(*this) -= p; }
  modint operator*(const modint &p) const { return modint(*this) *= p; }
  modint operator/(const modint &p) const { return modint(*this) /= p; }
  bool operator==(const modint &p) const { return val == p.val; }
  bool operator!=(const modint &p) const { return val != p.val; }
  modint inverse() const {
    int a = val, b = mod, u = 1, v = 0, t;
    while (b > 0) {
      t = a / b;
      swap(a -= t * b, b), swap(u -= t * v, v);
    }
    return modint(u);
  }
  modint pow(ll n) const {
    if (n < 0) return inverse().pow(-n);
    assert(n >= 0);
    modint ret(1), mul(val);
    while (n > 0) {
      if (n & 1) ret *= mul;
      mul *= mul;
      n >>= 1;
    }
    return ret;
  }
  static constexpr int get_mod() { return mod; }
  // (n, r), r は 1 の 2^n 乗根
  static constexpr pair<int, int> ntt_info() {
    if (mod == 120586241) return {20, 74066978};
    if (mod == 167772161) return {25, 17};
    if (mod == 469762049) return {26, 30};
    if (mod == 754974721) return {24, 362};
    if (mod == 880803841) return {23, 211};
    if (mod == 943718401) return {22, 663003469};
    if (mod == 998244353) return {23, 31};
    if (mod == 1004535809) return {21, 582313106};
    if (mod == 1012924417) return {21, 368093570};
    return {-1, -1};
  }
  static constexpr bool can_ntt() { return ntt_info().fi != -1; }
};

#ifdef FASTIO
template <int mod>
void rd(modint<mod> &x) {
  fastio::rd(x.val);
  x.val %= mod;
  // assert(0 <= x.val && x.val < mod);
}
template <int mod>
void wt(modint<mod> x) {
  fastio::wt(x.val);
}
#endif

using modint107 = modint<1000000007>;
using modint998 = modint<998244353>;
#line 2 "mod/mod_inv.hpp"

// long でも大丈夫

// (val * x - 1) が mod の倍数になるようにする

// 特に mod=0 なら x=0 が満たす

ll mod_inv(ll val, ll mod) {
  if (mod == 0) return 0;
  mod = abs(mod);
  val %= mod;
  if (val < 0) val += mod;
  ll a = val, b = mod, u = 1, v = 0, t;
  while (b > 0) {
    t = a / b;
    swap(a -= t * b, b), swap(u -= t * v, v);
  }
  if (u < 0) u += mod;
  return u;
}
#line 2 "mod/crt3.hpp"

constexpr u32 mod_pow_constexpr(u64 a, u64 n, u32 mod) {
  a %= mod;
  u64 res = 1;
  FOR(32) {
    if (n & 1) res = res * a % mod;
    a = a * a % mod, n /= 2;
  }
  return res;
}

template <typename T, u32 p0, u32 p1>
T CRT2(u64 a0, u64 a1) {
  static_assert(p0 < p1);
  static constexpr u64 x0_1 = mod_pow_constexpr(p0, p1 - 2, p1);
  u64 c = (a1 - a0 + p1) * x0_1 % p1;
  return a0 + c * p0;
}

template <typename T, u32 p0, u32 p1, u32 p2>
T CRT3(u64 a0, u64 a1, u64 a2) {
  static_assert(p0 < p1 && p1 < p2);
  static constexpr u64 x1 = mod_pow_constexpr(p0, p1 - 2, p1);
  static constexpr u64 x2 = mod_pow_constexpr(u64(p0) * p1 % p2, p2 - 2, p2);
  static constexpr u64 p01 = u64(p0) * p1;
  u64 c = (a1 - a0 + p1) * x1 % p1;
  u64 ans_1 = a0 + c * p0;
  c = (a2 - ans_1 % p2 + p2) * x2 % p2;
  return T(ans_1) + T(c) * T(p01);
}

template <typename T, u32 p0, u32 p1, u32 p2, u32 p3>
T CRT4(u64 a0, u64 a1, u64 a2, u64 a3) {
  static_assert(p0 < p1 && p1 < p2 && p2 < p3);
  static constexpr u64 x1 = mod_pow_constexpr(p0, p1 - 2, p1);
  static constexpr u64 x2 = mod_pow_constexpr(u64(p0) * p1 % p2, p2 - 2, p2);
  static constexpr u64 x3 = mod_pow_constexpr(u64(p0) * p1 % p3 * p2 % p3, p3 - 2, p3);
  static constexpr u64 p01 = u64(p0) * p1;
  u64 c = (a1 - a0 + p1) * x1 % p1;
  u64 ans_1 = a0 + c * p0;
  c = (a2 - ans_1 % p2 + p2) * x2 % p2;
  u128 ans_2 = ans_1 + c * static_cast<u128>(p01);
  c = (a3 - ans_2 % p3 + p3) * x3 % p3;
  return T(ans_2) + T(c) * T(p01) * T(p2);
}

template <typename T, u32 p0, u32 p1, u32 p2, u32 p3, u32 p4>
T CRT5(u64 a0, u64 a1, u64 a2, u64 a3, u64 a4) {
  static_assert(p0 < p1 && p1 < p2 && p2 < p3 && p3 < p4);
  static constexpr u64 x1 = mod_pow_constexpr(p0, p1 - 2, p1);
  static constexpr u64 x2 = mod_pow_constexpr(u64(p0) * p1 % p2, p2 - 2, p2);
  static constexpr u64 x3 = mod_pow_constexpr(u64(p0) * p1 % p3 * p2 % p3, p3 - 2, p3);
  static constexpr u64 x4 = mod_pow_constexpr(u64(p0) * p1 % p4 * p2 % p4 * p3 % p4, p4 - 2, p4);
  static constexpr u64 p01 = u64(p0) * p1;
  static constexpr u64 p23 = u64(p2) * p3;
  u64 c = (a1 - a0 + p1) * x1 % p1;
  u64 ans_1 = a0 + c * p0;
  c = (a2 - ans_1 % p2 + p2) * x2 % p2;
  u128 ans_2 = ans_1 + c * static_cast<u128>(p01);
  c = static_cast<u64>(a3 - ans_2 % p3 + p3) * x3 % p3;
  u128 ans_3 = ans_2 + static_cast<u128>(c * p2) * p01;
  c = static_cast<u64>(a4 - ans_3 % p4 + p4) * x4 % p4;
  return T(ans_3) + T(c) * T(p01) * T(p23);
}
#line 2 "poly/convolution_naive.hpp"

template <class T, typename enable_if<!has_mod<T>::value>::type* = nullptr>
vc<T> convolution_naive(const vc<T>& a, const vc<T>& b) {
  int n = int(a.size()), m = int(b.size());
  if (n > m) return convolution_naive<T>(b, a);
  if (n == 0) return {};
  vector<T> ans(n + m - 1);
  FOR(i, n) FOR(j, m) ans[i + j] += a[i] * b[j];
  return ans;
}

template <class T, typename enable_if<has_mod<T>::value>::type* = nullptr>
vc<T> convolution_naive(const vc<T>& a, const vc<T>& b) {
  int n = int(a.size()), m = int(b.size());
  if (n > m) return convolution_naive<T>(b, a);
  if (n == 0) return {};
  vc<T> ans(n + m - 1);
  if (n <= 16 && (T::get_mod() < (1 << 30))) {
    for (int k = 0; k < n + m - 1; ++k) {
      int s = max(0, k - m + 1);
      int t = min(n, k + 1);
      u64 sm = 0;
      for (int i = s; i < t; ++i) { sm += u64(a[i].val) * (b[k - i].val); }
      ans[k] = sm;
    }
  } else {
    for (int k = 0; k < n + m - 1; ++k) {
      int s = max(0, k - m + 1);
      int t = min(n, k + 1);
      u128 sm = 0;
      for (int i = s; i < t; ++i) { sm += u64(a[i].val) * (b[k - i].val); }
      ans[k] = T::raw(sm % T::get_mod());
    }
  }
  return ans;
}
#line 2 "poly/convolution_karatsuba.hpp"

// 任意の環でできる
template <typename T>
vc<T> convolution_karatsuba(const vc<T>& f, const vc<T>& g) {
  const int thresh = 30;
  if (min(len(f), len(g)) <= thresh) return convolution_naive(f, g);
  int n = max(len(f), len(g));
  int m = ceil(n, 2);
  vc<T> f1, f2, g1, g2;
  if (len(f) < m) f1 = f;
  if (len(f) >= m) f1 = {f.begin(), f.begin() + m};
  if (len(f) >= m) f2 = {f.begin() + m, f.end()};
  if (len(g) < m) g1 = g;
  if (len(g) >= m) g1 = {g.begin(), g.begin() + m};
  if (len(g) >= m) g2 = {g.begin() + m, g.end()};
  vc<T> a = convolution_karatsuba(f1, g1);
  vc<T> b = convolution_karatsuba(f2, g2);
  FOR(i, len(f2)) f1[i] += f2[i];
  FOR(i, len(g2)) g1[i] += g2[i];
  vc<T> c = convolution_karatsuba(f1, g1);
  vc<T> F(len(f) + len(g) - 1);
  assert(2 * m + len(b) <= len(F));
  FOR(i, len(a)) F[i] += a[i], c[i] -= a[i];
  FOR(i, len(b)) F[2 * m + i] += b[i], c[i] -= b[i];
  if (c.back() == T(0)) c.pop_back();
  FOR(i, len(c)) if (c[i] != T(0)) F[m + i] += c[i];
  return F;
}
#line 2 "poly/ntt.hpp"

template <class mint>
void ntt(vector<mint>& a, bool inverse) {
  assert(mint::can_ntt());
  const int rank2 = mint::ntt_info().fi;
  const int mod = mint::get_mod();
  static array<mint, 30> root, iroot;
  static array<mint, 30> rate2, irate2;
  static array<mint, 30> rate3, irate3;

  assert(rank2 != -1 && len(a) <= (1 << max(0, rank2)));

  static bool prepared = 0;
  if (!prepared) {
    prepared = 1;
    root[rank2] = mint::ntt_info().se;
    iroot[rank2] = mint(1) / root[rank2];
    FOR_R(i, rank2) {
      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];
    }
    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];
    }
  }

  int n = int(a.size());
  int h = topbit(n);
  assert(n == 1 << h);
  if (!inverse) {
    int len = 0;
    while (len < h) {
      if (h - len == 1) {
        int p = 1 << (h - len - 1);
        mint rot = 1;
        FOR(s, 1 << len) {
          int offset = s << (h - len);
          FOR(i, p) {
            auto l = a[i + offset];
            auto r = a[i + offset + p] * rot;
            a[i + offset] = l + r;
            a[i + offset + p] = l - r;
          }
          rot *= rate2[topbit(~s & -~s)];
        }
        len++;
      } else {
        int p = 1 << (h - len - 2);
        mint rot = 1, imag = 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++) {
            u64 mod2 = u64(mod) * mod;
            u64 a0 = a[i + offset].val;
            u64 a1 = u64(a[i + offset + p].val) * rot.val;
            u64 a2 = u64(a[i + offset + 2 * p].val) * rot2.val;
            u64 a3 = u64(a[i + offset + 3 * p].val) * rot3.val;
            u64 a1na3imag = (a1 + mod2 - a3) % mod * imag.val;
            u64 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);
          }
          rot *= rate3[topbit(~s & -~s)];
        }
        len += 2;
      }
    }
  } else {
    mint coef = mint(1) / mint(len(a));
    FOR(i, len(a)) a[i] *= coef;
    int len = h;
    while (len) {
      if (len == 1) {
        int p = 1 << (h - len);
        mint irot = 1;
        FOR(s, 1 << (len - 1)) {
          int offset = s << (h - len + 1);
          FOR(i, p) {
            u64 l = a[i + offset].val;
            u64 r = a[i + offset + p].val;
            a[i + offset] = l + r;
            a[i + offset + p] = (mod + l - r) * irot.val;
          }
          irot *= irate2[topbit(~s & -~s)];
        }
        len--;
      } else {
        int p = 1 << (h - len);
        mint irot = 1, iimag = iroot[2];
        FOR(s, (1 << (len - 2))) {
          mint irot2 = irot * irot;
          mint irot3 = irot2 * irot;
          int offset = s << (h - len + 2);
          for (int i = 0; i < p; i++) {
            u64 a0 = a[i + offset + 0 * p].val;
            u64 a1 = a[i + offset + 1 * p].val;
            u64 a2 = a[i + offset + 2 * p].val;
            u64 a3 = a[i + offset + 3 * p].val;
            u64 x = (mod + a2 - a3) * iimag.val % mod;
            a[i + offset] = a0 + a1 + a2 + a3;
            a[i + offset + 1 * p] = (a0 + mod - a1 + x) * irot.val;
            a[i + offset + 2 * p] = (a0 + a1 + 2 * mod - a2 - a3) * irot2.val;
            a[i + offset + 3 * p] = (a0 + 2 * mod - a1 - x) * irot3.val;
          }
          irot *= irate3[topbit(~s & -~s)];
        }
        len -= 2;
      }
    }
  }
}
#line 8 "poly/convolution.hpp"

template <class mint>
vector<mint> convolution_ntt(vector<mint> a, vector<mint> b) {
  if (a.empty() || b.empty()) return {};
  int n = int(a.size()), m = int(b.size());
  int sz = 1;
  while (sz < n + m - 1) sz *= 2;

  // sz = 2^k のときの高速化。分割統治的なやつで損しまくるので。

  if ((n + m - 3) <= sz / 2) {
    auto a_last = a.back(), b_last = b.back();
    a.pop_back(), b.pop_back();
    auto c = convolution(a, b);
    c.resize(n + m - 1);
    c[n + m - 2] = a_last * b_last;
    FOR(i, len(a)) c[i + len(b)] += a[i] * b_last;
    FOR(i, len(b)) c[i + len(a)] += b[i] * a_last;
    return c;
  }

  a.resize(sz), b.resize(sz);
  bool same = a == b;
  ntt(a, 0);
  if (same) {
    b = a;
  } else {
    ntt(b, 0);
  }
  FOR(i, sz) a[i] *= b[i];
  ntt(a, 1);
  a.resize(n + m - 1);
  return a;
}

template <typename mint>
vector<mint> convolution_garner(const vector<mint>& a, const vector<mint>& b) {
  int n = len(a), m = len(b);
  if (!n || !m) return {};
  static constexpr int p0 = 167772161;
  static constexpr int p1 = 469762049;
  static constexpr int p2 = 754974721;
  using mint0 = modint<p0>;
  using mint1 = modint<p1>;
  using mint2 = modint<p2>;
  vc<mint0> a0(n), b0(m);
  vc<mint1> a1(n), b1(m);
  vc<mint2> a2(n), b2(m);
  FOR(i, n) a0[i] = a[i].val, a1[i] = a[i].val, a2[i] = a[i].val;
  FOR(i, m) b0[i] = b[i].val, b1[i] = b[i].val, b2[i] = b[i].val;
  auto c0 = convolution_ntt<mint0>(a0, b0);
  auto c1 = convolution_ntt<mint1>(a1, b1);
  auto c2 = convolution_ntt<mint2>(a2, b2);
  vc<mint> c(len(c0));
  FOR(i, n + m - 1) { c[i] = CRT3<mint, p0, p1, p2>(c0[i].val, c1[i].val, c2[i].val); }
  return c;
}

vector<ll> convolution(vector<ll> a, vector<ll> b) {
  int n = len(a), m = len(b);
  if (!n || !m) return {};
  if (min(n, m) <= 2500) return convolution_naive(a, b);

  ll mi_a = MIN(a), mi_b = MIN(b);
  for (auto& x: a) x -= mi_a;
  for (auto& x: b) x -= mi_b;
  assert(MAX(a) * MAX(b) <= 1e18);

  auto Ac = cumsum<ll>(a), Bc = cumsum<ll>(b);
  vi res(n + m - 1);
  for (int k = 0; k < n + m - 1; ++k) {
    int s = max(0, k - m + 1);
    int t = min(n, k + 1);
    res[k] += (t - s) * mi_a * mi_b;
    res[k] += mi_a * (Bc[k - s + 1] - Bc[k - t + 1]);
    res[k] += mi_b * (Ac[t] - Ac[s]);
  }

  static constexpr u32 MOD1 = 1004535809;
  static constexpr u32 MOD2 = 1012924417;
  using mint1 = modint<MOD1>;
  using mint2 = modint<MOD2>;

  vc<mint1> a1(n), b1(m);
  vc<mint2> a2(n), b2(m);
  FOR(i, n) a1[i] = a[i], a2[i] = a[i];
  FOR(i, m) b1[i] = b[i], b2[i] = b[i];

  auto c1 = convolution_ntt<mint1>(a1, b1);
  auto c2 = convolution_ntt<mint2>(a2, b2);

  FOR(i, n + m - 1) { res[i] += CRT2<u64, MOD1, MOD2>(c1[i].val, c2[i].val); }
  return res;
}

template <typename mint>
vc<mint> convolution(const vc<mint>& a, const vc<mint>& b) {
  int n = len(a), m = len(b);
  if (!n || !m) return {};
  if (mint::can_ntt()) {
    if (min(n, m) <= 50) return convolution_karatsuba<mint>(a, b);
    return convolution_ntt(a, b);
  }
  if (min(n, m) <= 200) return convolution_karatsuba<mint>(a, b);
  return convolution_garner(a, b);
}
#line 4 "poly/fps_inv.hpp"

template <typename mint>
vc<mint> fps_inv_sparse(const vc<mint>& f) {
  int N = len(f);
  vc<pair<int, mint>> dat;
  FOR(i, 1, N) if (f[i] != mint(0)) dat.eb(i, f[i]);
  vc<mint> g(N);
  mint g0 = mint(1) / f[0];
  g[0] = g0;
  FOR(n, 1, N) {
    mint rhs = 0;
    for (auto&& [k, fk]: dat) {
      if (k > n) break;
      rhs -= fk * g[n - k];
    }
    g[n] = rhs * g0;
  }
  return g;
}

template <typename mint>
vc<mint> fps_inv_dense_ntt(const vc<mint>& F) {
  vc<mint> G = {mint(1) / F[0]};
  ll N = len(F), n = 1;
  G.reserve(N);
  while (n < N) {
    vc<mint> f(2 * n), g(2 * n);
    FOR(i, min(N, 2 * n)) f[i] = F[i];
    FOR(i, n) g[i] = G[i];
    ntt(f, false), ntt(g, false);
    FOR(i, 2 * n) f[i] *= g[i];
    ntt(f, true);
    FOR(i, n) f[i] = 0;
    ntt(f, false);
    FOR(i, 2 * n) f[i] *= g[i];
    ntt(f, true);
    FOR(i, n, min(N, 2 * n)) G.eb(-f[i]);
    n *= 2;
  }
  return G;
}

template <typename mint>
vc<mint> fps_inv_dense(const vc<mint>& F) {
  if (mint::can_ntt()) return fps_inv_dense_ntt(F);
  const int N = len(F);
  vc<mint> R = {mint(1) / F[0]};
  vc<mint> p;
  int m = 1;
  while (m < N) {
    p = convolution(R, R);
    p.resize(m + m);
    vc<mint> f = {F.begin(), F.begin() + min(m + m, N)};
    p = convolution(p, f);
    R.resize(m + m);
    FOR(i, m + m) R[i] = R[i] + R[i] - p[i];
    m += m;
  }
  R.resize(N);
  return R;
}

template <typename mint>
vc<mint> fps_inv(const vc<mint>& f) {
  assert(f[0] != mint(0));
  int n = count_terms(f);
  int t = (mint::can_ntt() ? 160 : 820);
  return (n <= t ? fps_inv_sparse<mint>(f) : fps_inv_dense<mint>(f));
}
#line 2 "seq/geometric_sequence_sum_formula.hpp"

// r != 1 とする. sum_{i=0}^{n-1}i^kr^i = c + r^n f(n) と書ける. return : c, f
template <typename mint>
pair<mint, vc<mint>> geometric_sequence_sum_formula(mint r, int K) {
  assert(r != 1);
  vc<mint> F(K + 1);
  F[0] = 1;
  FOR(i, K + 1) F[i] -= r * fact_inv<mint>(i);
  F = fps_inv<mint>(F);
  mint c = F[K] * fact<mint>(K);
  reverse(all(F));
  FOR(k, K + 1) F[k] *= -fact_inv<mint>(k) * fact<mint>(K);
  return {c, F};
}
#line 2 "convex/line_min_function.hpp"

#line 2 "geo/convex_hull.hpp"

#line 2 "geo/base.hpp"
template <typename T>
struct Point {
  T x, y;

  Point() : x(0), y(0) {}

  template <typename A, typename B>
  Point(A x, B y) : x(x), y(y) {}

  template <typename A, typename B>
  Point(pair<A, B> p) : x(p.fi), y(p.se) {}

  Point operator+=(const Point p) {
    x += p.x, y += p.y;
    return *this;
  }
  Point operator-=(const Point p) {
    x -= p.x, y -= p.y;
    return *this;
  }
  Point operator+(Point p) const { return {x + p.x, y + p.y}; }
  Point operator-(Point p) const { return {x - p.x, y - p.y}; }
  bool operator==(Point p) const { return x == p.x && y == p.y; }
  bool operator!=(Point p) const { return x != p.x || y != p.y; }
  Point operator-() const { return {-x, -y}; }
  Point operator*(T t) const { return {x * t, y * t}; }
  Point operator/(T t) const { return {x / t, y / t}; }

  bool operator<(Point p) const {
    if (x != p.x) return x < p.x;
    return y < p.y;
  }
  T dot(const Point& other) const { return x * other.x + y * other.y; }
  T det(const Point& other) const { return x * other.y - y * other.x; }

  double norm() { return sqrtl(x * x + y * y); }
  double angle() { return atan2(y, x); }

  Point rotate(double theta) {
    static_assert(!is_integral<T>::value);
    double c = cos(theta), s = sin(theta);
    return Point{c * x - s * y, s * x + c * y};
  }
  Point rot90(bool ccw) { return (ccw ? Point{-y, x} : Point{y, -x}); }
};

#ifdef FASTIO
template <typename T>
void rd(Point<T>& p) {
  fastio::rd(p.x), fastio::rd(p.y);
}
template <typename T>
void wt(Point<T>& p) {
  fastio::wt(p.x);
  fastio::wt(' ');
  fastio::wt(p.y);
}
#endif

// A -> B -> C と進むときに、左に曲がるならば +1、右に曲がるならば -1
template <typename T>
int ccw(Point<T> A, Point<T> B, Point<T> C) {
  T x = (B - A).det(C - A);
  if (x > 0) return 1;
  if (x < 0) return -1;
  return 0;
}

template <typename REAL, typename T, typename U>
REAL dist(Point<T> A, Point<U> B) {
  REAL dx = REAL(A.x) - REAL(B.x);
  REAL dy = REAL(A.y) - REAL(B.y);
  return sqrt(dx * dx + dy * dy);
}

// ax+by+c
template <typename T>
struct Line {
  T a, b, c;

  Line(T a, T b, T c) : a(a), b(b), c(c) {}
  Line(Point<T> A, Point<T> B) { a = A.y - B.y, b = B.x - A.x, c = A.x * B.y - A.y * B.x; }
  Line(T x1, T y1, T x2, T y2) : Line(Point<T>(x1, y1), Point<T>(x2, y2)) {}

  template <typename U>
  U eval(Point<U> P) {
    return a * P.x + b * P.y + c;
  }

  template <typename U>
  T eval(U x, U y) {
    return a * x + b * y + c;
  }

  // 同じ直線が同じ a,b,c で表現されるようにする
  void normalize() {
    static_assert(is_same_v<T, int> || is_same_v<T, long long>);
    T g = gcd(gcd(abs(a), abs(b)), abs(c));
    a /= g, b /= g, c /= g;
    if (b < 0) { a = -a, b = -b, c = -c; }
    if (b == 0 && a < 0) { a = -a, b = -b, c = -c; }
  }

  bool is_parallel(Line other) { return a * other.b - b * other.a == 0; }
  bool is_orthogonal(Line other) { return a * other.a + b * other.b == 0; }
};

template <typename T>
struct Segment {
  Point<T> A, B;

  Segment(Point<T> A, Point<T> B) : A(A), B(B) {}
  Segment(T x1, T y1, T x2, T y2) : Segment(Point<T>(x1, y1), Point<T>(x2, y2)) {}

  bool contain(Point<T> C) {
    T det = (C - A).det(B - A);
    if (det != 0) return 0;
    return (C - A).dot(B - A) >= 0 && (C - B).dot(A - B) >= 0;
  }

  Line<T> to_Line() { return Line(A, B); }
};

template <typename REAL>
struct Circle {
  Point<REAL> O;
  REAL r;
  Circle(Point<REAL> O, REAL r) : O(O), r(r) {}
  Circle(REAL x, REAL y, REAL r) : O(x, y), r(r) {}
  template <typename T>
  bool contain(Point<T> p) {
    REAL dx = p.x - O.x, dy = p.y - O.y;
    return dx * dx + dy * dy <= r * r;
  }
};
#line 4 "geo/convex_hull.hpp"

// allow_180=true で同一座標点があるとこわれる
// full なら I[0] が sorted で min になる
template <typename T, bool allow_180 = false>
vector<int> ConvexHull(vector<Point<T>>& XY, string mode = "full", bool sorted = false) {
  assert(mode == "full" || mode == "lower" || mode == "upper");
  ll N = XY.size();
  if (N == 1) return {0};
  if (N == 2) {
    if (XY[0] < XY[1]) return {0, 1};
    if (XY[1] < XY[0]) return {1, 0};
    return {0};
  }
  vc<int> I(N);
  if (sorted) {
    FOR(i, N) I[i] = i;
  } else {
    I = argsort(XY);
  }
  if constexpr (allow_180) { FOR(i, N - 1) assert(XY[i] != XY[i + 1]); }

  auto check = [&](ll i, ll j, ll k) -> bool {
    T det = (XY[j] - XY[i]).det(XY[k] - XY[i]);
    if constexpr (allow_180) return det >= 0;
    return det > T(0);
  };

  auto calc = [&]() {
    vector<int> P;
    for (auto&& k: I) {
      while (P.size() > 1) {
        auto i = P[P.size() - 2];
        auto j = P[P.size() - 1];
        if (check(i, j, k)) break;
        P.pop_back();
      }
      P.eb(k);
    }
    return P;
  };

  vc<int> P;
  if (mode == "full" || mode == "lower") {
    vc<int> Q = calc();
    P.insert(P.end(), all(Q));
  }
  if (mode == "full" || mode == "upper") {
    if (!P.empty()) P.pop_back();
    reverse(all(I));
    vc<int> Q = calc();
    P.insert(P.end(), all(Q));
  }
  if (mode == "upper") reverse(all(P));
  while (len(P) >= 2 && XY[P[0]] == XY[P.back()]) P.pop_back();
  return P;
}
#line 4 "convex/line_min_function.hpp"

// 1 次関数の max を [L,R,a,b] の列として出力
// https://qoj.ac/contest/1576/problem/8505
template <typename Re, typename T>
vc<tuple<Re, Re, Re, Re>> line_min_function_real(vc<pair<T, T>> LINE) {
  assert(!LINE.empty());
  using P = Point<T>;
  vc<P> point;
  for (auto& [x, y]: LINE) point.eb(P(x, y));
  auto I = ConvexHull(point, "lower");
  point = rearrange(point, I);
  int N = len(point);
  if (N >= 2 && point[N - 1].x == point[N - 2].x) { POP(point), --N; }
  reverse(all(point)); // 傾きは大きい方から
  Re l = -infty<Re>;
  vc<tuple<Re, Re, Re, Re>> ANS;
  FOR(i, N) {
    Re r = infty<Re>;
    auto [a, b] = point[i];
    if (i + 1 < N) {
      auto [c, d] = point[i + 1];
      if (a == c) continue;
      assert(a > c);
      r = Re(d - b) / (a - c);
      chmax(r, l), chmin(r, infty<Re>);
    }
    if (l < r) ANS.eb(l, r, a, b), l = r;
  }
  return ANS;
}

// 1 次関数の max を [L,R,a,b] の列として出力
template <typename Re, typename T>
vc<tuple<Re, Re, Re, Re>> line_max_function_real(vc<pair<T, T>> LINE) {
  assert(!LINE.empty());
  for (auto& [a, b]: LINE) a = -a, b = -b;
  auto ANS = line_min_function_real<Re, T>(LINE);
  for (auto& [l, r, a, b]: ANS) a = -a, b = -b;
  return ANS;
}

// LINE(a,b,c): y=(ax+b)/c, 評価点は整数
// 1 次関数の min を [L,R,a,b,c] の列として出力
// オーバーフロー安全
vc<tuple<ll, ll, ll, ll, ll>> line_min_function_rational(vc<tuple<ll, ll, ll>> LINE, ll L, ll R) {
  // 傾き降順
  sort(all(LINE), [&](auto& L, auto& R) -> bool {
    auto& [a1, b1, c1] = L;
    auto& [a2, b2, c2] = R;
    return i128(a1) * c2 > i128(a2) * c1;
  });
  vc<tuple<ll, ll, ll, ll, ll>> ANS;
  for (auto& [a2, b2, c2]: LINE) {
    while (1) {
      if (ANS.empty()) {
        ANS.eb(L, R, a2, b2, c2);
        break;
      }
      auto& [L1, R1, a1, b1, c1] = ANS.back();
      i128 s = i128(c2) * a1 - i128(a2) * c1; // >= 0
      i128 t = i128(b2) * c1 - i128(b1) * c2;
      if (s == 0) {
        // 平行なので小さい方だけを残す
        if (t >= 0) break;
        ANS.pop_back();
        if (len(ANS)) get<1>(ANS.back()) = R;
        continue;
      }
      i128 x = ceil<i128>(t, s);
      // x 以上で 2 の方が下に来る
      if (x <= L1) {
        ANS.pop_back();
        continue;
      }
      if (x < R) {
        R1 = x;
        ANS.eb(x, R, a2, b2, c2);
        break;
      } else {
        break;
      }
    }
  }
  return ANS;
}

// LINE(a,b,c): y=(ax+b)/c, 評価点は整数
// 1 次関数の max を [L,R,a,b,c] の列として出力
// オーバーフロー安全
vc<tuple<ll, ll, ll, ll, ll>> line_max_function_rational(vc<tuple<ll, ll, ll>> LINE, ll L, ll R) {
  for (auto& [a, b, c]: LINE) a = -a, b = -b;
  auto ANS = line_min_function_rational(LINE, L, R);
  for (auto& [L, R, a, b, c]: ANS) a = -a, b = -b;
  return ANS;
}

// LINE(a,b): y=ax+b, 評価点は整数
// 1 次関数の min を [L,R,a,b] の列として出力
// オーバーフロー安全
vc<tuple<ll, ll, ll, ll>> line_min_function_integer(vc<pair<ll, ll>> LINE, ll L, ll R) {
  // 傾き降順
  sort(all(LINE), [&](auto& L, auto& R) -> bool {
    auto& [a1, b1] = L;
    auto& [a2, b2] = R;
    return a1 > a2;
  });
  vc<tuple<ll, ll, ll, ll>> ANS;
  for (auto& [a2, b2]: LINE) {
    while (1) {
      if (ANS.empty()) {
        ANS.eb(L, R, a2, b2);
        break;
      }
      auto& [L1, R1, a1, b1] = ANS.back();
      if (a1 == a2) {
        if (b1 <= b2) break;
        ANS.pop_back();
        if (len(ANS)) get<1>(ANS.back()) = R;
        continue;
      }
      ll x = ceil<ll>(b2 - b1, a1 - a2);
      // x 以上で 2 の方が下に来る
      if (x <= L1) {
        ANS.pop_back();
        continue;
      }
      if (x < R) {
        R1 = x;
        ANS.eb(x, R, a2, b2);
        break;
      } else {
        break;
      }
    }
  }
  return ANS;
}

// LINE(a,b,c): y=(ax+b)/c, 評価点は整数
// 1 次関数の min を [L,R,a,b,c] の列として出力
// c>0, (ax+b)c がオーバーフローしない,
vc<tuple<ll, ll, ll, ll>> line_max_function_integer(vc<pair<ll, ll>> LINE, ll L, ll R) {
  for (auto& [a, b]: LINE) a = -a, b = -b;
  auto ANS = line_min_function_integer(LINE, L, R);
  for (auto& [L, R, a, b]: ANS) a = -a, b = -b;
  return ANS;
}

// (L,R,func) の下側と上側をマージするときなどに使う用
template <typename T>
vc<tuple<T, T, T, T, T, T>> merge_46(vc<tuple<T, T, T, T>> A, vc<tuple<T, T, T, T>> B) {
  vc<tuple<T, T, T, T, T, T>> ANS;
  reverse(all(A));
  reverse(all(B));
  while (len(A) && len(B)) {
    auto& [l1, r1, a1, b1] = A.back();
    auto& [l2, r2, a2, b2] = B.back();
    assert(l1 == l2);
    T r = min(r1, r2);
    ANS.eb(l1, r, a1, b1, a2, b2);
    l1 = r, l2 = r;
    if (r1 == r) POP(A);
    if (r2 == r) POP(B);
  };
  return ANS;
}

// (L,R,func) の下側と上側をマージするときなどに使う用
// f(L,R,a1,b1,a2,b2)
template <typename T, typename F>
void merge_46(const vc<tuple<T, T, T, T>>& A, const vc<tuple<T, T, T, T>>& B, F f) {
  int i = 0, j = 0;
  while (i < len(A) && j < len(B)) {
    auto& [l1, r1, a1, b1] = A[i];
    auto& [l2, r2, a2, b2] = B[j];
    T l = max(l1, l2), r = min(r1, r2);
    if (l < r) f(l, r, a1, b1, a2, b2);
    (r1 < r2 ? i : j)++;
  }
}

// (L,R,func) の下側と上側をマージするときなどに使う用
// f(L,R,a1,b1,a2,b2)
template <typename T, typename F>
void merge_58(const vc<tuple<T, T, T, T, T>>& A, const vc<tuple<T, T, T, T, T>>& B, F f) {
  int i = 0, j = 0;
  while (i < len(A) && j < len(B)) {
    auto& [l1, r1, a1, b1, c1] = A[i];
    auto& [l2, r2, a2, b2, c2] = B[j];
    T l = max(l1, l2), r = min(r1, r2);
    if (l < r) f(l, r, a1, b1, c1, a2, b2, c2);
    (r1 < r2 ? i : j)++;
  }
}
#line 2 "mod/floor_sum_of_linear_polynomial.hpp"

#line 2 "alg/monoid_pow.hpp"

// chat gpt
template <typename U, typename Arg1, typename Arg2>
struct has_power_method {
private:
  // ヘルパー関数の実装
  template <typename V, typename A1, typename A2>
  static auto check(int)
      -> decltype(std::declval<V>().power(std::declval<A1>(),
                                          std::declval<A2>()),
                  std::true_type{});
  template <typename, typename, typename>
  static auto check(...) -> std::false_type;

public:
  // メソッドの有無を表す型
  static constexpr bool value = decltype(check<U, Arg1, Arg2>(0))::value;
};

template <typename Monoid>
typename Monoid::X monoid_pow(typename Monoid::X x, ll exp) {
  using X = typename Monoid::X;
  if constexpr (has_power_method<Monoid, X, ll>::value) {
    return Monoid::power(x, exp);
  } else {
    assert(exp >= 0);
    X res = Monoid::unit();
    while (exp) {
      if (exp & 1) res = Monoid::op(res, x);
      x = Monoid::op(x, x);
      exp >>= 1;
    }
    return res;
  }
}
#line 3 "mod/floor_monoid_product.hpp"

// https://yukicoder.me/submissions/883884
// https://qoj.ac/contest/1411/problem/7620
// U は範囲内で ax+b がオーバーフローしない程度
// yyy x yyyy x ... yyy x yyy (x を N 個)
// k 個目の x までに floor(ak+b,m) 個の y がある
// my<=ax+b における lattice path における辺の列と見なせる
template <typename Monoid, typename X, typename U>
X floor_monoid_product(X x, X y, U N, U a, U b, U m) {
  U c = (a * N + b) / m;
  X pre = Monoid::unit(), suf = Monoid::unit();
  while (1) {
    const U p = a / m, q = b / m;
    a %= m, b %= m;
    x = Monoid::op(x, monoid_pow<Monoid>(y, p));
    pre = Monoid::op(pre, monoid_pow<Monoid>(y, q));
    c -= (p * N + q);
    if (c == 0) break;
    const U d = (m * c - b - 1) / a + 1;
    suf = Monoid::op(y, Monoid::op(monoid_pow<Monoid>(x, N - d), suf));
    b = m - b - 1 + a, N = c - 1, c = d;
    swap(m, a), swap(x, y);
  }
  x = monoid_pow<Monoid>(x, N);
  return Monoid::op(Monoid::op(pre, x), suf);
}
#line 2 "alg/monoid/monoid_for_floor_sum.hpp"
// sum i^k1floor^k2: floor path で (x,y) から x 方向に進むときに x^k1y^k2 を足す
template <typename T, int K1, int K2>
struct Monoid_for_floor_sum {
  using ARR = array<array<T, K2 + 1>, K1 + 1>;
  struct Data {
    ARR dp;
    T dx, dy;
  };

  using value_type = Data;
  using X = value_type;
  static X op(X a, X b) {
    static constexpr int n = max(K1, K2);
    static T comb[n + 1][n + 1];
    if (comb[0][0] != T(1)) {
      comb[0][0] = T(1);
      FOR(i, n) FOR(j, i + 1) { comb[i + 1][j] += comb[i][j], comb[i + 1][j + 1] += comb[i][j]; }
    }

    array<T, K1 + 1> pow_x;
    array<T, K2 + 1> pow_y;
    pow_x[0] = 1, pow_y[0] = 1;
    FOR(i, K1) pow_x[i + 1] = pow_x[i] * a.dx;
    FOR(i, K2) pow_y[i + 1] = pow_y[i] * a.dy;

    // +dy
    FOR(i, K1 + 1) {
      FOR_R(j, K2 + 1) {
        T x = b.dp[i][j];
        FOR(k, j + 1, K2 + 1) b.dp[i][k] += comb[k][j] * pow_y[k - j] * x;
      }
    }
    // +dx
    FOR(j, K2 + 1) {
      FOR_R(i, K1 + 1) { FOR(k, i, K1 + 1) a.dp[k][j] += comb[k][i] * pow_x[k - i] * b.dp[i][j]; }
    }

    a.dx += b.dx, a.dy += b.dy;
    return a;
  }

  static X to_x() {
    X x = unit();
    x.dp[0][0] = 1, x.dx = 1;
    return x;
  }
  static X to_y() {
    X x = unit();
    x.dy = 1;
    return x;
  }
  static constexpr X unit() { return {ARR{}, T(0), T(0)}; }
  static constexpr bool commute = 0;
};
#line 5 "mod/floor_sum_of_linear_polynomial.hpp"

// 全部非負, T は答, U は ax+b がオーバーフローしない
template <typename T, int K1, int K2, typename U>
array<array<T, K2 + 1>, K1 + 1> floor_sum_of_linear_polynomial_nonnegative(U N, U a, U b, U mod) {
  static_assert(is_same_v<U, u64> || is_same_v<U, u128>);
  assert(a == 0 || N < (U(-1) - b) / a);
  using Mono = Monoid_for_floor_sum<T, K1, K2>;
  auto x = floor_monoid_product<Mono>(Mono::to_x(), Mono::to_y(), N, a, b, mod);
  return x.dp;
};

// sum_{L<=x<R} x^i floor(ax+b,mod)^j
// a+bx が I, U でオーバーフローしない
template <typename T, int K1, int K2, typename I>
array<array<T, K2 + 1>, K1 + 1> floor_sum_of_linear_polynomial(I L, I R, I a, I b, I mod) {
  static_assert(is_same_v<I, ll> || is_same_v<I, i128>);
  assert(L <= R && mod > 0);
  if (a < 0) {
    auto ANS = floor_sum_of_linear_polynomial<T, K1, K2, I>(-R + 1, -L + 1, -a, b, mod);
    FOR(i, K1 + 1) {
      if (i % 2 == 1) { FOR(j, K2 + 1) ANS[i][j] = -ANS[i][j]; }
    }
    return ANS;
  }
  assert(a >= 0);
  I ADD_X = L;
  I N = R - L;
  b += a * L;
  I ADD_Y = floor<I>(b, mod);
  b -= ADD_Y * mod;
  assert(a >= 0 && b >= 0);

  using Mono = Monoid_for_floor_sum<T, K1, K2>;
  using Data = typename Mono::Data;
  using U = std::conditional_t<is_same_v<I, ll>, u64, u128>;
  Data A = floor_monoid_product<Mono, Data, U>(Mono::to_x(), Mono::to_y(), N, a, b, mod);
  Data offset = Mono::unit();
  offset.dx = T(ADD_X), offset.dy = T(ADD_Y);
  A = Mono::op(offset, A);
  return A.dp;
};
#line 2 "poly/fps_div.hpp"

#line 5 "poly/fps_div.hpp"

// f/g. f の長さで出力される.
template <typename mint, bool SPARSE = false>
vc<mint> fps_div(vc<mint> f, vc<mint> g) {
  if (SPARSE || count_terms(g) < 200) return fps_div_sparse(f, g);
  int n = len(f);
  g.resize(n);
  g = fps_inv<mint>(g);
  f = convolution(f, g);
  f.resize(n);
  return f;
}

// f/g ただし g は sparse
template <typename mint>
vc<mint> fps_div_sparse(vc<mint> f, vc<mint>& g) {
  if (g[0] != mint(1)) {
    mint cf = g[0].inverse();
    for (auto&& x: f) x *= cf;
    for (auto&& x: g) x *= cf;
  }

  vc<pair<int, mint>> dat;
  FOR(i, 1, len(g)) if (g[i] != mint(0)) dat.eb(i, -g[i]);
  FOR(i, len(f)) {
    for (auto&& [j, x]: dat) {
      if (i >= j) f[i] += x * f[i - j];
    }
  }
  return f;
}
#line 2 "nt/primetable.hpp"

template <typename T = int>
vc<T> primetable(int LIM) {
  ++LIM;
  const int S = 32768;
  static int done = 2;
  static vc<T> primes = {2}, sieve(S + 1);

  if (done < LIM) {
    done = LIM;

    primes = {2}, sieve.assign(S + 1, 0);
    const int R = LIM / 2;
    primes.reserve(int(LIM / log(LIM) * 1.1));
    vc<pair<int, int>> cp;
    for (int i = 3; i <= S; i += 2) {
      if (!sieve[i]) {
        cp.eb(i, i * i / 2);
        for (int j = i * i; j <= S; j += 2 * i) sieve[j] = 1;
      }
    }
    for (int L = 1; L <= R; L += S) {
      array<bool, S> block{};
      for (auto& [p, idx]: cp)
        for (int i = idx; i < S + L; idx = (i += p)) block[i - L] = 1;
      FOR(i, min(S, R - L)) if (!block[i]) primes.eb((L + i) * 2 + 1);
    }
  }
  int k = LB(primes, LIM + 1);
  return {primes.begin(), primes.begin() + k};
}
#line 3 "mod/powertable.hpp"

// a^0, ..., a^N

template <typename mint>
vc<mint> powertable_1(mint a, ll N) {
  // table of a^i

  vc<mint> f(N + 1, 1);
  FOR(i, N) f[i + 1] = a * f[i];
  return f;
}

// 0^e, ..., N^e

template <typename mint>
vc<mint> powertable_2(ll e, ll N) {
  auto primes = primetable(N);
  vc<mint> f(N + 1, 1);
  f[0] = mint(0).pow(e);
  for (auto&& p: primes) {
    if (p > N) break;
    mint xp = mint(p).pow(e);
    ll pp = p;
    while (pp <= N) {
      ll i = pp;
      while (i <= N) {
        f[i] *= xp;
        i += pp;
      }
      pp *= p;
    }
  }
  return f;
}
#line 3 "seq/famous/bernoulli.hpp"

template <typename mint>
vc<mint> bernoulli_number(int N) {
  int n = N / 2;
  vc<mint> F(n + 1), G(n + 1);
  mint pow = 1;
  FOR(i, n + 1) {
    F[i] = fact_inv<mint>(2 * i) * pow;
    G[i] = fact_inv<mint>(2 * i + 1) * pow;
    pow *= inv<mint>(4);
  }
  F = fps_div<mint>(F, G);
  vc<mint> B(N + 1);
  if (1 <= N) B[1] = -inv<mint>(2);
  FOR(i, n + 1) B[2 * i] = F[i] * fact<mint>(2 * i);
  return B;
}

template <typename mint>
mint single_bernoulli(int n) {
  // https://atcoder.jp/contests/xmascon23/tasks/xmascon23_e
  if (n == 0) return 1;
  if (n == 1) return -inv<mint>(2);
  /*
  B_n = [x^n/n!] x / (exp(x)-1) = F(1-e^x)
  F(x) = 1+(1/2)x+(1/3)x^2+...
  これを x^n で打ち切る
  F(x) = 1+(1/2)x+(1/3)x^2+...+(1/n+1)x^n, G(x) = F(1-x)
  (xF(x)) d/dx = 1-x^{n+1}/1-x
  ((1-x)G(x)) -d/dx = 1-(1-x)^{n+1}/x = H(x)
  */
  vc<mint> G(n + 2);
  mint sm = 0;
  FOR(i, 1, n + 2) {
    mint c = C<mint>(n + 1, i);
    mint h = (i % 2 == 0 ? c : -c);
    // H(x) = ... gx^{i-1}
    G[i] = h * inv<mint>(i);
    sm += inv<mint>(i);
  }
  G[0] = sm;
  FOR(i, n) G[i + 1] += G[i];
  vc<mint> pow = powertable_2<mint>(n, n);
  mint ans = 0;
  FOR(i, n + 1) { ans += pow[i] * G[i]; }
  return ans;
}
#line 2 "seq/famous/faulhaber.hpp"

// sum_[1,n]i^k=f(n)
template <typename mint>
vc<mint> faulhaber_formula(int k) {
  vc<mint> F = bernoulli_number<mint>(k + 1);
  if (1 <= k) F[1] = inv<mint>(2);
  reverse(all(F));
  F[0] = 0;
  FOR(r, k + 1) { F[k - r + 1] *= fact<mint>(k) * fact_inv<mint>(r) * fact_inv<mint>(k + 1 - r); }
  return F;
}

// sum_[1,n]i^k=f(n)
template <typename mint>
vvc<mint> faulhaber_formula_2d(int n) {
  vc<mint> B = bernoulli_number<mint>(n);
  if (1 <= n) B[1] = inv<mint>(2);
  vvc<mint> ANS(n + 1);
  FOR(k, n + 1) {
    ANS[k].resize(k + 2);
    FOR(j, k + 1) ANS[k][k + 1 - j] = inv<mint>(k + 1) * C<mint>(k + 1, j) * B[j];
  }
  return ANS;
}
#line 4 "convex/lattice_point_sum_polynomial.hpp"

// ax+by<=c という半平面たち. 非有界は 0 埋め.
// 格子点 (x,y) に対して x^iy^j の sum. i<=K1, j<=K2
template <typename mint, int K1, int K2>
array<array<mint, K2 + 1>, K1 + 1> lattice_point_sum_polynomial(vc<tuple<ll, ll, ll>> LINE) {
  ll L = -infty<ll>, R = infty<ll>;
  vc<tuple<ll, ll, ll>> LINE1, LINE2;
  for (auto& [a, b, c]: LINE) {
    if (b == 0) {
      assert(a != 0);
      if (a > 0) { chmin(R, floor<ll>(c, a) + 1); }
      elif (a < 0) { chmax(L, ceil<ll>(-c, -a)); }
    } else {
      if (b > 0) { LINE2.eb(-a, c, b); }
      if (b < 0) { LINE1.eb(a, -c, -b); }
    }
  }
  if (L >= R) { return {}; }
  if (LINE1.empty() || LINE2.empty()) return {};

  auto LOWER = line_max_function_rational(LINE1, L, R);
  auto UPPER = line_min_function_rational(LINE2, L, R);

  array<array<mint, K2 + 2>, K1 + 1> S;
  FOR(i, K1 + 1) FOR(j, K2 + 1) S[i][j] = 0;

  bool bad = 0;
  auto wk = [&](ll L, ll R, ll a1, ll b1, ll c1, ll a2, ll b2, ll c2) -> void {
    // 交点 t/s
    i128 s = i128(a2) * c1 - i128(a1) * c2;
    i128 t = i128(b1) * c2 - i128(b2) * c1;
    if (s == 0) {
      if (t > 0) return;
    }
    if (s > 0) { chmax(L, ceil<i128>(t, s)); }
    if (s < 0) { chmin(R, floor<i128>(-t, -s) + 1); }
    if (L >= R) return;
    if (L == -infty<ll> || R == infty<ll>) {
      bad = 1;
      return;
    }
    auto ADD = floor_sum_of_linear_polynomial<mint, K1, K2 + 1, ll>(L, R, a2, b2, c2);
    auto SUB = floor_sum_of_linear_polynomial<mint, K1, K2 + 1, ll>(L, R, a1, b1 - 1, c1);
    FOR(i, K1 + 1) FOR(j, K2 + 2) S[i][j] += ADD[i][j] - SUB[i][j];
  };

  merge_58(LOWER, UPPER, wk);
  array<array<mint, K2 + 1>, K1 + 1> ANS{};
  if (bad) return ANS;

  static vvc<mint> CF;
  if (CF.empty()) { CF = faulhaber_formula_2d<mint>(K2); }
  FOR(i, K1 + 1) {
    FOR(j, K2 + 1) {
      FOR(k, j + 2) { ANS[i][j] += CF[j][k] * S[i][k]; }
    }
  }
  return ANS;
}
#line 1 "alg/monoid/monoid_for_floor_sum_pq.hpp"

// floor path で (x,y) から x 方向に進むときに p^xq^yx^k1y^k2 を足す
template <typename T, int K1, int K2>
struct Monoid_for_floor_sum_pq {
  using ARR = array<array<T, K2 + 1>, K1 + 1>;
  struct Data {
    ARR dp;
    T dx, dy;
    T prod;
  };

  static pair<T, T> &get_pq() {
    static pair<T, T> pq = {T(1), T(1)};
    return pq;
  }

  static void set_pq(T p, T q) { get_pq() = {p, q}; }

  using value_type = Data;
  using X = value_type;
  static X op(X a, X b) {
    static constexpr int n = max(K1, K2);
    static T comb[n + 1][n + 1];
    if (comb[0][0] != T(1)) {
      comb[0][0] = T(1);
      FOR(i, n) FOR(j, i + 1) { comb[i + 1][j] += comb[i][j], comb[i + 1][j + 1] += comb[i][j]; }
    }

    array<T, K1 + 1> pow_x;
    array<T, K2 + 1> pow_y;
    pow_x[0] = 1, pow_y[0] = 1;
    FOR(i, K1) pow_x[i + 1] = pow_x[i] * a.dx;
    FOR(i, K2) pow_y[i + 1] = pow_y[i] * a.dy;

    FOR(i, K1 + 1) FOR(j, K2 + 1) { b.dp[i][j] *= a.prod; }

    // +dy
    FOR(i, K1 + 1) {
      FOR_R(j, K2 + 1) {
        T x = b.dp[i][j];
        FOR(k, j + 1, K2 + 1) b.dp[i][k] += comb[k][j] * pow_y[k - j] * x;
      }
    }

    // +dx
    FOR(j, K2 + 1) {
      FOR_R(i, K1 + 1) { FOR(k, i, K1 + 1) a.dp[k][j] += comb[k][i] * pow_x[k - i] * b.dp[i][j]; }
    }

    a.dx += b.dx, a.dy += b.dy, a.prod *= b.prod;
    return a;
  }

  static X to_x() {
    X x = unit();
    x.dp[0][0] = 1, x.dx = 1, x.prod = get_pq().fi;
    return x;
  }
  static X to_y() {
    X x = unit();
    x.dy = 1, x.prod = get_pq().se;
    return x;
  }
  static constexpr X unit() { return {ARR{}, T(0), T(0), T(1)}; }
  static constexpr bool commute = 0;
};
#line 3 "mod/floor_sum_of_linear_polynomial_pq.hpp"

// 全部非負, T は答, U は ax+b がオーバーフローしない
template <typename T, int K1, int K2, typename U>
array<array<T, K2 + 1>, K1 + 1> floor_sum_of_linear_polynomial_nonnegative_pq(T p, T q, U N, U a, U b, U mod) {
  static_assert(is_same_v<U, u64> || is_same_v<U, u128>);
  assert(a == 0 || N < (U(-1) - b) / a);
  using Mono = Monoid_for_floor_sum_pq<T, K1, K2>;
  Mono::set_pq(p, q);
  auto x = floor_monoid_product<Mono>(Mono::to_x(), Mono::to_y(), N, a, b, mod);
  return x.dp;
};

// sum_{L<=x<R} x^i floor(ax+b,mod)^j
// a+bx が I, U でオーバーフローしない
template <typename T, int K1, int K2, typename I>
array<array<T, K2 + 1>, K1 + 1> floor_sum_of_linear_polynomial_pq(T p, T q, I L, I R, I a, I b, I mod) {
  static_assert(is_same_v<I, ll> || is_same_v<I, i128>);
  assert(L <= R && mod > 0);
  using Mono = Monoid_for_floor_sum_pq<T, K1, K2>;
  Mono::set_pq(p, q);
  if (a < 0) {
    auto ANS = floor_sum_of_linear_polynomial_pq<T, K1, K2, I>(p.inverse(), q, -R + 1, -L + 1, -a, b, mod);
    FOR(i, K1 + 1) {
      if (i % 2 == 1) { FOR(j, K2 + 1) ANS[i][j] = -ANS[i][j]; }
    }
    return ANS;
  }
  assert(a >= 0);
  I ADD_X = L;
  I N = R - L;
  b += a * L;
  I ADD_Y = floor<I>(b, mod);
  b -= ADD_Y * mod;
  assert(a >= 0 && b >= 0);

  using Mono = Monoid_for_floor_sum_pq<T, K1, K2>;
  using Data = typename Mono::Data;
  using U = std::conditional_t<is_same_v<I, ll>, u64, u128>;
  Data A = floor_monoid_product<Mono, Data, U>(Mono::to_x(), Mono::to_y(), N, a, b, mod);
  Data offset = Mono::unit();
  offset.dx = T(ADD_X), offset.dy = T(ADD_Y);
  A = Mono::op(offset, A);
  T mul = p.pow(ADD_X) * q.pow(ADD_Y);
  FOR(i, K1 + 1) FOR(j, K2 + 1) A.dp[i][j] *= mul;
  return A.dp;
};
#line 7 "convex/lattice_point_sum_polynomial_pq.hpp"

// 格子点 (x,y) に対して p^xq^yx^iy^j の sum. i<=K, j<=L
template <typename mint, int K1, int K2>
array<array<mint, K2 + 1>, K1 + 1> lattice_point_sum_polynomial_pq(mint p, mint q, vc<tuple<ll, ll, ll>> LINE) {
  assert(p != 0 && q != 0);
  if (p == 1 && q == 1) return lattice_point_sum_polynomial<mint, K1, K2>(LINE);
  if (q == 1) {
    for (auto& [a, b, c]: LINE) { swap(a, b); }
    auto tmp = lattice_point_sum_polynomial_pq<mint, K2, K1>(q, p, LINE);
    array<array<mint, K2 + 1>, K1 + 1> ANS{};
    FOR(i, K1 + 1) FOR(j, K2 + 1) ANS[i][j] = tmp[j][i];
    return ANS;
  }
  assert(q != 1);
  ll L = -infty<ll>, R = infty<ll>;
  vc<tuple<ll, ll, ll>> LINE1, LINE2;
  for (auto& [a, b, c]: LINE) {
    if (b == 0) {
      assert(a != 0);
      if (a > 0) { chmin(R, floor<ll>(c, a) + 1); }
      elif (a < 0) { chmax(L, ceil<ll>(-c, -a)); }
    } else {
      if (b > 0) { LINE2.eb(-a, c, b); }
      if (b < 0) { LINE1.eb(a, -c, -b); }
    }
  }
  if (L >= R) { return {}; }
  if (LINE1.empty() || LINE2.empty()) return {};

  auto LOWER = line_max_function_rational(LINE1, L, R);
  auto UPPER = line_min_function_rational(LINE2, L, R);

  array<array<mint, K2 + 1>, K1 + 1> S{};
  array<mint, K1 + 1> T{};

  bool bad = 0;
  auto wk = [&](ll L, ll R, ll a1, ll b1, ll c1, ll a2, ll b2, ll c2) -> void {
    if (bad) return;
    // 交点 t/s
    i128 s = i128(a2) * c1 - i128(a1) * c2;
    i128 t = i128(b1) * c2 - i128(b2) * c1;
    if (s == 0) {
      if (t > 0) return;
    }
    if (s > 0) { chmax(L, ceil<i128>(t, s)); }
    if (s < 0) { chmin(R, floor<i128>(-t, -s) + 1); }
    if (L >= R) return;
    if (L == -infty<ll> || R == infty<ll>) {
      bad = 1;
      return;
    }
    b2 = b2 + c2, b1 = b1 - 1 + c1;
    auto ADD_S = floor_sum_of_linear_polynomial_pq<mint, K1, K2, ll>(p, q, L, R, a2, b2, c2);
    auto SUB_S = floor_sum_of_linear_polynomial_pq<mint, K1, K2, ll>(p, q, L, R, a1, b1, c1);
    auto ADD_T = floor_sum_of_linear_polynomial_pq<mint, K1, 0, ll>(p, 1, L, R, a2, b2, c2);
    auto SUB_T = floor_sum_of_linear_polynomial_pq<mint, K1, 0, ll>(p, 1, L, R, a1, b1, c1);
    FOR(i, K1 + 1) FOR(j, K2 + 1) S[i][j] += ADD_S[i][j] - SUB_S[i][j];
    FOR(i, K1 + 1) T[i] += ADD_T[i][0] - SUB_T[i][0];
  };

  merge_58(LOWER, UPPER, wk);
  array<array<mint, K2 + 1>, K1 + 1> ANS{};
  if (bad) return ANS;

  FOR(k, K2 + 1) {
    auto [c, f] = geometric_sequence_sum_formula(q, k);
    FOR(i, K1 + 1) {
      ANS[i][k] += c * T[i];
      FOR(j, k + 1) { ANS[i][k] += f[j] * S[i][j]; }
    }
  }

  return ANS;
}
Back to top page