浮動小数点の罠?

PHP FreeBSD tech

以前のエントリで FreeBSD の PHP では log(64) の結果 float(6) を整数にキャストするとなぜか int(5) になるということを書きました。今回はその続き。
まず PHP の log 関数は ext/standard/math.c に実装されており、引数が二つの場合は

RETURN_DOUBLE(log(Z_DVAL_PP(num)) / log(Z_DVAL_PP(base)));

となっています。\log_m{n} = \frac{\log_o{n}}{\log_o{m}} の公式通りですね。


これを基に C でいろいろ試したところ、どうやら x86 系プロセッサの浮動小数点演算の癖 (バグか仕様かは知らない) が原因のような印象を受けました。
単に FPU だけでなく OS やコンパイラも含めての問題だとは思いますが。
また、検証の過程で Mac OS X の math.h には log2(double) があるけど FreeBSD には無いことを知ったり、任意の底で対数をとるとき、直接演算するのと、いったん分子と分母を変数に格納したときとで結果が異なる場合があることが分かりました。


前回のまとめとしては、知らず知らずのうちに地雷を踏んでしまう可能性があるので、ループ回数の判定に浮動小数点数は一切絡ませないほうが良いということになります。
もしかして常識なんでしょうか？


以下、テスト用プログラム (メモあり) をまるごとコピペ。

/*
Mac OS X:
cc -Wall -DHAVE_LOG2 log2.c -o log2
FreeBSD:
cc -Wall log2.c -lm -o log2

php/ext/standard/math.c の log(num, base) の実装:
RETURN_DOUBLE(log(Z_DVAL_PP(num)) / log(Z_DVAL_PP(base)));

zsh で簡易テスト
x=1; repeat 10; do; ./log2 $((2**$x)); x=$(($x+1)); echo; done

FreeBSD では my_log と my_log_v2 の結果が異なることがあった
真数が 2**(3**(2**n)), 2**(13**(2**n)) のときに現象を確認

Intel Mac では log2 と my_log(_v2) の結果が異なることがあった
真数が 2**(3**(2**n)), 2**(7**(2**n)), 2**(13**(2**n)) のときに現象を確認

結果が異なるとき、期待されない側の値は仮数部の最小の0でないビット以下が
反転して近似値になっている (e.g. 1.1 (1.5) -> 1.0111.. (1.4999..))
このため、int にキャストするとき小数点以下が切り捨てられる模様

ちなみに PowerPC Mac では問題なし -> x86 の浮動小数点演算器が原因?
*/
#include <stdio.h>
#include <string.h>
#include <math.h>

#define DEFAULT_NUMBER 64

double my_log(double, double);
double my_log_v2(double, double);
void dump_double(double);

int main(const int argc, const char **argv)
{
	int ret;
	double num, loga;
	const double base = 2;

	ret = 0;
	num = DEFAULT_NUMBER;

	if (argc > 1) {
		double param;
		if (sscanf(argv[1], "%lf", &param) && param > 0) {
			num = param;
		} else {
			ret = 1;
		}
	}

#ifdef HAVE_LOG2
	loga = log2(num);
	printf("log2(%lf) = %lf (", num, loga);
	dump_double(loga);
	printf(") -> %d\n", (int)loga);
#endif

	loga = my_log(num, base);
	printf("my_log(%lf, %lf) = %lf (", num, base, loga);
	dump_double(loga);
	printf(") -> %d\n", (int)loga);

	loga = my_log_v2(num, base);
	printf("my_log_v2(%lf, %lf) = %lf (", num, base, loga);
	dump_double(loga);
	printf(") -> %d\n", (int)loga);

	return ret;
}

double my_log(double num, double base)
{
	return log(num) / log(base);
}

double my_log_v2(double num, double base)
{
	double n, b;

	n = log(num);
	b = log(base);

	return n / b;
}

void dump_double(double d)
{
	unsigned char buf[50];
	size_t i, len;

	len = sizeof(double);
	memcpy(&(buf[0]), &d, len);
	for (i = 0; i < len; i++) {
		printf("%02x", buf[i]);
	}
}