3

一般的な仮定は、1 / x * x == 1. 一般的な IEEE 754 準拠のハードウェアでこれを破る最小の正の整数は?

乗法逆数の仮定が失敗すると、下手に書かれた有理数演算は機能しなくなります。C や C++ を含む多くの言語はデフォルトで浮動小数点数をゼロへの丸めを使用して整数に変換するため、小さなエラーでも整数の結果が 1 ずれることがあります。

簡単なテスト プログラムでは、さまざまな結果が得られます。

#include <iostream>

int main () {
    {
        double n;
        for ( n = 2; 1 / n * n == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
        for ( ; (int) ( 1 / n * n ) == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
    }
    {
        float n;
        for ( n = 2; 1 / n * n == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
        for ( ; (int) ( 1 / n * n ) == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
    }
}

GCC 4.3.4 を使用した ideone.com での結果は次のとおりです

41 (5.42101e-20)
45 (5.42101e-20)
41 (5.42101e-20)
45 (5.42101e-20)

GCC 4.5.1 を使用しても同じ結果が得られますが、エラー マージンは正確にゼロであると報告されています。

私のマシン (GCC 4.7.2 または Clang 4.1) では、結果は次のようになります。

49 (1.11022e-16)
49 (1.11022e-16)
41 (5.96046e-08)
41 (5.96046e-08)

これは、--fast-mathオプションに関係なくです。使うと-mfpmath=387驚くほど

41 (5.42101e-20)
41 (5.42101e-20)
41 (5.42101e-20)
41 (5.42101e-20)

値 5×10 -20は、64 ビット仮数に対応するイプシロン、つまり Intel 80 ビット拡張精度を使用した内部計算を暗示しているようです。

これは、FPU ハードウェアに大きく依存しているようです。テストに適した信頼できる値はありますか?

注: 言語標準やコンパイラが浮動小数点数システムについて何を保証しているかは気にしませんが、一般的なプログラミング システムには多くの意味のある保証があるとは思いません。数値と実世界のコンピューターとの相互作用について疑問に思っています。

4

2 に答える 2

4

倍精度の場合:

1/41 = 0x1.8F9C18F9C18FAP-6、および41*0x1.8F9C18F9C18FAP-6 = 0x1.000000000000028、1/45 = 0x1.6C16C16C16C17p-6、45*0x1.6.16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16C16 、これは 1 に丸められます。

でも、

1/49 = 0x1.4e5e0a72f0539p-6、および 49*0x1.4e5e0a72f0539p-6 = 0x0.fffffffffffffa4、これは 0x0.ffffffffffffff8 = 0x1.ffffffffffffff0p-1 に丸められます

ただし、49 には逆数があります。0x1.4e5e0a72f053ap-6 です。

より一般的には、f が [1, 2) の浮動小数点数の場合、f は逆数になります。通常の偶数への丸め演算では、数値が [1 - 2 -54 , 1 + 2 -53 ] にある場合、数値は 1 に丸められます。1/f に最も近い double、たとえば d は、1/f から 2 -54未満離れていることに注意してください。d > 1/f の場合、ゴールデンです。1 < f*d < f*(1/f+2 -54 ) <= 1 + 2 -54 * f < 1 + 2 -53なので、f*d は 1 に丸められます。d < 1/f の場合、f *d は 1 - 2 -53に丸められる場合があります。そうであれば、f*d は [1 - 2 -53 , 1 - 2 -54 ) にあります。e = 2 -53 + dを取る場合、e*f > 1 および e*f = d*f + 2 -53 *f < 1 - 2 -53+ 2 -52 = 1 + 2 -53、これも 1 に丸められます。

EDIT : 2 つの連続する double 間のストライドが 2 倍ずれているため、上記の推論は間違っています。逆数を持たない double の例は 0x1.ffffffbfffffe です。0x1.0000002000001p-1 は小さすぎますが、0x1.0000002000002p-1 は大きすぎます。逆数を持たない整数の最小の例は 237 です。この数値は小さすぎますが、その次の double は大きすぎます。

于 2012-11-26T07:41:53.497 に答える
0

この問題は、C++ の整数への変換方法の選択に関連しているようです。

比較のための Ada バージョンを次に示します。32 ビット、64 ビット、および 80 ビットの浮動小数点数をテストします (7、15、および 18 桁を要求するか、最初の 2 桁に組み込み型を使用してください)。

最初に結果とメモ、以下のコード。

$ gnatmake fp_torture.adb
gcc -c fp_torture.adb
gnatbind -x fp_torture.ali
gnatlink fp_torture.ali
$ ./fp_torture
 41 ( 5.96046E-08)
Error representing float  2.14748E+09 as integer
 49 ( 1.11022302462516E-16)
 2147483647 ( 0.00000000000000E+00)
 41 ( 5.42101086242752217E-20)
 2147483647 ( 0.00000000000000000E+00)
$

ご覧のとおり、浮動小数点計算は C++ の障害点を再現し、387 個の 80 ビット浮動小数点数の使用を確認します。ただし、(1 に非常に近い数値) を整数に戻すと、比較は機能します。

これを見て、C++ の例に適切な丸めを追加すると、比較が機能するようになります。MAX_INT に終了条件を追加すると、「double n」が機能します。

n のインクリメントに失敗すると、"float n" にポイントが来る++nので、反復子は反復を停止しますが、それは別の問題です!

以下の Ada バージョンはジェネリックを作成するので、任意の float 型でインスタンス化できます。(2^31 - 1 が 32 ビット float に変換され、バック オーバーフローが発生するため、例外ハンドラーが必要です...)

with Ada.Text_IO;   
use Ada.Text_IO;

procedure FP_Torture is

    generic
       type Float_Type is digits <>;
    procedure Test_FP;

    procedure Test_FP is
       F : Float_Type;
    begin
       -- for ( n = 2; 1 / n * n == 1; ++ n ) ;
       for i in 2 .. Natural'Last loop
          F := Float_Type(i);
          exit when 1.0 / F * F /= 1.0;
       end loop;
       Put_Line(natural'image(natural(F)) & " (" 
               & Float_Type'image(1.0 - (1.0 / F * F)) & ")");

       -- for ( ; (int) ( 1 / n * n ) == 1; ++ n ) ;
       for i in 1 .. Natural'Last  loop
          F := Float_Type(i);
          exit when natural(1.0 / F * F) /= 1;
       end loop;
       Put_Line(Natural'image(Natural(F)) & " (" 
               & Float_Type'image(1.0 - (1.0 / F * F)) & ")");
    exception
       when Constraint_Error => 
           Put_Line("Error representing float " & Float_Type'image(F) 
                    & " as integer");
    end;

    type Big_Float is digits 18;

    procedure Test7 is new Test_FP(Float);
    procedure Test15 is new Test_FP(Long_Float);
    procedure Test18 is new Test_FP(Big_Float);

begin
    Test7;
    Test15;
    Test18;
end FP_Torture;
于 2012-11-26T16:24:36.847 に答える