binomial_mod.hpp¶
预处理固定合数模数后回答多次组合数 \(\binom nk\)；用素数幂分解和 CRT 处理不可逆阶乘。
\[
\displaystyle \binom{n}{k} \bmod m
\]
Complexity: Time: O(m) preprocessing and O(sum_{p|m} log_p n) per table query; O(min(k,n-k) omega(m)) for the one-shot function. Space: O(m) for the table and O(omega(m)) for the one-shot function.
AC 记录：binomial_coefficient。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(m) preprocessing and O(sum_{p|m} log_p n) per table
/// query; O(min(k,n-k) omega(m)) for the one-shot function.
/// Space: O(m) for the table and O(omega(m)) for the one-shot function.
#include "noya/factorize.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <vector>
namespace noya {
namespace binomial_mod_internal {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128;
inline u64 multiply(u64 a, u64 b, u64 mod) { return u64(u128(a) * b % mod); }
inline u64 power(u64 val, u64 exp, u64 mod) {
u64 res = 1 % mod;
while (exp > 0) {
if (exp & 1) {
res = multiply(res, val, mod);
}
val = multiply(val, val, mod);
exp >>= 1;
}
return res;
}
inline u64 inverse(u64 val, u64 mod) {
i128 olr = val;
i128 r = mod;
i128 os = 1;
i128 s = 0;
while (r != 0) {
i128 quo = olr / r;
i128 nr = olr - quo * r;
olr = r;
r = nr;
i128 ns = os - quo * s;
os = s;
s = ns;
}
assert(olr == 1);
os %= i128(mod);
if (os < 0) {
os += mod;
}
return u64(os);
}
inline u64 prime_power_residue(u64 n, u64 k, u64 p, int exp, u64 mod) {
std::int64_t vp = 0;
u64 uni = 1 % mod;
for (u64 idx = 1; idx <= k; idx++) {
u64 num = n - k + idx;
u64 den = idx;
while (num % p == 0) {
num /= p;
vp++;
}
while (den % p == 0) {
den /= p;
vp--;
}
uni = multiply(uni, num % mod, mod);
uni = multiply(uni, inverse(den % mod, mod), mod);
}
assert(vp >= 0);
if (vp >= exp) {
return 0;
}
return multiply(uni, power(p, vp, mod), mod);
}
struct prime_power_binomial {
u64 p = 0;
int exp = 0;
u64 mod = 1;
std::vector<std::uint32_t> pre;
prime_power_binomial(u64 p_, int ex_) : p(p_), exp(ex_) {
for (int i = 0; i < exp; i++) {
mod *= p;
}
pre.resize(std::size_t(mod) + 1);
pre[0] = 1;
for (u64 i = 1; i <= mod; i++) {
pre[i] = pre[i - 1];
if (i % p != 0) {
pre[i] = std::uint32_t(multiply(pre[i], i, mod));
}
}
}
u64 valuation_factorial(u64 n) const {
u64 res = 0;
while (n > 0) {
n /= p;
res += n;
}
return res;
}
u64 unit_factorial(u64 n) const {
u64 res = 1;
while (n > 0) {
if ((n / mod) & 1) {
res = multiply(res, pre[mod], mod);
}
res = multiply(res, pre[n % mod], mod);
n /= p;
}
return res;
}
u64 choose(u64 n, u64 k) const {
if (k > n) {
return 0;
}
u64 vp = valuation_factorial(n) - valuation_factorial(k) -
valuation_factorial(n - k);
if (vp >= u64(exp)) {
return 0;
}
u64 uni = unit_factorial(n);
uni = multiply(uni, inverse(unit_factorial(k), mod), mod);
uni = multiply(uni, inverse(unit_factorial(n - k), mod), mod);
return multiply(uni, power(p, vp, mod), mod);
}
};
} // namespace binomial_mod_internal
/// @brief Preprocess one fixed composite modulus for many binomial queries.
/// Each prime-power factor stores prefix products with multiples of its prime
/// removed. Splitting n! into complete residue blocks and recursively stripping
/// one prime from every multiple yields its unit part and p-adic valuation in
/// logarithmic time. The prime-power answers are joined by CRT.
class binomial_mod_table {
public:
explicit binomial_mod_table(std::uint32_t mod) : md(mod) {
assert(mod >= 1);
if (mod == 1) {
return;
}
for (auto [p, exp] : factorize(mod)) {
fs.emplace_back(p, exp);
}
}
std::uint32_t modulus() const { return md; }
std::uint32_t choose(std::uint64_t n, std::uint64_t k) const {
using namespace binomial_mod_internal;
if (md == 1 || k > n) {
return 0;
}
u64 res = 0;
u64 cm = 1;
for (const auto &fct : fs) {
u64 rsd = fct.choose(n, k);
u64 cur = res % fct.mod;
u64 dif = (rsd + fct.mod - cur) % fct.mod;
u64 scl = multiply(dif, inverse(cm % fct.mod, fct.mod), fct.mod);
res += cm * scl;
cm *= fct.mod;
res %= cm;
}
return std::uint32_t(res);
}
private:
std::uint32_t md;
std::vector<binomial_mod_internal::prime_power_binomial> fs;
};
/// @brief Compute C(n,k) modulo an arbitrary 64-bit positive modulus in
/// O(min(k,n-k) times the number of prime factors of modulus).
inline std::uint64_t binomial_mod(std::uint64_t n, std::uint64_t k,
std::uint64_t mod) {
using namespace binomial_mod_internal;
assert(mod >= 1);
if (mod == 1 || k > n) {
return 0;
}
k = std::min(k, n - k);
u64 res = 0;
u64 cm = 1;
for (auto [p, exp] : factorize(mod)) {
u64 pp = 1;
for (int cnt = 0; cnt < exp; cnt++) {
pp *= p;
}
u64 rsd = prime_power_residue(n, k, p, exp, pp);
u64 cur = res % pp;
u64 dif = u64((u128(rsd) + pp - cur) % pp);
u64 scl = multiply(dif, inverse(cm % pp, pp), pp);
res = u64(u128(res) + u128(cm) * scl);
cm *= pp;
res %= cm;
}
return res;
}
} // namespace noya
#ifndef NOYA_BINOMIAL_MOD_HPP
#define NOYA_BINOMIAL_MOD_HPP 1
/// @complexity Time: O(m) preprocessing and O(sum_{p|m} log_p n) per table
/// query; O(min(k,n-k) omega(m)) for the one-shot function.
/// Space: O(m) for the table and O(omega(m)) for the one-shot function.
#include "noya/factorize.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <vector>
namespace noya {
namespace binomial_mod_internal {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128;
inline u64 multiply(u64 a, u64 b, u64 mod) { return u64(u128(a) * b % mod); }
inline u64 power(u64 val, u64 exp, u64 mod) {
u64 res = 1 % mod;
while (exp > 0) {
if (exp & 1) {
res = multiply(res, val, mod);
}
val = multiply(val, val, mod);
exp >>= 1;
}
return res;
}
inline u64 inverse(u64 val, u64 mod) {
i128 olr = val;
i128 r = mod;
i128 os = 1;
i128 s = 0;
while (r != 0) {
i128 quo = olr / r;
i128 nr = olr - quo * r;
olr = r;
r = nr;
i128 ns = os - quo * s;
os = s;
s = ns;
}
assert(olr == 1);
os %= i128(mod);
if (os < 0) {
os += mod;
}
return u64(os);
}
inline u64 prime_power_residue(u64 n, u64 k, u64 p, int exp, u64 mod) {
std::int64_t vp = 0;
u64 uni = 1 % mod;
for (u64 idx = 1; idx <= k; idx++) {
u64 num = n - k + idx;
u64 den = idx;
while (num % p == 0) {
num /= p;
vp++;
}
while (den % p == 0) {
den /= p;
vp--;
}
uni = multiply(uni, num % mod, mod);
uni = multiply(uni, inverse(den % mod, mod), mod);
}
assert(vp >= 0);
if (vp >= exp) {
return 0;
}
return multiply(uni, power(p, vp, mod), mod);
}
struct prime_power_binomial {
u64 p = 0;
int exp = 0;
u64 mod = 1;
std::vector<std::uint32_t> pre;
prime_power_binomial(u64 p_, int ex_) : p(p_), exp(ex_) {
for (int i = 0; i < exp; i++) {
mod *= p;
}
pre.resize(std::size_t(mod) + 1);
pre[0] = 1;
for (u64 i = 1; i <= mod; i++) {
pre[i] = pre[i - 1];
if (i % p != 0) {
pre[i] = std::uint32_t(multiply(pre[i], i, mod));
}
}
}
u64 valuation_factorial(u64 n) const {
u64 res = 0;
while (n > 0) {
n /= p;
res += n;
}
return res;
}
u64 unit_factorial(u64 n) const {
u64 res = 1;
while (n > 0) {
if ((n / mod) & 1) {
res = multiply(res, pre[mod], mod);
}
res = multiply(res, pre[n % mod], mod);
n /= p;
}
return res;
}
u64 choose(u64 n, u64 k) const {
if (k > n) {
return 0;
}
u64 vp = valuation_factorial(n) - valuation_factorial(k) -
valuation_factorial(n - k);
if (vp >= u64(exp)) {
return 0;
}
u64 uni = unit_factorial(n);
uni = multiply(uni, inverse(unit_factorial(k), mod), mod);
uni = multiply(uni, inverse(unit_factorial(n - k), mod), mod);
return multiply(uni, power(p, vp, mod), mod);
}
};
} // namespace binomial_mod_internal
/// @brief Preprocess one fixed composite modulus for many binomial queries.
/// Each prime-power factor stores prefix products with multiples of its prime
/// removed. Splitting n! into complete residue blocks and recursively stripping
/// one prime from every multiple yields its unit part and p-adic valuation in
/// logarithmic time. The prime-power answers are joined by CRT.
class binomial_mod_table {
public:
explicit binomial_mod_table(std::uint32_t mod) : md(mod) {
assert(mod >= 1);
if (mod == 1) {
return;
}
for (auto [p, exp] : factorize(mod)) {
fs.emplace_back(p, exp);
}
}
std::uint32_t modulus() const { return md; }
std::uint32_t choose(std::uint64_t n, std::uint64_t k) const {
using namespace binomial_mod_internal;
if (md == 1 || k > n) {
return 0;
}
u64 res = 0;
u64 cm = 1;
for (const auto &fct : fs) {
u64 rsd = fct.choose(n, k);
u64 cur = res % fct.mod;
u64 dif = (rsd + fct.mod - cur) % fct.mod;
u64 scl = multiply(dif, inverse(cm % fct.mod, fct.mod), fct.mod);
res += cm * scl;
cm *= fct.mod;
res %= cm;
}
return std::uint32_t(res);
}
private:
std::uint32_t md;
std::vector<binomial_mod_internal::prime_power_binomial> fs;
};
/// @brief Compute C(n,k) modulo an arbitrary 64-bit positive modulus in
/// O(min(k,n-k) times the number of prime factors of modulus).
inline std::uint64_t binomial_mod(std::uint64_t n, std::uint64_t k,
std::uint64_t mod) {
using namespace binomial_mod_internal;
assert(mod >= 1);
if (mod == 1 || k > n) {
return 0;
}
k = std::min(k, n - k);
u64 res = 0;
u64 cm = 1;
for (auto [p, exp] : factorize(mod)) {
u64 pp = 1;
for (int cnt = 0; cnt < exp; cnt++) {
pp *= p;
}
u64 rsd = prime_power_residue(n, k, p, exp, pp);
u64 cur = res % pp;
u64 dif = u64((u128(rsd) + pp - cur) % pp);
u64 scl = multiply(dif, inverse(cm % pp, pp), pp);
res = u64(u128(res) + u128(cm) * scl);
cm *= pp;
res %= cm;
}
return res;
}
} // namespace noya
#endif // NOYA_BINOMIAL_MOD_HPP
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <utility>
#include <vector>
/// @complexity Time: O(m) preprocessing and O(sum_{p|m} log_p n) per table
/// query; O(min(k,n-k) omega(m)) for the one-shot function.
/// Space: O(m) for the table and O(omega(m)) for the one-shot function.
/// @complexity Time: O(log^3 n) primality testing; Pollard-rho factorization is expected about O(n^(1/4)).
/// Space: O(log n) recursion and factors.
namespace noya {
namespace factorize_internal {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
inline u64 multiply_mod(u64 a, u64 b, u64 mod) {
return u64(u128(a) * b % mod);
}
inline u64 power_mod(u64 a, u64 exp, u64 mod) {
u64 res = 1;
while (exp > 0) {
if (exp & 1) {
res = multiply_mod(res, a, mod);
}
a = multiply_mod(a, a, mod);
exp >>= 1;
}
return res;
}
inline bool miller_rabin(u64 n) {
if (n < 2) {
return false;
}
for (u64 p :
std::array<u64, 12>{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
if (n % p == 0) {
return n == p;
}
}
int shf = __builtin_ctzll(n - 1);
u64 odd = (n - 1) >> shf;
for (u64 bas :
std::array<u64, 7>{2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
if (bas % n == 0) {
continue;
}
u64 val = power_mod(bas % n, odd, n);
if (val == 1 || val == n - 1) {
continue;
}
bool cmp = true;
for (int i = 1; i < shf; i++) {
val = multiply_mod(val, val, n);
if (val == n - 1) {
cmp = false;
break;
}
}
if (cmp) {
return false;
}
}
return true;
}
inline u64 splitmix64(u64 &st) {
u64 z = (st += 0x9e3779b97f4a7c15ULL);
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
return z ^ (z >> 31);
}
inline u64 pollard_rho(u64 n) {
if (n % 2 == 0) {
return 2;
}
if (n % 3 == 0) {
return 3;
}
static u64 st = 0x123456789abcdef0ULL;
while (true) {
u64 y = splitmix64(st) % (n - 1) + 1;
u64 c = splitmix64(st) % (n - 1) + 1;
constexpr u64 blk = 128;
u64 g = 1;
u64 r = 1;
u64 q = 1;
u64 x = 0;
u64 sy = 0;
auto nxt = [&](u64 val) {
return u64((u128(multiply_mod(val, val, n)) + c) % n);
};
while (g == 1) {
x = y;
for (u64 i = 0; i < r; i++) {
y = nxt(y);
}
for (u64 off = 0; off < r && g == 1; off += blk) {
sy = y;
for (u64 i = 0; i < std::min(blk, r - off); i++) {
y = nxt(y);
u64 dif = x > y ? x - y : y - x;
q = multiply_mod(q, dif, n);
}
g = std::gcd(q, n);
}
r <<= 1;
}
if (g == n) {
do {
sy = nxt(sy);
u64 dif = x > sy ? x - sy : sy - x;
g = std::gcd(dif, n);
} while (g == 1);
}
if (g != n) {
return g;
}
}
}
inline void collect_factors(u64 n, std::vector<u64> &res) {
if (n == 1) {
return;
}
if (miller_rabin(n)) {
res.push_back(n);
return;
}
u64 fct = pollard_rho(n);
collect_factors(fct, res);
collect_factors(n / fct, res);
}
} // namespace factorize_internal
/// @brief Deterministic Miller-Rabin primality test for unsigned 64-bit
/// integers.
inline bool is_prime(std::uint64_t n) {
return factorize_internal::miller_rabin(n);
}
/// @brief Return the prime factors of n with multiplicity in increasing order.
inline std::vector<std::uint64_t> prime_factors(std::uint64_t n) {
assert(n >= 1);
std::vector<std::uint64_t> res;
factorize_internal::collect_factors(n, res);
std::sort(res.begin(), res.end());
return res;
}
/// @brief Return the prime factorization of n as (prime, exponent) pairs.
inline std::vector<std::pair<std::uint64_t, int>> factorize(std::uint64_t n) {
std::vector<std::pair<std::uint64_t, int>> res;
for (std::uint64_t p : prime_factors(n)) {
if (res.empty() || res.back().first != p) {
res.emplace_back(p, 1);
} else {
res.back().second++;
}
}
return res;
}
} // namespace noya
namespace noya {
namespace binomial_mod_internal {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128;
inline u64 multiply(u64 a, u64 b, u64 mod) { return u64(u128(a) * b % mod); }
inline u64 power(u64 val, u64 exp, u64 mod) {
u64 res = 1 % mod;
while (exp > 0) {
if (exp & 1) {
res = multiply(res, val, mod);
}
val = multiply(val, val, mod);
exp >>= 1;
}
return res;
}
inline u64 inverse(u64 val, u64 mod) {
i128 olr = val;
i128 r = mod;
i128 os = 1;
i128 s = 0;
while (r != 0) {
i128 quo = olr / r;
i128 nr = olr - quo * r;
olr = r;
r = nr;
i128 ns = os - quo * s;
os = s;
s = ns;
}
assert(olr == 1);
os %= i128(mod);
if (os < 0) {
os += mod;
}
return u64(os);
}
inline u64 prime_power_residue(u64 n, u64 k, u64 p, int exp, u64 mod) {
std::int64_t vp = 0;
u64 uni = 1 % mod;
for (u64 idx = 1; idx <= k; idx++) {
u64 num = n - k + idx;
u64 den = idx;
while (num % p == 0) {
num /= p;
vp++;
}
while (den % p == 0) {
den /= p;
vp--;
}
uni = multiply(uni, num % mod, mod);
uni = multiply(uni, inverse(den % mod, mod), mod);
}
assert(vp >= 0);
if (vp >= exp) {
return 0;
}
return multiply(uni, power(p, vp, mod), mod);
}
struct prime_power_binomial {
u64 p = 0;
int exp = 0;
u64 mod = 1;
std::vector<std::uint32_t> pre;
prime_power_binomial(u64 p_, int ex_) : p(p_), exp(ex_) {
for (int i = 0; i < exp; i++) {
mod *= p;
}
pre.resize(std::size_t(mod) + 1);
pre[0] = 1;
for (u64 i = 1; i <= mod; i++) {
pre[i] = pre[i - 1];
if (i % p != 0) {
pre[i] = std::uint32_t(multiply(pre[i], i, mod));
}
}
}
u64 valuation_factorial(u64 n) const {
u64 res = 0;
while (n > 0) {
n /= p;
res += n;
}
return res;
}
u64 unit_factorial(u64 n) const {
u64 res = 1;
while (n > 0) {
if ((n / mod) & 1) {
res = multiply(res, pre[mod], mod);
}
res = multiply(res, pre[n % mod], mod);
n /= p;
}
return res;
}
u64 choose(u64 n, u64 k) const {
if (k > n) {
return 0;
}
u64 vp = valuation_factorial(n) - valuation_factorial(k) -
valuation_factorial(n - k);
if (vp >= u64(exp)) {
return 0;
}
u64 uni = unit_factorial(n);
uni = multiply(uni, inverse(unit_factorial(k), mod), mod);
uni = multiply(uni, inverse(unit_factorial(n - k), mod), mod);
return multiply(uni, power(p, vp, mod), mod);
}
};
} // namespace binomial_mod_internal
/// @brief Preprocess one fixed composite modulus for many binomial queries.
/// Each prime-power factor stores prefix products with multiples of its prime
/// removed. Splitting n! into complete residue blocks and recursively stripping
/// one prime from every multiple yields its unit part and p-adic valuation in
/// logarithmic time. The prime-power answers are joined by CRT.
class binomial_mod_table {
public:
explicit binomial_mod_table(std::uint32_t mod) : md(mod) {
assert(mod >= 1);
if (mod == 1) {
return;
}
for (auto [p, exp] : factorize(mod)) {
fs.emplace_back(p, exp);
}
}
std::uint32_t modulus() const { return md; }
std::uint32_t choose(std::uint64_t n, std::uint64_t k) const {
using namespace binomial_mod_internal;
if (md == 1 || k > n) {
return 0;
}
u64 res = 0;
u64 cm = 1;
for (const auto &fct : fs) {
u64 rsd = fct.choose(n, k);
u64 cur = res % fct.mod;
u64 dif = (rsd + fct.mod - cur) % fct.mod;
u64 scl = multiply(dif, inverse(cm % fct.mod, fct.mod), fct.mod);
res += cm * scl;
cm *= fct.mod;
res %= cm;
}
return std::uint32_t(res);
}
private:
std::uint32_t md;
std::vector<binomial_mod_internal::prime_power_binomial> fs;
};
/// @brief Compute C(n,k) modulo an arbitrary 64-bit positive modulus in
/// O(min(k,n-k) times the number of prime factors of modulus).
inline std::uint64_t binomial_mod(std::uint64_t n, std::uint64_t k,
std::uint64_t mod) {
using namespace binomial_mod_internal;
assert(mod >= 1);
if (mod == 1 || k > n) {
return 0;
}
k = std::min(k, n - k);
u64 res = 0;
u64 cm = 1;
for (auto [p, exp] : factorize(mod)) {
u64 pp = 1;
for (int cnt = 0; cnt < exp; cnt++) {
pp *= p;
}
u64 rsd = prime_power_residue(n, k, p, exp, pp);
u64 cur = res % pp;
u64 dif = u64((u128(rsd) + pp - cur) % pp);
u64 scl = multiply(dif, inverse(cm % pp, pp), pp);
res = u64(u128(res) + u128(cm) * scl);
cm *= pp;
res %= cm;
}
return res;
}
} // namespace noya