Skip to content

Commit

Permalink
Merge pull request #91 from rainbou-kpr/89-combination
Browse files Browse the repository at this point in the history
crt, combinatorics
  • Loading branch information
KowerKoint authored Sep 27, 2023
2 parents 1a9d6b1 + 55bb8b0 commit 3532649
Show file tree
Hide file tree
Showing 7 changed files with 425 additions and 3 deletions.
10 changes: 7 additions & 3 deletions .verify-helper/timestamps.remote.json
Original file line number Diff line number Diff line change
Expand Up @@ -24,9 +24,10 @@
"test/atcoder-edpc-g.test.cpp": "2023-06-06 14:52:29 +0900",
"test/atcoder-past202012-n.test.cpp": "2023-08-01 18:34:30 +0900",
"test/atcoder-past202012-n.test.py": "2023-09-07 14:26:13 +0900",
"test/yosupo-convolution-mod-1000000007.test.cpp": "2023-09-03 11:57:52 +0900",
"test/yosupo-convolution-mod-2-64.test.cpp": "2023-09-03 11:57:52 +0900",
"test/yosupo-convolution-mod.test.cpp": "2023-09-03 11:57:52 +0900",
"test/yosupo-binomial-coefficient.test.cpp": "2023-09-16 00:25:06 +0900",
"test/yosupo-convolution-mod-1000000007.test.cpp": "2023-09-16 00:07:15 +0900",
"test/yosupo-convolution-mod-2-64.test.cpp": "2023-09-16 00:07:15 +0900",
"test/yosupo-convolution-mod.test.cpp": "2023-09-16 00:07:15 +0900",
"test/yosupo-enumerate-palindromes.test.cpp": "2023-05-31 16:27:00 +0900",
"test/yosupo-enumerate-quotients.test.cpp": "2023-05-03 12:22:21 +0900",
"test/yosupo-enumerate-quotients.test.py": "2023-05-23 13:25:17 +0900",
Expand All @@ -47,6 +48,9 @@
"test/yosupo-unionfind.test.cpp": "2023-04-28 12:54:27 +0900",
"test/yosupo-unionfind.test.py": "2023-05-26 10:38:14 +0900",
"test/yosupo-zalgorithm.test.cpp": "2023-06-15 15:22:46 +0900",
"test/yukicoder-117.test.cpp": "2023-09-16 00:11:43 +0900",
"test/yukicoder-186.test.cpp": "2023-09-16 00:07:15 +0900",
"test/yukicoder-187.test.cpp": "2023-09-16 00:07:15 +0900",
"test/yukicoder-2292.test.cpp": "2023-08-01 17:59:54 +0900",
"test/yukicoder-rotate-enlarge_1.test.cpp": "2023-09-27 09:30:35 +0900",
"test/yukicoder-rotate-enlarge_2.test.cpp": "2023-09-27 09:30:35 +0900",
Expand Down
228 changes: 228 additions & 0 deletions cpp/combinatorics.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,228 @@
#pragma once

#include <vector>
#include "number-theory.hpp"

