ds/offline_query/coeffient_query_2d.hpp

• View this file on GitHub
• Last update: 2023-11-01 01:33:38+09:00
• Include: #include "ds/offline_query/coeffient_query_2d.hpp"

Depends on

• alg/monoid/add.hpp
• ds/fenwicktree/fenwicktree.hpp

Required by

• ds/offline_query/rectangle_add_rectangle_sum.hpp

Verified with

• test/library_checker/datastructure/rectangle_sum_cf2d.test.cpp
• test/library_checker/datastructure/static_rectangle_add_rectangle_sum.test.cpp
• test/yukicoder/1490.test.cpp

Code

#include "ds/fenwicktree/fenwicktree.hpp"

// A, B：定数
// Sparse Laurent Polynomial f(x,y) を与える
// [x^py^q] f(x,y)/(1-x)^A(1-y)^B をたくさん求める
// O(AB N logN) 時間
template <int A, int B, typename T, typename XY>
struct Coefficient_Query_2D {
  struct Mono {
    using value_type = array<T, A * B>;
    using X = value_type;
    static X op(X x, X y) {
      FOR(i, A * B) x[i] += y[i];
      return x;
    }
    static constexpr X unit() { return X{}; }
    static constexpr bool commute = 1;
  };

  vc<tuple<XY, XY, T>> F;
  vc<pair<XY, XY>> QUERY;

  Coefficient_Query_2D() {}
  void add_query(XY x, XY y, T c) { F.eb(x, y, c); }
  void sum_query(XY p, XY q) { QUERY.eb(p, q); }

  // div_fact：最後に (A-1)!(B-1)! で割るかどうか。ふつうは割る。
  vc<T> calc(bool div_fact = true) {
    // 加算する点の x について座圧
    sort(all(F),
         [&](auto& a, auto& b) -> bool { return get<0>(a) < get<0>(b); });
    vc<XY> keyX;
    keyX.reserve(len(F));
    for (auto&& [a, b, c]: F) {
      if (keyX.empty() || keyX.back() != a) keyX.eb(a);
      a = len(keyX) - 1;
    }
    keyX.shrink_to_fit();
    // y 昇順にクエリ処理する
    const int Q = len(QUERY);
    vc<int> I(Q);
    iota(all(I), 0);
    sort(all(I),
         [&](auto& a, auto& b) -> bool { return QUERY[a].se < QUERY[b].se; });
    sort(all(F),
         [&](auto& a, auto& b) -> bool { return get<1>(a) < get<1>(b); });
    FenwickTree<Mono> bit(len(keyX));
    vc<T> res(Q);
    int ptr = 0;
    for (auto&& qid: I) {
      auto [p, q] = QUERY[qid];
      // y <= q となる F の加算
      while (ptr < len(F) && get<1>(F[ptr]) <= q) {
        auto& [ia, b, w] = F[ptr++];
        XY a = keyX[ia];
        // w(p-a+1)...(p-a+A-1)(q-b+1)...(q-b+B-1) を p,q の多項式として
        vc<T> f(A), g(B);
        f[0] = w, g[0] = 1;
        FOR(i, A - 1) { FOR_R(j, i + 1) f[j + 1] += f[j] * T(-a + 1 + i); }
        FOR(i, B - 1) { FOR_R(j, i + 1) g[j + 1] += g[j] * T(-b + 1 + i); }
        reverse(all(f)), reverse(all(g));
        array<T, A * B> G{};
        FOR(i, A) FOR(j, B) G[B * i + j] = f[i] * g[j];
        bit.add(ia, G);
      }
      auto SM = bit.sum(UB(keyX, p));
      T sm = 0, pow_p = 1;
      FOR(i, A) {
        T prod = pow_p;
        FOR(j, B) { sm += prod * SM[B * i + j], prod *= T(q); }
        pow_p *= T(p);
      }
      res[qid] = sm;
    }
    if (div_fact && (A >= 3 || B >= 3)) {
      T cf = T(1);
      FOR(a, 1, A) cf *= T(a);
      FOR(b, 1, B) cf *= T(b);
      for (auto&& x: res) x /= cf;
    }
    return res;
  }
};

#line 2 "alg/monoid/add.hpp"

template <typename X>
struct Monoid_Add {
  using value_type = X;
  static constexpr X op(const X &x, const X &y) noexcept { return x + y; }
  static constexpr X inverse(const X &x) noexcept { return -x; }
  static constexpr X power(const X &x, ll n) noexcept { return X(n) * x; }
  static constexpr X unit() { return X(0); }
  static constexpr bool commute = true;
};
#line 3 "ds/fenwicktree/fenwicktree.hpp"

template <typename Monoid>
struct FenwickTree {
  using G = Monoid;
  using E = typename G::value_type;
  int n;
  vector<E> dat;
  E total;

  FenwickTree() {}
  FenwickTree(int n) { build(n); }
  template <typename F>
  FenwickTree(int n, F f) {
    build(n, f);
  }
  FenwickTree(const vc<E>& v) { build(v); }

  void build(int m) {
    n = m;
    dat.assign(m, G::unit());
    total = G::unit();
  }
  void build(const vc<E>& v) {
    build(len(v), [&](int i) -> E { return v[i]; });
  }
  template <typename F>
  void build(int m, F f) {
    n = m;
    dat.clear();
    dat.reserve(n);
    total = G::unit();
    FOR(i, n) { dat.eb(f(i)); }
    for (int i = 1; i <= n; ++i) {
      int j = i + (i & -i);
      if (j <= n) dat[j - 1] = G::op(dat[i - 1], dat[j - 1]);
    }
    total = prefix_sum(m);
  }

