特定の数値にしたときのみ結果がnanになる

Question

前提・実現したいこと

c言語で数値解析を行っています。プログラムは完成していて、実際に動かすことができているのですが、一部の結果だけnanになってしまいます。

発生している問題・エラーメッセージ

分割数を変化させると結果がnanになる

該当のソースコード

/*
    数値積分（台形面積近似）
    sqrt(4-x^2)の0から2までの面積（1/4円の面積）を求める
*/
int main()
{
    //分割数・刻み幅
    double m = 550000000;                         //分割数
    
    double h = (1.0-0.0)/(double)m;     //0-2までをm分割したときの1分割分の幅
    //xの値と和（面積計算）の値
    double x,sum;
    //ループ変数
    int i;

    //初期化
    sum = h;        //h/2*(f(0)+f(2))
    x = h;          //f(0+h)から開始

    //面積計算
    for(i=0;i<(m-1);i++){
        //面積和の計算
        sum += h*(sqrt(1.0-x*x));
        //途中経過出力
//      printf("%4d %f %f %f\n",i,x,(4.0-(x+h)*(x+h)),sum);
        //xを刻み幅分ステップアップ
        x += h;
    }

    //結果
    printf("%.15lf\n",sum*4);
}

試したこと

分割数mの変更
m=1000000000,m=100000000の時は実行結果が表示されましたが、m=550000000の時にnanになりました。

補足情報（FW/ツールのバージョンなど）

あまり詳しくないのでわかる範囲で
PC:mac
実行:vim version 8.1.2292

Comments

[Milly]

テキストエディタ(vim) は内容に関係していないため、vim タグは外した方がよいかと思います。
プログラムの話題では代わりにコンパイラ等のバージョンを書くとよいでしょう。(「gcc 7.5.0」等)

Answer 1

誤差の積み重ねによりxが1を超えます。
すると1-x*xが負になるのでsqrt(1-x*x)がnanになります。
xにhを足し込んでいくのではなく、iから求めるようにしてみてください

Comments

[退会済みユーザー]

@質問者様
これだけだと難しいかもしれないと思ったのでちょっとdoubleの中身を表示するコード書きました。
中身が2進数の小数を表しているので、10進数の小数は代入しただけで誤差が発生します。
今回は10進数の小数を代入している部分はありませんが、hに値を代入するところで誤差が発生します。
このhの加算時に誤差が累積していくのですが、今回のように桁上りを繰り返すような計算の場合、その誤差が単純な掛け算にならず、後半になるほど誤差が増えてしまいます。
毎回乗算していけば誤差も単純な掛け算になるのでマシになりますよという話です。

加算: x += h;
乗算: x = h * (i + 2);
ご参考までに。

#include <math.h>
#include <stdio.h>

union double_union {
double d;
unsigned long long u;
};

void print_binary(unsigned long long value, int nod) {
int i;
for (i = nod - 1; i >= 0; --i) {
int digit = 1 << i;
printf("%1d", ((value & (1ull << i)) != 0) ? 1 : 0);
}
}

/**
* @brief
* 倍精度実数の表示
*
* https://en.wikipedia.org/wiki/Double-precision_floating-point_format
*
* @param x 表示したい数
*/
void print_double(double x) {
union double_union du;
du.d = x;
int sign = (int)(du.u >> (64 - 1));
int exponent = (int)((du.u >> (64 - 1 - 11)) & ((1ull << 11) - 1));
unsigned long long fraction = du.u & ((1ull << 52) - 1);
printf("x(double): %e\n", du.d);
printf("x(hex): 0x%016llx\n", du.u);
printf("sign: %01x\n", sign);
printf("exponent: %d\n", exponent);
printf("fraction(hex): 0x%013llx\n", fraction);
printf("fraction(bin): ");
print_binary(fraction, 52);
printf("\n");
printf("value: %c%1d.", sign ? '-' : ' ', (exponent == 0) ? 0 : 1);
print_binary(fraction, 52);
printf(" x 2^(%d)\n", exponent - 1023 + (exponent == 0));
}

/*
数値積分（台形面積近似）
sqrt(4-x^2)の0から2までの面積（1/4円の面積）を求める
*/
int main()
{
//分割数・刻み幅
double m = 550000000; //分割数

double h = (1.0-0.0)/(double)m; //0-2までをm分割したときの1分割分の幅

// ちょっと表示
print_double(h);

//xの値と和（面積計算）の値
double x,sum;
//ループ変数
int i;

//初期化
sum = h; //h/2*(f(0)+f(2))
x = h; //f(0+h)から開始

//面積計算
for(i=0;i<(m-1);i++){
//面積和の計算
sum += h*(sqrt(1.0-x*x));
//途中経過出力
// printf("%4d %f %f %f\n",i,x,(4.0-(x+h)*(x+h)),sum);

// ちょっと表示
if ((i+1) == (m-1)) {
print_double(x);
}

//xを刻み幅分ステップアップ
x += h;
}

//結果
printf("%.15lf\n",sum*4);
}

Answer 2

まずは、doubleの表現範囲を調べてみれば。
そのうえで、その数値での計算結果を当てはめていけば、その数値のdouble上の表現がどうなるかわかるかと思います