Pixel Pedals of Tomakomai

北海道苫小牧市出身の初老の日常

トップ > 技術 > 小数第四位までの小数を10000倍して整数に直せるのか

小数第四位までの小数を10000倍して整数に直せるのか

技術

数値が小数第四位まで与えられるとする。例えば、 1234.5678 のような形式である。

それらの数値を使って計算をしたいのだが、何も考えないと計算機上ではこれらの数値は double 型などで保持され、誤差が発生することは有名だろう。

>>> 1001.0000 * 1.1000 == 1101.1000
False
>>> 1001.0000 * 1.1000
1101.1000000000001

誤差が発生するのは嫌なので、入力を 10000 倍して整数として扱いたくなる。例えば、 1234.5678 は 12345678 に変換したい。

しかし、この操作にも注意が必要である。 https://qiita.com/mod_poppo/items/910b5fb9303baf864bf7 で解説されている通り、 10000 倍したり、 int を取ったりするだけではうまくいかない。

>>> 1.2347 * 10000
12346.999999999998
>>> int(1.2347 * 10000)
12346

ということで、四捨五入することになる。

>>> int(1.2347 * 10000 + .5)
12347

ところで、 この四捨五入する方法は果たして万能なんだろうか 。答えから言うと、残念ながら万能ではなく、桁数が大きくなるとうまくいかなくなる。

>>> int(274877906944.1832 * 10000 + .5)
2748779069441833

int と .5 を使ったナイーブな実装だからではという人もいるかもしれないが、 round を使っても別の値でうまく行かない。

>>> round(274877906944.1832 * 10000)  # うまくいった！
2748779069441832
>>> int(274877906944.1835 * 10000 + .5)  # うまくいった！
2748779069441835
>>> round(274877906944.1835 * 10000)  # 駄目じゃん・・・
2748779069441834

この差は浮動小数点を「丸める」という仕様がそもそも人によって変わるということに起因している1。以下の 8 年ほど前のエントリに詳しい。

http://blog.practical-scheme.net/shiro/20131229-flonum-rounding

こんな感じのスクリプトで調べると、整数部が 274877906943 では発生しないようだ。

for i in range(1, 10000):
    x = f"274877906943.{i:04d}"
    ans = 274877906943 * 10000 + i
    if int(float(x) * 10000 + .5) != ans:
        print(f"int NG: {x}")
    if round(float(x) * 10000) != ans:
        print(f"round NG: {x}")

この現象について、もう少し数学的に考えてみよう。 double のフォーマットは色々なところに書かれていると思うが、今回は以下を参考にした。

https://www.k-cube.co.jp/wakaba/server/floating_point.html

誤差の伝搬による限界

URL に書かれている内容によると、 double の仮数部は 52bit あり、最高位の隠しビットを含めると 53bit で表現されていると言える。 52bit のうち整数部に n bit を使うとすると、小数に使えるのは 52 - n bit であり、一番細かい桁は 2^ {- 52 + n} ということになる。真の値から一番近い表現に丸められるとすると、誤差は最小位の桁の半分となるため 2 ^{-53 + n} となる。

10000 倍をしたときに、この誤差が \frac{10^ {-4}}{2} 未満であれば round したときに意図する値に丸められると言える。つまり 2 ^{-53 + n} \lt \frac{10^ {-4}}{2} が条件であり、解くと n \lt 52 - 4 \frac{\log 10}{\log 2} \lt 38.7... となる。

つまり、整数 bit が 39 以上だと問題が発生することがわかる。しかし、先の数 274877906944 は 2^ {39} = 549755813888 よりも小さく、これだけではうまく説明がつかない。

掛け算後の桁あふれによる限界

先程の計算では整数部が n bit であるとして計算したが、 10000 倍した後には整数部の bit 数は増えるはずで、これにより小数 bit が少なくなり、値を丸めるために誤差が発生すると考えられる。簡単のため 2 、 \log _ 2 10000 桁ズレるものとしよう。これにより、小数 bit の最小桁は 52 - n - \log _ 2 10000 となり、この桁に丸めるために \frac{2^ {- 52 + n + \log _ 2 10000}}{2} = \frac{10000 \cdot 2^ {-52 + n}}{2} の誤差が新たに加わる。

先の誤差と合計すると、 10000 \cdot 2 ^{-53 + n} + \frac{10000 \cdot 2^ {-52 + n}}{2} = 10000 \cdot 2^ {-52 + n} であり、ここから先程と同様に不等式を作って解くと n \lt 51 - 4 \frac{\log 10}{\log 2} \lt 37.7... となり、整数 bit が 38 以上だと問題が起こることがわかる。このときの値は 2^ {38} = 274877906944 で、 10000 倍して round してもうまくいかない例で紹介した値と一致している。

