4

<math.h>sqrt 関数について質問です。

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

int main(void) {
    long n = 100000000;
    long l = 94941695;
    double d = 94941695;
    printf("sqrt [long]  : %.10f\n", sqrt(n * n - l * l));
    printf("sqrt [double]: %.10f\n", sqrt(n * n - d * d));
}

を実行すると、次の結果が得られました:

$ ./a.out
sqrt [long]  : 31401823.9999999851
sqrt [double]: 31401824.0000000000

何故結果が異なるのでしょうか。
なお、10 ^ 16 - 94941695 ^ 2 = 986074550526975で、
  √(986074550526975) = 3.14018239999999840... × 10 ^ 7
でした。よろしくお願いします。

2 件の回答 2

3

演算精度のせいで、sqrt(n * n - l * l)に与えられる引数の値とsqrt(n * n - d * d)に与えられる引数の値が異なるためです。

(なお、longの表す桁数は処理系によって異なりますが、ここでは出力例から、longは64ビット整数型を表すものとします。)

n * n - l * lを引数とする場合、nおよびlは共にlong型であり、また、途中計算のn * nl * lの結果もlong型の範囲に収まるため、計算結果は正確にlong型の986074550526975となり、その値が暗黙の型変換でdoubleに変換された後にsqrtが呼ばれます。

それに対して、n * n - d * dを引数とする場合、n * nは一旦longで計算され、正確にlong型の10000000000000000となりますが、引き算の相手のd * ddouble型なので、その値はdouble型に変換されます。

ここで肝心なのは、double型の精度は10進換算で15〜16桁であり、1000000000000000094941695 * 94941695の計算結果として期待される9013925449473025もその限界を超えているため、doubleではそれらの値を正確には表せない のです。

言葉による説明だけは分かりにくいので、途中結果を見てみましょう。

int main(int argc, const char * argv[]) {
    long n = 100000000;
    long l = 94941695;
    double d = 94941695;
    printf("n * n [long] : %ld\n", n * n);
    printf("l * l [long] : %ld\n", l * l);
    printf("n * n - l * l: %.10f\n", (double)(n * n - l * l));
    printf("sqrt [long]  : %.10f\n", sqrt(n * n - l * l));
    printf("n * n[double]: %.10f\n", (double)(n * n));
    printf("d * d[double]: %.10f\n", d * d);
    printf("n * n - d * d: %.10f\n", (double)(n * n - d * d));
    printf("sqrt [double]: %.10f\n", sqrt(n * n - d * d));
    return 0;
}

結果(<-以下は注釈):

n * n [long] : 10000000000000000
l * l [long] : 9013925449473025
n * n - l * l: 986074550526975.0000000000
sqrt [long]  : 31401823.9999999851
n * n[double]: 10000000000000000.0000000000  <- `print`で表示した結果は正しく見えるが1の位まで正確には表現できない
d * d[double]: 9013925449473024.0000000000   <- `print`で表示した結果にも誤差が表れている
n * n - d * d: 986074550526976.0000000000    <- `986074550526975`じゃない!
sqrt [double]: 31401824.0000000000

(doubleは内部では2進数で表現されているので、10進数との変換の際に誤差が入ったり、逆に誤差が相殺されて見掛け上正しい値に見えてしまうことがあるので、10進表記で全てを理解しようとするのには無理があるんですが、ここでは見掛けにはっきり違いが表れているので詳細には立ち入らないことにします。)

と言うわけで、doubleで計算を行う場合、途中計算まで含めて「doubleの精度で正しく表せるかどうか」を気にしないと、整数で正確な計算をしたのと同じ結果にはならないことに注意しないといけません。

12
  • 細かいですが、10,000,000,000,000,000=‭23'86F2'6FC1'0000‬については下16bitが0なので必要精度40bitで、これ自体はdoubleで正確に表せます。
    – sayuri
    2020年5月28日 23:21
  • @sayuri さん、「正確に」と言う言葉の問題かもしれませんが、1000000000000000010000000000000001との違いを表現できないのに「正確に表せます」と言ってしまうのには問題があるように思われます。ただそれを言うなら「誤差を含む」と言う言い方も良くありませんね。そこは訂正します。
    – OOPer
    2020年5月28日 23:25
  • @OOPer さん、そのように正確さを定義すると、たとえば 0 と 0.00000000...0001 の違いを表現できないことを理由に「0 を正確に表せない」ことになってしまいませんか?
    – nekketsuuu
    2020年5月28日 23:32
  • @nekketsuuu さん、0.00000000...0001を問題にする用途なら「0 を正確に表せない」と言うのはその通りですよ。今は整数演算との誤差の話をしているので、1の位まで正確に計算結果に表せるかどうかが問題にされていると考えるべきです。
    – OOPer
    2020年5月28日 23:38
  • うーん、あまり納得できていません。おそらく「値を正確に表せない」という部分にひっかかっているのだと思います。値それ自体は正確に表せているものの、何かしら演算をすると(演算相手が整数でも)正確じゃなくなる場合があるという状況を的確に表現できていないと感じています。「double ではこの値と演算すると誤差が出る可能性があります」くらいなら分かります。
    – nekketsuuu
    2020年5月28日 23:48
1

sqrt ではなく、その前段の n * n - d * d の時点で丸め誤差が生まれています。

printf("%.10f\n", double(n * n - l * l));
printf("%.10f\n", n * n - d * d);

を実行すると

986074550526975.0000000000
986074550526976.0000000000

と出力され、下 1 桁が異なります。

IEEE 754 の倍精度では、2 進数で 53 桁まで情報を保持することができます(一番上の整数部分の桁を含む)。これは 10 進数だと約 15.95 桁で、今 17 桁 ひく 16 桁の計算をしようとしているのでどこかで丸めが起きる可能性があります。


では、ちゃんとビット表現まで確認してみましょう。実は d * d の時点で丸め誤差が起きています。簡単に言うと、

long n = 100000000;
long l = 94941695;
double d = 94941695;
printf("%ld\n", l * l);
printf("%.10f\n", d * d);

9013925449473025
9013925449473024.0000000000

と出力します。この部分をビット表現で確認すると、

// 上は long、下は double です。double は仮数部だけ見てください。
l * l =    0000000000100000000001100001111000010000011100101001110000000001
d * d = 0 10000110100 0000000000110000111100001000001110010100111000000000
//                   ↑ この桁から比較してください。
//        ↑ 指数部
//      ↑ 符号部

というように一番下の 1 が落ちてしまっていることが分かります。

このビット表現は下のプログラムで確かめました。

#include <stdio.h>

typedef union {
    double d;
    long l;
} value;

void print_binary(long n, int digit)
{
    if (digit <= 64) {
        print_binary(n >> 1, digit + 1);

        printf((n & 1) ? "1" : "0");
        if (digit == 53 || digit == 64) {
            printf(" ");
        }
    }
    return;
}

void print_binary_long(long l) {
    print_binary(l, 1);
    printf("\n");
    return;
}

void print_binary_double(double d) {
    value v;
    v.d = d;
    print_binary(v.l, 1);
    printf("\n");
    return;
}

int main(void) {
    long l = 94941695;
    double d = 94941695;

    printf("l2 =  ");
    print_binary_long(l * l);

    printf("d2 = ");
    print_binary_double(d * d);

    return 0;
}

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。