hnwの日記

PHPのround関数を読み解く (1)丸め桁数が大きすぎ・小さすぎる場合

PHPのround関数はPHP5.3.0で一新されましたが、その挙動は複雑です。以前の記事「PHP5.3.0alpha3のround関数の実装がPHP5.2.6と変わった」では典型的な挙動を紹介しましたが、実際はもう少し細かい場合分けがあります。僕自身の整理も兼ねて、RFC*1およびPHPのCソースコードと対応づけながら何回かに分けて紹介していきます。


まずは、第2引数(精度)が大きすぎまたは小さすぎる場合の処理を紹介します。例えば次のコードを考えてみましょう。

<?php
ini_set("precision",20);
var_dump(round(5.3e-25,24)); //float(9.999999999999999237E-25)
var_dump(1e-24); //float(9.999999999999999237E-25)


0.0000…ときて小数点以下第25位が5で第26位が3であるような数を小数点以下第24位で丸めるという例で、理想的な答えとしては1e-24になります。1e-24はIEEE754でピッタリ表現できる数ではありませんから、ピッタリ表現できる数の中で1e-24に最も近い数になれば正解だと言えます。実際、5.3.0以降のPHPでは上のように2つの結果が一致しますので、理想的な結果になっているとわかります。


ところで、round関数を素直に実装してしまうと、上のような結果にはなりません。x>0の場合、round(x)の素直な実装はceil(x+0.5)であり、精度指定がある場合は乗算や除算を行って整数丸めになるような桁数に調整して丸めることになります。


つまり、今回の例で言うと5.3e-25に1e24を掛けて0.5くらいの数を作り、これを丸めた数1.0を1e24で割る、というのが素直な実装による処理です。ところが、1e24は浮動小数点数で正確に表せないため、この最終段階の除算で誤差が生じてしまうことがあります。実際、素直な実装をしているPythonで同じ計算を行うと次のような結果になります。

#!/usr/bin/python
print "%.18e" % round(5.3e-25,24) # 1.000000000000000107e-24
print "%.18e" % 1e-24 # 9.999999999999999237e-25


このような問題に対処するため、PHPはround関数の第2引数が大きすぎたり小さすぎたりする場合には、この最終段での割り算の代わりに10進小数の文字列(この場合なら"1.0e-24")を作り、浮動小数点数文字列から最も近い浮動小数点数を得る内部関数zend_strtod*2を利用しています。


これがPythonのバグかというと非常に微妙だと思います。PHPの方が正しい値を返すのは確かですが、大抵のプログラマにとってどうでもいい場所に分岐を増やさないで欲しいなあ、というのが僕個人の感想です。


とはいえ、精度指定つきの丸めにはこんな罠があるんだ、というのは面白い知識だと思います。

RFC上の記述

この部分について、RFCの記述を引用しておきます。例として出している数5.3e-24は、おそらく5.3e-25の間違いでしょう。

Special handling for large places difference


This was taken from the 2004 proposal: If the numer of places are very large, then 10^places is very large, too, and cannot be represented in an exact manner anymore. This will cause inaccuracies with the final division, as explained earlier.


The solution for that is that the rounded double is converted to a string and e-places is added to that string. Take, for example, the number 5.3e-24 which you may want to round to 24 places precision (which it of course already is). After rounding, the float value is 1.0 and that has to be divided by 1e24. But 1e24 is too large to be exactly represented, so instead a string 1.0e-24 is generated and passed through strtod(). strtod() on the other hand will make sure that the nearest double representation for that number is chosen. This of course has a performance penalty, but anybody wanting to round such small numbers will probaby be willing to pay for it.


This change will make sure that very small (< 1e-22) or large (> 1e22) numbers will be rounded correctly to the given precision.

対応するCソースコード

以下はPHP5.3.6の実装です。第2引数の精度指定が23以上か-23以下の場合にelseの方の分岐を通り、計算結果を一度文字列にしているのがわかります。

        /* see if it makes sense to use simple division to round the value */
        if (abs(places) < 23) {
                if (places > 0) {
                        tmp_value = tmp_value / f1;
                } else {
                        tmp_value = tmp_value * f1;
                }
        } else {
                /* Simple division can't be used since that will cause wrong results.
                   Instead, the number is converted to a string and back again using
                   strtod(). strtod() will return the nearest possible FP value for
                   that string. */

                /* 40 Bytes should be more than enough for this format string. The
                   float won't be larger than 1e15 anyway. But just in case, use
                   snprintf() and make sure the buffer is zero-terminated */
                char buf[40];
                snprintf(buf, 39, "%15fe%d", tmp_value, -places);
                buf[39] = '\0';
                tmp_value = zend_strtod(buf, NULL);
                /* couldn't convert to string and back */
                if (!zend_finite(tmp_value) || zend_isnan(tmp_value)) {
                        tmp_value = value;
                }
        }


(8/2追記)ちなみに、なぜ1e23が境界線になっているかというと、1e22まではdoubleで誤差無く表現できるからです。これは不等式「5^22<2^53<5^23」で検算できます(2は何回掛けても指数部にしか影響を与えないので、10ではなく5の累乗で計算することになります)。

*1:https://wiki.php.net/rfc/rounding

*2:William D Clinger, "How to Read Floating Point Numbers Accurately", ACM SIGPLAN '90, pp.92--101, 1990.に基づいた実装