gcc 丸めモードと最適化オプション

% gcc -v
...
gcc バージョン 4.2.4 (Ubuntu 4.2.4-1ubuntu4)
% gcc hoge.c -lm && ./a.out
0: 4.00099999999999980105e+01, 4.00100000000000051159e+01
1: 4.00099999999999980105e+01, 4.00100000000000051159e+01
% gcc -O1 hoge.c -lm && ./a.out
0: 4.00099999999999980105e+01, 4.00099999999999980105e+01
1: 4.00099999999999980105e+01, 4.00100000000000051159e+01
% gcc -O2 hoge.c -lm && ./a.out
0: 4.00099999999999980105e+01, 4.00099999999999980105e+01
1: 4.00099999999999980105e+01, 4.00100000000000051159e+01
% gcc -O3 hoge.c -lm && ./a.out
0: 4.00099999999999980105e+01, 4.00099999999999980105e+01
1: 4.00099999999999980105e+01, 4.00099999999999980105e+01

数値演算の結果がおかしなことになってしまった.
hoge.c でやっていることは 40 と 0.01 を上丸めと下丸めで足し合わせて 40.01 を含む区間を求める.

丸めモードの指定は fesetround を使用して実現しているが, 上丸めがうまくいかない.
ほかの値でやると下丸めがうまくいかないこともある.

add0 と add1 で値が変わるのも理解できない.

以下使用したソースコード (hoge.c)

#include <stdio.h>
#include <fenv.h>
#define UP      fesetround(FE_UPWARD)
#define DOWN    fesetround(FE_DOWNWARD)
#define NEAR    fesetround(FE_TONEAREST)


typedef struct {
        double sup;
        double inf;
} intv;


#define ADD(x,y,z)      DOWN; (z)->inf=(x)->inf+(y)->inf; UP; (z)->sup=(x)->sup+(y)->sup;


intv add0(intv x, intv y)
{
        intv z;
        ADD(&x,&y,&z);
        NEAR;
        return z;
}

intv add1(intv x, intv y)
{
        intv z;
        ADD(&x,&y,&z);
        if (z.inf <= z.sup) { NEAR;}
        return z;
}

int main()
{
        intv d1, d2, d3;
        d1.inf = d1.sup = 40;
        d2.inf = d2.sup = 0.01;
        d3 = add0(d1, d2);
        printf("0: %1.20e, %1.20e\n", d3.inf, d3.sup);
        d3 = add1(d1, d2);
        printf("1: %1.20e, %1.20e\n", d3.inf, d3.sup);
        return 0;
}