うまく整数にできない例

実際に何が起こっているか、もう少し見てみよう。 Haskell では 274877906944.0007 を 10000 倍してうまく整数化することができない。

Prelude Data.Ratio> round $ 274877906944.0007 * 10000
2748779069440006

Rational は有理数で正確にリテラルの値を表現できるので、これとの差をとってみると、この時点で真の値との誤差が 0.000029 程度あることがわかる。リテラルを見るだけでは想像がつかない、割と大きな値なのではないだろうか。

Prelude Data.Ratio> x = 274877906944.0007 :: Rational
Prelude Data.Ratio> y = 274877906944.0007 :: Double
Prelude Data.Ratio> d = toRational y - x :: Rational
Prelude Data.Ratio> fromRational d
-2.861328125e-5

ただ、ここまでであれば 10000 倍して小数点を四捨五入すれば元の数に戻るはずである。 10000 倍後の誤差は、悲しいことにちょうど 0.5 となる。

Prelude Data.Ratio> x = 2748779069440007 :: Rational
Prelude Data.Ratio> y = 274877906944.0007 * 10000 :: Double
Prelude Data.Ratio> d = toRational y - x :: Rational
Prelude Data.Ratio> fromRational d
-0.5
Prelude Data.Ratio> d
(-1) % 2

もう少し状況を見てみよう。 Double で表された元の値を見ると、このような有理数になっている。

Prelude Data.Ratio> y = 274877906944.0007 :: Double
Prelude Data.Ratio> toRational y
4503599627370507 % 16384

手元に Haskell の Decimal パッケージがなかったので横着して Python で試してみると、この値は 10 進数で以下の通り。先程確認したとおり、真の値より 0.000029 少ない値だ。

>>> Decimal(4503599627370507) / Decimal(16384)
Decimal('274877906944.00067138671875')

これを 10000 倍するのは人間であれば簡単で、 2748779069440006.7138671875 である。ここまでで発生しているのが「誤差の伝播による限界」で計算したものであり、 10000 倍されて真の値と比べると 0.29 の誤差が発生することにある。

さて、このまま四捨五入できればよかったのだが、 この値は 2^ {51} と 2^ {52} の間にある。

Prelude Data.Ratio> (2^51, 2^52)
(2251799813685248,4503599627370496)

よって、整数のために仮数の bit が 51 必要であり、小数に当てられた bit は 1 bit しかない。そのため、この値を 2748779069440006.5 か 2748779069440007.0 に丸めなければならない。近い方を選択すると前者になり、 round の偶数丸めの仕様から 2748779069440006 と違う値に変換されてしまう。ここで生じた誤差は「掛け算後の桁あふれによる限界」で述べたもので、誤差がマイナス方向に 2 度蓄積することで希望する値に丸められていないことがわかった。

うまく整数にできる例

逆に整数部が 274877906943 のときは何が起きているのだろうか。 274877906943.9944 を整数にするときの処理を詳しく見てみよう。

まず、 Double にした段階での誤差を見ると 0.000015 であり、先程より僅かに小さい。

Prelude Data.Ratio> x = 274877906943.9944 :: Rational
Prelude Data.Ratio> y = 274877906943.9944 :: Double
Prelude Data.Ratio> x - toRational y
39 % 2560000
>>> Decimal(39) / Decimal(2560000)
Decimal('0.000015234375')

現在整数 bit は 37 bit であり、誤差は、

Prelude Data.Ratio> 2 ** (-15.0) / 2.0
1.52587890625e-5

未満であるから、現在の整数の bit 数では最大に近い誤差である（わざとそういう数を選んだ）。

さて、この Double の値の正確な数字を見ておこう。

Prelude Data.Ratio> toRational y
1125899906842601 % 4096
>>> Decimal(1125899906842601) / Decimal(4096)
Decimal('274877906943.994384765625')

10000 倍すると、 2748779069439943.84765625 である。整数部 2748779069439943 は先程と同様に 51bit を要する数であり、小数に充てられる bit 数は 1 bit のみである。これを round することになるが、先程と違うのは元の誤差が 0.0000153 未満であることがわかっているので 10000 倍しても 0.153 であり、 0.25 よりも小さいので正確に元の数字の方に丸められるということである。 実際、今回は 2748779069439943.5 か 2748779069439944.0 に丸めなければならないが、 2748779069439944.0 へ正確に丸められることになる。

  1. shiroさんのエントリを紹介したかったので書いてしまったけど、そもそも四捨五入と偶数丸めで仕様が根本的に違うものを比べているので当たり前ではあった。↩

  2. 本来は動かす桁数が小数になることはないので、値の大きさによって場合分けして整数桁ずらすべき。↩