Pythonの四捨五入で四苦八苦した話。（その3）

前回のあらすじ：Rubyの四捨五入すごい。PyPIにあったパッケージどれもダメ。

Rubyでの四捨五入のアルゴリズム

さて、そんなわけで、Pythonがどうしようもなかったので、なぜか誤差なく四捨五入できてたRubyでの実装を調べてみた。

Rubyでの四捨五入は、numeric.cのflo_round()で実装されている。

flo_round()ではいくつかの難しいケースへの対処も入っていて、通常処理の本体はマクロROUND_CALL()から呼び出されるround_half_up()。これを解読してPythonで書き直すと、次のようになる：

def myround3(number, ndigits=0):
    shift_amount = math.pow(10, ndigits)

    if number >= 0:
        shifted = math.floor(number * shift_amount)
        bound = (shifted + 0.5) / shift_amount
        if number >= bound:
            shifted += 1
    else:
        shifted = math.ceil(number * shift_amount)
        bound = (shifted - 0.5) / shift_amount
        if number <= bound:
            shifted -= 1

    return shifted / shift_amount

これを実際に試してみると、次のような感じ：

myround3(0.5)  # => 1.0, OK
myround3(1.5)  # => 2.0, OK
myround3(-0.5)  # => -1.0, OK
myround3(-1.5)  # => -2.0, OK
myround3(12345.67, -1)  # => 12350.0, OK
myround3(-12345.67, -1)  # => -12350.0, OK
myround3(1.255, 2)  # => 1.26, OK
myround3(-1.255, 2)  # => -1.26, OK

基本的な四捨五入はもちろん、鬼門だった1.255を小数第三位で四捨五入するのも正しくできてるのが分かる。

では、なぜこのアルゴリズムだとうまくいくのか？

ポイントは、シフト前の状態で判断しているということ。

これまでのアルゴリズムだと、1.25xxxをシフトして125.xxxにし、そこに0.5を足して126を超えるかどうかで切り上げ/切り捨てを判断してた。つまり、シフト後の状態で判断をしていたことになる。けど、これだとシフトで誤差が露呈してしまい（1.255をシフトすると125.499...になる）、うまくいかなかった。

一方で、Rubyでのアルゴリズムは1.25xxxを四捨五入するときの境界値（この場合は1.255）を作り出し、これとの比較をすることで切り上げ/切り捨てを判断している。つまり、シフト前の状態のままで判断をおこなっている。こうすると何が嬉しいかというと、比較のときに元の値をそのまま使えるので、計算での誤差がそこに乗らない。

もちろん、境界値を作る計算では誤差が乗るんだけど、面白いのは、誤差が乗ったとしてもその誤差は元の数値と同程度にしかならないということ。この非対称性が面白い。

これは具体的に見てみるとよく分かる：

1.255 * 100  # => 125.49999999999999, 誤差が露呈する
125.5 / 100  # => 1.255, 誤差が乗ってない（ように見える；内部では1.255と同じ誤差が乗ってる）

このように、シフトには誤差が乗るけど、アンシフトには誤差が乗らないという性質をうまく使ったのがRubyでの実装となっていた。いやー、賢い・・・

さらに厳しいケースへ

このRubyのリポジトリでは単体テストも用意されてて、RSpecで書かれたspec/ruby/core/float/round_spec.rbおよびtest-unitで書かれたtest/ruby/test_float.rb, test/ruby/test_numeric.rbなどがある。

これらの中にはバグチケットの番号とともに書かれた厳しいテストケースなどもある。

たとえば、「42.0を小数点以下308桁になるように四捨五入したら42.0になること」みたいな。これをさっきのmyround3()でやってみると、こうなる：

myround3(42.0, 308)  # => OverflowError: cannot convert float infinity to integer

これは、10の308乗もの数を42.0に掛けるので、浮動小数点数で表現できる桁を超えてしまって（オーバーフロー）、math.floor()が切り捨てできないよってなるのが原因。こうなる場合、四捨五入する桁はずっと小さい桁になってるから元の数値に影響を与えないことになるので、元の数値をそのまま返すことになる。そして、オーバーフローが発生するかどうかは、ざっくり次のようなロジックで書ける：

def will_overflow(number, ndigits):
    _, num_binexp = math.frexp(number)
    num_exp = num_binexp / 3.322259  # log_2(10) = 3.322259
    return (num_exp + ndigits) >= sys.float_info.dig

sys.float_info.digは浮動小数点数を10進数で表現したときに信頼できる桁数を表す定数。つまり、この桁数を $K$ とすると、 $10^K$ 以上の数は表現できないことになる。そして、元の数値が $10^m$ と表現でき、シフトする量が $10^n$ なのだとすれば、これを掛けた値は $10^m \times 10^n = 10^{m+n}$ となる。そこで、 $10^{m+n} \ge 10^K$ 、すなわち $m + n \ge K$ となったら、オーバーフローという感じ。（なお、math.frexpでは2進数での桁数が取れるので、それを10進数の桁数に直す計算が最初に入ってる）

逆に、「42.0を2³⁰桁で四捨五入したら0.0になること」というのも。これをさっきのmyround3()でやってみると、こうなる：

myround3(42.0, -(2**30))  # => ZeroDivisionError: float division by zero

これは逆に浮動小数点数で表現できる一番小さい幅よりも小さい値になってしまい（アンダーフロー）、値が0に落ちてしまってゼロ除算を発生させてしまったのが原因。これも同様にアンダーフローが発生するか（(num_exp + ndigits) < 0か）を確認し、発生していたら、それは極端に大きな桁で四捨五入していることになるので、0を返せばいい。

さらに、さきほどのsys.float_info.digに関連して、これより大きな桁については信頼できる数字ではなくなってしまうため、次のようなことが起こる：

for i in range(1, 32):
    print(f"[{i:02}] {myround3(0.5, i)}")
# 以下の出力になる
# [01] 0.5
# ...（省略）
# [15] 0.5
# [16] 0.5000000000000001
# [17] 0.5
# ...（省略）
# [31] 0.5

実行したマシンではsys.float_info.digが15なんだけど、それを超えた16桁のときに誤差が入り込んでしまってる。

こうなるともう浮動小数点数では対処しきれなくて、Rubyでは有理数を使った計算に切り替えて対処している。 Pythonの場合は10進数を10進数のまま扱えるようにするdecimalを使うといい。

パッケージへ

とまぁそんな感じで、たかが四捨五入、されど四捨五入。ちゃんとやろうとすると、かなり大変なことが分かる。

なので、パッケージを作って単体テストもRubyから移植し、マトモで堅牢な四捨五入を使えるようにした。それがsaneroundパッケージ。

四捨五入で四苦八苦しないためにも、ぜひ使ってほしい。

今日はここまで！

いものやま。

雑多な知識の寄せ集め

Pythonの四捨五入で四苦八苦した話。（その3）

Rubyでの四捨五入のアルゴリズム

さらに厳しいケースへ

パッケージへ