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

桁落ちと情報落ち

ここで今更(前回の説明で盛大に使っている)ではありますが「桁落ち」と「情報落ち」について確認します。

桁落ち

有効数字に制限のある計算において、たとえば、

141.4214 - 141.4213 = 0.0001

のように、値が何桁も一致する数の減算を行った結果、有効数字が(場合によってはほとんど)残らなくなるのが「桁落ち」です。

情報落ち

一方、

123.456 + 0.123456 = 123.579

のように、値が何桁も離れた数の加減算の結果、絶対値が小さい方の数の下の桁の情報が無視されるのが「情報落ち」です。小さい方の数があまりに小さいと、やはり完全に無視されることになります。
セオリーとしては、桁落ちは必ず避け、計算の最後で最小限の情報落ちが起こる(意味的に考えて、無視してよい誤差であることが多い)ようにするのが一般的です。

2進による mod

前回説明した、「より正確な mod」を、実際に計算機で計算できるコードにすると、次のようになります。コード全体を示した後、順番に説明してゆきます。ここでは Ruby を使っていますが、普通のプログラミング言語であれば、どれでも同じようにして実装できます。

precise_divrem

Ruby の divmod メソッドのように、整数の範囲での商と、剰余を求めるメソッドです。参考文献の p.112 にある Pascal のコードから、変数名 p を避けるなどの変更を加えて Ruby に書きかえたものです。
2進の、引きっ放し法による除算アルゴリズムのほぼそのままですが、以下簡単に解説します。
x から y をいくつ取ることができるか、という取り尽くしによる除算です。
まず、y の半分を z とします。そして、絶対値が z より小さくなるまで、x をどんどん半分にしてゆきます(何回、半分にしたかを覚えておきます)。
この後、正の場合と負の場合がある処理になりますが、まず、正の場合で説明します。
絶対値が z より大きくなるまで x を倍にし、z (すなわち y/2) より大きくなったら x から y を引きます。この時 x は y より小さいので、x は負になります。
負の場合も同様に、絶対値が z より大きくなるまで x を倍にしてゆきます。絶対値が z より大きくなったら、x に y を足します。これで正に戻ります。
ここの手続きを図にすると、次のようになります。

	           y を引く
	     x <------------ x

        y を足す 
   x ------------> x

-y	-z	0	z	y
|	|	|	|	|
+BBBBBBB+-------+-------+AAAAAAA+

-z~0~z の範囲に x がある場合は、x を倍にします。A の範囲に入ったら y を引き、B の範囲に入ったら y を足します。
最初に x を半分にしたぶんを、後半の x を倍にする操作で、取り戻したら終わりです。この時 n に、x からいくつ y を取ることができるか、すなわち整数の範囲の商が得られていて、x にはその残り、つまり剰余が得られます。
y の加算や減算で、「情報落ちは一切起こらず、桁落ちは起きている」ということに注目してください。ここで情報落ちが起きていないために、後の計算で回復させることができるような、有効数字より下の桁の情報が残るわけです。一方で桁落ちにより、下の桁に 0 が補われますが、これは、x が、存在しているよりも下の桁には 0 が続いている正確な数だとすれば、正しい操作だということになります。このように、通常のセオリーと逆に、「情報落ちを避け、コントロールされた桁落ちを起こす」のが、この場合は正解なのです。
なお、mod ではなく rem という名前を使っていますが、これは、他の多くの浮動小数点数の mod(特に IEEE 754 の mod)と違い、負の剰余を返す場合があるためです(当然、商として、負の剰余に整合した値を返します)。IEEE 754 が定義する reminder に似ているので rem という名前にしましたが、丁度半分の時に常に負の剰余を返す点が、IEEE 754 の reminder とも異なります(参考 http://stackoverflow.com/questions/158120/math-ieeeremainder-returns-negative-results-why )。

実際に計算してみる

以上で、必要な特性を持つ divmod が得られましたので、実際に計算に使ってみます。ここでは、1000000000.0 `mod` 2π を求めてみました。
まず、

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

IEEE 754 の倍精度を前提としていますが、補正用の、Math::PI の打ち切られた残りの値です。浮動小数点数の十六進による記法を使っています。
続いて、

def my_divmod x, y
  n = (x / y).floor.to_i
  return n, x - n.to_f * y
end

比較用の「手抜き」divmod です。

puts "bigdecimal"
bd_x = BigDecimal("#{x}000000000_0000000000")
bd_pi = BigDecimal("3.1415926535_8979323846_2643383279")
bd_2 = BigDecimal("2.0000000000_0000000000_0000000000")
puts (bd_x % (bd_2 * bd_pi))
puts ""

正しい結果が得られているか確認するために、BigDecimal による多倍長計算を使います。
手元の環境(FreeBSD8、ruby 1.9.3)での出力結果は次のようになりました。

bigdecimal
0.577395423501385169409807203806E0

precise_divrem
0.5773954235013852

my_divmod
0.5773954001662309

Float#divmod
0.5773954235013852

precise_divrem を使って、正しい結果が得られています。当然ながら手抜き divmod ではまともな結果が得られていませんが、注目すべきは処理系組み込みの Float#divmod でも正しい結果が得られていることでしょう。
rubyソースコードを見ると( http://svn.ruby-lang.org/cgi-bin/viewvc.cgi/trunk/numeric.c?annotate=40437#l868 )プラットフォームに fmod があればそちらを、無ければ mod = x - z * y というようになってますから、プラットフォームの fmod により正しい結果が得られているようだ、ということになります。
IEEE 754 の規格票を持ってないので、そちらの正確なことはわかりませんが、)man fmod や POSIXhttp://pubs.opengroup.org/onlinepubs/007904975/functions/fmod.html )では x - i*y を返す、というところまでしか規定されていませんから、これは実装が独自(?)にやっていることと言えます。FreeBSD の libm は Sun 作の FDLIBM( http://www.netlib.org/fdlibm/readme )がベースですから、そちらのソースを見てみますと( http://www.netlib.org/fdlibm/e_fmod.c / http://svnweb.freebsd.org/base/head/lib/msun/src/e_fmod.c?view=markup )かなりややこしいことをやっていますので、おそらく本稿で示したような正確な値を返すよう実装されているのだと思います。
また、余談になりますが、FDLIBM の三角関数ルーチン( http://www.netlib.org/fdlibm/k_sin.c )では、2 個目の引数として、この補正用の値を直接受け取るように作られています( __kernel_sin(double x, double y, int iy) の y がそれ)。

参考文献

近似式のプログラミング

近似式のプログラミング

"Hacker's Delight" より

I think that I shall never envision
An op unlovely as division.

An op whose answer must be guessed
And then, through multiply, assessed;

An op for which we dearly pay,
In cycles wasted every day.


Division code is often hairy;
Long division's downright scary.

The proofs can overtax your brain,
The ceiling and floor may drive you insane.

Good code to divide takes a Knuthian hero,
But even God can't divide by zero!