数値計算ガイド

付録 D

浮動小数点演算について

---

注 - この付録は、1991 年 3 月発行の "Computing Surveys" に掲載された "Every Computer Scientist Should Know About Floating-Point Arithmetic" 稿 (David Goldberg著) を再編集し、著作権を有する Association for Computing Machinery 社 (Copyright 1991) の許可のもとに、印刷しなおしたものです。

---

概要

浮動小数点演算は、概して難解な問題として受け取られることがあります。コンピュータシステムの内部には、浮動小数点演算がいたるところに存在する点を考えれば、この認識には多少、現実とのずれがあるようです。たとえば、ほとんどのコンピュータ言語に、浮動小数点のデータ型が使用されています。また PC からスーパーコンピュータに至るまであらゆるコンピュータで浮動小数点アクセラレータが使用されています。さらに浮動小数点アルゴリズムをコンパイルするために、コンパイラが呼び出されることがよくあります。現実には、オーバーフローなどの浮動小数点例外に対応していないオペレーティングシステムなど皆無であるとも言えます。この付録では、コンピュータシステムの設計担当者に対して直接的な影響を与える浮動小数点のさまざまな側面について説明します。具体的には、浮動小数点表現の背景、丸め誤差、IEEE 浮動小数点標準について説明します。最後にコンピュータ設計者が浮動小数点をどのようにサポートできるかを示す例をいくつか紹介します。

カテゴリ/サブジェクト記述子 (一次 ) C.0 [コンピュータシステムの編成] ：一般 － 命令セットの設計；D.3.4 [プログラミング言語] ：プロセッサ － コンパイラ、最適化；G.1.0 [数値解析] ：一般 － コンピュータ演算、エラー解析、数値アルゴリズム (二次 )

D.2.1 [ソフトウェアエンジニアリング] ：要件/仕様 － 言語；D.3.4 [プログラミング言語] ：正式な定義と理論 － セマンティクス；D.4.1 [オペレーティングシステム] ：プロセス管理 － 同期化

一般用語：アルゴリズム、設計、言語

その他のキーワードと表現：非正規化数、例外、浮動小数点、浮動小数点標準、段階的アンダーフロー、保護桁、NaN、オーバーフロー、相対誤差、丸め誤差、丸めモード、ulp、アンダーフロー

はじめに

コンピュータシステムの設計作業を進める上で、浮動小数点演算に関する情報が必要になることがよくあります。しかし現実には、浮動小数点演算に関する詳しい資料は意外なほど少ないものです。この問題を扱った数少ない書籍の1つである "Floating-Point Computation" (Pat Sterbenz 著) も、すでに廃刊になっています。この付録では、コンピュータシステムの設計に大きくかかわる浮動小数点演算の各側面について詳しく説明します。この付録は大きく 3 部に分けられます。「丸め誤差」の部では、四則演算の基本的な操作に各種の丸めモードを適用した場合の影響について説明します。また ulp と相対誤差という 2 種類の丸め誤差の測定方法に関する背景情報を示します。2 部では、多くのハードウェアベンダー各社によって急速な勢いで採用されつつある IEEE 浮動小数点標準について説明します。IEEE 標準には、基本的な演算に関する丸めの方法が定義されています。IEEE 標準については、「丸め誤差」に示す事項にもとづいて説明します。具体的には命令セットの設計、最適化コンパイラ、および例外処理について説明します。

この付録で浮動小数点について説明するときは、必ず、その説明が妥当であるという正当な理由も併せて示してあります。これは、その理由に関する説明自体が、単純な計算だけでは済ますことができない場合がよくあるからです。また説明の主旨と直接的には関係のない内容については「詳細」という項に独立して示してあります。「詳細」は読み飛ばしてもかまいません。また、この付録には特に、定理に関する「証明」が数多く紹介してあります。「証明」の終わりの箇所は■の記号を付けて表わします。また「証明」を省略した場合は、■の記号を定理の文の直後に示します。

丸め誤差

多数の実数を有限個のビットの中へ無限に圧縮していくには、近似値の表現が必要になります。整数には無限の数がありますが、プログラム中では通常、整数演算の結果は 32 ビットごとに格納されます。一方、任意の固定数ビットで、実数を伴う演算を行う場合は、どれほど多くのビットを使用しても、結果の値を正確には表現できなくなってきます。したがって、浮動小数点演算の結果に、丸め操作を加えることによって、有限の範囲内に収めて表現しなければなりません。この丸め操作のときに生じる誤差は、浮動小数点演算から切り離すことのできない特性と言えます。丸め誤差の測定方法については、「相対誤差と ulp」を参照してください。

そもそも、通常の浮動小数点演算で丸め誤差の問題を避けることができないのだとすれば、基本的な演算処理で多少の丸め誤差が出ることに何か問題があるのでしょうか。この根本的な疑問こそ、この項で扱う主要な問題にほかなりません。「保護桁」では保護桁について説明しています。保護桁を確保することにより、2 つの近い値を減算するときの誤差を低減することができます。保護桁は本来、IBM が 1968 年に、System/360 アーキテクチャの倍精度フォーマットとして保護用の桁を追加した (単精度フォーマットにはすでに保護桁が設定されていました) ことに対応して、既存の全マシンを改訂したことがきっかけで重要視されるようになりました。保護桁の効用を示す例を 2 つ紹介します。

IEEE 標準の仕様は、単なる保護桁の必要性を要求するだけのものではありません。四則演算のほか、平方根計算のアルゴリズムが用意され、各アルゴリズムと同じ結果を生成できる実装が必要になります。したがって、あるマシンから別のマシンにプログラムを移植した場合でも、両方のマシンでが IEEE 標準がサポートされていれば、基本演算により各ビットごとに、まったく同じ結果が生成されることになります。「正確な丸め操作」には、この正確な仕様にもとづいた別の例を示します。

浮動小数点フォーマット

これまでにも、実数を表現するさまざまな形式が提案されてきました。このうち最も一般的な形式は浮動小数点表現 ¹です。浮動小数点表現では、基数 β (つねに偶数とします) と精度 p を使用します。β = 10 で、p = 3 とすると、0.1 という値は1.00 × 10^-1 と表わします。また β = 2 で、p = 24 の場合、小数点数 0.1 は正確に表現することはできませんが、およそ 1.10011001100110011001101 × 2^-4 となります。一般に、浮動小数点数は ± d.dd... d × β^e と表わされます (ここで d.dd... d は「有意桁」²と呼ばれ、桁数は p となります)。さらに正確に言うと、� d₀ . d₁ d₂ � d_p-1 × β^e は、次の数を表わします。

(1)

浮動小数点数という用語は、上記のような形式で表現できる実数のことを指します。浮動小数点表現では、最大許容指数 e_max と最小許容指数 e_min という 2 つのパラメータを使用します。有意桁がβp 個で、許容指数が e_max - e_min + 1 とすれば、浮動小数点数は次のようにコーディングできます。

最後の + 1 は符号ビットを表わします。ここでコーディングそのものは、重要ではありません。実数が浮動小数点数としては正確に表現されない理由が 2 つあります。最も一般的なケースは小数点数 0.1 の例です。これは有限の 10 進小数点数として表記されますが、2 進数では無限の反復表現になります。したがって、β= 2 の場合、0.1 という数値は正確には 2 つの浮動小数点数の間に存在しており、いずれの値でも正しくは表現できないことになります。

また、別の状況として実数が有効範囲を超える場合、すなわち絶対値がβ × 以上、あるいは 1.0 × 以下になっているケースが考られます。この付録では、最初の問題に限って扱います。なお、有効範囲を超える数値については、「無限大」、および 「非正規化数」の各項で説明します。

浮動小数点表現は、必ずしも一意の表現にはなりません。
例えば、0.01 × 10¹ と 1.00 × 10^-1 はいずれも 0.1 を表わします。先行桁がゼロ以外 (上記の公式 (1) で d₀ ≠ 0 の場合) であれば、表現は「正規化されている」と言います。浮動小数点数 1.00 × 10^-1 は正規化されますが、0.01 × 10¹ は正規化されません。β=2、p = 3、e_min = -1、e_max = 2 の場合、図 D-1 に示すとおり、計 16 個の正規化浮動小数点数が存在します。図の太線は有意桁が 1.00 となる値を表わします。浮動小数点表現を正規化すると、一意の表現が得られます。しかし、この制限によってゼロを表現することができなくなります。ゼロを表現する一般的な方法として、1.0 × を使用します。これにより、非負の実数を表わす番号が浮動小数点表現の辞書編成順に対応しているという事実が保持されます³。指数を k ビットフィールドに格納すると、1 ビット分は 0 を表わすために予約する必要があるので、2^k - 1 の値だけが指数用として確保されることになります。

なお、浮動小数点数の × は表記の一部として使用されるもので、浮動小数点の乗算記号とは異なります。× 記号の意味は、文脈から簡単判断することができます。
例えば、(2.5 × 10^-3) × (4.0 × 10²) という式では、浮動小数点の乗算を 1 回だけ行うという意味になります。

図 D-1 正規化数 β= 2、p = 3、e_min = -1、e_max = 2

相対誤差と ulp

浮動小数点演算では本来、丸め誤差を避けることができないので、この誤差を測定する方法を確立しておくことが重要です。仮に β=10、p = 3という浮動小数点フォーマットの例を考えてみます (この付録では、同じ例を使用します)。浮動小数点演算の結果が 3.12 × 10^-2 の場合で、解を無限の精度に対して計算すると 0.0314 になり、最終桁に 2 ユニットの誤差があることがわかります。同じように、実数 0.0314159 を 3.14 × 10^-2 として表わすと、最終桁の誤差は 0.159 となります。一般に浮動小数点数 d.d...d × β^e を使用して z を表わした場合、最終桁の誤差は｜d.d...d - (z/β^e)｜β^p-1 ユニットとなります^{4 5}。ulpという用語は、"units in the last place" (最終桁のユニット) を略したものです。計算の結果が、正確な値にメモリーの近い浮動小数点値となる場合、0.5 ulpも誤差があることになります。浮動小数点値と実数の値の誤差を表わす場合に、相対誤差の近似値を求めることもできます。相対誤差とは、浮動小数点数と実数の差を実数値で割った結果のことです。例えば、3.14159 を 3.14 × 10⁰ によって近似表現した場合の相対誤差は 0.00159/3.14159=0.0005 となります。

0.5 ulp に相当する相対誤差を計算する場合、実際の値に最も近い値を d.dd...dd × β^e で表わすと、0.00...00β′× β^e もの誤差が出ることがあります (ここで β′は β/2 の桁を表わし、浮動小数点数の有意桁の誤差が p ユニット、また誤差の有意桁 p ユニットはゼロになります)。この誤差は ((β/2)β^-p) × β^e です。d.dd...dd × β^e の形式の数の絶対誤差はすべて同じになる一方、値はすべてβ^e から β × β^e の範囲にあるので、相対誤差は ((β/2)β^-p) × β^e/β^e と ((β/2)β^-p) × β^e/β^e+1 の範囲に存在します。

(2)

特に、0.5 ulp に相当する相対誤差は、β を係数として可変となります。この係数のことを「ワブル」と言います。ε= (β/2)β^-p を上式 (2) の上限に設定した場合、実数値が最も近い浮動小数点数に丸められると、相対誤差は必ず ε によって境界が設定されることになります。このことを「マシン・イプシロン」と言います。

上の例で、相対誤差は 0.00159/3.14159 = 0.0005 となります。このような小さい数を避けるために、相対誤差は通常、係数にεを乗じた数として表わします
(この場合、ε = (β/2)β^-p= 5(10)^-3 = 0.005 となります) 。したがって、相対誤差は
(0.00159/3.14159)/0.005) ε= 0.1ε と表現されます。

ulp と相対誤差の違いを説明するために、実値 x = 12.35 を例に取ります。この値は
= 1.24 × 10¹ により近似値として表わすことができます。誤差は 0.5 ulp で、相対誤差は 0.8εです。次に 8 を計算します。正確な値は 8 x = 98.8 で、計算値は
8 = 9.92 × 10¹ です。誤差は 4.0 ulp ですが、相対誤差は依然として 0.8εです。すなわち相対誤差は変わらない一方、ulp 誤差の測定値は 8 倍大きくなっています。一般に、基数が β の場合、ulp による相対誤差は最大 β の係数分だけ変動する可能性があります。逆に、上式 (2) に示すとおり、0.5 ulp という固定誤差によりβ 分だけ変動する相対誤差が生じることになります。

丸め誤差を表わす最も自然な方法は ulp を使用する形式です。例えば、最も近い浮動小数点数に丸めるという操作は、誤差が 0.5 ulp 以下になるということです。ただし、各種の公式が原因で起きる丸め誤差を分析する場合は、相対誤差の方が正確です。これについては 226ページで詳しく分析します。εは、最も近い浮動小数点値より、ワブル係数β 分だけ多くなる可能性があるので、公式による誤差の計算は、β の小さいマシンではそれだけ緊密になります。

丸め誤差の大きさだけが問題の場合、ulpとε の差は多くてもβ 以内に抑えられるので、入れ換えて使用してもかまいません。例えば、浮動小数点数に n ulp だけ誤差がある場合、汚染桁の数が logβ n となるという意味です。計算上の相対誤差が nε の場合は、次のように表わされます。

contaminated digits log_β n (3)

保護桁

2 つの浮動小数点数の差を計算する 1 つの方法として、差そのものの値を正確に求めて、これを最も近い値に丸めることができます。この方法は、オペランドのサイズが大きくなると、その分、コストがかかるようになります。p = 3 とすると 2.15 × 10¹² - 1.25 × 10^-5 は次のように計算できます。

x = 2.15 × 10¹²
y = 0.0000000000000000125 × 10¹²
x - y = 2.1499999999999999875 × 10¹²

これにより 2.15 × 10¹² に丸められます。浮動小数点対応のハードウェアは通常、これらの桁をすべて占有するのではなく、固定された桁数を対象として演算を行います。確保された桁数が p の場合、これより小さいオペランドが右にシフトされると、各桁は単に破棄されます (丸められません)。2.15 × 10^12 - 1.25  × 10^-5 は次のようになります。

x = 2.15 × 10¹²
y = 0.00 × 10¹²
x - y = 2.15 × 10¹²

解は、差を正確に計算した場合とまったく同じになり、この結果が丸められます。別の例を考えてみます。10.1 - 9.93 は次のようになります。

x = 1.01 × 10¹
y = 0.99 × 10¹
x - y = 0.02 × 10¹

正確な解は 0.17 なので、計算上の差は 30 ulp となり、全部の桁に誤差が生じています。なぜ、これほど大きな誤差が出るのでしょう。

定理 1

パラメータ b と p を使用した浮動小数点フォーマットで p 桁の差を計算した場合、結果の相対誤差は最大、b - 1 になる。

証明

式 x - y で β - 1 の相対誤差は、x = 1.00...0 、y = .ρρ...ρ (ρ = β - 1) の場合に生じます。ここで y の桁数は p です (すべて ρ)。正確な差は x - y = β^-p です。ただし p 桁だけで解答を求めると、y の最下位桁がシフトオフするので、計算上の差は β^-p+1 です。したがって、誤差は β^-p - β^-p+1 = β^-p (β - 1) で、相対誤差は β^-p(β - 1)/β^-p = β - 1 です。■

β = 2 の場合、相対誤差は結果と同じくらい大きくなることがあります。またβ = 10 の場合、9 倍になることがあります。言い換えると、β = 2 の場合、上記の式 (3) は汚染桁の数が log₂(1/ε) = log₂(2^p) = p になることを示します。すなわち結果の p 桁はすべて誤っていることになります。この状況を解決するために、保護桁として 1 桁分を追加してみます。すなわち小さい数字を p + 1 桁に切り捨てて、減算の結果を p 桁に丸めます。この保護桁を使用すると、前の例は次のようになります。

x = 1.010 × 10¹
y = 0.993 × 10¹
x - y = 0.017 × 10¹

また解も正しくなります。110 - 8.59 の場合、1 桁の保護桁を設けると、結果の相対誤差は、εより大きくなることがあります。

x = 1.10 × 10²
y = 0.085 × 10²
x - y = 1.015 × 10²

これを丸めると 102 になります。一方、正しい解は 101.41 で、相対誤差は 0.006 となり、ε = 0.005 より大きくなります。一般には、結果の相対誤差は εよりわずかに大きくなります。正確に表わすと次のようになります。

定理 2

x と y が β と p のパラメータを使用したフォーマットになっている場合に、 p + 1 桁 (1 保護桁) で減算を行うと、結果の相対誤差は 2ε より少なくなる。

この定理については 「丸め誤差」で証明します。なお、x と y は正負のいずれであってもかまわないので、上記の定理では加算を使用しています。

相殺

この最後の項の結論は、保護桁がなければ、2 つの近い値で減算をした場合に生じる相対誤差は、きわめて大きくなる可能性がある、ということです。すなわち減算 (負符号を使用した加算 ) を伴う式を評価した結果、相対誤差が大きくなりすぎて、桁がすべて無意味になる可能性があります (定理1)。互いに近い値を減算すると、オペランドの最上位桁がそれぞれ一致して、相互に相殺されます。この相殺には悪性と良性の 2 種類があります。

悪性の相殺は、オペランドが丸め誤差の影響を受ける場合に生じます。例えば、二次方程式の根の公式の中にある b² - 4ac の場合です。b² と 4ac は浮動小数点の乗算の結果となるので、いずれも丸め誤差の影響を受けます。仮にこれらの値が、最も近い浮動小数点数に丸められ、精度が 0.5 ulp 以内になっているとします。減算を行うときの相殺によって正確な桁が多数消失し、丸め誤差によって汚染された桁だけが残されることがあります。したがって、差の中に ulp の大きい誤差が含まれる可能性があります。例えば、b = 3.34、a = 1.22、c = 2.28 であると仮定します。 b² - 4ac の正確な値は、0.0292 です。しかし、 b² は 11.2 に、また 4ac は 11.1 にそれぞれ丸められるので、最終的な解は 0.1 となり、仮に 11.2 - 11.1 が 0·1 であったとしても、70 ulpもの誤差が出ます ⁶。減算自体が誤差の原因になったのはなく、以前の乗算で導入された誤差が明らかになったということです。

一方、良性の相殺は、既知の数量を減算するときに起きます。x と y に丸め誤差がない場合、保護桁を使用して減算を行うと、定理 2 により x - y の差の相対誤差は非常に小さく (2ε 以下) なります。

悪性の相殺の原因となる公式は編成を変更することにより、問題を解決できる場合があります。次の二次方程式の根の公式の例を見てみます。

(4)

b² が の場合、 には相殺が行われないので、 となります。ただし、2 つの公式のうちのもう一方の加算 (減算) が原因となって悪性の相殺が行われます。この問題を避けるには、次のように r₁ の分子と分母を乗算して (r₂ の場合も同様に)、

次の式を求めます。

(5)

で、b > 0 の場合、公式 (4) を使用して r₁ を計算すると、相殺が行われます。したがって、r₁ の計算には公式 (5)、また r₂ の計算には公式 (4) を使用するようにしてください。一方、b < 0 の場合は、r₁ の計算に公式 (4)、また r₂ の計算に公式 (5) を使用してください。

x² - y² の式も悪性の相殺を伴う公式です。(x - y)(x + y) として評価する方が正確です⁷。二次方程式の根の公式とは異なり、この式の方にも減算を伴いますが、丸め誤差が介在しない良性の相殺なので悪性ではありません。定理 2 により、x - y の相対誤差は多くても 2ε に抑えられます。これは、x + y にもあてはまります。相対誤差の少ない数量を乗算しても、積の相対誤差は小さくなります (「丸め誤差」を参照) 。

正確な値と計算値の混乱を避けるために、次のような表記を使用します。まず x - y は、x と y の正確な差を表わすものとします。また x から y を減算した計算値を表わす場合は、x  y と表記します。同様に、 、 、 は順に加算、乗算、除算の各計算値を表わします。すべて大文字の場合、特定の関数の計算値 (例えば LN(x)、SQRT(x) など) を表わします。小文字と一般的な数学記号は、値そのもの (例えば ln(x)、 など) を表わすものとします。

x² - y² の近似値を表わす場合、(x y) (x y) という表記はきわめて効果的ですが、x と y の浮動小数点数自体が と の真の値の近似値になることがあります。例えば、 と が、2 進表記では正確には表わせない 10 進数である場合があります。この場合、x y が x - y の近似値を表わすにしても、真の式 - に比べ、相対誤差は大きくなるので、x² - y² に対する (x + y)(x - y) の利点は、それほど大きくはなくなります。(x + y)(x - y) と x² - y² の計算負荷はほとんど同じなので、この場合に限り (x + y)(x - y) の方が有利ということになります。ただし一般には、入力が近似値ではないことが多いので、悪性の相殺を良性の相殺に置き換える意味はありません。とは言え、データが正しくない場合でも、相殺を全面的に排除 (二次方程式の根の公式のように) できれば、意味があります。この付録では、あるアルゴリズムに対する浮動小数点の入力が正しく、結果も可能な限り正確に計算されることを前提とします。

x² - y² の式は、(x - y)(x + y) として書き換えると、悪性の相殺が良性の相殺に置き換えられるので、より正確になります。

次に、悪性の相殺を良性の相殺に置き換えることのできる公式の例を紹介します。

三角形の面積は、次のように 3 辺 a、b、c の長さで直接表わすことができます。

(6)

三角形が非常に平たい形状である (a b + c) と仮定します。s a となるので、式 (6) の項 (s - a) では、2 つの近い値 (一方に丸め誤差が含まれる可能性があります) の減算を行います。例えば a = 9.0、b = c = 4.53 の場合、正しい値は 9.03 で A は 2.342...となります。計算値 s (9.05) に 2 ulp の誤差がありますが、A の計算値は 3.04 なので誤差は 70 ulp となります。

平坦な三角形の場合でも、正確な結果が得られるように、公式 (6) を次のように書き換えることができます (Kahan、1986)。

(7)

a、b、c は a ≧ b ≧ c の条件を満たさない場合は、それぞれの名前を変更してから、公式 (7) を使用します。 (6) と (7) の右辺が代数的に同一であることは比較的簡単にチェックできます。a、b、c の値により、面積 2.35 が求められます。この誤差は 1 ulp で、最初の公式よりはるかに正確であることがわかります。

この例の場合、公式 (7) は公式 (6) よりはるかに正確であるので、一般に公式 (7) が有効であると言えます。

定理 3

公式 (7) を使用して三角形の面積を求めたときに起きる丸め誤差は、保護桁で減算を行い、 e ≦ 0.005 で、平方根が 1/2 ulp 以下の範囲で計算される限り、11ε以内に抑えられる。

e < 0.005という条件は、事実上あらゆる浮動小数点システムで満たされます。例えば、β = 2 の場合、 p ≧ 8 であれば、必ず e < 0.005 になり、β = 10 であれば p ≧ 3 で十分ということになります。

式の相対誤差に関する定理 3 の記述では、式が浮動小数点演算によって計算されることがわかります。実際に、相対誤差は式に関するものです。

SQRT((a (b c)) (c (a b)) (c (a b)) (a (b c))) 4 (8)

公式 (8) の構造は複雑に入り組んでいるので、定理の記述では、E を丸付きの表記で表わすのではなく、E の計算値として扱います。

誤差の範囲は悲観的すぎる場合があります。上記のような数値の例の場合、公式 (7) の計算値は 2.35 で、正しい値である 2.34216 と比較すると、相対誤差は 0.7εとなります。これは、11εよりはるかに小さい値です。誤差の範囲を計算する主な目的は、正確な境界を設定することではなく、公式に数値の問題が存在しないことを確認することにあります。

