二項係数の和（オンライン） (math/combinatorics/online-binomial-sum.hpp)

View this file on GitHub
View document part on GitHub
Last update: 2026-07-11 16:54:42+09:00
Include: #include "math/combinatorics/online-binomial-sum.hpp"

概要

整数 $n,m$ と重み $r$ に対する二項係数の prefix sum を $\displaystyle F(n,m)=\sum_{i=0}^{n-1}r^i\binom{m}{i}$ とおく。
半開区間の左端、右端をそれぞれ $l,u$ とする。 $\displaystyle \sum_{i=l}^{u-1}r^i\binom{m}{i}$ をオンラインで求める。
$\binom{m}{i}=0\;(i>m)$ として扱う。
$n$ と $m$ をそれぞれバケットに分け、バケット境界と各軸の上端をサンプル座標とする。両軸のサンプル座標の直積で $F(n,m)$ を前計算する。
クエリでは、マンハッタン距離が最小であるサンプル点から $n$ 方向と $m$ 方向に遷移する。
バケットサイズを指定できる。

以下、クエリで与えられる $m$ の最大値を $M$ 、バケットサイズを $B$ とする。

$r=0$ の場合は閉形式を用いる。
$r=-1$ の場合は $\displaystyle F(n,m)=(-1)^{n-1}\binom{m-1}{n-1}$ を用い、 $r+1$ による除算を行わない。ただし、 $n=0$ または $m=0$ の場合は別に処理する。
$r$ が $0$ でも $-1$ でもない場合、 $n$ 方向の遷移には境界項 $\displaystyle r^n\binom{m}{n}$ を用いる。 $m$ 方向の遷移にはパスカルの三角形の等式から得られる $\displaystyle F(n,m+1)=(r+1)F(n,m)-r^n\binom{m}{n-1}$ を用いる。
バケットサイズ $B$ を指定しない場合、 $r=0$ であれば $B=1$ とし、既定のバケットサイズを計算しない。
バケットサイズ $B$ を指定せず、 $r\ne 0$ の場合、 $B^2>M$ を満たす最小の $2$ の冪が選ばれる。

使い方

max_m を $M$ 、bucket_size を $B$ とおく。

OnlineBinomialSum(int max_m, T r, int bucket_size)
- $0\le m\le M$ のクエリに対する前計算を行う。
- r は重みである。
- 引数の bucket_size をバケットサイズとして前計算を行う。
- 前提: $M\ge 0,\;B>0$ 。
- 前提: T は素数 $p$ を法とする体の型であり、整数からの構築、四則演算、等値比較を持つ。
- 前提: std::numeric_limits<T>::is_integer は false である。
- 前提: $M<p$ 。
- 備考: 合成数を法とする型は対象としない。
OnlineBinomialSum(int max_m, T r = T(1))
- $0\le m\le M$ のクエリに対する前計算を行う。
- r は重みであり、省略時は $1$ である。
- $r=0$ の場合、バケットサイズ $B$ は $1$ であり、既定のバケットサイズを計算しない。
- $r\ne 0$ の場合、バケットサイズ $B$ は、 $B^2>M$ を満たす最小の $2$ の冪が選ばれる。
- 前提: $M\ge 0$ 。
- 前提: T は素数 $p$ を法とする体の型であり、整数からの構築、四則演算、等値比較を持つ。
- 前提: std::numeric_limits<T>::is_integer は false である。
- 前提: $M<p$ 。
- 備考: 合成数を法とする型は対象としない。
T binom_prefix_sum(int n, int m) const
- $\displaystyle \sum_{i=0}^{n-1}r^i\binom{m}{i}$ を返す。
- 前提: $n\ge 0,\;0\le m\le M$ 。
- 備考: $n>m+1$ の場合も assert 違反にせず、全体の和を返す。
T binom_sum(int l, int u, int m) const
- $\displaystyle \sum_{i=l}^{u-1}r^i\binom{m}{i}$ を返す。
- 前提: $0\le l\le u,\;0\le m\le M$ 。
- 備考: $u>m+1$ または $l>m$ の場合も assert 違反にしない。

計算量

max_m を $M$ 、bucket_size を $B$ とおく。

