67100


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

[code]

    #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));
    }
    
[/code]

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

[code]

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

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


CC BY-SA 4.0
2020-05-28T20:22:48.190
67098
2020-05-28T23:35:31.817
2020-05-28T23:00:42.553
19110
37961
question
c 浮動小数点数
longとdoubleでsqrtの戻り値が異なる


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

[code]

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

を実行すると

[code]

    986074550526975.0000000000
    986074550526976.0000000000
    
[/code]

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

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

* * *

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

[code]

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

は

[code]

    9013925449473025
    9013925449473024.0000000000
    
[/code]

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

[code]

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

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

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

[code]

    #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;
    }
    
[/code]


CC BY-SA 4.0
2020-05-28T22:58:50.593
67099
2020-05-28T23:34:42.550
2020-05-28T23:34:42.550
19110
19110
67098
answer


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

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

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

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

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

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

[code]

    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;
    }
    
[/code]

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

[code]

    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
    
[/code]

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

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


CC BY-SA 4.0
2020-05-28T23:02:52.450
67100
2020-05-28T23:35:31.817
2020-05-28T23:35:31.817
13972
13972
67098
answer