Menu

記事

概要

英語版
このページの英語版を見る

機械のすぐそば:D言語の浮動小数点演算子

目次
  1. はじめに
  2. 浮動小数点数の謎を解く
  3. 条件数
    1. 加算と減算
    2. 乗算と除算
    3. 累乗と対数
  4. NaNペイロード
  5. IEEE丸めモード
  6. IEEE 例外ステータスフラグ
  7. 浮動小数点数と「純粋なnothrow」
  8. 結論
  9. 参考文献および関連情報

はじめに

ドン・クラッグストン著

コンピュータはもともと、数学を実行するための装置として考案された。初期のコンピュータは、ほとんどの時間を方程式の解法に費やしていた。現在では、工学や科学のコミュニティはコンピューティングの世界のごく一部を占めるに過ぎないが、その昔の素晴らしい遺産が残されている。現在では、ほとんどのコンピュータが数学計算を正確かつ極めて高速に実行するための優れたハードウェアを搭載している。残念ながら、ほとんどのプログラミング言語では、プログラマーがこのハードウェアを十分に活用することが難しい。さらに大きな問題は、ドキュメントが不足していることである。多くの数学プログラマーにとっても、浮動小数点演算の側面は依然として謎に包まれたままである。

システムプログラミング言語として、プログラミング言語Dはプログラマーとコンパイラー、プログラマーとマシン間の障壁を取り除くことを試みている。この哲学は浮動小数点演算のサポートにおいて特に顕著である。ハードウェアを正確に理解することの重要性を示す個人的な逸話がある。

初めて浮動小数点演算で悪夢のような問題に遭遇したのは、100回に1回くらいの割合でハングアップするC++プログラムだった。最終的に、問題は時折終了しなくなるwhileループにあることが判明した。コードの要点はリスト1に示す。

double q[8];
...
int x = 0;
while (x < 8)
{
    if (q[x] >= 0) return true;
    if (q[x] < 0) ++x;
}
return false;

当初、私はこの無害に見えるループがなぜ失敗するのか全く理解できなかった。しかし最終的に、qが適切に初期化されていなかったことが判明した。q[7]にはランダムな文字列が含まれていた。時折、その文字列はすべてのビットが設定されており、q[7]はNaN(非数)であった。NaNはコンパイラのドキュメントには記載されておらず、私が発見した唯一の情報はインテルのアセンブリ命令セットのドキュメントに記載されていた。NaNを含む比較はすべて偽となるため、q[7]は>= 0でも< 0でもなく、プログラムが停止してしまう。この不愉快な発見をするまで、私はNaNの存在すら知らなかった。浮動小数点数型はすべて正、負、またはゼロのいずれかであると誤って思い込んでいた私は、まさに「愚か者の楽園」に暮らしていたのだ。

D言語では、私の経験はまったく異なるものになっていたことだろう。浮動小数点数の「奇妙な」機能は、この言語ではより可視性が高く、数値プログラマーの教育効果も向上する。 未初期化の浮動小数点数型変数はコンパイラによってNaNに初期化されるため、問題のループは断続的にではなく、毎回失敗することになる。 D言語の数値プログラマーは通常、浮動小数点数の「無効」例外を有効にした状態でプログラムを実行する。そのような状況下では、プログラムが未初期化変数にアクセスした時点でハードウェア例外が発生し、デバッガーが呼び出される。 浮動小数点演算の「奇妙な」機能に簡単にアクセスできることで、より教育を受けたプログラマーが生まれ、混乱が減り、デバッグが迅速化され、テストが改善され、そして願わくば、より信頼性が高く正確な数値プログラムが作成されることになる。 本記事では、"プログラミング言語 D"における浮動小数点演算のサポートについて、その概要を説明する。

浮動小数点数の謎を解く

Dは、すべての組み込み浮動小数点型がIEEE 754算術に準拠し、動作が完全に予測可能であることを保証している(これは、すべてのプラットフォームで同一の結果が得られることとは異なることに注意)。IEEE 754-2008は、浮動小数点算術に関するIEEE 754規格の最新改訂版である。Dは、754-2008への完全準拠に向けて進歩している。

現在DでサポートされているIEEE標準浮動小数点型は、float とdouble である。さらに、Dはreal 型もサポートしており、これはCPUがサポートしている場合は「IEEE 80ビット拡張」となり、そうでない場合はdouble と同じである。今後、754-2008の新しい型であるquadruple 、decimal64 、decimal128 が追加される予定である。

