cp-library-cpp

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

View the Project on GitHub NotLeonian/cp-library-cpp

:heavy_check_mark: Binomial Coefficient Sum
(library/math/binomial_coefficient_sum.hpp)

Binomial Coefficient Sum

Depends on

Verified with

Code

#ifndef SUISEN_SUM_OF_BINOM
#define SUISEN_SUM_OF_BINOM

#include "library/algorithm/mo.hpp"
#include "library/math/pow_mods.hpp"
#include "library/math/factorial.hpp"

#include <tuple>

namespace suisen {
    // {n,m} in qs: Sum[i=0,n] r^i Binom[m,i]
    template <typename Mint>
    std::vector<Mint> binom_sum_prefix(const std::vector<std::pair<int, int>>& qs, const Mint& r = 1) {
        const int q = qs.size();
        std::vector<Mint> res(qs.size());

        if (r == -1) {
            factorial<Mint> fac;
            for (int i = 0; i < q; ++i) {
                auto [n, m] = qs[i];
                n = std::min(n, m);
                if (n < 0) {
                    res[i] = 0;
                } else if (m == 0) {
                    // n = m = 0
                    res[i] = 1;
                } else {
                    res[i] = (n % 2 == 1 ? -1 : +1) * fac.binom(m - 1, n);
                }
            }
            return res;
        }

        std::vector<int> qid;
        std::vector<std::pair<int, int>> ranges;
        qid.reserve(q), ranges.reserve(q);
        int max_m = 1;
        for (int i = 0; i < q; ++i) {
            auto [n, m] = qs[i];
            n = std::min(n, m);
            if (n < 0) {
                res[i] = 0;
            } else if (m == 0) {
                // n = m = 0
                res[i] = 1;
            } else {
                qid.push_back(i);
                ranges.emplace_back(n, m);
                max_m = std::max(max_m, m);
            }
        }
        pow_mods<Mint> pow_r(r, max_m + 1);
        factorial<Mint> fac(max_m + 1);
        const Mint inv_r1 = (r + 1).inv();

        Mo mo(max_m, std::move(ranges));
        Mint sum = 1;
        const std::vector<Mint> res_mo = mo.solve(
            [&] {
                return sum;
            },
            [&](const int n) { // Add Left
                const int m = mo.get_right();
                // (n + 1, m) -> (n, m)
                sum -= pow_r[n + 1] * fac.binom(m, n + 1);
            },
            [&](const int n) { // Del Left
                const int m = mo.get_right();
                // (n, m) -> (n + 1, m)
                sum += pow_r[n + 1] * fac.binom(m, n + 1);
            },
            [&](const int m) { // Add Right
                const int n = mo.get_left();
                // (n, m) -> (n, m + 1)
                sum = (r + 1) * sum - pow_r[n + 1] * fac.binom(m, n);
            },
            [&](const int m) { // Del Right
                const int n = mo.get_left();
                // (n, m + 1) -> (n, m)
                sum = (pow_r[n + 1] * fac.binom(m, n) + sum) * inv_r1;
            }
        );
        for (int i = 0; i < int(qid.size()); ++i) res[qid[i]] = res_mo[i];
        return res;
    }

    // {a,b,m} in qs: Sum[i=a,b] r^i Binom[m,i]
    template <typename Mint>
    std::vector<Mint> binom_sum(const std::vector<std::tuple<int, int, int>>& qs, const Mint& r = 1) {
        const int q = qs.size();
        std::vector<std::pair<int, int>> ls(q), rs(q);
        for (int i = 0; i < q; ++i) {
            const auto [a, b, m] = qs[i];
            ls[i] = { a - 1, m };
            rs[i] = { b, m };
        }
        std::vector<Mint> suml = binom_sum_prefix(ls, r);
        std::vector<Mint> sumr = binom_sum_prefix(rs, r);
        std::vector<Mint> sum(q);
        for (int i = 0; i < q; ++i) sum[i] = sumr[i] - suml[i];
        return sum;
    }
} // namespace suisen



