4

コードは次のとおりです。

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

static double const x = 665857;
static double const y = 470832;

int main(){
    double z = x*x*x*x -(y*y*y*y*4+y*y*4);
    printf("%f \n",z);
    return 0;
}

不思議なことに（私にとって）このコードは、GCC 4.6を使用して32ビットマシンで（または私の場合のように64ビットマシンでは-m32フラグを使用して）コンパイルすると「0.0」を出力します。浮動小数点演算について私が知っている限り、それらをオーバーフロー/アンダーフローしたり、精度を失ったりする可能性がありますが... 0？どのように？

前もって感謝します。

4

4 に答える 4

7

問題は、数値がオーバーフローすることではありません。問題は、double には、減算の 2 つのオペランドを区別するのに十分な精度がないことです。

の値x*x*x*xは 196573006004558194713601 です。

の値y*y*y*y*4+y*y*4は 196573006004558194713600 です。

これらの数値は 78 ビットで、最後のビットだけが異なります。倍精度数は 53 ビットしかありません。他の数値は 53 ビットのみに丸められます。

あなたの場合、2 つのオペランドは同じ数値に丸められるため、それらの差は 0 です。

z の式を少し書き直すと、さらに奇妙なことが起こります。

double z = x * x * x * x - ((y * y + 1) * y * y * 4);

この変更により、33554432 が得られます。なんで？中間結果が丸められる方法により、右側のオペランドの最後のビットが異なるためです。最後のビットの値は 2^(78-53)=2^25 です。

于 2012-05-08T19:56:41.550 に答える
6

任意精度の整数で式を評価する:

Prelude> 665857^4 - 4*(470832^4 + 470832^2)
1

通常、 adoubleの精度は 53 ビットのみで、中間結果の精度は 78 ビットであるため、結果を正確に計算するには精度が不十分であり、丸められ、ある時点で最後のビットが忘れられます。

于 2012-05-08T19:55:21.577 に答える
4

コードに浮動小数点のオーバーフローまたはアンダーフローはありません。2 つの量は 1.96573006 × 10^23 のオーダーで、ほぼ に収まりますdouble。あなたの例は、壊滅的なキャンセルを単に示しています.2つの近い量を差し引くと、結果の相対的な精度が恐ろしくなります。

http://en.wikipedia.org/wiki/Loss_of_significanceを参照

于 2012-05-08T19:55:19.187 に答える
3

これは、IEEE 754 が浮動小数点数を正規化された形式で表現する方法の結果です。float または double またはその他の IEEE 754 準拠の表現が次のように格納されます。

1.xxxxxxxxxxxxxxxxxxx * 2^exp

ここxxxxxxxxxxxxxxxxxxxで、 は仮数の小数部であるため、仮数自体は常に [1, 2) の範囲内にあります。常に 1 である整数部分は表現に格納されません。xビット数は精度を定義します。の場合は 52 ビットですdouble。指数はオフセット形式で格納されます (値を取得するには 1023 を減算する必要があります) が、現在は関係ありません。

64 ビット IEEE 754 の 665857^4 は次のとおりです。

0 10001001100 (1)0100110100000001111100111011101010000101110010100010
+ exponent    mantissa

(最初のビットは符号ビットです: 0 = 正、1 - 負; 括弧内のビットは実際には保存されません)

80 ビット x86 拡張精度では、次のようになります。

0 10001001100    (1)0100110100000001111100111011101010000101110010100010
0 100000001001100 1 010011010000000111110011101110101000010111001010000111000111011

(ここで、整数部分は明示的に表現の一部です。IEEE 754 からの逸脱です。わかりやすくするために仮数を揃えました)

64 ビット IEEE 754 および 80 ビット x86 拡張精度の 4*470832^4 は次のとおりです。

0 10001001100    (1)0100110100000001111100111011101001111111010101100111
0 100000001001100 1 010011010000000111110011101110100111111101010110011100100010000

64 ビット IEEE 754 および 80 ビット x86 拡張精度の 4*470832^2 は次のとおりです。

0 10000100110    (1)1001110011101010100101010100100000000000000000000000
0 100000000100110 1 100111001110101010010101010010000000000000000000000000000000000

最後の 2 つの数値を合計するときの手順は次のとおりです。小さい値の指数は大きい値の指数と一致するように調整され、仮数は値を維持するために右にシフトされます。2 つの指数は 38 異なるため、小さい方の仮数が 38 ビット右にシフトされます。

調整された 64 ビット IEEE 754 および 80 ビット x86 拡張精度の 470832^2*4:

 this bit came from 1.xxxx ------------------------------v
0 10001001100    (0)0000000000000000000000000000000000000110011100111010|1010
0 100000001001100 0 0000000000000000000000000000000000000110011100111010101001010101

これで、両方の量の指数が同じになり、仮数を合計できます。

0 10001001100 (1)0100110100000001111100111011101001111111010101100111|0010
0 10001001100 (0)0000000000000000000000000000000000000110011100111010|1010
--------------------------------------------------------------------------
0 10001001100 (1)0100110100000001111100111011101010000101110010100001|1100

内部で合計が 80 ビットのより高い精度で行われるため、バーの右側に 80 ビット精度のビットの一部を残しました。

次に、64 ビット + 80 ビット表現のいくつかのビットで減算を実行しましょう。

minuend    0 10001001100 (1)0100110100000001111100111011101010000101110010100001|1100
subtrahend 0 10001001100 (1)0100110100000001111100111011101010000101110010100001|1100
-------------------------------------------------------------------------------------
difference 0 10001001100 (0)0000000000000000000000000000000000000000000000000000|0000

純粋な0！完全な 80 ビットで計算を実行すると、再び純粋な 0 が得られます。

ここでの本当の問題は、指数が 2^77 の 64 ビット精度で 1.0 を表現できないことです。仮数には 77 ビットの精度がありません。これは 80 ビット精度にも当てはまります。仮数部には 63 ビットしかなく、指数が 2^77 の場合、1.0 を表すのに必要なビットよりも 14 ビット少なくなります。

それでおしまい！それは、数学の授業で教えられたように何も機能しない科学計算の素晴らしい世界です...

于 2012-05-08T21:44:46.130 に答える