浮動小数点演算について知っておくべきこと

遠い時代、IT業界にとっては、これは前世紀の70年代で、数学者（プログラマーは以前呼ばれていました）はドンキホーテのように、当時の小さな風車のサイズのコンピューターとの不平等な戦いで戦いました。 深刻なタスクが設定されました：軌道からの画像からの海洋の敵の潜水艦の検索、長距離ミサイルの弾道の計算など。 それらを解決するために、コンピューターは実数で動作する必要があります。実数は、ご存じのように連続体であり、メモリは有限です。 したがって、この連続体をゼロと1の有限セットにマッピングする必要があります。 速度、サイズ、精度の妥協点を求めて、科学者は浮動小数点数（またはブルジョアなら浮動小数点）を提案しました。



何らかの理由で、対応するデータ型がすべてのプログラミング言語に存在することを考えると、浮動小数点演算はコンピューターサイエンスのエキゾチックな分野と見なされます。 正直に言うと、CPUとGPUで同じ問題を解決しながら、私自身はコンピューター演算をそれほど重視していませんでしたが、結果は異なります。 この領域の隠されたコーナーには、非常に奇妙で奇妙な現象があります。非可換性および非連想性の算術演算、符号付きゼロ、不等号の差がゼロなどです。 この氷山の根は数学に深く入り込んでおり、切り口の下では表面にあるものだけを概説しようとします。

1.基本

整数のセットは無限ですが、特定の問題を解決するときに生じる整数を表すために、このようなビット数をいつでも取得できます。 実数のセットは無限であるだけでなく連続的でもあるため、いくらいくら取るにしても、正確な表現を持たない数字に遭遇することは避けられません。 浮動小数点数は、実数を表示する方法の1つです。これは、精度と許容値の範囲の間の妥協点です。



浮動小数点数は、通常、符号、 指数 、順序、および仮数に分割された個別の数字のセットで構成されます。 順序と仮数は、符号とともに、次の形式で浮動小数点数の表現を与える整数です。







数学的には、次のように書かれています。



（-1） ^s ×M×B ^E 、ここでsは符号、B-base、Eは次数、Mは仮数です。



基数は、数字の数体系を決定します。 基数B = 2の浮動小数点数（2進表現）が丸め誤差に最も耐性があることが数学的に証明されているため、実際には、基数2と、あまり一般的ではない10のみが発生します。浮動小数点は次のようになります。



（-1） ^s ×M×2 ^E



仮数と順序とは何ですか？ 仮数は、実数の最上位ビットを表す固定長の整数です。 仮数が3ビット（| M | = 3）で構成されているとします。 たとえば、バイナリシステムでは101 _2に等しい数値「5」を使用します。 最上位ビットは2 ² = 4に対応し、中央（これはゼロです）は2 ¹ = 2で、最下位ビットは2 ⁰ = 1です。 順序は、最高ランクのベース（2）の次数です。 この例では、E = 2です。 このような数値は、たとえば「1.01e + 2」など、いわゆる「科学的」標準形式で記述すると便利です。 仮数が3文字で構成され、順序が2であることはすぐにわかります。



仮数の同じ3ビットを使用して小数を取得するとします。 たとえば、E = 1の場合、これを行うことができます。 その後、私たちの数は等しくなります



1.01e + 1 = 1×2 ¹ + 0×2 ⁰ + 1×2 ^-1 = 2 + 0.5 = 2.5



ここでは、E = 1なので、最初の桁の2のべき乗（小数点の前に来る）は「1」です。 右側（小数点以下）にある他の2つの放電は、2 ^E-1と2 ^E-2 （それぞれ2 ⁰と2 ^-1 ）の寄与を提供します。 明らかに、Eを調整することにより、同じ数を異なる方法で表すことができます。 仮数の長さ| M | = 4の例を考えてみましょう。 番号「2」は次のように表すことができます。



2 = 10（バイナリ）= 1.000e + 1 = 0.100e + 2 = 0.010e + 3。 （E = 1、E = 2、E = 3）



