cp-algo/math/combinatorics.hpp

• View this file on GitHub
• Last update: 2025-12-30 00:27:16+01:00
• Include: #include "cp-algo/math/combinatorics.hpp"

Depends on

• cp-algo/math/common.hpp
• cp-algo/util/big_alloc.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)
• Chromatic Polynomial (verify/math/chromatic_poly.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 "../util/big_alloc.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 big_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 big_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 big_vector<T> prec_low(1 << 16);
        static big_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>
    big_vector<T> bulk_invs(auto const& args) {
        big_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>
#include <bit>
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;
        }
        auto ans = x;
        for(int j = std::bit_width<uint64_t>(n) - 2; ~j; j--) {
            ans = op(ans, ans);
            if((n >> j) & 1) {
                ans = op(ans, x);
            }
        }
        return ans;
    }
    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 1 "cp-algo/util/big_alloc.hpp"



#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#include <vector>
#include <string>
#include <cstddef>
#include <iostream>
#include <generator>
#include <forward_list>

// Single macro to detect POSIX platforms (Linux, Unix, macOS)
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
#  define CP_ALGO_USE_MMAP 1
#  include <sys/mman.h>
#else
#  define CP_ALGO_USE_MMAP 0
#endif

namespace cp_algo {
    template <typename T, size_t Align = 32>
    class big_alloc {
        static_assert( Align >= alignof(void*), "Align must be at least pointer-size");
        static_assert(std::popcount(Align) == 1, "Align must be a power of two");
    public:
        using value_type = T;
        template <class U> struct rebind { using other = big_alloc<U, Align>; };
        constexpr bool operator==(const big_alloc&) const = default;
        constexpr bool operator!=(const big_alloc&) const = default;

        big_alloc() noexcept = default;
        template <typename U, std::size_t A>
        big_alloc(const big_alloc<U, A>&) noexcept {}

        [[nodiscard]] T* allocate(std::size_t n) {
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align = std::max<std::size_t>(alignof(T),  Align);
#if CP_ALGO_USE_MMAP
            if (padded >= MEGABYTE) {
                void* raw = mmap(nullptr, padded,
                                PROT_READ | PROT_WRITE,
                                MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
                madvise(raw, padded, MADV_HUGEPAGE);
                madvise(raw, padded, MADV_POPULATE_WRITE);
                return static_cast<T*>(raw);
            }
#endif
            return static_cast<T*>(::operator new(padded, std::align_val_t(align)));
        }

        void deallocate(T* p, std::size_t n) noexcept {
            if (!p) return;
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align  = std::max<std::size_t>(alignof(T),  Align);
    #if CP_ALGO_USE_MMAP
            if (padded >= MEGABYTE) { munmap(p, padded); return; }
    #endif
            ::operator delete(p, padded, std::align_val_t(align));
        }

    private:
        static constexpr std::size_t MEGABYTE = 1 << 20;
        static constexpr std::size_t round_up(std::size_t x) noexcept {
            return (x + Align - 1) / Align * Align;
        }
    };

    template<typename T> using big_vector = std::vector<T, big_alloc<T>>;
    template<typename T> using big_basic_string = std::basic_string<T, std::char_traits<T>, big_alloc<T>>;
    template<typename T> using big_deque = std::deque<T, big_alloc<T>>;
    template<typename T> using big_stack = std::stack<T, big_deque<T>>;
    template<typename T> using big_queue = std::queue<T, big_deque<T>>;
    template<typename T> using big_priority_queue = std::priority_queue<T, big_vector<T>>;
    template<typename T> using big_forward_list = std::forward_list<T, big_alloc<T>>;
    using big_string = big_basic_string<char>;

    template<typename Key, typename Value, typename Compare = std::less<Key>>
    using big_map = std::map<Key, Value, Compare, big_alloc<std::pair<const Key, Value>>>;
    template<typename T, typename Compare = std::less<T>>
    using big_multiset = std::multiset<T, Compare, big_alloc<T>>;
    template<typename T, typename Compare = std::less<T>>
    using big_set = std::set<T, Compare, big_alloc<T>>;
    template<typename Ref, typename V = void>

    using big_generator = std::generator<Ref, V, big_alloc<std::byte>>;
}

// Deduction guide to make elements_of with big_generator default to big_alloc
namespace std::ranges {
    template<typename Ref, typename V>
    elements_of(cp_algo::big_generator<Ref, V>&&) -> elements_of<cp_algo::big_generator<Ref, V>&&, cp_algo::big_alloc<std::byte>>;
}


#line 6 "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 big_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 big_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 big_vector<T> prec_low(1 << 16);
        static big_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>
    big_vector<T> bulk_invs(auto const& args) {
        big_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);
        }
    }
}

