yukicoder No.590 Replacement
実装が結構大変。
https://yukicoder.me/problems/no/590
解法
順列なので辺 $ i \to A_i$ のグラフを考えよう。このようなグラフはサイクルの集まりになることに注意したい。
\( (x,y) \to (i,i) \) と考えるのではなく、逆向きの操作を考えて、\( (i,i) \to (x,y) \) となる対 \( (x,y) \) の数を数えることにする。
もし \( (i,i) \to (j,j) \) となるような \(j\) が存在しなければ、\( (i,i) \to (x,y) \) となる対 \((x,y)\) の個数はサイクル長の LCM である。個数が $n$ ならコストは $n(n-1)/2$ であるから、個数を求めることが本質である。存在するときが難しい。
まず問題になるのは \( (i,i) \) と \( (j,j) \) の連結性判定である。よって次は連結性判定について考える。
辺 $A_i \to i$ で構成されたグラフを $A$、辺 $B_i \to i$ で構成されたグラフを $B$ とする。グラフ $G$ において、頂点 $i$ が属する連結成分の代表頂点を $[i]_G$ と書くことにすると、少なくとも連結であるためには $([i]_A,[i]_B) = ([j]_A,[j]_B)$ である必要がある。$(i,i)$ と $(j,j)$ が連結である条件はもう一つある。各サイクルに 0 から順に通し番号を割り振る。$(i)_G$ を頂点 $i$ に振られた番号だとすると $(i)_B-(i)_A \equiv (j)_B-(j)_A \pmod{\gcd(n, m)}$ が必要である。ここで $n$, $ m $ はそれぞれ頂点 $i$, $j$ が属する連結成分の頂点数である。理由は図を見よ。
$(A,A)\to(B,B)\to(C,C)\to(A,A)\to\cdots$ となり、ラベル差が不変になっていることが分かる。またラベル差が等しい組は、一方から一方へ必ず到達できることも確認できる(大きなグループでみると互いに素であるため)。
条件をまとめると、
- $([i]_A,[i]_B) = ([j]_A,[j]_B)$
- $(i)_B-(i)_A \equiv (j)_B-(j)_A \pmod{\gcd(n, m)}$
となる。まずこの条件を用いて $(i,i)$ を分類する。
次に $(i,i) \to (j,j)$ に何手で到達できるかを求める。いま $(x_1,x_1), (x_2,x_2), \ldots, (x_s, x_s)$ が連結だったとする。中国剰余定理により $(x_i, x_i)$ に対応する $z_i$ が存在し、$z_i$ でソートすることで、$(x_k,x_k)$ の次の頂点対が分かる。さらに $z$ の差を見ることで距離も同時に分かる。以上より解くことが出来た。
#include <iostream> #include <algorithm> #include <vector> #include <map> #include <set> using namespace std; const long long mod = 1e9 + 7; long long gcd(long long x, long long y) { if (y == 0) return x; return gcd(y, x % y); } long long lcm(long long x, long long y) { return x / gcd(x, y) * y; } vector<vector<int>> cycle(vector<int> a) { const int n = a.size(); vector<vector<int>> ret(n); vector<bool> vis(n); for (int i = 0; i < n; i++) { if (vis[i]) continue; for (int k = i; !vis[k]; k = a[k]) { ret[i].push_back(k); vis[k] = true; } } return ret; } long long modinv(long long a, long long m) { long long b = m, x = 1, y = 0; while (b != 0) { long long q = a / b; a -= b * q; x -= y * q; std::swap(a, b); std::swap(x, y); } return x < 0 ? x + m : x; } long long crt(long long a1, long long m1, long long a2, long long m2) { long long v = (a2 - a1) * modinv(m1, m2) % m2; if (v < 0) v += m2; return a1 + v * m1; } map<long long, int> factors(long long n) { map<long long, int> ret; for (int i = 2; i * i <= n; i++) { while (n % i == 0) { ret[i]++; n /= i; } } if (n != 1) ret[n] = 1; return ret; } int power(int a, int b) { int ret = 1; for (int i = 0; i < b; i++) { ret *= a; } return ret; } long long crt_g(long long a1, long long m1, long long a2, long long m2) { auto f1 = factors(m1); auto f2 = factors(m2); int M1 = 1; int M2 = 1; set<int> st; for (auto kv : f1) st.insert(kv.first); for (auto kv : f2) st.insert(kv.first); for (int k : st) { if (f1[k] >= f2[k]) { M1 *= power(k, f1[k]); } else { M2 *= power(k, f2[k]); } } return crt(a1 % M1, M1, a2 % M2, M2); } int main() { int n; cin >> n; vector<int> a(n), b(n); for (int i = 0; i < n; i++) { scanf("%d", &a[i]); a[i]--; } for (int i = 0; i < n; i++) { scanf("%d", &b[i]); b[i]--; } vector<vector<int>> cycleA = cycle(a); vector<vector<int>> cycleB = cycle(b); vector<int> rootA(n), rootB(n); vector<int> xa(n), xb(n); for (int i = 0; i < n; i++) { for (int j = 0; j < cycleA[i].size(); j++) { xa[cycleA[i][j]] = j; rootA[cycleA[i][j]] = i; } for (int j = 0; j < cycleB[i].size(); j++) { xb[cycleB[i][j]] = j; rootB[cycleB[i][j]] = i; } } map<pair<int, int>, vector<int>> mp; for (int i = 0; i < n; i++) { mp[make_pair(rootA[i], rootB[i])].push_back(i); } long long ans = 0; for (auto kv : mp) { int A = kv.first.first; int B = kv.first.second; int lenA = cycleA[A].size(); int lenB = cycleB[B].size(); long long G = gcd(lenA, lenB); long long L = lcm(lenA, lenB); map<int, vector<int>> mp2; for (int i : kv.second) { int d = (xb[i] - xa[i]) % G; if (d < 0) d += G; mp2[d].push_back(i); } for (auto kv2 : mp2) { vector<long long> xs; for (int i : kv2.second) { int x = xa[i]; int y = xb[i] - kv2.first; if (y < 0) y += lenB; xs.push_back(crt_g(x, lenA, y, lenB)); } sort(xs.begin(), xs.end()); xs.push_back(xs[0] + L); for (int i = 0; i + 1 < xs.size(); i++) { long long d = (xs[i + 1] - xs[i]) % mod; ans += d * (d + mod - 1) % mod * ((mod + 1) / 2); ans %= mod; } } } cout << ans << endl; }