同じ番号には複数の表現があることに注意してください。 これは、機器にとっては便利ではありません。 数値を比較するとき、およびそれらに対して算術演算を実行するときは、表現の多様性を考慮する必要があります。 さらに、これは、表現の数が有限であり、繰り返しはまったく表現できる数の数を減らすため、経済的ではありません。 したがって、すでに非常に最初のマシンで、彼らはトリックを使用し始め、仮数部の最初のビットを常に正にしました。 このようなプレゼンテーションは、 normalizedと呼ばれていました 。







暗黙的なユニットをメモリに保存する必要がないため、これにより1ビットが節約され、数値の一意の表現が提供されます。 この例では、「2」には単一の正規化表現（「1.000e + 1」）があり、仮数は「000」としてメモリに保存されます。 上位ユニットは暗黙的に暗示されます。 しかし、数値の正規化された表現では、新しい問題が発生します-この形式でゼロを表現することは不可能です。



厳密に言えば、正規化された数値の形式は次のとおりです。



（-1） ^s ×1.M×2 ^E.



問題解決の質は、浮動小数点数の表現の選択に大きく依存します。 このようなビューを標準化する問題にスムーズに取り組みました。

2.少しの歴史

60年代および70年代には、浮動小数点数、丸め方法、および算術演算を表す単一の標準がありませんでした。 その結果、プログラムは非常にポータブルではありませんでした。 しかし、さらに大きな問題は、さまざまなコンピューターに「奇妙な」要素があり、プログラムでそれらを認識して考慮する必要があることでした。 たとえば、2つの異なる数値の差はゼロを返しました。 その結果、式「X = Y」と「XY = 0」が競合しました。 職人は、問題を避けるために乗算と除算の操作の前に割り当て「X =（XX）+ X」を行うなど、非常にトリッキーなトリックでこの問題を回避しました。



浮動小数点数の表現のための単一の標準を作成するというイニシアチブは、1976年にIntelが8086およびi432の新しいコプロセッサ用の「最適な」算術演算を開発しようとする試みと一致しました。 この分野の科学クジラ、教授。 ジョン・パーマーとウィリアム・カーハン。 後者は彼のインタビューで、Intelが算術演算を開発していた真剣さにより、他の企業が団結して標準化プロセスを開始することを余儀なくしたと述べました。



誰もが真剣でした。なぜなら、あなたのアーキテクチャを促進し、それを標準にすることは非常に有益だからです。 彼らの提案は、モトローラのZilogのNational SuperconductorのDECによって提示されました。 メインフレームのメーカーであるCrayとIBMは、傍観者の目から見ていた。 Intelはもちろん、新しい算術演算も導入しました。 提案された仕様の作成者は、William Kahan、Jerome Kunen、Harold Stoneであり、彼らの提案はすぐに「KCS」と呼ばれました。



すぐに、DECのVAXとIntelの「KCS」の2つを除いて、すべてのオファーが拒否されました。 VAX仕様ははるかに単純で、PDP-11コンピューターに既に実装されており、そのパフォーマンスを最大限に引き出す方法が明確でした。 一方、KCSには、「特別な」数字や「非正規化された」数字など、多くの便利な機能が含まれていました（詳細は以下）。



「KCS」では、すべての算術アルゴリズムが厳密に定義されており、実装では結果がそれらに一致する必要があります。 これにより、この仕様のフレームワーク内で厳密な計算を表示できます。 数学者が数値的方法で問題を解決するために使用し、解の特性を証明した場合、これらの特性がプログラムに保存されるという保証はありませんでした。 KCS演算の厳密さにより、浮動小数点演算に基づいた定理を証明することができました。



DECは、仕様を標準にするためにあらゆることを行っています。 彼女は、KCS算術演算が原則としてDECと同じパフォーマンスを達成できなかったという評判の高い科学者の支持を求めました。 皮肉なことに、Intelは仕様を生産的にする方法を知っていましたが、これらのトリックは企業秘密でした。 Intelが秘密を認めず、秘密の一部を明らかにしなかった場合、DECの猛攻撃を封じ込めることはできなかったでしょう。



標準化の戦いの詳細については、 Kahan教授のインタビューを参照してください。浮動小数点数の表現がどのように見えるかを見ていきます。

