2013/07/26

実数部が0にならない話

• 関数電卓再び

「-2の1.5乗」は純虚数になるはずで、HP 35sでは実際そうなるのだけど、 プログラミング言語処理系では実数部に誤差が乗ってしまう、という話。 Gaucheもそうなる。

gosh> (expt -2 1.5)
-5.195564742431613e-16-2.82842712474619i

(expt x y)は、結果が「明らかに実数」でない場合は次のように計算している。

    x^y == e^{y * log(x)}
        == e^{y * log(|x|)} * e^{y * arg(x) * i}

xが負の実数ならarg(x) = πなので、

    x^y == e^{y*log(-x)} * e^{πyi}
        == e^{y*log(-x)}*cos(πy) + e^{y*log(-x)}*sin(πy)i

y = 1.5の時実部が0にならない理由は、1.5πをdoubleで正確に表現できないから cos(1.5π)が0にならないってことなんだけど、なんかそれもださい気がする。 1.5自体はdoubleで正確に表現できるんだし。

sinpi(x) = sin(x*π) といった関数があれば、 xがキリの良い数(1/2の倍数)の時に正確な値が出せそうだ。 浮動小数点数演算でもsin(0), cos(0)は常に正確なので、 1/2πの倍数が0にくるように入力を適切に変換してやればよい。 そんな感じで作ってみたのがこれ。

gosh> (for-each (^x (print (%sinpi x))) (lrange -2 2.001 1/8))
0.0
0.3826834323650898
0.7071067811865475
0.9238795325112867
1.0
0.9238795325112867
0.7071067811865475
0.3826834323650898
0.0
-0.3826834323650898
-0.7071067811865475
-0.9238795325112867
-1.0
-0.9238795325112867
-0.7071067811865475
-0.3826834323650898
0.0
0.3826834323650898
0.7071067811865475
0.9238795325112867
1.0
0.9238795325112867
0.7071067811865475
0.3826834323650898
0.0
-0.3826834323650898
-0.7071067811865475
-0.9238795325112867
-1.0
-0.9238795325112867
-0.7071067811865475
-0.3826834323650898
0.0
#<undef>

これを使うと実部の誤差はこの通り消える。

gosh> (expt -2 1.5)
0.0-2.82842712474619i

確か以前にもsin(nπ)あたりが綺麗に出ないのに残念な気分に なった覚えがあるんだけど、少なくともGaucheのソースでは 他にsinpiとかが使える場所は無い感じなので、試しに入れてはみたが キープするかどうかは微妙。しばらく様子見。

Tag: Gauche