Improve this page Github へのログインが必要です。 簡単な修正は、ここから fork、オンライン編集、pull request ができます。 大きな修正については、 通常の clone で行って下さい。 Page wiki 関連するWikiページを参照・編集 English このページの英語版（原文）

マシンの中のリアル : Dの浮動小数点数

はじめに

by Don Clugston

コンピュータは元々、数学を扱うために作られた装置でした。最初期のコンピュータはほとんど全ての時間を数式を解くことに費やしています。現在でこそ科学計算はコンピュータの対象分野の小さな一部分に過ぎませんが、ここには過去からの素晴らしい遺産が確かに残っています。現代のほぼ全てのコンピュータは、数学的な計算を正確にそして非常に高速に行うための優れたハードウェアを備えているのです。しかし、残念なことに、ほとんどのプログラミング言語ではプログラマがこのハードウェアの機能をフル活用するのは難しくなっています。そしてさらに大きな問題は、ドキュメントの欠如です。数学的な処理を得意とするプログラマであっても、浮動小数点算術の様々な側面が謎のまま脇に置かれている場合が見受けられます。

システムズプログラミング言語として、D言語は、プログラマとコンパイラ、そしてプログラマとマシンの間の障壁を全て取り除くべく試みられています。この哲学は、特に浮動小数点の扱いにおいて顕著です。ここで私の個人的な逸話をお話しすることで、ハードウェアに関する正確な理解がいかに重要であるか、説明できればと思います。

私の最初の「浮動小数点の悪夢」は、100回に1回かそのくらいの実行でだけハングするC++プログラムでした。問題を追いかけたところ、最終的に、ある while が無限ループしていることを突き止めました。そのコードの重要な部分を抜粋すると、次のようになります: 

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 に触れていたのはインテルのアセンブリ命令セット説明書だけでした! 私のプログラムがダメになったのは、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-bit 拡張' をサポートしていればそれになり、そうでなければ double. と同じになります。将来的には、754-2008 から quadruple, decimal64, decimal128 という新たな型が加わるでしょう。

これらの型の特性は言語組み込みの プロパティ を使って簡単にアクセスできます。例えば、float.max は float が保持できる最大の数を表します。float.mant_dig は仮数部に保持できる桁数（ビット数）です。

D での計算を理解するには、IEEE 浮動小数点計算法の基礎を理解する必要があります。根本的な点として、浮動小数点数は、無限個ある実数をわずかなバイト数に落とし込んだものです。 IEEE 32-bit 浮動小数点ではわずか40億の数しか区別できません。そんな痛ましいほど小さな表現空間であっても、IEEE 浮動小数点数は数学的実数を扱えているんだという幻想を維持することに著しい成功を収めています。とはいえ、その幻想が叩き壊された時のことを理解することも重要です。

ほとんどの問題は、表現可能な数字の分布から引き起こされます。IEEE の数直線は数学的数直線とは全く異なるものです。


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

IEEE 数直線の半分が -1 から 1 の間にあることに注目してください。float で表現可能な数のうち10億個が0から0.5の間にある一方、0.5から1の間には800万個しかありません。これは計算の正確さに対する重大な含みを持ちます。ゼロ付近では実効精度は非常に高くなります。このあといくつかの例においてこの優位性を示すために -1 から 1 までの数を別に扱います。

特別な数にも注目してください。±∞ があります。また、±float.min_normal と 0 の間の数、いわゆる「非正規化数」は、精度を殺して表現されます。さらに、+0と-0という 2種類の ゼロとして表現されます。そして最後は "NaN" (非数)で、意味をなさない値であり、最初の例にあるような厄介事を起こします。