$r$ が $0$ でも $-1$ でもない場合、

コンストラクタ: 時間 $O(M^2/B+M)$ 、空間 $O((M/B+1)^2+M)$
binom_prefix_sum(n, m): 時間 $O(B)$
binom_sum(l, u, m): 時間 $O(B)$

である。各軸の隣接するサンプル座標の差は高々 $B$ であるため、最も近いサンプル点からクエリ点までの遷移回数は高々 $B$ である。

$r=-1$ の場合、

コンストラクタ: 時間 $O(M)$ 、空間 $O(M)$
binom_prefix_sum(n, m): 時間 $O(1)$
binom_sum(l, u, m): 時間 $O(1)$

である。

$r=0$ の場合、全ての処理は時間 $O(1)$ 、空間 $O(1)$ である。

コンストラクタで $B$ を指定せず、 $r$ が $0$ でも $-1$ でもない場合は $B=O(\sqrt M)$ である。このとき、

コンストラクタ: 時間 $O(M\sqrt M)$ 、空間 $O(M)$
binom_prefix_sum(n, m): 時間 $O(\sqrt M)$
binom_sum(l, u, m): 時間 $O(\sqrt M)$

である。

Verified with

Code

#ifndef MATH_COMBINATORICS_ONLINE_BINOMIAL_SUM_HPP
#define MATH_COMBINATORICS_ONLINE_BINOMIAL_SUM_HPP

// Σ_{i=l}^{u-1} r^i binom(m,i) をオンラインで求める。
// 0 <= l <= u と 0 <= m <= max_m を仮定する。
// n と m のバケット境界および上端の直積で累積和をサンプルし、
// クエリの点に 2 次元的に最も近いサンプル点から復元する。
// T は素数を法とする体の型である。
// std::numeric_limits<T>::is_integer が false であり、
// 法 p について max_m < p であることを仮定する。
// r = 0 では、バケットサイズを指定しない場合もその計算を行わない。
// 前計算とクエリは時間計算量、空間計算量ともに O(1)。
// r = -1 では交代二項和の閉形式を用いる。
// B をバケットサイズとして、r が 0, -1 でない場合の時間計算量は
// 前計算 O(max_m^2 / B + max_m)、クエリ O(B)。
// 空間計算量は O((max_m / B + 1)^2 + max_m)。

#include <cassert>
#include <limits>
#include <vector>

