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;
}

感想

グリーンの定理とか久々に使った(使わなくても直感的に式は出せる気もする、どうなんだろう)。

AOJ 2569 - Putter

解法

線分の通り方を n! 通り全探索する。
ある線分を通る時に、その線分で点を鏡写しにしていけば、直線がそれらの線分を貫くかどうか、という問題になる。
この判定は角度(厳密にはベクトル)で行う。
始点から、ある線分を通るような範囲をベクトルで表すことにする。
つまり、始点から線分の端点へのベクトル2つを、より反時計回り方向にある方とそうでない方で管理する。
すべての範囲が共通部分を持てばOK。
そのために、反時計回りにある方のベクトルだけを考えて、それらの中でもっとも時計回り方向にあるものを考える(時計回りの方は逆)。
こうして得られた2つのベクトルが、反時計回りのほうがちゃんと反時計回りの方にあるならOK。
ここまでの計算は全部外積だけでできる。

ただし、反時計回りのやつを時計回りに、時計回りのやつを反時計回りの方に寄せていった結果、最終的に寄せた大きさが180度を超えて条件合致してしまうケースがある。
なので、最も(反)時計まわりにあるものを考えるというよりは、ベクトルの範囲を1つずつ見ていって、常に共通部分を持つもの、とする。
convex polygon なのでこれでちゃんと判定可能(今条件合致していて、次に鏡写しにしたときのベクトルの範囲を考えたときに、範囲が大きく変動して、偶然条件合致と判定されないということ)。

ソースコード

#include <bits/stdc++.h>
using namespace std;

using ld = double;
using point = std::complex<ld>;
using polygon = std::vector<point>;

struct line {
    line() : a(0, 0), b(0, 0) {}
    line(point a_, point b_) : a(a_), b(b_) {}
    point a, b;
};

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);
}

point proj(line const& l, point const& p) {
    ld t = dot(p - l.a, l.a - l.b) / std::norm(l.a - l.b);
    return l.a + t * (l.a - l.b);
}

point reflect(line const& l, point const& p) {
    return p + (proj(l, p) - p) * 2.0;
}

point read_point() {
    ld x, y; cin >> x >> y;
    return point{x, y};
}

int main() {
    int n;
    while(cin >> n, n) {
        const auto sp = read_point();
        polygon ps;
        for(int i = 0; i < n; ++i) {
            ps.emplace_back(read_point());
        }

        vector<int> ord(n);
        iota(begin(ord), end(ord), 0);
        int ans = 0;
        do {
            auto ts = ps;
            vector<point> lvs, rvs;
            for(auto i : ord) {
                const auto a = ts[i], b = ts[(i + 1) % n];
                lvs.push_back(a + (b - a) / abs(b - a) * 1e-6 - sp);
                rvs.push_back(b + (a - b) / abs(b - a) * 1e-6 - sp);
                if(cross(rvs.back(), lvs.back()) < 0) swap(lvs.back(), rvs.back());
                for(int j = 0; j < n; ++j) {
                    if(j == i || (i + 1) % n == j) continue;
                    ts[j] = reflect(line{a, b}, ts[j]);
                }
            }
            auto lv = lvs[0], rv = rvs[0];
            bool check = true;
            for(int i = 1; i < n; ++i) {
                if(cross(lv, lvs[i]) <= 0) swap(lv, lvs[i]);
                if(cross(rv, rvs[i]) >= 0) swap(rv, rvs[i]);
                check &= cross(rv, lv) >= 0;
            }
            ans += check;
        } while(next_permutation(begin(ord), end(ord)));

        cout << ans << endl;
    }
}

感想

(反)時計回り側のベクトルを全部統合してから判定してて最初死んでいた。
多分本番だと無限に悩むやつ。そして多分通らないまま終わる。悲しいね。

AOJ 2164 - Revenge of the Round Table

解法

回転して同じものはとりあえず重複して数える。
というのも、もし周期 s の並べ方があったとしたら、回転して同じとみなせるものは s 通りだとわかるからだ。

まず、先頭を A に固定して、ループを考慮しない数え上げをする。
次に、ループを考慮して数え上げるが、これは先頭と末尾の A の個数を決め打ちして、間に B から始まって B で終わるものを考えればOK.
これで先頭 A でループも考慮したものが得られて、先頭を B にしたものも同じなので 2 倍する。

さっきので求まったテーブルを dp[i] とすると、周期が i で長さも i であるような列の総数 dp2[i] は
 dp2[i] = dp[i] - (i のすべての約数 j に対して dp2[j] の総和)
となる。

求める答えは、n のすべての約数について dp2[i] / i の総和となる。
計算量は O(N^2)

ソースコード

#include <bits/stdc++.h>
using namespace std;

using ll = long long;

constexpr int mod = 1000003;

ll inv(ll n) {
    ll a = n, b = mod, u = 1, v = 0;
    while(b) {
        ll t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= mod;
    if(u < 0) u += mod;
    return u;
}

int main() {
    int n, k;
    while(cin >> n >> k, n) {
        if(n == 1) {
            cout << 2 << endl;
            continue;
        }
        ll ans = 0;
        if(k >= n) {
            ans += 2;
            k = n - 1;
        }

        vector<vector<vector<int>>> dp1(n + 1, vector<vector<int>>(k + 1, vector<int>(2)));
        dp1[1][1][0] = 1; // start with A
        for(int i = 1; i < n; ++i) {
            for(int j = 1; j <= k; ++j) {
                for(int pre = 0; pre < 2; ++pre) {
                    if(j + 1 <= k) {
                        (dp1[i + 1][j + 1][pre] += dp1[i][j][pre]) %= mod;
                    }
                    (dp1[i + 1][1][!pre] += dp1[i][j][pre]) %= mod;
                }
            }
        }

        vector<ll> dp2(n + 1), sum(n + 1);
        for(int i = 1; i <= n; ++i) {
            for(int j = 0; j <= k; ++j) {
                (dp2[i] += dp1[i][j][1]) %= mod; // ends with B
                (sum[i] += dp1[i][j][0]) %= mod;
            }
            for(int loop = 2; loop <= min(k, i); ++loop) {
                (dp2[i] += 1LL * sum[i - loop] * (loop - 1)) %= mod;
            }
            (dp2[i] += dp2[i]) %= mod; // consider starting with B
        }

        vector<int> dp3(n + 1);
        for(int i = 1; i <= n; ++i) {
            dp3[i] = dp2[i];
            for(int j = 1; j < i; ++j) {
                if(i % j == 0) {
                    (dp3[i] += mod - dp3[j]) %= mod;
                }
            }
        }
        for(int i = 1; i <= n; ++i) {
            if(n % i != 0) continue;
            (ans += dp3[i] * inv(i)) %= mod;
        }

        cout << ans << endl;
    }
}

感想

n == 1 がコーナーケースになってる実装だから良くないなあ。
実装が悪い。