#endif // SUISEN_SUM_OF_BINOM
#line 1 "library/math/binomial_coefficient_sum.hpp"



#line 1 "library/algorithm/mo.hpp"



#include <algorithm>
#include <cmath>
#include <numeric>
#include <vector>

namespace suisen {
    struct Mo {
        Mo() = default;
        Mo(const int n, const std::vector<std::pair<int, int>> &queries) : n(n), q(queries.size()), b(bucket_size(n, q)), qs(queries), ord(q) {
            std::iota(ord.begin(), ord.end(), 0);
            std::sort(
                ord.begin(), ord.end(),
                [&, this](int i, int j) {
                    const auto &[li, ri] = qs[i];
                    const auto &[lj, rj] = qs[j];
                    const int bi = li / b, bj = lj / b;
                    if (bi != bj) return bi < bj;
                    if (ri != rj) return bi & 1 ? ri > rj : ri < rj;
                    return li < lj;
                }
            );
        }

        // getter methods used in updating functions: AddL, DelL, etc.
        auto get_left()  const { return l; }
        auto get_right() const { return r; }
        auto get_range() const { return std::make_pair(l, r); }
        auto get_query_id() const { return query_id; }

        /**
         * [Parameters]
         * Eval : () -> T : return the current answer
         * AddL : int -> any (discarded) : add    `l` to   the current range [l + 1, r)
         * DelL : int -> any (discarded) : delete `l` from the current range [l, r)
         * AddR : int -> any (discarded) : add    `r` to   the current range [l, r)
         * DelR : int -> any (discarded) : delete `r` from the current range [l, r + 1)
         * 
         * [Note]
         * starting from the range [0, 0).
         */
        template <typename Eval, typename AddL, typename DelL, typename AddR, typename DelR>
        auto solve(Eval eval, AddL add_l, DelL del_l, AddR add_r, DelR del_r) {
            l = 0, r = 0;
            std::vector<decltype(eval())> res(q);
            for (int qi : ord) {
                const auto &[nl, nr] = qs[query_id = qi];
                while (r < nr) add_r(r), ++r;
                while (l > nl) --l, add_l(l);
                while (r > nr) --r, del_r(r);
                while (l < nl) del_l(l), ++l;
                res[qi] = eval();
            }
            return res;
        }
    
        /**
         * [Parameters]
         * Eval : () -> T : return the current answer
         * Add : int -> any (discarded) : add    `i` to   the current range [i + 1, r) or [l, i)
         * Del : int -> any (discarded) : delete `i` from the current range [i, r) or [l, i + 1)
         * 
         * [Note]
         * starting from the range [0, 0).
         */
        template <typename Eval, typename Add, typename Del>
        auto solve(Eval eval, Add add, Del del) {
            return solve(eval, add, del, add, del);
        }

    private:
        int n, q, b;
        int query_id = -1;
        std::vector<std::pair<int, int>> qs;
        std::vector<int> ord;
        int l = 0, r = 0;

        static int bucket_size(int n, int q) {
            return std::max(1, int(::sqrt(3) * n / ::sqrt(std::max(1, 2 * q))));
        }
    };
} // namespace suisen


#line 1 "library/math/pow_mods.hpp"



#line 5 "library/math/pow_mods.hpp"

namespace suisen {
    template <int base_as_int, typename mint>
    struct static_pow_mods {
        static_pow_mods() = default;
        static_pow_mods(int n) { ensure(n); }
        const mint& operator[](int i) const {
            ensure(i);
            return pows[i];
        }
        static void ensure(int n) {
            int sz = pows.size();
            if (sz > n) return;
            pows.resize(n + 1);
            for (int i = sz; i <= n; ++i) pows[i] = base * pows[i - 1];
        }
    private:
        static inline std::vector<mint> pows { 1 };
        static inline mint base = base_as_int;
        static constexpr int mod = mint::mod();
    };

