0から始める計算幾何学第02回近似的アルゴリズム

厳密な解を出すことをとっとと諦めて、少しだけ誤差を含んだ解を簡単に出してみよう。
今回は、そんな近似的アルゴリズムのお話です。

今回取り上げる問題は、昨年2月に行われたオンラインコンテストで出題された「Airport」です。現在はこちらから見ることができます。

概要

xy 平面上に n 個の点がある。ここに直線を引き、各点と直線の平均距離を求めることを考える。
ただし、全ての点が線から見て同じ側にあるように線を引かなければならない (点と線は重なってもよい)。
平均距離が最も小さくなるように線を引いたとき、その平均距離はいくらか。

この問題を読んで、どこに線を引けばいいかすぐに分かる人はいいのですが、あいにく私はそれほど頭がよくありません。もっと楽してズルしていただきな方法を探すことになります。
例えば、こんな方法はどうでしょうか。

適当な場所に線を引いてみる

平均距離が最小となる直線は、少なくとも1個の点と重なっていることは明らかです。
なので、まず傾き θ の直線を下から近づけていき、どこかの点と重なったところで止めて平均距離を測ります。
この操作を全ての θ について行い、求まった平均距離の中で最小のものを見つければ、それが答えとなります。

ですがもちろん θ は実数ですから、無限個ある全ての θ を調べることなどできるはずがありません。
ここで妥協を差し挟みましょう。0°～360°の間で等間隔に1000個の角度を取り、その中で平均距離の最小を見つけるのです。

そんな乱暴なことをして大丈夫なのかと思うかもしれませんが、とりあえずコードを書いてみましょう。

const int MAX_HOME_NUM = 10000; // 座標の最大個数
const int DIVISION = 1000;      // 分割数
main (){
Point home[MAX_HOME_NUM];
int homenum;
// ここで入力処理
　double avedist=getdistsum(home,homenum,0,2*M_PI)/homenum;
　printf ("%.3lf\n",avedist);　
}
double getdistsum(Point *home,int hnum,double anglefrom,double anglesize){
// (anglefrom)°から(anglefrom+anglesize)°を(DIVISION)個に分割し、
// 平均距離の最小を返します
Point rotate[MAX_HOME_NUM];
double distance[DIVISION];
const double theta=anglesize/DIVISION;
int i;
int min=0;
for (i=0;i<DIVISION;i++){
// theta°だけ回転させて
rotation(home,rotate,anglefrom+theta*i,hnum);
// 平均距離を求める（後述）
distance[i]=calcdistance(rotate,hnum);
if (distance[i] < distance[min])
min=i;
}
return distance[min];
}
void rotation(Point *before,Point *after, double theta,int pnum){
// 回転させます
int i;
for (i=0;i<pnum;i++){
after[i].x=before[i].x*cos(theta)-before[i].y*sin(theta);
after[i].y=before[i].x*sin(theta)+before[i].y*cos(theta);
}
return;
}
double calcdistance(Point *home,int hnum){
// 回転させた後の座標に対して、平均距離を求めます
int i;
int min;
double distance;
　
min=0;
distance=home[0].y;
for (i=1;i<hnum;i++){
distance+=home[i].y;
if (home[i].y<home[min].y){
min=i;
}
}
return distance-home[min].y*hnum;
}

平均距離の測り方は難しそうに見えますが、座標を時計回りにθだけ回転させれば (つまり、直線をx軸と平行にすれば) 簡単に求められます。

さて、このプログラムに問題文のSample Inputを入れてみると、なんと4つのテストケースのうち3つに対して正しく答えを出します。
ただ、2つ目のテストケースに対しては正答 0.000 に対し 18.500 と、大きく異なる値が出てしまいました。いったい何が原因なのでしょう。

いろいろと足りてない

これはアルゴリズムが間違っているのではありません。単に 1000 という分割数が少なすぎただけだったのです。その証拠に、分割数を 1,000,000 にしてみると、0.016 という正答により近づいた値を出力します。
このまま分割数を上げていけば、正答に辿り着くかのように見えます。

が、そうは問屋が卸しません。確かに出力はより正確になったけれども、今度は計算時間が膨大になってしまったのです。
手元の機械では、分割数 1,000,000 として Sample Input に対して出力するのにかかる時間が約15秒。この問題の時間制限を (おそらく) 既に超過しています。
計算時間を短くする。精度を上げる。この2つの問題を同時に解決するために、もう一捻り必要です。

分割、そして分割

闇雲に分割数を上げても上手くいかないことは分かりました。それでは、「答えがありそうなところだけ」細かく分割していくのはどうでしょうか。
つまり、こういうことになります。

分割数 1000 として、平均距離の最小値を求める。
最小値があった角度の前後 (360/1000)°をさらに 1000 分割して、平均距離の最小値を求める。
さらにその前後を 1000 分割し......

：
：
：

そして、答えの精度が高くなってきたら計算を打ち切るのです。
while ループでもいいのですが、今回は再帰で実装してみました。
以下、変更のある関数のみ書き出しました。

const double PRECISION = 1e-6;　　　//　精度
double getdistsum(Point *home,int hnum,double anglefrom,double anglesize){
　Point rotate[MAX_HOME_NUM];
　double distance[DIVISION];
　const double theta=anglesize/DIVISION;
　int i;
　int min=0;
　for (i=0;i<DIVISION;i++){
　　rotation(home,rotate,anglefrom+theta*i,hnum);
　　distance[i]=calcdistance(rotate,hnum);
　　if (distance[i] < distance[min])
　　　min=i;
　}
　// ここから変更点
　if (distance[(min+1)%DIVISION]-distance[min] < PRECISION)
　　return distance[min];
　else
　　return getdistsum(home,hnum,anglefrom+theta*(min-1),theta*2);
　// ここまで変更点
}

こうすると、短時間で見事正答を出してくれるのです。
この問題がオンラインジャッジに追加されたことを知ったのがつい先ほどなのでまだAcceptは受け取っていませんが、おそらくこれでOKに違いありません。
それを数学的に証明するにはもう少しステップを踏む必要がありますが (1000は十分大きな分割数なのか (もしかしたら不十分かも) など)、ここでは言及しないでおきます。

難解な理論に頼らなくても、このように上手く計算力を使ってやれば、簡単に答えが出ることもありますよ。というのが今回のテーマでした。

致命的欠点

しかし、この方法も万能というわけではありません。類似の問題である Secret in Shadow も一見このように解けそうですが、実は出来ないのです。
原因は、10e-10 というべらぼうに高い出力の精度にあります。math ヘッダの sin() 関数や cos() 関数は、この要求に応えてくれるほど精度が高くないのです。
別の方法を探すか、あるいはいっそ自分で超高精度の三角関数を実装するか。すると多倍長浮動小数点数の自作も必要になってきそうです。あまりに険しい山ですが、登ってみるのも一興かもしれません。

担当 : 田山 (define よりも const が好き)