3.今日の浮動小数点数のプレゼンテーション

KCS開発者が勝ち、現在、その子孫はIEEE754標準に組み込まれています。 その中の浮動小数点数は、次のように符号（s）、仮数（M）、および順序（E）の形式で表示されます。



（-1） ^s ×1.M×2 ^E



発言。 新しいIEE754-2008標準では、2を底とする数値に加えて、10を底とする数値、いわゆる10進浮動小数点数があります。



Wikipediaで見つけられる過剰な情報で読者を混乱させないために、単精度（浮動）の1つのデータ型のみを考慮します。 半精度、倍精度、拡張精度の数値は同じ機能を持ちますが、順序と仮数の範囲が異なります。 単精度数（浮動/単精度）では、順序は8ビットで構成され、仮数は23です。有効な順序はE-127として定義されます。 たとえば、0.15625という数値はメモリに次のように書き込まれます。





ウィキペディアからの図。



この例では：

• 符号s = 0（正の数）
• 注文E = 01111100 ₂ -127 ₁₀ = -3
• 仮数M = 1.01 ₂ （最初の単位は明示的ではありません）
• その結果、F = 1.01 ₂ e-3 = 2 ^-3 +2 ^-5 = 0.125 + 0.03125 = 0.15625

3.1特別な番号：ゼロ、無限、および不確実性

IEEE754では、数値「0」はE = E _min -1（単一の場合は-127）に等しい次数とゼロ仮数の値で表されます。 独立した数値としてゼロを導入すると（正規化された表現ではゼロを表現できないため）、算術の多くの奇妙な問題を回避できました。 また、ゼロの操作は個別に処理する必要がありますが、通常は通常の数値よりも高速に実行されます。



また、IEEE754は特別な番号の表現を提供し、その操作は例外を引き起こします。 これらの数値には、無限大（±∞）と不確実性（NaN）が含まれます。 これらの数値により、オーバーフロー時に適切な値を返すことができます。 無限大は、E = E _max +1の次数とゼロ仮数を持つ数値として表されます。 オーバーフローおよびゼロ以外の数値をゼロで除算することにより、無限大を取得できます。 開発者は、配当と除数が特定の数になる傾向がある場合、限界の存在に基づいて部門の無限大を決定しました。 したがって、c / 0 ==±∞（たとえば、3/0 = +∞、および-3 / 0 =-∞）。これは、被除数が定数になり、除数がゼロになる場合、制限が無限になるためです。 0/0では、制限は存在しないため、結果は不確実になります。



不確実性またはNaN（非数値から）は、算術演算が常に何らかの意味のない値を返すことができるように考案された表現です。 IEEE754では、NaNはE = E _max +1の数値として表され、仮数はゼロではありません。 NaNを使用した操作は、NaNを返します。 必要に応じて、プログラムが解釈できる仮数部に情報を書き込むことができます。 これは規格では指定されておらず、仮数はほとんどの場合無視されます。



NaNを入手するにはどうすればよいですか？ 次のいずれかの方法で：

• ∞+（-∞）
• 0×∞
• 0/0、∞/∞
• sqrt（x）、x <0

定義上、NaN≠NaNです。したがって、変数の値を確認するには、それを自分で比較するだけです。

なぜゼロ記号（または+0対-0）

好奇心reader盛な読者は、浮動小数点数の記述表現に符号のみが異なる2つのゼロがあることに気づいたと思います。 したがって、3・（+0）= + 0、および3・（-0）=-0です。 ただし、+ 0 = -0を比較する場合。 標準では、符号は意図的に保持されたため、乗算または除算の際にオーバーフローまたは有意性の喪失の結果として無限またはゼロになる式は、依然として最も正しい結果を提供できます。 たとえば、ゼロに符号がなかった場合、1 / =∞と1 /-∞は0であるため、x =±∞の場合、式1 /（1 / x）= xは真になりません。



別の例：

（+∞/ 0）+∞= +∞、一方で（+∞/ -0）+∞= NaN



