AOJ 1301 Malfatti Circles

問題文

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

解法

まず，一つの円を決め打ちします．
すると，残りの2円の位置は，二分探索などをして求めることが出来ます．
ここで，残りの2円が離れすぎている場合，最初に決め打ちした円が小さすぎるということになります．逆に，近すぎる場合（交わっている場合），最初に決め打ちした円が大きすぎるということになります．
なので，最初の円を二分探索でちょうどいいサイズにすると，答えが求まります．
二分探索の上限の円の中心位置は，三角形の重心にしておかないとめんどくさいことになりそうな気がします．

ソースコード

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

using ld = long double;

constexpr ld eps = 1e-8;

// ライブラリ部分は長いので略．

point rot(point p, ld ang) {
    ld x = real(p) * cos(ang) - imag(p) * sin(ang);
    ld y = real(p) * sin(ang) + imag(p) * cos(ang);
    return point(x, y);
}

circle calc(point cp, ld r, point const& limit, line l1, line l2) {
    point lb = l1.a, ub = limit;
    ld rad = 0;
    for(int i = 0; i < 60; ++i) {
        point m = 0.5l * (lb + ub);
        rad = dist_lp(l1, m);

        if(abs(cp - m) + eps < r + rad) {
            ub = m;
        } else {
            lb = m;
        }
    }
    return circle(lb, rad);
}

int main() {
    int x1, y1, x2, y2, x3, y3;
    while(cin >> x1 >> y1 >> x2 >> y2 >> x3 >> y3) {
        if(x1 == 0 && y1 == 0 && x2 == 0 && y2 == 0 && x3 == 0 && y3 == 0) {
            break;
        }
        point p1(x1, y1), p2(x2, y2), p3(x3, y3);
        point vec1 = rot(p2 - p1, (arg(p3 - p1) + arg(p2 - p1)) / 2 - arg(p2 - p1));
        point vec2 = rot(p1 - p3, (arg(p1 - p3) + arg(p2 - p3)) / 2 - arg(p1 - p3));
        point g = is_ll(line(p1, p1 + vec1), line(p3, p3 + vec2));

        point lb = p1, ub = g; 
        ld r1, r2, r3;
        for(int i = 0; i < 60; ++i) {
            point m = (ub + lb) * 0.5l;
            r1 = dist_lp(line(p1, p2), m);

            auto c2 = calc(m, r1, g, line(p2, p1), line(p2, p3));
            auto c3 = calc(m, r1, g, line(p3, p1), line(p3, p2));
            r2 = c2.r;
            r3 = c3.r;
            if(abs(c2.p - c3.p) + eps < c2.r + c3.r) {
                lb = m;
            } else {
                ub = m;
            }
        }

        cout << setprecision(10) << fixed;
        cout << r1 << ' ' << r2 << ' ' << r3 << endl;
    }
}

感想

多分かなり遅い実装をしていて，ループ回数をそれぞれ100ずつにするとTLEギリギリ通るぐらいでした（座標のサイズ的に100回も回す必要はないんですけどね…）．
まあ通ればなんでもいいよね．速度的には一方の二分探索をやめて，2次方程式で直接求めるほうが早いでしょう．