/**
* @brief 組み合わせ
*/
template <typename Modint>
class Combination {
static std::vector<Modint> fact, inv_fact;

public:
/**
* @brief n までの階乗とその逆元を前計算する
* @param n
* @note O(n) 必要になったら呼び出されるが、予め大きなnに対して呼び出しておくことで逆元の直接計算を減らせる
*/
inline static void extend(int n) {
int m = fact.size();
if (n < m) return;
fact.resize(n + 1);
inv_fact.resize(n + 1);
for (int i = m; i <= n; ++i) {
fact[i] = fact[i - 1] * i;
}
inv_fact[n] = fact[n].inv();
for (int i = n; i > m; --i) {
inv_fact[i - 1] = inv_fact[i] * i;
}
}

/**
* @brief n の階乗を返す
* @param n
* @return n!
* @note extend(n), O(1)
*/
inline static Modint factorial(int n) {
extend(n);
return fact[n];
}
/**
* @brief n の階乗の逆元を返す
* @param n
* @return n!^-1
* @note extend(n), O(1)
*/
inline static Modint inverse_factorial(int n) {
extend(n);
return inv_fact[n];
}
/**
* @brief n の逆元を返す
* @param n
* @return n^-1
* @note extend(n), O(1)
*/
inline static Modint inverse(int n) {
extend(n);
return inv_fact[n] * fact[n - 1];
}

/**
* @brief nPr を返す
* @param n
* @param r
* @return nPr
* @note extend(n), O(1)
*/
inline static Modint P(int n, int r) {
if (r < 0 || n < r) return 0;
extend(n);
return fact[n] * inv_fact[n - r];
}
/**
* @brief nCr を返す
* @param n
* @param r
* @return nCr
* @note extend(n), O(1)
*/
inline static Modint C(int n, int r) {
if (r < 0 || n < r) return 0;
extend(n);
return fact[n] * inv_fact[r] * inv_fact[n - r];
}
/**
* @brief nHr を返す
* @param n
* @param r
* @return nHr
* @note extend(n+r-1), O(1)
*/
inline static Modint H(int n, int r) {
if (n < 0 || r < 0) return 0;
if (n == 0 && r == 0) return 1;
return C(n + r - 1, r);
}

/**
* @brief nPr を定義どおり計算する
* @param n
* @param r
* @return nPr
* @note O(r)
*/
inline static Modint P_loop(long long n, int r) {
if (r < 0 || n < r) return 0;
Modint res = 1;
for (int i = 0; i < r; ++i) {
res *= n - i;
}
return res;
}
/**
* @brief nCr を定義どおり計算する
* @param n
* @param r
* @return nCr
* @note O(min(r, n-r))
*/
inline static Modint C_loop(long long n, long long r) {
if (r < 0 || n < r) return 0;
if(r > n - r) r = n - r;
extend(r);
return P_loop(n, r) * inv_fact[r];
}
/**
* @brief nHr を定義どおり計算する
* @param n
* @param r
* @return nHr
* @note O(r)
*/
inline static Modint H_loop(long long n, long long r) {
if (n < 0 || r < 0) return 0;
if (n == 0 && r == 0) return 1;
return C_loop(n + r - 1, r);
}

/**
* @brief nCr を Lucas の定理を用いて計算する
* @param n
* @param r
* @return nCr
* @note expand(Mod), O(log(r))
*/
inline static Modint C_lucas(long long n, long long r) {
if (r < 0 || n < r) return 0;
if (r == 0 || n == r) return 1;
Modint res = 1;
while(r > 0) {
int ni = n % Modint::mod(), ri = r % Modint::mod();
if (ni < ri) return 0;
res *= C(ni, ri);
n /= Modint::mod();
r /= Modint::mod();
}
return res;
}
};

template <typename Modint>
std::vector<Modint> Combination<Modint>::fact{1, 1};
template <typename Modint>
std::vector<Modint> Combination<Modint>::inv_fact{1, 1};

/**
* @brief mod p^q での二項係数を求める構造
* @note 前計算O(p^q) 参考: https://nyaannyaan.github.io/library/modulo/arbitrary-mod-binomial.hpp
*/
struct CombinationPQ {
int p, q;
int pq;
std::vector<int> fact_p, inv_fact_p;
int delta;
CombinationPQ(int p, int q) : p(p), q(q) {
pq = 1;
for(int i = 0; i < q; i++) pq *= p;
fact_p.resize(pq);
fact_p[0] = 1;
for(int i = 1; i < pq; i++) {
if(i % p == 0) fact_p[i] = fact_p[i - 1];
else fact_p[i] = (long long)fact_p[i - 1] * i % pq;
}
inv_fact_p.resize(pq);
inv_fact_p[pq - 1] = modinv(fact_p[pq - 1], pq);
for(int i = pq - 1; i > 0; i--) {
if(i % p == 0) inv_fact_p[i - 1] = inv_fact_p[i];
else inv_fact_p[i - 1] = (long long)inv_fact_p[i] * i % pq;
}
if(p == 2 && q >= 3) delta = 1;
else delta = -1;
}
/**
* @brief nCr mod p^q を返す
* @param n long long
* @param r long long
* @return nCr mod p^q
* @note O(log(n))
*/
int C(long long n, long long r) {
if(r < 0 || n < r) return 0;
long long m = n - r;
int ans = 1;
std::vector<int> epsilon;
while(n > 0) {
ans = (long long)ans * fact_p[n % pq] % pq;
ans = (long long)ans * inv_fact_p[m % pq] % pq;
ans = (long long)ans * inv_fact_p[r % pq] % pq;
n /= p;
m /= p;
r /= p;
epsilon.push_back(n - m - r);
}
if(delta == -1 && epsilon.size() >= q && accumulate(epsilon.begin()+q-1, epsilon.end(), 0) % 2 == 1) ans = pq - ans;
if(ans == pq) ans = 0;
int e = accumulate(epsilon.begin(), epsilon.end(), 0);
if(e >= q) ans = 0;
else {
for(int i = 0; i < e; i++) {
ans = (long long)ans * p % pq;
}
}
return ans;
}
};
84 changes: 84 additions & 0 deletions cpp/number-theory.hpp
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#pragma once

#include <numeric>
#include <vector>
#include "modint.hpp"