template <class T> struct OnlineBinomialSum {
    static_assert(!std::numeric_limits<T>::is_integer,
                  "std::numeric_limits<T>::is_integer must be false.");

  public:
    int max_m;
    int bucket_size;
    T r;
    T r_plus_one;
    bool r_is_zero;
    bool r_is_minus_one;
    T r_plus_one_inverse;
    std::vector<T> factorial;
    std::vector<T> inverse_factorial;
    std::vector<T> integer_inverse;
    std::vector<T> power_r;
    std::vector<int> sample_n_list;
    std::vector<int> sample_m_list;
    std::vector<T> sample_sum_table;

    explicit OnlineBinomialSum(int max_m, T r, int bucket_size)
        : max_m(max_m), bucket_size(bucket_size), r(r), r_plus_one(r + T(1)),
          r_is_zero(r == T()), r_is_minus_one(r_plus_one == T()),
          r_plus_one_inverse(T()) {
        assert(max_m >= 0);
        assert(bucket_size > 0);

        if (r_is_zero) {
            return;
        }

        factorial.assign(max_m + 1, T(1));
        for (int i = 1; i <= max_m; ++i) {
            factorial[i] = factorial[i - 1] * T(i);
        }

        inverse_factorial.assign(max_m + 1, T(1));
        inverse_factorial[max_m] = T(1) / factorial[max_m];
        for (int i = max_m; i >= 1; --i) {
            inverse_factorial[i - 1] = inverse_factorial[i] * T(i);
        }

        if (r_is_minus_one) {
            return;
        }

        r_plus_one_inverse = T(1) / r_plus_one;

        integer_inverse.assign(max_m + 1, T());
        for (int i = 1; i <= max_m; ++i) {
            integer_inverse[i] = factorial[i - 1] * inverse_factorial[i];
        }

        power_r.assign(max_m + 2, T());
        power_r[0] = T(1);
        for (int i = 0; i <= max_m; ++i) {
            power_r[i + 1] = power_r[i] * r;
        }

        sample_n_list = make_sample_list(max_m + 1);
        sample_m_list = make_sample_list(max_m);
        sample_sum_table.assign(sample_n_list.size() * sample_m_list.size(),
                                T());

        const int sample_n_count = static_cast<int>(sample_n_list.size());
        const int sample_m_count = static_cast<int>(sample_m_list.size());
        for (int sample_m_index = 0; sample_m_index < sample_m_count;
             ++sample_m_index) {
            const int sample_m = sample_m_list[sample_m_index];
            T sum = T();
            T term = T(1);
            int current_n = 0;

            for (int sample_n_index = 0; sample_n_index < sample_n_count;
                 ++sample_n_index) {
                const int sample_n = sample_n_list[sample_n_index];
                while (current_n < sample_n && current_n <= sample_m) {
                    sum += term;
                    if (current_n < sample_m) {
                        term *= r;
                        term *= T(sample_m - current_n);
                        term *= integer_inverse[current_n + 1];
                    }
                    ++current_n;
                }
                sample_sum_table[sample_m_index * sample_n_count +
                                 sample_n_index] = sum;
            }
        }
    }

    explicit OnlineBinomialSum(int max_m, T r = T(1))
        : OnlineBinomialSum(max_m, r,
                            r == T() ? 1 : default_bucket_size(max_m)) {}

    T binom_prefix_sum(int n, int m) const {
        assert(n >= 0);
        assert(m >= 0);
        assert(m <= max_m);

        return binom_prefix_sum_unchecked(n, m);
    }

    T binom_sum(int l, int u, int m) const {
        assert(l >= 0);
        assert(l <= u);
        assert(m >= 0);
        assert(m <= max_m);

        return binom_prefix_sum_unchecked(u, m) -
               binom_prefix_sum_unchecked(l, m);
    }

  private:
    T binom_prefix_sum_unchecked(int n, int m) const {
        if (n == 0) {
            return T();
        }
        if (n > m) {
            n = m + 1;
        }
        if (r_is_zero) {
            return T(1);
        }
        if (r_is_minus_one) {
            if (m == 0) {
                return T(1);
            }
            if (n > m) {
                return T();
            }

            T ans = binomial(m - 1, n - 1);
            if ((n - 1) % 2 == 1) {
                ans = T() - ans;
            }
            return ans;
        }

        const int sample_n_index = nearest_sample_index(sample_n_list, n);
        const int sample_m_index = nearest_sample_index(sample_m_list, m);
        const int sample_n_count = static_cast<int>(sample_n_list.size());
        int current_n = sample_n_list[sample_n_index];
        int current_m = sample_m_list[sample_m_index];
        T sum =
            sample_sum_table[sample_m_index * sample_n_count + sample_n_index];

        while (current_n < n) {
            if (current_n <= current_m) {
                sum += power_r[current_n] * binomial(current_m, current_n);
            }
            ++current_n;
        }

        while (current_n > n) {
            --current_n;
            if (current_n <= current_m) {
                sum -= power_r[current_n] * binomial(current_m, current_n);
            }
        }

        while (current_m < m) {
            sum *= r_plus_one;
            if (current_n - 1 <= current_m) {
                sum -= power_r[current_n] * binomial(current_m, current_n - 1);
            }
            ++current_m;
        }

        while (current_m > m) {
            --current_m;
            sum += power_r[current_n] * binomial(current_m, current_n - 1);
            sum *= r_plus_one_inverse;
        }

        return sum;
    }

    T binomial(int n, int k) const {
        return factorial[n] * inverse_factorial[k] * inverse_factorial[n - k];
    }

    std::vector<int> make_sample_list(int limit) const {
        const int full_bucket_count = limit / bucket_size;
        std::vector<int> sample_list;
        sample_list.reserve(full_bucket_count + 2);
        for (int index = 0; index <= full_bucket_count; ++index) {
            sample_list.push_back(index * bucket_size);
        }
        if (sample_list.back() != limit) {
            sample_list.push_back(limit);
        }

        return sample_list;
    }

    int nearest_sample_index(const std::vector<int> &sample_list,
                             int value) const {
        int index = value / bucket_size;
        const int sample_count = static_cast<int>(sample_list.size());
        if (index + 1 >= sample_count) {
            return index;
        }
        if (value - sample_list[index] <= sample_list[index + 1] - value) {
            return index;
        }

        return index + 1;
    }

    static int default_bucket_size(int max_m) {
        assert(max_m >= 0);

        int bucket_size = 1;
        while (static_cast<long long>(bucket_size) * bucket_size <= max_m) {
            bucket_size *= 2;
        }

        return bucket_size;
    }
};

