2014年4月20日日曜日

消費税8%問題とceil関数の罠

1円が合わないのはなぜ? 


最近、消費税が8%に上がり、お買い物が楽しくなった今日この頃。チョット気になる問題が僕を襲いました。なんとceil関数の罠にドップリはまり、2、3日、酒浸りになりました。

ことの発端が、この計算
<?php
$v = ceil(15000 * 1.08);
print($v . "\n");
?>
cuomo@localhost ~ $ php sample.php 
16201
ちょ、まて、俺が望んでるのは、16200だぞ!

これに気づきませんでした...

「1円!違ってるけどどうなってんのっ!、いっぱい1円違うじゃないのぉーーー (怒)」

言われるはずです、早速直して対応したのですが、なんで?って思いますよね...
調べましょう、Gentooで......


そもそもの間違い


ceilのせいでは....ないんですよ、ceilはドキュメントに書いてあるとおりに動作しています、手で計算した結果を鵜呑みにして、さらに和をかけて思いこみの沙汰が酷過ぎ、単なる浅はかさがこの事態を生んだ模様です 笑...。

調査開始


まずPHPのceil関数をしらべてみたら、ext/standard/math.cでceil関数を呼んでいる、ceil関数はgentooですとglibcパッケージ(libm)に入っていますので、そちらを確認してみましょう。

<pre>
PHP_FUNCTION(ceil)
{
    zval **value;

    if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "Z", &value) == FAILURE) {
        return;
    }
    convert_scalar_to_number_ex(value);

    if (Z_TYPE_PP(value) == IS_DOUBLE) {
        RETURN_DOUBLE(ceil(Z_DVAL_PP(value))); /* ceil関数を呼び出し */
    } else if (Z_TYPE_PP(value) == IS_LONG) {
        convert_to_double_ex(value);
        RETURN_DOUBLE(Z_DVAL_PP(value));
    }
    RETURN_FALSE;
}
</pre>


ですね、これを元にceilを利用したサンプル関数を作ってみる

実際はどんなことしてるんだ


glibcのコードをソースコードレベルで解析したいので、面倒ですがglibcをdebug版でコンパイルしなおす。まぁ、emergeするだけなので楽勝ですが、出きるまでビールでも呑んでましょう。
サンプルコードとgdbを利用して、libmのceil関数へブレークポイントを仕掛けて、実行してみる。


#include <stdio.h>
#include <math.h>

double subceil(int);

int main(void)
{
  double v;
  v = subceil(15000);
  printf("%f\n", v);
  return 0;
}

double subceil(int v)
{
  return ceil(v * 1.08);
}
cuomo@localhost ~ $  gcc -g ceil.c -L/usr/lib64/debug/lib64 -lm -lc -O0 -o c.out
これを、gdbで止めてみると../sysdeps/x86_64/fpu/multiarch/s_ceil.S:37で停止するのですがSSE4のroundsd命令で 一撃でやられちゃってよくわからないので(自分がSSEをよくわかってない...笑)、正確ではないのですが違う関数ceillで再度確認。
Breakpoint 1, __ceil_sse41 () at ../sysdeps/x86_64/fpu/multiarch/s_ceil.S:37
37              roundsd $2, %xmm0, %xmm0
(gdb) l
32      END(__ceil)
33      weak_alias (__ceil, ceil)
34
35
36      ENTRY(__ceil_sse41)
37              roundsd $2, %xmm0, %xmm0 /* xmm0レジスタ一撃ですわ */
38              ret
39      END(__ceil_sse41)

自分的にわかりにくかったので、ceillでもう一度
#include <stdio.h>
#include <math.h>

double subceil(int);

int main(void)
{
  double v;
  v = subceil(15000);
  printf("%f\n", v);
  return 0;
}

double subceil(int v)
{
  return ceill(v * 1.08); /* ceil()をceill()に修正 */
}

もう一度、ビルドし直し、gdbでブレークしてみる

cuomo@localhost ~/Code/binaries/CProgs/ceil $ gdb ./c.out
GNU gdb (Gentoo 7.6.2 p1) 7.6.2
Copyright (C) 2013 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-pc-linux-gnu".
For bug reporting instructions, please see:
<http://bugs.gentoo.org/>...
Reading symbols from /home/cuomo/Code/binaries/CProgs/ceil/c.out...done.
(gdb) b ceill
Breakpoint 1 at 0x400520
(gdb) r
Starting program: /home/cuomo/Code/binaries/CProgs/ceil/./c.out
warning: no loadable sections found in added symbol-file system-supplied DSO at 0x7ffff7ffa000
warning: Could not load shared library symbols for linux-vdso.so.1.
Do you need "set solib-search-path" or "set sysroot"?

