AOJ 2373 HullMarathon

問題文

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2373

解法

自分は数式で捉えて、ラグランジュの未定乗数法で解いた。
ラグランジュの未定乗数法については以下のサイトがわかりやすい（ただし厳密な証明はない）。
mathtrain.jp

与えられた問題を定式化する。
今、r[0], r[1], ..., r[n - 1] をこの順に反時計回りに使うとして、それぞれの間の成す角を θ[0], ..., θ[n - 1] とする。
すると、以下の問題になる。

maximize sum(r[i] * r[(i + 1) % n] * sin(θ[i]) / 2) (for i = 0, ..., n - 1)
s.t. sum(θ[i]) = 2π (for i = 0, ..., n - 1)

ラグランジュの未定乗数法の通りに方程式を立てると、
r[0] * r[1] / (r[i] * r[(i + 1) % n]) * cos(θ[0]) = cos(θ[i])　　(forall i = 0, ..., n - 1)
を満たすようなθを求める問題になる。これは二分探索でできる。

凸包パートが必要そうに見えるが、全探索するとサボれる。
r のうちどれを使うか最初に決めて、それらだけで next_permutation を行うと、単に三角形の面積の和を求めるだけでよくなる。
ただし、next_permutation をするときに先頭を固定してしわないと定数倍で落ちるので注意（最初を固定して定数倍落とすのが嬉しい場面はたまにあるので覚えて損はない）。

ソースコード

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

using ld = long double;

const ld pi = acos(-1.0);
constexpr ld eps = 1e-9;

int main() {
    int n; cin >> n;
    vector<ld> rr(n);
    for(int i = 0; i < n; ++i) {
        cin >> rr[i];
    }
    sort(begin(rr), end(rr));

    ld ans = 0;
    for(int S = 1; S < (1 << n); ++S) {
        vector<ld> r;
        for(int i = 0; i < n; ++i) {
            if(S & (1 << i)) {
                r.push_back(rr[i]);
            }
        }
        const int m = r.size();
        if(m < 3) continue;
        do {
            ld tans = 0, ang_sum = 0;
            auto check = [&] (ld theta) {
                ang_sum = theta;
                tans = 0.5 * r[0] * r[1] * sin(theta);
                for(int i = 1; i < m; ++i) {
                    const auto ang = acos(max(-1.0L, min(1.0L, r[0] * r[1] / r[i] / r[(i + 1) % m] * cos(theta))));
                    tans += 0.5 * r[i] * r[(i + 1) % m] * sin(ang);
                    ang_sum += ang;
                }
                return ang_sum < 2 * pi;
            };
            ld lb = 0, ub = pi;
            for(int lp = 0; lp < 100; ++lp) {
                const auto mid = (lb + ub) / 2;
                (check(mid) ? lb : ub) = mid;
            }
            if(abs(ang_sum - 2 * pi) < eps) {
                ans = max(ans, tans);
            }
        } while(next_permutation(begin(r) + 1, end(r)));
    }
    cout << fixed << setprecision(10) << ans << endl;
}