この場合、無限大はNaNよりどのように優れていますか？ NaNが算術式に現れる場合、式全体の結果は常にNaNになるという事実。 式で無限が検出された場合、結果はゼロ、無限、または通常の浮動小数点数になります。 たとえば、1 /∞= 0。

3.3非正規化数

非正規非正規化（非正規）数とは、簡単な例を考えてみましょう。 仮数の長さ| M | = 2ビット（+ 1ビットの正規化）と-1≤E≤2のオーダーの値の範囲を持つ正規化表現を考えてみましょう。 この場合、16個の数字を取得します。







大きなストロークは、1.00の仮数を持つ数字を示します。 ゼロから最も近い数値（0-0.5）までの距離は、この数値から次の数値（0.5-0.625）までの距離よりも大きいことがわかります。 つまり、0.5から1までの任意の2つの数値の差は、これらの数値が等しくなくても0になります。 さらに悪いことに、1より大きい数の差は0.5と0の間のギャップに入ります。たとえば、「1.5-1.25 = 0」です（図を参照）。



すべてのプログラムが「ゼロに近い穴」に陥るわけではありません。 70年代の統計によると、平均して、各コンピューターは月に一度この問題に遭遇しました。 コンピューターが大衆化したことを考えると、KCS開発者はこの問題をハードウェアレベルで解決するのに十分なほど深刻だと考えました。 彼らが提案した解決策は次のとおりでした。 E = E _min -1（フロートの場合は「-127」）および仮数がゼロの場合、数値はゼロに等しいと見なされます。 仮数がゼロでない場合、数値は非ゼロとみなされ、その順序はE = E _minであると想定され、仮数の暗黙的な上位ビットはゼロに設定されます。 このような数値は非正規化と呼ばれます。



厳密に言えば、浮動小数点数は次のようになります。



（-1）E min≤E≤Emaxの場合、 ^s ×1.M×2 ^E （正規化数）



（-1） ^s × ^0。E = E _min -1の場合、M×2 ^Emin 。 （非正規化数）



例に戻りましょう。 E _min = -1 新しい順序値E = -2を導入します。この場合、数値は非正規化されます。 その結果、数値の新しい表現が得られます。







0から0.5の間隔は非正規化された数値で埋められ、上記の0の例（0.5-0.25および1.5-1.25）で失敗しないようにすることができます。 これにより、ゼロに近い数値の丸め誤差に対してビューがより堅牢になりました。



しかし、プロセッサで数値の非正規化表現を使用する贅沢は無料では提供されません。 そのような数はすべての算術演算で異なる方法で処理する必要があるという事実により、そのような算術での作業を効果的にすることは困難です。 これにより、プロセッサでのALUの実装にさらに困難が生じます。 非正規化された数値は非常に便利ですが、万能薬ではなく、ゼロへの丸めを追跡する必要があります。 したがって、この機能は標準の開発において障害となり、最も強い抵抗に会いました。

3.4 IEEE754の番号シーケンス

IEEE754形式で数値を表現する驚くべき機能の1つは、順序と仮数が、実行される整数{n}のシーケンスを形成するように次々に配置されることです。



n <n + 1⇒F（n）<F（n + 1）。ここで、F（n）は整数nからビットを順序と仮数に分割することにより形成される浮動小数点数です。



したがって、正の浮動小数点数を取得して整数に変換し、「1」を追加すると、この算術で表現可能な次の数が得られます。 Cでは、これは次のように実行できます。

float a=0.5; int n = *((int*) &a); float b = *((float*) &(++n)); printf(" %e  : %e,  (%e)\n", a, b, ba);
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

このコードは、32ビット整数のアーキテクチャでのみ機能します。

4.浮動小数点演算の落とし穴

今すぐ練習します。 プログラミング時には特に注意する必要がある浮動小数点演算の機能を検討してください。

4.1丸め

特に倍精度を使用している場合は、最新の浮動小数点演算の丸め誤差によるエラーに対処することは困難です。 IEEE754規格の丸め規則では、算術演算の結果は、正確な値で実行され、この形式で表される最も近い数値に丸められるかのようになります。 これにはALUからの追加の努力が必要であり、一部のコンパイラオプション（gccの「-ffast-math」など）はこの動作を無効にする場合があります。 IEEE754の丸めの機能：

