点の集合を包含する球
石畑

今回は，幾何の問題を取り上げる．3 次元空間内の位置が与えられた複数の点すべてを包含する最小の球を求める問題である．2001 年 11 月の函館大会で問題H として出題された．

■幾何の問題に関する一般論

清（明治大学理工学部）
ishihata@cs.meiji.ac.jp

うな式の計算を一生懸命やってからでないと，プログラムを書き始めることはできない．
実は，3 点の座標を (x1, y1), (x2, y2), (x3, y3) とすると，三角形の面積は次の式によって簡単に求まる．
S = �  (y3y1)(x2x1)(y2y1)(x3x1) �

これなどは，知らなければとうてい導けない式だろプログラミングコンテストの問題は，数多くの分野う．座標に関する代数計算を必死でやれば得られるこから出題される．これは，いろいろな得意分野を持つとは確かだが，限られた時間のコンテストでそんなこ学生たちを公平に競わせるという観点から，意図的にとをしている暇はないはずである．

行われていることであろう．数値計算，組合せ探索，上の式の絶対値記号の中の部分，(y3y1)(x2x1)

グラフ，テキスト処理，… などさまざまなテーマの問(y2y1)(x3x1) には，三角形の面積計算以外の場面で題が出題される．

も利用価値がある．特に，点の位置関係を判定する場これら数多く存在する分野の中でも，幾何は特に重合に利用できる点が重要である．この式の値は，3 点要な分野の 1 つである．どのコンテストでも 1 問や 2が反時計回り順に与えられているなら正，時計回りな問は必ず出題されている．実際に数えたことはないが，ら負になる．同じことを別の言い方で言うと，点 1 におそらく問題全体の 20% 近くは幾何の問題ではないか

立って点 2 の方向を見ているとき，点 3 が左側にあれと思われる．

ば正，右側にあれば負になる．たとえば，点の集合の一口に幾何の問題といっても，問題の傾向はもちろんさまざまである．図形の種類だけでも，多角形を扱凸包 (convex hull) を求める問題などで利用すると便利な式である．

うもの，円を扱うもの，点の集合を扱うもの，… とい結論として，幾何の問題を能率的に解くためには，ろいろである．2 次元もあれば，3 次元もある．面積をこの式を始めとする幾何計算の常識をいくつか身に付求める問題，特定の性質を満たす経路を求める問題，けておくべきである．コンテストでは参考書の持込み何かの値の最大化または最小化を行う問題，… と計算が許されているから，公式集に頼るという考え方も成の目的も数え切れない．

立するが，その場合でも必要とする式をすぐに引き出

幾何の問題を解くときは，プログラミングの技法をせるよう練習しておく必要がある．

考える前に，対象となる図形に関する数学的な考察を幾何問題のもう 1 つの特徴として，さまざまな解法行わなければならないことが多い．たとえば，多角形が存在し得るということを挙げておきたい．いろいろの面積を求める問題を考えてみよう．いろいろな解きな発想に基づいて，複数の方向からの攻略が可能な問方があるが，多角形を三角形に分割した上で，三角形題が多い．初等幾何の問題が補助線一発でスマートにの面積の和として求める方法を思いつくことは難しく解けるのと同じで，難しく見えた問題が，考え方を少ない．しかし，3 点の x-y 座標が与えられているとしし変えるだけで，いとも簡単に解けることもある．プて，三角形の面積を求めるにはどうしたらよいだろう

ログラミング問題としての面白さという点では一番でか．3 点の座標に対して面積の公式（底辺×高さ÷ 2）あろう．

を直接適用できるはずもない．高校の数学でやったよ

今回取り上げる問題も，常識的な解法のほかに，複IPSJ Magazine Vol.43 No.9 Sep. 2002

1009

数の攻略法が存在する．そのうちの 1 つは，いかにも切れ味の鋭い見事な解法である．

■問題：点の集合を包含する球■代数方程式を解いて球の中心を求める方法
幾何の問題の例にもれず，この問題も数学的な考察なしではどこから手を付けたらよいかすら分からない．

