cp-algo/math/combinatorics.hpp

• View this file on GitHub
• Last update: 2025-05-27 12:29:19+02:00
• Include: #include "cp-algo/math/combinatorics.hpp"

Depends on

• cp-algo/math/common.hpp

Required by

• cp-algo/linalg/frobenius.hpp
• cp-algo/math/factorials.hpp
• cp-algo/math/poly.hpp

Verified with

• Binomial Coefficient (Prime Mod) (verify/combi/binom.test.cpp)
• Many Factorials (verify/combi/many_facts.test.cpp)
• Many Factorials (dynamic Mod) (verify/combi/many_facts_dynamic.test.cpp)
• Characteristic Polynomial (verify/linalg/characteristic.test.cpp)
• Counting Eulerian Circuits (verify/linalg/euler_circs.test.cpp)
• Pow of Matrix (Frobenius) (verify/linalg/pow_fast.test.cpp)
• Counting Spanning Trees (Directed) (verify/linalg/spanning_directed.test.cpp)
• Counting Spanning Trees (Undirected) (verify/linalg/spanning_undirected.test.cpp)
• Bell Number (verify/poly/bell.test.cpp)
• Division of Polynomials (verify/poly/div.test.cpp)
• Exp of Power Series (verify/poly/exp.test.cpp)
• Find Linear Recurrence (verify/poly/find_linrec.test.cpp)
• Polynomial Interpolation (verify/poly/inter.test.cpp)
• Polynomial Interpolation (Geometric Sequence) (verify/poly/inter_geo.test.cpp)
• Inv of Power Series (verify/poly/inv.test.cpp)
• Inv of Polynomials (verify/poly/inv_mod.test.cpp)
• Consecutive Terms of Linear Recursion (verify/poly/linrec_consecutive.test.cpp)
• Kth term of Linearly Recurrent Sequence (verify/poly/linrec_single.test.cpp)
• Log of Power Series (verify/poly/log.test.cpp)
• Multipoint Evaluation (verify/poly/multieval.test.cpp)
• Multipoint Evaluation (Geometric Sequence) (verify/poly/multieval_geo.test.cpp)
• Conversion to Newton Basis (verify/poly/newton.test.cpp)
• Pow of Power Series (verify/poly/pow.test.cpp)
• Product of Polynomial Sequence (verify/poly/prod_sequence.test.cpp)
• Polynomial Root Finding (verify/poly/roots.test.cpp)
• Shift of Sampling Points of Polynomial (verify/poly/sampling.test.cpp)
• Sqrt of Power Series (verify/poly/sqrt.test.cpp)
• $\sum\limits_{i=0}^\infty r^i i^d$ (verify/poly/sum_exp_poly.test.cpp)
• Polynomial Taylor Shift (verify/poly/taylor.test.cpp)

Code

#ifndef CP_ALGO_MATH_COMBINATORICS_HPP
#define CP_ALGO_MATH_COMBINATORICS_HPP
#include "../math/common.hpp"
#include <cassert>
#include <ranges>
namespace cp_algo::math {
    // fact/rfact/small_inv are caching
    // Beware of usage with dynamic mod
    template<typename T>
    T fact(auto n) {
        static std::vector<T> F(maxn);
        static bool init = false;
        if(!init) {
            F[0] = T(1);
            for(int i = 1; i < maxn; i++) {
                F[i] = F[i - 1] * T(i);
            }
            init = true;
        }
        return F[n];
    }
    // Only works for modint types
    template<typename T>
    T rfact(auto n) {
        static std::vector<T> F(maxn);
        static bool init = false;
        if(!init) {
            int t = (int)std::min<int64_t>(T::mod(), maxn) - 1;
            F[t] = T(1) / fact<T>(t);
            for(int i = t - 1; i >= 0; i--) {
                F[i] = F[i + 1] * T(i + 1);
            }
            init = true;
        }
        return F[n];
    }
    template<typename T, int base>
    T pow_fixed(int n) {
        static std::vector<T> prec_low(1 << 16);
        static std::vector<T> prec_high(1 << 16);
        static bool init = false;
        if(!init) {
            init = true;
            prec_low[0] = prec_high[0] = T(1);
            T step_low = T(base);
            T step_high = bpow(T(base), 1 << 16);
            for(int i = 1; i < (1 << 16); i++) {
                prec_low[i] = prec_low[i - 1] * step_low;
                prec_high[i] = prec_high[i - 1] * step_high;
            }
        }
        return prec_low[n & 0xFFFF] * prec_high[n >> 16];
    }
    template<typename T>
    std::vector<T> bulk_invs(auto const& args) {
        std::vector<T> res(std::size(args), args[0]);
        for(size_t i = 1; i < std::size(args); i++) {
            res[i] = res[i - 1] * args[i];
        }
        auto all_invs = T(1) / res.back();
        for(size_t i = std::size(args) - 1; i > 0; i--) {
            res[i] = all_invs * res[i - 1];
            all_invs *= args[i];
        }
        res[0] = all_invs;
        return res;
    }
    template<typename T>
    T small_inv(auto n) {
        static auto F = bulk_invs<T>(std::views::iota(1, maxn));
        return F[n - 1];
    }
    template<typename T>
    T binom_large(T n, auto r) {
        assert(r < maxn);
        T ans = 1;
        for(decltype(r) i = 0; i < r; i++) {
            ans = ans * T(n - i) * small_inv<T>(i + 1);
        }
        return ans;
    }
    template<typename T>
    T binom(auto n, auto r) {
        if(r < 0 || r > n) {
            return T(0);
        } else if(n >= maxn) {
            return binom_large(T(n), r);
        } else {
            return fact<T>(n) * rfact<T>(r) * rfact<T>(n - r);
        }
    }
}
#endif // CP_ALGO_MATH_COMBINATORICS_HPP

