[0] 扱う問題
問題解説も同時にするので、問題を先に見たい人はここを開いて、適当に検索してください。
扱う問題一覧
ABC 336 G
Codeforces 925 G
AGC 051 D
ARC 146 E
Codeforces Hello 2024 E
ARC 117 E
yukicoder No. 310
[1] 概要
BEST theorem (BEST 定理) とはオイラーグラフのオイラー路の数え上げに非常に有効な定理です。
この定理を用いることで解くことができる競プロの問題に対する解法と、その実装例を紹介します。
[2] BEST 定理とは
グラフの $G$ 全ての辺を一回ずつ通る閉路をオイラー閉路と言います。$G$ がオイラー閉路を持つ必要条件は全ての頂点について入次数と出次数が同じことと孤立点を除いて連結であることの両方を満たすことです。この条件を満たすグラフをオイラーグラフと言います。
全ての辺が区別できる連結な有向オイラーグラフ $G$ について、以下が成り立ちます。
$G$ に含まれるオイラー閉路の個数は次の式で表せる。
$$\displaystyle c(G,u)\prod_{u\in V} (\text{outdeg}(u) - 1)!$$
ここで、 $c(G,u)$ を $u$ を根とした全域木(全ての辺が $u$ の方向を向いている)の数、$V$ を $G$ の頂点集合、$\text{outdeg}(u)$ を $u$ の出次数とします。また、連結な有向オイラーグラフ $G$ に対する $c(G,u)$ は $u$ をどの頂点としても値が変わらないことが知られています。
以下から引用しました。証明等も以下の解説を読んでください。
これを用いることで、全ての辺を一回ずつ通る walk の数え上げもできます。始点と終点に辺をつないでオイラー閉路の数え上げに帰着すれば良いからです。
コンテスト中に問題を解く上で重要なことは、辺の通る回数の指定された walk の数え上げは解ける可能性が高くて、頂点の通る回数の指定された walk の数え上げは解ける可能性が低いということです。
頂点の通る回数の指定されたときはなんとかして辺の通る回数を指定するように変換した方が良い可能性が高いです。
$(\text{outdeg}(u) - 1)!$ に関しては入力に対して線形で解けることが多いので基本的にここの計算が難しいことは少ないです。$c(G,u)$ は特殊な形をしているときに高速に計算することができます。
[3] 準備
modint や 階乗を事前に前計算しておくライブラリがあると便利です。
検索すればいくらでも出てくると思いますが、自分が今回使ったものを貼っておきます。
modint
struct mint{ static constexpr int m = 998244353; int x; mint() : x(0){} mint(long long x_):x(x_ % m){if (x < 0) x += m;} int val(){return x;} mint &operator+=(mint b){if ((x += b.x) >= m) x -= m; return *this;} mint &operator-=(mint b){if ((x -= b.x) < 0) x += m; return *this;} mint &operator*=(mint b){x= (long long)(x) * b.x % m; return *this;} mint pow(long long e) const { mint r = 1,b =*this; while (e){ if (e & 1) r *= b; b *= b; e >>= 1; } return r; } mint inv(){return pow(m - 2);} mint &operator/=(mint b){return *this *= b.pow(m - 2);} friend mint operator+(mint a, mint b){return a += b;} friend mint operator-(mint a, mint b){return a -= b;} friend mint operator/(mint a, mint b){return a /= b;} friend mint operator*(mint a, mint b){return a *= b;} friend bool operator==(mint a, mint b){return a.x == b.x;} friend bool operator!=(mint a, mint b){return a.x != b.x;} };
Binomial
namespace po167{ template<class T> struct Binomial{ std::vector<T> fact_vec, fact_inv_vec; void extend(int m = -1){ int n = fact_vec.size(); if (m == -1) m = n * 2; if (n >= m) return; fact_vec.resize(m); fact_inv_vec.resize(m); for (int i = n; i < m; i++){ fact_vec[i] = fact_vec[i - 1] * T(i); } fact_inv_vec[m - 1] = T(1) / fact_vec[m - 1]; for (int i = m - 1; i > n; i--){ fact_inv_vec[i - 1] = fact_inv_vec[i] * T(i); } } Binomial(int MAX = 2){ fact_vec.resize(1, T(1)); fact_inv_vec.resize(1, T(1)); extend(MAX + 1); } T fact(int i){ if (i < 0) return 0; while (int(fact_vec.size()) <= i) extend(); return fact_vec[i]; } T invfact(int i){ if (i < 0) return 0; while (int(fact_inv_vec.size()) <= i) extend(); return fact_inv_vec[i]; } T C(int a, int b){ if (a < b || b < 0) return 0; return fact(a) * invfact(b) * invfact(a - b); } T invC(int a, int b){ if (a < b || b < 0) return 0; return fact(b) * fact(a - b) *invfact(a); } T P(int a, int b){ if (a < b || b < 0) return 0; return fact(a) * invfact(a - b); } T inv(int a){ if (a < 0) return inv(-a) * T(-1); if (a == 0) return 1; return fact(a - 1) * invfact(a); } }; }
[4] 次数が少ないグラフ
有向行列木定理が存在します。これは行列木定理とほとんど同じです。行列木定理を知らない人は例えばこのサイトを見てください。
$L_{i,j} = G$ の頂点 $i$ から頂点 $j$ に向かう辺の数に $-1$ をかけたもの
$L_{i,i} = G$ の頂点 $i$ の自己ループを除いた出次数
上記のように定義された行列 $L$ に対して以下が成り立つ。
$$c(u,G) = \text{det}L_{u}$$
ただし、 $L_{u}$ は $L$ から $u$ 行目 $u$ 列目を取り除いたもの
$N\times N$ の行列 $A$ に対する $\text{det} A$ は $O(N^{3})$ で計算できます。
以下は Library Checker へのリンクです
det を 入力で与えられる正整数で割ったあまりを求める問題
今回は mod が問題ごとに固定された素数であるようなものしかないので、上が解ければ十分です。
この後の問題を解くために、有向グラフ $G$ が与えられたとき、そのオイラー路やオイラーパスを数え上げるライブラリを作りました。使い方等は以下の問題例を参考にしてみてください。中身の説明も暇なときにやるかもやらないかも。
オイラー 路/パス 数え上げ
namespace po167{ template<class S> std::vector<S> count_outdeg(std::vector<std::vector<S>> &G){ int N = G.size(); std::vector<S> outdeg(N); for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) outdeg[i] += G[i][j]; return outdeg; } template<class S> std::vector<S> count_indeg(std::vector<std::vector<S>> &G){ int N = G.size(); std::vector<S> indeg(N); for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) indeg[j] += G[i][j]; return indeg; } // O(|M|^{3}) template<class T> T det_matrix(std::vector<std::vector<T>> M){ int N = M.size(); if (N == 0) return 1; if (N == 1) return M[0][0]; if (N == 2) return M[0][0] * M[1][1] - M[0][1] * M[1][0]; T res = 1; for (int i = 0; i < N; i++){ for (int j = i; j < N; j++){ if (M[j][i] != 0){ if (j != i){ swap(M[i], M[j]); res *= -1; } break; } } if (M[i][i] == 0) return 0; res *= M[i][i]; if (i + 1 == N) break; T v = 1 / M[i][i]; for (int j = i + 1; j < N; j++){ T t = M[j][i] * v; for (int k = i; k < N; k++){ M[j][k] -= M[i][k] * t; } } } return res; } template<class T,class S> std::vector<std::vector<T>> Directed_Matrix_tree_Theorem_sub(std::vector<std::vector<S>> &G, int u = 0){ int N = G.size(); std::vector L(N - 1, std::vector<T>(N - 1, 0)); std::vector<S> outdeg(N); for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) if (i != j) outdeg[i] += G[i][j]; for (int i = 0; i < N; i++) for (int j = 0; j < N; j++){ if (i == u || j == u) continue; int a = i, b = j; if (u < i) a--; if (u < j) b--; if (i == j) L[a][b] = outdeg[i]; else L[a][b] -= G[i][j]; } return L; } template<class T,class S> T Directed_Matrix_tree_Theorem(std::vector<std::vector<S>> &G, int u = 0){ if (int(G.size()) == 1) return T(1); return det_matrix(Directed_Matrix_tree_Theorem_sub<T>(G, u)); } // s.t // forall i,j // 0 <= G[i][j] // if not connected // use remove iso template<class T,class S> std::pair<T, std::vector<std::vector<T>>> Count_Euler_Circuit_sub(std::vector<std::vector<S>> G, std::vector<T> &fact_base, bool fact_in = true){ int N = G.size(); if (N == 0) return {1, {{1}}}; std::vector<S> outdeg = count_outdeg(G); std::vector<S> indeg = count_indeg(G); // indeg == outdeg ? for (int i = 0; i < N; i++) if (indeg[i] != outdeg[i]) return {0, {{1}}}; // connected ? std::vector<bool> seen(N); std::vector<int> order={0}; seen[0] = 1; for (int i = 0; i < N; i++){ if (i == int(order.size())) return {0, {{1}}}; for (int j = 0; j < N; j++){ if (G[order[i]][j] != 0 && !seen[j]){ seen[j] = 1; order.push_back(j); } } } T res = 1; if (fact_in) for (int i = 0; i < N; i++) res *= fact_base[outdeg[i] - 1]; return {res ,Directed_Matrix_tree_Theorem_sub<T, S>(G)}; } template<class T,class S> T Count_Euler_Circuit(std::vector<std::vector<S>> G, std::vector<T> &fact_base, bool fact_in = true){ auto tmp = Count_Euler_Circuit_sub(G, fact_base, fact_in); if (tmp.first == 0) return T(0); return tmp.first * det_matrix(tmp.second); } template<class S> std::vector<std::vector<S>> remove_isolated_vertex( std::vector<std::vector<S>> G ){ int N = G.size(); std::vector<int> seen(N, -1); for (int i = 0; i < N; i++) for (int j = 0; j < N; j++){ if (G[i][j] != 0) seen[i] = 1, seen[j] = 1; } int ind = 0; for (int i = 0; i < N; i++) if (seen[i] == 1) seen[i] = ind++; std::vector res(ind, std::vector<S>(ind)); for (int i = 0; i < N; i++) if (seen[i] != -1){ for (int j = 0; j < N; j++) if (seen[j] != -1){ res[seen[i]][seen[j]] = G[i][j]; } } return res; } // s.t // forall i,j // 0 <= G[i][j] // if not connected // use remove iso template<class T,class S> std::pair<T, std::vector<std::vector<T>>> Count_Eulerian_Trail_sub(std::vector<std::vector<S>> G,std::vector<T> &fact_base, bool fact_in = true){ bool ok = 1; int N = G.size(); if (N == 0) return {1, {{1}}}; std::vector<S> outdeg = count_outdeg(G), indeg = count_indeg(G); S sum = 0; for (int i = 0; i < N; i++){ sum += outdeg[i]; } int st = -1, ed = -1; for (int i = 0; i < N; i++) if (outdeg[i] != indeg[i]){ if (std::abs(outdeg[i] - indeg[i]) > 1){ ok = 0; break; } if (outdeg[i] > indeg[i]){ if (ed != -1) ok = 0; ed = i; } else { if (st != -1) ok = 0; st = i; } } if (!ok) return {0, {{1}}}; if ((st == -1) ^ (ed == -1)) return {0, {{1}}}; if (st == -1){ auto tmp = Count_Euler_Circuit_sub(G, fact_base, fact_in); return {T(sum) * tmp.first, tmp.second}; } G[st][ed]++; return Count_Euler_Circuit_sub(G, fact_base, fact_in); } template<class T,class S> T Count_Eulerian_Trail(std::vector<std::vector<S>> G,std::vector<T> &fact_base, bool fact_in = true){ auto tmp = Count_Eulerian_Trail_sub<T, S>(G, fact_base, fact_in); if (tmp.first == 0) return 0; return tmp.first * det_matrix(tmp.second); } }
与えるグラフ $G$ は $N\times N$ の二次元配列であることが想定されていて、$N = 0$ のときは $1$ が返ってくるようになっているはずです。
また、長さが $\text{outdeg}$ 以上の配列で $i$ 番目の要素が $i$ の階乗であるようなものを引数に入れる必要があります。これが嫌なのであれば、その配列を入れる部分に適当な mint の配列を入れ、その後に $0$ を入れる必要があります。 すると、 $c(G,u)$ の値のみ返ってきます。
ABC 336 G - 16 Integers
問題
長さ $16$ で和が $N$ の非負整数列 $X=(X_{0},X_{1},\dots X_{15})$ が与えられます。
長さ $3+N$ の数列 $A$ で、以下の条件を全て満たすものの場合の数を $998244353$ で割ったあまりを求めてください。
- $A$ の任意の要素は $0$ もしくは $1$
- $A_{i} \times 8 + A_{i + 1} \times 4 + A_{i + 2} \times 2 + A_{i + 3} = j$ を満たす $i$ の数は $X_{j}$ 個
$1\leq N\leq 10^{6}$
解法
$(i,j,k,l)$ から $(j,k,l,m)$ へ辺をはると、頂点の通る回数が指定された walk の数え上げになって、解くことが難しいです。$(i,j,k,l)$ が存在するということは、連続した $3$ 文字を見て $(i,j,k)$ から $(j,k,l)$ になる部分があるということです。
そこで、 $(i,j,k,l)$ を頂点にするのではなく辺だと見ると、$(i,j,k)$ から $(j,k,l)$ へ結ぶ辺であると言えます。
つまり、以下の $8$ 頂点のグラフのオイラー路を数え上げれば良いです。
- $0\leq i < 16$ について、 頂点 $\lfloor i\rfloor$ から 頂点 $i \pmod{8}$ への辺が $X_{i}$ 本ある
実装例
long long solve_abc336_g(std::vector<int> X,int bit = 4){ po167::Binomial<mint> C(1001001); if (bit == 1){ return (C.fact(X[0] + X[1]) * C.invfact(X[0]) * C.invfact(X[1])).val(); } int len = (1 << (bit - 1)); std::vector G(len, std::vector<int>(len, 0)); for(int i = 0; i < len * 2; i++){ G[i / 2][i & (len - 1)] = X[i]; } mint ans = po167::Count_Eulerian_Trail(po167::remove_isolated_vertex(G), C.fact_vec); for(auto x: X) ans *= C.invfact(x); return ans.val(); }
提出例
CF Round 925 - One-Dimensional Puzzle
問題
以下の図形がそれぞれ $c_{i}$ 個あります。
これらを全て横一列に凸凹が合うように並べる方法を $998244353$ で割ったあまりを求めてください。
$1\leq \text{testcase}\leq 2\times 10^{5}$
$0\leq c_{i}\leq 10^{6}$
解法
先ほどの ABC の問題の bit 数が減っただけです。
実装例
po167::Binomial<mint> table(4004004); long long solve_cf_1931_g(std::vector<int> C){ std::vector G(2, std::vector<int> (2)); G[0][1] = C[0]; G[1][0] = C[1]; G[1][1] = C[2]; G[0][0] = C[3]; if (C[0] + C[1] + C[2] + C[3] == 0) return 1; mint ans = po167::Count_Eulerian_Trail(po167::remove_isolated_vertex(G), table.fact_vec); for (auto x: C) ans *= table.invfact(x); return ans.val(); }
先ほどの ABC336 の問題を解くのに使った関数をそのまま適用した提出例
AGC 051 D - C4
問題
以下の無向グラフにおいて、 $S$ から $S$ への walk で、辺の通る回数がそれぞれ $a,b,c,d$ 回であるようなものの場合の数を $998244353$ で割ったあまりを求めてください。
$0\leq a,b,c,d\leq 5\times 10^{5}$
解法
まず、$a,b,c,d$ の偶奇が一致していないものは $0$ 回です。
それぞれの辺についての時計回りに通る回数と反時計回りに通る回数の差は等しいです。
よって、$1$ つの辺について、時計回りに通る回数を決めると、全ての辺について、時計回りに通る回数と反時計回りに通る回数が決まります。
ある辺の時計回りに通るを決め打って、それぞれの回数について、有向グラフのオイラー閉路の数え上げができれば答えが求まります。
実装例
long long solve_agc051_d(std::vector<int> A){ const int LIM = 500500; po167::Binomial<mint> C(LIM); int M = A[0]; for(int i = 0; i < 4; i++){ if((A[i] + A[(i + 1) % 4]) & 1){ return 0; } M = std::min(M, A[i]); } mint ans = 0; for(int d = -((M + 1) / 2); d <= M / 2; d++){ std::vector G(4, std::vector<int>(4)); mint tmp = 1; for(int i = 0; i < 4; i++){ int a = i; int b = (i + 1) % 4; G[a][b] = (A[i] + 1) / 2 + d; G[b][a] = A[i] / 2 - d; tmp *= C.invfact(G[a][b]); tmp *= C.invfact(G[b][a]); } ans += tmp * po167::Count_Euler_Circuit(G, C.fact_vec); } ans *= ((A[0] + A[3]) / 2); return ans.val(); }
提出例
[5] 線グラフ
$N$ 頂点のグラフで、$i$ から $i+1$ への辺と $i+1$ から $i$ への辺がそれぞれ $A_{i}$ 本あるようなグラフの有効全域木の場合の数は $\prod_{i=1}^{N-1} A_{i}$ で表せます。また、頂点の index の差が $1$ でない辺が $1$ 本あるようなオイラーグラフにおいても、その特殊な辺の始点を $u$ として $c(G,u)$ を計算するのは簡単です。
この線グラフはある整数から始めて $+1,-1$ を繰り返すような設定の問題から生成することができます。
ARC 146 E - Simple Speed
問題
長さ $N$ の正整数列 $A$ が与えられます。
以下の条件を満たす長さ $\sum A$ の正整数列 $B$ の場合の数を $998244353$ で割ったあまりを求めてください。
- $B$ の中にある $i$ の数は $A_{i}$ 個
- $B$ の隣り合う要素の差の絶対値は全て $1$
$1\leq N\leq 2\times 10^{5}$
$1\leq A_{i}\leq 2\times 10^{5}$
解法
このように、終点から始点への辺を追加することで、オイラー路の数え上げへと帰着できます。
$i$ の小さい順に dp をして求められます。このとき情報として、$i$ から $i+1$ への辺は何本あるのかと、 $i+1$ から $i$ への辺は何本あるのかと、始点、終点それぞれについて、右側左側どちらにあるのか、というものを持っておく必要があります。
この状態の量が実は線形になるので、この dp は $O(N)$ で可能です。
実装例
int solve_arc146_e(int N, std::vector<int> A){ po167::Binomial<mint> B(200200); std::map<std::pair<int, int>, mint> dp; dp[{0, 0}] = 1; for (int i = 0; i < N; i++){ std::map<std::pair<int, int>, mint> n_dp; for (auto tmp: dp){ int l_indeg = tmp.first.first; int l_outdeg = l_indeg; int t = tmp.first.second; if (t == 1) l_outdeg--; if (t == 2) l_outdeg++; if (l_indeg < 0 || l_outdeg < 0) continue; if (i != 0){ // 連結性を保ち続ける if (l_indeg == 0 && t % 3 == 0) continue; // 同じ辺を区別しないようにする tmp.second *= B.invfact(l_indeg); tmp.second *= B.invfact(l_outdeg); } for (int k = 0; k < 4; k++) if ((k & t) == 0 ){ int r_indeg = A[i] - l_indeg; int r_outdeg = A[i] - l_outdeg; if(k & 1) r_indeg--; if(k & 2) r_outdeg--; mint val = tmp.second; // 全域木の数だけ乗算する if ((k & 2) == 0){ if (t & 2) val *= l_outdeg; else val *= r_outdeg; } // 自分の出次は相手の入次 if (val != 0) n_dp[{r_outdeg, k + t}] += val; } } std::swap(dp, n_dp); } mint ans = dp[{0, 3}]; // prod (outdeg - 1)! for (int i = 0; i < N; i++){ ans *= B.fact(A[i] - 1); } return ans.val(); }
提出例
CF Hello 2024 E - Counting Prefixes
問題
長さ $N$ の $-1,+1$ からなる数列 $A=(A_{1}\dots A_{N})$ に対して、 $\displaystyle S_{i} = \sum_{j = 1}^{i} A_{j}$ とし、$(S_{1},\dots S_{N})$ を昇順ソートしたものを $f(A)$ とします。
広義単調増加な長さ $N$ の数列 $p$ が与えられるので、 $f(A)=p$ となる $A$ の場合の数を $998244353$ で割ったあまりを求めてください。
$N\leq 5000$
解法
先ほどの ARC の問題の始点が $0$ で固定されているだけです。
始点の更新だけ決められた位置で行えば、同じように解くことができます。
実装例
int solve_cf_1919_E(int n, std::vector<int> p){ po167::Binomial<mint> B(n + 1); p.push_back(0); std::sort(p.begin(), p.end()); std::vector<int> A = {1}; int Z = -p[0]; for (int i = 1; i <= n; i++){ if(p[i - 1] == p[i]){ A[int(A.size()) - 1]++; } else if(p[i - 1] + 1 == p[i]){ A.push_back(1); } else{ return 0; } } std::map<std::tuple<int, int, int>, mint> dp; dp[{0, 0, 0}] = 1; for (int i = 0; i < int(A.size()); i++){ std::map<std::tuple<int, int, int>, mint> n_dp; for (auto tmp: dp){ auto[l_indeg, l_outdeg, t] = tmp.first; mint val = tmp.second; if (l_indeg < 0 || l_outdeg < 0) continue; if (i != 0){ if(l_indeg == 0 && l_outdeg == 0) continue; val *= B.invfact(l_indeg); val *= B.invfact(l_outdeg); } for (int k = 0; k < 2; k++) if (k + t != 2){ int r_indeg = A[i] - l_indeg; int r_outdeg = A[i] - l_outdeg; mint n_val = val; if (i == Z){ r_outdeg--; } else if (i < Z){ n_val *= r_outdeg; } else{ n_val *= l_outdeg; } if (k == 1){ r_indeg--; } n_dp[{r_outdeg, r_indeg, k + t}] += n_val; } } std::swap(n_dp, dp); } mint ans = dp[{0, 0, 1}]; for(auto x: A) ans *= B.fact(x - 1); return ans.val(); }
提出例
Submission #247368236 - Codeforces
ARC 117 E - Zero-Sum Ranges 2
問題
$N,K$ が与えられます。
長さ $2N$ の数列 $A$ で以下の条件を全て満たすものの場合の数を求めてください。
$A$ は $N$ 個の $+1$ と $N$ 個の $-1$ からなる
$A_{l} + A_{l+1} + \cdots A_{r} = 0$ となる $(l,r)$ の組の数がちょうど $K$ 個
$1\leq N\leq 30$ $1\leq K\leq N^{2}$
解法
頂点 $0$ から始めて、頂点 $0$ で終わる長さ $2N$ の walk で、頂点 $i$ を通過した回数を $a_{i}$ 回としたとき、$\sum \dfrac{a_{i}(a_{i}-1)}{2} =K$ となる walk の場合の数を数えれば良いです。
端から頂点を繋いでいって、最後頂点 $0$ で繋げるような dp をします。
$dp[i][j][k]$ を、今まで $i$ 本辺を使って、 直前が $j$ 本で、 $k$ 個 zero-sum になる範囲があるものの場合の数に関する値とします。この場合の数はどのように更新すれば良いのでしょうか。
例えば写真下部について、このときの場合の数について考えます。
- $c(G,u)$ に関する寄与 : $a\times b\times c$
- $\text{outdeg}$ に関する寄与 : $(a-1)!\times (a+b-1)! \times (b+c-1)!$
- 辺を区別しないことに関する寄与 : $\dfrac{1}{(a!)^{2}(b!)^{2}(c!)^{2}}$
これらの積は以下のように表現できます。
$$\dfrac{(a+0-1)!}{(0!)((a-1)!)}\times \dfrac{(a+b-1)!}{(a!)((b-1)!)}\times\dfrac{(b+c-1)!}{(b!)((c-1)!)}\times \dfrac{1}{c!}$$
よって、 $dp[i][j][k]$ には、この積に $j$ の階乗をかけたものを持っておきます。そして、この状態に $1$ つの頂点と $l$ 本の辺を加える操作は以下の dp の更新に対応づけられます。
$$dp[i+l][l][k + \frac{(j+l)(j+l-{1})}{2}] += dp[i][j][k]\times \frac{(j+l-{1})!}{(j!)((l-{1})!)}$$
このように更新することで、除算を発生させることなく、常に整数で状態数を持つことができます。
最後に条件を満たすような二つの組をマージし、それらの積の和を取れば良いです。
二項係数を事前に計算しておけば、遷移は $O(1)$ で行えます。dp の空間計算量が $O(N^{4})$ で、 遷移が高々 $N$ 通りなので、時間計算量は $O(N^{5})$ です。ただし、無駄な状態が多いので、実際の実行時間は短いです。
$\dfrac{60!}{30!30!}\leq 10^{18}$ であることから、答えは long long におさまります。よって、modint で求めて復元などをしなくても良いです。
実装例
long long solve_arc117_e(int N, int K){ std::vector C(N + 1, std::vector<long long>(N + 1)); C[0][0] = 1; for (int i = 0; i < N; i++) for (int j = 0; j <= i; j++){ C[i + 1][j] += C[i][j]; C[i + 1][j + 1] += C[i][j]; } // 次数の総和、一個前の次数、(l, r) の組み合わせの数 std::vector dp(N + 1, std::vector(N + 1, std::vector<long long>(K + 1))); dp[0][0][0] = 1; for (int sum = 0; sum < N; sum ++){ for (int fr = 0; fr <= sum; fr++){ for (int count = 0;count <= K; count++) if (dp[sum][fr][count] != 0){ for (int k = 1; k + sum <= N; k++){ long long val = dp[sum][fr][count]; val *= C[fr + k - 1][k - 1]; int n_sum = sum + k; int n_count = count + (fr + k) * (fr + k - 1) / 2; if(n_count > K) break; dp[n_sum][k][n_count] += val; } } } } long long ans = 0; for (int sum = 0; sum <= N; sum++){ for (int l = 0; l <= sum; l++){ for (int count = 0; count < K; count++){ for (int r = 0; sum + r <= N; r++){ int r_c = K - count - (l + r) * (l + r + 1) / 2; if (r_c < 0) break; ans += dp[sum][l][count] * dp[N - sum][r][r_c] * C[l + r][l]; } } } } return ans; }
提出例
[6] スパースな行列
Library Checker に以下のような問題があります。
スパースな行列の行列式を mod 998244353 で求める問題
Black Box Linear Algebra というアルゴリズムを用いてスパースな行列に対する行列式を高速に求めることができます。
詳細はこちら
Nyaan さんによる実装例はこちら
上記を参考にした自分の Black Box Linear Algebra の実装例
namespace po167{ std::random_device po_seed_gen; std::mt19937 po_random_gen(po_seed_gen()); long long randLong(long long l = 0, long long r = (1ll << 62)){ return std::uniform_int_distribution<long long>(l, r - 1)(po_random_gen); } template<class T> std::vector<T> rand_vec(int sz, long long l =0, long long r = (1ll << 62)){ std::vector<T> res(sz); for(auto &x: res) x = T(randLong(l, r)); return res; } //https://ikatakos.com/pot/programming_algorithm/number_theory/barlekamp_massey template<class T> std::vector<T> Berlekamp_Massey(std::vector<T> s){ std::vector<T> C={1}, B={1}; int L = 0; int m = 1; T b = 1; for (int n = 0; n < int(s.size()); n++){ T d = 0; for (int j = 0; j < L + 1; j++){ d += C[j] * s[n - j]; } if (d == 0){ m++; continue; } T iv = d / b; if (2 * L <= n){ auto tC = C; while (int(C.size()) <= n + 1 - L) C.push_back(0); for (int i = 0; i < int(B.size()); i++) C[i + m] -= iv * B[i]; B = tC; L = n + 1 - L; b = d; m = 1; } else{ for(int i = 0; i < int(B.size()); i++){ C[i + m] -= iv * B[i]; } m++; } } return C; } // https://nyaannyaan.github.io/library/matrix/black-box-linear-algebra.hpp template<class T> T det_sparse_matrix(std::vector<std::vector<T>> A, T mj = 0){ int n = A.size(); std::vector<T> D; struct pos_mat{ int x; int y; T val; }; std::vector<pos_mat> p; for (int i = 0; i < n; i++){ for(int j = 0; j < n; j++) if (A[i][j] != mj) { p.push_back({i, j, A[i][j] - mj}); } } while (true){ while (true){ D = rand_vec<T>(n); bool ok = 1; for (auto x: D) if (x == 0) ok = 0; if (ok) break; } std::vector<pos_mat> AD = p; for (int i = 0; i < int(AD.size()); i++) AD[i].val *= D[AD[i].y]; std::vector<T> u = rand_vec<T>(n), v = rand_vec<T>(n); std::vector<T> b(n); std::vector<T> a(2 * n + 1); b = u; for (int i = 0; i < 2 * n + 1; i++){ T sum = 0; for (int j = 0; j < n; j++){ sum += u[j] * D[j]; a[i] += u[j] * v[j]; } sum *= mj; for (int j = 0; j < n; j++){ u[j] = sum; } for (pos_mat tmp: AD){ u[tmp.x] += tmp.val * b[tmp.y]; } b = u; } auto mp = Berlekamp_Massey(a); if (mp.back() == 0) return 0; if (int(mp.size()) != n + 1) continue; T res = mp.back(); if (n & 1) res *= -1; T tmp = 1; for (auto d: D) tmp *= d; return res / tmp; } } }
Library Checker への提出例
行列 $A$ についてランダムな対角行列 $D$ をとって $A\leftarrow AD$ と置換した後、任意のベクトル $u$ に対して $u A$ を求める時間計算量を $T$ としたとき、 $O(N^{2} + NT)$ で求められます。
以上より、非零な要素の数を $K$ として $O(N^{2}+NK)$ です。
やり方は上記の記事等を見てください。 BM(Berlekamp_Massey) 法が関係しています。
また、$O(N^{3})$ の解法をそのままに、行列を事前にシャッフルすることで高速に求める方針もあります。
yukicoder No.310 - 2文字しりとり
問題
$N$ 頂点で全ての $(i,j)$ に対して、 $i$ から $j$ へちょうど $1$ 本辺がはられている $N^{2}$ 辺の有効グラフから、入力で与えられる $K$ 本の辺を取り除いた有効グラフについて、オイラー路の場合の数を $10^{9}+7$ で割ったあまりを求めてください。
$1\leq N\leq 4000$ $1\leq K\leq 4000$
解法
ほとんどの頂点間に $1$ 本の辺がはられているため、 $c(G,u)$ を計算するために行列式を求めたい行列のうち、 $-1$ でない要素は高々 $K + N$ 個です。
よって、Black Box Linear Algebra を用いて、 $O(N(K+N))$ で求めることができます。
実装例
int solve_yuki_0310(int N, int M, std::vector<int> A, std::vector<int> B){ if (N * N == M) return 1; std::vector G(N ,std::vector<int>(N, 1)); for (int i = 0; i < M; i++){ G[A[i] - 1][B[i] - 1] = 0; } po167::Binomial<mint> bin(N + 100); auto tmp = po167::Count_Eulerian_Trail_sub(po167::remove_isolated_vertex(G), bin.fact_vec); return (tmp.first * po167::det_sparse_matrix<mint>(tmp.second , -1)).val(); }
提出例
終わり
最初に紹介した ABC の問題の解説を見るまで BEST 定理を自分は知らなかったのですが、どうやらかなり前から名が知られていたみたいです。
せっかくなのでライブラリ化しましたが、ここまで様々なバージョンが出ていると、もう出題されないのではという気がしてしまいました。