3 次元空間の中に n 個の点があって，これらの点の位置が与えられている．これらの点すべてを含む最小図形的な考察を少し行えば，次のような結論に到達することは難しくない．の球を求めるプログラムを書くことが問題である．球が点を含むという場合，点は球面の上に乗っていてもよいし，球の内部に位置するのでもよい．
点の個数 n は，4 以上 30 以下である．点の位置は，y，z のいずれもごく普通の x-y-z 座標で与えられる．x ，n 個の点すべてを包含する最小の球は次のいずれかである．
(1) n 個のうちの 4 個の点で作られる四面体に外接する球

0 以上 100 以下の浮動小数点値である．また，点と点(2) n 個のうちの 3 個の点で作られる三角形の外接円をは少なくとも 0.01 以上離れていると規定されている．

大円とする球（三角形の外接円をその直径を軸とし点の位置に関するこれらの制約は，プログラムを書くて回転して得られる球）上で必ずしも必要ではないが，極端に数値的な性質の(3) n 個のうちの 2 個の点を直径の両端とする球

悪いデータはないと分かっているだけで安心感が生まれる．逐次近似を行う場合の初期値の設定や，誤差の見積りにも好都合である．

この結論の正しさを確認する推論は次のようなものである．説明や図解を簡単にするために，ここではプログラムは，求めた球の半径を出力しなければ2 次元平面上の最小包含円を求める問題を対象にとっならない．半径が求まっていれば，球の中心の位置もて説明する．問題の本質は変わらないので，十分に類分かっているはずだが，こちらの出力は要求されてい推できるはずである．

ない．

2 次元の場合は，ケース (1) がなくなる．ケース (2)半径の出力は，小数点以下 5 桁までと規定されていと (3) のいずれかであることを示せばよい．

る．また，真の値からの誤差が 0.00001 以下でなければならないとも要求されている．これは，実数値を答(a) n 個の点すべてを包含する円の円周上に乗っているとする問題でよく使われる誤差の規定法である．小数点が 0 個だったとする（図 -1(a)）．この場合，円の点以下 6 桁目の四捨五入を間違える程度までは許すと半径を小さくしていくことによって，より小さい円いう意味を含んでいる．10 進で 7 桁程度の精度があれが得られる．したがって，この円が最小であるはずばよいのだから，普通にプログラムを書けば，誤差のはない．円の縮小は，n 個の点のいずれかにぶつかせいで間違う恐れはほとんどない．特別に誤差が大きるまで続けることができる．

くなりやすい問題でもない限り，あまり誤差のことを気にする必要はないだろう．
以下では，入力した点の個数と位置が次のような変数に記録されていると仮定する．

(b) 円周上の点が 1 個だったとする（図 -1(b)）．この場

合，その点を円周上に乗せたまま，円の中心をその点に近づけることによって，より小さい円が得られる．この場合も，この円が最小であるはずはない．
(c) 円周上の点が 2 個だったとする（図 -1(c)）．2 点が直径の両端であれば，それが最小である．そうでな

typedef struct { double x, y, z; } pos;
int n;
pos point[30];ければ，この 2 点を乗せたまま，円の中心を 2 点のpos は，点の 3 次元座標値を記録するための型であれる．円の中心はもともと 2 点の垂直 2 等分線の上る．変数 n に点の個数，配列 point の point[0] 〜point[n-1] に各点の位置を記録する．

中点に近づけることによって，より小さい円が得らにあったことに注意．
(d) 円周上の点が 3 個以上だったとする．この場合，円はそれらの点のうちの 3 個で作られる三角形の外接円になっている．

1010 43巻9号

情報処理

2002年9月

（a）

（b）

（c）

図 -1

結局，(a) と (b) の場合は最小ではなく，(c) が最小

これよりオーダが低い．

なら (3) のケースに，(d) が最小なら (2) のケースに該O(n5) の計算量では遅すぎる感じがするが，この問当することが分かる．以上で，(2) と (3) のケースを尽題の場合は，この程度でも十分間に合う．n の最大がくせば最小包含円が求まることを説明できた．

30 だからである．このアルゴリズムでも正解として元の問題のプログラムの作り方に戻ろう．上の推論

許すよう n の上限を小さく定めてあるのだと想像できさえできれば，次のような単純明快なプログラムに到る．一般に，幾何の問題の場合，アルゴリズムの計算達することは難しくない．

