浮動小数点数の(より)正確な mod (その1)

浮動小数点演算で mod なんか計算するのか? と思うかもしれませんが、たとえば sin(10000.0) のような計算をする時、よくできている浮動小数点演算の処理系では、10000.0 `mod` 2π を計算してから、近似式で結果を求めています(細かく言うと、最終的には -π/4~π/4 に落とし込む)。さもなくば、単にでたらめな結果を返すか、エラーとなるでしょう。0 から遠く離れても正確に近似するような近似式は、無いからです。
昔話をすると、8 ビットパソコン時代に、学校の文化祭などでパソコンに計算させて結果をグラフィカルに表示する題材として、「バイオリズム」なるものが流行ったものでした。誕生日を基点として、一週間前後の周期で気分曲線などがサインカーブを描くとかいったもので、果たしてそんな長期にわたって、数学関数のように正確に気分が浮き沈みするものかどうか、考えてみれば甚だ怪しいシロモノでした。果たして当時の Basic がどこまで正確な計算をしていたのかも微妙なところです。
話を元に戻して、この時、10000.0 `mod` 2π について、単純な計算では、実は正確さが落ちます。以下では、より正確な値を求める方法を解説します。やはり、よくできている浮動小数点演算の処理系では、ここで述べているような正確さを向上させる計算を内部で行っているはずです。
まず、10000.0 は近似値ではなく、正確にその値だとします。つまり、意味的には 10000.000... とどこまでも 0 が続くものと考えます。また、実際に計算機に実装する際は、倍精度などの二進浮動小数になりますが、ここでは説明のため、十進で有効数字 7 桁の浮動小数で計算します。
十進で 2π の値は、6.283185307...です。有効数字 7 桁ですから、8 桁目の 3 以下を四捨五入で切り捨てて、6.283185 です(あとの計算が理解しやすいように、切り捨てになる桁を選びました。実際には切り上げでも計算できます)。
まず、10000.0 `mod` 6.283185 の計算にも注意が必要です。最初に、よくある、除算結果の整数部分を乗じて引く方法でやってみます。

10000.0 ÷ 6.283185 = 1591.5495...

ですから、

10000.0 - ( 1591.0 × 6.283185 )
= 10000.0 - 9996.547
= 3.453000

となります。ここで、(1)引くべき値を求めるための乗算をした結果、桁数が伸びて、有効数字の制限のために正確な下の桁の値が失われたこと (2)近い数の減算をした結果、桁落ちが起きて、結果の下のほうの桁が全く無意味なこと、の 2 点に注意してください。実際に、0.001 程度の誤差があり、有効数字が 4 桁もありません。
この計算を次のようにすると、これよりは正確な結果が得られます。

10000.0 - 6283.185 = 3716.815
3716.815 - ( 5.0 × 628.3185 )
= 3716.815 - ( 4.0 × 628.3185 ) - 628.3185
= 3716.815 - 2513.274 - 628.3185
= 1203.541 - 628.3185
= 575.2225
575.2225 - ( 9.0 × 62.83185 )
= 575.2225 - ( 8.0 × 62.83185 ) - 62.83185
= 575.2225 - 502.6548 - 62.83185
= 72.56770 - 62.83185
= 9.735850
9.735850 - 6.283185
= 3.452665

姑息な変形ですが、有効数字を拡張して計算したりはしていないのに(末尾の桁がたまたま 5 だったので、偶然うまくいっているという面もありますが)情報が失われていないことに注意してください。
実は、計算途中で桁落ちそのものは起きているので、絶対的な有効数字としてはたいして変わっていないのですが、この末尾の桁には、ある情報が残っています。
最初の、2π の真の値から切り捨てた値のうちの、(ここでの有効数字のぶんである)上位 7 桁 0.0000003071796 を考えます。2π の最初の近似値 6.283185 = a, 0.0000003071796 = b とすると、10000.0 `mod` 2π の、より正確な値は 10000.0 `mod` ( a + b ) ということになります。これを以下のように変形してゆくと、

10000.0 - ( 1591 * ( a + b ) )
= 10000.0 - ( 1591 * a ) - ( 1591 * b )

10000.0 - ( 1591 * a ) の値として、先ほど正確に計算した値を使い、

3.452665 - ( 1591.0 * 0.0000003071796 )
= 3.452665 - 0.0004887227
= 3.452176

このようにして補正することで、有効数字一杯までの正しい結果を取り戻すことができました(たとえば Ruby で、10000.0 % ( 2.0 * Math::PI ) を計算して、比較してみてください)。
次回はこれを2進でどのようにおこなうかを示し、計算機における計算で試してみます。

参考文献

近似式のプログラミング

近似式のプログラミング