NaNはなぜ存在するのでしょうか? それは価値のある役割を負っています。浮動小数点演算から 未定義のふるまいを根絶する のです。 これは浮動小数点を完全に予測可能なものにします。int 型で 3/0 をしようとすれば、ハードウェアゼロ除算トラップハンドラが呼び出されて、おそらくあなたのプログラムは停止するでしょう。一方、浮動小数点の 3.0/0.0 は ∞ になります。また、数値オーバーフロー (たとえば real.max*2) も ∞ を生みます。 アプリケーションによりますが、∞ は完全に正当な結果であるかもしれませんが、より典型的には、それはエラーを示します。0.0 / 0.0 などの無意味な操作はNaNを生みます。$(しかしあなたのプログラムは制御不能にはなりません。)一見すると、無限とNaNは不要に見えるかもしれません――なぜ整数型のように単にエラー扱いしないのでしょうか？つまるところ、ゼロ除算を避けるのは簡単で、単に除算の前に毎回分母がゼロかチェックするだけです。本当の困難はオーバーフローがもたらすのです。乗算でオーバーフローがで起こるかどうか前もって決定するのは非常に困難でしょう。

非正規化数は、ある種のずれを防ぎ、"x - y == 0 と x == y は等値である"といった重要な関係を保存するのに必要です。

∞ はオーバーフローによって生まれるので、+∞ と -∞ の両方が必要です。また、+0 と -0 の両方が、"x>0 のとき 1/(1/x) > 0 である"といった恒等式を守るために必要です。しかしながら、その他ほとんど場合には、+0と-0の間に違いはありません。

これらの ‘特別な値’ は普通それほど効率的でないことには注意すべきです。たとえばx86マシンでは、NaN、無限大、非正規化数を含む乗算は普通の数の演算の20～50倍も遅くなります。もしあなたの数値演算コードが予想外に遅いのであれば、あなたはうっかりこれらの特別な値をたくさん作っている可能性があります。後に述べる浮動小数点例外トラップを有効にしておけば、これを簡単に検証できます。

マシンがしていることを覆い隠してしまう最大の要因の1つは、二進数とdecimalの間の変換です。結果表示の際 "%a" フォーマットを使うことでこれを取り除くことができます。これはすこぶる有益なデバッグツールであり、浮動小数点アルゴリズム開発時に飛び切り役に立つ補助になるでしょう。 0x1.23Ap+6 のような形式の16進浮動小数点フォーマットも、ソースコード中であなたの入力したデータが正確にあなたの意図したものであることを保障します。

量子化された浮動小数点型の特質

表現可能な値が限られているという事実は、数学的実数ではありえないいくつかの操作を提供することになります。ある数 x が与えられたとき、 nextUp(x) は x より大きい次の表現可能な数を与えます。 nextDown(x) は x より小さい次の表現可能な数を与えます。

数値解析の専門家は、しばしばある種のエラーを記述するのに「最終計算単位/最小桁の単位 ("units in the last place"; ulp)」という、非常に微妙なニュアンスで、しかしぞんざいに扱われがちな用語を使います。 [注: 最も形式的な定義は [J.-M. Muller, "On the definition of ulp(x)",INRIA Technical Report 5504 (2005).] にあります: 実数 x が型 F の隣り合う有限浮動小数点数 a と b の間にあり、かつどちらとも等しくないとき、ulp(x) = abs(b-a) であり、そうでなければ ulp(x) = x*F.epsilon である。加えて ulp(NaN) は NaN であり、ulp(±F.infinity) = ±F.max*F.epsilon である。] 私としては、もっとシンプルな定義が好みです: つまり、「二つの数xとyの間のulp単位での差は、xからyに到達するまでに呼び出す必要があるnextUp()やnextDown()の回数である」。 [注: xやyが浮動小数点数ではない実数ならば、この値は整数ではなくなります。] Dのライブラリ関数 feqrel(x, y) は、x と y の間で共通している上位ビットの個数を数えます。これは精度落ち問題をチェックする簡単な方法です。

浮動小数点数は量子化された値であるという性質は、次のような興味深い結果をもたらします。

条件数 (condition number)

nextUpを使うことで、簡単に条件数の概算値を計算できます。

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

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

この例が示すのは、非常に大きな数 x での1ビットの誤差は、exp(x)になると12ビットの誤差にまで拡大されてしまうことです！ 誤差は次第に大きくなり、realの最大値付近では約6000倍にまで拡大されます。このことを、そのような値xの条件数は6000である、といいます。