量を小さくすることはあまり重要でなく，解の求め方を図形的な考察によって構成することができればよし• あらゆる点の組合せに対して (1) ，(2) ，(3) で決まとしているものが多いようである．

る球を片端から求める．それらの球の中で，他のすこの方式のプログラムはかなり長いものになるので，べての点を含んでいて，かつ最小のものを見つける．

実際のプログラムを示すことは省略する．注意すべき点を 1 つだけ述べておこう．

実際には，(1) ，(2) ，(3) の計算は包含球の中心位4 点がたまたま同一平面上にあった場合，この 4 点置の候補を与えるだけと割り切ってしまう方が簡明でで決まる四面体の外接球の中心は求まらない．不定まあろう．何はともあれ，それらの位置を中心とする包たは不能の方程式を解くかたちになる．具体的には，含球を作ってみて，半径最小のものを見つければ，そ計算の途中で除数 0 の割り算が生じるはずである．同れが求める最小包含球である．中心位置が与えられて一直線上の 3 点があった場合も同様である．

いれば，包含球を求める操作は簡単である．n 個の点の中で最も遠い点までの距離を半径とすればよい．
(1) で決まる球の中心の座標は，4 点からの距離が等しいという代数方程式を立てて，これを解けば求まる．

プログラムには，これらの場合を検出する仕掛けを入れておかなければならない．これらの場合が本来の問題の解であることはないので，検出したら，その場合を単に無視するだけでよい．

未知数が x ，y ，z の 3 つ，方程式も 3 つの 3 元連立一次最近の計算機は，浮動小数点型の値として無限大や方程式である．(2) の方が複雑（3 点と同一平面上にあNaN を扱える．これを利用すれば，同一平面上の 4 点るという式が必要）だが，こちらも 3 元連立一次方程などを特別扱いせずに済む可能性がある．不定または式を解くことによって求められる．

不能の方程式の場合は，球の中心として何でも（無限(1) の球は，全部で nC4 個ある．これらの球の中心遠でもデタラメな位置でも）かまわないから，適当なそれぞれについて，n 個の点までの距離を計算する必点の座標を返してくれればよい．その点を中心とする要がある．球の中心を求める作業の回数が O(n ) あり，包含球は最小のものより大きいはずだから，単に無視距離の計算にそれぞれ O(n) の計算量がかかるので，されるだけである．筆者は，このような考え方に基づ全体の計算量は O(n ) になる．(2) や (3) の計算量は，いてプログラムを書いたのだが，これはうまくいかな

45IPSJ Magazine Vol.43 No.9 Sep. 2002

1011

（a）（b）
図 -2

かった．数値計算ライブラリーの一部に NaN に十分に対応していないものがあるらしく，包含球の半径が0.0 になってしまった．やはり特別扱いは必要である．
複数個の点を包含する円や球を求める問題は，計算• プログラムでは，球の中心を含む可能性のある複数個の立方体を常に把握する．
• 最初に記憶する立方体は x ，y ，z 方向とも最小 0 ，最大 100 のもの 1 個だけである．

幾何学の分野で十分に研究し尽くされた有名問題の• 各立方体の中心から n 個の点までの距離の最大値を1 つである．ここに述べたのと同様の解法をアルゴリ求める．その点を中心として n 個の点を包含する球ズム的に改良すると O(n) にまで計算量を小さくできを作ったときの半径である．これを r とする．立方ることが知られている．文献として次の 2 つを挙げて体の 1 辺の長さを t とすると，立方体の中に中心をおく．

置いてできる包含球の半径が r √3t2 を下回ることはない．√3t2 は，立方体の中心から 8 つの頂点

1) ドバーグ他（浅野訳）: コンピュータジオメトリ，近代科学社 (2000).

までの距離である．
• 計算の進行に伴って，多数の立方体についての r を2) Welzl, E.: Smallest Enclosing Disks ( balls and次々に求めていくが，それらの r の値の最小値を記ellipsoids ) , in H.Maurer ( ed. ) New Results and憶する．これを rmin と呼ぶことにしよう．ある立方New Trends in Computer Science, Lecture Notes in体について，r √3t2 が rmin を超えた場合は，そのComputer Science, Vol.555, pp.23-33 (1988).