#ifndef CP_ALGO_MATH_COMBINATORICS_HPP
#define CP_ALGO_MATH_COMBINATORICS_HPP
#include "../math/common.hpp"
#include "../util/big_alloc.hpp"
#include <cassert>
#include <ranges>
namespace cp_algo::math{template<typename T>T fact(auto n){static big_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];}template<typename T>T rfact(auto n){static big_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 big_vector<T>prec_low(1<<16);static big_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>big_vector<T>bulk_invs(auto const&args){big_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

#line 1 "cp-algo/math/combinatorics.hpp"
#line 1 "cp-algo/math/common.hpp"
#include <functional>
#include <cstdint>
#include <cassert>
#include <bit>
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;auto bpow(auto const&x,auto n,auto const&one,auto op){if(n==0){return one;}auto ans=x;for(int j=std::bit_width<uint64_t>(n)-2;~j;j--){ans=op(ans,ans);if((n>>j)&1){ans=op(ans,x);}}return ans;}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 1 "cp-algo/util/big_alloc.hpp"
#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#include <vector>
#include <string>
#include <cstddef>
#include <iostream>
#include <generator>
#include <forward_list>
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
#  define CP_ALGO_USE_MMAP 1
#  include <sys/mman.h>
#else
#  define CP_ALGO_USE_MMAP 0
#endif
namespace cp_algo{template<typename T,size_t Align=32>class big_alloc{static_assert(Align>=alignof(void*),"Align must be at least pointer-size");static_assert(std::popcount(Align)==1,"Align must be a power of two");public:using value_type=T;template<class U>struct rebind{using other=big_alloc<U,Align>;};constexpr bool operator==(const big_alloc&)const=default;constexpr bool operator!=(const big_alloc&)const=default;big_alloc()noexcept=default;template<typename U,std::size_t A>big_alloc(const big_alloc<U,A>&)noexcept{}[[nodiscard]]T*allocate(std::size_t n){std::size_t padded=round_up(n*sizeof(T));std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){void*raw=mmap(nullptr,padded,PROT_READ|PROT_WRITE,MAP_PRIVATE|MAP_ANONYMOUS,-1,0);madvise(raw,padded,MADV_HUGEPAGE);madvise(raw,padded,MADV_POPULATE_WRITE);return static_cast<T*>(raw);}
#endif
return static_cast<T*>(::operator new(padded,std::align_val_t(align)));}void deallocate(T*p,std::size_t n)noexcept{if(!p)return;std::size_t padded=round_up(n*sizeof(T));std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){munmap(p,padded);return;}
#endif
::operator delete(p,padded,std::align_val_t(align));}private:static constexpr std::size_t MEGABYTE=1<<20;static constexpr std::size_t round_up(std::size_t x)noexcept{return(x+Align-1)/Align*Align;}};template<typename T>using big_vector=std::vector<T,big_alloc<T>>;template<typename T>using big_basic_string=std::basic_string<T,std::char_traits<T>,big_alloc<T>>;template<typename T>using big_deque=std::deque<T,big_alloc<T>>;template<typename T>using big_stack=std::stack<T,big_deque<T>>;template<typename T>using big_queue=std::queue<T,big_deque<T>>;template<typename T>using big_priority_queue=std::priority_queue<T,big_vector<T>>;template<typename T>using big_forward_list=std::forward_list<T,big_alloc<T>>;using big_string=big_basic_string<char>;template<typename Key,typename Value,typename Compare=std::less<Key>>using big_map=std::map<Key,Value,Compare,big_alloc<std::pair<const Key,Value>>>;template<typename T,typename Compare=std::less<T>>using big_multiset=std::multiset<T,Compare,big_alloc<T>>;template<typename T,typename Compare=std::less<T>>using big_set=std::set<T,Compare,big_alloc<T>>;template<typename Ref,typename V=void>using big_generator=std::generator<Ref,V,big_alloc<std::byte>>;}namespace std::ranges{template<typename Ref,typename V>elements_of(cp_algo::big_generator<Ref,V>&&)->elements_of<cp_algo::big_generator<Ref,V>&&,cp_algo::big_alloc<std::byte>>;}
#line 6 "cp-algo/math/combinatorics.hpp"
#include <ranges>
namespace cp_algo::math{template<typename T>T fact(auto n){static big_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];}template<typename T>T rfact(auto n){static big_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 big_vector<T>prec_low(1<<16);static big_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>big_vector<T>bulk_invs(auto const&args){big_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