Hatena::ブログ(Diary)

hnwの日記 このページをアンテナに追加 RSSフィード

[プロフィール]
 | 

2007年6月3日(日) PHPのround関数の謎が少し解けた このエントリーを含むブックマーク このエントリーのブックマークコメント

2週間以上前の記事「PHPの奇妙なround関数」がすごいことになっていますね。最近書き始めたばかりの日記にこんなに人が来るなんて、有名人の集客力は流石だなあ、などと感心しています。


その集客力のおかげかもしれませんが、FreeBSDMac OS Xだと挙動が違うよ、というコメントを頂きました。実際にFreeBSDで試してみたところ、確かにLinuxと異なる、いわばマトモな挙動です。その原因がわかりました、というのが本稿の概要です。僕がモタモタ記事を書いている間に理由がわかっちゃった人も居るかとは思いますし、より詳細なところまで把握した人も居そうですが、僕なりに現時点でわかったことを書いてみます。


前回の記事で、PHP_ROUND_FUZZという定数が「少なくとも僕の手元の環境では」0.50000000001と定義されている、と書きました。この詳細を説明すると、configureスクリプトで下記のコードを実行して、exit codeが1なら0.50000000001に、0なら0.5に定義しています。

#include <math.h>
/* keep this out-of-line to prevent use of gcc inline floor() */
double somefn(double n) {
  return floor(n*pow(10,2) + 0.5);
}
int main() {
  return somefn(0.045)/10.0 != 0.5;
}

僕の手元では0.50000000001になっている、というだけでも面白い事実だと思ったのですが、configure次第だということを書かなかったのはフェアでは無かったかもしれません。これでミスリードされた人も居たかと思います。ごめんなさい。


一方で、僕はこのコードの実行結果は浮動小数点数の四則演算とfloor(3)の挙動にのみ依存していて、少なくともx86系のCPUであれば僕の環境と同じ挙動を示すはずだと考えていました。このコードの挙動がFreeBSDLinuxで違うなどというのは全く想像しませんでしたが、実際のところx86FreeBSD上でconfigureするとPHP_ROUND_FUZZは0.5になります。


どちらもCPUx86系なのに実行結果が違う理由は、gccコンパイル結果が違うためです。Linuxgcc -O2で上のコードをコンパイルすると、上記ソースコード中のコメントとは裏腹に、n*pow(10,2)+0.5の計算結果はレジスタに乗ったままでインライン展開されたfloorに評価されてしまいます。つまり、80bitで評価されてしまいます。その結果、精度が十分過ぎるために0.045*100+0.5が5未満となってしまい、configureでのテストに失敗してしまいます。一方FreeBSDコンパイルすると期待通りfloor(3)が呼ばれますので、一旦計算結果はメモリに記録されて*1、64bitに丸められてから評価されます。結果、0.045*100+0.5が5.0となり、configureでのテストは成功します。


補足しておくと、x86系の浮動小数点数レジスタは80bitあります。doubleの64bitの値として受け取った値を、途中経過では80bitの精度で計算して、必要なときに64bitに丸めているわけです。gcc最適化オプションに-ffloat-storeというのがありますが、これをつけてコンパイルするとインライン展開されたfloorに入る前にfstplしてからfldlで取り出すようになります。無駄なメモリアクセスをして精度が落ちる、つまり丸めが起こるのでLinuxでもFreeBSDと同じ結果が得られるようになります。