立方体を捨てる．最小包含球の中心がその立方体の中に存在する可能性がないからである．

■空間を分割する方法• 捨てられずに残った立方体は，x ，y ，z の 3 つの方向に2 等分して，合計 8 つの小立方体に分ける．そして，幾何の問題は，正確に解こうとすると大変だが，近似的な解法なら簡単ということが多い．この問題の場合も，これが当てはまる．

同様の操作を繰り返す．
• 立方体の大きさが誤差の限界を下回るまで小さくなったら，アルゴリズムを停止する．

幾何問題の近似解法としてよく使われる戦略に，空間を立方体の集まりとして把握し，それを順次分割し

図 -2 に，この 方 式を 2 次 元 で 図 解し たものを 示ていくものがある．立体図形の表現に使うデータ構造

す．図で×印を付けた点が最小包含円の中心である．に octree（8 分木）と呼ばれるものがあるが，これと同

図 -2(a) では，平面全体を 16 の正方形に分けてあるが，様の方法である．

このうち候補として生き残っているのは 5 つである（斜この問題の場合は，次のような考え方でこの方法を適用できる．

1012 43巻9号

線を付けたもの）．これらの正方形をそれぞれ 4 分割（2 次元だから 8 分割でなく 4 分割）したものが図 -2(b)情報処理

2002年9月

である．この時点では，6 つの正方形が候補として生

continue;if (k >= Max_Candidate)goto end_expand;new_candidate[k].p = p;
new_candidate[k].r = r;k++;

き残っている．生き残りの正方形の分割を続けていけば，最小包含円の中心が存在し得る領域の面積をどんどん小さくすることができる．
この方式を実際のプログラムとして書いたものを示}}
end_expand:count = k;for (i = 0; i < count; i++)candidate[i] = new_candidate[i];sort_candidate(count);}
return (min_sofar);す．初めに，使用する型や大域変数の宣言である．
typedef struct { pos p; double r; } can;can candidate[Max_Candidate];
can new_candidate[Max_Candidate];
int

sign[8][3] = {{ 1, 1, 1 },{ 1, 1, -1 },{ 1, -1, 1 },{ 1, -1, -1 },{ -1, 1, 1 },{ -1, 1, -1 },{ -1, -1, 1 },{ -1, -1, -1 } };}while 文は，立方体の大きさを半分にする操作の繰返しを表す．for 文が二重になっているが，このうちの外側の方で候補立方体を順次調べていき，内側で 1 つの立方体を 8 つに分割する操作を行っている．立方体

can は，候補立方体を記録するための型である．メを分割したら，すぐにそれを候補として残すか捨てる

ンバ p は立方体の中心位置の座標，r はその点を中心かを判断している．compute_radius(p) は，点 p をとしたときの包含球の半径を表す．

中心とした包含球を作ったときの半径を求める関数で

配列 candidate と new_candidate に候補立方ある．

体を記憶する．筆者のプログラムでは，配列の長さ最初の方法に比べて，少なくともプログラミングはMax_Candidate を 100000 とした．配列 sign は，ずっと楽である．こんな簡単な方法でも，大抵の場合8 等分を簡単に行うための補助配列である．

は答が正しく求まる．
このプログラムの欠点は，候補立方体の数が多くな

