cp-algo/math/combinatorics.hpp

• View this file on GitHub
• Last update: 2026-10-16 14:16:08+02: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/fps.hpp
• cp-algo/math/laurent.hpp
• cp-algo/math/poly.hpp
• cp-algo/math/poly/calculus.hpp
• cp-algo/math/poly/chirpz.hpp
• cp-algo/math/poly/compose.hpp
• cp-algo/math/poly/eval.hpp
• cp-algo/math/poly/powmod.hpp
• cp-algo/math/poly/series.hpp
• cp-algo/math/poly/series/exp.hpp
• cp-algo/math/poly/series/log.hpp
• cp-algo/math/poly/series/pow.hpp
• cp-algo/math/poly/sparse.hpp
• cp-algo/math/poly/sparse/exp.hpp
• cp-algo/math/poly/sparse/impl.hpp
• cp-algo/math/poly/sparse/inv.hpp
• cp-algo/math/poly/sparse/log.hpp
• cp-algo/math/poly/sparse/pow.hpp
• cp-algo/math/poly/sparse/sqrt.hpp
• cp-algo/math/poly/transform.hpp
• tests/fps.cpp
• tests/frobenius.cpp
• tests/poly.cpp
• tests/poly_series.cpp
• tests/poly_sparse.cpp
• tests/power.cpp

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)
• Lazy FPS: Exp of Power Series (verify/fps/exp.test.cpp)
• Lazy FPS: Inv of Power Series (verify/fps/inv.test.cpp)
• Lazy FPS: Log of Power Series (verify/fps/log.test.cpp)
• Lazy FPS: Product with a Known Polynomial (verify/fps/mul.test.cpp)
• Lazy FPS: Pow of Power Series (verify/fps/pow.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 Number (verify/math/chromatic_number.test.cpp)
• Chromatic Polynomial (verify/math/chromatic_poly.test.cpp)
• Bell Number (verify/poly/bell.test.cpp)
• Bernoulli Number (verify/poly/bernoulli.test.cpp)
• Composition of Formal Power Series (verify/poly/compose.test.cpp)
• Exp of Power Series (verify/poly/exp.test.cpp)
• Polynomial Interpolation (verify/poly/inter.test.cpp)
• Polynomial Interpolation (Geometric Sequence) (verify/poly/inter_geo.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)
• Prefix Sum of Polynomial (verify/poly/prefix_sum.test.cpp)
• Polynomial Root Finding (verify/poly/roots.test.cpp)
• Shift of Sampling Points of Polynomial (verify/poly/sampling.test.cpp)
• Exp of Power Series (Sparse) (verify/poly/sparse/exp.test.cpp)
• Inv of Power Series (Sparse) (verify/poly/sparse/inv.test.cpp)
• Log of Power Series (Sparse) (verify/poly/sparse/log.test.cpp)
• Pow of Power Series (Sparse) (verify/poly/sparse/pow.test.cpp)
• Sqrt of Power Series (Sparse) (verify/poly/sparse/sqrt.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 const auto F = [] {
            big_vector<T> F(maxn);
            F[0] = T(1);
            for(int i = 1; i < maxn; i++) {
                F[i] = F[i - 1] * T(i);
            }
            return F;
        }();
        return F[n];
    }
    // Only works for modint types
    template<typename T>
    T rfact(auto n) {
        static const auto F = [] {
            big_vector<T> F(maxn);
            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);
            }
            return F;
        }();
        return F[n];
    }
    template<typename T, int base>
    T pow_fixed(int n) {
        static const auto prec = [] {
            big_vector<T> prec_low(1 << 16), prec_high(1 << 16);
            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 std::pair{std::move(prec_low), std::move(prec_high)};
        }();
        return prec.first[n & 0xFFFF] * prec.second[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>
#include <vector>
#include <algorithm>
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

    // Nonnegative 64-bit exponents, with an associative operation and its identity.
    // Windows >1 precompute odd powers only when that saves operations.
    template<int window = 1>
    auto bpow(auto const& x, auto n, auto const& one, auto op) {
        static_assert(window >= 1 && window <= 6);
        if constexpr(window > 1) {
            if(n == 0) {return one;}
            int bits = std::bit_width(uint64_t(n));
            auto low_bit = [&](int high) {
                int low = std::max(0, high - window + 1);
                while(!((n >> low) & 1)) {low++;}
                return low;
            };
            int first = low_bit(bits - 1);
            int cost = (1 << (window - 1)) + first;
            for(int j = first - 1; j >= 0;) {
                if(!((n >> j) & 1)) {j--;}
                else {cost++; j = low_bit(j) - 1;}
            }
            // Do not pay for the table when binary powering uses fewer operations.
            if(cost >= bits + std::popcount(uint64_t(n)) - 2) {return bpow<1>(x, n, one, op);}
            using T = std::decay_t<decltype(x)>;
            std::vector<T> odd;
            odd.reserve(1 << (window - 1));
            odd.push_back(x);
            auto square = op(x, x);
            while(odd.size() < size_t(1 << (window - 1))) {odd.push_back(op(odd.back(), square));}
            auto ans = odd[(n >> first) / 2];
            for(int j = first - 1; j >= 0;) {
                if(!((n >> j) & 1)) {ans = op(ans, ans); j--;}
                else {
                    int low = low_bit(j), length = j - low + 1;
                    auto digit = (n >> low) & ((1u << length) - 1);
                    for(int i = 0; i < length; i++) {ans = op(ans, ans);}
                    ans = op(ans, odd[digit / 2]);
                    j = low - 1;
                }
            }
            return ans;
        } else {
            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;
        }
    }
    template<int window = 1>
    auto bpow(auto x, auto n, auto ans) {
        return bpow<window>(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>
#line 10 "cp-algo/util/big_alloc.hpp"
#include <string>
#include <cstddef>
#line 13 "cp-algo/util/big_alloc.hpp"
#include <iostream>
#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::has_single_bit(Align), "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;

        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);
                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>>;
}


#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 const auto F = [] {
            big_vector<T> F(maxn);
            F[0] = T(1);
            for(int i = 1; i < maxn; i++) {
                F[i] = F[i - 1] * T(i);
            }
            return F;
        }();
        return F[n];
    }
    // Only works for modint types
    template<typename T>
    T rfact(auto n) {
        static const auto F = [] {
            big_vector<T> F(maxn);
            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);
            }
            return F;
        }();
        return F[n];
    }
    template<typename T, int base>
    T pow_fixed(int n) {
        static const auto prec = [] {
            big_vector<T> prec_low(1 << 16), prec_high(1 << 16);
            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 std::pair{std::move(prec_low), std::move(prec_high)};
        }();
        return prec.first[n & 0xFFFF] * prec.second[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 const auto F= []{ big_vector<T>F(maxn); F[0] =T(1); for(int i=1;i<maxn;i++){ F[i] =F[i-1] *T(i); } return F; }(); return F[n]; } template<typename T> T rfact(auto n){ static const auto F= []{ big_vector<T>F(maxn); 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); } return F; }(); return F[n]; } template<typename T,int base> T pow_fixed(int n){ static const auto prec= []{ big_vector<T>prec_low(1<<16),prec_high(1<<16); 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 std::pair{std::move(prec_low),std::move(prec_high)}; }(); return prec.first[n&0xFFFF] *prec.second[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>
#include <vector>
#include <algorithm>
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; template<int window=1> auto bpow(auto const&x,auto n,auto const&one,auto op){ static_assert(window>=1&&window<=6); if constexpr(window>1){ if(n==0){return one;} int bits=std::bit_width(uint64_t(n)); auto low_bit= [&](int high){ int low=std::max(0,high-window+1); while(!((n>>low)&1)){low++;} return low; }; int first=low_bit(bits-1); int cost=(1<<(window-1))+first; for(int j=first-1;j>=0;){ if(!((n>>j)&1)){j--;} else{cost++;j=low_bit(j)-1;} } if(cost>=bits+std::popcount(uint64_t(n))-2){return bpow<1>(x,n,one,op);} using T=std::decay_t<decltype(x)>; std::vector<T>odd; odd.reserve(1<<(window-1)); odd.push_back(x); auto square=op(x,x); while(odd.size()<size_t(1<<(window-1))){odd.push_back(op(odd.back(),square));} auto ans=odd[(n>>first)/2]; for(int j=first-1;j>=0;){ if(!((n>>j)&1)){ans=op(ans,ans);j--;} else{ int low=low_bit(j),length=j-low+1; auto digit=(n>>low)&((1u<<length)-1); for(int i=0;i<length;i++){ans=op(ans,ans);} ans=op(ans,odd[digit/2]); j=low-1; } } return ans; }else{ 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; } } template<int window=1> auto bpow(auto x,auto n,auto ans){ return bpow<window>(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>
#line 10 "cp-algo/util/big_alloc.hpp"
#include <string>
#include <cstddef>
#line 13 "cp-algo/util/big_alloc.hpp"
#include <iostream>
#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::has_single_bit(Align),"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; 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); 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>>; }
#line 6 "cp-algo/math/combinatorics.hpp"
#include <ranges>
namespace cp_algo::math{ template<typename T> T fact(auto n){ static const auto F= []{ big_vector<T>F(maxn); F[0] =T(1); for(int i=1;i<maxn;i++){ F[i] =F[i-1] *T(i); } return F; }(); return F[n]; } template<typename T> T rfact(auto n){ static const auto F= []{ big_vector<T>F(maxn); 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); } return F; }(); return F[n]; } template<typename T,int base> T pow_fixed(int n){ static const auto prec= []{ big_vector<T>prec_low(1<<16),prec_high(1<<16); 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 std::pair{std::move(prec_low),std::move(prec_high)}; }(); return prec.first[n&0xFFFF] *prec.second[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