2013/12/29

2進小数の10進桁での丸め

• RubyとPythonとC#のround関数のバグっぽい挙動について

5.015を小数点以下2桁のところで丸めたらいくつになるべきか、という話。 小数に2進浮動小数点数を使っている場合、5.015は正確に表現できず、 もっとも近い表現可能な値は5.015よりわずかに小さな値となる。 とすれば、丸めた結果は近い方である5.01になるべきだろうか、 それとも元の表記から考えられる5.02になるべきだろうか。

この問題には、唯一の正解はなさそうだ。 そもそも5.01もしくは5.02も2進浮動小数点数で正確に表現できないので、 結果も浮動小数点数で得るとすれば、これは「丸め」操作ではないのだ (2進で考えた場合、桁数は減っていない)。 「丸め」操作なら、「二重丸めはしてはいけない」という原則があるけれど、 そもそも丸めでないのでこの原則が適用できるかどうかも曖昧である。

ただ、モデルを定めれば、そのモデルでのあるべき解というのは考えられると思う。

以下、混乱を避けるために、正確な数を [5.015]と表記する。 一方、(NaN, 無限大を除く)2進浮動小数点数はそれぞれ正確な数値と対応しているけれど、 10進表記で書くと長大になるので、通常「曖昧さの無い範囲で最短の表現」が取られる。 それを{5.015}と書こう。どちらか解釈が確定してない場合は普通に5.015と書く。

つまり{5.015}というのは[5.015]にもっとも近い2進浮動小数点数で、 倍精度であればその正確な値は [5646388032815759/1125899906842624] である。 (この有理数は、倍精度浮動小数点数の定義から計算できる。Gaucheの場合、 0.9.3.3以前であれば (inexact->exact 5.015)、 現在の開発版であれば (real->rational 5.015 0 0 #f) とする。)

ではいくつかモデルを考えてみよう。

1. 浮動小数点数は、その2進数値が表現する正確な値に対応していると考えるモデル

5.015と言っている数は厳密に {5.015} = [5646388032815759/1125899906842624] を示している、 もしくは{5.015}を代表値として±ε/2の誤差を含んでいる、と考える立場。 数値計算ではこっちが普通だと思う。

• 1a. {5.015}は[5.015] より少し小さい。 従って[5.02]ではなく[5.01]に丸めるべき、と考える。
  [5.01]も2進浮動小数点数で正確には表現できないので、 もっとも近い浮動小数点数 {5.01} を答えとして採用する。 これはdoubleなら[2820379266640773/562949953421312]に等しく、 [5.01]より少し小さい。

• 1b. 「近い方へ丸める」には近さを判定しなければならないが、[5.01]と {5.015} の距離は2進浮動小数点演算では正確に計算できない。 1a.の方法は、距離の判定で近似演算を行い、さらに結果を近似しているので、 二重丸めの臭いがする。
  そこで、比較演算に使う先もあらかじめ正確値を使ってみよう。 つまり、{5.015}と{5.01}、{5.02}をそれぞれ比べるのだ。

  gosh> (real->rational 5.01 0 0 #f)
  2820379266640773/562949953421312
  gosh> (real->rational 5.02 0 0 #f)
  1413004383087493/281474976710656
  gosh> (real->rational 5.015 0 0 #f)
  5646388032815759/1125899906842624
  gosh> (- 5646388032815759/1125899906842624 2820379266640773/562949953421312)
  5629499534213/1125899906842624
  gosh> (- 5646388032815759/1125899906842624 1413004383087493/281474976710656)
  -5629499534213/1125899906842624
  

  この通り、{5.015}は{5.01}、{5.02}のちょうど中間にある。 「偶数丸め」原則を結果の10進表記にあてはめて考えるなら、{5.02}が答えとして採用される。

2. 浮動小数点数は、その10進表記が本来表したかった数であって、正確に表現できない場合は便宜上もっとも近い値で近似している、とするモデル

つまり、5.015というのは[5.015]という意図なんだけど、 計算機では表現できないから仕方なく{5.015}を代わりに使ってる、と考える。

こちらは日常の感覚に近い。我々は10進表記を見たら、「その数字は真の値とは違う、仮の表記」 と考えるよりは、「その数字が真の値、計算機内部の表現が仮の表現」と考えるだろう。

ややこしいが、この違いは誤差を考慮しても現れる。下図において、5.015という表記を 見た時に日常的にとらえる感覚が上、計算機に適した解釈が下である。 (誤差分布をガウシアンみたいに書いたけど、多くの場合は一様分布の方がふさわしいかも)。

この場合、5.015は本来[5.015]なんだから、小数点数以下2桁で丸めたら 難しいことを言わずに結果の表記が5.02になって欲しい。(その5.02が内部で どう表現されるかは計算機が良きにはからえ)。

このモデルで計算する場合は、計算機内部で{5.015}を[5.015]に 変換したうえで、10進数としての処理を行う必要がある。もちろん[5.015]は 2進浮動小数点数では表現できないので、例えば[5015/1000]として扱うのだ。

ちとややこしいのはこの変換が表記に左右されることで、 例えばもともとの数値が5.01534と表示されていたら人間は自然に[5.01534]である と認識するだろう (ユーザが有効数字をimplicitに想定するということ)。 なので計算機は「10進表記で最適な表示をした場合の桁数を算出して、それが表す 正確な数を求める」という操作をしないとならない。

この点では、「一度文字列にして操作する」のようなナイーブな実装が案外的を射ていたりする。

このモデルでは、日常的なアプリケーションでユーザの直感に反する振る舞いは出てこないだろう。

ただ、10進表記に直す時に一回近似が起きてて、そこからさらに丸めるってのが 二重丸めくさくて怪しい感じはする。

Gaucheで実装するとしたらどうしよう?

「10進表記の小数点以下何桁で丸める」というのは計算途中ではなく ほぼ常にユーザへの呈示の時に使われると思うので、 ユーザにとって混乱の少ない表記になるのが望ましいだろう。

しかし、Gaucheの浮動小数点数の扱いは原則として1.のモデルを取っている。 つまり浮動小数点数を正確に解釈する場合は、定義に従ってその示す値を使う、ということだ。 なのでここだけモデルを変えるのも気持ち悪い。

いやまてよ、そもそもこの操作は「10進表記に直す」時にしか問題にならないのだから、 double, int -> double であるような「10進n桁で丸める」関数を考えるのが間違いなのかもしれない。結果が文字列で得られれば良いのだから、 double, int -> stringという関数にすべき、つまり round関数のバリエーションとするのは誤りで (だってdoubleのドメインでは丸めじゃないもの)、number->stringのオプションとして実装すべきかも。

number->stringで桁数を指定するのは、Gaucheが使ってる Burger&Dybvigの論文にちゃんと誤差最小にする方法が出てるので、 素直にそれを実装するのが正解かな。

(追記2013/12/31 02:33:28 UTC):いやまてよ、Burger&Dybvig使って小数点以下2桁で打ちきると、{5.015}は[5.01]に近いから誤差最小で5.01になっちゃうのかな? 要調査。 もしそうなら、必要に応じて1bのアルゴリズムを別途用意すべきかも。

過去の浮動小数点数に関する話題:

• 実数部が0にならない話
• 拡張浮動小数点数の扱い
• またまた浮動小数点数話
• 浮動小数点数遊び
• 浮動小数点数

Tags: Programming, Gauche