昨夜はどうもありがとうございました．解決策が2つみつかりました．

>>668 などで指摘されている通り，double型(64ビット)の変数がレジスタ
(80ビット)に割り当てられた際に，余分な桁がまぎれこみ，計算結果に
違いが出たようです．

私が思いついたのは次の2つで，結局，2を選択しました．

1. g++のコンパイルの際に-ffloat-storeオプションをつける．
これにより，double型(64ビット)変数のレジスタ(80ビット)への割り
当てを抑制できました．

2. double型の変数をlong double型で宣言しなおす．
もともとlong double型は80ビットですので，レジスタに割り当てても
　 不要な桁が入りこまないようになります．

ちなみに，計算速度は，1，2で，大きな差はありませんでした．