Ellipse circumference

楽しみにしていた Ellipse circumference が公開になった。まったくかなわなかった notogawa さんの C の解。

float x,a,b;
main(n){
    for(;~scanf("%f%f",&a,&b);printf("%.f\n",x))
        for(n=4e3;--n;)x=a*6.2832+(n-2+.75/n)/n*(1-b*b/a/a)*x;
}

これ、すごいな。2 回目の(外側)ループ以降、x が初期化されていない。なんでだ？と思ったが、内側ループが n=4e3 と高次の項から低次へ向かう方向で計算している。つまり、x の最初の値、つまり、前の計算結果(printfで出力しているもの)は、誤差範囲に葬り去っているということだ！まじかぁ。すごいな。
＃＃＃とりあえずここまで。仕事中なので、残りはあとで書こう。。。

2/2追記：

notogawa さんご自身がすでに詳しい解説をされていた。⇒ http://d.hatena.ne.jp/notogawa/20110201/1296557658
C で、n を 8e3 から 2 ずつ減らす方法でさらに短くならないか？とやってみたが、結果は同じ 119B だった。まぁ、この程度の式変形はすでにケアされているか。。。

float x,a,b;main(n){for(;~scanf("%f%f",&a,&b);printf("%.f\n",x))for(n=8e3;n-=2;)x=a*6.2832+(n-4+3./n)/n*(1-b*b/a/a)*x;}

追記：

っと思ったら 1B 縮んだ。