良性の相殺に変更できる最後の例は (1 + x)ⁿ (ここで ) です。この式は金融計算でよく使用されます。例えば、年率 6 % の銀行口座に毎日 100 ドルずつ預金した場合の複利計算の例を考えてみます。n = 365 で、i = 0.06 とすると、1 年後の貯蓄額は

100 ドルとなります。

β = 2、 p = 24 として計算すると結果は $37615.45 となり、実際の金額の $37614.05 とは $1.40 の差が出ます。この問題の理由は簡単です。式 1 + i/n では、0.0001643836 に 1 を加算するので、i/n の下位ビットが失われるということです。この丸め誤差は、
1 + i/n に n 乗すると、増幅します。

(1 + i/n)ⁿ という複雑な式は e^{nln(1 + i/n)} と書き換えることができます。この場合の問題は、x が小さい場合の ln (1 + x) の計算です。1 つの方法として ln (1 + x) x の近似値を使用することができます。この場合、支払いは $37617.26 となり、誤差は $3.21 なので、単純な公式よりも精度が落ちます。ただし、定理 4 に示すとおり、(1 + x) をさらに正確に計算する方法があります (HP、1982年)。この公式により、$37614.07 という結果が得られ、誤差がわずか 2 セント以内になります。

定理 4 では、LN(x) により ln (x) の近似値が 1/2 ulp 以内の誤差で求められることを前提とします。これにより、x の値が小さいときに、1 x が、x の下位ビットにある情報を喪失するために、LN (1 x) により ln (1 + x) の近似値が得られない問題が解決されます。すなわち、 のときに、ln (1 + x) の計算値は実際の値の近似値にはならないという問題です。

定理 4

ln (1+ x) を次の公式で計算した場合、

保護桁を使用して減算を行い、e < 0.1、または ln を 1/2 ulp 以内の誤差で計算すると、0 ≦ x < の場合に相対誤差は 5e 以内となる。

この公式は、x がどのような値であっても適用されますが、 の場合、すなわち本来の公式 ln (1 + x) で悪性の相殺が起きる状況に限り意味があります。公式は一見、意味不明のように思われますが、これが正しく動作する理由は簡単です。ln (1 + x) を

と書き換えてみます。左辺の要素はそのまま計算できます

が、右辺の要素 μ(x) = ln(1 + x)/x は、x に 1 を加算するときに大きな丸め誤差が生じます。しかし、ln(1 + x) x のため μ はほとんど定数です。したがって、x をわずかに変更するだけでも、大きな誤差を伴います。

すなわち、 の場合、 ( ) を計算すると、xμ(x) = ln(1 + x) のほぼ近似値が得られます。 と + 1 を正しく計算できる の値があるでしょうか。1 + はちょうど 1 x になるので、 = (1 x) 1という値です。

この項の結論として、既知の近い値を減算するとき (良性の相殺を伴う) は、保護桁により精度が保証されるということが言えます。不正確な結果をもたらす公式は、良性の相殺を使用して、数値の精度を上げることにより、書き換えることができます。ただし、この方法は保護桁を使用して減算を行う場合に限り利用できます。また保護桁を使用しても、加算器のビットが 1 つ多くなるだけなので、それほどの犠牲も伴いません。例えば 54 ビットの倍精度加算器の場合、余分なコストは 2 % 程度で済みます。この程度のコストによって、三角形の面積を求める公式 (6) や式 ln (1 + x) などの多数のアルゴリズムを実行する利点が得られます。最近のコンピュータには通常、保護桁が用意されているので、これに対応できないシステムは、ほとんどありません (Cray システムなどを除きます)。

正確な丸め操作

保護桁を使用して浮動小数点演算を行うと、正確に計算した後で最も近い浮動小数点数に丸めた場合よりも、精度が落ちます。後者の計算方法を「正確な丸め操作」と言います ⁸。定理 2 の前に紹介した例は、1 桁分の保護桁を確保したからといって、必ずしも正確な結果が得られるとは限らないことを示しています。前の項には、正しい動作を保証するために、保護桁を使用するアルゴリズムの例をいくつか紹介しました。ここでは、正確な丸めを行うためのアルゴリズムの例を示します。

ここまでの段階では、まだ丸めの定義を示していません。丸めの操作は、切り下げと切り上げの区切りをどこで付けるかという点 (例えば、12.5 を 12 とするか 13 とするかという問題) を除き、難しくありません。1 つの方法として、10 までの値を四捨五入して、0、1、2、3、4 を切り捨て、5、6、7、8、9 を切り上げるという考え方があります。したがって、12.5は 13 として扱います。DEC の VAX ではこのような四捨五入が採用されています。また、5 で終わる数は 2 つの丸め方法のちょうど境目にあるという意味で、丸めを行う回数全体の半分を切り上げ、半分を切り下げるという考え方もあります。この 50% ずつ 2 分するという方法には、丸めの結果が少なくても、最下位の桁が偶数になっていなければならないという条件を伴います。したがって、12.5 は、2 が偶数であるため、13 ではなく、12 に切り捨てられます。切り上げと偶数への切り捨てのどちらがよいかという問題について、Reiser と Knuth は偶数への切り捨ての方が望ましいと指摘しています (1975)。

定理 5

x と y を浮動小数点数として、x₀ = x、 x₁ = (x₀ y) y、...、x_n =(x_n-1 y) y と定義する。 と を偶数への切り捨てにより同じように丸めると、すべての n について x_n = x、また n ≧ 1 のすべてについて x_n = x₁ となる。■

この結果を分類する上で、x = 1.00、y = -0.555 として β = 10、 p =3 の場合を考えてみます。これを切り上げる場合、x₀ y = 1.56、x₁ = 1.56 0.555 = 1.01、
x₁ y = 1.01 0.555 = 1.57というシーケンスになり、x_n = 9.45 (n ≦ 845) になる⁹ まで、x_n の連続値が 0.01 ごとに増えていきます。偶数に切り捨てる場合、x_n はつねに 1.00 になります。この例は、切り上げのルールを使用すると、計算結果がしだいに上方にずれていき、一方、偶数への切り捨てを使用すると、このずれが生じないことを示しています。この後の説明では、偶数への切り捨てを一貫して使用することにします。

多重精度の演算では、正確な丸めが 1 回だけ行われます。精度を上げるには、2 つの基本的なアプローチがあります。1 つは、きわめて大きい有意桁を使用した浮動小数点数です。この有意桁はワード配列に格納され、アセンブリ言語でこれらの数を操作するためのルーチンのコーディングが行われます。もう 1 つのアプローチは、通常の浮動小数点数の配列としてさらに高い精度の浮動小数点数を表わす方法です。この場合、無限の精度で配列要素を追加することにより、高い精度の浮動小数点数を復元します。ここでは、2 つ目のアプローチについて説明します。浮動小数点数の配列を使用する利点は、高水準言語により移植可能な形でコーディングできることにあります。この反面、正確な丸めの演算が必要になります。

このシステムにおける乗算の重要な点は、各加数が x y と同じ精度を持つように、x y の積を和として表わすことにあります。これは、x と y を分割することによって行われます。x = x_h + x_l 、y = y_h + y_l と書くと、実際の積は xy = x_hyh + x_h y_l + x_l y_h + x_l y_l となります。 x と y に p ビットの有意桁がある場合 [p/2] ビットで x_l、x_h、y_h、y_l を表わすことができれば、加数の有意桁も p ビットになります。 p が偶数であれば、簡単に分割できます。x₀.x_l...x_{p -1}の数は、x₀.x_l...x_{p/2 - 1} と 0.0 ... 0x_p/2 ... x_{p - 1} の和として表わすことができます。一方、p が奇数の場合、このように簡単には分割できません。ただし、負数を使用して、1 ビット分を追加できます。例えば、β = 2、 p = 5、x = 0.10111 の場合、x は x_h= 0.11、および x_l = -0.00001として分割できます。Dekker (1971) によると、分割方法は簡単ですが、1 桁の保護桁以外にも余分に桁が必要です。

定理 6

β > 2 であっても p が偶数であるという制限のもとに、p を浮動小数点精度として、浮動小数点数を正確に丸めるものとする。このとき、k = [p/2] が精度の半分 (切り上げ) で、m = β^k + 1 とすると、x は x = x_h + x_l (ただし
x_h = (m x) (m x x)、x_l = x x_h で、各 x_i は [p /2] ビット精度で表現可能とする) として分割できる。

この定理が実際の例にどのように適用されるかを見るために、β = 10、 p = 4、b = 3.476、a = 3.463、c = 3.479 であると仮定します。b² - ac を最も近い浮動小数点数に丸めると、0·03480になり、b b = 12.08、a c = 12.05 になります。したがって、b² - ac の計算値は 0·03となります。これには 480 ulp の誤差があります。定理 6 にもとづいて b = 3.5 - 0.024、a = 3.5 - 0.037、c = 3.5 - 0.021 とすると、
b² は 3.5² - 2 × 3.5 × 0.024 + 0.024² となります。各加数はすべて正確なので、
b² = 12.25 - 0.168 + 0.000576 になります (この時点で和の評価はされていません)。同じように、ac = 3.52 - (3.5 × 0.037 + 3.5 × 0.021) + 0.037 × 0.021 = 12·25 - 0.2030 + 0.000777 になります。

最後に、2 つの連続項ごとに減算すると、b² - ac の近似値が
0 0.0350 0.000201= 0.03480 として求められます。これは、正確に丸めた場合の結果と一致します。定理 6 で正確な丸めが必要になることを確認するために、p = 3、β = 2、x = 7、さらに m = 5、mx = 35、m x = 32 の例を考えてみます。1 桁の保護桁を使用して減算を行うと、(m x ) x = 28 となります。したがって、x_h = 4、x_l = 3 となるので、x_l は [p /2] = 1 ビットでは表わせません。

正確な丸めの最後の例として、m を 10 で除算する場合を考えてみます。この結果、一般には m/10 に一致しない浮動小数点数が得られます。β = 2 の場合、正確な丸めを使用して、m/10 と 10 を掛けると、不思議なことに m がリストアされます。実際には、さらに一般的な事実があてはまります (Kahan)。証明は詳細にわたるものですが、このような詳しい説明まで関心のない方は、「IEEE 標準」まで読み飛ばしてもかまいません。

定理 7

β = 2 のとき、m と n が |m|<2^{p - 1} の整数で、n に n = 2ⁱ + 2^j という特殊な形式がある場合、浮動小数点演算を正確に丸めるとすれば、 (m n ) n ＝ m となる。

証明

2 のべき乗でスケーリングするのは、有意桁ではなく、単に指数が変更されるだけなので、悪影響はありません。q = m/n の場合、2^{p - 1} ≦ n < 2^p となるように n をスケーリングし、

< q < 1 となるように m をスケーリングします。したがって 2^{p - 2} < m < 2^p となります。m には p 個の有意ビットがあるので、2 進小数点の右に最大 1 ビットが必要です。m の符号を変更しても悪影響はないので、q > 0とします。

= m n の場合、定理を証明するには、次を示す必要があります。

(9)

これは、m の 2 進小数点の右側に最大 1 ビットが使用されているので、n により m に丸められます。

中間の場合の扱いについては、｜n - m｜= の場合に、スケーリングしていない初期の m は｜m｜< 2^{p - 1}、また下位ビットは 0 になるので、スケーリング後の m の下位ビットも 0 になります。したがって、中間の値は m に丸められます。
q = .q₁q₂ ...、 = .q₁q₂ ... q _p 1 であると仮定します。|n - m|の近似値を求めるには
まず、｜ - q｜=｜N/2^{p + 1} - m/n｜ (Nは奇数の整数) を計算します。n = 2ⁱ + 2^j、および 2^{p - 1} ≦ n < 2^p であるから、特定の k ≦ p - 2 については n = 2^{p - 1} + 2^k でなければなりません。したがって、次のようになります。

分子は整数です。N は奇数なので、実際には奇数の整数となります。したがって、
｜ - q｜≧ 1/ (n2 p+1-k ) となります。q < (q > の場合も同様) と仮定¹⁰すると、
n < mで、次のようになります。

｜m-n ｜= m-n = n(q- ) = n(q-( -2-p-1) ) ≦

= (2^p-1+2k)2^-p-1+2^-p-1+k =

これにより、 (9) が成り立ち、定理が証明されることになります¹¹。■

定理は、2ⁱ + 2^j をβⁱ + β^j に置き換える限り、任意の基数 β について真になります。ただし β が大きくなるに従ってβⁱ + β^j の分母の差が大きくなっていきます。

ここで、本付録の最初に提示した「基本的な演算処理で多少の丸め誤差が出ることに何か問題があるか？」という根本的な疑問に立ち返ることにします。この答えは「大いに問題はある」です。その理由は、正確な基本演算では相対誤差が低く抑えられるという意味で、公式の「妥当性」を証明できるからです。

各公式にわずかな相対誤差を伴うという意味で、公式が正しいことを証明できるという理由によって、丸め誤差の問題はあるということになります。「相殺」では、保護桁を使用するアルゴリズムによっては、上記の意味で正しい結果が得られるものがあると説明しました。ただし、これらの公式への入力が、不正確な測定にもとづく数値である場合、定理 3 と 4 で言う境界は無意味になってきます。この理由は、x と y が、ある測定値の近似値である場合、良性の相殺 x - y が悪性に変わるからです。しかし、正確な計算は、定理 6 と 7 で示した正確な関係を実証できるという意味で、データが不正確な場合でも有効であることに違いはありません。各浮動小数点変数がそれぞれ、実際の値の近似値にすぎない場合でも有効であるということです。

IEEE 標準

浮動小数点演算には、IEEE 754 と IEEE 854 という 2 つの標準があります。IEEE 754 は、単精度の場合 β = 2、 p = 24、また倍精度の場合 p = 53 を使用する 2 進数標準です (IEEE 1987)。この標準により、単精度と倍精度の場合の正確なビット・レイアウトが指定されます。IEEE 854 では、β = 2、または β = 10 のいずれかを使用します。IEEE 754 と異なり、IEEE 854 の場合は浮動小数点数をビットに符号化する形式は指定されません (Codyほか 1984)。p に対する特定の値は必要ありませんが、単精度と倍精度に認められる p の値が制限されます。なお IEEE 標準という用語は、以上 2 つの標準に共通する特性を表わす場合に使用します。

本項では、IEEE 標準の概要について説明します。以降の各項では、この標準の各特性と、これが採用された理由を示します。なお IEEE 標準が最も優れた浮動小数点標準であるかどうかを分析するのは本資料の主旨ではないので、あくまで標準としてそのまま受け入れ、その基本概念を説明することにします。詳細については、IEEE 標準自体の資料を参照してください (IEEE 1987、Cody ほか、1984)。

フォーマットと操作

基数

IEEE 854 で β = 10 が使用されている理由は明らかです。基数 10 は、日常的に親しみやすい数であるからです。β = 10 は、特に計算結果を 10 進数で表わすような電卓などに適しています。

IEEE 854 で、基数が 10 以外の場合に必ず 2 を使用する理由はいくつかあります。「相対誤差と ulp」で、理由の 1 つを説明しました。β を 2 にすると、相対誤差として計算した場合、0.5 ulp の丸め誤差が β の倍数で変動するので、誤差の分析結果がはるかに厳しくなり、相対誤差にもとづく誤差はほとんどの場合、単純に分析できるという理由です。また、基数を大きくした場合の有効な精度に関する理由も考えられます。β = 16、 p = 1 と、β = 2、 p = 4 の場合を比較してみます。いずれのシステムも有意桁は 4 ビットとします。例えば、15/8 を計算します。β = 2 の場合、15 は 1.111 × 23 と、また 15/8 は 1.111 × 20 と表現されます。したがって、15/8 はそのまま正確に表わされます。一方、β = 16 の場合、15 は F × 160 (F は 15 に対応する 16 進表現) と表わします。しかし、15/8 は 1 × 160となり、正しいビットは 1 ビットしかありません。一般に、基数 16 を使用すると最大 3 ビットを失うことがあるので、 p の 16 進桁の精度は、有効精度が 2 進の 4p 分ではなく、4p -3 にまで低下することがあります。β の値が大きくなるに従って、こうした問題が出てくるとすれば、IBM は System/370 に、なぜ β = 16 を採用したのでしょう。その真意は IBM 以外に知るよしもありませんが、2 つの理由が考えられます。まず、指数の範囲が広がるという点です。System/370 の単精度では、β = 16 と p = 6 が使用されています。したがって、有意桁には 24 ビットが必要になります。24 ビットは 32 ビットに収まるので、指数用の 7 ビットと符号用の 1 ビットが残ります。したがって、表現可能な数字の範囲は 16-2 6から 162 6 = 22 8 になります。β = 2 で同じ指数範囲を確保するには、指数分に 9 ビットが必要となり、有意桁にはわずか 22 ビットしか使用できません。一方、β = 16 とすると、前述のとおり、有効精度が 4p - 3 = 21ビットにまで低下するという問題があります。これも、β = 2 で 1 ビット分、精度を上げる (本項で後述) ことができる場合は、さらに悪くなります。したがって、β = 2 のマシンには、β = 16 のマシンの 21 ～ 24 ビットに相当する 23 ビット分の精度が確保されます。

IBM が β = 16 を採用したもう 1 つの理由はシフトに関するものです。2 つの浮動小数点数を加算した場合に、各指数が異なっていると、基数の位置を合わせるために、いずれかの有意桁をシフトしなければならないことがあります。これが原因でパフォーマンスが低下することがあります。β = 16、p = 1 のシステムでは、1～15 までの数はすべて指数が同じになるので、( ) = 105 のどの組み合わせを加算してもシフトは必要ありません。一方、β= 2、p = 4 のシステムでは、これらの数の指数の範囲は 0～3 なので、105 ペアのうちの 70 についてはシフトが必要になります。

最近のハードウェアでは、オペランドのサブセットに対するシフトを避けても、パフォーマンス上、それほど大きな利点は得られないので、ワブルの少ない β= 2 の方が有利な基数であると言えます。β= 2 のもう 1 つの利点は、有意桁の余分なビットが得られる点です¹²。浮動小数点数は必ず正規化されるので、有意桁の最上位ビットは必ず 1 になり、これに対応するビットを無駄にするのは意味がありません。この手法を効果的に利用したフォーマットのことを「隠れビット」と言います。「浮動小数点フォーマット」では、0 に対する特殊な規約が必要であることを説明しました。そこで示した方法は、e_min - 1 の指数と、すべてゼロの有意桁により、
ではなく、0 を表わすというものです。

IEEE 754 の単精度は、符号用の 1 ビット、指数用の 8 ビット、および有意桁用の 23 ビットを使用して、32 ビットをエンコードします。ただし、隠れビットを使用するので、23 ビットだけを使用してエンコードした場合でも、有意桁は 24 ビット (p =24) になります。

精度

IEEE 標準では、単精度、倍精度、拡張単精度、および拡張倍精度という 4 種類のフォーマットを定義しています。IEEE754 の場合、単精度と倍精度はほぼ、大部分の浮動小数点ハードウェアで要求される機能に対応することができます。単精度では 32 ビットのシングルワード、また倍精度では 2 つ連続した 32 ビットワードが占有されます。拡張精度は、精度を多少上げ、指数の範囲を拡張したフォーマットです
(表 D-1) 。

表 D-1   IEEE 754フォーマットのパラメータ

パラメータ	フォーマット
	単精度	拡張単精度	倍精度	拡張倍精度
p	24	≧32	53	≧64
emax	+127	≧1023	+1023	> 16383
emin	-126	≦-1022	-1022	≦ -16382
指数の幅 (ビット)	8	≦11	11	≧15
フォーマットの幅 (ビット)	32	≧43	64	≧79

IEEE 標準は、拡張精度によって追加する余分なビット数の下限だけを指定します。拡張倍精度フォーマットの最低許容限度は、上の表によると 79 ビットになっていますが、便宜上 80 ビットフォーマットと呼ばれることがあります。この理由は、拡張精度のハードウェア実装では一般に隠れビットは使用しないので、79 ビットではなく、80 ビットを使用するからです。¹³

IEEE 標準では、拡張精度に最も重点が置かれるので、倍精度に関する推奨は何もありません。代わりに、サポートされる最も広い基本フォーマットに対応する拡張フォーマットをインプリメントするように強く推奨されています。

拡張精度が採用された理由は、電卓によるところがあります。電卓は通常、10 桁を使用しますが、内部的には 13 桁を使用します。13 桁のうち、10 桁だけを表示することにより、電卓は一見、指数、コサインなどの関数を 10 桁までの精度で計算する「ブラックボックス」を使用しているように思われます。指数、対数、コサインなどの関数を 10 桁までの精度で比較的効率よく計算できる電卓では、余分な桁がいくつか必要になります。結局、500 ユニット程度の誤差で対数の近似値を計算する有理式を見つけるのはそれほど難しくはありません。したがって、13 桁を使用して計算することにより、10 桁で正しい解が得られるようになります。この余分な 3 ビットを隠れビットとして使用することにより、オペレータは単純なモデルによって電卓を操作することができます。

IEEE 標準の拡張精度も同じような作用を持っています。これにより、各ライブラリは単精度 (または倍精度) により 0.5 ulp 程度の誤差の範囲で効率的に計算できるようになります。したがって、ユーザーは単純なモデルによって各ライブラリを使用することができます。例えば、乗算や対数の呼び出しなどに限らず、基本的な操作によっても、0.5 ulp 以内の誤差で正確に値が返されるようになります。ただし拡張精度を使用するときは、これを使用していることをユーザーにも明らかにすることが重要です。例えば、電卓で表示する値の内部表現が、表示形式と同じ精度で丸められている場合、以降の操作の結果は隠れビットに依存するので、ユーザーには予測できなくなることがあります。

拡張精度についてさらに詳しく説明するために、IEEE 754 の単精度表現と 10 進数を変換するときの問題について考えてみます。可能であれば、10 進数を読み戻す時点で、単精度数も復元できるように、十分な桁数を確保して単精度数を出力するのが理想的です。単精度の 2 進数を復元するのに、9 桁の 10 進数があれば十分であることがわかっています (「2 進数から 10 進数への変換」を参照してください)。10 進数を一意の 2 進表現に戻す場合、1 ulp 程度の小さい丸め誤差があっても誤った解になるので、致命的となります。拡張精度が効率的なアルゴリズムに不可欠となるような状況を次に説明します。単精度の場合、10 進数を単精度の 2 進数に変換するには、きわめて簡単な方法があります。まず、整数 N として 9 桁の 10 進数を読み取り、小数点は無視します。

表 D-1 から、p ≧ 32 の場合、10⁹ < 2³² 4.3 × 10⁹ となるので、N は単精度拡張フォーマットで正確に表現できることがわかります。次に N をスケーリングするのに必要な 10P のべき乗を求めます。これは、10 進数の指数と、その時点までに無視された小数点の位置を考慮した組み合わせになります。10^|P| を計算します。|P| ≦ 13 であれば、10¹³ = 2¹³ 5¹³ で、5¹³ < 2³² となるので、これも正確に表現することができます。最後に、N と 10|P| を乗算 (または p < 0 の場合は除算) します。この最後の計算が正しく行われれば、最も近い 2 進数が求められます。最後の乗算 (除算) を正確に行う方法については、「2 進数から 10 進数への変換」で詳しく説明しています。したがって、|P| ≦ 13 の場合、単精度拡張フォーマットを使用することにより、9 桁の 10 進数が最も近い 2 進数に変換 (すなわち正確な丸め) されることになります。また |P| > 13 の場合、単精度拡張フォーマットでは上のアルゴリズムに十分対応できないので、必ずしも正確に 2 進数が計算できるとは限りません。ただし Coonen (1984) により、2 進数から 10 進数への変換と、その逆の操作で、もとの 2 進数が復元されることが実証されています。