#endif

#line 1 "math/combinatorics/online-binomial-sum.hpp"



// Σ_{i=l}^{u-1} r^i binom(m,i) をオンラインで求める。
// 0 <= l <= u と 0 <= m <= max_m を仮定する。
// n と m のバケット境界および上端の直積で累積和をサンプルし、
// クエリの点に 2 次元的に最も近いサンプル点から復元する。
// T は素数を法とする体の型である。
// std::numeric_limits<T>::is_integer が false であり、
// 法 p について max_m < p であることを仮定する。
// r = 0 では、バケットサイズを指定しない場合もその計算を行わない。
// 前計算とクエリは時間計算量、空間計算量ともに O(1)。
// r = -1 では交代二項和の閉形式を用いる。
// B をバケットサイズとして、r が 0, -1 でない場合の時間計算量は
// 前計算 O(max_m^2 / B + max_m)、クエリ O(B)。
// 空間計算量は O((max_m / B + 1)^2 + max_m)。

#include <cassert>
#include <limits>
#include <vector>

template <class T> struct OnlineBinomialSum {
    static_assert(!std::numeric_limits<T>::is_integer,
                  "std::numeric_limits<T>::is_integer must be false.");

  public:
    int max_m;
    int bucket_size;
    T r;
    T r_plus_one;
    bool r_is_zero;
    bool r_is_minus_one;
    T r_plus_one_inverse;
    std::vector<T> factorial;
    std::vector<T> inverse_factorial;
    std::vector<T> integer_inverse;
    std::vector<T> power_r;
    std::vector<int> sample_n_list;
    std::vector<int> sample_m_list;
    std::vector<T> sample_sum_table;

    explicit OnlineBinomialSum(int max_m, T r, int bucket_size)
        : max_m(max_m), bucket_size(bucket_size), r(r), r_plus_one(r + T(1)),
          r_is_zero(r == T()), r_is_minus_one(r_plus_one == T()),
          r_plus_one_inverse(T()) {
        assert(max_m >= 0);
        assert(bucket_size > 0);

        if (r_is_zero) {
            return;
        }

        factorial.assign(max_m + 1, T(1));
        for (int i = 1; i <= max_m; ++i) {
            factorial[i] = factorial[i - 1] * T(i);
        }

        inverse_factorial.assign(max_m + 1, T(1));
        inverse_factorial[max_m] = T(1) / factorial[max_m];
        for (int i = max_m; i >= 1; --i) {
            inverse_factorial[i - 1] = inverse_factorial[i] * T(i);
        }

        if (r_is_minus_one) {
            return;
        }

        r_plus_one_inverse = T(1) / r_plus_one;

        integer_inverse.assign(max_m + 1, T());
        for (int i = 1; i <= max_m; ++i) {
            integer_inverse[i] = factorial[i - 1] * inverse_factorial[i];
        }

        power_r.assign(max_m + 2, T());
        power_r[0] = T(1);
        for (int i = 0; i <= max_m; ++i) {
            power_r[i + 1] = power_r[i] * r;
        }

        sample_n_list = make_sample_list(max_m + 1);
        sample_m_list = make_sample_list(max_m);
        sample_sum_table.assign(sample_n_list.size() * sample_m_list.size(),
                                T());

        const int sample_n_count = static_cast<int>(sample_n_list.size());
        const int sample_m_count = static_cast<int>(sample_m_list.size());
        for (int sample_m_index = 0; sample_m_index < sample_m_count;
             ++sample_m_index) {
            const int sample_m = sample_m_list[sample_m_index];
            T sum = T();
            T term = T(1);
            int current_n = 0;

            for (int sample_n_index = 0; sample_n_index < sample_n_count;
                 ++sample_n_index) {
                const int sample_n = sample_n_list[sample_n_index];
                while (current_n < sample_n && current_n <= sample_m) {
                    sum += term;
                    if (current_n < sample_m) {
                        term *= r;
                        term *= T(sample_m - current_n);
                        term *= integer_inverse[current_n + 1];
                    }
                    ++current_n;
                }
                sample_sum_table[sample_m_index * sample_n_count +
                                 sample_n_index] = sum;
            }
        }
    }

    explicit OnlineBinomialSum(int max_m, T r = T(1))
        : OnlineBinomialSum(max_m, r,
                            r == T() ? 1 : default_bucket_size(max_m)) {}

    T binom_prefix_sum(int n, int m) const {
        assert(n >= 0);
        assert(m >= 0);
        assert(m <= max_m);

        return binom_prefix_sum_unchecked(n, m);
    }

    T binom_sum(int l, int u, int m) const {
        assert(l >= 0);
        assert(l <= u);
        assert(m >= 0);
        assert(m <= max_m);

        return binom_prefix_sum_unchecked(u, m) -
               binom_prefix_sum_unchecked(l, m);
    }

  private:
    T binom_prefix_sum_unchecked(int n, int m) const {
        if (n == 0) {
            return T();
        }
        if (n > m) {
            n = m + 1;
        }
        if (r_is_zero) {
            return T(1);
        }
        if (r_is_minus_one) {
            if (m == 0) {
                return T(1);
            }
            if (n > m) {
                return T();
            }

            T ans = binomial(m - 1, n - 1);
            if ((n - 1) % 2 == 1) {
                ans = T() - ans;
            }
            return ans;
        }

        const int sample_n_index = nearest_sample_index(sample_n_list, n);
        const int sample_m_index = nearest_sample_index(sample_m_list, m);
        const int sample_n_count = static_cast<int>(sample_n_list.size());
        int current_n = sample_n_list[sample_n_index];
        int current_m = sample_m_list[sample_m_index];
        T sum =
            sample_sum_table[sample_m_index * sample_n_count + sample_n_index];

        while (current_n < n) {
            if (current_n <= current_m) {
                sum += power_r[current_n] * binomial(current_m, current_n);
            }
            ++current_n;
        }

        while (current_n > n) {
            --current_n;
            if (current_n <= current_m) {
                sum -= power_r[current_n] * binomial(current_m, current_n);
            }
        }

        while (current_m < m) {
            sum *= r_plus_one;
            if (current_n - 1 <= current_m) {
                sum -= power_r[current_n] * binomial(current_m, current_n - 1);
            }
            ++current_m;
        }

        while (current_m > m) {
            --current_m;
            sum += power_r[current_n] * binomial(current_m, current_n - 1);
            sum *= r_plus_one_inverse;
        }

        return sum;
    }

    T binomial(int n, int k) const {
        return factorial[n] * inverse_factorial[k] * inverse_factorial[n - k];
    }

    std::vector<int> make_sample_list(int limit) const {
        const int full_bucket_count = limit / bucket_size;
        std::vector<int> sample_list;
        sample_list.reserve(full_bucket_count + 2);
        for (int index = 0; index <= full_bucket_count; ++index) {
            sample_list.push_back(index * bucket_size);
        }
        if (sample_list.back() != limit) {
            sample_list.push_back(limit);
        }

        return sample_list;
    }

    int nearest_sample_index(const std::vector<int> &sample_list,
                             int value) const {
        int index = value / bucket_size;
        const int sample_count = static_cast<int>(sample_list.size());
        if (index + 1 >= sample_count) {
            return index;
        }
        if (value - sample_list[index] <= sample_list[index + 1] - value) {
            return index;
        }

        return index + 1;
    }

    static int default_bucket_size(int max_m) {
        assert(max_m >= 0);

        int bucket_size = 1;
        while (static_cast<long long>(bucket_size) * bucket_size <= max_m) {
            bucket_size *= 2;
        }

        return bucket_size;
    }
};