src/jiangly/math/13C-Q-Binomial.hpp
Code
/** 生成函数（Binomial 任意模数二项式）
* 2023-08-22: https://codeforces.com/contest/896/submission/219861532
**/
std::vector<std::pair<int, int>> factorize(int n) {
std::vector<std::pair<int, int>> factors;
for (int i = 2; static_cast<long long>(i) * i <= n; i++) {
if (n % i == 0) {
int t = 0;
for (; n % i == 0; n /= i)
++t;
factors.emplace_back(i, t);
}
}
if (n > 1)
factors.emplace_back(n, 1);
return factors;
}
constexpr int power(int base, i64 exp) {
int res = 1;
for (; exp > 0; base *= base, exp /= 2) {
if (exp % 2 == 1) {
res *= base;
}
}
return res;
}
constexpr int power(int base, i64 exp, int mod) {
int res = 1 % mod;
for (; exp > 0; base = 1LL * base * base % mod, exp /= 2) {
if (exp % 2 == 1) {
res = 1LL * res * base % mod;
}
}
return res;
}
int inverse(int a, int m) {
int g = m, r = a, x = 0, y = 1;
while (r != 0) {
int q = g / r;
g %= r;
std::swap(g, r);
x -= q * y;
std::swap(x, y);
}
return x < 0 ? x + m : x;
}
int solveModuloEquations(const std::vector<std::pair<int, int>> &e) {
int m = 1;
for (std::size_t i = 0; i < e.size(); i++) {
m *= e[i].first;
}
int res = 0;
for (std::size_t i = 0; i < e.size(); i++) {
int p = e[i].first;
res = (res + 1LL * e[i].second * (m / p) * inverse(m / p, p)) % m;
}
return res;
}
constexpr int N = 1E5;
class Binomial {
const int mod;
private:
const std::vector<std::pair<int, int>> factors;
std::vector<int> pk;
std::vector<std::vector<int>> prod;
static constexpr i64 exponent(i64 n, int p) {
i64 res = 0;
for (n /= p; n > 0; n /= p) {
res += n;
}
return res;
}
int product(i64 n, std::size_t i) {
int res = 1;
int p = factors[i].first;
for (; n > 0; n /= p) {
res = 1LL * res * power(prod[i].back(), n / pk[i], pk[i]) % pk[i] * prod[i][n % pk[i]] % pk[i];
}
return res;
}
public:
Binomial(int mod) : mod(mod), factors(factorize(mod)) {
pk.resize(factors.size());
prod.resize(factors.size());
for (std::size_t i = 0; i < factors.size(); i++) {
int p = factors[i].first;
int k = factors[i].second;
pk[i] = power(p, k);
prod[i].resize(std::min(N + 1, pk[i]));
prod[i][0] = 1;
for (int j = 1; j < prod[i].size(); j++) {
if (j % p == 0) {
prod[i][j] = prod[i][j - 1];
} else {
prod[i][j] = 1LL * prod[i][j - 1] * j % pk[i];
}
}
}
}
int operator()(i64 n, i64 m) {
if (n < m || m < 0) {
return 0;
}
std::vector<std::pair<int, int>> ans(factors.size());
for (int i = 0; i < factors.size(); i++) {
int p = factors[i].first;
int k = factors[i].second;
int e = exponent(n, p) - exponent(m, p) - exponent(n - m, p);
if (e >= k) {
ans[i] = std::make_pair(pk[i], 0);
} else {
int pn = product(n, i);
int pm = product(m, i);
int pd = product(n - m, i);
int res = 1LL * pn * inverse(pm, pk[i]) % pk[i] * inverse(pd, pk[i]) % pk[i] * power(p, e) % pk[i];
ans[i] = std::make_pair(pk[i], res);
}
}
return solveModuloEquations(ans);
}
};
#line 1 "src/jiangly/math/13C-Q-Binomial.hpp"
/** 生成函数（Binomial 任意模数二项式）
* 2023-08-22: https://codeforces.com/contest/896/submission/219861532
**/
std::vector<std::pair<int, int>> factorize(int n) {
std::vector<std::pair<int, int>> factors;
for (int i = 2; static_cast<long long>(i) * i <= n; i++) {
if (n % i == 0) {
int t = 0;
for (; n % i == 0; n /= i)
++t;
factors.emplace_back(i, t);
}
}
if (n > 1)
factors.emplace_back(n, 1);
return factors;
}
constexpr int power(int base, i64 exp) {
int res = 1;
for (; exp > 0; base *= base, exp /= 2) {
if (exp % 2 == 1) {
res *= base;
}
}
return res;
}
constexpr int power(int base, i64 exp, int mod) {
int res = 1 % mod;
for (; exp > 0; base = 1LL * base * base % mod, exp /= 2) {
if (exp % 2 == 1) {
res = 1LL * res * base % mod;
}
}
return res;
}
int inverse(int a, int m) {
int g = m, r = a, x = 0, y = 1;
while (r != 0) {
int q = g / r;
g %= r;
std::swap(g, r);
x -= q * y;
std::swap(x, y);
}
return x < 0 ? x + m : x;
}
int solveModuloEquations(const std::vector<std::pair<int, int>> &e) {
int m = 1;
for (std::size_t i = 0; i < e.size(); i++) {
m *= e[i].first;
}
int res = 0;
for (std::size_t i = 0; i < e.size(); i++) {
int p = e[i].first;
res = (res + 1LL * e[i].second * (m / p) * inverse(m / p, p)) % m;
}
return res;
}
constexpr int N = 1E5;
class Binomial {
const int mod;
private:
const std::vector<std::pair<int, int>> factors;
std::vector<int> pk;
std::vector<std::vector<int>> prod;
static constexpr i64 exponent(i64 n, int p) {
i64 res = 0;
for (n /= p; n > 0; n /= p) {
res += n;
}
return res;
}
int product(i64 n, std::size_t i) {
int res = 1;
int p = factors[i].first;
for (; n > 0; n /= p) {
res = 1LL * res * power(prod[i].back(), n / pk[i], pk[i]) % pk[i] * prod[i][n % pk[i]] % pk[i];
}
return res;
}
public:
Binomial(int mod) : mod(mod), factors(factorize(mod)) {
pk.resize(factors.size());
prod.resize(factors.size());
for (std::size_t i = 0; i < factors.size(); i++) {
int p = factors[i].first;
int k = factors[i].second;
pk[i] = power(p, k);
prod[i].resize(std::min(N + 1, pk[i]));
prod[i][0] = 1;
for (int j = 1; j < prod[i].size(); j++) {
if (j % p == 0) {
prod[i][j] = prod[i][j - 1];
} else {
prod[i][j] = 1LL * prod[i][j - 1] * j % pk[i];
}
}
}
}
int operator()(i64 n, i64 m) {
if (n < m || m < 0) {
return 0;
}
std::vector<std::pair<int, int>> ans(factors.size());
for (int i = 0; i < factors.size(); i++) {
int p = factors[i].first;
int k = factors[i].second;
int e = exponent(n, p) - exponent(m, p) - exponent(n - m, p);
if (e >= k) {
ans[i] = std::make_pair(pk[i], 0);
} else {
int pn = product(n, i);
int pm = product(m, i);
int pd = product(n - m, i);
int res = 1LL * pn * inverse(pm, pk[i]) % pk[i] * inverse(pd, pk[i]) % pk[i] * power(p, e) % pk[i];
ans[i] = std::make_pair(pk[i], res);
}
}
return solveModuloEquations(ans);
}
};
Back to top page