longとdoubleでsqrtの戻り値が異なる


`<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  
でした。よろしくお願いします。




演算精度のせいで、`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`の精度で正しく表せるかどうか」を気にしないと、整数で正確な計算をしたのと同じ結果にはならないことに注意しないといけません。




演算精度のせいで、`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`の精度で正しく表せるかどうか」を気にしないと、整数で正確な計算をしたのと同じ結果にはならないことに注意しないといけません。

