doubleの==比較は思ったとおりにならない時ある

浮動小数点の比較の際の注意(訂正版)
※ 本補足資料の最初の版で、浮動小数点の絶対値を計算する関数を abs としていま
したが、fabs の誤りでした。
C 言語で浮動小数点を使った数同士の比較を行うとき、整数の比較と同じ感覚で次のよ
うに書いたプログラムがよく見られます。

double D;
// D に関する面積や導関数などの計算
if(D == 0){ // なにか処理 }

このように浮動小数点同士の比較を==で行うことはあまりお勧めできません。演習が
提供している C 言語の資料にもこのようなプログラム例があります*1が、これは解説の便
宜のためであって本当はあまりよくありません。

C 言語で浮動小数点同士の比較を行うときは、==を使わずに次のようにするべきです。

#include <math.h>
// ヘッダファイル math.h から絶対値を求める関数 fabs() の定義を読み込む
double D;
double epsilon = 1e-10;
// 代わりに、十分小さな数を以下のように定数として定義しておいてもよい。
// #define EPSILON 1e-10
// D に関する面積や導関数などの計算
if(fabs(D) < epsilon ){ // なにか処理 }

つまり、2つの浮動小数点 x と y が等しいかどうか判定するには、それらの差の絶対値x − y が十分小さいかどうかで判断します。(1e-10 は 1.0 × 10−10 のことです。C 言
語では MeN と書いて、M × 10N を表現します。)
計算機内部では浮動小数点は 64bit(0 または 1 の 64 桁の列) で表現しますが、それらは
あくまでも近似値に過ぎません。例えば、円周率 π = 3.141592 · · · は無限に続く小数点
のどこかから先を切り捨てたものになります。
浮動小数点の計算を繰り返すと誤差が発生し、多くの場合小数点の末尾付近の桁は本当
の数値とは異なるものとなります。したがって、数学的には一致するはずの2つの数値が
コンピュータ上の計算では完全には一致しないということが頻繁に発生します。このよう
な理由で、64bit が全て一致、すなわち小数点の並びが完全に一致することが要求される
== による比較は浮動小数点同士の比較方法としては不適切です。*2
浮動小数点の計算によって発生する誤差には他にも桁落ちというよく知られた問題があ
り、C 言語の演習資料にも書かれています。1/3 などの 2 以外の因数を持つ数による割り
算、sqrt などの算術関数で見られます。
浮動小数点の 0 判定を行う必要がある典型的なパターンとして、その浮動小数点が割り
算の分母に来る場合があります。もちろん、本当に 0 で割ってしまうのは避けなければな
りませんが、正確に 0 でなくてもそれに非常に近い数で割るのも好ましくありません。な
ぜなら、極端な桁数の変化が起こってしまうからです。
倍精度浮動小数点は 15 桁程度の精度を持ちますが、これは 1×100 から見ると 1×10−15
以下は実質的に 0 だということです。倍精度浮動小数点のアンダーフロー限界(表現でき
る最少の絶対値)は 10−300 くらいですが、これが分母にくると今度は 10300 くらいの値
となり、その値を使った計算にはもはや信頼性がありません。オーバーフローの危険性も
あります。
閾値の epsilon をどれくらいにするのがよいかというのは計算内容にも影響されるの
で一概には言えません。今回の演習では、倍精度実数の精度が 15 桁くらいだということ
と、(特に何も指定しない場合) printf が 7 桁くらいまでしか表示しないということから
ϵ = 10−10 くらいの設定でよいでしょう。
なお、この辺りについてもっと詳しく知りたい人は、演習室に配備してある参考書
「ニューメリカルレシピ・イン・シー — C 言語による数値計算のレシピ」の 1.3 節および
5.5 節あたりを読んでみてください。

ネタ元