AOJ埋め18

AOJに飽きてきたので今日はコンピュータ・ジオメトリを部室から借りてちょっと読んでいた．やっぱり難しいトピックが多い．
2章の平面走査はなんとなく分かったつもりになれたけどやはりまだあまり分かってないのでもう1回読み返したい．

今日やったのは1題だけ．積分の問題．

2088 Spirograph

hypotrochoidsという曲線の長さを求める問題．
何それ… という感じだったのですぐにググった．実在する曲線らしい．へー．
ある特定のときには楕円になるらしいけど，楕円は周の長さを解析的に求めるのが難しいはずなので，どうやら曲線を折れ線で近似するなり数値積分を実行するなりして値を求める問題らしい．
誰も解いてないっぽいけど，要求精度がとても緩いしなんとかなるんでは? と思って書く．
和を求める際に注意しないといけないとおもったのは，たとえば

double res=0;
for (int j=0; j<M; ++j) res+=f(j);
printf("%.9f\n",res);

のようにするとまずそうだということ．
なぜならこれだと解が大きな値になる場合，resはループに応じてどんどん大きくなるけど，そうなるとres+=f(j);のところで， resはとても大きいけど f(j)はとても小さい ということが起きてしまう．こうなると和の足し算で誤差の積み重ねがどんどん大きくなって最終的には誤差が大きくなってしまいうる．
doubleの精度は15桁程度あるから上の実装でも実は良かったりするのかもしれないけど，不安だったので和を積む際には下のように実装して誤差が重ならないようにした．

//lowerからupperまでの和を求める
double solve(int lower,int upper) {
    if (lower==upper-1) return f(lower);
    int mid=(lower+upper)/2;
    return solve(lower,mid)+solve(mid,upper);
}

..

printf("%.9f",solve(0,M));

で，問題は関数f(j)の部分および関数の周期である．これはWikipediaにある内コトロイドのパラメータ表示に対してsqrt{(dx/dθ)^2+(dy/dθ)^2}を積み重ねればいい．微分が面倒なら，単に折れ線の長さを求めてもいいのかもしれない．
周期は，R=0のときは2π(ただの円)，R>0のときは2nπ (nはP-Q/Qを既約表示したときの分母) とかになる．
問題を解いている最中，実際に和をとるのは0→2πでよくて，対称性からそれをn倍すれば元の値になるはず… とずっと考えていたのだけど，この考えは普通に間違っていて，それだと正しい値は出ない．ちゃんと素直に0→*2nπまで積分しないといけない．

実装すべきコード量自体はとても小さいけどデバッグがしにくくて骨の折れる問題だった．