AOJ 1158 - ICPC: Intelligent Congruent Partition of Chocolate
解法
\(_{36}C_{18}\) は流石に間に合わない。
しかし、2つの領域を合同かつ連結になるように選んでいけば、探索空間はそんなに大きくないように思える（実際大きくない）。
なので、今探索中の領域1に連結なチョコレートを付け加えたときに、領域2にも対応する位置のチョコレートを使うと確定させるような探索が良さそう。
そのためにまず、領域2が領域1に対してどれだけ回転しているか＆裏返しかを決めておく。
このもとで、一番左上のチョコレートを領域1として使うと確定させ、そのチョコレートが領域2においてどこに対応するかを全部試す。
すると、チョコレートを領域1に付け加えたとき、それが領域2のどこに対応するかがわかるようになる。
領域1にチョコレートを付け加えるときに、連結になるようにしておけば自動的に領域2も連結になるのでOK。
途中でバッティングしたり範囲外に出たりせず全部使い切れれば YES とすればよい。
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = long long; constexpr int dy[4] = {0, 1, 0, -1}; constexpr int dx[4] = {1, 0, -1, 0}; int main() { int w, h; while(cin >> w >> h, w) { vector<vector<int>> r(h, vector<int>(w)), idx(h, vector<int>(w, -1)); vector<int> xs, ys; for(int i = 0; i < h; ++i) { for(int j = 0; j < w; ++j) { cin >> r[i][j]; if(r[i][j] == 1) { idx[i][j] = ys.size(); ys.push_back(i); xs.push_back(j); } } } const int n = xs.size(); if(n & 1) { cout << "NO" << endl; continue; } auto in_range = [&] (int y, int x) { return 0 <= y && y < h && 0 <= x && x < w; }; auto calc_dy = [] (int mir, int rot, int s_dy, int s_dx) { if(rot == 0) return s_dy; else if(rot == 1) return (mir ? -s_dx : s_dx); else if(rot == 2) return -s_dy; else return (mir ? s_dx : -s_dx); }; auto calc_dx = [] (int mir, int rot, int s_dy, int s_dx) { if(rot == 0) return (mir ? -s_dx : s_dx); else if(rot == 1) return s_dy; else if(rot == 2) return (mir ? s_dx : -s_dx); else return -s_dy; }; auto solve = [&] () { for(int j = 1; j < n; ++j) { for(int rot = 0; rot < 4; ++rot) { for(int mir = 0; mir < 2; ++mir) { unordered_set<ll> vis; stack<pair<ll, ll>> st; vis.insert(1); st.emplace(1, 1LL << j); while(!st.empty()) { const ll used1 = st.top().first, used2 = st.top().second; if(__builtin_popcountll(used1) * 2 == n) { return true; } st.pop(); for(int k = 0; k < n; ++k) { if(!(used1 & (1LL << k))) continue; for(int d = 0; d < 4; ++d) { const int ny1 = ys[k] + dy[d], nx1 = xs[k] + dx[d]; const int ny2 = ys[j] + calc_dy(mir, rot, ny1 - ys[0], nx1 - xs[0]); const int nx2 = xs[j] + calc_dx(mir, rot, ny1 - ys[0], nx1 - xs[0]); if(!in_range(ny1, nx1) || idx[ny1][nx1] == -1) continue; if(!in_range(ny2, nx2) || idx[ny2][nx2] == -1) continue; const int id1 = idx[ny1][nx1], id2 = idx[ny2][nx2]; if((used1 | used2) & (1LL << id1) || (used1 | used2) & (1LL << id2) || id1 == id2) continue; if(vis.count(used1 | (1LL << id1))) continue; vis.insert(used1 | (1LL << id1)); st.emplace(used1 | (1LL << id1), used2 | (1LL << id2)); } } } } } } return false; }; cout << (solve() ? "YES" : "NO") << endl; } }
感想
国内予選ならではって感じの問題だった。
AOJ 2682 - Polygon Guards
解法
答えは 9 以下なので全探索で間に合う（えぇ…）。
前処理で、ある点に配置したときにどこが見えるようになるかを bit で管理すると判定が O(1) になる。
視線の線分が多角形内部に全部入っているかどうかの判定は適当にやるとハマる。
まず、その線分と多角形の交点を全部求める。
線分の両端から線分の内側にちょっと進んだところが多角形の内部にあるかをまず判定。
交点については、線分の方向ベクトルの2方向にちょっと進んだところが多角形内部にあるかで判定。
この2つをすべてクリアした場合のみ、その点が見える。
AOJ でも 2sec ぐらいかかったので本番で出たらもうちょっと高速化しないとダメかも。
ソースコード
#include <bits/stdc++.h> using namespace std; // ライブラリ部分は省略 int dfs(polygon const& ps, int i, ll cur_vis, ll used, vector<ll> const& vis) { const int n = ps.size(); if(cur_vis == (1LL << n) - 1) return __builtin_popcountll(used); if(i == n) return 9; int res = dfs(ps, i + 1, cur_vis, used, vis); if(__builtin_popcountll(used) < 8) { res = min(res, dfs(ps, i + 1, cur_vis | vis[i], used | (1LL << i), vis)); } return res; } int main() { int n; cin >> n; polygon ps(n); for(int i = 0; i < n; ++i) { int x, y; cin >> x >> y; ps[i] = point(x, y); } vector<ll> vis(n); for(int i = 0; i < n; ++i) { vis[i] = 1LL << i; for(int j = 0; j < n; ++j) { if(i == j) continue; segment s(ps[i], ps[j]); bool check = true; for(int k = 0; k < n; ++k) { segment t(ps[k], ps[(k + 1) % n]); if(!isis_ss(s, t)) continue; const auto delta = (ps[j] - ps[i]) / abs(ps[i] - ps[j]) * 0.001; for(auto p : is_ss(s, t)) { if(abs(s.a - p) > eps && abs(s.b - p) > eps) { check &= is_in_polygon(ps, p + delta) != 2 && is_in_polygon(ps, p - delta) != 2; } } check &= is_in_polygon(ps, ps[i] + delta) != 2; check &= is_in_polygon(ps, ps[j] - delta) != 2; } if(check) { vis[i] |= 1LL << j; } } } cout << dfs(ps, 0, 0, 0, vis) << endl; }
感想
この問題みたいな多角形だと、n / 4 以下でできることが示せるらしい。
一般の多角形なら n / 3 とのこと。
たしかに多角形を三角形 or 四角形で分割することを考えるとそんな気もする。
AOJ 2724 - Laser Cutter
解法
結論を言えば、最小重み二部マッチングに帰着できる。
線分の交点と端点が頂点になった有向グラフを考えると、考えるべき問題は「いくつか辺を追加することで、最小重みのオイラーグラフを作る」になる。
辺を追加するとしたら、線分の端点同士（終点 -> 始点) を考えるだけで良い。
実際、端点以外での交点は、入次数と出次数が等しくなるようにしかならないためである。
したがって、端点同士でマッチングを取れば良いことがわかる。
（一応言っておくと、交点を経由するのは回り道でしかないのでありえない。）
最小重み二部マッチングは最小費用流やハンガリアン法で解けて、O(N^3logN) か O(N^3) となる。
最小費用流で解くときは、重みが整数値じゃないので、ポテンシャルの判定式に eps をつけないとバグるから気をつけること。
ソースコード
#include <bits/stdc++.h> using namespace std; using pii = pair<int, int>; constexpr double eps = 1e-10; struct edge { int to, rev, cap; double cost; edge(int t, int c, double ct, int r) : to(t), rev(r), cap(c), cost(ct) {} }; using graph = vector<vector<edge>>; void add_edge(graph& g, int from, int to, int cap, double cost) { g[from].emplace_back(to, cap, cost, g[to].size()); g[to].emplace_back(from, 0, -cost, g[from].size() - 1); } double min_cost_flow(graph& g, int s, int t, int f) { using P = pair<double, int>; const double inf = 1e18; double res = 0; vector<double > h(g.size()), dist(g.size()); vector<int> prevv(g.size()), preve(g.size()); while(f > 0) { priority_queue<P, vector<P>, greater<>> que; fill(begin(dist), end(dist), inf); dist[s] = 0; que.emplace(0, s); while(!que.empty()) { const auto cur_d = que.top().first; const int v = que.top().second; que.pop(); if(dist[v] < cur_d) continue; for(int i = 0; i < (int)g[v].size(); ++i) { auto& e = g[v][i]; if(e.cap > 0 && dist[e.to] > dist[v] + e.cost + h[v] - h[e.to] + eps) { dist[e.to] = dist[v] + e.cost + h[v] - h[e.to]; prevv[e.to] = v; preve[e.to] = i; que.emplace(dist[e.to], e.to); } } } if(dist[t] == inf) return -1; for(int v = 0; v < (int)g.size(); ++v) { h[v] += dist[v]; } auto d = f; for(int v = t; v != s; v = prevv[v]) { d = min(d, g[prevv[v]][preve[v]].cap); } f -= d; res += d * h[t]; for(int v = t; v != s; v = prevv[v]) { auto& e = g[prevv[v]][preve[v]]; e.cap -= d; g[v][e.rev].cap += d; } } return res; } int main() { int n; cin >> n; int ix, iy; cin >> ix >> iy; vector<int> sx(n), sy(n), gx(n), gy(n); for(int i = 0; i < n; ++i) { cin >> sx[i] >> sy[i] >> gx[i] >> gy[i]; } double ans = 0; for(int i = 0; i < n; ++i) { ans += hypot(gx[i] - sx[i], gy[i] - sy[i]); } graph g(n * 2 + 2); const int src = n * 2, sink = n * 2 + 1; for(int i = 0; i < n; ++i) { for(int j = 0; j < n; ++j) { add_edge(g, i, j + n, 1, hypot(sx[j] - gx[i], sy[j] - gy[i])); } add_edge(g, src, i, 1, 0); add_edge(g, i + n, sink, 1, 0); } ans += min_cost_flow(g, src, sink, n); cout << fixed << setprecision(10) << ans << endl; }
ICPC 模擬国内予選 2019
はじめに
2019/Practice/模擬国内予選 - ACM-ICPC Japanese Alumni Group
の参加記です。
kazuma, nakano と SleepingDragon として出ました。
コンテスト内容
- A を読む、やるだけなので投げる、WAになる
- 2回目のケースで1回目のテストケース結果提出してた（は？）
- C を読む、こういうのは僕は嫌いなので無視
- D を読む、rotate の位置全探索して編集距離やるだけだと思った
- 若干バグったけどサンプルが合うので投げる、WA
- nakano に前消すときは rotate のコスト減るの忘れてない？と指摘される（忘れてた
- 直しても WA 、絶望
- 15分～20分程度悩んだ結果、サンプルケースの結果提出してたことに気づく（は？？？）
- ICPC は初めてか？肩の力抜けよ
- よく知らないがこの間に B, C, E が通っていた
- G は聞いたらやるだけだったので書く
- バグった
- 終了1分後にサンプルが合った（テストケース見たけど通ってた、悲しいね
5完で全体 13 位、これで6完できないのはダメだね。
本番じゃなくてよかった（本番だったら普通に予選敗退しそう
反省・感想
- 謎が発生したらすぐに相談しような
- サンプルの結果を投げるのをやめろ（提出2回もミスってるの初心者すぎる
- Heno_World おめでとうございます、想像以上に強かった
AOJ 2691 - Cost Performance Flow
解法
流量1ずつ流していったときの最小費用流のコストのグラフは、下に凸な折れ線グラフになる。
求める値は、このグラフと点 \((M, 0)\) との距離であり、グラフの形からこの距離も下に凸な関数となる。
なので、折れ線の各線分に対して \((M, 0)\) との距離関数を求め、極値を計算する。
今流量を \(f_i\) だけ流しており、その時の最小コストが \(S_i\) であるとする。
この状態から、追加で 1 流すとしたときのコストを \(c_i\) とする。
\(f_i\) の代わりに \(f_i + \Delta ~~(0 \leq \Delta \leq 1)\) 流すとすると、その時のコスパは $$(M - f_i - \Delta)^2 + (S_i + c_i\Delta)^2$$ である。
これを微分して極値を得る \(\Delta\) を求めると $$\Delta = \frac{M - f_i - S_i c_i}{c_i^2 + 1}$$ となる。
これが 0 以上 1 以下なら、候補になる。
これらの候補と折れ線の端点の中で一番小さいコストが求める値である。
ところで、適当に実装すると（有理数部分で）オーバーフローするため、頑張ってオーバーフローしないようにするか、int128 を使ってごまかそう。
ソースコード
#include <bits/stdc++.h> using namespace std; using ll = __int128; // めっちゃ長いのでライブラリ部分は省略 int main() { int n, m; cin >> n >> m; capacity_weighted_graph<int, int> g(n); int s, t; cin >> s >> t; s--, t--; for(int i = 0; i < m; ++i) { int a, b, u, c; cin >> a >> b >> u >> c; add_edge(g, a - 1, b - 1, u, c); } const ll M = [&] { auto tmp = g; return max_flow(tmp, s, t); }(); rational<ll> ans{M * M, 1}; auto calc_cost = [&] (ll fi, ll ci, ll ci_sum, rational<ll> delta) { const auto a = M - (fi + delta), b = ci_sum + ci * delta; return a * a + b * b; }; ll ci_sum = 0; for(int i = 0; i < M; ++i) { const ll ci = min_cost_flow(g, s, t, 1); if(0 <= M - i - ci * ci_sum && M - i - ci * ci_sum <= ci * ci + 1) { auto delta = rational<ll>{M - i - ci * ci_sum, ci * ci + 1}; delta.reduce(); ans = min(ans, calc_cost(i, ci, ci_sum, delta)); } ci_sum += ci; ans = min(ans, calc_cost(i + 1, 0, ci_sum, rational<ll>{0, 1})); } cout << (long long)ans.numerator() << "/" << (long long)ans.denominator() << endl; }
感想
オーバーフローが罠すぎる（なんか雑に gcd 取るだけだったらオーバーフローした）。
AOJ 2202 - Canal: Water Going Up and Down
解法
全長が小さいので DP する。
reach[i][j] := 船 i が位置 j に到達する最短時刻
これはあくまでも「到達」時刻であり、その地点を出発できる時刻ではない（門の上下などで停泊する必要があるため）。
あとはこのテーブルの更新と、出発できる時刻を同時に求めていけばOK。
計算量は O(MK)
ソースコード
#include <bits/stdc++.h> using namespace std; constexpr double inf = 1e20; int main() { int n, m, k; while(cin >> n >> m >> k, n) { const int max_len = n + k + 10; vector<double> t1(n), t2(n), V(m); vector<int> idx(max_len, -1); vector<double> gate_time(n); for(int i = 0; i < n; ++i) { int X, L, F, D, UD; cin >> X >> L >> F >> D >> UD; t1[i] = (double)L / D; // east -> west t2[i] = (double)L / F; // west -> east if(UD == 1) { swap(t1[i], t2[i]); gate_time[i] = t1[i]; } idx[X] = i; } for(auto& v : V) cin >> v; vector<vector<double>> reach(m, vector<double>(max_len)); for(int i = 0; i < m; ++i) { double cur = 0; for(int j = 0; j < max_len; ++j) { if(j == 0) { cur = i / V[i]; } else { cur += 1 / V[i]; } if(i != 0 && j + 1 < max_len) { cur = max(cur, reach[i - 1][j + 1]); } reach[i][j] = cur; if(idx[j] != -1) { const int id = idx[j]; const double enter = max(cur, gate_time[id]); cur = enter + t2[id]; gate_time[id] = enter + t1[id] + t2[id]; } } } cout << fixed << setprecision(10) << reach[m - 1][k] << endl; } }
感想
到達時刻と出発可能時刻で式を書き間違えていて破滅した。こういうときに限ってサンプルが全部通るっていうね。
AOJ 2385 - Shelter
解法
こんなの積分するしかないでしょという気分になるので立式する。
とりあえずボロノイ図を考えると、凸多角形領域 \(D\) とその内部の1点に対する問題に帰着できる。
内部の点を \((a, b)\) とする。
この上で式を立てると $$\int\!\!\!\int_D (x - a)^2 + (y - b)^2 dxdy$$ となる。
グリーンの定理より $$\oint_{\partial D} -\frac{1}{3}(y - b)^3dx + \frac{1}{3}(x - a)^3dy$$ を求めればよい。
多角形上の線分の端点を \((x_1, y_1), (x_2, y_2)\) とすれば、\(0 \leq t \leq 1\) を用いて $$x = x_1 + (x_2 - x_1)t$$ $$y = y_1 + (y_2 - y_1)t$$ とおけて、求める周回積分は $$-\frac{1}{3}\int_0^1 \{(y_2 - y_1)t + (y_1 - b)\}^3(x_2 - x_1)dt \\ + \frac{1}{3}\int_0^1 \{(x_2 - x_1)t + (x_1 - a)\}^3(y_2 - y_1)dt$$ となる（厳密にはすべての線分に対してこれを求めて総和を取る）。
これを頑張って手計算して答えを得る。最後に全体の面積で割るのを忘れないように。
ボロノイ図を作るのに \(O(n^3)\) かかるが、\(n \leq 100\) なので間に合う。
ソースコード
#include <bits/stdc++.h> using namespace std; using ld = long double; using point = std::complex<ld>; using polygon = std::vector<point>; constexpr ld eps = 1e-10; ld dot(point const& a, point const& b) { return std::real(std::conj(a) * b); } ld cross(point const& a, point const& b) { return std::imag(std::conj(a) * b); } int ccw(point a, point b, point c) { b -= a; c -= a; if(cross(b, c) > eps) return 1; // a -> b -> c : counterclockwise if(cross(b, c) < -eps) return -1; // a -> b -> c : clockwise if(dot(b, c) < 0) return 2; // c -> a -> b : line if(std::norm(b) < std::norm(c)) return -2; // a -> b -> c : line return 0; // a -> c -> b : line } struct line { line() : a(0, 0), b(0, 0) {} line(point a, point b) : a(a), b(b) {} point a, b; }; point is_ll(line s, line t) { point sv = s.b - s.a, tv = t.b - t.a; assert(cross(sv, tv) != 0); return s.a + sv * cross(tv, t.a - s.a) / cross(tv, sv); } // left side polygon convex_cut(polygon const& p, line l) { const int N = p.size(); polygon res; for(int i = 0; i < N; ++i) { auto a = p[i], b = p[(i + 1) % N]; if(ccw(l.a, l.b, a) != -1) { res.push_back(a); } if(ccw(l.a, l.b, a) * ccw(l.a, l.b, b) < 0) { if(cross(a - b, l.a - l.b) == 0) continue; // cut line が辺に覆いかぶさる res.push_back(is_ll(line(a, b), l)); } } return res; } ld area(polygon const& p) { const int N = p.size(); ld res = 0; for(int i = 0; i < N; ++i) { res += cross(p[i], p[(i + 1) % N]); } return res / 2; } line separate_line(point const& p1, point const& p2) { const auto mid = (p1 + p2) * 0.5L; return line{mid, mid + (mid - p1) * point(0, 1)}; } int main() { int m, n; cin >> m >> n; polygon ps(m), shelter(n); for(int i = 0; i < m; ++i) { int x, y; cin >> x >> y; ps[i] = point(x, y); } for(int i = 0; i < n; ++i) { int x, y; cin >> x >> y; shelter[i] = point(x, y); } ld ans = 0; for(int i = 0; i < n; ++i) { auto range = ps; for(int j = 0; j < n; ++j) { if(i == j) continue; auto l = separate_line(shelter[i], shelter[j]); if(cross(shelter[i] - l.a, l.b - l.a) > 0) { swap(l.a, l.b); } range = convex_cut(range, move(l)); } const int sz = range.size(); const auto a = real(shelter[i]), b = imag(shelter[i]); for(int j = 0; j < sz; ++j) { const ld x1 = real(range[j]), y1 = imag(range[j]); const ld x2 = real(range[(j + 1) % sz]), y2 = imag(range[(j + 1) % sz]); const ld t1 = (x2 - x1) * (pow(y2 - y1, 3) / 4 + pow(y2 - y1, 2) * (y1 - b) + 1.5 * (y2 - y1) * pow(y1 - b, 2) + pow(y1 - b, 3)); const ld t2 = (y2 - y1) * (pow(x2 - x1, 3) / 4 + pow(x2 - x1, 2) * (x1 - a) + 1.5 * (x2 - x1) * pow(x1 - a, 2) + pow(x1 - a, 3)); ans += t2 - t1; } } ans /= 3 * area(ps); cout << fixed << setprecision(10) << ans << endl; }
感想
グリーンの定理とか久々に使った（使わなくても直感的に式は出せる気もする、どうなんだろう）。