• 標準で最も近いものへの丸めは、以前のようには行われません。 0.5を1（切り上げ）に丸めると、丸め誤差が無限大に増加する一連の操作があることが数学的に示されています。 したがって、IEEE754は偶数に丸め規則を適用します。 したがって、12.5は12に、13.5は14に丸められます。
• 浮動小数点演算の丸めに関する最も危険な操作は、減算です。 近い数字を減算すると、有効数字が失われる可能性があります。
  
  相対誤差が大幅に増加する可能性があります。
• 多くの一般的な数式について、数学者は丸め誤差を大幅に減らすことができる特別な形式を開発しました。 たとえば、式「x ² -y ² 」の計算は、式「（xy）（x + y）」を使用して最適に計算されます。

4.2非連想算術演算

浮動小数点演算では、ルール（a * b）* c = a *（b * c）は算術演算には適用されません。 例えば



（10 ²⁰ +1）-10 ²⁰ = 0≠（10 ²⁰ -10 ²⁰ ）+ 1 = 1



数値を合計するプログラムがあるとします。

 double s = 0.0; for (int i=0; i<n; i++) s = s + t[i];
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

一部のコンパイラは、デフォルトで、複数のALUを同時に使用するようにコードを書き換えることができます（nを2で割ると仮定します）。

 double sa[2], s; sa[0]=sa[1]=0.0; for (int i=0; i<n/2; i++) {    sa[0]=sa[0]+t[i*2+0];    sa[1]=sa[1]+t[i*2+1]; } S=sa[0]+sa[1];
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

加算の操作は連想的ではないため、これらの2つのプログラムは異なる結果を生成する可能性があります。

4.3数値定数

すべての10進数に2進浮動小数点表現があるわけではないことに注意してください。 たとえば、数値「0.2」は単精度で「0.200000003」として表されます。 したがって、「0.2 + 0.2≈0.4」。 別の絶対誤差

ケースは高くないかもしれませんが、サイクルでそのような定数を使用すると、累積エラーを取得できます。

4.4最低2つの値の選択

最小値を選択する必要がある2つの値があるとします。 Cでは、これは次のいずれかの方法で実行できます。

1. x <y？ x：y
2. x <= y？ x：y
3. x> y？ y：x
4. x> = y？ y：x

多くの場合、コンパイラは1つのプロセッサ命令で実行されるため、これらを同等と見なし、常に最初のオプションを使用します。 ただし、±0とNaNを考慮すると、これらの演算はどのような意味でも同等ではありません。

x	y	x <y？ x：y	x <= y？ x：y	x> y？ y：x	x> = y？ y：x
+0	-0	-0	+0	+0	-0
ナン	1	1	1	ナン	ナン

4.5数値の比較

浮動小数点数を処理するときに非常に一般的なエラーは、同等性をチェックするときに発生します。 例えば

 float fValue = 0.2; if (fValue == 0.2) DoStuff();
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ここでの間違いは、0.2は正確なバイナリ表現を持たないことです。次に、0.2は倍精度の定数であり、変数fValueは単一であり、この比較の動作についての保証はありません。



最良の、しかしまだ間違った方法は、差を許容絶対誤差と比較することです：

 if (fabs(fValue – fExpected) < 0.0001) DoStuff(); // fValue=fExpected?
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

このアプローチの欠点は、この数値自体の増加に伴い、数値を表すエラーが増加することです。 したがって、プログラムが「10000」を予期している場合、上記の等式は最も近い隣接番号（10000,000977）については満たされません。 これは、プログラムが単精度から倍精度への変換を行う場合に特に当てはまります。