Breakpoint 1, ceill () at ../sysdeps/x86_64/fpu/s_ceill.S:12
7
8       #include <machine/asm.h>
9
10
11      ENTRY(__ceill)
12              fldt    8(%rsp)                 /*Floating Point Load */ 
13
14              fstcw   -4(%rsp)                /* store fpu control word */
15
16              /* We use here %edx although only the low 1 bits are defined.
17                 But none of the operations should care and they are faster
18                 than the 16 bit operations.  */
19              movl    $0x0800,%edx            /* コントロールレジスタのRCフィールドを切り上げにセット */
20              orl     -4(%rsp),%edx
21              andl    $0xfbff,%edx
22              movl    %edx,-8(%rsp)
23              fldcw   -8(%rsp)                /* コントロールレジスタを設定 */
24
25              frndint                         /* 切上処理の実行 */
26
27              fldcw   -4(%rsp)                /* コントロールレジスタを元に戻す */
28
29              ret
30      END (__ceill)
31      weak_alias (__ceill, ceill)
(gdb)
多分こちらはx87の浮動小数点計算処理でSSE命令より余計な部分が多いのですが、何をやっているのかわかりやすいのでこちらで確認。 x87の細かい使い方は専門の方へお願いするとして、数バイトstep実行させた後のST(0)の値、ここへ15000*1.08の計算結果が格納されている。この時点でお分かりですよね、そう小数点の部分が0ではないのです、これをfrndintでやってしまえばそれはそうなりますよって話 笑...
(gdb) info all-reg
...
st0            16200.000000000001818989403545856476 (raw 0x400cfd20000000000800)
...

ここからはつまらない話


x86拡張倍精度で80bitの形式で保存されていて、16進数で0x400cfd20000000000800という値を保存している、で、さらにこれを2進数80bitで表現すると...
0 100000000001100 1 111110100100000000000000000000000000000000000000000100000000000


  • 符号部(s) 0(1bit)
  • 指数部(e) 100000000001100(15bit バイアス値 16383)
  • 整数部(p) 1(1bit)
  • 仮数部(f)  111110100100000000000000000000000000000000000000000100000000000(64bit)
でわけられ、下の式で計算される
(-1)^s * 2^(e-16383) * p.f
※fは2進数を小数に変換したもの
これに習って計算してみると
8192 * 1.9775390625 ~= 16200.00000000000181... ※1.9775390625 ~= 1 + (1/2 + 1/4 + 1/8 + 1/16 + 1/32 + 1/128 + 1/1024 + 1/4503599627370496) こういうことに内部的になっているため、frndintを実行した場合に1を加えた数値になります。
実際の計算を終えた後の、浮動小数点レジスタの値はこの様な感じになってました。
st(0)の値がしっかり16201になってますね。
(gdb) info all-reg
...
...
es             0x0      0
fs             0x0      0
gs             0x0      0
st0            16201    (raw 0x400cfd24000000000000)
st1            0        (raw 0x00000000000000000000)
st2            0        (raw 0x00000000000000000000)
st3            0        (raw 0x00000000000000000000)
st4            0        (raw 0x00000000000000000000)
st5            0        (raw 0x00000000000000000000)
st6            0        (raw 0x00000000000000000000)
st7            0        (raw 0x00000000000000000000)
fctrl          0xb7f    2943
fstat          0x3a20   14880
ftag           0x3fff   16383
fiseg          0x7fff   32767
fioff          0xf7b5b610       -139086320
...
...

まとめると


通常の計算結果で端数が出ない浮動小数点の計算でも、機械の内部ではぴったり収まらないケースがあるということ。小数点は意外とムズカシイ... なので、解決方法として
  • 仕様としてなるべく切捨てを採用するよう上司に涙ながらに訴える
  • ceil(round($v, x))で小数点以下のどこで切り上げるか明確にする
  • 切上げの場合は、安易に掛け算とceilを利用しないで、足し算、引き算などで頑張る(やり方は分かりません)
  • Haskellのような分数と遅延評価を使える言語で頑張る(みんなにいやがられる)

ちなみにHaskell
Prelude> 15000 + (15000 * 8/100 :: Rational)
16200 % 1

みなさんceil関数をこんな感じで使ってる所ってない?..なんか気持ち悪いですね....ヒヤリハットです

最後に、消費税をもう今後上げない...無理でしょうけど

0 件のコメント:

コメントを投稿