double solve_by_partition(void){
int count, k, i, d;
pos p;
double size, r, root3, min_sofar;candidate[0].p.x = 50.0;
candidate[0].p.y = 50.0;
candidate[0].p.z = 50.0;
count = 1;
size = 100.0;
min_sofar = 1.0e+50;
root3 = sqrt(3.0);
while (size > 0.000002) {size = size/2.0;
k = 0;for (i = 0; i < count; i++) {for (d = 0; d < 8; d++) {p = candidate[i].p;
p.x += sign[d][0]*size/2.0;
p.y += sign[d][1]*size/2.0;
p.z += sign[d][2]*size/2.0;
r = compute_radius(p);if (r < min_sofar)
min_sofar = r;
if (r > min_sofar+size/2.0*root3)って，記憶域不足になる可能性があることである．最悪のケースは，2 点を直径の両端とする球が答である
場合に生じる．2 点を結ぶ線分の中点から，この線分と直角の方向に少し移動した点は，2 点までの距離が真の中心とあまり変わらない．このような点を含む立方体は，r √3t2 の打切りではなかなか捨てられずに残る．線分の中点を中心とする円盤上にこのような立方体が多数存在するので，候補立方体の個数はすぐに数万のオーダを超え，億に達することもある．
候補の数がこのように爆発しても，元の問題が解けさえすればよい．これを可能にするために，上のプロ
グラムでは，各ステップの最後の時点で，候補として残っている立方体をよさそうな順にソートしている．
つまり，立方体の中心を中心としたときの包含球の半径（candidate[i].r の値）が小さいものほど配列の前の方に位置するように並べ替えをしている．これを担当するのが関数 sort_candidate( count ) である．
分割の結果として，立方体の数が多くなりすぎた場
IPSJ Magazine Vol.43 No.9 Sep. 20021013

合は，その時点で分割操作を打ち切って，配列の後の方にある立方体を分割する操作は行わない．真の解は，

��おそらく配列の前の方にあるはずなので，打切りの前に処理が終わっているだろうと期待してのことである．
このほか，ここには示さないが，筆者のプログラム

����

��
��
��には候補数の爆発そのものを避けるための工夫も組み込んである．2 つの点を直径の両端とする球が包含球になっているかどうかを調べる．このチェック自体は簡単である．包含球になっていれば，それが答なので，立方体分割の作業を始めずに直ちに結果を返す．

����

図 -3

上に述べたとおり，候補の数が最も多くなるのは，直径の両端の 2 点だけが最小包含球の決定に貢献する場合である．この場合だけ特別扱いすれば，候補数爆うまくいかない．これは，現在位置の周囲で，より包発の可能性を大きく減らすことができる．

含球の半径の小さい点を探して，そちらに移動するともちろん，この工夫だけで十分とは言えない．候補

いう方法である．半径が小さくなる方向を見つけるこ

数の爆発が起こるデータはいくらでも作れる．いったとが大変なので，この方法はほとんど成立しない．場ん爆発が起こったら，何しろ計算の途中でデータを捨合によっては，きわめて小さな角度の方向に移動しなててしまっているのだから，必ず正しい答が求まるといと半径を小さくできないことがある．

は保証できない．しかし，コンテストで使われた審判データ（爆発を起こす例も含まれている）で試してみた

次に示すプログラムでは，現在位置の移動の規則として次のようなものを採用している．

ところ，候補の立方体の数が 400 あれば，すべての例について正解が得られた．候補の数を数万の程度にま• 各ステップでは，n 個の点の中で現在位置から最もで増やせば，間違える恐れはほとんどないと考えてよ

遠い点の方向に向かって適当な距離だけ移動する．

さそうである．

• 移動距離は，ステップを繰り返すに従って徐々に小計算スピードは，あまり速いとは言えない．候補数さくしていく．

を 100000 として，審判データで走らせると，コンテストの時間制限（3 分）に引っ掛かるかどうか微妙なとだまされたように感じるかもしれないが，このようころである．候補数をこの半分くらいに減らせば安全な単純なルールによる移動で，確かに最小包含球の中であろう．

心が求まる．図 -3 に 2 次元の例を示すので，これで納得してもらいたい．図で，S1 〜 4 は与えられた n 個の■球の中心に近づく操作を繰り返す方法

点を表す．P0 〜 3 が現在位置の推移である．P0 から各点までの距離を計算すると S1 が最遠方なので，S1 の今度は，同じ近似解法でも，上のプログラムとはま方向へ向かって動き，P1 に移動する．ここで再び各点ったく違った発想に基づく方法を示そう．任意の位置までの距離を計算すると S3 が遠いので，その方向にあから出発して，小さな移動を繰り返す（現在位置を次々る P2 へ移動する．以下同様である．1 回当たりの移動と変えていく）ことによって，最小包含球の中心に順距離を徐々に小さくしていけば，求める最小包含球の次近づいていくという方法である．空間内の 1 点を求中心（×印）に近づいていくことが見て取れると思う．

めるための逐次近似法として最も自然な考え方であろう．たとえば Newton 法を使って高次方程式の根を求めるのと同様の方法である．
ここまでは誰でも思いつくのだが，ハードルは高い．
ポイントは，どのような規則に従って移動する先の位置を決めるかということである．最も単純な考え方としては，山登り法（hill climbing）があるが，これでは1014 43巻9号

情報処理

2002年9月

double solve_by_movement(void){
int k, i, t;
pos p;double move, max;
p.x = 50.0;
p.y = 50.0;

p.z = 50.0;
move = 0.5;
while (move > 1.0e-8) {for (t = 0; t < 100; t++) {max = 0.0;
for (i = 0; i < n; i++)if (distance(p, point[i]) > max) {max = distance(p, point[i]);k = i;}
p.x += (point[k].x-p.x)*move;
p.y += (point[k].y-p.y)*move;
p.z += (point[k].z-p.z)*move;}
move = move/2.0;}
return (max);}

変数 move は，最も遠い点までの距離に対する移動タに依存して，悪い性質を示すということもなさそうである．この問題の解法として決定版と言ってよいだ
ろう．
初めに述べたとおり，幾何問題には複数の解法が成立することが多い．包含球問題は，その好例である．
幾何学的な考察に基づく厳密解法，空間を立方体に分けることによる近似解法，求める点に徐々に近づいていく近似解法と，バラエティに富む解法が可能であった．プログラミングの技法としてもそれぞれ面白い点を含んでいて，挑戦しがいのある問題だと思う．

■付言
Program Promenade の執筆者 4 人のうち，筆者だけがプログラミングコンテストの現役の審判である．

距離の比を表す．初めは，最も遠い点に向かって，そコンテストの審判団は，出題する問題それぞれにつの距離の半分（0.5 倍）だけ近づく．これを 100 回繰りいていろいろなプログラムを準備する．正解として想返した後，比率を 0.25 にして，また 100 回繰り返す．

定するプログラム，アルゴリズムが悪くて遅すぎるプ

以下同様で，移動距離の比を半分半分と小さくしてい

ログラム，場合分けが不十分で正解にならないプログ

きながら，それぞれ 100 回ずつ移動を繰り返している．

ラム，等々である．これらのプログラムに基づいて，0.5 とか 100 回とかは，いい加減に決めた値で，こどんなデータを準備すればプログラムの正しさやスピれが最善だという根拠は何もない．収束の具合に関すードを検査できるか検討する．この作業の過程で，複る考察を加えることによって，もっと少ない繰返し回数の審判から，いろいろな面白いプログラムのアイデ数で解に到達することは可能だと思われる．しかし，ィアが提示されることが多い．

ここに示した方法でも十分に速いし，確実に答が得ら今回示したプログラムは，審判団の他のメンバが提れる．細かな改良を加えることにあまり意味はないだ

案した解法を筆者なりの記法で書き直したものである．

ろう．

本当は，提案者である審判の実名を記して謝辞を述べこの方法で最小包含球の中心を求められることの証るべきなのだが，審判団の内部事情を明かすことは好明（現在位置と最小包含球の中心の間の距離が次第にましくないと思われる．ここでは，これらの事情を記小さくなっていく様子を示すことになるだろう）は結録するだけにとどめたい．

構大変である．しかし，直観的な理解でよければ，難（平成 14 年 8 月 13 日受付）しいことはない．その鍵になるのは，元の問題を一種の最適化問題として定式化することである．
問題で要求しているのは minc maxi distance(pi, c) を求めることである．点 pi と空間上の任意の点 c の距
離を distance(pi, c) として，n 個の点の中の最大のものを選ぶと，これが包含球の半径になる．この最大値（半径）を最小にするような点 c が求める球の中心である．現在位置を最も遠い点の方向に移動すれば，maxi distance(pi, c) が小さくなる．この操作を繰り返すことによって最小値に到達できることも，図形的に考えれば，ほぼ明らかである．
この方法は，プログラムの行数が最も少ないし，確実に答が得られる．実行時間もきわめて短い．デー
IPSJ Magazine Vol.43 No.9 Sep. 2002

1015