浮動小数点演算について
浮動小数点ルーチンの使用 への移動
浮動小数点数を取り扱うには、データの内部表現を理解しておく必要があります。プログラムを作成する際には、次のような有限精度の問題点を認識している必要があります。
- 整数値(整数型)の場合は、オーバーフローが起こり得ることを考慮する必要があります。
- 浮動小数点値(単精度、倍精度、拡張精度)の場合は、精度の損失が起こり得ることを考慮する必要があります。
目次
有限精度のもたらす影響
浮動小数点数の精度の損失は、以下のような複数の原因で生じるおそれがあります。
- 浮動小数点リテラル(たとえば 0.1)を浮動小数点変数に代入する場合、目的とする値を誤差なしで格納できるだけのビット数がその浮動小数点変数にない可能性があります。その浮動小数点リテラルには、非常に多くのビット数が必要かもしれません。あるいは、無限精度の表現を行うために無限個のビットさえ必要かもしれません。
uses System.SysUtils; var X: Single; Y: Double; begin X := 0.1; Y := 0.1; Writeln('X =', X); Writeln('Y =', Y); Writeln(Format('Exponent X: %x', [X.Exp])); Writeln(Format('Mantissa Y: %x', [Y.Mantissa])); ReadLn; end.
- コンソール出力は次のとおりです。
X = 1.00000001490116E-0001 Y = 1.00000000000000E-0001 Exponent Y: 3FB Mantissa Y: 1999999999999A
- 上記のコードから、Single 精度表現の 9 桁目にエラーが確認できます。 Double 精度表現にもエラーがあります。 このエラーを表すために、Double 精度値の、生の Exponent(指数) と Mantissa(仮数) を使用してみましょう。 16 進数
$3FB
の 指数には、10 進数表記1019
があります。 Double 数の内部表現には、1023 に等しい バイアス があるため、指数 = 1019 - 1023 = -4
となります。 16 進数$1999999999999A
の 仮数には、2進数表現11001100110011001100110011001100110011001100110011010
があります。 つまり、Y の 2 進数表現は1.1001100110011001100110011001100110011001100110011010*2-4
または0.00011001100110011001100110011001100110011001100110011010*2-4
となります。 この数は、Double 精度の近似値である点に注意してください。 厳密な0.1
の数は、無限循環小数0.0(0011)
で表します。
- ただし、こうして生じる可能性のある最大誤差は 0.5 ulp であることに注意することが大切です。
- 浮動小数点演算を実行する場合、各ステップ(演算)で、それ固有の誤差が生じるおそれがあります。その理由は、一部の演算については、計算結果を完全な精度で格納できないからです。たとえば、2 つの数値(S1 ビットと S2 ビット)を乗算する場合(これは整数型の場合にも浮動小数点型の場合にも当てはまります)、結果を完全な精度で表現するには S1 + S2 ビットが必要になります。
- 演算で生じる誤差の "量" はプロセッサ モデルと演算の種類によって異なります。加法演算で生じる誤差は比較的小さくなります。乗法演算で生じる誤差は比較的大きくなります。
浮動小数点数の精度の損失(誤差)は計算を通じて伝播されますが、それでも正しいアルゴリズムを設計するのはプログラマの役割であることを理解することが重要です。
浮動小数点変数は、2 のべき乗のスケールを持つ整数変数と見なすことができます。浮動小数点変数に極値を代入すると、スケールは自動的に調整されます。浮動小数点変数が決してオーバーフローしないような印象を持つことがあるのは、そのためです。そして、実際にそのとおりなのですが、他方では、危険な問題があります。つまり、浮動小数点変数はかなりの誤差を蓄積したり、非正規化数となるおそれがあります。
より大きいデータ型の使用
整数のオーバーフローや浮動小数点数の精度低下の問題(一般に有限精度の影響)を解決する最も簡単な方法は、同類(整数か浮動小数点数)で容量のより大きいデータ型を使用することです。つまり、ShortInt でオーバーフローする場合、LongInt、FixedInt、または Int64 に容易に切り替えることができます。同様に、単精度浮動小数点型で精度が不十分な場合、倍精度浮動小数点型に切り替えることができます。しかし、考慮すべき事項が次のように 2 つあります。
- 記憶容量がより大きいデータ型でもまだ不十分な場合があります。
- 記憶容量がより大きいデータ型の場合は、より多くのメモリが必要で、演算に必要な CPU サイクルも増える可能性があります。
制御設定
32 ビット プラットフォームでは、x87 FPU 制御ワード(CW)の 2 ビット分が丸めモードの指定に割り当てられます。詳細については、『Intel® 64 and IA-32 Architectures Software Developer's Manual Volume 1: Basic Architecture > 8.1.5.3 Rounding Control Field(64/IA-32 インテル® アーキテクチャ・ソフトウェア・デベロッパーズ・マニュアル 上巻:基本アーキテクチャ 8.1.5.3 項「丸め制御フィールド」)』を参照してください。64 ビット プログラムの場合は、SSE 制御レジスタ MSXCSR で丸めモードが指定されます。これらは、System.Math.SetRoundMode を使って変更できます。
浮動小数点変数を操作する一部の RTL 関数は、FPU 丸めモードの影響を受ける可能性があります。RTL ルーチンの結果が FPU 制御ワードに基づいて実際にどう変わるかは、実装されるアルゴリズムしだいです。結果をターゲット型に合わせるのに丸めが必要なあらゆる演算に、丸めは影響を及ぼします。たとえば、浮動小数点乗算はほとんど常に丸めを伴います。関数が多数の浮動小数点乗算から成る場合は、丸めモードの影響を強く受けます。
丸めモードは、区間演算を実装するのに使用されることがあります。つまり、大まかに言えば、同じアルゴリズムを切り上げモードで実行したあと、切り捨てモードで同じことを繰り返し、次に、2 つの結果の違いを調べます。これで、丸めと精度不足で生じる可能性がある誤差がわかります。
使用例
財務計算
IEEE 浮動小数点は財務計算には妥当でない可能性があります。精度の要件が通常非常に厳しいからです。整数型(プリミティブ整数か Currency)または BCD 型の使用を検討しなければなりません。
Data.FmtBcd ユニットは BCD 演算のサポートを提供します。BCD 形式には、各 10 進数字(基数 10 の数字)が 4 ビット メモリ(ニブル)でコード化される、という重要な特徴があります。
次のコードは、TBcd 値を便宜上バリアントとして使用する方法を示しています。
Delphi の場合:
var
X: Variant;
begin
X := VarFMTBcdCreate('0.1', 16 { Precision }, 2 { Scale });
Writeln(VarToStr(X));
// ...
C++ の場合:
#include <Variants.hpp>
#include <FMTBcd.hpp>
int _tmain(int argc, _TCHAR* argv[]) {
Variant x = VarFMTBcdCreate("0.1", 16 /* Precision */, 2 /* Scale */);
printf("%ls", VarToStr(x).c_str());
// ...
コンソール出力は次のとおりです。
0.1
上記のコードでは、BCD 変数のおかげでテキスト形式から数値形式への変換が完ぺきであることがわかります。
財務計算には、Currency 型を使用できます。Currency 型は本質的に 10000 でスケールされた整数です(10 で割り切れることになります)。Currency 変数には 10進数字を 4 つ格納でき、この上限を超えるものはすべて丸められます。
Delphi の場合:
var
X, Y: Currency;
begin
X := 0.1;
Y := 0.00001;
Writeln(X);
Writeln(Y);
// ...
C++ の場合:
#include <System.hpp>
int _tmain(int argc, _TCHAR* argv[]) {
Currency x, y;
x = 0.1;
y = 0.00001;
printf("%2.14E\n", (double)x);
printf("%2.14E\n", (double)y);
// ...
コンソール出力は次のとおりです。
1.00000000000000E-0001 0.00000000000000E+0000
C++ による Currency の実装については $BDS\include\rtl\syscurr.h
を参照してください。
Currency の区間は [-922337203685477.5807; 922337203685477.5807] です。
物理(科学技術)計算/シミュレーション
科学技術計算には、一般に、相当な計算が必要であり、浮動小数点精度を上げることは望ましくない可能性があります。拡張精度演算はあまりサポートされていません(Delphi for x64 を参照)。
SSE を使用する場合は、SSE レジスタには 2 つの倍精度変数または 4 つの単精度変数を格納できることを覚えておいてください。したがって、倍精度演算よりも多くの単精度演算を並行実行できます。
非常に興味深く有用なアプローチは、単精度浮動小数点を使用しつつ、累積誤差を定期的に減らすことです。多くのアプリケーションでは、精度の小さい損失に耐えることができ、そのずれを何とかして解消することが重要です。そのような実装の例には、四元数を使用した 3D 空間回転があります。『Physically Based Modeling > Rigid Body Dynamics (ONLINE SIGGRAPH 2001 COURSE NOTES)』(SIGGRAPH 2001 オンライン コース資料: 物理ベース モデリング|剛体力学)を参照してください。
デジタル信号処理
すべての DSP 変数には、一般に、誤差が含まれています。データの精度と計算量の間の最適なトレードオフを検討しなければなりません。
全般的なまとめ
"見えているものが実態ではない"
ソース コードに 10 進数字で書かれている浮動小数点数と、画面に表示されている浮動小数点数は、メモリ内のデータとはおそらく異なります。コンソールでの表示がメモリ内のデータを正確に表していると決めてかからないでください。10 進から 2 進への変換(およびその逆変換)は、どのような場合にも完ぺきには行えません。
整数型変数、BCD 型変数、Currency 型変数を使用すると、IEEE 浮動小数点表現の誤差を避けることができます。
データ フローの理解
x86 上の Delphi では、単精度(Single 型)浮動小数点式の中間結果は常に Extended 型として暗黙に格納されます。
単精度(Single 型)浮動小数点値のみ関与するすべての x64 算術演算および算術式では、デフォルトで、中間結果を倍精度(Double 型)浮動小数点値として格納することにより、高い精度を保っています。その結果、そのような演算は、オペランドが倍精度と明示される場合よりも低速になります(コンパイル済みコードでは、演算ごとに単精度の値が倍精度に変換されるため)。実行速度が一番の関心事である場合は、問題のコードに {$EXCESSPRECISION OFF} 指令を付けて、中間の倍精度値を使用しないようにすることができます。そうでない場合は、デフォルトの指令({$EXCESSPRECISION ON})のままにして、結果値の精度を向上させることをお勧めします。{$EXCESSPRECISION OFF} 指令は、x64 のターゲット CPU の場合にのみ効果があります。
C++ では、浮動小数点リテラルは、単精度か倍精度のどちらかの浮動小数点数を表すことができます。そのどちらであるかは、f
サフィックスの有無で決まります。C++ で書かれた浮動小数点リテラルの末尾に f
が付加されている場合、コンパイラは単精度を選びます。中間値の精度を理解するには、公開されている ISO 標準を参照してください。
浮動小数点演算は結合的でない可能性あり
演算子ごとに誤差が生じるため、計算の実行順序が重要な意味を持ちます。
『CERT C Secure Coding Standards, Recommendation FLP01-C』(CERT C セキュア コーディング標準、勧告 FLP01-C)を参照してください。
浮動小数点例外
浮動小数点演算は、floating-point overflow、0 による除算、非正規化値、NaNs の生成、他の不正な浮動小数点演算の実行のような、不正な状況を引き起こす可能性があります。 大抵は、そのような状況は、浮動小数点例外を発生させます。 System.Math ユニットは次のものを提供します:
- 発生し得る浮動小数点例外の集合。
- 現在のプラットフォームで現在アクティブになっている浮動小数点例外を取得するメソッド。
- 現在マスクされている例外を現在の浮動小数点演算ハードウェアから取得するメソッドと、現在の浮動小数点演算ハードウェアでマスクされる例外を設定するメソッド。
浮動小数点演算ハードウェアによって、浮動小数点例外を制御する手段は異なります。
- Intel 32 ビット プロセッサでは、その手段は FPU 制御ワードです。
- Intel 64 ビット プロセッサでは、その手段は SSE 制御ワードです。
- ARM アーキテクチャでは、浮動小数点例外をサポートしていません。そのため、ARM アーキテクチャでは、すべての浮動小数点例外が常にマスクされています。
関連項目
- 浮動小数点数制御ルーチン
- 浮動小数点数の丸めルーチン
- 浮動小数点比較ルーチン
- 実数型の内部表現
- 数値型の内部表現
- 浮動小数点例外
- System.Round
- System.Math.SetRoundMode
- 浮動小数点の精度の制御(x64 用 Delphi コンパイラ指令 {$EXCESSPRECISION})
- Extended 型の互換性({$EXTENDEDCOMPATIBILITY} Delphi)
- SetRoundMode(Delphi)コード例
- How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?(William Kahan 著: 浮動小数点計算における丸めのやみくもな評価はいかに無意味か)(英語版)
- What Every Computer Scientist Should Know About Floating-Point Arithmetic(David Goldberg 著: すべてのコンピュータ科学者に必須の浮動小数点演算の知識)(英語版)
- Numerical stability(Wikipedia.org: 数値的安定性)
- CERT C Coding Standards, Floating Point (FLP)(CERT C コーディング標準、浮動小数点(FLP))
- Floating-point to fixed-point conversion for DSP(Google プロジェクト ホスティング: DSP を対象とした浮動小数点から固定小数点への変換)(英語版)
- Floating-Point Rounding(Google プロジェクト ホスティング: 浮動小数点丸め)(英語版)