double-doubleでマシンイプシロンは0となる。

昨夜山下さんからdouble-doubleのマシンイプシロンはどういう数なのかということを問われた。なんと0である。従ってマシンイプシロンの用途としては不向きだ、ということである。以下少し説明する。

浮動小数点フォーマットにおけるマシンイプシロンeを
1+e = 1となる数とする。だいたいIEEE 754 doubleでは1.1e-16となる。
求めるときには、IEEE 754 doubleでは、基数が2なので、
#include <stdio.h>
int main(void){
double one = 1.;
double eps = 1.;
double one_plus_eps;
while (1) {
one_plus_eps = one + eps;
printf ("eps: %.16e 1+eps: %.16e\n", eps, one_plus_eps);
if (one == one_plus_eps) {
break;
}
eps = eps / 2.0;
}
return 0;
}
とすればよい。しかしいわゆるdouble-double(いわゆる倍々精度 cf. Hida and Baily's QD)では、これを定義するのが難しい。普通に考えると、2.2e-32程度ではなかろうかと予想して、以下のプログラムを走らせると、

#include <stdio.h>
#include <qd/dd_real.h>

int main(void){
dd_real one = 1.;
dd_real eps = 1.;
dd_real one_plus_eps;
while (1) {
one_plus_eps = one + eps;
printf ("eps: %.16e 1+eps: %.16e\n", eps.x[0], one_plus_eps.x[0]);
if (one == one_plus_eps) {
break;
}
eps = eps / 2.0;
}
return 0;
}
..
eps: 7.9050503334599447e-323 1+eps: 1.0000000000000000e+00
eps: 3.9525251667299724e-323 1+eps: 1.0000000000000000e+00
eps: 1.9762625833649862e-323 1+eps: 1.0000000000000000e+00
eps: 9.8813129168249309e-324 1+eps: 1.0000000000000000e+00
eps: 4.9406564584124654e-324 1+eps: 1.0000000000000000e+00
eps: 0.0000000000000000e+00 1+eps: 1.0000000000000000e+00
...となる。なんと"0"と、全然違った、しかも全く役に立たない値になる。

これは、double-doubleなので、二つのIEEE 754 doubleをもって、一つの実数を表すため、epsはIEEE 754の(denormalizedな)最小の数である、4.9e-324をもってしても、1+eps=1とはならないからである。足し算において切り捨てが行われないからである。つまりdouble-doubleでは、いかなるゼロより大きな値を持ってしても、1+eps=1となるのが存在しないのである。