Day 094-Q5 | Welzl's Algorithm（最小包含円）

問題

2次元平面上に $N$ 個の点が与えられる。これら全ての点を内部または境界上に含む円のうち、半径が最小のもの（最小包含円）を求め、中心の座標と半径を出力せよ。

入力形式

N
x_1 y_1
:
x_N y_N

制約

$1 \le N \le 10^5$

$-10^4 \le x_i,y_i \le 10^4$（実数）

誤差 $10^{-6}$ 以内で正解（特別ジャッジ）

入出力例

入力例1

出力例1

2.0000000000 1.5000000000 2.5000000000

直角三角形の最小包含円は斜辺を直径とする円（タレスの定理）。斜辺$(4,0)$-$(0,3)$、中点$(2,1.5)$、半径$2.5$。

入力例2

2
0 0
2 0

出力例2

1.0000000000 0.0000000000 1.0000000000

概念図

ヒント（段階的開示）

ヒント1: 方向性

全3点組の外接円から最大を取る$O(N^4)$全探索では間に合わない。最小包含円は高々3点（境界上の点）だけで一意に決まるという性質を利用し、点を1点ずつ追加しながら円を更新する方法を考えよ。

ヒント2: アプローチ

点をシャッフルしてから処理する。現在の円に含まれない点$p_i$が見つかったら「最小包含円の境界に$p_i$が必ず含まれる」性質を使い、$p_i$を固定して残りの点で円を作り直す（Welzl's algorithmの反復版）。境界2点なら中点を中心とする円、3点なら外接円。乱択で期待$O(N)$を達成。

ヒント3: 誘導（コード骨格）

def min_enclosing_circle(points):
    random.shuffle(points)
    c, r = points[0], 0.0
    for i in range(1, len(points)):
        if dist(points[i], c) > r + EPS:
            c, r = points[i], 0.0
            for j in range(i):
                if dist(points[j], c) > r + EPS:
                    c = midpoint(points[i], points[j])
                    r = dist(points[i], c)
                    for k in range(j):
                        if dist(points[k], c) > r + EPS:
                            c = circumcenter(points[i], points[j], points[k])
                            r = dist(points[i], c)
    return c, r

3重ループに見えるが内側ループの実行回数は乱択により期待定数回に抑えられ、全体で期待$O(N)$になる。

模範解答 (Python)

import sys
import math
import random


def dist(a, b):
    return math.hypot(a[0] - b[0], a[1] - b[1])


def circumcenter(a, b, c):
    ax, ay = a
    bx, by = b
    cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-10:
        pts = [a, b, c]
        best_pair = max(
            ((pts[i], pts[j]) for i in range(3) for j in range(i + 1, 3)),
            key=lambda pq: dist(pq[0], pq[1]),
        )
        return ((best_pair[0][0] + best_pair[1][0]) / 2, (best_pair[0][1] + best_pair[1][1]) / 2)
    ux = ((ax**2 + ay**2) * (by - cy) + (bx**2 + by**2) * (cy - ay) + (cx**2 + cy**2) * (ay - by)) / d
    uy = ((ax**2 + ay**2) * (cx - bx) + (bx**2 + by**2) * (ax - cx) + (cx**2 + cy**2) * (bx - ax)) / d
    return (ux, uy)


def min_enclosing_circle(points):
    pts = points[:]
    random.shuffle(pts)
    n = len(pts)
    EPS = 1e-7

    if n == 1:
        return pts[0], 0.0

    c = pts[0]
    r = 0.0
    for i in range(1, n):
        if dist(pts[i], c) > r + EPS:
            c, r = pts[i], 0.0
            for j in range(i):
                if dist(pts[j], c) > r + EPS:
                    c = ((pts[i][0] + pts[j][0]) / 2, (pts[i][1] + pts[j][1]) / 2)
                    r = dist(pts[i], c)
                    for k in range(j):
                        if dist(pts[k], c) > r + EPS:
                            c = circumcenter(pts[i], pts[j], pts[k])
                            r = dist(pts[i], c)
    return c, r


def main():
    data = sys.stdin.buffer.read().split()
    n = int(data[0])
    pts = []
    idx = 1
    for _ in range(n):
        x = float(data[idx]); y = float(data[idx + 1])
        idx += 2
        pts.append((x, y))

    c, r = min_enclosing_circle(pts)
    sys.stdout.write(f"{c[0]:.10f} {c[1]:.10f} {r:.10f}\n")


main()

計算量: 乱択インクリメンタル法により期待 $O(N)$（事前シャッフルが本質）。

Step-by-Step 解説

1シャッフルと初期円
点をランダム順に並び替え、最初の1点を中心・半径0の円として開始する。

2外側ループ
現在の円に含まれない点$p_i$が現れたら「最小包含円の境界には必ず$p_i$が含まれる」性質(Welzlの補題)を使い円を作り直す。

3中間ループ（2点円）
それでも入らない点$p_j$が見つかったら境界は$p_i,p_j$の両方を含むことが確定し、線分を直径とする円を暫定解にする。

4内側ループ（3点円=外接円）
さらに入らない点$p_k$が見つかれば境界は3点に確定し、外接円を計算する。

5縮退ケースの処理
3点がほぼ一直線上のとき外接円公式は不安定になるため、最も離れた2点を直径とする円にフォールバックする（数学的に正しい）。

よくあるミス

ミス	原因	正しい書き方
点をシャッフルせず処理し最悪$O(N^3)$になる	乱択の意義を理解していない	必ず`random.shuffle`してから処理
3点がほぼ一直線上のときゼロ除算エラー	縮退ケースの考慮漏れ	`abs(d)<EPS`でフォールバック処理
境界上の点を「含まれない」と誤判定	誤差処理不足	`dist > r+EPS`のように許容誤差を入れる
出力精度が低く誤差超過	桁数指定を怠る	`.10f`など十分な桁数で出力

次のステップ

発展: 3次元の最小包含球への拡張（同様の乱択法がそのまま使える）
発展: 動的に点が追加されるオンライン最小包含円（Core-set近似）
次回予告: Master Levelローテーション継続（次回セッションで新テーマ）

自己評価

完全理解おおむね理解要復習

Day 094-Q5 — Welzl's Algorithm（最小包含円）

問題