倍精度がサポートされていれば、上記のアルゴリズムは拡張単精度フォーマットではなく倍精度で実行されます。ただし、倍精度フォーマットを 17 桁の 10 進数に変換し、再度戻す場合には拡張倍精度フォーマットが必要になります。

指数

指数には負と正の 2 つの関係が考えられるので、符号を表現するための特定の方法を選択しなければなりません。符号付きの数は、符号と数量を組み合わせる方法と、2 の補数で表現する方法があります。符号/数量の形式は、IEEE フォーマットによる有意桁の符号を表わす場合に使用されます。すなわち 1 ビットを符号用として確保し、残りのビットで数量を表わします。また 2 の補数で表わす方法は整数演算のときによく使用されます。この方法では、[ -2^p-1 , 2^p-1 - 1] の範囲にある数を 2^p の剰余に相当する非負の最小値で表わすことができます。

IEEE 2 進数標準で指数を表現する場合は、いずれの方法も使用しません。代わりに「バイアス表現」を使用します。指数を 8 ビットで格納する単精度の場合、バイアスは 127 となります (倍精度の場合は 1023 です)。 が符号なしの整数を表わす指数ビットである場合、浮動小数点数の指数は - 127 になります。これはバイアス指数 と区別するために、「バイアスなしの指数」と呼ばれることがあります。

表 D-1 を見ると、単精度の場合は e_max = 127、e_min = -126 であることがわかります。|e_min |< e_max となるのは、最小値 の逆数がオーバーフローしないという理由によるものです。最大値の逆数がアンダーフローすることはありますが、これはオーバーフローの場合ほど深刻な問題ではありません。「基数」では、0 を表わすのに e_min - 1を使用すると説明しました。また 「特殊な数」では、e_max + 1 について説明します。IEEE 単精度の場合、これはバイアス指数の範囲が
e_min - 1 = 127 から e_max + 1 = 128 にあるということになります。一方、バイアスなしの指数の範囲は、0 から 255 までとなり、8 ビットで表現できる非負の数に相当します。

操作

IEEE 標準では、四則演算の結果が正確に丸められます。すなわち結果を正確に計算した後、最も近い浮動小数点数に丸める (偶数への丸め) ということです。「保護桁」では、正確な差や 2 つの浮動小数点数の和を計算するのは、指数が大きく異なるときにコスト高になると指摘しました。さらに、相対誤差を最小限に抑制しながら差を計算する便利な方法も紹介しました。ただし、1 桁の保護桁を使用して計算したからといって、必ずしも、正確な結果を計算した上で丸めを求めた場合と同等の結果が得られるとは限りません。

そこで 2 桁目の保護桁と 3 桁目のスティッキビットを使用することにより、1 桁の保護桁の場合より多少のコストは伴うものの、正確な計算の後に丸めを行なった場合と同じ結果を得ることができます (Goldberg 1990)。この方法により、標準をより効率的に実現できます。

算術演算の結果を完全に指定することにより、ソフトウェアの移植性が改善されます。あるマシンから別のマシンにプログラムを移植するときに、両方のマシンで IEEE 標準がサポートされていれば、その間に生成される中間結果はすべてソフトウェアバグに原因があり、演算上の差によるものでないことが確信できます。また正確な指定によって、浮動小数点に関する推論も行いやすくなります。浮動小数点に関する証明は、さまざまな演算処理が原因で起きる各種の問題に対処する必要がない場合でさえ、難しいものです。整数プログラムが正しいことを証明できるのと同じように、浮動小数点プログラムの場合も、証明する内容が計算結果の丸め誤差によって特定の条件が満たされるというものであっても、証明することができます。定理 4 は、このような証明の例を示したものです。こうした証明は、推論する操作の内容が正確に指定できれば、はるかに容易になります。いったん、あるアルゴリズムが IEEE 演算に適合していることがわかれば、IEEE 標準をサポートするマシンである限り、いずれであっても正しく動作します。

Brown (1981) は、一般的な浮動小数点ハードウェアに見られる原理を発表しました。ただし、このシステムでの証明が必ずしも、「相殺」や 「正確な丸め操作」に示すアルゴリズムの妥当性を裏付けることにはなりません。これらのアルゴリズムの特性がすべてのハードウェアでサポートされているとは限らないからです。さらに、Brown の原理は、単に正確な計算と丸め操作に関する単純な定義よりも、一段と複雑になっています。したがって、Brown の原理をもとに定理を導くのは、操作が単に正確な丸めであることを前提として証明する場合よりも難しくなります。

浮動小数点標準によって対応する操作の内容については、全面的な同意が確立されているわけではありません。+、-、×、/ の四則演算以外に、IEEE 標準では平方根、剰余、整数/浮動小数点の変換の各操作を正確に丸めるように規定されています。また、内部フォーマットと10 進数の変換も正確に丸めなければなりません (ただし、きわめて大きい数の場合は除きます)。Kulisch と Miranker (1986) は、正確に指定された操作のリストに内積を追加することを提案しました。IEEE 演算にもとづいて内積を計算すると、最終的な解が大きく異なることがあると指摘されています。

例えば、和は内積の特殊なケースなので、 ((2 × 10^-30 + 10³⁰) - 10³⁰ ) -10^-30 の和は
10^-30と同じになります。ただし IEEE 対応のマシンでは、計算結果が -10^-30 となります。高速の乗算器を実現するより少ないハードウェアコストで、誤差を 1 ulp 以内に抑えて内積を計算することができます (Kirchner/Kulish、1987)。^{14 15}

標準として規定される操作はすべて、10 進数と 2 進数の変換処理の場合を除き、正確に丸めを行う必要があります。これは、変換処理を除き、すべての操作について正確に丸められる効率的なアルゴリズムが実証されているからです。なお変換処理の場合は、最も効果的であるとされるアルゴリズムでも、正確な丸めによる計算に比べ、多少精度が落ちます (Coonen、1984) 。

IEEE 標準では、「テーブル・メーカのジレンマ」があるために、超越関数は必ずしも正確に丸める必要はありません。これを説明するために、4 つの小数点の位置で指数関数のテーブルを作成する場合を考えてみます。exp (1.626) = 5.0835 とします。これを丸めると、5.083、または 5.084 のいずれになるのでしょう。exp (1.626) を慎重に計算すると、5.08350 となります。さらに 5.083500、次に 5.0835000 になります。exp が超越関数である場合、exp (1.626) が 5.083500 ...0ddd、または 5.0834999 ...9ddd のいずれであるかを判断できるまでに相当時間がかかる可能性があります。したがって、超越関数について、無限大の精度で計算した後、その結果を丸めた場合と同じような精度を求めるのは現実的ではありません。別のアプローチとして、超越関数をアルゴリズム的に指定することもできます。ただし、すべてのハードウェアアーキテクチャに適用できる 1 つのアルゴリズムが存在するわけではありません。最近のマシンで超越関数を計算する場合に、有理近似値、CORDIC¹⁶、および大型テーブルという 3 つの手法を使用できます。いずれの手法もそれぞれ独自のクラスのハードウェアに適用されるもので、現在のところ、各種のハードウェアに幅広く対応できる単一のアルゴリズムはありません。

特殊な数

浮動小数点ハードウェアによっては、各ビットパターンが有効な浮動小数点数を表わす場合があります。IBM System/370 は、その一例です。一方、VAX^TM では、「予約オペランド」という特殊な数を表わす特定のビットパターンが確保されています。この考え方は CDC 6600 に由来するもので、INDEFINITE と INFINITY という特殊な数のビットパターンが使用されます。

IEEE 標準には、この伝統を継承して、NaN (Not a Number ) と無限大という概念があります。特殊な数がサポートされていないと、負数の平方根などの例外状況が起きた場合、計算をアボートする以外に対処する方法がありません。IBM System/370 FORTRAN では、-4 のような負数の平方根の計算に出会うと、デフォルトの動作としてエラーメッセージを出力します。各ビットパターンは、有効な数を表わすので、平方根の戻り値は浮動小数点数でなければなりません。System/370 FORTRAN の場合 が戻されます。一方、IEEE 演算の場合は、NaN が戻されます。

IEEE 標準により、次の特殊な値が指定されます (表 D-2 を参照してください)。±0、非正規化数、± 、および NaN (次の項で説明するとおり、NaN が存在します) です。これらの特殊な値は、e_max + 1 または e_min - 1 の指数によってすべてエンコードされます (なお、0 には e_min - 1 の指数があることはすでに指摘したとおりです)。

表 D-2   IEEE 754 の特殊な値  

指数	関数	表示
e = emin - 1	f = 0	±0
e = emin - 1	f ≠ 0
emin ≦ e ≦emax	--	1.f × 2e
e = emax + 1	f = 0
e = emax + 1	f ≠ 0	NaN

NaNs

これまで、0/0 や などの計算は、その時点で計算の停止の原因となる回復不能なエラーとして扱われてきました。ただし、こうした状況でも計算をそのまま続行することに意味があるケースも考えられます。たとえば、ある関数 f のゼロを求めるサブルーチン、zero(f)の例を考えてみます。従来、ゼロ検出サブルーチンが検索を行う期間として、関数が定義された間隔 [a,b] をユーザーが指定しなければなりませんでした。このサブルーチンのことを zero(f,a,b) と言います。さらに便利なゼロ検出サブルーチンでは、この追加情報を入力する必要はありません。この汎用的なゼロ検出サブルーチンは、単に関数の入力だけが要求され、領域の指定などは無意味となるような電卓などに適用できます。しかし、実際にはほとんどのゼロ検出サブルーチンで領域の指定が必要になります。その理由は簡単です。ゼロ検出サブルーチンはいくつかの値で関数 f を調べることにより処理を行います。f の領域を超えた値があると、f に対するコードが 0/0 または になり、計算が停止して、ゼロ検出プロセスをアボートしてしまうことがあるからです。

この問題は、NaN という特殊な値を導入し、0/0 や の式の計算では、処理が停止する代わりに NaN を生成するように指定すれば、解決します。表 D-3 には、NaN の原因となる状況をいくつか示してあります。zero(f) で f の領域を超えた値があると、f のコードにより NaN が返され、ゼロ検出サブルーチンの処理が続行します。すなわち、zero(f)は誤った推測をしたことについて「とがめられる」ことはありません。この例をもとに、NaN と通常の浮動小数点数を組み合わせた場合の結果を簡単に確認することができます。f の最後の文は return (-b + sqrt(d))/(2*a) となります。d < 0 の場合、f は NaN を返します。d < 0 なので、NaN と別の数の和が NaN の場合は、sqrt(d) は NaN、-b + sqrt(d) も NaN になります。同様に、除算のいずれかのオペランドが NaN の場合、商も NaN になります。一般に、浮動小数点計算で NaN が要求された場合は、結果も NaN になります。

表 D-3   NaN が生成される操作

操作	NaN が生成される原因
+	+ (- )
x	0 x
/	0/0, /
REM	x REM 0, REM y
	(x < 0 のとき)

ユーザーが領域を入力する必要のないゼロ・ソルバを作成する別のアプローチとして、シグナルを使用することもできます。ゼロ検出サブルーチンは、浮動小数点例外に対するシグナルハンドラをインストールすることができます。f が領域外と評価され、例外が発生すると、制御はゼロ・ソルバに渡されます。このアプローチの問題点は、各言語ごとにシグナルを処理する方法が異なる (その方法がない場合もあります) ので、移植性がまったくないという点です。

IEEE 754 では、NaN は指数 e_max + 1 とゼロ以外の有意桁を伴う浮動小数点数として表わすことがあります。実装ごとに、システム依存の情報を有意桁に入れることができます。したがって、一意の NaN が存在するわけではなく、一連の NaN で全体が構成されることになります。NaN と通常の浮動小数点数を組み合わせると、結果は NaN オペランドと同じになります。したがって、長い計算の結果が NaN の場合、有意桁にあるシステム依存の情報は、最初の NaN が生成された時点で生成される情報となります。実際に、最後の文には注意が必要です。両方のオペランドが NaN の場合、結果はいずれか一方の NaN になるので、必ずしも最初に生成された NaN であるとは限りません。

無限大

NaN を使用して、0/0 や のような式が検出されたときに計算を続行させる場合と同じように、無限大によってオーバーフローが起きたときに処理を続けることができます。これは単に最大許容数を戻すよりも、はるかに安全な方法です。
例えば、β= 10、 p = 3、e_max = 98 のときの を考えてみます。x = 3 × 10⁷⁰、y = 4 × 10⁷⁰ の場合、x² はオーバーフローし、9.99 × 10⁹⁸ に置き換えられます。同じように、y²、および x² + y² はそれぞれ順にオーバーフローして、9.99 × 10⁹⁸ に置き換えられます。最終的には、 となり、正解の 5 × 10⁷⁰ とはほど遠い結果となります。IEEE 演算の場合、x² の結果は、y² 、x² + y² 、 と同じように、 となります。最終的な結果は となり、これは正しい解から大きくはずれた通常の浮動小数点数を戻すより、はるかに安全であると言えます¹⁷。

0 を 0 で除算しても、NaN になります。ゼロ以外の数を 0 で割ると、 が返されます (1/0 = 、-1/0 = - )。これを区別する理由は次のとおりです。x がある限界に接近するに従い、f(x) → 0、g(x) → 0 になる場合、f(x)/g(x) の値は任意になります。例えば、f(x) = sin x、g(x) = x の場合、x → 0 だと f(x)/g(x) → 1 ですが、f(x) = 1- cos x の場合、f(x)/g(x) → 0 になります。0/0 を、2 つの小さい値の商を制限する条件として捉えると、0/0 は任意になります。したがって、IEEE 標準で、0/0 は NaN になります。ただし、c > 0 の場合、f(x) → c、g(x) → 0 であれば、あらゆる分析関数 f と g について f(x)/g(x) →± となります。小さい x に対して g(x) < 0 の場合、f(x)/g(x) →- となります。これ以外の場合は、限界が + になります。したがって、IEEE 標準では、c ≠ 0である限り、c/0 =± が定義されます。 の符号は、通常の c と 0 の符号に依存するので、-10/0 = - 、また -10/-0 = + となります。オーバーフローが原因で生成される と、ゼロによる除算が原因で起きる との違いは、ステータスフラグをチェックすれば区別できます (詳細については、「フラグ」で説明していま
す ) 。前者の場合はオーバーフロー・フラグが、また後者の場合はゼロによる除算のフラグがそれぞれセットされます。

オペランドとして無限大を生成する操作の結果を決定するルールは簡単です。
無限大を有限数 x に置き換え、その限界を x → に設定します。
したがって、 なので、3/ = 0 になります。

同様に、4 - = - 、 = になります。限界がない場合、結果は NaN になるので、 / も NaN になります (表 D-3に追加の例を示しています)。これは、0/0 が NaN になると結論する場合の推論と一致します。

ある部分式が NaN と評価されると、式全体の値も NaN になります。ただし、± の場合、1/ = 0 のようなルールがあるので、式の値は通常の浮動小数点数になる場合があります。ここでは、無限大の演算についてルールを適用する実際の例を示します。
関数 x/(x² + 1) を計算する場合を考えてみます。

これは、x が より大きくなるとオーバーフローするので、よい公式ではありませんが、無限大の計算により、1/x の近似値ではなく 0 が生成されるので、誤った解が得られます。ただし x/(x² + 1) は 1/(x + x^-1) と書き換えることができます。この書き換えにより、早い段階でオーバーフローが起きることはなくなり、無限大の演算により、x = 0 : 1/(0 + 0^-1) = 1/(0 + ) = 1/ = 0 のときに正しい値が得られます。無限大の演算がないと、1/(x + x^-1) の式では、x = 0 かどうかをテストする必要があります。この場合、余分な命令が追加されるほか、パイプラインが破壊される可能性があります。この例は、無限大の演算によって、特殊なケースをチェックする必要がなくなるという一般的な事例を示したものです。ただし公式については、無限大時に、疑似的な動作 (x/(x² + 1) のように) を示すことがないように、慎重に調べる必要があります。

符号付きゼロ

ゼロは、指数 e_min -1 とゼロの有意桁によって表わすことができます。符号ビットは 2 つの値を取ることができるので、ゼロにも +0と -0 の 2 種類があります。+0 と -0 を比較する場合に、区別してしまうと、if (x=0)のような簡単なテストでも、x の符号によって予測できない結果になることがあります。したがって、IEEE 標準では、-0 < +0 ではなく、-0 = +0 となるように定義されています。ゼロの符号をつねに無視することも可能ですが、IEEE 標準では無視されません。乗算、または除算で符号付きのゼロを伴う場合、解の符号を計算するときには通常の符号が適用されます。したがって、3·(+0) = +0、および +0/-3 = -0になります。またゼロに符号がない場合は、x =± の場合、1/(1/x) = x は成り立たなくなります。その理由は 1/- と 1/+ はいずれも 0、1/0 は + となり、符号情報が失われるからです。1/(1/x) = x の意味を保持するには、一方の無限大だけを使用します。ただしこの場合も、オーバーフローした数の符号が失われるという結果は変わりません。

符号付きゼロを使用するもう 1 つの例は、アンダーフローと、log のような 0 で不連続となる関数に関するものです。IEEE 演算では、x < 0 のときに、log 0 = - 、log x を NaN と定義するのが当然になっています。ゼロにアンダーフローした小さい負数を表わす x の例を考えてみます。符号付きのゼロのために、x は負数になるので、log は NaN を戻すことができます。ただし、符号なしのゼロの場合、log 関数はアンダーフローした負数と 0 とを区別できなくなるので、- を返します。ゼロで不連続になるもう 1 つの関数の例として、数値の符号を返す signum 関数があります。

符号付きゼロの最も興味のある例は、複素数演算です。

単純な例として、 という方程式を考えてみます。これは、z � 0の場合に真になります。

z = -1 の場合は、 、および となります。

したがって、 となります。この問題は、平方根が複数の値となり、複素数平面全体に連続する値を選択できないという事実に由来するものです。ただし、すべて負の実数で構成されるブランチカットを考慮の対象から外せば、平方根は連続になります。これでも、-x + i0 (ただし x > 0) の形式の負の実数について、どう扱うかという問題は残ります。符号付きゼロにより、この問題を完全な形で解決することができます。
x + i(+0) の形式の数は符号と 、またブランチカットのもう一方にある x + i(-0) の形式の数は、逆の符号と になります。実際に、 を計算する公式により、次の結果が得られます。

となるので、z = -1 + i0 であれば、
1/z = 1/(-1 + i0) = [(-1 - i0)] / [(-1 + i0)(-1 - i0)] = (-1 - i0)/((-1)2 - 02 ) = -1 + i(-0)
となり、
、一方 となります。
したがって、IEEE 演算では、すべての z の意味が保持されます。さらに高度な例が Kahan によって紹介されています (1987)。+0 と -0 を区別することに利点はありますが、わかりにくくなることがあります。例えば、符号付きゼロにより、x = +0、および y = -0 のときに偽となる x = y ⇔ 1/x = 1/y の関係が壊れます。ただし、IEEE 委員会では、符号付きゼロを使用した場合の利点の方が欠点を上回っていると判断しています。

非正規化数

β= 10、 p = 3、e_min = -98 で正規化された浮動小数点数の例を考えてみます。
x = 6.87 × 10^-97 と y = 6.81 × 10^-97 は、最小の浮動小数点数 1.00 × 10^-98 の 10 倍以上大きいので、通常の浮動小数点数のように思われます。しかし、これらの浮動小数点数には、x ≠ y! であっても、x y = 0 になるという変わった特性があります。
これは、x - y = 0.06 × 10^-97 = 6.0 × 10^-99 が、正規化数としては小さすぎて表現できないため、ゼロにフラッシュしなければならないという理由によるものです。この特性を保持することは、どの程度重要なのでしょう。

x = y ⇔ x - y = 0 ? (10)

if (x ≠ y) then z = 1/(x-y)というコード・フラグメントを記述し、後で、見かけ上のゼロ除算によってプログラムをアボートさせることも簡単にできます。しかし、このようなバグ検出の方法は、心理的にも時間的にも大きな負担になります。さらに高度な方法として、コンピュータサイエンスの参考書には、現在のところ大規模なプログラムの妥当性を検証するのは現実的ではないものの、これを実証することを念頭にプログラムを設計すると、効率的なコードができあがると指摘されています。例えば、不変式を取り入れるのは、仮にこれが証明の一部としては利用できない場合でも、効果的な手段です。浮動小数点コードによって、通常のコードの場合と同じように、信頼できる証明可能な事実が得られるようになります。
例えば、公式 (6) を分析する場合、x/2 < y < 2x x y = x - y という事実がわかれば非常に役立ちます。同じように、公式 (10) が真であることがわかれば、信頼できる浮動小数点コードも記述しやすくなります。大部分の数にあてはまるからという前提は、何も証明することにはなりません。

IEEE 標準では、非正規化¹⁸数を使用して公式 (10) 、およびその他の関係を証明します。非正規化数はよく議論の対象となる部分であり、754 で認可されるまでに長い期間がかかったのは、おそらくこの理由によるものと思われます。パフォーマンスの高いハードウェアはほとんどの場合、IEEE 互換であることがすなわち、非正規化数を直接サポートすることとは限りません。非正規化数を使用、あるいは生成する時点でハードウェアがトラップされ、ソフトウェアに IEEE 標準のシミュレーションを一任する形式になっています¹⁹。非正規化数の背景にある考えは、Goldberg (1967) が提案したもので、きわめて単純です。指数の e_min の場合、有意桁は正規化する必要がないので、β= 10、 p = 3、e_min = -98 であれば、0.98 × 10^-98 も浮動小数点数になるので、1.00 × 10^-98 は最小の浮動小数点数ではなくなります。

β= 2 の場合に隠れビットを使用すると、e_min の指数を伴う数は、暗黙に先行するビットがあるために必ず有意桁が 1.0 以上になるので、多少の支障が出てきます。この解決方法は、0 を表わす場合のそれと同じで、表 D-2 にまとめてあります。e_min 指数を使用して、非正規化数を表わすことができます。すなわち有意桁フィールドのビットが b ₁, b ₂, ... b p _-1で、指数の値が e の場合、表現する数は 1.b ₁b ₂ ... b_{p -1} × 2^e となり、e = e_min -1 の場合は 0·b ₁b ₂ ... b p -1 × 2^e+1 となります。非正規化数には e_min -1 ではなく、e_min の指数を伴うので、指数の +1 は必須になります。

本項の最初に示したβ= 10、 p =3、e_min = -98、x = 6.87 × 10^-97、y = 6.81 × 10^-97 の例を再び取り上げます。非正規化数の場合、x - y はゼロへのフラッシュは行われず、代わりに 0·6 × 10^-98 という非正規化数で表わされます。この動作のことを段階的な「アンダーフロー」と言います。段階的アンダーフローを使用すれば、公式 (10) がつねに真になることを簡単に検証することができます。

　　　 　　

図 D-2 段階的アンダーフローとゼロフラッシュの比較