これらの型の特性は、プロパティを介して簡単に言語からアクセスできる。例えば、float.max は浮動小数点数に格納できる最大値であり、float.mant_dig は仮数部に格納される桁数(ビット)である。

D言語で数学を理解するには、IEEE浮動小数点演算の基本的な理解が必要である。基本的に、これは無限に存在する実数を少数のバイトにマッピングするものである。IEEE 32ビット浮動小数点数として表現できるのは、わずか40億個の異なる数値である。このような哀れなほど小さな表現空間であっても、IEEE浮動小数点演算は数学的な実数が使用されているかのような錯覚を維持するのに非常に優れた働きをするが、その錯覚が崩れる場合を理解しておくことは重要である。

ほとんどの問題は、これらの表現可能な数値の分布から生じる。IEEE数値軸は数学的な数値軸とは全く異なる。

     +     +-----------+------------+    ..   +    ..    +----------+----------+     +       #
-infinity -float.max  -1  -float.min_normal   0   float.min_normal  1  float.max infinity  NaN

IEEE 数直線は、-1 から +1 の間に半分が位置していることに注目してほしい。 0 から 0.5 の間には 10 億の表現可能な浮動小数点数があるが、0.5 から 1 の間には 800 万しかない。 これは精度に関して重要な意味を持つ。ゼロ付近では有効な精度が非常に高いのだ。 このことを利用して、-1 から +1 の範囲の数値を個別に扱ういくつかの例を提示する。

また、特別な数値にも注目してほしい。&plusmn;&infin;; つまり、&plusmn;float.min_normalと0の間の、いわゆる「サブノーマル」は、縮小された精度で表現される。さらに、+0と-0の2つのゼロが存在すること、そして最後に「NaN」(「Not-a-Number」)というナンセンスな値がある。NaNは、リスト1で多くの問題を引き起こした。

なぜ NaN が存在するのか? 浮動小数点演算における未定義の動作を排除するという重要な役割がある。 これにより、浮動小数点演算は完全に予測可能となる。int 型では、3/0 の演算によりゼロによるハードウェア除算トラップハンドラが呼び出され、プログラムが終了する可能性があるが、浮動小数点の除算 3.0/0.0 は &infin; となる。 数値のオーバーフロー(例:real.max*2 )も &infin; となる。アプリケーションによっては、∞は完全に妥当な結果である場合もあるが、通常はエラーを示す。0.0 / 0.0 のような無意味な演算は NaN を生み出すが、プログラムの制御は失われない。一見したところ、∞や NaN は不要のように思える。整数の場合と同様に、エラーとしてしまえばよいではないか。結局のところ、ゼロによる除算を避けるのは簡単である。すべての除算の前に分母がゼロでないかチェックすればよいだけだ。真の難しさはオーバーフローから生じる。掛け算でオーバーフローが起こるかどうかを事前に判断するのは極めて難しい。

ある種の異常を回避し、「x - y == 0 であるのは、x == y の場合のみ」といった重要な関係を維持するには、サブノーマルが必要である。

∞はオーバーフローによって生じる可能性があるため、+∞と-∞の両方が必要となる。また、x と1/(1/x) が常に同じ符号ビットを持つようにするためには、+0と-0の両方が必要となる。しかし、ほとんどの他のケースでは、+0と-0に違いはない。

これらの「特別な値」は通常、効率的ではないことに留意すべきである。例えば、x86マシンでは、NaN、無限大、非正規数を含む乗算は通常の数値の演算よりも20倍から50倍遅くなる可能性がある。数値コードが予想外に遅い場合、これらの特殊な値を意図せず多数作成している可能性がある。後述する浮動小数点例外トラップを有効にすると、これを確認する簡単な方法となる。

マシンが何をしているのかを不明瞭にする最大の要因のひとつは、2進数と10進数の変換にある。結果を表示する際に"%a" フォーマットを使用することで、この要因を排除できる。これは非常に貴重なデバッグツールであり、浮動小数点アルゴリズムを開発する際には非常に役立つ。0x1.23Ap+6 16進浮動小数点フォーマットは、ソースコードでも使用でき、入力データが意図した通りのものであることを確認できる。

浮動小数点数の量子化の性質

取り得る値が限られているという事実により、数学的な実数では不可能な演算が可能になる。ある数値 x が与えられた場合、nextUp(x) は x より大きな次の表現可能な数値を返す。nextDown(x) は x より小さな次の表現可能な数値を返す。