    template <typename mint>
    struct pow_mods {
        pow_mods() = default;
        pow_mods(mint base, int n) : base(base) { ensure(n); }
        const mint& operator[](int i) const {
            ensure(i);
            return pows[i];
        }
        void ensure(int n) const {
            int sz = pows.size();
            if (sz > n) return;
            pows.resize(n + 1);
            for (int i = sz; i <= n; ++i) pows[i] = base * pows[i - 1];
        }
    private:
        mutable std::vector<mint> pows { 1 };
        mint base;
        static constexpr int mod = mint::mod();
    };
}


#line 1 "library/math/factorial.hpp"



#include <cassert>
#line 6 "library/math/factorial.hpp"

namespace suisen {
    // 引数として与える値に対して、法が十分大きいことを仮定する
    template <typename T, typename U = T>
    struct factorial {
        factorial() = default;
        factorial(int n) { ensure(n); }

        static void ensure(const int n) {
            int sz = _fac.size();
            if (n + 1 <= sz) return;
            int new_size = std::max(n + 1, sz * 2);
            _fac.resize(new_size), _fac_inv.resize(new_size);
            for (int i = sz; i < new_size; ++i) _fac[i] = _fac[i - 1] * i;
            _fac_inv[new_size - 1] = U(1) / _fac[new_size - 1];
            for (int i = new_size - 1; i > sz; --i) _fac_inv[i - 1] = _fac_inv[i] * i;
        }

        T fac(const int i) {
            ensure(i);
            return _fac[i];
        }
        T operator()(int i) {
            return fac(i);
        }
        U fac_inv(const int i) {
            ensure(i);
            return _fac_inv[i];
        }
        // i の逆数
        // i = 0 の場合は assert 違反となる
        U inv(const int i) {
            assert(i > 0);
            ensure(i);
            return _fac_inv[i] * _fac[i - 1];
        }
        U binom(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[r] * _fac_inv[n - r];
        }
        // binom(n, r) の逆数
        // binom(n, r) = 0 の場合は assert 違反となる
        U binom_inv(const int n, const int r) {
            assert(r >= 0 and n >= r);
            ensure(n);
            return _fac_inv[n] * _fac[r] * _fac[n - r];
        }
        // n 種類から重複を許して r 個選ぶ場合の数
        // x_1+x_2+...+x_n=r（x_i は非負整数）となる x の個数でもある
        // multichoose(n, r) = binom(n + r - 1, r)
        U multichoose(const int n, const int r) {
            if (n < 0 or r < 0) return 0;
            return r > 0 ? binom(n + r - 1, r) : U(1);
        }
        // n 種類から重複を許して r 個選ぶ場合の数 multichoose(n, r) の逆数
        // x_1+x_2+...+x_n=r（x_i は非負整数）となる x の個数の逆数でもある
        // multichoose(n, r) = binom(n + r - 1, r)
        // multichoose(n, r) = 0 の場合は assert 違反となる
        U multichoose_inv(const int n, const int r) {
            assert(n >= 0 and r >= 0);
            return r > 0 ? binom_inv(n + r - 1, r) : U(1);
        }
        template <typename ...Ds, std::enable_if_t<std::conjunction_v<std::is_integral<Ds>...>, std::nullptr_t> = nullptr>
        U polynom(const int n, const Ds& ...ds) {
            if (n < 0) return 0;
            ensure(n);
            int sumd = 0;
            U res = _fac[n];
            for (int d : { ds... }) {
                if (d < 0 or d > n) return 0;
                sumd += d;
                res *= _fac_inv[d];
            }
            if (sumd > n) return 0;
            res *= _fac_inv[n - sumd];
            return res;
        }
        U perm(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[n - r];
        }
        // perm(n, r) の逆数
        // perm(n, r) = 0 の場合は assert 違反となる
        U perm_inv(const int n, const int r) {
            assert(r >= 0 and n >= r);
            ensure(n);
            return _fac_inv[n] * _fac[n - r];
        }
    private:
        static std::vector<T> _fac;
        static std::vector<U> _fac_inv;
    };
    template <typename T, typename U>
    std::vector<T> factorial<T, U>::_fac{ 1 };
    template <typename T, typename U>
    std::vector<U> factorial<T, U>::_fac_inv{ 1 };
} // namespace suisen


