tan(355/226)

(追記: 「正確さ」と書くべき場所の多くで「精度」となっていたのを修正しました)
ウェブサイト『関数電卓マニアの部屋』( http://teamcoil.sp.u-tokai.ac.jp/calculator/index.html )を見ていて知ったのだが、torture test という関数電卓の限界的なスペックを調査しているサイトがあり、そのテストの #1 が trig accuracy というもので、accuracy of tan(355/226) すなわち tan(355/226) の正確さをチェックしている( http://www.voidware.com/calcs/torture_trigs.htm )。
ここで 355/226 は、円周率πの近似として知られる 355/113 の半分で、π/2 よりわずかに大きい数である。このため tan(355/226) は絶対値がそこそこ大きい負の値になる。
具体的には、torture test のサイトが示している値では、-7497258.18532558711290507183... という値である。WolframAlpha による計算でも確かめられる( http://www.wolframalpha.com/input/?i=tan%28355%2F226%29 )。
これを ruby で Math.tan(355.0/226.0) のように計算しても、-7497258.179140373 と、かなり低い正確さの結果しか得られない。
IEEE 754の倍精度で、可能な限り正しい値が得られたならば、-7497258.185325587 となるはずである。以下、この誤差の原因を追跡する。
まず 355.0/226.0 は 1.5707964601769913 という値になる。浮動小数点数の16進表示では、0x1.921fb78121fb8p+0 である。
この前後での Math.tan の結果を見てみる。

Math.tan(Float("0x1.921fb78121fb7p+0")) #=> -7497258.191621251
Math.tan(Float("0x1.921fb78121fb8p+0")) #=> -7497258.179140373
Math.tan(Float("0x1.921fb78121fb9p+0")) #=> -7497258.166659494

割り切れない除算のため誤差があるはずだが、この付近の値では、引数で表現可能なぎりぎりの値の違いに対しても関数値が大きく変わるため、うまく計算できていないようだということがわかる。

つまり(検索するとちらほら見つかるのだが)この計算を「三角関数の正確さ」のチェックだと考えてはいけない、ということである。ここで示した私の例の場合、(仮にtanが完璧だったとしても)それに与える引数の時点で精度が足りていないために、不正確さはその時点で発生している。

tan(x) == cot(π/2 - x), cot(x) == 1/tan(x) より、
tan(x) == 1/tan(π/2 - x) であるので、これで計算してみる。
まず、π/2 - 355.0/226.0 は、普通に計算しただけだと、桁落ちのためにたいした正確さが出ない。
ruby で、次のような定数やメソッドを定義して、それらを使う。
定数

PI_D = Float("0x1.1a62633145c06p-53")

これは、π(真値) == Math::PI + PI_D であるような値。
次のメソッド divmod2 は、

divmod2 x, y #=> a, b

このように呼び出すと、Float で表現できる限りにおいて a = x / y を計算したうえで、a と、その余り b を返す。すなわち(仮に、より正確に計算できたなら)x == y * a + b となる(やっつけなので、ゼロ除算対策などはしていない)。

これらを使って π/2 - 355.0/226.0 を計算すると、

a, b = divmod2 355.0, 226.0
x = (Math::PI*0.5 - a) + (PI_D*0.5 - b/226.0)

となる。この x を使って、

p 1.0 / Math.tan(x)

とすると -7497258.185325587 という、満足できる結果が得られる。