PHPのround関数の謎が少し解けた

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:というか、スタックに積まれて