適切な比較手順を選択するのは難しく、興味のある読者です。 ブルース・ドーソンの記事を参照してください。 浮動小数点数と整数変数への変換を比較することを提案しています。 これは、移植性のある方法ではありませんが、最高です。

 bool AlmostEqual2sComplement(float A, float B, int maxUlps) {    // maxUlps        ,     // NaN          assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);    int aInt = *(int*)&A;    //    aInt,  ,         if (aInt < 0) aInt = 0x80000000 - aInt;    //aInt &= 0x7fffffff; //(.   Vayun)    //   bInt    int bInt = *(int*)&B;    if (bInt < 0) bInt = 0x80000000 - bInt;    /*aInt &= 0x7fffffff;*/    unsigned int intDiff = abs(aInt - bInt); /*(.   Vayun)*/    if (intDiff <= maxUlps)        return true;    return false; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

このプログラムでは、maxUlps（Units-In-Last-Placeから）は、チェック値と期待値の間にある浮動小数点数の最大数です。 この変数のもう1つの意味は、比較対象の数字の2進数（最下位から始まる）の数を逃すことです。 たとえば、maxUlps = 16は、下位4ビット（log ₂ 16）が一致しない可能性があることを意味し、数値は依然として等しいと見なされます。 さらに、数値10000と比較すると、絶対誤差は0.0146に等しくなり、0.001と比較すると、誤差は0.00000001（10 ^-8 ）未満になります。

5. IEE754のサポートの完全性を確認します

プロセッサがIEEE754標準に完全に準拠している場合、標準データ型（Cのfloat / doubleなど）を使用するプログラムは、異なるコンピューターで同じ結果を生成すると思いますか？ あなたは間違っています。 移植性とコンプライアンスは、コンパイラと最適化オプションの影響を受けます。 William KahanはC（Fortranのバージョンがあります）でプログラムを作成しました。これにより、IEEE754のリンク「アーキテクチャ+コンパイラ+オプション」かどうかを確認できます。 「浮動小数点パラノイア」と呼ばれ、そのソースコードはダウンロードできます 。 同様のプログラムがGPUで利用可能です。 そのため、たとえば、Intelコンパイラ（icc）はデフォルトで「緩和された」IEEE754モデルを使用するため、すべてのテストが実行されるわけではありません。 オプション「-fp-model precise」を使用すると、標準に厳密に準拠してプログラムをコンパイルできます。 GCCコンパイラには「-ffast-math」オプションがあり、これを使用するとIEEE754の不一致が発生します。

おわりに

最後に、有益な物語。 GPUでテストプロジェクトに取り組んでいたとき、1つのプログラムのシリアルバージョンとパラレルバージョンがありました。 ランタイムを比較すると、300倍の加速が得られたため非常に満足しています。 しかし後になって、GPUでの計算は「崩れ」、NaNに変わり、GPUでの計算は通常の数値よりも高速であることが判明しました。 もう1つ興味深いのは、GPUエミュレーター（CPU上）の同じプログラムが正しい結果を生成したが、GPU自体ではなかったことです。 その後、このGPUがIEEE754標準を完全にサポートしておらず、直接的なアプローチが機能しなかったことが問題であることが判明しました。



現在、浮動小数点演算はほぼ完璧です。 ほとんどの場合、素朴なアプローチが機能し、すべての機能を考慮していないプログラムは正しい結果をもたらします。また、記載されている落とし穴はエキゾチックなケースのみに関係します。 しかし、常に警戒しなければなりません。コンピューター数学などの問題では、熊手を踏むのは簡単です。



PS重要な点についてはuqlockに感謝します。 お金を浮動小数点数として保存することはできません。 この場合、重要なカテゴリを区別できません。 プログラミング言語に固定小数点データ型がない場合、状況から抜け出し、お金を整数として保存できます。これは、ペニー（場合によってはペニーの一部）を意味します。



PPSユーザーに誤植とエラーを発見しました： gribozavr 、 kurokikaze 、 Cenness 、 TheShock 、 perl_demon 、 GordTremor 、 fader44 、 DraculaDis 、 icc 、 f0rbidik 、 Harkonnen 、 AlexanderYastrebov 、 Vayun 、 EvilsInterrupt ！

文学

1. IEE754標準の開発に関するWilliam Kahanへのインタビュー 。
2. すべてのコンピューター科学者が浮動小数点演算について知っておくべきこと、David Goldberg-数学計算の本。
3. 浮動小数点数の比較、ブルース・ダウソン 。