Binomial coefficient
(Mylib/Combinatorics/binomial_coefficient.cpp)
Operations
-
ext_lucas(p, q)
-
operator(n, k)
- $\binom{n}{k} \bmod p ^ q$
-
binomial_coefficient(m)
Requirements
Notes
Problems
References
Depends on
Verified with
Code
#pragma once
#include "Mylib/Number/Mod/mod_inv.cpp"
#include "Mylib/Number/Mod/mod_pow.cpp"
#include "Mylib/Number/chinese_remainder_algorithm.cpp"
#include "Mylib/Number/extended_gcd.cpp"
namespace haar_lib {
class ext_lucas {
std::vector<int64_t> prod_, inv_;
int p_, q_;
int64_t m_;
public:
ext_lucas(int p, int q) : p_(p), q_(q), m_(1) {
for (int i = 0; i < q; ++i) m_ *= p_;
prod_.assign(m_, 1);
inv_.assign(m_, 1);
for (int i = 1; i < m_; ++i) {
prod_[i] = prod_[i - 1] * (i % p_ == 0 ? 1 : i) % m_;
}
inv_[m_ - 1] = mod_inv(prod_[m_ - 1], m_);
for (int i = m_ - 1; i > 0; --i) {
inv_[i - 1] = inv_[i] * (i % p_ == 0 ? 1 : i) % m_;
}
}
int64_t operator()(int64_t n, int64_t k) const {
int64_t r = n - k;
int64_t e = 0;
int64_t eq = 0;
int64_t ret = 1;
for (int i = 0;;) {
if (n == 0) { break; }
(ret *= prod_[n % m_]) %= m_;
(ret *= inv_[k % m_]) %= m_;
(ret *= inv_[r % m_]) %= m_;
n /= p_;
k /= p_;
r /= p_;
e += n - k - r;
if (e >= q_) return 0;
i += 1;
if (i >= q_) eq += n - k - r;
}
if ((p_ != 2 or q_ < 3) and eq % 2 == 1) ret = m_ - ret;
(ret *= mod_pow(p_, e, m_)) %= m_;
return ret;
}
};
class binomial_coefficient {
std::vector<std::pair<int, int>> m_primes;
std::vector<ext_lucas> lu;
std::vector<int64_t> ms;
public:
binomial_coefficient(int m) {
for (int64_t i = 2LL; i * i <= m; ++i) {
if (m % i == 0) {
int t = 1, c = 0;
while (m % i == 0) {
m /= i;
++c;
t *= i;
}
m_primes.emplace_back(i, c);
ms.push_back(t);
}
}
if (m != 1) {
m_primes.emplace_back(m, 1);
ms.push_back(m);
}
for (auto [p, q] : m_primes) {
lu.push_back(ext_lucas(p, q));
}
}
int64_t operator()(int64_t n, int64_t k) const {
if (n < k or n < 0 or k < 0) return 0;
std::vector<int64_t> bs;
for (auto &a : lu) { bs.push_back(a(n, k)); }
return chinese_remainder_algorithm(bs, ms).value().first;
}
};
} // namespace haar_lib
#line 2 "Mylib/Number/Mod/mod_inv.cpp"
#include <cstdint>
#include <utility>
namespace haar_lib {
constexpr int64_t mod_inv(int64_t a, int64_t m) {
int64_t b = m, u = 1, v = 0;
while (b) {
int64_t t = a / b;
a -= t * b;
a = a ^ b;
b = a ^ b;
a = a ^ b;
u -= t * v;
u = u ^ v;
v = u ^ v;
u = u ^ v;
}
u %= m;
if (u < 0) u += m;
return u;
}
} // namespace haar_lib
#line 3 "Mylib/Number/Mod/mod_pow.cpp"
namespace haar_lib {
constexpr int64_t mod_pow(int64_t n, int64_t p, int64_t m) {
int64_t ret = 1;
while (p > 0) {
if (p & 1) (ret *= n) %= m;
(n *= n) %= m;
p >>= 1;
}
return ret;
}
} // namespace haar_lib
#line 2 "Mylib/Number/chinese_remainder_algorithm.cpp"
#include <cassert>
#include <optional>
#include <vector>
#line 2 "Mylib/Number/extended_gcd.cpp"
#include <tuple>
namespace haar_lib {
auto ext_gcd(int64_t a, int64_t b) -> std::tuple<
int64_t, // gcd
int64_t, // p
int64_t // q
> {
if (b == 0) return std::make_tuple(a, 1, 0);
const auto [d, q, p] = ext_gcd(b, (a + b) % b);
return std::make_tuple(d, p, q - a / b * p);
}
} // namespace haar_lib
#line 6 "Mylib/Number/chinese_remainder_algorithm.cpp"
namespace haar_lib {
std::optional<std::pair<int64_t, int64_t>> chinese_remainder_algorithm(
int64_t b1, int64_t m1,
int64_t b2, int64_t m2) {
const auto [d, p, q] = ext_gcd(m1, m2);
if ((b2 - b1) % d != 0) return std::nullopt;
const int64_t m = m1 * m2 / d;
const int64_t t = ((b2 - b1) * p / d) % (m2 / d);
const int64_t r = (b1 + m1 * t + m) % m;
return {{r, m}};
}
std::optional<std::pair<int64_t, int64_t>> chinese_remainder_algorithm(
const std::vector<int64_t> &bs,
const std::vector<int64_t> &ms) {
assert(bs.size() == ms.size());
int64_t R = 0, M = 1;
for (int i = 0; i < (int) bs.size(); ++i) {
const auto res = chinese_remainder_algorithm(R, M, bs[i], ms[i]);
if (not res) return std::nullopt;
const auto [r, m] = *res;
R = r;
M = m;
}
return {{R, M}};
}
} // namespace haar_lib
#line 6 "Mylib/Combinatorics/binomial_coefficient.cpp"
namespace haar_lib {
class ext_lucas {
std::vector<int64_t> prod_, inv_;
int p_, q_;
int64_t m_;
public:
ext_lucas(int p, int q) : p_(p), q_(q), m_(1) {
for (int i = 0; i < q; ++i) m_ *= p_;
prod_.assign(m_, 1);
inv_.assign(m_, 1);
for (int i = 1; i < m_; ++i) {
prod_[i] = prod_[i - 1] * (i % p_ == 0 ? 1 : i) % m_;
}
inv_[m_ - 1] = mod_inv(prod_[m_ - 1], m_);
for (int i = m_ - 1; i > 0; --i) {
inv_[i - 1] = inv_[i] * (i % p_ == 0 ? 1 : i) % m_;
}
}
int64_t operator()(int64_t n, int64_t k) const {
int64_t r = n - k;
int64_t e = 0;
int64_t eq = 0;
int64_t ret = 1;
for (int i = 0;;) {
if (n == 0) { break; }
(ret *= prod_[n % m_]) %= m_;
(ret *= inv_[k % m_]) %= m_;
(ret *= inv_[r % m_]) %= m_;
n /= p_;
k /= p_;
r /= p_;
e += n - k - r;
if (e >= q_) return 0;
i += 1;
if (i >= q_) eq += n - k - r;
}
if ((p_ != 2 or q_ < 3) and eq % 2 == 1) ret = m_ - ret;
(ret *= mod_pow(p_, e, m_)) %= m_;
return ret;
}
};
class binomial_coefficient {
std::vector<std::pair<int, int>> m_primes;
std::vector<ext_lucas> lu;
std::vector<int64_t> ms;
public:
binomial_coefficient(int m) {
for (int64_t i = 2LL; i * i <= m; ++i) {
if (m % i == 0) {
int t = 1, c = 0;
while (m % i == 0) {
m /= i;
++c;
t *= i;
}
m_primes.emplace_back(i, c);
ms.push_back(t);
}
}
if (m != 1) {
m_primes.emplace_back(m, 1);
ms.push_back(m);
}
for (auto [p, q] : m_primes) {
lu.push_back(ext_lucas(p, q));
}
}
int64_t operator()(int64_t n, int64_t k) const {
if (n < k or n < 0 or k < 0) return 0;
std::vector<int64_t> bs;
for (auto &a : lu) { bs.push_back(a(n, k)); }
return chinese_remainder_algorithm(bs, ms).value().first;
}
};
} // namespace haar_lib
Back to top page
Mod pow
(Mylib/Number/Mod/mod_pow.cpp)