#line 1 "cp-algo/math/combinatorics.hpp"


#line 1 "cp-algo/math/common.hpp"


#include <functional>
#include <cstdint>
#include <cassert>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
    const int maxn = CP_ALGO_MAXN;
#else
    const int maxn = 1 << 19;
#endif
    const int magic = 64; // threshold for sizes to run the naive algo

    auto bpow(auto const& x, auto n, auto const& one, auto op) {
        if(n == 0) {
            return one;
        } else {
            auto t = bpow(x, n / 2, one, op);
            t = op(t, t);
            if(n % 2) {
                t = op(t, x);
            }
            return t;
        }
    }
    auto bpow(auto x, auto n, auto ans) {
        return bpow(x, n, ans, std::multiplies{});
    }
    template<typename T>
    T bpow(T const& x, auto n) {
        return bpow(x, n, T(1));
    }
    inline constexpr auto inv2(auto x) {
        assert(x % 2);
        std::make_unsigned_t<decltype(x)> y = 1;
        while(y * x != 1) {
            y *= 2 - x * y;
        }
        return y;
    }
}

#line 5 "cp-algo/math/combinatorics.hpp"
#include <ranges>
namespace cp_algo::math {
    // fact/rfact/small_inv are caching
    // Beware of usage with dynamic mod
    template<typename T>
    T fact(auto n) {
        static std::vector<T> F(maxn);
        static bool init = false;
        if(!init) {
            F[0] = T(1);
            for(int i = 1; i < maxn; i++) {
                F[i] = F[i - 1] * T(i);
            }
            init = true;
        }
        return F[n];
    }
    // Only works for modint types
    template<typename T>
    T rfact(auto n) {
        static std::vector<T> F(maxn);
        static bool init = false;
        if(!init) {
            int t = (int)std::min<int64_t>(T::mod(), maxn) - 1;
            F[t] = T(1) / fact<T>(t);
            for(int i = t - 1; i >= 0; i--) {
                F[i] = F[i + 1] * T(i + 1);
            }
            init = true;
        }
        return F[n];
    }
    template<typename T, int base>
    T pow_fixed(int n) {
        static std::vector<T> prec_low(1 << 16);
        static std::vector<T> prec_high(1 << 16);
        static bool init = false;
        if(!init) {
            init = true;
            prec_low[0] = prec_high[0] = T(1);
            T step_low = T(base);
            T step_high = bpow(T(base), 1 << 16);
            for(int i = 1; i < (1 << 16); i++) {
                prec_low[i] = prec_low[i - 1] * step_low;
                prec_high[i] = prec_high[i - 1] * step_high;
            }
        }
        return prec_low[n & 0xFFFF] * prec_high[n >> 16];
    }
    template<typename T>
    std::vector<T> bulk_invs(auto const& args) {
        std::vector<T> res(std::size(args), args[0]);
        for(size_t i = 1; i < std::size(args); i++) {
            res[i] = res[i - 1] * args[i];
        }
        auto all_invs = T(1) / res.back();
        for(size_t i = std::size(args) - 1; i > 0; i--) {
            res[i] = all_invs * res[i - 1];
            all_invs *= args[i];
        }
        res[0] = all_invs;
        return res;
    }
    template<typename T>
    T small_inv(auto n) {
        static auto F = bulk_invs<T>(std::views::iota(1, maxn));
        return F[n - 1];
    }
    template<typename T>
    T binom_large(T n, auto r) {
        assert(r < maxn);
        T ans = 1;
        for(decltype(r) i = 0; i < r; i++) {
            ans = ans * T(n - i) * small_inv<T>(i + 1);
        }
        return ans;
    }
    template<typename T>
    T binom(auto n, auto r) {
        if(r < 0 || r > n) {
            return T(0);
        } else if(n >= maxn) {
            return binom_large(T(n), r);
        } else {
            return fact<T>(n) * rfact<T>(r) * rfact<T>(n - r);
        }
    }
}

Back to top page