数値解析者はしばしばエラーを「最後の桁の単位」(ulp)という驚くほど微妙な用語で表現するが、これはしばしば無頓着に使用される。[脚注: 最も正式な定義は[J.-M. Muller, "On the definition of ulp(x)",INRIA Technical Report 5504 (2005)]で見られる。x が、有限の連続する浮動小数点数型Fの2つの数値aとbの間の実数であり、そのいずれとも等しくない場合、ulp(x)=abs(b-a)となる。それ以外の場合、ulp(x) =x*F.epsilon となる。さらに、ulp(NaN)はNaNであり、ulp(&plusmn;F.infinity) = &plusmn;F.max*F.epsilon となる。 私は、もっとずっとシンプルな定義を好む。2つの数値 x と y の間の ulp の差は、x から y に移動するために nextUp() または nextDown() を呼び出す必要がある回数である。 [脚注:x または y が浮動小数点数ではなく実数である場合、これは整数にはならない。] Dライブラリ関数feqrel(x, y) は、xとyの間に等しいビット数を返す。これは、精度の損失問題をチェックする簡単な方法である。

浮動小数点数の量子化には、興味深い結果をもたらす性質がある。

条件数

nextUp を使用すれば、条件数を簡単に概算できる。

real x = 0x1.1p13L;
real u = nextUp(x);

int bitslost = feqrel(x, u) - feqrel(exp(x), exp(u));

これは、xの値が非常に大きい場合、xの1ビットのエラーがexp(x)の12ビットの精度を破壊することを示している。 エラーは最後の桁で約6000単位増加している。したがって、このxの値における条件数は6000である。

float、double、realのセマンティクス

市場を支配するx86マシンでは、浮動小数点演算は伝統的に8087数学コプロセッサの派生品で実行されてきた。これらの「x87」浮動小数点ユニットは、IEEE754演算を初めて実装したものである。SSE2命令セットはx86-64プロセッサ用の代替手段であるが、浮動小数点32ビットx86マシンではx87が唯一のポータブルなオプションである(32ビットAMDプロセッサではSSE2はサポートされていない)。

x87は、他のほとんどの浮動小数点ユニットと比較すると異例である。x87は80ビットのオペランドのみをサポートしており、以降「real80」と呼ばれる。double およびfloat のオペランドはすべてまず80ビットに変換され、すべての算術演算は80ビット精度で実行され、必要に応じて結果は64ビットまたは32ビット精度に縮小される。つまり、結果は最大64ビット演算しかサポートしないマシンよりもはるかに正確になる可能性がある。しかし、ポータブルコードの記述には課題もある。 (脚注:x87では仮数部の長さをdouble やfloat と同じに減らすことができるが、指数部はreal80のままなので、非正規数では異なる結果が得られる。double を正確にエミュレートするには、浮動小数点演算コードの速度が1桁遅くなる)。

x87ファミリーを除くと、80ビット浮動小数点演算をサポートしているのは、モトローラ68K(ただし、コールドファイアはサポートしていない)とItaniumプロセッサのみである。

同様の問題は、FMA(融合乗算累積)命令にも関連しており、これはPowerPC、Itanium、Sparc、Cellなど、ますます多くのプロセッサで利用可能になっている。このようなプロセッサでは、x*y + z などの式を評価する際、x*y は通常の2倍の精度で実行される。そうでなければ完全な精度の損失を招くような計算も、代わりに正確に計算される。 高度なシステムプログラミング言語の課題は、すべてのプラットフォーム上で予測可能な動作を提供し、かつ利用可能なハードウェアを十分に活用できる抽象化を作成することである。

この状況に対するDの取り組みは、以下の観察から生じている。

  1. すべてのプロセッサで同一の動作を確実に実現するには、パフォーマンス面で非常にコストがかかる。特に、x87では致命的である。
  2. 非常に多くのプログラムは、特定のプロセッサでしか動作しない。移植性を確保するために、より精度の高いハードウェアの使用を否定するのは残念なことである。
  3. 移植性と高精度の要件は同時に必要とされることはない。double の精度が不十分である場合、一部のプロセッサのみ精度を上げることは助けにならない。
  4. 言語は特定のプロセッサの特定の機能に縛られるべきではない。

重要な設計目標は、使用するプロセッサに関係なく、double 型のみをサポートするシステムで得られる精度よりも悪くならないようなコードを記述できることである。

脚注:real は、Javaプログラミング言語のボルネオ提案における「固有」に近い[Ref Borneo]。

x*y + z*w を評価することを検討する。x, y, z とw はdoubleである。

  1. double r1 = x * y + z * w;
  2. double a = x * y; double r2 = a + z * w;
  3. 実数 b = x * y; 実数 r3 = b + z * w;

最適化の過程で、(2)と(3)は(1)に変換される可能性があるが、これは実装に依存する。 (2)の場合、余分な丸め処理が導入されるため、特に問題となる。

「単純な」CPUでは、r1==r2==r3である。この値をr0と呼ぶ。 PowerPCでは、r2==r3であるが、FMAを使用できるため、r1は他の値よりも正確である可能性がある。 x86では、r1==r3であり、PowerPCの場合ほどではないが、r0よりも正確である可能性がある。 しかし、r2はr0よりも正確でない可能性がある。

中間値としてreal を使用することで、double のみをサポートする単純なCPUよりも結果が悪くなることは決してないことが保証される。

組み込み型の特性

基本的な浮動小数点の特性は、epsilon 、min_normal 、max である。6つの整数特性は、これら3つのlog2またはlog10に単純に置き換えられる。

float double real80 倍精度 decimal64 decimal128
イプシロン 0x1p-23 0x1p-52 0x1p-63 0x1p-112 1e-16 (1p-54) 1e-34 (1p-113)
[min_normal 0x1p-126 0x1p-1022 0x1p-16382 0x1p-16382 1e-383 1e-6143
...最大) 0x1p+128 0x1p+1024 0x1p+16384 0x1p+16384 1e+385 1e+6145
バイナリプロパティ
mant_dig 24 53 64 113 53 112
min_exp -125 -1021 -16381 -16381 &nbsp;
最大経験値 +128 +1024 +16384 +16384 &nbsp;
小数点以下のプロパティ
6 15 18 33 16 34
min_10_exp -37 -307 -4932 -4932 -382 -6142
max_10_exp +38 +308 +4932 +4932 385 +6145

コンパイル時に異なるCPUに適応するコードを記述する場合は、static if をmant_dig プロパティとともに使用する。例えば、static if (real.mant_dig==64) は80ビットの実数が使用可能である場合に真となる。 バイナリ型では、dig プロパティは有効な10進数の最小桁数のみを与える。表現可能なすべての数値に一意の10進表現を確保するには、2桁の追加が必要である。同様に、10進数では、mant_dig は有効な2進数の下限となる。

浮動小数点数型に関する便利な関係式F 、ここでx とy は型であるF

加算と減算

乗算と除算

2進数型に対する累乗と対数

NaNペイロード

IEEE 754規格によると、「ペイロード」はNaNの仮数部に格納することができる。このペイロードには、それがどのように、またはなぜ生成されたかに関する情報が含まれている可能性がある。歴史的に、この潜在的に強力な機能を利用したプログラミング言語はほとんどない。D言語では、このペイロードは正の整数で構成される。

決してポインタを整数ペイロードとしてNaN内に格納してはならない。ガベージコレクタがそれを発見できなくなるからだ!

IEEE丸めモード

丸めモードはスコープ内で制御される。丸めモードはそのスコープの終了時に以前の状態に復元される。 4つの丸めモードを設定できる。デフォルトのモードである「最も近い値に丸める」は統計的には最も正確であるが、直感的ではない。同値の場合、結果は偶数に丸められる。

丸めモード rndint(4.5) rndint(5.5) rndint(-4.5) 注釈:
最も近い 4 6 -4 偶数に切り上げる
切り捨て 4 5 -5
切り上げ 5 6 -4
ゼロに丸める 4 5 -4

丸めモードを変更する理由はほとんどない。 切り上げモードと切り捨てモードは、特に間隔演算を迅速に実装できるようにするために作成されたもので、特定のライブラリでは不可欠であるが、それ以外ではほとんど使用されない。 ゼロに丸めるモードは、浮動小数点数を整数に変換する際に使用される。モードの切り替えは特にIntelマシンでは遅いので、cast(int) の動作を正確に複製するために、ゼロに丸めるモードに切り替えると便利である。

丸めモードを変更する理由として一般的に挙げられるのは、数値の安定性を確認するという理由だけである。丸めモードを変更すると計算結果が大幅に異なる場合は、丸め誤差が発生していることを示す明確な兆候である。

IEEE 例外ステータスフラグ

すべてのIEEE準拠プロセッサには、プログラムが認識する必要がある「奇妙な」事象が発生したことを示す特別なステータスビットが含まれている。例えば、ieeeFlags.divideByZero はゼロで割ることによって無限大が生成されたかどうかを示す。これらは「粘着性」ビットであり、一度設定されると、明示的にクリアされるまで設定された状態が維持される。計算の最後に一度だけこれをチェックするだけで、ほとんど失敗することのない何千もの比較を回避できる可能性がある。

以下は、検出可能な奇妙なもののリストである:

無効
これは、NaN が生成された場合に設定される。これは、∞ - ∞、∞ * 0、0 * ∞、0/0、∞/∞、∞%∞、x%0 、または任意の数値x の場合に発生する可能性がある。また、sqrt(-1) などのいくつかの他の演算でも NaN が生成される可能性がある。「シグナリングNaN」にアクセスした際にも無効な状態が設定され、未初期化変数の使用が示される。これはほとんどの場合、プログラミングエラーを示している。
オーバーフロー
2つの数値の合計がreal.max よりも大きくなるほど大きな数値を加算または乗算した結果、無限大(&infin;)が生成された場合に設定される。これは、ほとんどの場合、結果が不正であることを示しており、是正措置が必要である。
ゼロによる除算
ゼロで割ることによって &plusmn;&infin; が生成された場合に設定される。これは通常、プログラミングエラーを示すが、常にそうなるとは限らない。ゼロで割っても、正しい結果を返す計算の型もある。 (例えば、x == 0 の場合、1/(1+ 1/x) == 0 )。 ゼロに近い非常に小さな数値で割った場合も、無限の結果が生成されるが、divisionByZero フラグではなく、オーバーフローフラグが設定される。
アンダーフロー
これは、2つの数値が減算または除算され、その結果が小数点以下桁数が不足しているために精度が失われた場合に発生する。極端なアンダーフローが発生すると、結果はゼロになる。アンダーフローが問題を引き起こすことはほとんどなく、通常は無視できる。
不正確
これは丸め処理が行われたことを示す。ほとんどの浮動小数点演算でこのフラグが設定される!これは、数値解析の草創期に用いられた難解なテクニックをサポートするためにハードウェアに組み込まれたものと思われる。これは常に無視してよい。

浮動小数点トラップは、上記のカテゴリーのいずれかで有効にすることができる。有効にすると、ハードウェア例外が発生する。 これは、非常に有用なデバッグ支援となる。 より高度な使用方法として、まだどのプラットフォームでもサポートされていないものがあるが(!)、ハードウェア例外ハンドラとして使用するネストされた関数を提供することが挙げられる。これは、オーバーフローおよびアンダーフローの例外に対して最も有用である。

浮動小数点数と「純粋な nothrow」

浮動小数点演算は、どんな些細なものであっても浮動小数点丸め状態の影響を受け、スティッキーフラグに書き込まれる。 ステータスフラグと制御状態は、このように「隠れた変数」であり、潜在的にすべてのpure 関数に影響を与える可能性がある。また、浮動小数点トラップが有効になっている場合、浮動小数点演算はすべてハードウェア例外を生成することができる。 Dは、pure およびnothrow 関数が呼び出された場合でも、限られた状況下で浮動小数点制御モードと例外フラグを使用できるようにする機能を提供する。

[TODO: 私は2つの提案をしたが、まだウォルターを説得できていない!]。

結論

Dは汎用プログラミング言語であり、多くの高度な概念をサポートしているが、最新の浮動小数点ハードウェアのほぼすべての機能に直接かつ便利にアクセスできる。このため、堅牢で高性能な数値コードの開発に最適な言語となっている。また、マシンに対する深い理解を促す言語でもあり、革新や新しいアルゴリズムの開発に適した土壌となっている。

参考文献および関連情報

  1. 「コンピュータ科学者が知っておくべき浮動小数点演算」
  2. 「浮動小数点演算の老大家へのインタビュー:ウィリアム・カーハン氏へのチャールズ・セベランス氏によるインタビュー」
  3. N. Brisebarre、J-M Muller、S.K. Raina、「除数が事前に分かっている場合の正確な丸め浮動小数点除算の高速化」、IEEE Transactions on Computers、第53巻、1069-1072ページ(2004年)。
  4. 「ボルネオの言語」