float, double, real の意味論

市場を席巻しているx86機では、浮動小数点数演算は伝統的に8087コプロセッサの子孫によって実行されてきました。この "x87" 浮動小数ユニットは、IEEE754 算術を実装した最初のプロセッサです。x86-64 プロセッサでは SSE2 命令セットが代替となっていますが、今日でも、x87 は32bit x86マシンでポータブルに浮動小数演算を実現する手段として残っています (AMD の 32bit プロセッサはSSE2に対応していません)

x87 は、他の多くの浮動小数演算ユニットと比べるとかなり変わっています。80-bit 浮動小数点 だけ をサポートするので "real80" と呼ばれています。double や float の値は演算の前に必ず 80-bit に変換され、80-bit の精度で全ての演算が実行されます。その後で、必要ならば演算結果を 64-bit や 32-bit 精度に縮小します。この実装の意味するところは、最大64-bitの演算子かサポートしていないようなマシンと比べると、精度がかなり高くなる可能性があると言うことです。一方で、可搬性のあるコードを書こうとしたとき、これは枷ともなっています。 (注: x87 は仮数部を double や float と同じに下げる機能は提供していますが、指数部は real80 のままで、非正規化数に対して異なる計算結果になります。厳密に double の計算をシミュレートしようとすると、浮動小数点処理のコードが格段に遅くなってしまいます)

x87 ファミリの他では、Motorola 68K (ColdFire は違います) と Itanium プロセッサが 80-bit 浮動小数をサポートしています。

似たような問題が FMA (fused multiply and accumulate) 命令にも関連してきます。この種の命令をサポートするプロセッサは PowerPC, Itanium, Sparc, Cell, と増えてきています。このようなプロセッサでは、x*y + z のような式を計算する際には x*y の部分は通常の倍の精度で計算されます。本来はフルの精度ロスが起きるような計算であっても、正確に計算されることがあります。 高水準システムズプログラミング言語にとって、全てのプラットフォームで予測可能な動作をするように巧く抽象化しつつ同時にハードウェアの良さを活用するデザイン、というのは大きな挑戦です。

このテーマに関するD言語のアプローチは、以下の事実を踏まえて考えられています:

  1. 全てのプロセッサで完全に同じ動作をさせるのはパフォーマンスコストが物凄く高い。特に、x87 を考えた場合。
  2. 非常に多くのプログラムは特定のプロセッサのみで動かされる。実際には必要ない移植性のために、せっかく存在する高精度演算を犠牲にしてしまうとしたら不幸なことでしょう。
  3. 高精度性と移植性が同時に求められることはない。double の精度では不十分という状況で、特定のプロセッサでだけ精度が増えても仕方が無い。
  4. 言語の機能が特定のプロセッサに強く結びついてしまうのはよくない。

キーとなる設計目標は以下の通りです: 「どんなプロセッサで動かされるとしても、double 型だけをサポートするシステムで動かした時より精度が落ちることがない、そのようなコードが書けるべきである。」

