diff --git a/config.yml b/config.yml index 74d5a5f5c..7185ecb10 100644 --- a/config.yml +++ b/config.yml @@ -740,6 +740,9 @@ notebook: - ball_box_iil: 球与盒(球相同,盒相同,至少装一个球) code_ext: hpp test_ext: cpp + - factl_helper: 大数阶乘与阶乘逆辅助类 + code_ext: hpp + test_ext: cpp nt: - euler_phi_u32: Euler 函数(32 位) code_ext: hpp diff --git a/src/code/comb/binom.hpp b/src/code/comb/binom.hpp index 46dccb1e3..5a0ed007b 100644 --- a/src/code/comb/binom.hpp +++ b/src/code/comb/binom.hpp @@ -6,13 +6,13 @@ namespace tifa_libs::math { -template class fact = fact_helper> -struct binom : fact { - static CEXP auto mod() NE { return fact::mod(); } - CEXPE binom(u32 max_m = fact::DEFUALT_MAX) NE { this->reset(max_m + 1); } +template > +requires std::same_as +struct binom { + CEXPE binom(u32 max_m = fact::DEFUALT_MAX) NE { fact::ensure(max_m + 1); } // $\binom{m}{n}$ - CEXP mint mCn(uint_c auto m, uint_c auto n) CNE { return m < n ? 0 : mPn(m, n) * this->get_ifact(n); } + CEXP mint mCn(uint_c auto m, uint_c auto n) CNE { return m < n ? 0 : mPn(m, n) * fact::get_ifact(n); } // $\binom{m}{n}$ template CEXP mint mCn(T m, T n) CNE { return m < n || n < 0 ? 0 : mCn(to_uint_t(m), to_uint_t(n)); } @@ -20,11 +20,11 @@ struct binom : fact { template CEXP mint lucas(T m, T n) CNE { assert(mint::mod() > 1); - auto f = [this](auto &&f, auto m, auto n) NE -> mint { return n == 0 ? 1 : this->mCn(m % mod(), n % mod()) * f(f, m / mod(), n / mod()); }; + auto f = [this](auto &&f, auto m, auto n) NE -> mint { return n == 0 ? 1 : this->mCn(m % fact::mod(), n % fact::mod()) * f(f, m / fact::mod(), n / fact::mod()); }; return m < n || n < 0 ? 0 : f(f, to_uint_t(m), to_uint_t(n)); } // $\binom{m}{n} \cdot n!$ - CEXP mint mPn(uint_c auto m, uint_c auto n) CNE { return m < n ? 0 : this->get_fact(m) * this->get_ifact(m - n); } + CEXP mint mPn(uint_c auto m, uint_c auto n) CNE { return m < n ? 0 : fact::get_fact(m) * fact::get_ifact(m - n); } // $\binom{m}{n} \cdot n!$ template CEXP mint mPn(T m, T n) CNE { return m < n || n < 0 ? 0 : mPn(to_uint_t(m), to_uint_t(n)); } diff --git a/src/code/comb/fact_helper.hpp b/src/code/comb/fact_helper.hpp index f076e87ff..cde79a6b7 100644 --- a/src/code/comb/fact_helper.hpp +++ b/src/code/comb/fact_helper.hpp @@ -1,37 +1,46 @@ #ifndef TIFALIBS_COMB_FACT_HELPER #define TIFALIBS_COMB_FACT_HELPER -#include "gen_fact.hpp" -#include "gen_ifact.hpp" +#include "../util/traits.hpp" namespace tifa_libs::math { -template +template struct fact_helper { + using val_t = mint; static CEXP u32 DEFUALT_MAX = 10'000'001; - static CEXP u64 mod() NE { return mint::mod(); } - static inline vec fact, ifact; - - CEXP void reset(u32 n = DEFUALT_MAX) NE { - if (n = std::min((u32)mod(), n); n < fact.size()) { - fact.resize(n), ifact.resize(n); - return; - } - fact = gen_fact(n), ifact = gen_ifact(n); + static CEXP u64 mod() NE { return val_t::mod(); } + static inline vec fact, ifact; + + fact_helper() = delete; + + // ensure fact.size() >= sz + static CEXP void ensure(u32 sz = DEFUALT_MAX) NE { + if (sz = std::max(2_u32, std::min((u32)mod(), sz)); sz <= fact.size()) return; + u32 pre = (u32)fact.size(); + fact.resize(sz), ifact.resize(sz); + if (pre < 2) pre = 2, fact[0] = fact[1] = ifact[0] = ifact[1] = 1; + flt_ (u32, i, pre, sz) fact[i] = fact[i - 1] * i; + ifact.back() = fact.back().inv(); + for (u32 i = sz - 1; i > pre; --i) ifact[i - 1] = ifact[i] * i; } - CEXP mint get_fact(u64 n) CNE { + static CEXP val_t get_fact(u64 n) NE { if (n >= mod()) [[unlikely]] return 0; + if (fact.empty()) [[unlikely]] + ensure(); if (n < fact.size()) [[likely]] return fact[n]; - mint _ = fact.back() * n; + val_t _ = fact.back() * n; flt_ (u64, i, fact.size(), n) _ *= i; return _; } - CEXP mint get_ifact(u64 n) CNE { + static CEXP val_t get_ifact(u64 n) NE { if (n >= mod()) [[unlikely]] return 0; + if (fact.empty()) [[unlikely]] + ensure(); if (n < ifact.size()) [[likely]] return ifact[n]; return get_fact(n).inv(); diff --git a/src/code/comb/factl_helper.hpp b/src/code/comb/factl_helper.hpp new file mode 100644 index 000000000..10ced5ea9 --- /dev/null +++ b/src/code/comb/factl_helper.hpp @@ -0,0 +1,82 @@ +#ifndef TIFALIBS_COMB_FACTL_HELPER +#define TIFALIBS_COMB_FACTL_HELPER + +#include "../poly/ctsh_fps.hpp" +#include "fact_helper.hpp" + +namespace tifa_libs::math { + +template +struct factl_helper { + using val_t = TPN poly_t::val_t; + using base_t = fact_helper; + static CEXP u32 threshold = 2'000'001; + + factl_helper() = delete; + + static CEXP val_t get_fact(u64 n) NE { + if (n >= base_t::mod()) [[unlikely]] + return 0; + if (base_t::fact.empty()) [[unlikely]] + base_t::ensure(threshold); + if (n < threshold) return base_t::fact[n]; + return fact_(n); + } + static CEXP val_t get_ifact(u64 n) NE { + if (n >= base_t::mod()) [[unlikely]] + return 0; + if (base_t::fact.empty()) [[unlikely]] + base_t::ensure(threshold); + if (n < threshold) return base_t::ifact[n]; + return fact_(n).inv(); + } + + private: + static CEXP u32 LBSZ = 9; + static CEXP u32 SZ = 1 << LBSZ; + static inline u64 B = base_t::mod() >> LBSZ; + // f[i] = (i*B + 1) * ... * (i*B + B) + static inline poly_t f{1}; + static val_t fact_(u64 n) NE { + static bool inited = false; + if (!inited) { + inited = true; + f.reserve(SZ); + flt_ (u32, i, 0, LBSZ) { + poly_t g = ctsh_fps(f, val_t(1 << i), base_t::ifact, 3_u32 << i); + const auto get = [&](u32 j) { return j < (1_u32 << i) ? f[j] : g[j - (1 << i)]; }; + f.resize(2_u32 << i); + flt_ (u32, j, 0, 2_u32 << i) f[j] = get(2 * j) * get(2 * j + 1) * ((2 * j + 1) << i); + } + if (B > SZ) { + vec g = ctsh_fps(f, val_t(SZ), base_t::ifact, u32(B - SZ)); + std::ranges::move(g, std::back_inserter(f)); + } else f.resize(B); + flt_ (u32, i, 0, (u32)B) f[i] *= val_t(i + 1) * SZ; + f.insert(f.begin(), 1); + flt_ (u32, i, 0, (u32)B) f[i + 1] *= f[i]; + } + val_t res; + u64 q = n / SZ; + u32 r = n % SZ; + if (2 * r <= SZ) { + res = f[q]; + flt_ (u32, i, 0, r) res *= n - i; + } else if (q != B) { + res = f[q + 1]; + val_t den = 1; + flt_ (u32, i, 1, SZ - r + 1) den *= n + i; + res /= den; + } else { + res = base_t::mod() - 1; + val_t den = 1; + for (u64 i = base_t::mod() - 1; i > n; --i) den *= i; + res /= den; + } + return res; + } +}; + +} // namespace tifa_libs::math + +#endif \ No newline at end of file diff --git a/src/code/comb/gen_ball_box_ii.hpp b/src/code/comb/gen_ball_box_ii.hpp index 60733c1ea..68ee279f8 100644 --- a/src/code/comb/gen_ball_box_ii.hpp +++ b/src/code/comb/gen_ball_box_ii.hpp @@ -7,8 +7,8 @@ namespace tifa_libs::math { // f = \\prod_{i=1}^m 1/(1-x^i), deg(f) = n -template -CEXP poly gen_ball_box_ii(u32 m, u32 n, spnuu inv) NE { +template +CEXP poly gen_ball_box_ii(u32 m, u32 n, vec CR inv) NE { poly f(n + 1); flt_ (u32, i, 1, m + 1) flt_ (u32, k, 1, n / i + 1) f[i * k] += inv[k]; diff --git a/src/code/comb/gen_bernoulli.hpp b/src/code/comb/gen_bernoulli.hpp index eeb01ae44..be9870899 100644 --- a/src/code/comb/gen_bernoulli.hpp +++ b/src/code/comb/gen_bernoulli.hpp @@ -8,8 +8,8 @@ namespace tifa_libs::math { // bernoulli[i] = B_i, i=0,1,...,n -template mint> -CEXP poly gen_bernoulli(u32 n, vec CR fact, vec CR ifact) NE { +template mint> +CEXP poly gen_bernoulli(u32 n, vec CR fact, vec CR ifact) NE { if (!n) return poly{1}; poly b(n + 1); flt_ (u32, i, 0, n + 1) b[i] = ifact[i + 1]; diff --git a/src/code/comb/gen_stirling1_col.hpp b/src/code/comb/gen_stirling1_col.hpp index aa7519276..e0f6dc65a 100644 --- a/src/code/comb/gen_stirling1_col.hpp +++ b/src/code/comb/gen_stirling1_col.hpp @@ -9,8 +9,8 @@ namespace tifa_libs::math { // stirling1[i] = {i \\brack k}, i=0,1,...,n -template -CEXP poly gen_stirling1_col(u32 n, u32 k, spnuu fact, spnuu inv, spnuu invfact) NE { +template +CEXP poly gen_stirling1_col(u32 n, u32 k, vec CR fact, vec CR inv, vec CR invfact) NE { if (n < k) return poly(n + 1); poly f(n + 1); flt_ (u32, i, 1, n + 1) f[i] = inv[i]; diff --git a/src/code/comb/gen_stirling1_row.hpp b/src/code/comb/gen_stirling1_row.hpp index d8247174e..2ea7ee61c 100644 --- a/src/code/comb/gen_stirling1_row.hpp +++ b/src/code/comb/gen_stirling1_row.hpp @@ -9,8 +9,8 @@ namespace tifa_libs::math { // stirling1[i] = {n \\brace i}, i=0,1,...,n -template -CEXP poly gen_stirling1_row(u32 n, spnuu fact, spnuu ifact) NE { +template +CEXP poly gen_stirling1_row(u32 n, vec CR fact, vec CR ifact) NE { using mint = TPN poly::val_t; if (!n) return poly{1}; poly f{0, 1}; diff --git a/src/code/comb/gen_stirling2_col.hpp b/src/code/comb/gen_stirling2_col.hpp index 4c8b9ccdf..f7ad002de 100644 --- a/src/code/comb/gen_stirling2_col.hpp +++ b/src/code/comb/gen_stirling2_col.hpp @@ -8,8 +8,8 @@ namespace tifa_libs::math { // stirling2[i] = {i \\brack k}, i=0,1,...,n -template -CEXP poly gen_stirling2_col(u32 n, u32 k, spnuu fact, spnuu ifact) NE { +template +CEXP poly gen_stirling2_col(u32 n, u32 k, vec CR fact, vec CR ifact) NE { using mint = TPN poly::val_t; if (k > n) return poly(n + 1); auto g = [&](auto&& g, poly& f, u32 n) NE -> void { diff --git a/src/code/comb/qbinom.hpp b/src/code/comb/qbinom.hpp index 8b0337887..ce5112adf 100644 --- a/src/code/comb/qbinom.hpp +++ b/src/code/comb/qbinom.hpp @@ -6,12 +6,11 @@ namespace tifa_libs::math { -template class fact = fact_helper> +template > struct qbinom : binom { vec qfact, iqfact; - static CEXP u64 mod() NE { return mint::mod(); } - CEXPE qbinom(u32 q, u32 max_m = fact::DEFUALT_MAX) NE : binom(max_m), qfact(2) { + CEXPE qbinom(u32 q, u32 max_m = fact::DEFUALT_MAX) NE : binom(max_m), qfact(2) { assert(q), qfact[0] = qfact[1] = 1; mint x = 1; flt_ (u32, i, 2, max_m + 1) @@ -19,7 +18,6 @@ struct qbinom : binom { else break; flt_ (u32, i, 3, (u32)qfact.size()) qfact[i] *= qfact[i - 1]; iqfact = gen_invseq(qfact); - this->reset(max_m / ((u32)qfact.size() - 1) + 2); } // $\binom{m}{n}_q$ diff --git a/src/code/poly/comp_fpssps.hpp b/src/code/poly/comp_fpssps.hpp index c0f2b5273..0bf3161d5 100644 --- a/src/code/poly/comp_fpssps.hpp +++ b/src/code/poly/comp_fpssps.hpp @@ -9,8 +9,8 @@ namespace tifa_libs::math { // @return $f(g(x_0, \dots, x_{n-1}))$ -template -auto comp_fpssps(u32 n, poly f, vec g, spnuu fact, spnuu ifact) NE { +template +auto comp_fpssps(u32 n, poly f, vec g, vec CR fact, vec CR ifact) NE { using mint = TPN poly::val_t; static conv_subset ss; if (assert(n <= N); !f.size()) return vec(1 << n); diff --git a/src/code/poly/compinv_fps.hpp b/src/code/poly/compinv_fps.hpp index f0ecdce41..23558dcb3 100644 --- a/src/code/poly/compinv_fps.hpp +++ b/src/code/poly/compinv_fps.hpp @@ -10,8 +10,8 @@ namespace tifa_libs::math { // @return g s.t. $g(f(x)) \equiv x \pmod{\deg(f)+1}$ -template