2週間以上前の記事「PHPの奇妙なround関数」がすごいことになっていますね。最近書き始めたばかりの日記にこんなに人が来るなんて、有名人の集客力は流石だなあ、などと感心しています。
その集客力のおかげかもしれませんが、FreeBSDとMac 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であれば僕の環境と同じ挙動を示すはずだと考えていました。このコードの挙動がFreeBSDとLinuxで違うなどというのは全く想像しませんでしたが、実際のところx86なFreeBSD上でconfigureするとPHP_ROUND_FUZZは0.5になります。
どちらもCPUはx86系なのに実行結果が違う理由は、gccのコンパイル結果が違うためです。Linuxでgcc -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:というか、スタックに積まれて