`-ffloat-store'
     Do not store floating point variables in registers, and inhibit
     other options that might change whether a floating point value is
     taken from a register or memory.

     This option prevents undesirable excess precision on machines such
     as the 68000 where the floating registers (of the 68881) keep more
     precision than a `double' is supposed to have.  Similarly for the
     x86 architecture.  For most programs, the excess precision does
     only good, but a few programs rely on the precise definition of
     IEEE floating point.  Use `-ffloat-store' for such programs, after
     modifying them to store all pertinent intermediate computations
     into variables.

で、誰のせいかという話に立ち戻ると、configureでのバグのように思います。前回は仕様かと思いましたが、少なくともLinuxにおいては開発陣の意図と異なる状態だろうと思います。上のコード中のコメントからして、インライン版のfloorが使われてしまうと本来やりたいテストとは異なってしまうのではないでしょうか。Linux上でも-ffloat-storeをつけてコンパイルするか、最適化を-O0にしてしまえばマトモに生まれ変わりますけど、細かいバグを妙な方法で回避しているだけで本質的ではないように思います。


PHP中の人が何を考えて上記configureおよび例の0.50000000001の定義をしているかは僕の妄想なので、あちこちで事実のように引用されると困ってしまうのですが、少なくとも前回の妄想は間違いであるように思えてきました。「浮動小数点数の丸めの精度が怪しい環境なんて相手にしていられないから、round関数もアバウトにしておけばいいんじゃね?」ということなのかもしれません。で、Linuxが怪しいと勘違いされちゃったという状況なんですかね。これも妄想なので、真実を知っている人は本当に教えてほしいです。


最後に、両者のアセンブリ出力を貼っておきます。まずLinuxでの-O2の結果から。somefn:からretまでを見れば十分だと思いますが、念のため全部貼っておきます。

        .long   1120403456
        .align 4
.LC1:
        .long   1056964608
        .text
        .p2align 2,,3
.globl somefn
        .type   somefn, @function
somefn:
        pushl   %ebp
        movl    %esp, %ebp
        subl    $16, %esp
        flds    .LC0
        fmull   8(%ebp)
        fadds   .LC1
#APP
        fnstcw -2(%ebp)
#NO_APP
        movw    -2(%ebp), %ax
        andb    $243, %ah
        orb     $4, %ah
        movw    %ax, -4(%ebp)
#APP
        fldcw -4(%ebp)
        frndint
        fldcw -2(%ebp)
#NO_APP
        fstpl   -16(%ebp)
        fldl    -16(%ebp)
        leave
        ret
        .size   somefn, .-somefn
        .section        .rodata.cst4
        .align 4
.LC4:
        .long   1092616192
        .align 4
.LC5:
        .long   1056964608
        .text
        .p2align 2,,3
.globl main
        .type   main, @function
main:
        pushl   %ebp
        movl    %esp, %ebp
        subl    $8, %esp
        andl    $-16, %esp
        subl    $16, %esp
        pushl   $1067911741
        pushl   $1889785610
        call    somefn
        fdivs   .LC4
        flds    .LC5
        fxch    %st(1)
        popl    %eax
        fucompp
        fnstsw  %ax
        andb    $69, %ah
        xorb    $64, %ah
        setne   %al
        popl    %edx
        movzbl  %al, %eax
        leave
        ret
        .size   main, .-main
        .section        .note.GNU-stack,"",@progbits
        .ident  "GCC: (GNU) 3.4.6 20060404 (Red Hat 3.4.6-3)"

次がFreeBSDの-O2です。

        .file   "round-test.c"
        .section        .rodata.cst4,"aM",@progbits,4
        .p2align 2
.LC0:
        .long   1120403456
        .p2align 2
.LC1:
        .long   1056964608
        .text
        .p2align 2,,3
.globl somefn
        .type   somefn, @function
somefn:
        pushl   %ebp
        movl    %esp, %ebp
        flds    .LC0
        fmull   8(%ebp)
        fadds   .LC1
        fstpl   8(%ebp)
        leave
        jmp     floor
        .size   somefn, .-somefn
        .section        .rodata.cst4
        .p2align 2
.LC4:
        .long   1092616192
        .p2align 2
.LC5:
        .long   1056964608
        .text
        .p2align 2,,3
.globl main
        .type   main, @function
main:
        pushl   %ebp
        movl    %esp, %ebp
        subl    $8, %esp
        andl    $-16, %esp
        subl    $24, %esp
        pushl   $1067911741
        pushl   $1889785610
        call    somefn
        fdivs   .LC4
        flds    .LC5
        fxch    %st(1)
        fucompp
        fnstsw  %ax
        andb    $69, %ah
        addl    $16, %esp
        xorb    $64, %ah
        setne   %al
        movzbl  %al, %eax
        leave
        ret
        .size   main, .-main
        .ident  "GCC: (GNU) 3.4.4 [FreeBSD] 20050518"

追記:書き忘れていましたが僕の手元の環境のLinux2.6+gcc 4.1.0でもインライン版のfloorが使われます。

*1:というか、スタックに積まれて

hnwhnw 2007/06/03 18:45 前の記事のトラックバック元を読んでみると、Linuxでもマトモな環境があるみたいですね。configureのオプションが違うのか、コンパイラの挙動が違うのか、そんな違いがあったりするんでしょうか。

hnwhnw 2007/06/03 20:03 さらに追記ですが、configure中のコメントにつられて「インライン展開」と書いてしまいましたが、別にlibmから展開してるわけじゃなくてgccが独自に実装しているコードでよすね。きっとプロの方は気になるんだろうと思いますが、今更なので放置しておきます。

hnwhnw 2007/06/03 20:43 その後また進展がありました。この話題はまだ続きそうです。期待せずにお待ちください。どれくらいの人が続きを読んでくれるのかはさておき。

tsekinetsekine 2007/06/05 13:58 百歩譲って 0.5 + alpha を使うにしても、(1 + DBL_EPSILON) / 2 にすべきだと思うけどねー。

まぁ、どう考えても百歩譲れないけど。

hnwhnw 2007/06/05 16:34 ども。定量的な意見が出たのは初めてのような気がしますし、その方が理系的には生産的な方向ですよね。大感謝です。ただ、それだと具体的にはround(0.5-DBL_EPSILON/2)とround(0.5-DBL_EPSILON/4)とround(1.5-DBL_EPSILON)の挙動がおそらく変わりますが、2以上になると打ち切り誤差があるので、0.5のときと同じ結果になるかと思います。このアプローチを採用するときにalphaをいくつにすればいいかは案外自明では無いかと思います。

hnwhnw 2007/06/05 17:33 いちいち書くまでもないですけど、round(0.5-DBL_EPSILON/4*3)もそうですね。また、自明でないと言いましたが、来た引数に応じて「1個となり」までの距離の半分をalphaにする、という考え方もあるかと思います。ただ、この関数は第二引数の値によってfloorの前に乗算が入ることがあり、これらの打ち切り誤差などを考えると毎回「1個となり」の数まで含めるだけでいいのか、「2個となり」も含めるのか、というのはそれなりに難しいような気がします。また、引数が大きすぎると「1個となり」までの距離がかなりあったりして、「1個となり」を必ず含めていいのかという問題もあるかと思います。

tsekinetsekine 2007/06/05 21:12 あんまりよく考えてませんでした(汗

scale された後の値 x に対して、
alpha = min(max(floor(x), 1) * DBL_EPSILON / 2, magic_number)
みたいな感じかですかね?時間があったらもっとちゃんと考えてみます。ただ、こんなことに時間を割くのは無駄な気もしますが…

sodasoda 2007/06/06 19:41 Linux でもマトモな環境…というのは、gcc のバージョンの違いにより、生成されるコードが違っているんでしょうね。
sumii さんのところ (http://d.hatena.ne.jp/sumii/20070604/p1) にもありますが、これは gcc の違いというよりは、むしろ OS の違いです。
FreeBSD (に限らず全ての *BSD 系や Solaris/x86 など) では、たとえ例に挙げられている Linux の gcc -O2 の場合と同様に、FPU レジスタ上だけで演算を行なったとしても、80bit ではなく 64bit で計算が行なわれます。
Linux だけは、x86 FPU の control word レジスタに、拡張精度 (80bit) で計算するような初期設定がされているのが、この違いの原因です。
結局、この PHP のケースは、Linux のこのデフォルト動作に対する対処方法が間違っているということではないでしょうか? むしろ、Linux の場合でも、常に 64bit で計算を行なうようにするのが良いように思います。というのは、FPU レジスタで閉じた計算であるかどうかによって結果が異なるというのは、このケースに限らず、いろいろびっくりする結果を生むことが多いからです。
初期化時に次のような処理を行なえば、Linux の場合も他の OS と同様になります。(このコードは Linux/x86 依存です。念のため…)
#include <fpu_control.h>
#define FPU_PRECISION_MASK (_FPU_EXTENDED|_FPU_DOUBLE|_FPU_SINGLE) /*XXX*/
fpu_control_t npxcw;
_FPU_GETCW(npxcw);
npxcw = (npxcw & ~FPU_PRECISION_MASK) | _FPU_DOUBLE;
_FPU_SETCW(npxcw);
configure スクリプトから抜きだしたテストコードの main の先頭に上記のコードを加えると、CentOS 4.4 上で -O2 を指定した場合も、FreeBSD と同じ結果になることを確認しました。

sodasoda 2007/06/06 19:49 > これは gcc の違いというよりは、むしろ OS の違いです。

ちょっと意味不明な書き方でした。
ここの「これ」というのは、「configure から呼び出されるテストコードの結果の違い」を指しています。
「メモリにストアしたか否かで結果が異なる」というデフォルト動作を選択しているのはgcc ではなく OS (Linux/i386) 側ですから。

sodasoda 2007/06/06 20:21 さらに補足です。

> これは gcc の違いというよりは、むしろ OS の違いです。

これはもちろん、Linuxのバージョンによって、FPUの control wordの初期値が違っているという意味ではありません。
知る限り、この初期値は、x86 FPUのデフォルト設定そのままであり、Linuxの誕生以来変わってないと思います。

hnwhnw 2007/06/06 22:40 ありがとうございます。勉強になります。実はsumiiさんの記事のsodaさんのコメントを今日の昼頃に読んでいたところでした。

私はFPUの精度の指定を行えるということ自体を知りませんでしたので、教えていただいたような対応は全く思いつきもしないものでした。変数をvolatile指定するような対応はどうかな、などと考えていたのですが、教えて頂いた方法の方が確実で本質的ですね。

一方で、私が今日書いたエントリで話題にしているのですが、Pythonでもおそらく同様の現象が報告されているものの、「それはそれであるべき姿だ」という反応のようです(Pythonプロジェクト内のどういう立場の人なのか、私は知らないのですが)。やりとりの後半、下記のURLですね。

http://mail.python.org/pipermail/python-list/2005-September/340741.html

そう言い切れる知識に私は心から感心しましたが、一方で挙動を合わせられるなら合わせる方が現実的な対応だろうとも考えたりして、何が望ましい対応なのか混乱していたところでした。


よろしければもう一点だけご意見を頂きたいのですが、この_FPU_SETCWを使うなどしてx86系のFPUの精度をIEEE754に合わせる(という理解で正しいですよね)手法が何故PythonでもPHPでも採用されていないのか、sodaさんはどうお考えでしょうか。この知識があまり知られていないせいでしょうか。それとも、これを敢えて使わない理由が考えられるのでしょうか。よろしくお願いします。

sodasoda 2007/06/07 03:19 うーん、私は必要がない限り、浮動小数点演算には極力近付かないようにしているので良くわかりません。^^;
この件を知っているのは、たまたま各 OS の FPU の control word の初期値を比べる機会があったからに過ぎません。
というわけで、以下は全くの素人の推測です。といっても、おそらく hnw さんが既に考えているのと同じ内容のような気がしますが…
Python の方の答は原則論通りですよね。double の演算に拡張精度を使うことによって、びっくりする答がでることがあるわけですが、逆に桁落ちを避けられることもあるでしょうし、また拡張精度を使うことによって入ることがある誤差は、浮動小数点演算では当然予期すべき誤差の範囲なわけですから、Linux のデフォルトが必ずしも間違いであるとは言えません。ですから、Python 側では何もしないというのは一つの見識だと思います。ググってみたところ、Python のチュートリアル (http://www.python.jp/doc/2.4/tut/node16.html) からリンクされている The Perils of Floating Point (http://www.lahey.com/float.htm) でも、Too Much Precision の項にこれと同様な原因で起きる問題が載っているようです。
PHPの方は、control wordの設定で、他のOSと動作を合わせられるということを知らなかったのではないかという気がしてます。わざわざ configure で調べているということは、この挙動がシステム依存だということまでは気づいているからでしょうし。その違いを吸収したいのであれば、違っている部分(FPUレジスタ上での演算精度)を合わせるのが筋ではないでしょうか。これが前述のようなコメントを書いた理由です。
わざわざ演算精度を下げるのもどうかという議論もあるかもしれません。しかし拡張精度を使わないことによって、今回の例や The Perils of Floating Point にある例のように かえって誤差が減ることもあるわけですし、個人的には、浮動小数点演算の実装の詳細に慣れてない人にとって自然な動作は、むしろ拡張精度を使わない方ではないかと思います。詳しい人のために、拡張精度を使うためのスイッチを用意しておくのは良いかもしれませんが、PHPレベルの変数はほぼ必ずメモリにストアされるでしょうから、そこまでする必要は、おそらくないと思います。
あと、「FPUの精度をIEEE754に合わせる」という表現については、実際にIEEE754にあたったわけではないので推測ですが、たぶん誤りだと思います。IEEE754自体は拡張精度を明示的に許していますし、歴史的には、IEEE754はx86系のFPUの仕様を元に定められたという経緯があったと記憶していますので。doubleの演算に対して拡張精度を用いることをIEEE754が本当に許しているかどうかは確認していませんが…

hnwhnw 2007/06/09 14:53 ふむふむなるほど…っていやいや、どんな素人ですか!それはさておき、コメントありがとうございます。お返事が遅くなってしまいましたが、IEEE754に関して自分なりの理解に辿りつくのに時間がかかっていました。そんなレベルでこんな記事を書いていいのかどうか悩んでしまいますが、もう少し頑張ってみようと思います。

さて、IEEE754がdoubleの演算に対して拡張精度を許しているかどうかですが、sodaさんのおっしゃる通りで、どうやら許しているようですね。単精度や倍精度の演算結果が拡張精度のレジスタ(など)に格納される場合には拡張精度で丸めてもよいが、ユーザー(=高級言語のコンパイラなど)の指定によって倍精度や単精度に丸められるような仕組みも用意しなさい、というのが規格として決まっていることのようです。あくまで僕の理解ですけど。

残りについては僕が調べて考えるべきところを全部言われてしまった感じですね。結局はソフトウェアを開発する側が浮動小数点数の演算結果についてのポータビリティを重視するかどうかの問題で、Pythonは精度を落とすのを嫌ったけれども、PHPはおそらくポータビリティを取った(けど失敗した)ということでしょうね。

おかげさまで、自分の中でようやくまとまりました。本当にありがとうございました。心から感謝いたします。

 | 
ページビュー
1950996