#line 7 "library/math/binomial_coefficient_sum.hpp"

#include <tuple>

namespace suisen {
    // {n,m} in qs: Sum[i=0,n] r^i Binom[m,i]
    template <typename Mint>
    std::vector<Mint> binom_sum_prefix(const std::vector<std::pair<int, int>>& qs, const Mint& r = 1) {
        const int q = qs.size();
        std::vector<Mint> res(qs.size());

        if (r == -1) {
            factorial<Mint> fac;
            for (int i = 0; i < q; ++i) {
                auto [n, m] = qs[i];
                n = std::min(n, m);
                if (n < 0) {
                    res[i] = 0;
                } else if (m == 0) {
                    // n = m = 0
                    res[i] = 1;
                } else {
                    res[i] = (n % 2 == 1 ? -1 : +1) * fac.binom(m - 1, n);
                }
            }
            return res;
        }

        std::vector<int> qid;
        std::vector<std::pair<int, int>> ranges;
        qid.reserve(q), ranges.reserve(q);
        int max_m = 1;
        for (int i = 0; i < q; ++i) {
            auto [n, m] = qs[i];
            n = std::min(n, m);
            if (n < 0) {
                res[i] = 0;
            } else if (m == 0) {
                // n = m = 0
                res[i] = 1;
            } else {
                qid.push_back(i);
                ranges.emplace_back(n, m);
                max_m = std::max(max_m, m);
            }
        }
        pow_mods<Mint> pow_r(r, max_m + 1);
        factorial<Mint> fac(max_m + 1);
        const Mint inv_r1 = (r + 1).inv();

        Mo mo(max_m, std::move(ranges));
        Mint sum = 1;
        const std::vector<Mint> res_mo = mo.solve(
            [&] {
                return sum;
            },
            [&](const int n) { // Add Left
                const int m = mo.get_right();
                // (n + 1, m) -> (n, m)
                sum -= pow_r[n + 1] * fac.binom(m, n + 1);
            },
            [&](const int n) { // Del Left
                const int m = mo.get_right();
                // (n, m) -> (n + 1, m)
                sum += pow_r[n + 1] * fac.binom(m, n + 1);
            },
            [&](const int m) { // Add Right
                const int n = mo.get_left();
                // (n, m) -> (n, m + 1)
                sum = (r + 1) * sum - pow_r[n + 1] * fac.binom(m, n);
            },
            [&](const int m) { // Del Right
                const int n = mo.get_left();
                // (n, m + 1) -> (n, m)
                sum = (pow_r[n + 1] * fac.binom(m, n) + sum) * inv_r1;
            }
        );
        for (int i = 0; i < int(qid.size()); ++i) res[qid[i]] = res_mo[i];
        return res;
    }

    // {a,b,m} in qs: Sum[i=a,b] r^i Binom[m,i]
    template <typename Mint>
    std::vector<Mint> binom_sum(const std::vector<std::tuple<int, int, int>>& qs, const Mint& r = 1) {
        const int q = qs.size();
        std::vector<std::pair<int, int>> ls(q), rs(q);
        for (int i = 0; i < q; ++i) {
            const auto [a, b, m] = qs[i];
            ls[i] = { a - 1, m };
            rs[i] = { b, m };
        }
        std::vector<Mint> suml = binom_sum_prefix(ls, r);
        std::vector<Mint> sumr = binom_sum_prefix(rs, r);
        std::vector<Mint> sum(q);
        for (int i = 0; i < q; ++i) sum[i] = sumr[i] - suml[i];
        return sum;
    }
} // namespace suisen
Back to top page