二項係数のprefix sumの多点評価
(modulo/multipoint-binomial-sum.hpp)

• View this file on GitHub
• Last update: 2023-05-22 22:29:25+09:00
• Include: #include "modulo/multipoint-binomial-sum.hpp"

Depends on

• Mo's algorithm (misc/mo.hpp)
• modulo/binomial.hpp

Verified with

• verify/verify-unit-test/multipoint-binomial-sum.test.cpp

Code

#pragma once

#include "../misc/mo.hpp"
#include "binomial.hpp"

template <typename mint>
vector<mint> multipoint_binomial_sum(const vector<pair<int, int>>& qs) {
  int N = 2;
  for (auto& p : qs) N = max(N, p.first);
  Binomial<mint> b(N + 1);
  int Q = qs.size();
  Mo mo(N, Q);
  for (auto& p : qs) {
    assert(p.second <= p.first);
    assert(p.first <= N);
    mo.insert(p.second, p.first);
  }
  vector<mint> ans(Q);
  mint cur = 1;
  int n = 0, m = 0;
  auto al = [&](int) { cur -= b.C(n, m--); };
  auto ar = [&](int) { cur += cur - b.C(n++, m); };
  auto el = [&](int) { cur += b.C(n, ++m); };
  auto er = [&](int) { cur = (cur + b.C(--n, m)) * b.inv(2); };
  auto q = [&](int i) { ans[i] = cur; };
  mo.run(al, ar, el, er, q);
  return ans;
}

/**
 * @brief 二項係数のprefix sumの多点評価
 */

#line 2 "modulo/multipoint-binomial-sum.hpp"

#line 2 "misc/mo.hpp"

struct Mo {
  int width;
  vector<int> left, right, order;

  Mo(int N, int Q) : order(Q) {
    width = max<int>(1, 1.0 * N / max<double>(1.0, sqrt(Q * 2.0 / 3.0)));
    iota(begin(order), end(order), 0);
  }

  void insert(int l, int r) { /* [l, r) */
    left.emplace_back(l);
    right.emplace_back(r);
  }

  template <typename AL, typename AR, typename DL, typename DR, typename REM>
  void run(const AL &add_left, const AR &add_right, const DL &delete_left,
           const DR &delete_right, const REM &rem) {
    assert(left.size() == order.size());
    sort(begin(order), end(order), [&](int a, int b) {
      int ablock = left[a] / width, bblock = left[b] / width;
      if (ablock != bblock) return ablock < bblock;
      if (ablock & 1) return right[a] < right[b];
      return right[a] > right[b];
    });
    int nl = 0, nr = 0;
    for (auto idx : order) {
      while (nl > left[idx]) add_left(--nl);
      while (nr < right[idx]) add_right(nr++);
      while (nl < left[idx]) delete_left(nl++);
      while (nr > right[idx]) delete_right(--nr);
      rem(idx);
    }
  }
};

/**
 * @brief Mo's algorithm
 * @docs docs/misc/mo.md
 */
#line 2 "modulo/binomial.hpp"

#include <cassert>
#include <type_traits>
#include <vector>
using namespace std;

// コンストラクタの MAX に 「C(n, r) や fac(n) でクエリを投げる最大の n 」
// を入れると倍速くらいになる
// mod を超えて前計算して 0 割りを踏むバグは対策済み
template <typename T>
struct Binomial {
  vector<T> f, g, h;
  Binomial(int MAX = 0) {
    assert(T::get_mod() != 0 && "Binomial<mint>()");
    f.resize(1, T{1});
    g.resize(1, T{1});
    h.resize(1, T{1});
    if (MAX > 0) extend(MAX + 1);
  }

  void extend(int m = -1) {
    int n = f.size();
    if (m == -1) m = n * 2;
    m = min<int>(m, T::get_mod());
    if (n >= m) return;
    f.resize(m);
    g.resize(m);
    h.resize(m);
    for (int i = n; i < m; i++) f[i] = f[i - 1] * T(i);
    g[m - 1] = f[m - 1].inverse();
    h[m - 1] = g[m - 1] * f[m - 2];
    for (int i = m - 2; i >= n; i--) {
      g[i] = g[i + 1] * T(i + 1);
      h[i] = g[i] * f[i - 1];
    }
  }

  T fac(int i) {
    if (i < 0) return T(0);
    while (i >= (int)f.size()) extend();
    return f[i];
  }

  T finv(int i) {
    if (i < 0) return T(0);
    while (i >= (int)g.size()) extend();
    return g[i];
  }

  T inv(int i) {
    if (i < 0) return -inv(-i);
    while (i >= (int)h.size()) extend();
    return h[i];
  }

  T C(int n, int r) {
    if (n < 0 || n < r || r < 0) return T(0);
    return fac(n) * finv(n - r) * finv(r);
  }

  inline T operator()(int n, int r) { return C(n, r); }

  template <typename I>
  T multinomial(const vector<I>& r) {
    static_assert(is_integral<I>::value == true);
    int n = 0;
    for (auto& x : r) {
      if (x < 0) return T(0);
      n += x;
    }
    T res = fac(n);
    for (auto& x : r) res *= finv(x);
    return res;
  }

  template <typename I>
  T operator()(const vector<I>& r) {
    return multinomial(r);
  }

  T C_naive(int n, int r) {
    if (n < 0 || n < r || r < 0) return T(0);
    T ret = T(1);
    r = min(r, n - r);
    for (int i = 1; i <= r; ++i) ret *= inv(i) * (n--);
    return ret;
  }

  T P(int n, int r) {
    if (n < 0 || n < r || r < 0) return T(0);
    return fac(n) * finv(n - r);
  }

  // [x^r] 1 / (1-x)^n
  T H(int n, int r) {
    if (n < 0 || r < 0) return T(0);
    return r == 0 ? 1 : C(n + r - 1, r);
  }
};
#line 5 "modulo/multipoint-binomial-sum.hpp"

template <typename mint>
vector<mint> multipoint_binomial_sum(const vector<pair<int, int>>& qs) {
  int N = 2;
  for (auto& p : qs) N = max(N, p.first);
  Binomial<mint> b(N + 1);
  int Q = qs.size();
  Mo mo(N, Q);
  for (auto& p : qs) {
    assert(p.second <= p.first);
    assert(p.first <= N);
    mo.insert(p.second, p.first);
  }
  vector<mint> ans(Q);
  mint cur = 1;
  int n = 0, m = 0;
  auto al = [&](int) { cur -= b.C(n, m--); };
  auto ar = [&](int) { cur += cur - b.C(n++, m); };
  auto el = [&](int) { cur += b.C(n, ++m); };
  auto er = [&](int) { cur = (cur + b.C(--n, m)) * b.inv(2); };
  auto q = [&](int i) { ans[i] = cur; };
  mo.run(al, ar, el, er, q);
  return ans;
}

/**
 * @brief 二項係数のprefix sumの多点評価
 */

Back to top page