数値計算で生じる小さなごみ

質問

ある問題を数値的に解こうと思いFortranでプログラムを組んでいるのですが、倍精度の数値計算で以下のような事が起きて困っています。
プログラムの中で
A=B-C*D (変数は全て倍精度実数変数)
のような代入文があるのですが、
write(*,*) B,C,D
とすると、0.250000000000000 0.500000000000000 0.500000000000000
と表示されるのですが、
write(*,*) B-C*D や write(*,*)A
では-1.79570984817912D-017 などと表示されます。(この数値に意味はないのでしょうが、同じシチュエーションでは常に同じ値が表示されます)
入力した値の精度外の「誤差」ですが、どこに原因のある問題でしょうか?

因に使っているコンパイラはintelのサイトで入手したフリーの評価版ifc ver.6 でそれにpentium4に最適化するオプションをつけてコンパイルしています

回答No.3

noname#11476 · ベストアンサー

他に間違っている部分がなくて、その通りとするとおかしいですね。
本当にB,C,Dに0.2,0.5の数字が入っていてそうなるのであればコンパイラのバグかもしれませんが、オーダーが大きすぎるし、pen4が実際には計算している筈なので余り考えられません。
一応、最適化オプションをはずして試しても同じですか？
あと、単精度ではどうですか？

0.25, 0.5 という数字は２進数では割り切れる数値なのでそのようなことが起きること自体が不思議です。

気になるのは一見0.25とか0.5に見えても、何らかの計算の結果の数値とすればそのような現象になるのは至極当然のことです。
というのも、表示するときには必ず丸められて表示されるので、0.2とか0.5に対してD-17だと完全に丸められる領域だからです。
差分をとったときにはもちろんクローズアップするのでわかります。

試しに、 0.5 - 0.5 が 0になることを確認して、 C - 0.5 を表示させるなどすれば、Dが本当に0.5なのかどうかがわかります。
もう一つは、一度単精度の変数に代入して、また倍精度の変数にいれるとか、あるいは関数を使って下１０桁で四捨五入処理するなどですね。

この手の話は良くあることです。

お礼

ありがとうございました。
>試しに、 0.5 - 0.5 が 0になることを確認して、 C-0.5
>を表示させるなどすれば、Dが本当に0.5なのかどうかがわかります。
教えていただいたことを実行してみたところ、0.25に関してゴミがそのままでてきました。
なるほどこういうことがあるんですね。勉強になりました。
>この手の話は良くあることです。
のようですね。うっかり周りにいた詳しい人にも尋ねた所、大笑いされてしまいました。
どうもお世話になりました。

回答No.4

noname#21649

一つの可能性は．
Ｂ＝０．２５
のように指定すると．Ｂは倍精度．０．２５は単精度と解釈して．０．２５００００（中略）１２３．．．なんて数値が代入されている可能性があります。
０．２５ｄ１
のような倍精度を指定して代入式を作ってみてください。
これが失敗した場合には．全部整数で代入式を作り演算によって倍精度にする方法があります。

インテルの場合．たしか．有効桁以下は０以外の数値が入っている場合が多いので．これが変な数値の原因となっていると考えられます。

お礼

ありがとうございます
>０．２５ｄ１
>のような倍精度を指定して代入式を作ってみてください。
ご指摘の通り代入式もいい加減でした。レベルがばれそうで恥ずかしいです。全部書き直しました。直接の原因ではなかったみたいですがこれから気をつけることにします。
>インテルの場合．たしか．有効桁以下は０以外の数値が入っている場合が多いので．これが変な数値の原因となっていると考えられます。
確かにそうでした。これからは気をつけたいです。

みなさんありがとうございました。この場をかりてお礼申し上げます。

回答No.2

aton

同じくFortranはよく知らないのですが…。
お使いになっているFortranではBCD (Binary Coded Decimal, 二進化十進数→参考URL) 型の実数は扱えないのでしょうか?
これが使えれば，誤差の問題は解決出来るかもしれません。

参考URL：

http://www.ffortune.net/comp/develop/data/kazu.htm

お礼

ありがとうございました。
はずかしいことに二進化十進数は初耳でした。でも確かにこうしないとだめですよね。私のもっているfortranの教科書には実数は単精度、倍精度などの分類はありますが内部での表現の指定方法はないです。無いのかもしれないです。
他のfortranのコンパイラでは問題の無いプログラムなんですが。

回答No.1

noname#202350

Fortranは分からないので恐らくなのですが・・・。
コンピュータで計算を行う場合、大体に置いて2進数で計算を行います。

ところが、小数点を含む計算を行う場合に問題が出てくるのです。
というのは、10進数での特定の数字が2進数では表せないのです。
（詳しい事は昔教えてもらったのですが忘れてしまいました。）

そのため、2進数に変換が出来ないので2進数の近似値で計算を行います。
これによって誤差が生じてくるのです。

一般的に通貨型とか言われているような型でしたら誤差が出ないと思います。
たとえばVisualBasicではCurrency型がそれにあたります。

お礼

御返事ありがとうございました。
なにも考えないで計算させてましたが中の事も考えないとだめですよね。
いろいろ考えてみましたが、２進数で循環小数になる場合などは確かに有限桁で切ると誤差の原因になりそうです。
あと、質問には書き忘れましたが、write(*,*)0.25-0.5*0.5は正確に0となります。変数に値が入っている時だけなのです。