POSIX::floor() を使うと計算が合わない

Perl で大きな数字を扱うと、計算結果があわないという相談を受けた。

- Perl の丸め誤差?

349347958500 * 466593284000 / 699889926000 という計算が 232898639000 にならずに 232898638999 になるという。

そもそもすごい桁数。3000億 * 4000億って、これ何の計算なの? と聞くと、金額計算のプログラムで、巨大数値入力テストでのエラーケースだという。なるほど、納得。

- ソースコードを確認

ちゃんとソースを持ってきてくれていたので、該当箇所を探してみる。あった。って、floor って関数を呼んでるのね。Perl というかこれが悪いんじゃないの? これってどの名前空間のメソッド? と思って先頭行の use を探してみるとそれらしきものがあった。

use POSIX qw(math_h)

POSIX って、POSIX 規格のモジュール? となると、これって POSIX の仕様?

- 実行環境を確認

テスト環境のターミナルを使わせてくれたので環境を確認。

Welcome to Linux 2.2.16.

valeria:~> uname -a
Linux valeria 2.2.16 #2 Mon Aug 14 23:19:26 JST 2000 i686 unknown

相談してきた人は管理者ではないので、ディストリビューションは不明。TurboLinux なのかな? それにしてもなかなかに古いなあ。まあ、特定環境で使うものだし、仕方がないでしょう。

valeria:~> perl -v

This is perl, version 5.005_02 built for i486-linux

Copyright 1987-1998, Larry Wall

Japanization patch 4 by Yasushi Saito, 1996

Modified by Hirofumi Watanabe, 1996-1998
jperl5.005_02-981225
EUC version

Perl も古い。5.005_03 じゃなくて 5.005_02 か。おや? Hirofumi Watanabe って Ruby の ML などでみかける わたなべ ひろふみ さん?

- コマンドラインから試す

Perl って比較的簡単にコマンドラインから試せるから好き。

valeria:~> perl -MPOSIX -le 'print floor(349347958500 * 466593284000 / 699889926000)'
232898638999

あー、確かに 232898639000 よりも少なくなるね。

valeria:~> perl -MPOSIX -le 'print 349347958500 * 466593284000 / 699889926000'
232898639000

floor() を使わなければ OK と。

valeria:~> perl -le 'print floor(349347958500 * 466593284000 / 699889926000)'
Undefined subroutine &main::floor called at -e line 1.

もちろん、標準の名前空間に floor() が無いことも確認。

調べてみると、POSIX::floor は C のライブラリを呼んでるだけということがわかった。POSIX 規格の関数を実装したライブラリだというなら、そうだろうなあ。結局そこの仕様を調べないと原因も対処も確定しない。

相談してきた人によると、どうするかは今後検討するとのこと。相談者はもともと別件でこのプログラムを修正していた。テストケースを増やしたところ、このエラーを見つけることができたとのこと。テスト重要。

- 追記

上記のメモをご覧になった方から、解説と Perl および C のサンプルコードをメールで頂いた。ありがとうございます。その結果理解したことをメモ。

1. 今回の計算は「丸め誤差」ではなく、オーバーフローが原因。
2. 今回の計算は、IEEE754 に則った倍精度浮動小数点数で表現できる範囲を
超えている。
3. Perl は倍精度浮動小数点数で計算している。
4. Perl であれば Math::BigFloat を使うとより高い精度で計算ができる。

この分野って過去に情報処理技術者試験のために学んだくらいで、普段全く使わない。1の補数とか2の補数とか、少ないビットで大きい数を表すとかね。

追記をしていても、自分の中でも消化し切れてない感じがある。「オーバーフローしている」のに、なぜ「正しい計算結果が出ている (ように見える)」のかがよくわからない。たまたまそういう結果になるような数値を使って計算したから? 「下二桁が00だから、その部分を除いて考えると倍精度浮動小数点数の表現範囲に収まる」と思いかけたけど、明らかに足りないし。少しずつ勉強していこう。

- 追記2

以下でも解説を頂いた。ありがとうございます。
http://www.kt.rim.or.jp/~kbk/zakkicho/zakkicho11.html#D20060 ...

これにしたがってprintの整形ルーチンが四捨五入して数値を丸めているので一見合っているかのように見えるだけ。それに対してfloorで丸めると切り捨てになるので、桁あふれが起きて生じた小数部分の分、答えが違うということになる。

メールで頂いた解説と上記 URL の説明を読み、実際に自分で試してみてやっと理解することができた。
オーバーフローの結果失われた数字が小さかったことと、四捨五入がオーバーフローして失われた分を偶然回復するような動きとなったために、結果として正しい計算ができたように見えたということだ。