図 D-2 に、非正規化数を示します。上の数線は、正規化した浮動小数点数を表わします。0 から最小の正規化数 1.0 × までに格差があることに注意してください。浮動小数点演算の結果がこの差の範囲に入る場合は、ゼロへのフラッシュが行われます。下の数線は、一連の浮動小数点数に非正規化数を加えたときに、どう処理されるかを示したものです。この格差が埋められ、計算結果が 1.0 × より少なければ、最も近い非正規化数で表わされます。数線に非正規化数を加算すると、隣接する浮動小数点数どうしの間隔が規則的に取られます。隣接する数どうしの間隔は同じ場合と、βの係数によって異なる場合とがあります。非正規化数がなければ、β^-p+1 から までの間隔が、βの係数による規則的な変化ではなく、β^p-1 の係数に従って、急激に変化することになります。このため、アンダーフローの限界に近い正規化数に対して大きな相対誤差を伴うアルゴリズムはほとんどの場合、段階的アンダーフローを使用すると、この範囲内で正しく動作するようになります。

段階的アンダーフローがなければ、x - y のような単純な式でも、上の
x = 6.87 × 10^-97、および y = 6.81 × 10^-97 の例のように、正規化された入力に対して非常に大きな相対誤差を伴うことになります。大きな相対誤差は、次の例 (Demmel、1984) に示すように、相殺がない場合でも起きる可能性があります。a + ib と c + id という 2 つの複素数を除算する例を考えてみます。次の公式には、分母 c + id のいずれかの要素が より大きいときに、最終結果が範囲内にあっても、オーバーフローするという問題があります。

� i

この商を計算する場合、次の Smith の公式を使用すると、より効果的です。

(11)

Smith の公式を (2・10^-98 + i10^-98 ) / (4・10^-98 + i(2・10^-98)) に適用すると、段階的アンダーフローによって 0.5 という正しい解が得られます。ゼロ・フラッシュの場合は、
0.4 になり、誤差は 100 ulp となります。これは、1.0 × を下限とする引数の誤差範囲を保証する非正規化数としては典型的な値です。

例外、フラグ、トラップハンドラ

IEEE 演算で、ゼロによる除算、あるいはオーバーフローといった例外が起きると、デフォルトとして結果を報告し、処理を続行します。一般的なデフォルトとして、0/0 と には NaN、また 1/0 とオーバーフローには を返します。前の項では、例外が起きた場合でも、以上のデフォルト値によって処理を続行する方が望ましい例をいくつか紹介しました。いずれかの例外が起きると、ステータスフラグもセットされます。IEEE 標準の仕様では、ユーザーがステータスフラグを読み取り/書き込みできる手段を提供しなければなりません。フラグは、いったんセットされると、明示的にクリアーしない限り、そのままセット状態になるという意味で「スティッキ」な特性があります。本来の無限大を表わす 1/0 と、オーバーフローとの違いを区別するには、フラグをテストする以外に方法はありません。

例外が起きたときに、処理を続行することが妥当でない場合もあります。「無限大」では、x/ (x² +1) の例を示しました。

x > の場合、分母は無限大になるので、最終的な解が 0 という、まったく誤った結果になります。この公式の場合、1/ (x + x^-1) と書き換えれば問題は解決しますが、式の書き換えによって必ずしも問題が解決しない場合もあります。IEEE 標準では、トラップハンドラをインストールできるような仕様が強く推奨されています。例外が起きると、フラグをセットする代わりにトラップハンドラが呼び出され、トラップハンドラが返す値を操作の結果として使用します。ステータスフラグをクリアーまたはセットするという判断は、トラップハンドラに一任されます。これ以外の場合、フラグの値は未定義にすることができます。

IEEE 標準では、例外をその種類によって、オーバーフロー、アンダーフロー、ゼロによる除算、無効な操作、および不正確な操作という 5 つのクラスに分類しています。例外の各クラスごとに、個別のステータスフラグが用意されています。最初の 3 つの例外の意味はその名が示すとおりです。無効な操作とは、表 D-3 に示す状況と、NaN を伴うあらゆる比較操作のことです。無効操作例外の原因となる操作では、デフォルトとして NaN が返されます。ただし逆は真ではありません。ある操作に対するオペランドのいずれかが NaN の場合、結果も NaN になりますが、操作の内容が 表 D-3 に示す条件のいずれかに一致しない限り、無効操作の例外は起きません²⁰。

表 D-4   IEEE 754の例外¹

例外	トラップが無効な場合の結果	トラップハンドラの引数
オーバーフロー	± または±xmax	round(x2-a)
アンダーフロー	0、 または denormal	round(x2a)
ゼロによる除算		オペランド
無効な操作	NaN	オペランド
不正確な操作	round(x)	round(x)

1 x は操作の正確な結果、単精度は a =192、倍精度は a =1536、xmax = 1.11 ...11 × となります。

β= 10、 p = 3 のシステムで、3.5 4.2 = 14.7 は正確ですが、3.5 4.3 = 15.0 は正確ではないので (正しくは 3.5・4.3 = 15.05)、不正確操作の例外が起きます。「2 進数から 10 進数への変換」では、不正確操作の例外を使用したアルゴリズムについて説明しています。また表 D-4 に、5 種類の例外の動作をまとめてあります。

不正確な操作の例外が頻繁に起きる場合は、実装上の問題があるものと思われます。浮動小数点ハードウェアに独自のフラグが用意されておらず、代わりにオペレーティングシステムに割り込んで浮動小数点例外を報告する形式の場合は、不正確な操作の例外が起きたときのコストはきわめて高くなります。この問題は、ソフトウェアが保持するステータスフラグを使用すれば、解決します。例外が初めて起きた時点で、該当するクラスのソフトウェアフラグをセットし、浮動小数点ハードウェアに、対応する例外クラスをマスクするように指示します。これ以降の例外はすべて、オペレーティングシステムへの介入を伴わずに、処理されるようになります。ユーザーがステータスフラグをリセットすると、ハードウェアマスクが再度イネーブルになります。

トラップハンドラ

トラップハンドラを使用する意味は、まず下位の互換性を保持することにあります。例外が起きた時点でアボートとなる旧バージョンのコードでは、プロセスを打ち切るためのトラップハンドラをインストールすることができます。
これは特に、do S until (x � 100)のようなループを使用するコードに有効になります。NaN は、<、≦、>、 ≧ 、= (≠ はなし) を使用した数に比較されると必ず偽を返すので、x が NaN になると、このコードは無限ループに入ります。

オーバーフローする可能性のある のような積を計算するときに使用されるトラップハンドラをさらに別の用途として利用することもできます。対数を使用して、exp を計算する方法です。このアプローチの問題は、精度に欠けることに加え、オーバーフローが起きない場合でも単純な 式よりもコストが大きくなる点にあります。またこれらの問題を避けるためのオーバーフロー/アンダーフロー・カウンティングというトラップハンドラを使用する方法もあります (Sterbenz、1974)。

この考え方は次のとおりです。まずゼロに初期化される汎用カウンタがあります。 が特定の k に対してオーバーフローすると必ず、トラップハンドラがカウンタの値を 1 だけ増やし、指数をラップしてオーバーフローの値を返します。IEEE 754 の単精度では、e_max = 127 となるので、pk = 1.45 × 2¹³⁰ であれば、オーバーフローを起こし、トラップハンドラが呼び出されます。トラップハンドラは指数をもとの範囲にラップし、pk を 1.45 × 2^-62 (下記参照) に変更します。同様に、pk がアンダーフローすると、カウンタの値が減り、負の指数が正の指数にラップされます。乗算がすべて終わった時点でカウンタがゼロの場合、最終的な積は pn になります。カウンタが正であれば積はオーバーフローし、負であればアンダーフローします。部分積がいずれも範囲を超えていない場合は、トラップハンドラは呼び出されず、計算により余分なコストが発生することはありません。オーバーフローまたはアンダーフローが起きた場合でも、pk はそれぞれ、全面的な精度の乗算を使用して P_{k - 1} から計算されているので、対数で計算した場合よりもはるかに正確です。Barnett (1987) は、オーバーフロー/アンダーフロー・カウントの全面精度によって、前出の表で該当する誤差を検出できる公式を紹介しました。

IEEE 754 の仕様では、オーバーフロー/アンダーフロー・トラップハンドラが呼び出された時点で、これがラップの結果を引き数として渡されます。オーバーフローに対するラップは、結果を無限の精度で計算し、これを 2a で除算した後、対応する精度に丸める操作と定義されます。アンダーフローの場合は、結果を 2a で乗算します。指数α は単精度の場合は 192、倍精度の場合は 1536 となります。これは、1.45 × 2¹³⁰ が 1.45 × 2^-62 に変換されるからです。

丸めモード

IEEE 標準では、それぞれの操作を正確に計算し、その結果が丸められる (2 進数から 10 進数への変換は除く) ので、正確でない結果が生成されると必ず、丸めが行われます。丸めとは本来、最も近い値で表わすということです。これ以外に標準には、0 への丸め、+ への丸め、- への丸めという 3 種類の丸めモードがあります。- への丸めを整数操作への変換とともに使用すると、変換処理は切り下げ関数として行われます。一方 + に丸めると、切り上げ関数になります。0 への丸め、または - への丸めが有効になっているときに、正の値がオーバーフローすると、デフォルトの結果が + ではなく、表現可能な最大値になるので、丸めモードはオーバーフローに影響を与えることになります。同様に、+ への丸め、または 0 への丸めが有効になっているときに負の値がオーバーフローすると、負の最大値が生成されます。

丸めモードの 1 つの応用例として、インタバルの計算があります (もう 1 つの応用例は、「2 進数から 10 進数への変換」で示します)。インタバルの計算を使用する場合、2 つの値 x、y の和をインタバル として扱います。ここで は x y を- に向かって丸め、 は x y を + に向かって丸めるという意味です。加算の正確な結果は、インタバル の範囲にあります。丸めモードがなければ、インタバルの計算は通常、 と (ここで はマシン・イプシロン) によって実現します²¹。

この結果、インタバルのサイズが多めに見積もられます。インタバル計算の結果はインタバルとして計算されるので、通常、操作に対する入力もインタバルで指定します。2 つのインタバル と を加算すると、結果は (ここで は丸めモードを - にセットした 、また は丸めモードを + にセットした ) になります。

インタバル計算を使用して浮動小数点演算を行うと、最終的な解は、正確な計算結果を含むインタバルとして得られます。これは、インタバルが広すぎる場合 (よく起きることです)、正しい解がそのインタバル内の任意の値になるので、それほど便利ではありません。インタバルの計算は、多重精度の浮動小数点パッケージと併用すれば、さらに意味を持ちます。まず、特定の精度 p とともに計算を行います。インタバル計算の結果により、最終的な解が不正確であることがわかれば、最終的なインタバルが適切なサイズになるまで段階的に精度を上げながら、何度か計算しなおします。

フラグ

IEEE標準では、いくつかのフラグとモードがサポートされています。すでに説明したとおり、5 種類の例外 (アンダーフロー、オーバーフロー、ゼロによる除算、無効な操作、不正確な操作) に対してステータスフラグが用意されています。また丸めモードには、最も近い値への丸め、+ への丸め、0 への丸め、- への丸めの 4 種類があります。それぞれの例外に対してイネーブル・モード・ビットを 1 つ確保するように強くお勧めします。本項では、各モードとフラグをどのように有効に利用できるかを示す例をいくつか紹介します。さらに高度な例については、「2 進数から 10 進数への変換」に示しています。

xⁿ (ここで n は整数) を計算するサブルーチンの例を考えてみます。n > 0 の場合は、次のような単純なルーチンが使用できます。