(注: real は、Java言語に対するBorneo提案の ‘indigenous’ に似ています [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. real b = x * y; double r3 = b + z * w;

N最適化によっては (2) と (3) は (1) に変換されてしまうことがありますが、これは実装依存です。 場合 (2) が、余計な丸め誤差を入れてしまっているという意味で特に問題です。

"シンプルな" CPU では、r1==r2==r3 となります。この値を r0 と呼ぶことにしましょう。 PowerPC では、r2==r3 ですが、r1 はFMAを使えるのでそれよりも高精度になることがあります。 x86 では r1==r3 で、この値は r0 よりも高精度になりますが、PowerPC の時ほどの高精度ではありません。 一方で、r2 は r0 よりも精度が下がります。

まとめると、real 型を中間結果の格納に使うことで、double だけをサポートするシンプルなCPUと比べて精度が下がってしまうことはない、と保証されます。

組み込み型のプロパティ

浮動小数点数の基本的なプロパティは epsilon, min_normal, maxです。他の6つの整数型プロパティは、単にこれらのlog2とlog10です。

  float double real80 quadruple decimal64 decimal128
epsilon 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
..max) 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    
max_exp +128 +1024 +16384 +16384    
十進表現に関するプロパティ
dig 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-bit real が使える環境で true になります。 二進型では、dig は、十進で表したときの 最小の 有効十進桁数を表します。表現可能な値が全て完全に表現できるような十進桁数が必要な場合は、さらに2桁伸ばす必要があります。同様に、十進型では、mant_dig 有効な二進桁数の下限を表します。

浮動小数点型 F の値 x と y の有用な関係

加算と減算

乗算と除算

べき乗と対数

NaN ペイロード（追加情報）

IEEE 754 標準によると、NaN の仮数部には ‘ペイロード’ と呼ばれる追加情報を格納できます。ここには、NaN が作られた方法や理由などを含めることができます。歴史を振り返ると、ほとんど全てのプログラミング言語はこの潜在的な可能性を秘めた機能を使っていませんでした。Dでは、このペイロードは正の整数として操作できます。

絶対に NaN の整数ペイロードとしてポインタを格納してはいけません。ガベージコレクタがそのポインタを見つけられなくなってしまいます！

NCEG 比較演算

おなじみの <, >, <=, >= といった比較演算子に加えて、D はさらに "NCEG" 演算子をサポートしています。そのうち４つは普通の比較演算子の否定形を意味しています。他に、<>, <>=, !<>, !<>= が提供されます。この8つの新しい演算子は、NaN が関わってきたときにのみ、普通の演算子と異なる動きをします。従って、多くの場合これらはあまり目立って使われません。主に、計算を始める前にNaNを除く目的で活用されます。 もっとも使いやすい関係式は、おそらく以下のものです:

y がコンパイル時定数 (たとえば 0) ならば、これは !isNaN(x) や isNaN(x) と同じ意味になります。 更に参考までに、x==x も !isNaN(x) と同じ意味になり、x!=x は isNaN(x) を意味します。 abs(x) !< x.infinity は isNaN(x) || isInfinity(x) です。 このような関係式は、コンパイル時関数で使えるという点でまず有用です。 そのほかの NCEG 演算子にはほとんど使い道が知られていません。

IEEE 丸めモード

丸めモードはスコープ単位で制御することが可能です。スコープの終わりでは元の状態が復元されます。 丸めモードには4種類あります。デフォルトは Round to nearest で、統計的にはもっとも正確ですが、もっとも直感的な把握が難しくもあります。ちょうど中間の値は、偶数側に丸められます。

丸めモード rndint(4.5) rndint(5.5) rndint(-4.5) 備考
Round to nearest 4 6 -4 ちょうど中央値ならば偶数側に丸めます
Round down 4 5 -5  
Round up 5 6 -4  
Round to zero 4 5 -4  

丸めモードをあえて変更する理由はほぼありません。 round-up と round-down モードは、区間演算を高速に実装する目的のために特別に作られたモードです。ある種のライブラリにはこの機能がたしかに必要ですが、他で使われる機会はあまりありません。 round-to-zero モードは浮動小数点数を整数にキャストする際に使われます。モード切替は（特にIntelマシンでは）遅いので、最初から round-to-zero モードにしておいて、計算のもっともコアなループで cast(int) と正確に同じ挙動を高速に使う、という使い道はあるかもしれません。

他に丸めモードを切り替える理由としてよく挙がるものはもう一つだけあって、それは、数値計算の安定性の簡単なチェック用途です。丸めモードを変えたら結果が大きく変わってしまったような場合、どこかで大きな丸め誤差が発生してしまっていることの明確な印となります。

IEEE 例外ステータスフラグ

IEEE準拠のプロセッサは必ず、プログラムに伝える必要があるかもしれない"マズい"ことが起こったことを示す特別なステータスビットを持っています。例えば、ieeeFlags.divideByZero はゼロ除算で無限大が作られたことを示します。これらは 'sticky' なフラグです: つまり、一度セットされたら、明示的にクリアされるまでセットされっぱなしになります。言い方を変えると、フラグは計算の最後に一度だけチェックすることで、ほとんどの場合failしないような何千回ものチェックを避けることができます。

以下が、検出できる"マズい"ことの一覧です:

invalid
これは、NaN が出来た時に常にセットされます。∞ - ∞, ∞ * 0, 0 * ∞, 0/0, ∞/∞, ∞%∞, 任意の数xに対するx%0 を計算したときに発生します。その他の演算、たとえば sqrt(-1) なども NaN を生成します。 invalid フラグは、未初期化変数へのアクセスを表す 'signalling NaN' へのアクセス時にも発生します。ほぼ全ての場合、これはプログラミングのミスを示しています。
overflow
∞ が2つの数の加算や乗算の結果生成された場合にセットされます。これは、演算結果が real.max より大きくなるときに発生し、ほとんどの場合、演算結果が正しくなくなっていることを表します。なにか補正する処理が必要です。
divisionByZero
±∞ がゼロ除算によって発生したときにセットされます。これは普通プログラミングのミスを示していますが、そうでないこともあります; ある種の計算では、たとえ途中にゼロ除算が発生しても正しく結果を求めることが出来ます。 (例えば x == 0 のときに 1/(1+ 1/x) == 0 になります)。 注意点として、非常に小さいほとんどゼロの値で割ったときも無限大が結果として返りますが、このときはdivisionByZeroではなくoverflowフラグがセットされます。
underflow
2つの数の減算または除算の結果、結果が非正規化数になって精度が失われたときにセットされます。酷い場合は、結果は完全にゼロになってしまいます。underflow が問題となることは実際にはほとんどなく、単に無視できる場合がほとんどです。
inexact
丸めが発生したことを示します。ほとんど全ての浮動小数点演算がこのフラグをセットします！ どうも、数値解析の初期の頃の超絶トリックをサポートするためにハードウェアに含まれたようですが、今日では常に無視して構いません。

浮動小数点エラーのトラップは、ここに列挙したどのカテゴリについても有効化できます。有効化した場合、ハードウェア例外が生成されるようになります。 これは計り知れないほどデバッグの助けとなります。 より発展的な使用法（まだどのプラットフォームでもサポートされていませんが(!)）としては、ネスト関数をハードウェア例外のハンドラとして提供するといった機能が考えられます。これは特に、overflow と underflow 例外に対して役に立つでしょう。

Floating point and ‘pure nothrow’

どんな浮動小数点演算も、一番トリビアルなものであっても、演算は丸めモードの状態に影響され、sticky flags への書き込みを行います。状態フラグや制御ステートはつまり '隠された変数' であって、全ての pure 関数の動作に影響する可能性を含んでいます; さらに、浮動小数点のトラップが有効にされれば、どんな浮動小数点演算も、ハードウェア例外を投げることができるようになります。 D は、浮動小数点の制御モードや例外フラグを限定された状況で使えるようにする機能を提供しており、これは pure や nothrow 関数が呼ばれている時であっても使用できます。

[TODO: これについて2つ提案を出しましたが、まだWalterを説得できていません！]

まとめ

D は汎用プログラミング言語で多くの高水準な機能をサポートしていますが、一方で、モダンな浮動小数点演算ハードウェアのほぼ全機能に対する直接的で便利なアクセスをも可能としています。これによって、Dは、ロバストで高性能な数値演算コードを書くための優れた言語となっています。また同時に、マシンに対する深い理解を促すような言語デザインは、新しいアルゴリズムの発想の元となる豊かな土壌でもあります。

参考文献

  1. "What Every Computer Scientist Should Know About Floating-Point Arithmetic"
  2. "An Interview with the Old Man of Floating-Point: Reminiscences elicited from William Kahan by Charles Severance"
  3. N. Brisebarre, J-M Muller, and S.K. Raina, "Accelerating Correctly Rounded Floating-Point Division when the Divisor Is Known in Advance", IEEE Trans. on Computers, Vol 53, pp 1069-1072 (2004).
  4. "The Borneo language"