Expand Down Expand Up @@ -43,6 +44,89 @@ long long modpow(long long a, long long n, long long MOD) {
return res;
}

/**
* @brief 2式の連立合同式を、mが互いに素になるように変形する
* @param r1 long long
* @param m1 long long
* @param r2 long long
* @param m2 long long
* @note 矛盾する場合、r1 = r2 = m1 = m2 = -1となる
*/
void coprimize_simulaneous_congruence_equation(long long& r1, long long& m1, long long& r2, long long& m2) {
long long g = std::gcd(m1, m2);
if((r2 - r1) % g != 0) {
r1 = r2 = m1 = m2 = -1;
return;
}
m1 /= g, m2 /= g;
long long gi = std::gcd(g, m1);
long long gj = g / gi;
do {
g = std::gcd(gi, gj);
gi *= g, gj /= g;
} while(g != 1);
m1 *= gi, m2 *= gj;
r1 %= m1, r2 %= m2;
}

/**
* @brief 連立合同式を解く
* @param r vector<long long> 余りの配列
* @param m vector<long long> modの配列
* @return std::pair<long long, long long> (解, LCM) 解なしのときは{-1, -1}
*/
std::pair<long long, long long> crt(const std::vector<long long>& r, const std::vector<long long>& m) {
assert(r.size() == m.size());
if(r.size() == 0) return {0, 1};
int n = (int)r.size();
long long m_lcm = m[0];
long long ans = r[0] % m[0];
for (int i = 1; i < n; i++) {
long long rr = r[i] % m[i], mm = m[i];
coprimize_simulaneous_congruence_equation(ans, m_lcm, rr, mm);
if(m_lcm == -1) return {-1, -1};
long long t = ((rr - ans) * modinv(m_lcm, mm)) % mm;
if(t < 0) t += mm;
ans += t * m_lcm;
m_lcm *= mm;
}
return {ans, m_lcm};
}

/**
* @brief 連立合同式の最小の非負整数解 % MODを求める
* @param r vector<long long> 余りの配列
* @param m vector<long long> modの配列
* @param MOD long long
* @return std::pair<long long, long long> (最小解 % MOD, LCM % MOD) 解なしのときは{-1, -1}
*/
std::pair<long long, long long> crt(const std::vector<long long>& r, const std::vector<long long>& m, long long MOD) {
assert(r.size() == m.size());
if(r.size() == 0) return {0, 1};
int n = (int)r.size();
std::vector<long long> r2 = r, m2 = m;
// mを互いに素にする
for(int i = 1; i < n; i++) {
for(int j = 0; j < i; j++) {
coprimize_simulaneous_congruence_equation(r2[i], m2[i], r2[j], m2[j]);
if(m2[i] == -1) return {-1, -1};
}
}

m2.push_back(MOD);
std::vector<long long> prod(n+1, 1); // m2[0] * ... * m2[i - 1] mod m2[i]
std::vector<long long> x(n+1, 0); // i番目までの解 mod m2[i]
for(int i = 0; i < n; i++) {
long long t = (r2[i] - x[i]) * modinv(prod[i], m2[i]) % m2[i];
if(t < 0) t += m2[i];
for(int j = i + 1; j <= n; j++) {
(x[j] += t * prod[j]) %= m2[j];
(prod[j] *= m2[i]) %= m2[j];
}
}
return {x[n], prod[n]};
}

/**
* @brief 畳み込み
*/
Expand Down
40 changes: 40 additions & 0 deletions test/yosupo-binomial-coefficient.test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
#define PROBLEM "https://judge.yosupo.jp/problem/binomial_coefficient"

#include "../cpp/combinatorics.hpp"
#include "../cpp/number-theory.hpp"

#include <iostream>

int main() {
int t, m; std::cin >> t >> m;
std::vector<std::pair<int, int>> factors;
{
int mm = m;
for(int i = 2; i * i <= mm; ++i) {
if(mm % i == 0) {
int cnt = 0;
while(mm % i == 0) {
mm /= i;
cnt++;
}
factors.emplace_back(i, cnt);
}
}
if(mm > 1) factors.emplace_back(mm, 1);
}
std::vector<CombinationPQ> combpq;
for(int i = 0; i < factors.size(); ++i) {
auto [p, e] = factors[i];
combpq.emplace_back(p, e);
}
while(t--) {
long long n, k; std::cin >> n >> k;
std::vector<long long> x, y;
for(int i = 0; i < factors.size(); ++i) {
auto [p, e] = factors[i];
x.push_back(combpq[i].C(n, k));
y.push_back(combpq[i].pq);
}
std::cout << crt(x, y).first << std::endl;
}
}
Loading

0 comments on commit 3532649

Please sign in to comment.