PositivePower(x,n) { while (n is even) { x = x*x n = n/2 } u = x while (true) { n = n/2 if (n==0) return u x = x*x if (n is odd) u = u*x }

一方、n < 0 の場合に xⁿ を計算するさらに正確な方法は、PositivePower (1/x, -n)ではなく、1/PositivePower (x, -n) を呼び出すことです。この理由は、最初の式の場合、乗算を行う n の値ごとに、1/x の除算で得られた丸め誤差が含まれることになるからです。2 つ目の式では、正確な値の計算 (x) にもとづいて、最後に行う除算で丸め誤差が 1 回だけ伴います。しかし、この方法にも欠点があります。PositivePower (x, -n) がアンダーフローになると、アンダーフロー・トラップハンドラが呼び出されるか、あるいはアンダーフロー・ステータスフラグがセットされます。この動作は、x^-n がアンダーフローした時点で、xⁿ はオーバーフローになるか、範囲内に収まるかのいずれかになるので、誤りです²²。しかし、IEEE 標準ではユーザーがすべてのフラグにアクセスできるので、サブルーチンにより簡単にこれを修正することができます。単に、オーバーフロー/アンダーフロー・イネーブル・ビットをオフにして、オーバーフロー/アンダーフロー・ステータスビットを保存するだけのことです。この後、1/PositivePower (x, -n) を計算します。オーバーフロー、アンダーフローのいずれのステータスビットもセットされなければ、トラップ・イネーブル・ビットとともにステータスビットがリストアされます。いずれかのステータスビットがセットされた場合は、フラグがリストアされ、PositivePower (1/x, -n)を使用して再計算します。この場合は、適切な例外が起きます。

フラグを使用するもう 1 つのケースは、arccos x = 2 arctan という公式で arccosを計算するときに起きます。arctan ( ) の評価が π/2 となる場合、arccos (-1) は、無限大の演算に従って 2・arctan ( ) = πに正しく評価されます。ただし、 (1 - x)/(1 + x) の計算により、arccos (-1) は例外にならなくても、ゼロによる除算例外のフラグがセットされます。この問題は簡単に解決できます。まずゼロによる除算フラグの値を保存してから、arccos を計算し、計算が終了した後、もとの値をリストアします。

システムの側面

コンピュータシステムのあらゆる側面を設計する場合に、浮動小数点に関する知識が必要になります。コンピュータアーキテクチャには通常、浮動小数点命令が組み込まれているので、コンパイラでこれらの浮動小数点命令を生成すると同時に、オペレーティングシステムでは、各浮動小数点命令に対して例外条件が発生した場合の処置を判断しなければなりません。一般に、数値解析の解説書を読んでも、基本的な内容はソフトウェアのユーザーや開発者を対象に書かれているので、コンピュータシステムの設計者にはそれほど参考にはなりません。設計上の方針が、どの程度予期せぬ方向に発展するかを示す実例として、次の BASIC プログラムを考えてみます。

q = 3.0/7.0 if q = 3.0/7.0 then print "Equal": else print "Not Equal"

Borland 社の Turbo Basic を IBM PC 上で実行しながら、"Not Equal" と出力するプログラムです。この例の分析は、次の項で行います。

なお、上記のような特殊な例は浮動小数点数と同列に扱うことはできないという指摘もあると思われます。しかし、いずれの例も特定の範囲 E の中に収まるものという意味で、同じであると解釈してください。これは、多くの答えになると同時に、多くの疑問も呈するので、万能的な解答にはならないのは言うまでもありません。E の値はどのようになるのでしょう。x < 0 と y > 0 の値が E の範囲内にある場合、それぞれ異なる符号があるのに、実際にこれが等価であると言えるでしょうか。さらに、a ～ b ⇔|a - b| < E というルールで定義される関係は、a ～ b と b ～ c が a ～ c を意味するものではないので、等価の関係とみなすことはできません。

命令セット

あるアルゴリズムで、正確な結果を生成するために、精度の高いショートバーストが必要になることがあります。

一例として、二次方程式の根の公式 ( ) /2a があげられます。234ページに説明するとおり、b² 4acのとき、丸め誤差により二次方程式の根の公式で求めた根の桁の半分が汚染される可能性があります。b² - 4ac の部分計算を倍精度で行うことにより、根の倍精度ビットの半分 (すなわち単精度ビットの全部) が保持されます。

a、b、c の各値が単精度フォーマットのときに、2 つの単精度数を取り、倍精度の結果を生成する乗算命令があれば、b² - 4ac を倍精度で計算するのは簡単です。2 つの p 桁数を正確に丸めるには、処理の進行に従ってビットが破棄される可能性はあるものの、積の 2p ビット全体を生成するために乗数が必要になります。したがって、単精度オペランドから倍精度の積を計算するハードウェアは通常、単精度の乗数の場合に比べコストが高く、また倍精度の乗数よりもコストが低くなります。これにもかかわらず、最近の命令セットはオペランドと同じ精度の結果を生成する命令だけを提供する傾向があります²³。

2 つの単精度オペランドを組み合わせて倍精度の結果を生成する命令が、二次方程式の根の公式に限って適用される場合、命令セットに追加しても意味がありません。しかし、この命令にはさまざまな使い方があります。次の一次方程式の系を求める問題を考えてみます。

a₁₁x₁ + a₁₂x₂ +・・・+ a_1nx_n= b₁
a₂₁x₁ + a₂₂x₂ +・・・+ a_2nx_n= b₂
� � �
a_n1x₁ + a_n2x₂ +・・・+ a_nnx_n= b_n

これは、次のように Ax = b として行列形式に書き換えることができます。

x⁽¹⁾ の解をガウスの消去などの方法で計算すると仮定します。結果の精度を「反復的改善」の方法で上げることができます。

最初の計算を行います。

ξ= Ax⁽¹⁾ - b (12)

次に系を計算します。

Ay = ξ (13)

x⁽¹⁾ が正確な計算である場合、ξは y と同様にゼロ・ベクタになります。一般に、ξと y の計算は丸め誤差を招きやすいので、Ay ξ Ax⁽¹⁾-b = A (x⁽¹⁾- x) となります (ただし x は真の解 (未知) です)。y x⁽¹⁾- x となるので、解の概算は次のようになります。

x⁽²⁾ = x⁽¹⁾ - y (14)

(12)、 (13)、 (14) のステップをそれぞれ x⁽¹⁾ を x⁽²⁾ に、また x⁽²⁾ を x⁽³⁾ に置き換えて繰り返すことができます。x⁽ⁱ⁺¹⁾ が x⁽i⁾ より正確であるという指摘は正式なものではありません。詳細については、Golub/Van Loan (1989) を参照してください。

反復的改善を実行すると、ξは、近似値ではあるものの正確には一致しない浮動小数点数の差に相当する要素を含むベクタとなるので、悪性の相殺の問題が出てきます。したがって、反復的改善は、ξ= Ax⁽¹⁾ - b を倍精度で計算しない限り、それほど意味がありません。これも、完全な倍精度の結果が必要とされるところで、2 つの単精度数 (A と x⁽¹⁾) の積を計算するケースに相当します。

要約すると、2 つの浮動小数点数を乗算し、オペランドの 2 倍の精度で積を返す命令は、浮動小数点命令セットに追加することができます。このことがコンパイラに与える影響について次に説明します。

言語とコンパイラ

コンパイラと浮動小数点の相互関係については、Farnum (1988) の資料に詳しく説明されています。本項の大部分は、この資料から抜粋したものです。

あいまいさ

言語の定義によって、プログラムに関する文が証明されるように言語の意味を正確に規定できるのが理想的です。このことは、言語の整数部分についてはあてはまるものの、浮動小数点に関する限り、言語の定義自体があいまいであることも事実です。おそらく、言語の設計者の意識の中に、浮動小数点には丸め誤差が付きものであるため、証明できるものは何もないという思い込みがあるように思われます。仮にこれが本当であるなら、これまでの説明は、この論拠の誤った認識について実証してきたことになります。本項では、言語定義に共通して見られるあいまいさについて、その対処方法などを中心に説明します。

驚くべきことに、言語によっては x が浮動小数点の変数 (例えば 3.0/10.0) である場合、10.0*x が現れるたびにすべて同じ値になるという明確な指定のないものがよくあります。例えば、Brown のモデルにもとづく Ada では、浮動小数点演算は Brown の原理だけを満たしていれば十分なので、式が多数の値を持っていてもかまわないことになります。浮動小数点をあいまいな形で捉えるという点は、浮動小数点の演算ごとに結果を正確に定義する必要のある IEEE モデルとはまったく異なるところです。
IEEE モデルでは、(3.0/10.0) *10.0 が 3 になることを証明 (定理 7) できますが、Brown のモデルではこれを証明できません。

大部分の言語定義に見られるもう 1 つのあいまいさは、オーバーフロー、アンダーフロー、その他の例外が起きたときの処理方法に関するものです。IEEE 標準では、例外の動作が正確に指定されるので、この標準をモデルとして使用した言語では、この点に関するあいまいさの問題を回避できます。

かっこの解釈の仕方もあいまいです。浮動小数点数は丸め誤差を伴うため、代数の結合法則は必ずしも浮動小数点数には適用できません。例えば、x = 10³⁰、y = -10³⁰、
z = 1 とした場合の (x+y)+z と x+(y+z) の解はまったく異なります。前者の場合は 1 で、後者は 0 になります。かっこを残すことの重要性は決して強調しすぎることはありません。定理 3、4 および 6 に示したアルゴリズムはすべて、かっこの位置に依存します。例えば、定理 6 の公式 x_h = mx - (mx - x) は、もしかっこがなければ x_h = x になり、アルゴリズム全体が無意味になります。かっこの意味を尊重する必要のない言語の定義は、浮動小数点の計算にはほとんど無意味です。

部分式の評価が正確に定義されない言語も数多くあります。仮に ds が倍精度、x と y が単精度であると仮定します。ds + x*yという式で、積は単精度と倍精度のどちらで計算されるのでしょう。また別の例で、x + m/n (ここで m と n はともに整数 ) の式は整数と浮動小数点数のいずれで計算されるのでしょう。この問題には 2 つの方法で対処できますが、いずれも完全な解決策にはなりません。1 つは、式の中の変数をすべて同じタイプとして扱う方法です。これは簡単な解決方法ですが、欠点もいくつかあります。まず、Pascal のような言語では、部分範囲のタイプによって部分範囲変数と整数の変数を併用できるので、単精度変数と倍精度変数の混用を禁止してもほとんど意味がありません。また、定数に関する問題もあります。0.1*x という式の場合、ほとんどの言語では 0.1 が単精度の定数と解釈されます。ここで仮に、プログラマが浮動小数点変数の宣言をすべて、単精度から倍精度に変更しようと考えたとします。0.1 をそのまま単精度定数として扱うと、コンパイル時にエラーになります。この場合、プログラマは浮動小数点定数をすべて探して、これを変更する必要があります。

2 つ目のアプローチは、式の混用を認めることにより、部分式の評価に関するルールを提供する方法です。これには参考になる例がいくつかあります。C 言語本来の定義では、浮動小数点式はすべて倍精度で計算しなければなりません (Kernighan/Ritchie、1978)。この条件は、本項の最初に紹介した例のような偏差をもたらす原因となります。
3.0/7.0 は倍精度で計算されますが、q が単精度変数の場合、商は単精度に丸めて格納されます。3/7 は 2 進の循環小数になるので、倍精度による計算値は、単精度で格納された値とは一致しなくなります。したがって、q = 3/7の計算は失敗します。これは、すべての式を最も高い精度で計算しても必ずしも効果的ではないことを示すものです。

もう 1 つの例は内積に関するものです。内積に多数の項が含まれる場合、和の丸め誤差が大きくなる可能性があります。この丸め誤差を減少させる方法として、倍精度の和を累積することができます (この方法については、「最適化プログラム」で詳しく説明しています)。d が倍精度変数で、x[] と y[] が単精度配列の場合、内積のループは d = d + x[i]*y[i] のようになります。単精度で乗算を行うと、倍精度変数に追加する直前に積が単精度に切り捨てられるので、倍精度による累積の利点が得られません。

前の 2 つの例に適用されるルールは、ある式のあらゆる変数中で最も高い精度で計算するというものです。そうすると、 q = 3.0/7.0 全体が単精度で計算され²⁴、真の論理値が与られます。一方、d = d + x[i]*y[i] は倍精度で計算されるので、倍精度の累積による利点が全面的に生かされます。ただし、このルールはあらゆるケースにそのまま適用できるわけではありません。dx と dy が倍精度変数の場合、y = x + single (dx-dy) の式には倍精度の変数が含まれますが、オペランドが両方とも単精度であることに加え、結果も単精度になるので、倍精度で和を求めても意味がありません。

さらに複雑な部分式に関する評価ルールは次のとおりです。まず、それぞれの操作に対し、仮の精度として対応するオペランドの最も高い精度を割り当てます。この割り当ては、式のツリー構造で葉の部分から根の方向に行います。次に、根から葉の方向に二次パスを設定します。このパスでは、各操作に最高の精度と、親から要求される精度を仮に割り当てます。q = 3.0/7.0 の場合、葉の部分は単精度になるので、操作はすべて単精度で実行されます。d = d + x[i]*y[i] の場合は、乗算の一時的な精度が単精度になりますが、二次パスでは親の操作から倍精度オペランドが要求されるので倍精度に昇格します。y = x + single (dx-dy) では、加算が単精度で実行されます。Farnum (1988) の資料では、このアルゴリズムのインプリメントがそれほど難しくないことが実証されています。

このルールの欠点は、部分式の評価が、中に組み込まれた式に依存する点にあります。これは、複雑なシーケンスの原因になることがあります。例えば、プログラムのデバッグ作業を行なっているときに、部分式の値を確認したい場合があります。この場合、プログラム中の部分式は、これが組み込まれた式に依存するので、単にデバッガに部分式を入力して評価するように指示することはできません。最後に、部分式には次の問題があります。10 進数定数から 2 進数定数への変換も 1 つの演算となるので、評価ルールも 10 進数定数の解釈に影響を与えるということです。これは 0.1 のように 2 進数では正確に表現できない定数の場合に特に重要になります。

その他、組み込まれている動作の 1 つとして言語に指数が含まれている場合にも、あいまいさが現れます。基本的な四則演算の場合とは異なり、指数の値は必ずしも明確には表現されない場合があります (Kahan/Coonen、1982)。** が指数関数の演算子である場合、(-3) **3 の値は -27 となります。ただし (-3.0) **3.0 などの場合は問題になります。** 演算子が整数のべき乗をチェックすると、(-3.0 ) **3.0 を
-3.0³ = -27 として計算します。一方、x^y = e^ylogx の公式を使用して実際の引数に ** を定義すると、log 関数によって結果が NaN (x < 0 のとき log(x) = NaN という定義にもとづく) になることがあります。これに FORTRAN の CLOG 関数を使用すると、ANSI FORTRAN の標準により、CLOG (-3.0) が iπ + log 3 (ANSI 1978) と定義されるので、解は -27 になります。Ada 言語では、整数のべき乗に限って指数を定義することにより、この問題を回避しています。一方、ANSI FORTRAN では、負数を実際の乗数でべき乗することはできません。

FORTRANの標準では、次のように規定されます。

「結果を数学的に定義できない算術演算は禁止される」

しかし、IEEE 標準に± という概念が導入されたことにより、「数学的に定義できない」という表現はまったく意味をなさなくなりました。1 つの定義として、「無限大」に示した方法を適用できます。たとえば、a^b の値を決定する場合に、
x → 0 に従ってそれぞれ f(x) → a、g(x) → b という特性を持つ非定数の解析関数 f と g の例を考えてみます。f(x)^g(x) はつねに同じ極限に近づくとすれば、これが a^b の値となります。この定義により、2^· = が成立します。1.0^· の場合、f(x) = 1、g(x) = 1/x のときに、限界は 1 に近づきますが、f(x) = 1 - x、g(x) = 1/x のときは限界が e^-1 となります。したがって、1.0^· は NaN となります。また 0 ⁰ の場合は、f(x)^g(x) = e^{g(x)log f(x)} となります。f と g は解析関数なので、0、f(x) = a ₁x¹ + a ₂x² + ...、および
g(x) = b ₁x¹ + b ₂x² + ... のとき値 0 を取ります。したがって、
lim_x � ₀g(x) log f(x) = lim_x � ₀x log (x(a ₁ + a ₂x + ...) ) = lim_x � ₀x log(a ₁x) = 0 になります。さらにすべての f と g に対して、f(x)g(x) → e⁰ = 1 となります。すなわち 0 ⁰ = 1 になります ^{25 26}。この定義を使用すると、すべての引数に対する指数関数の定義があいまいになります。特に (-3.0) **3.0 は -27 と定義されることになります。

IEEE 標準

先の「IEEE 標準」では、IEEE 標準の主な特長について説明しました。しかし IEEE 標準では、これらの特長をプログラミング言語から利用できるかについては何も保証されていません。したがって、この標準をサポートする浮動小数点ハードウェアと、C、Pascal、FORTRAN などのプログラミング言語の間に不一致が生じるのはよくあることです。IEEE の機能によっては、サブルーチン・ライブラリの呼び出しによってアクセスできるものもあります。例えば、IEEE 標準には、平方根は正確に丸めなければならないという条件があります。これに対応して、平方根関数がハードウェア上で直接インプリメントされていることがよくあります。この機能には、ライブラリにある平方根ルーチンを使用して簡単にアクセスすることができます。しかし、標準の特性によっては、サブルーチンでは簡単にインプリメントできないものもあります。例えば、コンピュータ言語ではほとんどの場合、浮動小数点のタイプは 2 種類程度しかサポートされません。一方、IEEE 標準では、4 種類の精度が用意されています (このうち、推奨されるのは単精度に拡張を加えた拡張単精度フォーマット、単精度、倍精度、拡張倍精度フォーマットです)。これ以外に、無限大の扱いという問題もあります。± を表わす定数は、サブルーチンから供給できます。ただし、一定の変数のイニシャライザとして定数式が要求される状況では、この定数が使用できなくなります。

さらに微妙な状況として、丸めモード、トラップ・イネーブル・ビット、トラップハンドラ、および例外フラグで構成される状態で計算を伴う操作があげられます。1 つのアプローチは、状態を読み取りまたは書き込みするためのサブルーチンを用意することです。さらに、新しい値を自動的にセットし、もとの値を返す操作を 1 回の操作で処理できる呼び出しも便利な場合があります。「フラグ」に説明したとおり、IEEE の状態を修正するための一般的なパターンは、ブロック、またはサブルーチンの範囲に限って、その状態を変更することです。これにより、プログラマ側には、ブロックの出口を逐一、探し出す負担と、状態が復元されたことを確認する作業が増えることになります。この場合、状態をブロックの範囲内に正しく設定できる言語は非常に便利です。Modeula-3 は、トラップハンドラにこの概念を適用した言語です (Nelson、1991)。

言語の中で IEEE 標準をインプリメントするときには、考慮すべき点がいくつかあります。あらゆる x について x - x = +0 となるので²⁷、(+0) - (+0) = +0 となります。
ただし、- (+0) = -0 となるので、-x は 0 - x として定義することはできません。NaN は、別の NaN の場合も含め、同じ値には決してならないので、x = x が必ずしも真にはならないという理由から、NaN を導入することは混乱の原因になります。IEEE で Isnan 関数を使用しないように推奨されている場合は、実際に、x ≠ x の式は NaN をテストする最も簡単な方法です。また NaN はほかの数との順序に関連性がないので、x ≦ y を非 x > y として定義することもできません。NaN を導入すると、浮動小数点数が部分的に順序付けられるので、<、=、>、または unordered のいずれかを返す compare 関数により、プログラマは比較処理を簡単に実行できるようになります。

IEEE 標準では、いずれかのオペランドが NaN の場合に、NaN を返す基本的な浮動小数点操作を定義してありますが、複合操作の場合には、この定義が必ずしも有効になるとは限りません。グラフをプロットするときに使用する倍率を求める場合、一連の値の最大値を計算することができます。この場合、最大値の操作で単に NaN を無視すれば、意味があります。

最後に、丸めの操作が問題になることがあります。IEEE 標準では、丸め操作が厳密に定義されており、これは丸めモードの現在の値に依存します。丸めモードは、型変換における暗黙の丸めや、言語中の明示的な round 関数と競合することがあります。すなわちIEEE の丸めを使用するプログラムでは、自然言語のプリミティブを使用できず、一方、言語のプリミティブは、しだいに普及の広がる IEEE マシンでインプリメントするには不十分ということになります。

最適化プログラム

コンパイラに関する解説書は、浮動小数点に関する問題を無視する傾向があるようです。例えば、Ahoら (1986) は、x/2.0 を x*0.5 に置き換えることにより、読者は
x/10.0 も 0.1*x に置き換えるものと解釈してしまいます。しかし、これらの式は、
0.1 を 2 進数で正確に表現することはできないので、2 進数マシンでは同じセマンティクスにはなりません。本書では、x*y - x*z を x*(y-z) と置き換えることをお勧めします。これは、y  z のときに、仮に 2 つの式がまったく異なる値を持っている場合でもあてはまります。最適化プログラムが言語の定義に違反してはならないという条件は、コードを最適化するときは、どのような代数単位元でも使用できるという記述を確かに正当化するものですが、その反面、浮動小数点のセマンティクスはそれほど重要ではないという印象も与えます。言語の標準でかっこの位置を尊重すべきかに関する指定の有無を問わず、すでに説明したとおり (x+y) + z と x +(y+z) の解がまったく異なることも考えられます。例えば、かっこの保持に深くかかわる問題は、次のコードに見ることができます。

eps = 1; do eps = 0.5*eps; while (eps + 1 > 1);

このコードは、マシン・イプシロンに関する概算を求めるものです。最適化コンパイラが eps + 1 > 1 ⇔ eps > 0 であることを認識すると、プログラムはまったく異なる動作を示すことになります。1 x が x (x e β^-p) を超えることがないように、最小の x を計算する代わりに、x/2 を 0 (x ) 方向に丸めを行う最大の x を計算します。この種の「最適化」を避けることは、最適化によって失われる便利なアルゴリズムを保持するという意味で重要になります。

積分、微分方程式の数値解といった数多くの問題は、多数の項を伴う和の計算を行うものです。加算操作を行うごとに、0.5 ulp もの誤差が介入する可能性があるので、数千もの項を伴う和の丸め誤差は非常に大きくなる場合があります。この問題を解決する簡単な方法は、部分的な加数を倍精度変数として格納し、各加算処理をそれぞれ倍精度で実行することです。計算を単精度で実行して、和を倍精度で実行することは、ほとんどのコンピュータシステムでサポートされています。ただし、計算がすでに倍精度で行われている場合に、精度を倍加するのはそれほど簡単ではありません。よく提案される 1 つの解決案として、各数値をソートして、小さい順に加算していく方法があります。しかし、和の精度を大幅に改善できる、はるかに効率的な方法があります。すなわち、次の方法です。

定理 8 (Kahan の加法公式)

次のアルゴリズムに従って、 を計算するものとする。

S = X[1]; C = 0; for j = 2 to N { Y = X[j] - C; T = S + Y; C = (T - S) - Y; S = T; }

S の計算値は、 (ここで ) になる。

固有の公式 を使用して和を計算すると、 (ただし |δ_j| < (n - j)e) になります。この結果は、 Kahan の加法公式の場合に比べると著しく改善されています。この単純な公式の中で、各加数は、ne もの摂動の影響は受けず、単に 2e で摂動されるだけです。詳細については、「加法の誤差」を参照してください。

浮動小数点演算が、代数法則に従っていることを前提とする最適化プログラムでは、結果が C = [T-S] - Y = [(S+Y)-S] - Y = 0 となり、アルゴリズムがまったく無意味になります。これらの例により、実数に適用される代数の単位元を浮動小数点変数を伴う式に使用する場合は、最適化プログラムを慎重に使用する必要があるということが言えます。

さらに、最適化プログラムは、定数を扱う場合にも浮動小数点コードのセマンティクスを変える可能性があります。1.0E-40*x の式では、暗黙のうちに 10 進数から 2 進数定数への変換が行われます。この定数は、2 進数で正確に表現することはできないので、不正確な操作の例外が起きます。また、式が単精度で評価された場合は、アンダーフロー・フラグがセットされます。定数が不正確である場合、2 進数への正確な変換は、IEEE 丸めモードの現在の値に依存します。したがって、コンパイル時に 1.0E-40 を 2 進数に変換する最適化プログラムにより、プログラムのセマンティクスが変更されることになります。ただし、27.5 のような定数は、最低の精度でも正確に表現できるので、コンパイル時に正しく変換されます。すなわち、つねに正確に表現されるので、例外も起きず、丸めモードの影響も受けないということです。コンパイル時に変換すべき定数には、例えば const pi = 3.14159265 のような定数宣言を使用する必要があります。

共通項の消去の場合も、最適化プログラムによって浮動小数点のセマンティクスを変えられる可能性があります。例えば、次のようなコードです。

C = A*B; RndMode = Up D = A*B;

A*B は一見、共通項のように思われますが、2 か所の評価位置でそれぞれ丸めモードが異なるので共通項としては扱いません。次の 3 つの例が考えられます。まず x = x は、論理定数 true で置き換えることはできません。これは、NaN が通常の浮動小数点数に比べて大小で表わすことはできないので、x が NaN のときに成立しない、
x = +0 の場合に -x = 0 - x が成立しない、また x < y は x ≧ y の逆ではない、という理由があるからです。

以上の例がある反面、浮動小数点コードに対して効果のある最適化もあります。まず、浮動小数点数に有効となる代数の単位元があります。IEEE 演算の例として x + y = y + x、 2 × x = x + x、1 × x = x、 0.5 × x = x/2 などがあります。ただし、これらの単純な単位元は、CDC や Cray のスーパーコンピュータ上では失敗します。そのほか、命令のスケジューリングやインライン・プロシージャ置換も、最適化によって効果の得られる例となります²⁸。

3 つ目の例として、dx = x*y (x と y はともに単精度の変数、dx は倍精度を表わします) の式を考えてみます。2 つの単精度数を乗算して倍精度数を出力するマシンでは、dx = x*yはそれぞれ該当する命令にマップされます。つまりオペランドを倍精度に変換してから、倍精度の乗算をするという一連の命令へのコンパイルは行われません。

コンパイラの設計者の中に、(x + y) + z から x + (y + z) への変換を禁止する制限は、移植不能なトリックを使用するプログラマにしか関連のないことなので、無意味であるという指摘もあります。おそらく、これは浮動小数点数が実数のモデルに従っているので、実数と同じ法則に従うべきであるという論拠にもとづくものです。実数のセマンティクスに伴う問題は、インプリメントするのにコストがかかるという点です。2 つの n ビット数を乗算すると、積は 2n ビットになります。

間隔の広い指数を伴う 2 つの n ビット数を加算するごとに、和のビット数は n + 指数間の間隔で表わされます。和は (e^max - e^min) + n ビットとなり、およそ 2・e^max + n ビットに相当します。何千回もの演算を行うアルゴリズム (例えば、線形代数の解など) では、多数の有意ビットを含む数を対象として操作を行うので、処理がきわめて低速になります。また、sin や cos などの超越関数の値は有理数にはならないので、これらのライブラリ関数をインプリメントするのはさらに難しくなります。LISP システムでは、正確な整数の演算が提供されていて、問題を解決しやすい場合があります。ただし、正確な浮動小数点演算はほとんど役に立ちません。

実際には、(x + y) + z ≠ x + (y + z) という事実を利用した便利なアルゴリズム (Kahan の加法公式など) が用意されており、次の範囲が有効である限り、正しく動作します。

a b = (a + b)(1 +δ)

この範囲は、市販されている大部分のハードウェアに適用されるので、数値プログラマがこうしたアルゴリズムを無視したり、コンパイラ設計者が浮動小数点変数に実数のセマンティクスが保持されることを前提として、これらのアルゴリズムを利用しないことは、賢明とは言えません。

例外処理

ここまでの説明は主に、システムが持つ精度の意味に関するものでした。トラップハンドラも、システムに関する興味深い問題を提起するものです。IEEE 標準では、ユーザーが 5 つのクラスの各例外に対するトラップハンドラを指定できるような仕様が強く推奨されています。「トラップハンドラ」では、ユーザーが定義できるトラップハンドラのアプリケーションをいくつか紹介しました。無効な操作やゼロによる除算が行われた場合、ハンドラには該当するオペランド、あるいは正確に丸めた結果を提供する必要があります。使用するプログラミング言語によってトラップハンドラは、プログラム中のその他の変数にもアクセスできる場合があります。すべての例外について、トラップハンドラは、実行された操作の内容、および意図された精度を識別できなければなりません。

IEEE 標準では、各操作が概念的にシリアル形式で構成され、何らかの例外が起きた時点で操作の内容とオペランドを識別できることを前提としています。パイプライン機能や多重演算ユニットをサポートしたマシンでは、例外が起きた場合に、単にトラップハンドラでプログラムカウンタをチェックさせるだけでは不十分になる場合もあります。すなわち操作の内容を正確に識別できるハードウェアが必要になるということです。

次のプログラム部分は、また別の問題を示したものです。

x = y*z; z = x*w; a = b + c; d = a/x;

2 つ目の乗算により例外が起きて、トラップハンドラが a の値を使用したいものと仮定します。加算と乗算を並列的に実行できるハードウェアでは、最適化プログラムによって、2 つ目の乗算を始める前に加算操作が移動されるので、この加算処理は最初の乗算と並行して実行することができます。したがって、2 つ目の乗算がトラップされても、a = b + c の実行が終わっているので、a の結果が一部変更されている可能性があります。あらゆる浮動小数点演算がトラップされる可能性があり、この結果、最適化をスケジュールする命令がすべて消去されることがあるため、コンパイラで、以上のような最適化を避けることはよい方法とは言えません。この問題は、トラップハンドラがプログラム中の変数に直接アクセスするのをすべて禁止することによって解決できます。代わりに、ハンドラには引数として、オペランドまたは結果を渡すことができます。

これでもなお、別の問題が残ります。次のコード部分の場合、2 つの命令が並列的に実行される可能性があります。

x = y*z; z = a + b;

乗算がトラップされると、その引数 z がすでに次の加算によって上書きされている可能性があります。特に、加算は乗算よりも高速で処理されるので、これがあてはまります。IEEE 標準をサポートするコンピュータシステムでは、z の値をハードウェアに保存するか、あるいは最初の時点でコンパイラにこのような状況を回避させるか、といった何らかの手段が必要になります。

Kahanは、トラップハンドラの代わりに、「前置換」という方法によってこうした問題を回避できると提案しています。この方法では、ある例外が起きた時点で、その結果として使用する値と例外の内容をユーザー自身が指定します。例えば、(sin x)/x を計算するコードで、ユーザーが x = 0 となることはほとんどないと判断した上で、x = 0 のテストを省略し、代わりに 0/0 トラップが起きた時点で、このケースに対応させることによってパフォーマンスを改善したいとします。IEEE トラップハンドラを使用する場合、ユーザーは、sin x/x の計算を始める前に、まず値 1 を返すハンドラを作成し、これをインストールしなければなりません。一方、前置換を利用すれば、無効な操作が起きた時点で、値 1 を使用するように指定するだけで済みます (例外が実際に起きる前に、使用する値を指定しておく必要があるからです)。トラップハンドラの場合、返される値はトラップが起きた時点で計算することができます。

前置換の利点は、ハードウェアの実装がわかりやすくなるという点です²⁹。例外のタイプが特定された時点で、これをもとに、必要な操作の結果を保存したテーブルを参照させることができます。前置換にはいくつかの利点がある一方、IEEE の標準がすでに広く普及していることが、ハードウェアメーカーによる受け入れを阻んでいるという問題があります。

詳細

ここまで、浮動小数点演算に関する主な特性について、さまざまな観点から検討してきました。この後の項では、浮動小数点が決して難題ではなく、これまでに指摘されてきた問題はすべて数学的に検証できるわかりやすい問題であるという点を実証していきます。この項は大きく 3 つに分けることができます。第 1 部では、誤差の分析に関する導入説明を行います。また 177ページに示した「丸め誤差」について詳しく説明します。第 2 部では、2 進数から 10 進数への変換について、「IEEE 標準」に関する補足説明を行います。第 3 部では、「システムの側面」で例示した「Kahanの加法公式」について詳しく説明します。

丸め誤差

「丸め誤差」の項では、保護桁を 1 桁確保すれば、加算と減算は必ず正確に実行されると説明しました (定理 2)。ここでは、その検証を行います。定理 2 は、加算に関するものと、減算に関するものの 2 つの部分で構成されます。減算に関する定理は次のとおりです。

定理 9

x と y が、パラメータβ と p を伴うフォーマットで正の浮動小数点数を表わす場合、p + 1 桁 (1 桁は保護桁) で減算を行うと、結果の相対丸め誤差は

e ≦ 2e 以下になる。

証明

x > y になるように、必要に応じて x と y を入れ換えます。x が x₀.x₁ ... xp-1 × b⁰ になるように x と y を基準化 (スケール) してもかまいません。y が y₀.y ₁ ...y_p-1'で表わされる場合、その差は正確です。また y が 0.y ₁ ... yp の場合、保護桁によって計算値の差は、丸め誤差が最大 e以内になるように、浮動小数点数に丸めた正確な差を表わすことになります。
一般に、y = 0.0 ... 0y_k+1 ... y_k+p、および は、p + 1 桁に切り捨てた y となります。
これにより、次の式が成り立ちます。

y - < (β - 1)(β^{-p - 1} + β^{-p - 2} +...+ β^{-p - k}) (15)

保護桁の定義から、x - y の計算値は x - を浮動小数点数に丸めた値、すなわち (x - ) + δ(丸め誤差 δ は次の式を満たします) になります。

|δ| ≦ (β/2)β^-p (16)

正確な差は、x - y なので、誤差は (x - y) - (x - + δ) = - y +δになります。この場合、次の 3 つのケースが考えられます。x - y ≧ 1 の場合、相対誤差は次の範囲になります。

≦β^-p [(β- 1)(β^-1 +...+β^-k)+β/2] <β^-p(1+β/2) (17)

x - < 1 の場合は、δ= 0 になります。x - y が取れる最小値は次のとおりです。

> (β- 1)(β^-1 +...+ β^-k), (ただしρ= β - 1)

したがって、この場合の相対誤差は次の範囲になります。

(18)

最後は、x - y < 1 のときに x - ≧ 1 となるケースです。これは、x - = 1 の場合 (すなわちδ= 0) に限って起きます。δ= 0 の場合は、公式 (18) が適用されるので、相対誤差の範囲は β-p < β-p (1 + β/2) となります。■

β= 2 の場合、範囲は正確に 2eとなり、この範囲は p → とした x = 1 + 2^2-p と
y = 2^1-p - 2^1-2p に対して得られます。同じ符号の数を加算する場合は、次の例に見るように、保護桁がなくても、正確な結果が得られます。

定理 10

x ≧ 0、y ≧ 0 の場合、x + y を計算するときの相対誤差は、保護桁がない場合でも、最大 2εに抑えられる。

証明

k の保護桁に対するアルゴリズムは、減算のアルゴリズムに似ています。x ≧ y の場合は、x と y の基数の位置が合うまで、y を右にシフトします。p + k の位置を超えた桁は無視されます。この p + k 桁の 2 つの数の和を求めます。次に p 桁に丸めます。

保護桁を使用しない場合の定理について検証します (一般的なケースは同じです)。x ≧ y ≧ 0、および x を d.dd... × β⁰ の形式に基準化 (スケール) することを前提とした一般論が失われることはありません。まず、キャリーアウトはないものとします。y の終わりからシフトオフした桁の値は β^-p+1 より少なく、和は 1 以上となるので、相対誤差は β^-p+1 /1 = 2e となります。

一方、キャリーアウトを伴う場合、シフトによる誤差は β-p⁺² の丸め誤差に追加する必要があります。

和は b 以上になるので、丸め誤差は

≦ 2ε

以下に抑えられます。■

これら 2 つの定理を組み合わせることにより、定理 2 を導くことができます。定理 2 は、1 回の操作を実行する場合に伴う相対誤差を求めるものです。x² -y² と (x + y) (x - y) の丸め誤差を比較する場合、多重操作に伴う相対誤差を考慮しなければなりません。x y の相対誤差は、δ₁= [(x y) - (x - y)]/(x - y) となり、|δ₁|≦ 2e の条件を満たします。また、次のように書くこともできます。

x y = (x - y) (1 + δ₁), |δ₁| ≦ 2e (19)

同様に、

x y = (x + y) (1 + δ₂), |δ₂|≦ 2e (20)

まず正確な積の計算と丸めによって乗算を行うとすると、相対誤差は最大 0.5 ulp となるので、浮動小数点数 u と v には次の式が成り立ちます。

u v = uv (1 + δ₃), |δ₃| ≦ e (21)

上の 3 つの方程式を組み合わせる (u = x y、および v = x y) と、次のようになります。

(x y) (x y) = (x - y) (1+δ₁) (x + y) (1 +δ₂) (1 +δ₃) (22)

したがって、 (x - y) (x + y) を計算するときに起きる相対誤差は次のようになります。

相対誤差は、δ₁+δ₂+δ₃+δ₁δ₂+δ₁δ₃+δ ₂δ ₃+δ ₁δ ₂δ ₃ になり、範囲は 5ε+ 8ε² となります。すなわち最大相対誤差は 5 (丸め誤差) となります (e が最小数で、e² はほぼ無視できます)。

(x x) (y y) を同じように分析しても相対誤差は小さい値にはなりません。これは、x と y という 2 つの近い値が x² - y² に代入されると、相対誤差は通常、非常に大きくなるためです。この場合、 (x y) (x y) で真になった分析を再び使用して、次の式を得ることができます。

(x x) (y y) = [x²(1 + δ₁) - y²(1 + δ₂)] (1 + δ₃)
= ((x² - y²) (1 + δ₁) + (δ₁ - δ₂)y2) (1 + δ₃)

x と y の値が近い場合、誤差の項 (δ₁-δ₂) y²は、最大 x² - y² と同じくらい大きくなることがあります。以上の計算によって、x² - y² より (x - y) (x + y) のほうがより正確であるということが実証されたことになります。

次に、三角形の面積に関する公式を分析します。公式 (7) を計算するときに起きる最大誤差を見積もるには、次の事実が必要になります。

定理 11

保護桁で減算を行うと、y/2 ≦ x ≦ 2y の場合に、x - y は正確に計算される。

証明

x と y に同じ指数がある場合、x y は正確に計算されます。これ以外の場合、定理の条件から、指数の差は最大 1 であると言えます。0 ≦ y ≦ x になるように、必要に応じて x と y を入れ換えて、x が x₀.x₁ ...x_{p - 1}、および y が 0.y ₁ ... y_p となるように x と y を基準化 (スケール) してもかまいません。x y を計算するアルゴリズムは、まず x - y を正確に計算して、次に浮動小数点に丸めるものです。差が 0.d ₁ ... d_p の形式の場合、この差の長さはすでに p 桁になっているので、丸め操作は必要ありません。x ≦ 2y となるので、x - y ≦ y、また y は 0.d₁ ... d_p となるので、x - y になります。■

β > 2 の場合、定理 11 を y/β≦ x ≦ βy に置き換えることはできません。y/2≦ x ≦ 2y の強い条件が必要になります。定理 10 の証明の直後に示した (x - y) (x + y) の誤差分析は、加算と減算の基本演算中に起きる相対誤差は小さい (定理 (19)と(20)) という事実を前提にしたものです。これは、最も一般的な誤差分析の例です。公式 (7) を分析するには、次の証明で示すとおり、定理 11 のような条件が必要になります。

定理 12

保護桁で減算を行う場合、a、b、c が三角形の 3 辺 (a ≧ b ≧ c) とすると、
(a +(b + c))(c -(a - b))(c +(a - b))(a +(b - c)) を計算するときの相対誤差は最大 16ε(ただし e < 0.005) になる。

証明

各項を 1 つずつ検証してみます。定理 10 より、b c = (b + c ) (1 +δ₁) (ただしδ₁ は相対誤差、|δ₁| ≦ 2ε) を導くことができます。最初の項の値は
(a (b c)) = (a +(b c)) (1 +δ₂) = (a +(b + c) (1 +δ₁)) (1 +δ₂) になり、次のようになります。

(a + b + c) (1 - 2ε)² ≦ [a + (b + c) (1 - 2ε)] � (1-2ε)
≦ a (b c)
≦ [a + (b + c) (1 + 2ε)] (1 + 2ε)
≦ (a + b + c) (1 + 2ε)²

すなわち、η₁ が 1 つなので、次の公式が成り立ちます。

(a (b c)) = (a + b + c) (1 + η₁)², |η₁| ≦ 2ε (24)

次の項には、a b に丸め誤差を含むことがあるので、c と a b に悪性の減算を伴う可能性があります。a、b、c は三角形の 3 辺 (a ≦ b + c ) をなし、この式を
c ≦ b ≦ a の順序と組み合わせると、a ≦ b + c ≦ 2b ≦ 2a となります。したがって、
a - b は定理 11 の条件を満たします。すなわち a - b = a b は正しいことになり、c (a - b) という減算は無害なので、定理 9 から次のように導くことができます。

(c (a b)) = (c - (a - b)) (1 +η₂), |η2| ≦ 2ε (25)

3 つ目の項は、2 つの正確な正数の和を計算するものなので、次のようになります。

(c (a b)) = (c + (a - b)) (1 + η₃), |η₃| ≦ 2ε (26)

最後の項は、定理 9 と 10 にもとづいて次のようになります。

(a (b c)) = (a + (b - c)) (1 + η₄)², |η₄| ≦ 2ε (27)

x y = xy (1+ξ) (ただし|ξ | ≦ε) となるように、乗算を正確に丸めると、公式 (24)、 (25)、 (26)、 (27) の組み合わせは次のようになります。

(a (b c)) (c (a b)) (c (a b)) (a (b c))
≦(a + (b + c)) (c - (a - b)) (c + (a - b)) (a + (b - c)) E

ただし、

E = (1 + η₁)² (1 +η₂) (1 + η₃) (1 +η₄)² (1 + ζ₁)(1 + ζ₂) (1 + ζ₃)

E の上限は、 (1 + 2ε)⁶(1 +ε)³ となり、範囲は 1 + 15ε+ O (ε²) となります。作成者によっては O (e2) を単に無視する場合もありますが、これは簡単に説明できます。
(1 + 2ε)⁶(1 +ε) ³ = 1 + 15ε+εR(ε)と記述すると、R(ε) は正の係数を持つ e の多項式となるので、εの漸増関数になります。R (0.005 ) = 0.505 なので、ε< 0.005 にはすべて R(ε)< 1 となり、E ≦ (1 + 2ε)⁶(1+ε)³ < 1 + 16εが成り立ちます。
E の下限を求める場合、1 - 15ε - εR(ε ) < E となるので、ε< 0.005 のとき、
1 - 16ε< (1 - 2ε)⁶(1 - ε)³ になります。これら 2 つの範囲を組み合わせると、
1 - 16ε< E < 1 + 16εになります。したがって、相対誤差は最大 16εです。■

定理 12 は、公式 (7) に悪性の相殺は起きないことを証明するものです。ここで公式 (7) が数値的に安定していることを証明するまでもありませんが、公式全体の有効範囲 (「相殺」に示す定理 3 で決まる範囲) を規定しておくのがよい方法です。

定理 3 の証明

次のように仮定します。

q = (a + (b + c)) (c - (a - b)) (c + (a - b)) (a + (b - c))

および

Q = (a (b c)) (c (a b)) (c (a b)) (a (b c))

定理 12 により、Q = q (1 +δ) (ただしδ≦16ε) が成り立ちます。
δ≦0.04/ (0.52) ² 0.15と仮定すれば、次の式は、簡単にチェックできます。

(28)

また、|δ| ≦ 16ε ≦ 16 (0.005 ) = 0.08 となるので、δが条件を満たすことになります。したがって、|δ₁| ≦ 0.52|δ| ≦ 8.5εの場合、 となります。

誤差を 0.5 ulp 以内として平方根を計算すると、 を計算するときの誤差は、|δ₂| ≦εの場合に、(1 +δ₁) (1 +δ₂) になります。β= 2 の場合、4 で除算をした場合でもこれ以外の誤差は導入されません。その他の場合、除算にはもう 1 つの因数
1 +δ₃ (|δ₃| ≦ε)が必要で、また、定理 12 の証明に示した方法によって、
(1 +δ₁) (1 +δ₂) (1 +δ₃) の最終的な誤差の範囲は 1 +δ₄ に依存します
(ただし |δ₄| ≦ 11ε)。■

定理 4 の直後に示した帰納的な説明をさらに正確に表現できるように、次の定理は μ(x) がどれほど定数に接近するかを示したものです。

定理 13

μ(x) = ln (1+ x)/x とすると、0 ≦ x ≦ の場合、 ≦ μ(x) ≦ 1 になり、この微分は

|μ'(x)| ≦ を満たすことになる。

証明

μ(x) = 1 - x/2 + x²/3 - ...は下降項の交代級数なので、x ≦ 1 の場合、
μ(x) ≧ 1 - x/2 ≧ となります。μの級数は入れ代わってμ(x)≦ 1となることは簡単に確認できます。また μ′(x) のテイラー級数も入れ代わります。

x ≦ 3/4 は下降項なので、- ≦μ′(x) ≦ - + 2x/3、または - ≦ μ' (x) ≦ 0 となるので、|μ′(x) | ≦ となります。■

定理 4 の証明

ln に対するテイラー級数、

は交代級数、0 < x -ln (1 + x) < x²/2 なので、ln (1 + x) を x で近似させるときの相対誤差は x/2 により範囲が決まります。1 x = 1 の場合、|x|<εとなるので、相対誤差の範囲は ε/2となります。

1 x ≠1 の場合、1 x = 1 + によって を定義します。0≦ x < 1 のとき、
(1 x) 1 = です。

除算と対数を ulp 以内の範囲で計算した場合、ln (1 + x)/((1 + x) - 1) の式の計算値は次のようになります。

(1 + ₁) (1 + ₂) = (1 + ₁) (1 + ₂) = µ( ) (1 + ₁) (1 + ₂) (29)

ここで、|δ₁| ≦ ε 、|δ₂| ≦ ε です。μ( ) の概算を求めるには、x から の間の特定のξに対して、次の平均値の定理を使用します。

μ( ) - μ(x) = ( - x)μ′(ξ) (30)

の定義から、| - x| ≦εとなり、これを定理 13 と組み合わせると、
|μ( ) - μ(x)|≦ ε/2、または |μ( )/μ(x) - 1| ≦ε/(2|μ(x)|)≦εとなります。すなわち |δ₃|≦εの場合にμ( ) = μ(x)(1 +δ ₃) となります。最後に x で乗算すると、最終的に δ₄ が導入されるので、x・ln (1 x)/((1 x) 1) は次のようになります。

ε< 0.1 であれば、 (1 +δ₁) (1 +δ₂) (1 +δ₃) (1 +δ₄) = 1 +δ (ただし|δ|≦ 5ε の場合) となります。■

公式 (19)、 (20)、 (21) を使用した誤差分析の興味深い例は、二次方程式の根の公式 に見ることができます。「相殺」で、方程式を書き換えることにより、±の操作が原因で起きる可能性のある相殺をいかに除去できるかを示しました。しかし、d = b ₂ - 4ac の計算時に起きる可能性のある相殺がもう 1 つあります。この相殺は、単に公式を入れ換えるだけでは取り除くことはできません。簡単に言うと、b² 4ac のとき、丸め誤差により、二次方程式の根の公式で計算する根の半分の桁まで汚染される可能性があります。ここで、略式の証明を示します (二次方程式の根の公式の誤差を除去するための別のアプローチは Kahan (1972) の資料に紹介されています)。

b² 4ac の場合、丸め誤差により、二次方程式の根の公式 で計算する根の桁の半分まで汚染される可能性がある。

証明：

(b b) (4a c) = (b² (1 +δ₁) - 4ac (1 +δ₂)) (1+δ₃) (ただし |δ_i| ≦ε) となります。³⁰

d = b² - 4ac を使用すると、 (d(1 +δ₁) - 4ac(δ₂ - δ₁)) (1 +δ₃) のように再書き込みされます。

この誤差のサイズを見積もるには、δi の第 2 項を無視します。この場合の絶対誤差は
d (δ₁+δ₃) - 4acδ₄ となります (ただし |δ₄| = |δ₁-δ₂|≦ 2εです)。 であるため、第 1 項 d (δ₁+δ₃) は無視できます。
第 2 項を見積もる場合は、ax² + bx + c = a (x - r₁) (x - r₂) という事実にもとづき、
ar₁r₂ = c となります。b² 4ac の場合、r_{1 r2} となるので、第 2 の項は
4acδ₄ 4a² r ₁δ₄² となります。

したがって、 の計算値は となります。

次の不等式

は、 (ただし ) であることを示すので、

a の絶対誤差は約 となり、根 r₁ r₂ の下位半分を破壊します。

δ₄ β^-p なので、 となり、絶対誤差は になります。

すなわち、平方根の計算では の計算を伴い、この式では ri の下位半分に対応する位置に意味のあるビットが格納されないため、ri の下位ビットは有意になりません。■

最後に、定理 6 を検証してみます。これは、「定理 14 と定理 8」の項で証明する事実にもとづく説明です。

定理 14

0 < k < p として、m = βk + 1 をセットし、浮動小数点演算が正確に計算されるものと仮定する。 (m x) (m x x) は有意桁 p - k に丸めた正確な x として計算される。さらに正確に言うと、最下位桁 k の左の位置に基数ポイントがあることを前提に、整数に丸めることにより、x が有意桁 x を使用して丸められる。

定理 6 の証明

定理 14 により、x_h は p - k =[p/2]に丸められます。キャリーアウトがなければ、x_h は [p/2] の有意桁で表わすことができます。キャリーアウトがある場合、
x = x₀.x₁ ... x_{p - 1} × β^e であれば、x_{p - k - 1} に 1 が加算されます。キャリーアウトが存在できるのは x_{p - k - 1} =β-1 の場合だけに限られ、x_h の下位桁は 1 + x_{p - k - 1}=0 となるので、xh は [p/2] 桁で表現することができます。

xl を扱えるように、 β^p-1 ≦ x ≦ βp -1 を満たす整数に x を基準化 (スケール) します。 (ただし は x の上位桁 p - k で、 は下位桁 k) とします。

考慮すべきケースが 3 つあります。 の場合、p - k 桁に x を丸める操作はチョッピング、および x_h = 、x_l = と同じになります。 の最大桁は k なので、p が偶数の場合、 の最大桁は k =[p/2] =[p/2]となります。

これ以外の場合は、b = 2、および < 2^k-1 は、有意桁 -1 ≦ [p/2] で表わすことができます。2 つ目のケースは、 > (β/2 ) β^k-1 の場合で、x_h を計算すると切り上げが行われるので、x_h= +β^k になり、x₁ = x - x_h = x - -β^k = -β^k になります。 の最大桁は k になるので、[p/2] で表現できます。最後の = (β/2 ) β^k-1 の場合、切り上げの有無によって、x_h= 、または + β^k となります。したがって、x_l は (β/2 )β^k-1、または (β/2 ) β^k-1-β^k = - β^k/2 になり、いずれも 1 桁で表わすことができます。■

定理 6 により、2 つの精度の数を正確な和として表現することができます。和を正確に表わすための公式があります。|x|≧ |y|のとき x + y = (x y ) + (x (x y) ) y (Dekker 1971、Knuth 1981、Theorem C 第 4.2.2 節)。ただし正確な丸め操作を使用した場合、この公式は、x = 0.99998、 y = 0.99997 の例に示すとおりβ= 2 の場合に限り真になり、β= 10 には適用されません。

2 進数から 10 進数への変換

単精度では、p = 24、2²⁴ < 10⁸ になるので、2 進数から 8 桁の 10 進数に変換しても、もとの 2 進数を復元できるものと期待しがちですが、実際にはそうでありません。

定理 15

2 進の IEEE 単精度数を最も近い 8 桁の 10 進数に変換した場合、必ずしも 10 進数からもとの 2 進数を一意に復元できるわけではない。ただし、9 桁の 10 進数を使用すれば、10 進数から最も近い 2 進数への変換により、もとの浮動小数点数が復元される。

証明

2 進の単精度数は半分開いた間隔 [10³, 2¹⁰) = [1000, 1024) に存在するので、2 進小数点の右側には 10 ビット、左側には 14 ビットがそれぞれ置かれます。したがって、この間隔の中に計 (2¹⁰ - 10³)2¹⁴ = 393,216 個の 2 進数が存在することになります。10 進数を 8 桁で表わすと、同じ間隔に (2¹⁰ -10³ ) 10⁴ = 240,000 個の 10 進数が存在します。240,000 個の 10 進数で 393,216 個の 2 進数を表現することができないのは明らかです。したがって、単精度の 2 進数を一意に表現するには、8 桁の 10 進数では不十分ということになります。

9 桁あれば十分であるという点は、各 2 進数間の間隔がつねに 10 進数どうしの間隔より大きいことを明らかにすれば、実証されます。これにより、各 10 進数 N ごとに、

[N - ulp, N + ulp] の間隔に最大 1 つの 2 進数が含まれることが保証されます。したがって、各 2 進数はそれぞれ一意の 10 進数に丸められる一方、10 進数は一意の 2 進数に変換されます。

2 進数どうしの間隔が 10 進数の間隔よりつねに大きいことを証明するために、 [10ⁿ, 10^{n + 1}] という間隔の場合を考えてみます。この間隔の場合、連続した 10 進数の間隔は、 10 ^{(n + 1) - 9} になります。 [10ⁿ, 2^m] で m が最小整数の場合、 10ⁿ < 2^m、
2 進数どうしの間隔は 2^{m - 24} となり、この間隔はしだいに広くなります。したがって、10^{(n + 1) - 9} < 2^{m - 24} であることを確認すれば十分ということになります。しかし実際には、10ⁿ < 2^m になるので、10^{(n + 1) - 9} = 10ⁿ10^-8 < 2^m10^-8 < 2^m2^-24 となります。■

同じ説明を倍精度の場合にも適用して、17 桁の 10 進数を使用すれば、倍精度数を一意に復元できることを証明できます。

2 進数と 10 進数の変換は、フラグの使用に関する例にも適用されます。「精度」で、10 進数の拡張から 2 進数を復元するには、10 進数から 2 進数への変換を正確に実行しなければならないと説明しました。この変換は、拡張単精度フォーマットの N と 10^|P| (いずれも p < 13 の場合は正確な数) を乗じた後、これを単精度に丸めることによって行います(p < 0 の場合は除算で、いずれの場合も同じ操作) 。N 廖�10^|P| の計算が正確でないことは言うまでもありません。拡張単精度から単精度に変換する正確な丸め操作を組み合わせたもの (N・10^|P|) であるからです。正確さを欠く可能性のある理由を見るために、β= 10、p = 2 (単精度の場合) 、また p = 3 (拡張単精度の場合) の単純なケースを考えてみます。積が 12.51 の場合、これは拡張単精度の乗算の 1 部として 12.5 に丸められます。単精度に丸めると、12 になります。しかし、積 12.51 を単精度に丸めると 13 になるので、解は正しくありません。この誤差は倍精度による丸めによって起きたものです。

IEEE フラグを使用することにより、この倍精度の丸め誤差を避けることができます。まず、不正確なフラグの現在の値を保存して、これをリセットします。丸めモードを「ゼロへの丸め」にセットします。この後、N・10^|P| を計算します。不正確フラグの新しい値を ixflag に格納して、丸めモードと不正確フラグをリストアします。ixflag が 0 であれば、N・10^|P| は正確になるので、(N・10^|P|) を最終ビットまで正確に丸めます。また ixflag が 1 であれば、「ゼロへの丸め」では必ず切り捨てが行われるので、桁が落ちている可能性があります。積の有意桁は、1.b₁...b₂₂b₂₃...b₃₁ のようになります。倍精度の丸め誤差は、b₂₃ ...b₃₁ = 10...0 の場合に起きる可能性があります。これらのケースを検証するには、単に b₃₁ に対し ixflag の論理 OR を実行します。これにより (N・10^|P|) の丸め操作は、すべてのケースについて正しく実行されるようになります。

加法の誤差

「最適化プログラム」では、長大な和を正確に計算するときの問題を指摘しました。この精度を改善するための最も簡単なアプローチは、精度を倍加することです。精度を倍加した結果、和の精度がどの程度改善されるかを概算する場合に、s1 = x₁、s₂ = s₁   x₂..., s_i = s_i - 1 x_i とします。次に、s_i = (1 + δ_i) (s_{i - 1} + x_i) (ただし|δ_i | ≦ε) として、δi の第 2 項を無視すると、次の式が得られます。

(31)

公式 (31) の最初の等式により、 の計算値は、x_j の摂動値に正確な加法を実行した場合と同じになります。最初の項 x₁ は nεによって摂動され、最後の項 x_n はεによって摂動されます。公式 (31) の 2 つ目の等式は、誤差の範囲が によって規定されることを示すものです。精度を倍加すると、εが二乗されます。IEEE 倍精度フォーマットで和を求めると、1/ε 10¹⁶ になるので、n の妥当な値について となります。
したがって、精度を倍加すると、nεの最大摂動が適用され、 に変更されます。これにより、Kahan の加法公式に関する 2εという誤差範囲 (定理 8) は、単精度よりもはるかに改善されるものの、倍精度の場合ほどの精度は得られないことになります。

Kahan の加法公式が有効になる理由をわかりやすく説明するために、次の図に従って考えてみます。

　

加法を実行するたびに、次のループに適用される係数 C が修正されます。まず、前のループで計算した修正値 C を X_j から減算して、正しく加数 Y を求めます。次にこの加数を現在の和に加算します。Y の下位ビット (すなわち Y_l ) は和には入れられません。次に T - S を計算して、Y の上位ビットを計算します。ここから Y を減じると、Y の下位ビットが復元されます。これは、図の最初の和で失われたビットに相当します。これらのビットは、次のループの修正要素となります。Knuth 資料 (1981) の 572 ページにある定理 8 の正式な証明は、「定理 14 と定理 8」に紹介されています。

まとめ

コンピュータシステムの設計者が、浮動小数点に関する部分を無視するのは珍しいことではありません。おそらく、これはコンピュータサイエンスの授業で浮動小数点について扱う機会がほとんど、あるいはまったくないことに起因しているものと思われます。このことは、浮動小数点が定量化できる問題ではなく、この問題を扱ったハードウェアやソフトウェアについてあれこれ騒ぎ立てる意味はないといった認識を広める原因にもなっています。

本資料では、浮動小数点について積極的に論理的な検討を加えることが可能であることを実証してきました。例えば、相殺を伴う浮動小数点アルゴリズムは、使用するハードウェアで保護桁を確保し、拡張精度がサポートされた逆操作の可能な 2 進数 /10 進数変換に関する効果的なアルゴリズムを用意すれば、相対誤差を最小限に抑えることができます。使用するコンピュータシステムで浮動小数点がサポートされていれば、信頼性のある浮動小数点ソフトウェアの構築作業が大幅に簡素化されます。本資料で紹介した 2 つの例 (保護桁と拡張精度) 以外に、「システムの側面」では、命令セットの設計からコンパイラの最適化に至るまで、浮動小数点を効率的にサポートするためのさまざまな例を紹介しました。

IEEE 浮動小数点標準が広く普及していることは、標準に準拠したコードの移植性が広がることを意味します。「IEEE 標準」では、IEEE 標準の機能を利用して、実践的な浮動小数点コードをどのように記述できるかを具体的に説明しました。

謝辞

本資料は、1988 年 5 月から 7 月、Sun Microsystems 社の David Hough 氏の協力のもとに Sun Microsystems 社で開催された W. Kahan 氏の講義を参考に作成したものです。本資料の草稿を見直しに並々ならぬご協力と貴重なコメントをいただいた Kahan 氏とXerox PARC の皆様、特に John Gilbert 氏には、深く感謝の意を表わすものです。また Paul Hilfinger 氏をはじめ、関係各者にも多大なご協力をいただいたことに感謝いたします。

参考資料

Aho, Alfred V., Sethi, R., and Ullman J. D. 1986. Compilers: Principles, Techniques and Tools, Addison-Wesley, Reading, MA.

ANSI 1978. American National Standard Programming Language FORTRAN, ANSI Standard
X3.9-1978, American National Standards Institute, New York, NY.

Barnett, David 1987. A Portable Floating-Point Environment, unpublished manuscript.

Brown, W. S. 1981. A Simple but Realistic Model of Floating-Point Computation, ACM Trans. on Math. Software 7(4), pp. 445-480.

Cody, W. J et. al. 1984. A Proposed Radix- and Word-length-independent Standard for Floating-point Arithmetic, IEEE Micro 4(4), pp. 86-100.

Cody, W. J. 1988. Floating-Point Standards -- Theory and Practice, in "Reliability in Computing: the role of interval methods in scientific computing", ed. by Ramon
E. Moore, pp. 99-107, Academic Press, Boston, MA.

Coonen, Jerome 1984. Contributions to a Proposed Standard for Binary Floating-Point Arithmetic, PhD Thesis, Univ. of California, Berkeley.

Dekker, T. J. 1971. A Floating-Point Technique for Extending the Available Precision,
Numer. Math. 18(3), pp. 224-242.

Demmel, James 1984. Underflow and the Reliability of Numerical Software,
SIAM J. Sci. Stat. Comput. 5(4), pp. 887-919.

Farnum, Charles 1988. Compiler Support for Floating-point Computation, Software-Practice and Experience, 18(7), pp. 701-709.

Forsythe, G. E. and Moler, C. B. 1967. Computer Solution of Linear Algebraic Systems,
Prentice-Hall, Englewood Cliffs, NJ.

Goldberg, I. Bennett 1967. 27 Bits Are Not Enough for 8-Digit Accuracy, Comm. of the
ACM. 10(2), pp 105-106.

Goldberg, David 1990. Computer Arithmetic, in "Computer Architecture: A Quantitative Approach", by David Patterson and John L. Hennessy, Appendix A, Morgan Kaufmann, Los Altos, CA.

Golub, Gene H. and Van Loan, Charles F. 1989. Matrix Computations, 2nd edition,
The Johns Hopkins University Press, Baltimore Maryland.

Graham, Ronald L. , Knuth, Donald E. and Patashnik, Oren. 1989. Concrete Mathematics, Addison-Wesley, Reading, MA, p.162.

Hewlett Packard 1982. HP-15C Advanced Functions Handbook.

IEEE 1987. IEEE Standard 754-1985 for Binary Floating-point Arithmetic, IEEE, (1985). Reprinted in SIGPLAN 22(2) pp. 9-25.

Kahan, W. 1972. A Survey Of Error Analysis, in Information Processing 71, Vol 2,
pp. 1214 - 1239 (Ljubljana, Yugoslavia), North Holland, Amsterdam.

Kahan, W. 1986. Calculating Area and Angle of a Needle-like Triangle, unpublished manuscript.

Kahan, W. 1987. Branch Cuts for Complex Elementary Functions, in "The State of the Art in Numerical Analysis", ed. by M.J.D. Powell and A. Iserles (Univ of Birmingham, England), Chapter 7, Oxford University Press, New York.

Kahan, W. 1988. Unpublished lectures given at Sun Microsystems, Mountain View, CA.

Kahan, W. and Coonen, Jerome T. 1982. The Near Orthogonality of Syntax, Semantics, and Diagnostics in Numerical Programming Environments, in "The Relationship Between Numerical Computation And Programming Languages", ed. by J. K. Reid, pp. 103-115,
North-Holland, Amsterdam.

Kahan, W. and LeBlanc, E. 1985. Anomalies in the IBM Acrith Package, Proc. 7th IEEE Symposium on Computer Arithmetic (Urbana, Illinois), pp. 322-331.

Kernighan, Brian W. and Ritchie, Dennis M. 1978. The C Programming Language,
Prentice-Hall, Englewood Cliffs, NJ.

Kirchner, R. and Kulisch, U. 1987. Arithmetic for Vector Processors, Proc. 8th IEEE Symposium on Computer Arithmetic (Como, Italy), pp. 256-269.

Knuth, Donald E., 1981. The Art of Computer Programming, Volume II, Second Edition, Addison-Wesley, Reading, MA.

Kulisch, U. W., and Miranker, W. L. 1986. The Arithmetic of the Digital Computer: A New Approach, SIAM Review 28(1), pp 1-36.

Matula, D. W. and Kornerup, P. 1985. Finite Precision Rational Arithmetic: Slash Number
Systems, IEEE Trans. on Comput. C-34(1), pp 3-18.

Nelson, G. 1991. Systems Programming With Modula-3, Prentice-Hall, Englewood Cliffs, NJ.

Reiser, John F. and Knuth, Donald E. 1975. Evading the Drift in Floating-point Addition, Information Processing Letters 3(3), pp 84-87.

Sterbenz, Pat H. 1974. Floating-Point Computation, Prentice-Hall, Englewood Cliffs, NJ.

Swartzlander, Earl E. and Alexopoulos, Aristides G. 1975. The Sign/Logarithm Number System, IEEE Trans. Comput. C-24(12), pp. 1238-1242.

Walther, J. S., 1971. A unified algorithm for elementary functions, Proceedings of the AFIP Spring Joint Computer Conf. 38, pp. 379-385.

定理 14 と定理 8

この項では、本書で扱っていない技術的な証明を 2 つ紹介します。

定理 14

0 < k < p として、m = β^k + 1 をセットし、浮動小数点演算が正確に計算されるものと仮定する。 (m x) (m x x) は有意桁 p - k に丸めた正確な x として計算される。さらに正確に言うと、最下位桁 k の左の位置に基数ポイントがあることを前提に、整数に丸めることにより、x が有意桁 x を使用して丸められる。

証明

証明の手順は、mx = β^kx + x の計算結果におけるキャリーアウトの有無によって 2 つのケースに分けられます。

キャリーアウトがない場合は、x が整数になるように基準化 (スケール) しても問題はありません。これにより mx = x + β^k x は次のようになります。

ここで x は、2 つの部分に区分されています。下位の k 桁は b、上位の p - k 桁は a と表わします。mx から m x を計算すると、下位の k 桁 (b で表わした桁) が丸められます。

m x = mx - x mod(β^k) + rβ^k (32)

r の値は、.bb...b が より大きければ 1、これ以外の場合は 0 になります。厳密に言うと次のようになります。

a.bb...b を a + 1 に丸める場合は r = 1、これ以外の場合は r = 0 (33)

次に m x - x = mx - x mod (βk) + rβk - x = βk(x + r) - x mod (βk) を計算します。次の図に、m x - x の計算が丸められる、つまり (m x) x であることを示します。上の段は βk (x + r) を表わします (ただし B は、最下位桁 b に r を加えた場合の桁を表わします)。
.bb...b < の場合は r = 0 となり、減算により B の桁から借りが生じますが、差が丸められるため最終的には丸めの差が上段と同じ、すなわち β^kx となります。

また .bb...b > の場合は、r = 1 となり、借りのために B から 1 が減算され結果は同じくβkxとなります。

最後に .bb...b = の場合を考えます。r = 0 の場合、B は偶数、Z は奇数となり、差が切り上げられるので、結果は β^kx となります。同様に r = 1 であれば、B は奇数、Z は偶数になり、差は切り捨てられます。結果は同じくβ^kx となります。以上をまとめると次のようになります。

(m x) x =β^kx (34)

(32) と (34) の方程式を組み合わせると、 (m x) - (m x x) = x - x mod (β^k) + ρ・β^k となります。これを計算すると、次の結果が得られます。

方程式 (33) で r を計算するルールは、a...ab...b を p - k 桁に丸める場合のルールと同じです。したがって、浮動小数点演算の精度で mx - (mx - x) を計算することは、x + β^k x でキャリーアウトなしのときに x を p - k 桁に丸める場合とまったく等しくなります。

x +β^kx でキャリーアウトが伴う場合、mx = β^kx + x は次のようになります。

したがって、m x = mx - x mod (β^k) + wβ^k (ただし Z < β/2 の場合に w = -Z) となります (w の正確な値は重要ではありません) 。次に m x - x = β^kx - x mod (β^k) + wβ^k になります。図解すると次のようになります。

これを丸めると (m x) x = β^kx + wβ^k- rβ^k (ただし、.bb...b > の場合、

または .bb...b = で b₀ = 1 の場合に r = 1) となります³¹。

また (m x) - (m x x) = mx - x mod (β^k) + wβ^k - (β^kx + wβ^k - rβ^k) = x - x mod(β^k) + rβ^k となります。この場合も、a...ab...b を p - k 桁に丸めるときに切り上げを伴うと、r = 1 になります。したがって、あらゆるケースについて定理 14 が証明されたことになります。■

定理 8 (Kahan の加法公式 )

を次のアルゴリズムに従って計算する。

S = X [1]; C = 0; for j = 2 to N { Y = X [j] - C; T = S + Y; C = (T - S) - Y; S = T; }

このとき S の計算値は、S =Σx_j (1 + δ_j) + O(Nε²) Σ |x_j| (ただし |δ_j|≦ 2ε) に等しくなる。

証明

まず、単純な公式 Σx_i の誤差をどのように概算するかを考えてみます。s ₁ = x₁、s_i = (1 +δ_i) (s_i-1 +xi) を導入します。和の計算値は、s_n となり、これは各項 x をδ_j のある式で乗じた値 x_i の和に相当します。x₁ の正確な係数は (1+δ₂)(1+δ₃) ... (1+δn) なので、番号の変更により x ₂ の係数は、(1+δ₃)(1+δ₄) ... (1+δn) にならなければなりません。定理 8 は、x₁ の係数が多少複雑になる点を除き、まったく同じように証明することができます。
s₀ = c₀ = 0 で、次のように仮定すると、

ただし、上記のギリシア文字はすべてεで範囲が決まります。s_k の x₁ の係数は極端な式ですが、s_k - c_k と c_k の x₁ の係数が計算しやすくなります。k = 1 のとき、次のようになります。

以上の式で仮に、x₁ の係数をそれぞれ C_k と S_k とすると、次のようになります。

Sk と Ck の公式が得られるように、i > 1 の場合の xi を伴う項をすべて無視して、Sk と Ck の定義を拡張します。結果は次のとおりです。

Sk と Ck はε² の位までしか計算されないので、これらの公式は次のように簡略化できます。

これらの式を使用すると、次のようになります。

一般に次のように帰納法でチェックすると簡単です。

最終的に sk の x₁ の係数が必要になります。この値を得るには、x_n + 1 = 0 とし、n + 1 という添え字の付いたギリシア文字をすべて 0 として s_n + 1 を計算します。
この結果、s_n + 1 = sn - cn となり、sn の x ₁ の係数は s_n + 1 の係数より小さくなります。この係数は、s_n = 1 +η₁ - γ₁ + (4n + 2)ε² = (1 + 2ε + O(nε²)) となります。■

IEEE 754 実装間の違い

---

注 - この節は、発行された論文の一部ではありません。この節は、IEEE 規格の一部を説明し、読者が論文から推量しやすい IEEE についての思い違いを正すために追加されたものです。この資料は David Goldberg によって記述されたものではありませんが、Goldberg 氏の許可を得てここに含められています。

---

前記の論文は、プログラマはプログラムの正確度を浮動小数点演算の特性に頼ることがあるため、浮動小数点演算は十分な注意を払って実装しなければならないことを示しています。IEEE 規格は注意深い実装を規定しており、正しく動作する便利なプログラムを作成してこの規格に準拠したシステムにのみ正確な結果を配布することは可能です。そのようなプログラムはすべての IEEE システムに移植できなければならないと読者は結論付けるかもしれません。実際、197ページの「あるマシンから別のマシンにプログラムを移植するときに、両方のマシンで IEEE 標準がサポートされていれば、その間に生成される中間結果はすべてソフトウェアバグに原因があり、演算上の差によるものでない」という所見が真実であるなら、移植性のあるソフトウェアの記述は比較的簡単でしょう。

残念ながら IEEE 規格は、同一のプログラムが IEEE 準拠のすべてのシステムでまったく同じ結果を出すとは保証していません。実際、さまざまな理由から、ほとんどのプログラムはシステムによって異なった結果を出します。一例を挙げると、ほとんどのプログラムは 10 進形式と 2 進形式の間で数値を変換させる必要がありますが、IEEE 規格はこの変換を行うための正確度を完全には指定していません。また、多くのプログラムはシステムライブラリによって供給される基本関数を使用しますが、IEEE 規格はこれらの関数についてはまったく規定していません。当然、ほとんどのプログラマはこれらの機能が IEEE 規格の範囲外であることを認識しています。

IEEE 規格に規定された数値形式と演算だけを使用するプログラムですら、システムによって異なる結果を出す場合があることに気づいていないプログラマも多いことでしょう。実際、IEEE 規格の執筆者は、異なる実装が異なる結果を生むことを認めました。それは、次に示す IEEE 754 規格の用語「宛先」の定義で明らかです 。「宛先は、ユーザーによって明示的に指定されるか、システムにより暗黙に供給される (たとえば、プロシージャのサブエクスプレッションや引数における中間結果)。言語によっては、中間計算の結果をユーザー制御の及ばない宛先に置くものもあるが、当規格は演算結果を宛先の形式とオペランドの値によって定義する。」(IEEE 754-1985、7 ページ) つまり IEEE 規格は各結果をその宛先の精度に正しく丸めることを規定していますが、その宛先の精度をユーザーのプログラムで決定するとは規定していません。したがって、各システムは独自の精度で結果を宛先に配布できるため、システムがすべて IEEE 規格に準拠していても、同一のプログラムで異なる結果が出てしまいます (ときどきこの程度が著しくなります)。

前記の論文に示されたいくつかの例は、浮動小数点演算の丸めに関する知識を前提にしています。このような例を頼みにするには、プログラマはプログラムがどのように解釈されるか (具体的に IEEE システムでは各算術演算の宛先の精度）を予測する必要があります。何と、IEEE 規格の「宛先」定義は、その狭間でプログラムがどのように解釈されるかを認識するプログラマの能力を揺るがしています。その結果、前記の例のいくつかは、高級言語内で一見移植可能なプログラムとして実装される場合、通常プログラマが予期するものとは異なる精度で宛先に結果を配布する IEEE システムでは正しく動作しない場合があります。正しく動作するものもあるかもしれませんが、それらの動作が平均的なプログラマの能力を超える場合があります。

この節では、IEEE 754 演算が通常使用する宛先形式の精度に基づいた IEEE 754 演算の既存の実装について分類します。そのあとで、論文の一部の例を取り上げ、プログラムが予期する精度よりも広い精度で結果を配布すると、予期された精度が使用されるときには正しいはずでも、間違った結果が生成される場合があることを示します。また、論文内の証明の 1 つを再考して、予期しない精度に取り組む (その精度がプログラムを無効にしないにせよ) ために必要な努力を示します。これらの例は、IEEE 規格のあらゆる規定にもかかわらず、IEEE 規格が認める異なる実装間の相違のために、動作が正確に予測でき、移植可能で効率的であるという数値ソフトウェアを記述できない場合があることを示しています。このようなソフトウェアを開発するには、IEEE 規格が認めている可変性を制限するとともに、プログラムが依存する浮動小数点意味論をプログラマが表現できるプログラミング言語と環境をまず作成する必要があります。

現在の IEEE 754 実装

IEEE 754 演算の現在の実装は、ハードウェア内でそれらがサポートする各種の浮動小数点形式の程度により 2 つのグループに分けることができます。Intel x86 ファミリのプロセッサで代表される「拡張ベース」システムは拡張倍精度形式をフルサポートしますが、単精度と倍精度については部分的にしかサポートしません。拡張ベースシステムは、単精度と倍精度でデータの読み込みと格納を行い、データを拡張倍精度形式との間ですばやく変換する命令を提供します。このシステムはまた、拡張倍精度形式でレジスタに保持される場合でも算術演算の結果が単精度または倍精度に丸められる、デフォルトとは別の特殊なモードを提供します。Motorola 68000 シリーズプロセッサは、結果をこれらのモードで単精度形式または倍精度形式の精度と範囲の両方に丸めます。Intel x86 および互換プロセッサは、単精度または倍精度形式の精度に結果を丸めますが、同じ範囲を拡張倍精度形式として保持します。ほとんどの RISC プロセッサを含む単精度および倍精度システムは、単精度形式と倍精度形式をフルサポートしますが、IEEE 準拠の拡張倍精度形式はサポートしません。IBM POWER アーキテクチャは単精度については部分的なサポートしか提供していませんが、この節では単精度 / 倍精度システムとして分類します。

拡張ベースシステムにおける演算動作が単精度 / 倍精度システムにおける動作とどう異なるかを確認するために、212ページの 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 を継承します。単精度 / 倍精度システムでは、使用する上で倍精度がもっとも効率的な形式であるため、この式は倍精度で評価されます。そのため、倍精度に正しく丸められた値 3.0/7.0 が q に代入されます。次の行で、式 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 を倍精度で評価することも拡張倍精度で評価することも考えられます (1 つの x86 コンパイラにおいて、最適化によりコンパイルするときにプログラムは「Equal」と出力し、デバッグのためにコンパイルするときは「Not Equal」と出力します)。拡張ベースシステム用のコンパイラの中には、丸め精度モードを自動的に変更し、レジスタ内に結果を生成する演算によって結果を (より広い範囲を使用して) 単精度または倍精度に丸めるものもあります。そのため、このようなシステムでは、そのソースコードを読み取り IEEE 754 演算の基本的な解釈を適用するだけでは、プログラムの動作を予測できません。私たちは、ハードウェアやコンパイラが IEEE 754 準拠の環境を提供しないからといって責めることはできません。ハードウェアは、求められるとおり、正確に丸められた結果をそれぞれの宛先に配布しています。コンパイラは、許可されているとおり、ユーザーの制御の及ばない宛先に中間結果を代入しています。

拡張ベースシステムにおける演算の落とし穴

世間一般の通念では、拡張ベースシステムは、単精度 / 倍精度システムで配布される結果より正確とは言えなくても、最低限正確な結果を生成しなければならないとされています。これは、拡張ベースシステムが常に少なくても単精度 / 倍精度システムと同じ正確度を提供しており、単精度 / 倍精度システムを超える正確度を提供する場合も多いためです。前記の C プログラムのような単純な例や以下で説明する例に基づいた比較的微妙なプログラムは、この通念がせいぜい理想論であることを示しています。単精度 / 倍精度システムには確かに移植できるような、見かけ上移植性のあるプログラムの中には、拡張ベースシステムでは明らかに不正確な結果を出すものがあります。これは、コンパイラとハードウェアが、時折プログラムが期待する以上の精度を提供しようとするためです。

現在の各種のプログラミング言語では、プログラムが期待する精度を指定することが困難です。「言語とコンパイラ」の節で触れているように、多くのプログラミング言語は、同じコンテキスト内で 10.0*x のような式が発生するごとにそれらが同じ値に評価されなければならないとは指定していません。Ada など一部の言語は、この点について IEEE 規格以前のさまざまな演算におけるバリエーションの影響を受けています。さらに新しい ANSI C のような言語は、標準準拠の拡張ベースシステムに影響を受けています。実際 ANSI C 標準は、コンパイラが浮動小数点式を通常その種類に対応づけられている精度よりも広い精度に評価することを明確に許可しています。その結果、式 10.0*x の値は、次のようなさまざまな要因によって異なります。-- その式が変数にただちに代入されるのか、それともより大きな式の一部として表示されるのか、式が比較に関係するかどうか、式が引数として関数に渡されるかどうか、また渡される場合、引数は値と参照のどちらとして渡されるか、現在の精度モード、プログラムのコンパイルに使用された最適化のレベル、プログラムのコンパイル時にコンパイラが使用した精度モードと式の評価方法など。

予測のつかない式の評価は、言語規格にすべての原因があるわけではありません。拡張ベースシステムは、可能なかぎり式が拡張精度レジスタで評価される場合にもっとも効率よく動作します。しかし、格納が必要な値は、要求されるもっとも狭い精度で格納されます。10.0*x があらゆる位置で同じ値に評価されるように言語で要求させると、これらのシステムでパフォーマンス低下というペナルティが避けられなくなります。残念ながら、構文上同じコンテキストで 10.0*x を別の値に評価することをこれらのシステムに許すと、それ自体のペナルティが正確な数値ソフトウェアを開発するプログラマにかかり、意図する意味論を表現するためにプログラムの構文に頼ることができなくなります。

実際のプログラムは、指定された式は常に同じ値に評価されるという仮定に頼っているのでしょうか。ln(1 + x) を計算する定理 4 で示されたアルゴリズムを思い出してください。次に、FORTRAN による記述を示します。

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 行目で log 関数に渡されるとき、コンパイラはその値をメモリーに格納して、値を単精度に丸めることができます。そのため x のサイズが、拡張精度で 1.0 + x が1.0 に丸められるほど小さくないが、単精度では 1.0 + x が 1.0 に丸められるくらい小さいという場合、log1p(x) が返す値は x ではなくゼロになり、相対誤差は 1 になります (5e よりもいくぶん大きい)。同様に、サブエクスプレッション 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 倍にあたるおよそ 2^-23 を返します。このような状況は、少なくとも 1 つのコンパイラで発生します。前記のコードが x86 システム用の Sun WorkShop Compiler 4.2.1 FORTRAN 77 コンパイラで -O 最適化フラグを使用してコンパイルされる場合、生成されたコードは 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 度めの丸めが結合を偶数に丸める規則によって決定される場合に起きます。この 2 番目の丸めが最初の丸めと同じ方向で行われる場合、最終的な丸め誤差は最終桁のユニットの半分を超えます。しかし、二重の丸めは倍精度演算にしか影響を与えません。初めに q ビットに丸められ、続いて p ビットに丸められるような、2 つの p ビット数値の合計、差、積、または商、あるいは p ビット数値の平方根は、q ≧ 2p + 2 の場合、結果があたかも 1 度だけで p ビットに丸められたかのように同じ値を生成することがわかります。拡張倍精度は十分広いため、単精度演算は二重の丸めを受けることがありません。

正確な丸めを前提とするアルゴリズムの中には、二重の丸めに対応できないものもあります。実際、正確な丸めを必要とせず IEEE 754 に準拠しないさまざまなマシンで正しく動作するようなアルゴリズムでも二重の丸めで動作しなくなる場合があります。これらの中でもっとも便利なのは、189ページで触れているシミュレートされた多重精度演算を実行するための移植性のあるアルゴリズムです。たとえば、浮動小数点数を上位と下位の部分に分割するための定理 6 で述べられた方法は、二重の丸め演算では正しく動作しません。倍精度数 2⁵² + 3 x 2²⁶ - 1 を、それぞれ最大 26 ビットで 2 つの部分に分割してみてください。それぞれの演算が倍精度に正しく丸めれるとき、上位の部分は 2⁵² + 2²⁷ となり下位の部分は 2²⁶ - 1 となりますが、それぞれの演算が初めに拡張倍精度に丸められ続いて倍精度に丸められるとき、上位部分として 2⁵² + 2²⁸、下位部分として -2²⁶ - 1 が生成されます。-2²⁶ - 1 は 27 ビットを必要とするため、その二乗は倍精度では正確に計算できません。もちろん、この数値の二乗を拡張倍精度で計算することは可能ですが、その結果生まれるアルゴリズムは単精度 / 倍精度システムには移植できません。また、多重精度の乗算アルゴリズムによるその後のステップは、すべての部分積が倍精度で計算されていることを前提とします。倍精度変数と拡張倍精度変数の組み合わせを正確に処理するには、実装の手間が大幅に増えます。

同様に、倍精度値の配列として表現される多重精度値を加える移植性のあるアルゴリズムは、二重に丸めを行う演算では失敗することがあります。このようなアルゴリズムは、通常Kahan の加法公式に類似した手法に頼っています。239ページに挙げられた、加法公式を簡単に示した図からわかるように、s と y が |s| ≧ |y| である浮動小数点変数であり、次の計算を実行する場合、ほとんどの演算において e は t の計算で発生した丸め誤差を正確に回復させます。

t = s + y; e = (s - t) + y;

しかしこの手法は、二重に丸められる演算では無効です。s = 2⁵² + 1、y = 1/2 - 2^-54 の場合、s + y は初めに拡張倍精度で 2⁵² + 3/2 に丸められ、さらに結合を偶数に丸める規則によって倍精度で 2⁵² + 2 に丸められます。したがって t の計算における最終的な丸め誤差は 1/2 + 2^-54 です。これは、倍精度で正確に表現することは不可能であり、そのため上記の式では正確に計算できません。この場合も、拡張倍精度のまま合計を計算することにより丸め誤差を回復できますが、そうなるとプログラムは最終出力を減らして倍精度に戻す作業を行う必要があり、二重の丸めがこのプロセスも悩ます可能性があります。このような理由から、これらの方法による多重精度演算のシミュレートを行う移植性のあるプログラムは多様なマシンで正しく効率的に動作しますが、拡張ベースシステムでは公表どおりには動作しません。

一見正確な丸めに依存しているようなアルゴリズムでも、実際は二重の丸めで正しく動作するものがあります。このようなケースでは、実装ではなくアルゴリズムが公表どおりに動作するかどうか検証するために、手間のかかる二重の丸めが行われます。このことを説明するため、定理 7 の変更例を次に証明します。

定理 7'

m と n が IEEE 754 倍精度で表現可能な |m| < 2⁵² の整数で、n が n = 2ⁱ + 2^j という特殊な式で表される場合、両方の浮動小数点演算が倍精度に正しく丸められるか、あるいは初めに拡張倍精度に丸められて続いて倍精度に丸められるとき、
(m n) n = m である。

証明

損失なしで m > 0 と仮定します。q = m n とします。2 の累乗でスケールし、
2⁵² ≦ m < 2⁵³、2⁵² ≦ q < 2⁵³ となる同等の設定を考えることができます。この場合、m と q の両方とも、最下位ビットがユニットの桁を占める整数 (つまり
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 ^{+ 2)} です。これは、|ne| ≦ 1/4 を意味します。m と次に大きな表現可能な数値の差は 1 であり、m と次に小さな表現可能な数値の差は m > 252 の場合 1、m = 2⁵² の場合 1/2 であることを思い出してください。このため、|ne| ≦ 1/4 の場合、m + ne は m に丸められます (m = 2⁵² で ne = -1/4 の場合でも、積は 結合を偶数に丸める規則によって m に丸められます)。同様に、n の式が 1 + 2^-k の場合、ne は 2^-k の整数倍で、|ne| < 1/2 + 2^-(k ^{+ 1)} です。これは、|ne| ≦ 1/2 を意味します。m は q よりも確実に大きいため、この場合 m = 2⁵² にはなりません。このため、m と隣接した表現可能な数値との差は ± 1 です。したがって、|ne| ≦ 1/2 の場合も m + ne は m に丸められます (|ne| = 1/2 の場合でも、m が偶数であるため 結合を偶数に丸める規則によって積は m に丸められます)。正確に丸められる演算の証明は以上です。

二重の丸めを行う演算では (実際に 2 度丸められた場合でも) q が正確に丸められた商となることがあるため、上記と同じく |e| < 1/2 です。この場合、q n が 2 度丸められるという事実を考慮するならば、上記段落の論証を利用できます。このことを説明するため、拡張倍精度形式は少なくとも 64 個の有効ビットを持つことを IEEE 規格が規定しており、そのため 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 > 2⁵² であることに注意してください)。

最後に、二重の丸めのために q が正しく丸められた商ではない場合を考えてみましょう。この場合、最悪の場合には |e| < 1/2 + 2^-(d ^{+ 1)} となります。この場合、d は拡張倍精度形式における余分なビットの数です (既存の拡張ベースシステムはすべて、ちょうど 64 個の有効ビットによる拡張倍精度形式をサポートしています。この形式では、
d = 64 - 53 = 11)。二重の丸めは、結合を偶数に丸める規則によって 2 度目の丸めが決定される場合に限り不正確な丸め結果を出すため、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/2 + 2^{-(d + 1)} です。どちらの場合も、積の最初の丸めにより m と最大 1/2 の差がある結果が配付され、前記の証明に従って 2 度目の丸めにより m に丸められます。■

前記の証明は、商が二重の丸めを引き起こす場合にのみ積は二重の丸めを引き起こし、その場合でも正しい結果に丸めることを示しています。前記の証明はまた、推論を二重の丸めの可能性を含めるところまで拡張することは、たった 2 つの浮動小数点演算しか含まないプログラムでも困難であることも示しています。もっと複雑なプログラムでは、二重の丸めの効果を体系的に説明することは通常不可能です。倍精度演算と拡張倍精度演算の一般的な組み合わせは言うまでもありません。

拡張精度に対するプログラミング言語サポート

前記の例は、拡張精度 per se が有害であると述べているのではありません。プログラマが拡張精度を選択して使用できる場合には、多くのプログラムが拡張精度の恩恵を得ることができます。残念ながら、現在のプログラミング言語は拡張精度をいつどのような方法で使用するべきかプログラマが指定するための十分な手段を提供していません。どのようなサポートが必要かを示すため、ここでは拡張精度の使用を制御する方法を検討します。

名目上の作業精度として倍精度を使用する移植性のあるプログラムでは、より広い精度の使用を制御する方法として次の 5 つを利用できます。

1. 拡張ベースのシステムでは、可能なかぎり拡張精度を使用してもっとも速いコードを生成するようにコンパイルする。ほとんどの数値ソフトウェアでは、「マシンイプシロン」が範囲設定する各演算の相対誤差以上の演算は必要としません。メモリー内のデータが倍精度で格納される場合、その精度における最大の相対丸め誤差として通常マシンイプシロンが使用されます。そのため、入力データは入力時に丸めがなされたと (正しくまたは誤って) 想定され、結果はそれらが格納されるときに同様に丸められます。したがって、拡張精度での中間結果の計算はより正確な結果を生成する場合もありますが、拡張精度は必須ではありません。この場合、感知されるほどプログラムの速度が低下しないときだけコンパイラで拡張精度を使用し、それ以外では倍精度を使用することをお勧めします。

2. 適度に高速で十分広い場合は倍精度よりも広い形式を使用し、それ以外ではほかの形式を使用する。計算の中には、拡張精度が使用できればもっと簡単に実行できるものがあります。しかし、そのような計算もいくらか余分な手間をかけるだけで、倍精度でも実行できます。倍精度値のベクトルのユークリッド正規形を計算することを想定します。要素の二乗を計算し、広い指数範囲を使用して IEEE 754 拡張倍精度形式でそれらの合計を累積することにより、実際の長さのベクトルに対する早まったアンダーフローまたはオーバーフローを自明な方法で防ぐことができます。拡張ベースシステムでは、これが正規形を計算するもっとも速い方法です。
   単精度 / 倍精度システムでは、拡張倍精度形式はソフトウェア内でエミュレートする必要があり (拡張倍精度形式がサポートされている場合)、アンダーフローまたはオーバーフローが発生していないか確認するために例外フラグをテストし、もし発生している場合には明示的なスケールにより計算を繰り返すために、このようなエミュレーションは単純に倍精度を使用するよりも大幅に速度が低下します。拡張精度のこのような使用をサポートするために、言語は、使用するメソッドをプログラムが選択できるように、適度に高速で利用できる中でもっとも広い形式を示すとともに、その形式が十分広いこと (たとえば倍精度よりも広い範囲であること) をプログラムが検証できるようにそれぞれの形式の精度と範囲を示す環境パラメータを提供する必要があります。

3. ソフトウェアでエミュレートする必要がある場合でも、倍精度より広い形式を使用する。ユークリッド正規形の例よりも複雑なプログラムでは、プログラマは 2 種類のプログラムを記述する必要を避け、代わりに速度が遅くても拡張精度に頼ることを望む場合があります。この場合も言語は、利用できる中で最も範囲の広い形式の範囲と精度をプログラムが確認できるように、環境パラメータを提供する必要があります。

4. 倍精度より広い精度を使用しない (範囲が拡張される場合があるが、倍精度形式の精度に正しく丸められる)。上記で説明しているいくつかの例のように、正しく丸められる倍精度演算に依存するようにいたって簡単に記述されるプログラムの場合、中間結果が倍精度よりも広い指数範囲のレジスタで計算できるとしても、言語は拡張精度が使用されないように指示する方法をプログラマに提供する必要があります。この方法で計算される中間結果でも、メモリーに格納されるときにアンダーフローする場合は二重の丸めを引き起こします。算術演算の結果が初めに 53 個の有効ビットに丸められ、非正規化される必要があるときに引き続いてより少ない有効ビットに再び丸められるような場合、最終結果は非正規化数に 1 度だけ丸める場合に取得された結果とは異なる場合があります。もちろん、このような二重の丸めが実際のプログラムに不利な影響を与える可能性はほとんどありません。

5. 倍精度形式の精度と範囲の両方に正しく丸められる。倍精度のこのきびしい強制は、倍精度形式の範囲と精度の両方の限度近くで、数値ソフトウェアまたは演算そのものをテストするプログラムに最も便利です。そのような注意深いテストプログラムは、移植可能な方法で記述するのは通常困難であり、特定の形式に結果を強制的に丸めるためにダミーのサブルーチンやほかのトリックを使用しなければならない場合にはさらに難しく、エラーも生じやすくなります。したがって、拡張ベースシステムを使用してあらゆる IEEE 754 実装に移植可能な強固なソフトウェアを開発するプログラマは、並外れた労力をかけることなく単精度 / 倍精度システムの演算をエミュレートできることの真価をすぐに認めるようになるでしょう。

これらの 5 つのオプションをすべてサポートするような言語は現在ありません。実際、拡張精度の使用を制御できる機能をプログラマに提供した言語はほとんどありません。顕著な例外の 1 つは、現在標準化の最終段階にある C 言語の最新リビジョン、ISO/IEC 98 99; 1999 プログラミング言語である C 標準 です。

現在の C 標準と同様、C99 標準 は、通常関連付けられているよりも広い形式で式を評価する実装を許可しています。しかし C99 標準は、式を評価する 3 つの方法のうちどれか 1 つを使用することを推奨しています。推奨されているこの 3 つの方法の区別は、式がより広い形式に「拡張」されるその程度にあります。C99 標準は、プリプロセッサマクロ 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_t は少なくとも float と同じ広さであり、double_t は少なくとも 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 標準が提供する手段で解決できます。式 1.0 + x が任意の型の変数に代入され、その変数が常に使用される場合、log1p 関数の C99 標準バージョンは正しく動作することが保証されています。しかし、倍精度数を上位と下位の部分に分割するための移植性のある効率的な C99 標準プログラムは、比較的困難です。double の式が倍精度に正しく丸められることが保証できない場合、どうすれば正確な位置で分割し、二重の丸めを避けることができるでしょうか。1 つの解決法として、double_t 型を使用すれば、単精度 / 倍精度システムでは倍精度に分割し、拡張ベースシステムでは拡張精度に分割できます。これで、どちらの場合も演算は正しく丸められます。定理 14 は、基礎となる演算の精度がわかればどのビット位置においても分割が可能であるとしています。この情報は、FLT_EVAL_METHOD と環境パラメータマクロで取得できます。次のフラグメントは、可能な実装例を示しています。

#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 標準 はこれを実施する可搬性のある方法を提供していません。浮動小数点をサポートするために C 標準に加える変更を指定した作業文書、Floating-Point C Edit の初期の草案は、丸め精度モードを備えたシステムにおける実装で、丸め精度の取得と設定を行う fegetprec および fesetprec 関数 (丸め方向の取得と設定を行う fegetround および fesetround 関数に類似したもの) を提供することを推奨していました。この推奨は、C99 標準草案に変更が加えられる前に除かれました。

また、異なる整数演算機能を持つ各種のシステム間の移植性をサポートするC99 標準のアプローチは、異なる浮動小数点アーキテクチャをサポートするより優れた方法を提示しています。それぞれの C99 標準実装は、その実装がサポートする整数型を定義する <inttypes.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 つが当てはまり、高速な型と正確な範囲の型を区別する必要はありません。しかし拡張ベースシステムでは、難しい選択にせまられます。拡張ベースシステムは、純粋な倍精度演算、純粋な拡張精度演算のどちらも、これらの組み合わせの場合と同じ効率ではサポートせず、プログラムによって異なる組み合わせを要求します。また、いつ拡張精度を使用するかという選択は、浮動小数点演算は「本来不正確」であり、したがって整数演算にふさわしくなく、その予測力もないとベンチマークから見なす傾向がある (時には数値アナリストによってこのように公然と語られることもあります) コンパイラライターに任せてはなりません。この選択はプログラマにゆだねるべきです。プログラマは、彼らの選択を十分表現できる言語を必要とします。

結論

前述の論評は、拡張ベースシステムを非難することがねらいではなく、その意図はいくつかの誤った議論 (その最たるものは、すべての IEEE 754 システムは同一プログラムに対しまったく同じ結果を出さなければならないというもの) を明らかにすることにあります。これまで拡張ベースシステムと単精度 / 倍精度システム間の違いに焦点を当ててきましたが、それぞれのグループにおける各種のシステム間にはさらに相違があります。たとえば、単精度 / 倍精度システムの中には、2 つの数を掛けて 3 つめの数を足す単一の命令で、最後の 1 つの丸めしか行わないものがあります。「乗算/加算の混合演算」と呼ばれるこの演算では、プログラムは単精度 / 倍精度システムごとに異なる結果を生成します。この演算では、拡張精度のように、プログラムが使用されるかどうか、そしていつ使用されるかによって、同じシステム上でも同じプログラムが異なる結果を生成する場合もあります (「乗算/加算の混合演算」は、分割せずに多重精度の乗算を実行するために非可搬的な方法で使用できますが、定理 6 の分割プロセスどおりにならない場合があります)。IEEE 規格が仮にこのような演算を予想しなかったとしても、中間の積はユーザー制御が及ばない正確に保持するだけの広さを持った「宛先」に配布され、最終の合計はその単精度または倍精度の宛先のサイズに合うように正しく丸められます。

IEEE 754 は一定のプログラムが配布すべき結果を正確に規定しているという考えは、やはり魅力的です。多くのプログラマは、プログラムをコンパイルするコンパイラやプログラムを実行するコンピュータにかかわりなくプログラムの動作を理解でき、プログラムが正しく動作することを証明できると信じたがります。この確信をサポートすることは、コンピュータシステムやプログラミング言語の設計者にとって大いに価値のある目標です。残念ながら、浮動小数点演算に関しては、この目標は実際、達成が不可能です。IEEE 規格の執筆者はこれを認識しており、この達成を試みてはいません。その結果、コンピュータ業界のほぼ全体がほとんどの IEEE 754 標準に準拠しているにもかかわらず、移植性のあるソフトウェアの開発を手掛けるプログラマは、予測不可能な浮動小数点演算に取り組み続けなければなりません。

プログラマが IEEE 754 の特徴を利用するのであれば、浮動小数点演算を予想可能なものにできるプログラミング言語が必要です。C99 標準では、FLT_EVAL_METHOD ごとにプログラムを複数作成しなければなりませんが、幾分予測性が向上しています。将来の言語では IEEE 754 意味論に依存する範囲を明白に示す構文により、単一のプログラムを記述すればすむかどうかはまだ明らかではありません。既存の拡張ベースシステムは、演算が一定のシステムでどのように行われるべきかはプログラマよりもコンパイラとハードウェアの方がよく認識できると私たちに決め込ませることにより、この見通しを脅かしています。この決め込みは 2 つ目の誤信です。計算結果に求められる正確度は、結果を生成するマシンではなく、結果から引き出される結論によってのみ決まります。それらの結論が何かを把握できるのは、プログラマ、コンパイラ、ハードウェアのうちせいぜいプログラマだけです。

¹ これ以外の表現として、浮動小数点スラッシュ、および符号付き対数があります (Matula/Kornerup 1985、Swartzlander/Alexopoulos 1975 )。

² この用語は、Forsythe/Moler が 1967 年、それまでの仮数に置き換わる表現として採用したものです。

³ ここでは、指数は有意桁の上位に格納することを前提とします。

⁴ z の値が β^e_max+1 を超えない場合、あるいは β^e_min に満たない場合を除きます。この範囲を超える数については、ここでは扱いません。

⁵ z' は、z に近似する浮動小数点数とします。|d.d...d - (z/β^e) |βp^-1 は |z'-z|/ulp(z') に等しくなります 。誤差を測定するさらに正確な公式は |z'-z|/ulp(z) です。--Ed

⁶ 70 ではなく 700 です。0.1 - 0.0292 = 0.0708 なので、ulp (0.0292 ) の誤差は 708 ulp となります。- Ed

⁷ (x - y)(x + y) の式では悪性の相殺は起きないものの または の場合は x² - y² よりもわずかに精度が落ちます。この場合、x² - y² の小さい方を計算するときに起きる丸め誤差が、最終的な減算に影響を与えることはないので、 (x - y)(x + y) の丸め誤差は 3、また x² - y² の丸め誤差は 2 になります。

⁸ 「正確な丸め操作」は英文で言う "exactly rounded operation" または "correctly rounded operation" のことです。

⁹ n = 845 の場合m、x_n = 9.45、x_n + 0.555 = 10.0、10,0 - 0.555 = 9.45 となります。したがって、n > 845 については x_n = x845 となります。

¹⁰ 進数の場合、q は に等しくありません。-Ed

¹¹ 2 以外の基数についても読者自身で確認してみてください。- Ed

¹² これは、Knuth (1981、211ページ) によると Konrad Zuse に考案されたとされていますが、Goldberg (1967) によって最初に発表されたようです。

¹³ Kahan によると、拡張精度には 64 ビットの有意桁が確保されます。これは、サイクル時間を上げずに、Intel 8087 上でキャリーを伝播できる最も幅広い精度であるからです。

¹⁴ Kahan、および LeBlanc (1985 ) が提案した基本演算の 1 つとして内積を加えることには反論も見られます。

¹⁵ Kirchner は次のように指摘しています。「1 クロック・サイクルあたり 1 つの部分積ずつ、1 ulp 以内の誤差によりハードウェア上で内積を計算することができる。」
¹⁶ CORDIC とは、"Coordinate Rotation Digital Computer"の頭字語で、これにより、ほとんどの処理がシフトや加算で構成される (乗算や除算がほとんどない計算) 超越関数を計算することができます (Walther 1971)。CORDIC により、ハードウェアを追加した場合でも速度に応じて乗算器アレイに比較できるようになります。d は Intel 8087 とMotorola 68881 の両方で使用されます。
¹⁷ これは、微妙な点です。IEEE 演算のデフォルトではオーバーフローした数が に丸められますが、デフォルトを変更することもできます (「丸めモード」を参照してください)。

¹⁸ これは、854 ではサブノーマル、754 では非正規と呼ばれます。

¹⁹ これが、標準化の最も複雑な側面です。頻繁にアンダーフローするプログラムは、ソフトウェアトラップを使用するハードウェアでは極端にパフォーマンスが低下します。
²⁰ 操作に「トラップ」する NaN が伴わない限り、無効操作の例外は起きません。詳細については、IEEE Std 754-1985 の 6.2 項を参照してください。-Ed

²¹ xと y の両方が負の場合、 が より大きくなる場合があります。

²² x < 1, n < 0 で x^-n が、アンダーフローの限界 2^{e min}よりわずかに小さい場合は、 となるので、IEEE の精度はすべて -e_min < e_max となることにより、オーバーフローは起きません。

²³ これはおそらく、設計者が浮動小数点命令の精度が実際の操作からは独立した直行関係の命令セットを好む傾向を反映したものです。乗算で特別なケースを認めることにより、この直行関係がなくなります。

²⁴ ここでは、3.0 が単精度の定数、3.0D0 が倍精度の定数であることを前提とします。

²⁵ 0 ⁰ = 1 という判断は、f が非定数であることを規定する制限に依存します。この制限がなければ、f を一意に 0 とすることにより、関数は、lim_x � ₀f(x)^g(x) の値として 0 を返します。したがって、0 ⁰は NaN と定義されることになります。

²⁶ 0 ⁰については、さまざまな議論の余地がありますが、最も説得力のある説明は "Concrete Mathematics" (Graham、Knuth、Patashnik共著) で、二項定理が成り立つには 0 ⁰ = 1 でなければならないと指摘しています。
²⁷ 丸めモードが - への丸めになっている場合は除きます (この場合は x - x = -0 です)。

²⁸ VAX の VMS 数学ライブラリでは、低速な CALLS や CALLG 命令の代わりに、該当するサブルーチンをコストのかからない形式で呼び出すという意味で、一種のインライン・プロシージャ置換を使用しています。

²⁹ 前置換の難しさは、ハードウェアによる直接的な実装が必要になる点、またソフトウェア上でインプリメントする場合は、連続的な浮動小数点のトラップを設定するところにあります。- Ed

³⁰ この略式の証明では、β= 2 を前提とするので、4 による乗算は正確に行われるため、δ_i は必要ありません。

³¹ (b^kx + wb^k) で b^kx の形式が保持される場合に限り、丸めにより、b^kx + wb^k - rb^k となります。-Ed