SCMの浮動小数点書き出しアルゴリズムは、Steele&White ("How to print floating-point numbers accurately, SIGPLAN '90, pp.112-123) に示されているFP^3アルゴリズムと基本的に同じである。 すなわち、書き出したい数値v (v ＞ 0) に対して、
( m - ε) × 10^e ≦ v ≦ ( m + ε) × 10^e 1 ≦ m ＜ 10 ε は上の式の成立させる最も小さい正の浮動小数点数
となるm, ε, e を求め、
1. m の整数部を書き出す。
2. ε ≦ mの小数部 ≦ 1.0 - ε なら
m ← (mの小数部)×10.0
ε ← ε × 10.0
として、1.から繰り返す
3. そうでなければ、 mの小数部 ＜ 0.5 なら そのまま終了
mの小数部 ＞ 0.5 なら最後の桁をひとつ増やす
(m の小数部＝0.5の時は丸め規則に従う)。
4. e ≠ 0 なら指数部を書き出す。
という具合である。
SCMはこれをdouble型の変数を使って実装している (src/scl.c, idbl2str)。 問題は、double型変数を使った場合、mのスケーリングによる誤差によって 最後の桁の方が怪しくなることだ。
次の2つの数値を書き出す場合を考えてみよう。
9007199254740797 × 2^-53 ≒ 9.999999999999784e-1 9007199254740796 × 2^-53 ≒ 9.999999999999782e-1
最初の m を計算するためには、仮数部を10倍する必要があるのだが…
9007199254740797 × 10 ＝ #x13ffffffffff862 9007199254740796 × 10 ＝ #x13ffffffffff858
これらの数値は少なくとも54ビットないと表現できない。 ところが、double型 (IEEE-754倍精度浮動小数点数)では仮数部は 53ビット(hidden bit含む)しかない。したがってこの計算をdoubleで 行うと最後のビットが丸められて、両者ともに
#x13ffffffffff86
になってしまう。このまま計算を進めると、どちらの数値も
9.999999999999783e-1
を印字してしまう。
すなわち、倍精度浮動小数点演算のみを使って倍精度浮動小数点数を 「正確に」書き出すのは原理的に不可能なのだ。
ではSCMはどのようにしてこの問題を回避しているのだろうか。 実は、SCMの表示ルーチンにこの問題を回避するコードは入っていないのだ。 にもかかわらず、SCMでは両者がちゃんと区別して印字される。
一体どうなっているのか不思議に思って計算の途中結果を 表示させてみると、さらに謎は深まる。 どちらの数値を与えても、出発点となる最初にスケールされたmの値(SCMのコードではf)は等しい。 にもかかわらず、それ以降の計算には差が出て来る。
秘密は、x86系の浮動小数点演算ユニットにある。 x86系の浮動小数点演算ユニットのレジスタは80bitの拡張浮動小数点数 (仮数部64bit)で、floatでもdoubleでも演算は一度80bitに直してから行われ、 結果が単精度または倍精度に丸められる。 ところが、ある種の最適化が行われた場合、 途中結果が浮動小数点レジスタに置いたままにされるため、80bitのまま 演算が進行する。
その値をスタック経由で関数に渡す時にはdoubleに丸められるので、 printfに渡して印字させるだけでは違いがわからなかったのだ。
この問題の微妙さは、次のようなコードで確かめることができる。
#include <stdio.h>
#include <math.h>
static double mult10(double m, double ten)
{
return m * ten;
}
int main(int argc, char **argv)
{
double m = 9.999999999999784e-1;
/*double m = 9.999999999999782e-1;*/
int d, i;
printf("m = %.20lg\n", m);
m = mult10(m, 10.0);
for (i = 0; i < 3; i++) {
printf("m = %.20lg\n", m);
d = m;
m -= d;
m *= 10.0;
}
return 0;
}
Linux (RH7.2) + gcc2.96, glibc 2.2.4 で、-O2でコンパイルしたものだと こうなった：
shiro:toccata% ./a.out # m = 9.999999999999784e-1 m = 0.99999999999997835065 # 初期値 m = 9.9999999999997832845 # 最初のスケール m = 9.9999999999978328447 # 次のスケール m = 9.9999999999783284466
shiro:toccata% ./a.out # m = 9.999999999999782e-1 m = 0.99999999999997823963 # 初期値 m = 9.9999999999997832845 # 最初のスケール m = 9.9999999999978239629 # 次のスケール m = 9.9999999999782396287
最初のスケールでmの値が両者で等しくなっているのに、次のスケールでは 異なっているのに注目。
この現象はコンパイラの最適化に依存する。この例でも、-O2をつけなければ 違いは発生しないし、mult10()の呼び出しをただの乗算に置き換えても 発生しない。また、FreeBSDなどではプロセスの浮動小数点演算モードを 明示的に拡張浮動小数点数を使うように設定してやらねばならないようである。
関連: Gauche:二重丸めの落とし穴
Tag: 数値