この『数値計算ガイド』の対象読者すべてにとって、David Goldberg 著 1991 年 3 月 Computing Surveys 発行の論説『浮動小数点演算について計算機科学者は何を知っておくべきか』が役立ちます (http://dl.acm.org/citation.cfm?id=103163 を参照)。
この概略は次のとおりです。
浮動小数点演算は、概して難解な問題と捉えられています。コンピュータシステム内部には、浮動小数点が至る所に存在することを考えると、これはむしろ意外なことです。ほとんどの言語で浮動小数点のデータ型が使用されています。PC からスーパーコンピュータまで、コンピュータには浮動小数点アクセラレータが使用されており、浮動小数点アルゴリズムをコンパイルするために、さまざまなコンパイラが随時呼び出されます。事実上、オペレーティングシステムはすべて、オーバーフローなどの浮動小数点例外に対応していると言えます。この論説では、コンピュータシステムの設計者に直接影響を与える浮動小数点の側面について説明します。浮動小数点表現の背景、丸め誤差、続いて IEEE 浮動小数点標準について説明し、最後にコンピュータシステム設計者が浮動小数点を適切にサポートする方法の例を紹介します。
この付録は、発行済みの Goldberg 氏の論説の一部ではありません。特定のポイントを明確にし、読者が論説から推測する可能性のある IEEE 標準に関する誤解を正すために追加されました。この資料は David Goldberg 氏による文書ではありませんが、氏の許可のもと記載しています。754-1985 と 854-1987 標準は、2 進と 10 進の両方の浮動小数点演算を指定した 754-2008 に置き換えられています。これは Goldberg 氏の論説には影響しません。
この付録では特に、IEEE 754 実装間の相違について説明しています。このトピックには、次のサブトピックがあります。
Goldberg 氏の論説では、プログラマがそのプログラムの正確さと精度を浮動小数点演算のプロパティーに依存する場合があるので、浮動小数点演算を慎重に実装する必要があると述べています。特に、IEEE 標準は慎重な実装が必要であり、この標準に従ったシステム上でのみ正しく機能し、正確な結果をもたらす有用なプログラムを作成することが可能です。一部の読者は、このようなプログラムがすべての IEEE 標準に移植できると判断する可能性もあります。実際、「プログラムが 2 つのマシン間で移され、どちらも IEEE 演算をサポートしている場合、中間結果が異なっていても、演算の差ではなく、ソフトウェアバグによるものに違いない」という見解が真であるならば、移植可能なソフトウェアの作成は簡単になります。
しかし、IEEE 標準では、準拠するすべてのシステムで同じプログラムが同一の結果をもたらすことを保証していません。ほとんどのプログラムは実際には、さまざまな理由により、別々のシステムで異なる結果を生成します。ひとつには、多くのプログラムは、システムライブラリに用意されている初等関数を使用しますが、標準はこれらの関数を完全には規定していないからです。1985 標準は、10 進形式と 2 進形式との間の数値変換を完全には規定しておらず、超越関数をまったく規定していませんでした。
ほとんどのプログラマは、IEEE 標準で定められた数値形式と演算だけを使用するプログラムでさえ、別々のシステムでは算出する結果が異なるという可能性を認識していません。実際、標準規格の作成者は、別々の実装では異なる結果を取得できるようにすることを意図していました。その意図は、IEEE 754 標準でのデスティネーションという用語の定義で明らかです。「デスティネーションは、ユーザーによって明示的に指定される場合も、システムによって暗黙的に与えられる場合 (たとえば、部分式の中間結果やプロシージャーの引数など) もあります。一部の言語には、デスティネーションにおける中間計算の結果をユーザーが制御できないようにするものがあります。それでもこの標準は、デスティネーションの形式とオペランドの値を単位として、演算の結果を定義します」(IEEE 754-1985、7 ページ)。つまり、IEEE 標準では、それぞれの結果が配置先のデスティネーションの精度に正しく丸められる必要がありますが、デスティネーションの精度をユーザーのプログラムで判断することを必要としていません。したがって、別々のシステムでは、デスティネーションの結果の精度が異なることがあるため、これらのシステムがすべて標準に準拠していても、同じプログラムで異なる結果が生成され、大幅に異なる結果になる場合もあります。
参照先の論説のいくつかの例は、浮動小数点演算を丸める方法に関する知識に依存しています。このような例を信頼するには、プログラマは、どのようにプログラムが解釈されるか、特に IEEE システムでは、各算術演算のデスティネーションの精度がどのようになるかについて予測できることが必要になります。しかし、IEEE 標準のデスティネーションの定義における盲点によって、どのようにプログラムが解釈されるかを認識するプログラマの能力が損なわれます。結果的に、Goldberg 氏の論説に示す例のいくつかは、高級言語で一見、移植可能なプログラムとして実装された場合、一般に、デスティネーションの結果がプログラマの予想と異なる精度になるような IEEE システムでは正しく機能しない可能性があります。ほかの例には機能するものもありますが、例が機能するかどうかの証明は、平均的なプログラマの能力の範疇外である可能性があります。
この付録では、IEEE 754 演算の既存の実装は、通常は実装で使用するデスティネーション形式の精度に基づいて分類されています。論説から、一部の例では、プログラムの予想より結果の精度が大きい場合、予想された精度が使用された場合には正しくなる結果でも、誤った結果が演算されることを示しています。論説では、想定外の精度によってプログラムが無効にならない場合でも、その精度を処理するために必要な作業について解説するため、1 つの論証が訂正されています。これらの例では、IEEE 標準のさまざまな規定にもかかわらず、この標準規格が異なる実装で許容している相違によって、動作を正確に予測できる移植可能で効率的な数値ソフトウェアの作成が妨げられている可能性を示唆しています。このようなソフトウェアを開発するには、プログラムが依存する浮動小数点演算セマンティックスを表現できるように、IEEE 標準が許容する多様性を制限するプログラミング言語と環境を最初に作成する必要があります。2008 バージョンの 1985 標準では、このようなプログラミング言語に関する推奨事項を定めています。
IEEE 754 演算の現在の実装は、ハードウェアで異なる浮動小数点形式のサポートの程度で区別される 2 つのグループに分けられます。Intel x86 ファミリのプロセッサに例示される拡張ベースシステムは、–xarch=386 オプションでコンパイルされており、拡張倍精度形式を完全にサポートしていますが、単精度と倍精度については一部しかサポートしていません。このシステムには、単精度と倍精度のデータをロードおよび格納する命令が用意されており、これらのデータをその場で拡張倍精度形式に、または拡張倍精度形式から変換します。また、算術演算の結果が、拡張倍精度形式でレジスタに保持されている場合でも、単精度または倍精度に丸められる特殊モード (デフォルトではない) も用意されています。Motorola 68000 シリーズプロセッサは、これらのモードで、単精度または倍精度形式の範囲に結果を丸めます。Intel x86 および互換プロセッサは、単精度または倍精度形式に結果を丸めますが、保持する範囲は拡張倍精度形式と同じです。ほとんどの RISC プロセッサを含む単精度/倍精度システムでは、単精度および倍精度形式を完全にサポートしていますが、IEEE 準拠の拡張倍精度形式はサポートしていません。x86 SSE2 拡張機能では、SSE2 レジスタで単精度/倍精度システムを提供しますが、拡張精度レジスタでは拡張精度を引き続きサポートします。
拡張ベースシステムと単精度/倍精度システムで計算の動作がどのように異なるかを確認するために、次のように、Goldberg 氏の論説の例「システムの諸側面」の C バージョンについて考えてみます。
int main() {
double q;
q = 3.0/7.0;
if (q == 3.0/7.0) printf("Equal\n");
else printf("Not Equal\n");
return 0;
}
定数 3.0 と 7.0 は倍精度浮動小数点数と解釈され、式 3.0/7.0 は double データ型を継承します。単精度/倍精度システムでは、倍精度のほうが効率的に使用できる形式のため、この式は倍精度で評価されます。したがって、q には正しく倍精度に丸められた値 3.0/7.0 が代入されます。次の行で、式 3.0/7.0 はふたたび倍精度で評価され、結果は、q に直前に代入した値に等しくなるので、プログラムは予想どおり「Equal」と出力します。
拡張ベースシステムでは、式 3.0/7.0 の型が double であっても、商はレジスタにおいて拡張倍精度形式で計算されるため、デフォルトモードで拡張倍精度に丸められます。しかし、結果の値が変数 q に代入されるときにとメモリーに格納され、q が double と宣言されているために値は倍精度に丸められます。次の行で、式 3.0/7.0 はふたたび拡張精度で評価されることがあり、q に格納された倍精度値とは異なる結果が生じ、この結果、プログラムは「Not equal」と出力します。
ほかの結果も考えられます。コンパイラは、2 行目の式 3.0/7.0 の値を格納し、丸めたあとで q と比較することも、格納せずに拡張精度で q をレジスタに保持することもできます。最適化コンパイラは、コンパイル時に式 3.0/7.0 を倍精度で、あるいは拡張倍精度で評価することが考えられます。ある x86 コンパイラにおいて、プログラムは、最適化でコンパイルするときには「Equal」を出力し、デバッグ用にコンパイルするときには「Not Equal」を出力します。また、拡張ベースシステムの一部のコンパイラには、丸め精度モードを自動的に変更し、より広い範囲を使用して、レジスタ内に結果を生成する演算によってその結果を単精度または倍精度に丸めるものもあります。そのため、これらのシステムでは、そのソースコードを読み取り、IEEE 754 演算の基本的な解釈を適用するだけでは、プログラムの動作を予測できません。IEEE 754 準拠環境を提供できない原因は、ハードウェアにもコンパイラにもありません。ハードウェアは、要求されているとおり各デスティネーションに丸めた結果を正しく配布しており、コンパイラは、許可されているとおりユーザーの制御の及ばないデスティネーションに中間結果を代入しています。
このセクションの残りでは、Oracle Solaris Studio 12 以前のリリースでデフォルトであった –xarch=386 でコンパイルされた x86 の動作について説明します。Oracle Solaris Studio 12.4 では、デフォルトは –xarch=sse2 です。
世間一般の通念では、拡張ベースシステムは、常に単精度/倍精度システムと少なくとも同じ精度、多くの場合はそれを超える精度を提供しているため、単精度/倍精度システムで配布される結果以上に正確ではないとしても、最低限は正確な結果を生成する必要があるとされています。次のセクションに示す単純な例、およびその例に基づく微妙なプログラムでは、この通念が単なる理想であることを示しています。単精度/倍精度システムには実際に移植できる、一見、移植性のあるプログラムの中には、コンパイラとハードウェアが重なり合って、プログラムの予想以上の精度を提供しようとするため、拡張ベースシステムで正しくない結果をもたらすものがあります。
現在のプログラミング言語では、プログラムが予想する精度を指定することが困難になっています。Goldberg 氏の論説の言語とコンパイラのセクションで述べているように、多くのプログラミング言語では、10.0*x などの式が同じコンテキストで複数出現したときに、それぞれが同じ値に評価される必要があるとは指定しません。この点に関しては、Ada などの一部の言語では、IEEE 標準以前のさまざまな演算のバリエーションの影響を受けています。これより新しい ANSI-C などの言語は、標準に準拠した拡張ベースシステムに影響を受けています。実際、ANSI C 標準では明示的に、コンパイラが浮動小数点式を、通常その型に関連付けられている精度よりも広い精度に評価することを許可しています。結果として、式 10.0*x の値は、次のようなさまざまな要因によって異なることがあります。式がただちに変数に代入されるのか、それともより大きな式のサブエクスプレッションとして表現されるのか、式が比較に関係するかどうか、式が引数として関数に渡されるかどうか、渡される場合、引数は値と参照のどちらとして渡されるのか。現在の精度モード。プログラムのコンパイルに使用された最適化のレベル、プログラムのコンパイル時にコンパイラが使用した精度モードと式の評価方法、などです。
予測のつかない式の評価のすべての原因が言語規格にあるわけではありません。拡張ベースシステムは、可能なかぎり式が拡張精度レジスタで評価されるときにもっとも効率的に実行しますが、格納する必要のある値は、要求されるもっとも短い精度で格納されます。10.0*x がどこでも同じ値に評価されるように言語に要求させると、これらのシステムのパフォーマンス低下を避けられなくなります。これらのシステムで、構文的に同等のコンテキストで 10.0*x を別の値に評価すること許可すると、それ自体のペナルティーが正確な数値ソフトウェアのプログラマに課され、意図した意味を表現するためにプログラムの構文に依存できなくなります。
次の例では、指定された式が常に同じ値に評価されるという仮定に、実際のプログラムが依存しているかどうかを調べます。Goldberg 氏の論説の定理 4 にある、Fortran で記述された ln(1 + x) を計算するアルゴリズムを参照してください。
real function log1p(x) real x if (1.0 + x .eq. 1.0) then log1p = x else log1p = log(1.0 + x) * x / ((1.0 + x) - 1.0) endif return
拡張ベースシステムでは、コンパイラは 3 行目の式 1.0 + x を拡張精度で評価し、その結果を 1.0 と比較します。ただし、同じ式が 6 行目のログ関数に渡される場合、コンパイラはその値をメモリーに格納して、単精度に丸めることができます。したがって、拡張精度では 1.0 + x が 1.0 に丸められるほど x は小さくないが、単精度では 1.0 + x が 1.0 に丸められるほどには小さい場合、log1p(x) で返される値は x ではなくゼロになり、相対誤差は 1 で 5ε よりもわずかに大きくなります。同様に、サブエクスプレッション 1.0 + x が繰り返し出現する 6 行目の式の残りが拡張精度で評価されるとします。この場合、x が小さいが、単精度で 1.0 + x が 1.0 に丸められるほど小さくはない場合、log1p(x) で返される値が正しい値を x と同じくらい上回り、この場合も相対誤差は 1 に近づきます。具体的な例としては、x を 2-24 + 2-47 とすると、x は、1.0 + x が、次に大きな数である 1 + 2-23 に丸められるような最小の単精度数になります。この場合、log(1.0 + x) はおよそ 2-23 になります。6 行目の式の分母は拡張精度で評価されるため、正確に計算され、x になります。したがって、log1p(x) はおよそ 2-23 を返し、これは正確な値のほぼ 2 倍になります。
これは実際に、少なくとも 1 つのコンパイラで発生します。前述のコードが、-O 最適化フラグを使用して x86 システム用の Sun WorkShop Compilers 4.2.1 Fortran 77 コンパイラでコンパイルされる場合、生成されたコードは説明のとおり 1.0 + x を計算します。その結果、この関数は log1p(1.0e-10) にはゼロを、log1p(5.97e-8) には 1.19209E-07 を配布します。
定理 4 のアルゴリズムが正しく動作するには、式 1.0 + x が出現するたびに同じ方法で評価される必要があります。このアルゴリズムは、1.0 + x が、あるときは拡張倍精度に評価され、別のときには単精度または倍精度に評価されるというような場合にのみ、拡張ベースシステムで失敗する可能性があります。log は、Fortran の総称組み込み関数なので、コンパイラは完全に拡張精度で式 1.0 + x を評価し、同じ精度でその対数を計算できますが、明らかに、コンパイラがそうすることは想定できません。ユーザー定義の関数を使用した場合の類似例についても想定できます。この場合、関数が単精度の結果を返したとしても、コンパイラはまだ拡張精度で引数を保持できますが、既存の Fortran コンパイラでこれを実行するものは、存在するとしてもごくわずかです。そのため、1.0 + x を変数に割り当てることによって、これを常に同じように評価するよう試行するかもしれません。しかし、変数 real を宣言した場合でも、ある変数が出現しても、これに拡張精度でレジスタに保持されている値を代入し、別の変数にはメモリーに格納されている単精代変数のある出現の代わりに、拡張精度でレジスタに保持されている値を使用し、別の変数には単精度でメモリーに格納されている値を代入するようなコンパイラによって失敗することがあります。そのために、拡張精度形式に対応する型の変数を宣言する必要があります。標準の FORTRAN 77 にはこの方法が提供されておらず、Fortran 95 にはさまざまな形式を記述するための SELECTED_REAL_KIND メカニズムが用意されていますが、拡張精度で変数を宣言されるように拡張精度で式を評価する実装を明示的には規定していません。つまり、標準 Fortran では、式 1.0 + x が実証を無効にするような方法で評価することを確実に防ぐように、このプログラムを記述する移植可能な方法はありません。
拡張ベースシステムでは、それぞれのサブエクスプレッションが格納され、そのために同じ精度に丸められている場合でも、正しく動作しない可能性がある例がほかにもあります。原因は二重丸めです。デフォルトの精度モードでは、拡張ベースシステムは、最初にそれぞれの結果を拡張倍精度に丸めます。この結果がその後倍精度に格納されると、ふたたび丸められます。これら 2 回の丸めが組み合わさることにより、最初の結果を正しく倍精度に丸めるときに得られるものとは異なる値が生成されることがあります。これは、拡張倍精度に丸められたときの結果が「中間の場合」である場合、つまり、2 つの倍精度数のちょうど中間になるため、2 回目の丸めが round-ties-to-even 規則によって決定される場合に起こる可能性があります。この 2 回目の丸めが 1 回目と同じ方向になる場合、最終的な丸め誤差は最後の 1 位の桁の半分を超えてしまいます。しかし、この二重丸めは、倍精度計算にしか影響しません。2 つの p ビット数の和、差、積、または商、あるいは p ビット数の平方根が、最初 q ビットに丸められ、続いて p ビットに丸められるような場合、q ≥ 2p + 2 であれば、結果は p ビットに 1 回だけ丸められた場合と同じになります。拡張倍精度は十分な広いため、単精度計算は二重丸めを受けることはありません。
正確な丸めを前提とするアルゴリズムの中には、二重丸めに対応しないものもあります。実際、正しい丸めを必要とせず、IEEE 754 に準拠していないさまざまなマシンで正しく機能するアルゴリズムでも、二重丸めで失敗することがあります。この中でもっとも便利なものは、Goldberg 氏の論説の定理 5 に述べられている、シミュレーションされた多倍精度演算を実行する、移植可能なアルゴリズムです。たとえば、定理 6 に述べられている、浮動小数点を上位部分と下位部分に分けるためのプロシージャ―は、二重丸め演算では正しく動作しません。倍精度数 252 + 3 × 226 – 1 を 2 つの部分に、それぞれ最大でも 26 ビットで分割してみます。それぞれの演算が正しく倍精度に丸められたら、上位部分は 252 + 227 で下位部分は 226 – 1 になりますが、それぞれの演算が最初に拡張倍精度に丸められ、続いて倍精度に丸められると、252 + 228 の上位部分と –226 – 1 の下位部分が生成されます。下位部分の数は 27 ビットを必要とするため、その 2 乗は倍精度では正確に計算できません。この数の 2 乗を拡張倍精度で計算することは可能ですが、その結果のアルゴリズムは単精度/倍精度システムには移植できなくなります。また、多倍精度乗算アルゴリズムの後半の手順では、すべての部分積が倍精度で計算されていることを前提としています。倍精度と拡張倍精度の変数の組み合わせを正しく処理するには、非常に実装コストがかかります。
同様に、倍精度数の配列として表される多倍精度数を追加するための移植可能なアルゴリズムは、二重丸め演算で失敗することがあります。これらのアルゴリズムは通常、カハンの総和公式に類似した手法に依存しています。総和公式の非公式な説明として、Goldberg 氏の論説の加算でのエラーのセクションでは次のように提案しています。s と y が |s| ≥ |y| である浮動小数点変数であり、次のように計算すれば、
t = s + y; e = (s - t) + y;
ほとんどの演算で、e が t の計算で生じた丸め誤差を正確に回復させます。ただし、この手法は二重丸め演算では機能しません。s = 252 + 1、y = 1/2 – 2-54 の場合、s + y が最初に拡張倍精度で 252 + 3/2 に丸められ、この値がさらに round-ties-to-even 規則によって倍精度で 252 + 2 に丸められます。したがって、t の計算での最終的な丸め誤差は 1/2 + 2-54 になり、これは倍精度では正しく表現できず、このため上記の式では正しく計算できません。この場合も、拡張倍精度で和を計算することにより丸め誤差を回復することは可能ですが、その場合、プログラムは最終出力を引き下げて倍精度に戻す作業が必要になり、二重丸めがこのプロセスにも影響することがあります。このような理由により、これらの方法で多倍精度演算をシミュレーションする移植可能なプログラムは、さまざまなマシンで正しく効率的に動作しますが、拡張ベースシステムでは公表されているようには動作しません。
最後に、一見正確な丸めに依存するようなアルゴリズムでも、実際は二重の丸めで正しく動作するものがあります。このような場合では、実装ではなく、アルゴリズムが公表どおりに動作するかどうかを検証するために、二重丸めを実行する必要があります。これを説明するため、定理 7 を別の例で証明します。
m と n が IEEE 754 の倍精度で表現可能な |m| < 252 の整数であり、n が特殊な式 n = 2i + 2j である場合、両方の浮動小数点演算が倍精度に正しく丸められるか、最初に拡張倍精度に丸められてから倍精度に丸められるとすると、(m ⊘n) ⊗ n = m になる。
損失なしの m > 0 と仮定します。q = m ⊘ n とします。2 の累乗でスケールする場合、252 ≤ m < 253 で q についても同様であり、その結果、m と q の両方が、その最下位ビットが 1 位の桁を占める整数になる (つまり、ulp(m) = ulp(q) = 1) 同等の設定を考えられます。スケール前に m < 252 と仮定したため、スケール後の m は偶数の整数になります。また、m と q のスケール後の値は m/2 < q < 2m を満たすため、対応する n の値は、m か q のどちらが大きいかに応じて次の 2 つの式のいずれかになります。q < m の場合、明らかに 1 < n < 2 であり、n は 2 つの 2 の累乗の和であるため、特定の k については n = 1 + 2-k であり、同様に q > m の場合は 1/2 < n < 1 であるため、n = 1/2 + 2-(k+1) になります。n は 2 つの 2 の累乗の和であるため、1 にもっとも近い値である n は、n = 1 + 2-52 です。m/(1 + 2-52) は、m より小さい 2 番目に小さな倍精度数より大きくないため、q = m にはなりません)。
q を計算するための丸め誤差を e と仮定し、その結果 q = m/n + e となり、計算された値 q ⊗ n は m + ne の (1 回または 2 回) 丸められた値になるとします。まず、それぞれの浮動小数点演算が倍精度に正しく丸められる場合を考えてみます。この場合、|e| < 1/2 になります。n の式が 1/2 + 2-(k+1) である場合、ne = nq – m は、2-(k+1) の整数の倍数であり、|ne| < 1/4 + 2-(k+1) になります。これは |ne| ≤ 1/4 であることを意味します。m と次に大きな表現可能な数の差は 1 であり、m と次に小さな表現可能な数の差は、m > 252 の場合は 1、m = 252 の場合は 1/2 であることを思い出してください。このため、|ne| ≤ 1/4 の場合、m + ne は m に丸められます。(m = 252 で ne = –1/4 の場合でも、積は、round-ties-to-even 規則により、m に丸められます)。同様に、n の式が 1 + 2-k である場合、ne は、2-k の整数倍であり、|ne| < 1/2 + 2-(k+1) です。これは |ne| ≤ 1/2 を意味します。m は厳密には q より大きいため、この場合、m = 252 にはなりません。そのため、m と隣接した表現可能な数値との差は ∓1 です。したがって、|ne| ≤ 1/2 の場合も、m + ne は m に丸められます。(|ne| = 1/2 の場合でも、m が偶数であるため、積は round-ties-to-even 規則で m に丸められます)。正確に丸められる演算の証明はこれで終了です。
二重丸め演算でも、q は (実際に 2 回丸められた場合でも) 正確に丸められた商になることがあるため、上記と同じく |e| < 1/2 です。この場合、q ⊗ n が 2 回丸められるという事実を考慮すると、前のパラグラフの引数に訴えることができます。このことを説明するために、IEEE 標準では拡張倍精度形式が少なくとも 64 の有効ビットを含むことを規定し、そのため、m ∓ 1/2 と m ∓ 1/4 の数は拡張倍精度で正確に表現できることに注意してください。したがって、n の式が 1/2 + 2-(k+1) であり、その結果 |ne| ≤ 1/4 になる場合、m + ne を拡張倍精度に丸めると、m と最大でも 1/4 の差がある結果が生成され、上記のように、この値は倍精度で m に丸められます。同様に、n の式が 1 + 2-k であり、その結果 |ne| ≤ 1/2 になる場合、m + ne を拡張倍精度に丸めた結果は、m と最大 1/2 の差がある結果が生成され、その値は倍精度では m に丸められます。この場合 m > 252 であることに注意してください。
最後に、二重丸めのために q が正しく丸められた商になっていない場合を考えてみます。このような場合、最悪の場合には |e| < 1/2 + 2-(d + 1) となり、ここでは d は拡張倍精度形式での余分なビット数です。既存の拡張ベースシステムはすべて、ちょうど 64 個の有効ビットを持つ拡張倍精度形式をサポートします。この形式では d = 64 – 53 = 11 になります。二重丸めは、2 回目の丸めが round-ties-to-even 規則で決定された場合にかぎって正しくない丸め結果を生成するため、q は偶数の整数である必要があります。したがって、n の式が 1/2 + 2-(k + 1) である場合は、ne = nq – m は 2-k の整数倍になり、|ne| < (1/2 + 2-(k + 1))(1/2 + 2-(d + 1)) = 1/4 + 2-(k + 2) + 2-(d + 2) + 2-(k + d + 2) になります。
k ≤ d の場合、これは |ne| ≤ 1/4 を意味します。k > d の場合、|ne| ≤ 1/4 + 2-(d + 2) になります。どちらの場合でも、積の最初の丸めにより、m と最大 1/4 の差のある結果が得られ、前に示した引数によって、2 回目の丸めで m に丸めます。同様に、n の式が 1 + 2-k である場合、ne は 2-(k - 1) の整数倍になり、|ne| < 1/2 + 2-(k + 1) + 2-(d + 1) + 2-(k + d + 1) になります。
k ≤ d の場合、これは |ne| ≤ 1/2 を意味します。k > d の場合、|ne| ≤ 1/4 + 2-(d + 1) になります。どちらの場合も、積の最初の丸めにより、m と最大 1/2 の差のある結果が得られ、再度、前に示した引数によって、2 回目の丸めで m に丸めます。
前述の証明では、商が二重丸めを引き起こす場合のみ積も二重丸めを引き起こし、その場合でも正しい結果に丸められることを示しています。この証明では、二重丸めの可能性を含めるように推論を拡張することは、浮動小数点演算が 2 つしかないプログラムでも困難になることも示しています。さらに複雑なプログラムでは、二重の丸めの影響について体系的に説明することは不可能であり、倍精度と拡張倍精度の演算の一般的な組み合わせは言うまでもありません。
前の例では、拡張精度自体が有害であることを示しているわけではありません。プログラマが拡張精度の使用を選択できる場合には、多くのプログラムが拡張精度のメリットを得ることができます。しかし、現在のプログラミング言語では、プログラマが拡張精度をいつどのような方法で使用するかを指定する十分な手段が用意されていません。どのようなサポートが必要かを示すため、拡張精度の使用を管理できる方法について検討してみます。
名目上の作業精度として倍精度を使用する移植可能なプログラムでは、より広い精度の使用を制御する 5 つの方法があります。
拡張ベースシステムでは、可能なかぎり拡張精度を使用して最速のコードを生成するようコンパイルする。ほとんどの数値ソフトウェアでは、「計算機イプシロン」が範囲を指定する各演算の相対誤差以上の演算は必要とされません。入力データは入力時に丸められていると想定され、その結果が格納されるときにも同様に丸められるため、メモリー内のデータが倍精度で格納されると、計算機イプシロンがその精度の最大の相対丸め誤差と見なされます。したがって、拡張精度の中間結果の計算では、より正確な結果を生成することもありますが、拡張精度は必須ではありません。この場合、感知できるほどプログラムの速度が低下しない場合にのみコンパイラで拡張精度を使用し、それ以外の場合は倍精度を使用することをお勧めします。
かなり高速で十分に広い場合は倍精度よりも広い形式を使用し、それ以外の場合はほかの形式を使用する。一部の計算では、拡張精度を使用できればより簡単に実行できるものもありますが、倍精度でも、いくらか取り組むだけで実行できるようになります。倍精度数のベクトルのユークリッドノルムの計算を考えてみます。要素の 2 乗を計算し、IEEE 754 拡張倍精度形式でその広い指数範囲を使用してそれらの和を累積すると、実際の長さのベクトルに対して早まったアンダーフローまたはオーバーフローを簡単に回避できます。拡張ベースシステムでは、これがノルムを計算するもっとも速い方法です。単精度/倍精度システムでは、拡張倍精度形式がサポートされている場合はソフトウェア内でエミュレートする必要がありますが、これは例外フラグをテストしてアンダーフローまたはオーバーフローが発生していないか確認し、発生した場合は明示的なスケールで計算を繰り返すため、単に倍精度を使用する場合より大幅に速度が低下します。言語では、拡張精度のこのような使用をサポートするため、使用する方法をプログラムが選択できるように適度に高速で利用可能なうちもっとも広い形式を示すとともに、その形式が十分広いこと (たとえば倍精度よりも広い範囲であること) をプログラムが検証できるようにする、それぞれの形式の精度と範囲を示す環境パラメータを提供する必要があります。
ソフトウェアでエミュレートする必要がある場合でも、倍精度より広い形式を使用する。ユークリッドノルムの例よりも複雑なプログラムの場合、プログラマは、2 つのバージョンのプログラムを作成する必要性を避け、代わりに、速度が遅くても拡張精度に頼ることがあります。この場合も、言語には、利用可能なうちもっとも範囲と精度の広い形式をプログラムが判断できるように、環境パラメータが用意されている必要があります。
倍精度より広い精度を使用しない (範囲が拡張された場合でも、結果は倍精度形式の精度に正しく丸められます)。前述のいくつかの例のように、正しく丸められる倍精度演算に依存するように非常に簡単に作成されるプログラムの場合、言語は、中間結果が倍精度より広い指数範囲のレジスタで計算できるとしても、拡張精度を使用しないように指示する方法をプログラマに提供する必要があります。この方法で計算された中間結果でも、メモリーに格納されるときにアンダーフローする場合は二重丸めを引き起こすことがあります。算術演算の結果が最初に 53 の有効ビットに丸められ、次に非正規化する必要があるときにそれより少ない有効ビットにふたたび丸められる場合、最終結果は、非正規化数に一度だけ丸めた場合に取得できる結果とは異なる場合があります。当然、このような二重丸めが、実際のプログラムに悪影響を及ぼすことはほとんどありません。
倍精度形式の精度と範囲の両方に正しく結果を丸める。この倍精度の厳密な強制は、倍精度形式の範囲と精度の両方の限度の近くで、数値ソフトウェアか演算自体をテストするプログラムにもっとも役立ちます。このような慎重を要するテストプログラムは、移植可能な方法で作成することは通常困難であり、特定の形式に結果を強制的に丸めるためにダミーサブルーチンやほかのトリックを使用する必要があるときにさらに難しく、エラーが起こりやすくなります。したがって、拡張ベースシステムを使用して、すべての IEEE 754 実装に移植可能である必要のある堅固なソフトウェアを開発するプログラマはすぐに、膨大な労力を必要とせず、単精度/倍精度システムの演算をエミュレートできることを高く評価するようになります。
現在、これら 5 つのオプションすべてをサポートしている言語はありません。実際、プログラマが拡張精度の使用を制御できる機能を提供した言語はほとんどありません。著名な例外の 1 つが、C 言語のメジャーリビジョンである ISO/IEC 9899:1999 プログラミング言語 - C 標準です。
C99 標準では、通常その型に関連付けられているよりも広い形式で式を評価する実装が可能ですが、式を評価する 3 つの方法のうち、どれか 1 つを使用することを推奨しています。3 つの推奨方法は、式がより広い形式に「変換」される程度によって区別され、実装では、プリプロセッサマクロ FLT_EVAL_METHOD を定義することによって使用する方法を特定することが推奨されています。FLT_EVAL_METHOD が 0 の場合、それぞれの式は、その型に対応する形式で評価されます。FLT_EVAL_METHOD が 1 の場合、float の式は、double に対応する形式に拡張されます。FLT_EVAL_METHOD が 2 の場合、float および double の式は long double に対応する形式に拡張されます。(実装では、式の評価方法が確定できないことを示すために FLT_EVAL_METHOD を –1 に設定することができます)。C99 標準では、<math.h> ヘッダーファイルで型 float_t と double_t を定義することも求められます。float と double は少なくとも同じ広さであり、それぞれ float および double の式を評価するために使用される型に一致するようになっています。たとえば、FLT_EVAL_METHOD が 2 の場合、float_t と double_t はどちらも long double です。最後に、C99 標準では、<float.h> ヘッダーファイルで、それぞれの浮動小数点型に対応する形式の範囲と精度を指定するプリプロセッサマクロを定義することが求められています。
C99 標準で求められる、または推奨されている機能の組み合わせは、上記の 5 つのオプションのすべてではなく一部をサポートしています。たとえば、実装が、long double 型を拡張倍精度形式にマッピングし、FLT_EVAL_METHOD を 2 と定義している場合、プログラマは、拡張精度が比較的高速であり、そのためユークリッドノルムの例のようなプログラムが、型 long double (または double_t) の中間変数を使用できると想定できます。一方、この同じ実装は、無名の式がメモリーに格納されるときでも (たとえば、コンパイラが浮動小数点レジスタをあふれさせる必要があるとき)、この式を拡張精度で保持する必要があり、その結果は、double に宣言された変数に割り当てた式の結果を格納し、倍精度に変換する必要があります。これは、式のレジスタに保持できた場合でも同様です。したがって、double と double_t のどちらの型も、現在の拡張ベースのハードウェアで最速のコードを生成するようにコンパイルすることはできません。
同様に、このセクションの例が示す問題は、すべてではなく、一部を C99 標準が提供する方法で解決できます。C99 標準バージョンの log1p 関数は、式 1.0 + x が (任意の型の) 変数に代入され、その変数が常に使用される場合に正しく動作することが保証されています。しかし、倍精度数を上位部分と下位部分に分割するための移植可能で効率的な C99 標準プログラムは、さらに困難です。double の式が倍精度に正しく丸められることを保証できない場合に、どうすれば正しい位置で分割し、二重丸めを回避できるでしょうか。1 つの解決方法として、double_t 型を使用して、単精度/倍精度システムでは倍精度に、拡張ベースシステムでは拡張精度に分割して、どちらの場合でも演算が正しく丸められるようできます。定理 14 では、ベースとなる演算の精度がわかっているならば、どのビット位置でも分割できるとしており、この情報は FLT_EVAL_METHOD と環境パラメータマクロで取得できます。
次のフラグメントは可能な実装例の 1 つを示しています。
#include <math.h>
#include <float.h>
#if (FLT_EVAL_METHOD==2)
#define PWR2 LDBL_MANT_DIG - (DBL_MANT_DIG/2)
#elif ((FLT_EVAL_METHOD==1) || (FLT_EVAL_METHOD==0))
#define PWR2 DBL_MANT_DIG - (DBL_MANT_DIG/2)
#else
#error FLT_EVAL_METHOD unknown!
#endif
...
double x, xh, xl;
double_t m;
m = scalbn(1.0, PWR2) + 1.0; // 2**PWR2 + 1
xh = (m * x) - ((m * x) - x);
xl = x - xh;
この解決方法を見つけるには、double の式が拡張精度で評価されることがあること、その後の二重丸め問題によってアルゴリズムが誤って動作する可能性があること、および定理 14 に従って拡張精度を代用できるということを認識している必要があります。よりわかりやすい解決方法として、それぞれの式が正しく倍精度に丸められるように指定することもできます。拡張ベースシステムでは、これには丸め精度モードの変更が必要なだけですが、C99 標準には、これを行う移植可能な方法は用意されていません。浮動小数点をサポートするために C90 標準に加える変更点を指定した作業文書である、浮動小数点 C 編集の初期ドラフトでは、丸め精度モードを使用したシステムでの実装で、丸め方向を取得し設定する fegetround および fesetround 関数に類似した、丸め精度を取得し設定する fegetprec および fesetprec 関数を提供することを推奨していました。この推奨は、C99 標準に変更が加えられる前に削除されました。
また、整数演算機能が異なるシステム間の移植可能性をサポートする C99 標準のアプローチは、異なる浮動小数点アーキテクチャーをサポートするより適切な方法を提示しています。C99 標準のそれぞれの実装には、実装がサポートする整数型を定義する <stdint.h> ヘッダーファイルが用意されており、そのサイズと効率性に応じて名前が付けられています。たとえば、int32_t は、ちょうど 32 ビット幅の整数型であり、int_fast16_t は、16 ビット幅以上の実装で最速の整数型であり、intmax_t はサポートされるもっとも広い整数型です。浮動小数点の型についても、類似のスキームを想像できます。たとえば、float53_t は、精度がちょうど 53 ビットだが、範囲がこれを超える可能性のある浮動小数点、float_fast24_t は、精度が 24 ビット以上である実装の最速の型、floatmax_t は、サポートするもっとも広い適度に高速な型の名前になります。高速型では、レジスタがあふれた結果として名前付きの変数の値が変更されてはならないという制約がありますが、拡張ベースシステム上のコンパイラでは、可能なかぎり最速のコードを生成できます。幅がちょうどの型では、拡張ベースシステム上のコンパイラは、前述の制限を条件として、丸め精度モードを指定の精度に丸めるように設定でき、より広い範囲を許可します。double_t は、IEEE 754 の倍精度形式の精度と範囲の両方を持つ型の名前として指定でき、厳密な倍精度評価を実現します。このようなスキームを、適宜名前を付けられた環境パラメータマクロとともに使用すると、前述の 5 つのオプションすべてを簡単にサポートし、プログラムで必要とされる浮動小数点演算セマンティックスをプログラマが簡単かつ明確に指定できるようになります。
拡張精度の言語サポートはそれほど複雑ででしょうか。単精度/倍精度システムでは、上記の 5 つのオプションのうち 4 つが該当し、高速な方と正確な幅の型を区別する必要はありません。ただし、拡張ベースシステムでは難しい選択を迫られます。拡張ベースシステムでは、純粋な倍精度計算、純粋な拡張精度計算のどちらも、2 つを組み合わせた場合ほど効率的にはサポートされず、プログラムに応じて異なる組み合わせが必要になります。また、拡張精度をいつ使用するかの選択は、コンパイラライターには任せないでください。彼らはベンチマークの結果、浮動小数点演算は「本質的に正確でなく」、整数演算には適しておらず、その予測能力もないと判断される傾向があり、また、数値アナリストが公然と語ることもあります。この選択はプログラマに委ねられる必要があり、プログラマは、自身の選択を表すことができる言語が必要になります。
前述の見解は、拡張ベースシステムを軽んじるためのものではなく、いくつかの誤った議論を明らかにすることを目的としており、その議論の最たるものは、すべての IEEE 754 システムが同じプログラムで同一の結果をもたらす必要があるというものでした。これまで拡張ベースシステムと単精度/倍精度システムとの違いに焦点を当てていましたが、これらの各系統内のシステム間にはさらに相違点があります。たとえば、一部の単精度/倍精度システムには、2 つの数を掛け合わせ、3 番目の数を加算し、最後に一度だけ丸めるという単一の命令があります。この演算は、融合型積和演算と呼ばれ、同じプログラムが単精度/倍精度システムごとに異なる結果を生成する場合があり、拡張精度のように、プログラムが使用されるかどうか、またいつ使用されるか応じて、同じシステム上でも同じプログラムが異なる結果を生成する場合もあります。融合型積和演算は、移植できない方法で分割せずに多倍精度乗算を実行する場合に使用できますが、定理 6 の分割プロセスが失敗する場合もあります。IEEE 標準でこのような演算を予期していなかったとしても、中間の積はユーザーの制御の及ばない、正確に保持できる十分な広さを持つ「デスティネーション」に配布され、最終的な和はその単精度または倍精度のデスティネーションに適合するように正しく丸められます。
それでもなお、IEEE 754 が一定のプログラムが配布する必要のある結果を正確に規定しているという考えは、非常に有用です。多くのプログラマは、プログラムをコンパイルするコンパイラやそれを実行するコンピュータに関係なく、プログラムの動作を理解し、それが正しく動作することを証明できると思い込みがちです。この見解をサポートすることは、コンピュータシステムやプログラミング言語の設計者にとっては、多くの点で意味のある目標になります。しかし、浮動小数点演算に関して言えば、この目標を達成することは事実上不可能です。IEEE 標準の作成者もこのことを理解しており、この実現を試みていません。その結果、コンピュータ業界のほぼ全体が IEEE 754 標準に準拠しているにもかかわらず、移植可能なソフトウェアのプログラマは、予測不可能な浮動小数点演算に取り組み続けなければなりません。
プログラマが IEEE 754 の機能を利用するのであれば、浮動小数点演算を予測可能にするプログラミング言語が必要になります。C99 標準では、プログラマは、FLT_EVAL_METHOD ごとに 1 つずつ、プログラムのバージョンを作成する必要がありますが、それなりに予測可能性が向上しています。将来の言語では、IEEE 754 セマンティックスに依存している範囲を明らかにする構文によって、単一のプログラムを作成できるかどうかはまだ明らかではありません。既存の拡張ベースシステムは、所定のシステムでどのように演算を実行するかについては、プログラマよりもコンパイラとハードウェアのほうがより認識できると想定させることによって、その展望を危うくしています。その想定が 2 番目の誤りであり、計算結果に求められる精度は、結果を生成したマシンではなく、そこから引き出される結論にのみ依存し、これらの結論が何であるかを知ることができるのは、プログラマ、コンパイラ、ハードウェアのうちせいぜいプログラマだけなのです。