  E prod_all() { return total; }
  E sum_all() { return total; }
  E sum(int k) { return prefix_sum(k); }
  E prod(int k) { return prefix_prod(k); }
  E prefix_sum(int k) { return prefix_prod(k); }
  E prefix_prod(int k) {
    chmin(k, n);
    E ret = G::unit();
    for (; k > 0; k -= k & -k) ret = G::op(ret, dat[k - 1]);
    return ret;
  }
  E sum(int L, int R) { return prod(L, R); }
  E prod(int L, int R) {
    chmax(L, 0), chmin(R, n);
    if (L == 0) return prefix_prod(R);
    assert(0 <= L && L <= R && R <= n);
    E pos = G::unit(), neg = G::unit();
    while (L < R) { pos = G::op(pos, dat[R - 1]), R -= R & -R; }
    while (R < L) { neg = G::op(neg, dat[L - 1]), L -= L & -L; }
    return G::op(pos, G::inverse(neg));
  }

  void add(int k, E x) { multiply(k, x); }
  void multiply(int k, E x) {
    static_assert(G::commute);
    total = G::op(total, x);
    for (++k; k <= n; k += k & -k) dat[k - 1] = G::op(dat[k - 1], x);
  }

  template <class F>
  int max_right(const F check) {
    assert(check(G::unit()));
    int i = 0;
    E s = G::unit();
    int k = 1;
    while (2 * k <= n) k *= 2;
    while (k) {
      if (i + k - 1 < len(dat)) {
        E t = G::op(s, dat[i + k - 1]);
        if (check(t)) { i += k, s = t; }
      }
      k >>= 1;
    }
    return i;
  }

  // check(i, x)
  template <class F>
  int max_right_with_index(const F check) {
    assert(check(0, G::unit()));
    int i = 0;
    E s = G::unit();
    int k = 1;
    while (2 * k <= n) k *= 2;
    while (k) {
      if (i + k - 1 < len(dat)) {
        E t = G::op(s, dat[i + k - 1]);
        if (check(i + k, t)) { i += k, s = t; }
      }
      k >>= 1;
    }
    return i;
  }

  int kth(E k) {
    return max_right([&k](E x) -> bool { return x <= k; });
  }
};
#line 2 "ds/offline_query/coeffient_query_2d.hpp"

// A, B：定数
// Sparse Laurent Polynomial f(x,y) を与える
// [x^py^q] f(x,y)/(1-x)^A(1-y)^B をたくさん求める
// O(AB N logN) 時間
template <int A, int B, typename T, typename XY>
struct Coefficient_Query_2D {
  struct Mono {
    using value_type = array<T, A * B>;
    using X = value_type;
    static X op(X x, X y) {
      FOR(i, A * B) x[i] += y[i];
      return x;
    }
    static constexpr X unit() { return X{}; }
    static constexpr bool commute = 1;
  };

  vc<tuple<XY, XY, T>> F;
  vc<pair<XY, XY>> QUERY;

  Coefficient_Query_2D() {}
  void add_query(XY x, XY y, T c) { F.eb(x, y, c); }
  void sum_query(XY p, XY q) { QUERY.eb(p, q); }

  // div_fact：最後に (A-1)!(B-1)! で割るかどうか。ふつうは割る。
  vc<T> calc(bool div_fact = true) {
    // 加算する点の x について座圧
    sort(all(F),
         [&](auto& a, auto& b) -> bool { return get<0>(a) < get<0>(b); });
    vc<XY> keyX;
    keyX.reserve(len(F));
    for (auto&& [a, b, c]: F) {
      if (keyX.empty() || keyX.back() != a) keyX.eb(a);
      a = len(keyX) - 1;
    }
    keyX.shrink_to_fit();
    // y 昇順にクエリ処理する
    const int Q = len(QUERY);
    vc<int> I(Q);
    iota(all(I), 0);
    sort(all(I),
         [&](auto& a, auto& b) -> bool { return QUERY[a].se < QUERY[b].se; });
    sort(all(F),
         [&](auto& a, auto& b) -> bool { return get<1>(a) < get<1>(b); });
    FenwickTree<Mono> bit(len(keyX));
    vc<T> res(Q);
    int ptr = 0;
    for (auto&& qid: I) {
      auto [p, q] = QUERY[qid];
      // y <= q となる F の加算
      while (ptr < len(F) && get<1>(F[ptr]) <= q) {
        auto& [ia, b, w] = F[ptr++];
        XY a = keyX[ia];
        // w(p-a+1)...(p-a+A-1)(q-b+1)...(q-b+B-1) を p,q の多項式として
        vc<T> f(A), g(B);
        f[0] = w, g[0] = 1;
        FOR(i, A - 1) { FOR_R(j, i + 1) f[j + 1] += f[j] * T(-a + 1 + i); }
        FOR(i, B - 1) { FOR_R(j, i + 1) g[j + 1] += g[j] * T(-b + 1 + i); }
        reverse(all(f)), reverse(all(g));
        array<T, A * B> G{};
        FOR(i, A) FOR(j, B) G[B * i + j] = f[i] * g[j];
        bit.add(ia, G);
      }
      auto SM = bit.sum(UB(keyX, p));
      T sm = 0, pow_p = 1;
      FOR(i, A) {
        T prod = pow_p;
        FOR(j, B) { sm += prod * SM[B * i + j], prod *= T(q); }
        pow_p *= T(p);
      }
      res[qid] = sm;
    }
    if (div_fact && (A >= 3 || B >= 3)) {
      T cf = T(1);
      FOR(a, 1, A) cf *= T(a);
      FOR(b, 1, B) cf *= T(b);
      for (auto&& x: res) x /= cf;
    }
    return res;
  }
};

Back to top page