四倍精度浮動小数点数

**四倍精度浮動小数点形式**（よんばいせいどふどうしょうすうてん、）は、浮動小数点数の形式の1つで、よく使われている通常の倍精度形式と比して、仮数部の長さが約2倍である。

Nicholas J. Higham『Accuracy and Stability of Numerical Algorithms』（2nd.e）が、カハンの言として引用している文によれば、

For now the 10-byte Extended format is a tolerable compromise between the value of extra-precise arithmetic and the price of implementing it to run fast; very soon two more bytes of precision will become tolerable, and ultimately a 16-byte format ... That kind of gradual evolution towards wider precision was already in view when IEEE Standard 754 for Floating-Point Arithmetic was framed.

（訳）これまでのところは、10バイトの拡張倍精度は、高精度演算の価値と高速に動作させるための実装コストとの間の妥協点としては許容できるものです。そうしてすぐにさらに2バイト分長い精度が妥当になるでしょう。そうして究極的には16バイトフォーマットもが…。この種の、より広い精度への段階的な進化はIEEE 754の規格が形成された時には既に想定されていました。

とのことである。（編注: 英語版Wikipediaが孫引きしている。英語版Wikipediaの引用元が示している引用元は「W. Kahan. Computer benchmarks versus accuracy. Draft manuscript, June 1994.」となっている）

IEEE 754-2008では、128ビットの二進フォーマットが公式に**binary128**として定義された。次に示す。

## IEEE 754 四倍精度二進浮動小数点フォーマット: binary128
IEEE 754-2008は四倍精度二進浮動小数点フォーマット**binary128**を以下のように定めている。
* 符号: 1
* 指数のビット数: 15
* 仮数のビット数: 113 (明示的には112)

これは十進換算で33桁〜36桁に相当する(十進小数→binary128→十進小数のラウンドトリップが可能な最大桁数が33桁、binary128→十進小数→binary128のラウンドトリップが可能な最大桁数が36桁)。

このフォーマットは、正規化数の場合は最上位の1が暗黙のうちに存在するものとして省略される（いわゆるケチ表現）。これにより113ビット(十進換算約34桁: \log_{10}(2^{113}) \approx 34.016)の精度の仮数部の112ビットのみがメモリフォーマット上に現れる。ビットの配置は次のようになっている。

### 指数部
他の精度の場合と同様に、指数はゲタ履き表現（IEEE 754ではexponent bias）で、オフセットは16383である。特別な場合（後述）についても基本的に他と同様である。

* Emin = 000116−3FFF16 = −16382
* Emax = 7FFE16−3FFF16 = 16383
* バイアス = 3FFF16 = 16383

ゲタ履き表現なので、実際の指数を得るには、ビットパターンによる表現からオフセットの16383を引く。

指数部が000016と7FFF16の場合は特別な場合で、それぞれ、ゼロ・非正規数、無限大・NaNの表現である。

### 値の範囲
最小の真に正の値(非正規化数)は2−16493 ≈ 10−4965であり、1ビットの精度である。
最小の正の正規化数は2−16382 ≈ 3.3621 × 10−4932であり、112ビットの精度がある。
最大の表現可能な値は216384 − 216272 ≈ 1.1897 × 104932である。

### 仮数部

TBD

### 例
ここでは十六進で表現する。符号、(バイアスされた)指数、仮数の、全てを含んでいる。

 3fff 0000 0000 0000 0000 0000 0000 0000   = 1
 c000 0000 0000 0000 0000 0000 0000 0000   = −2

 7ffe ffff ffff ffff ffff ffff ffff ffff   ≈  1.189731495357231765085759326628007 × 104932 (最大の四倍精度浮動小数点数)

 0000 0000 0000 0000 0000 0000 0000 0000   = 0
 8000 0000 0000 0000 0000 0000 0000 0000   = −0

 7fff 0000 0000 0000 0000 0000 0000 0000   = infinity
 ffff 0000 0000 0000 0000 0000 0000 0000   = −infinity

 3ffd 5555 5555 5555 5555 5555 5555 5555   ≈  1/3

デフォルトでは、1/3は倍精度と同様、切り捨てられる。これは仮数のビット数が奇数であるため、丸め位置以降のビット列が`0101...`となり、1/2ulpより小さいためである。

## サポート

四倍精度をサポートする処理系やアーキテクチャについて述べる。

利用者数が極めて少ないためか、利用環境は整っているとは言い難い現状がある。まず、処理系の独自拡張としての実装は存在するが、標準の言語仕様では四倍精度が定義されていない（独自拡張は、可搬性（portability）が必要な場合は利用できない）。マイクロプロセッサのハードウェアサポートは、携帯端末用は無論パーソナルコンピュータ用やサーバ用を含めほとんどが、倍精度までであり、ソフトウェアによる処理は一般に数倍以上遅い。

一部の処理系や、64ビットアーキテクチャの場合にC言語で `long double` が128ビット（`sizeof(long double)` が 16）のことがあるが、単に拡張倍精度フォーマットを格納するのに16バイトを使うだけであることもある。また、オプションで切り替える場合はバイナリ互換性がなくなることがありうるので注意を要する。独自拡張で（最新のGCCなど） `__float128` といった型が使えることもあるが、文字列表現との相互変換などのサポートが十分でないこともある。最近のGNU FORTRAN (gfortran) にはIEEE754の四倍精度のサポートがある。 

IBMのアーキテクチャには、次節で述べるdouble-doubleの原理による四倍精度の扱いをサポートしているものがある。また、IBM POWER9 （POWER ISA3.0）では、IEEE 754 の4倍精度FPと4倍精度十進FPをサポートした。

命令セットとしては用意があるがハードウェアでのサポートが無く、ソフトウェア実行になるというプラットフォームもある。たとえばSPARC-V9には四倍精度浮動小数点演算の命令があるが、2012年現在これをハードウェアで実装したSPARCの実機は無く、ソフトウェアで実行される。

専用計算機等での実装例の報告は少なくない。たとえばGRAPE-MPのように多倍長精度に特化した設計などもある。

（SSEでは128ビット長のレジスタを対象として浮動小数点演算を行うが、複数の単精度や倍精度演算を同時に行う機能のみで四倍精度演算のサポートはない。次節のdouble-doubleの実装に使用されることがある）

## double-double演算
通常の倍精度浮動小数点形式のペアを用いて、四倍精度による処理を擬似的に実現する手法がある。その手法を「double-double演算」などといい、またそれによって擬似的に実現されている四倍精度が「擬似四倍精度」などと言われることがある。「double-double演算」は固定精度の計算となるため、任意精度の計算と比較して高速に計算できる。

真の四倍精度計算とは異なり、仮数53ビットのIEEE倍精度数の1対を用いることで、double-double演算は少なくとも2×53=106ビットの仮数(もしくは可能性としては符号ビットの巧妙な扱いにより107ビット)を実現する計算手法である。これはIEEEのbinary128の113ビットの仮数よりも7ビット短いだけである。一方、指数部は11ビットのままであるので、表現できる値の大きさの範囲は基本的に倍精度と同じである。これは四倍精度の指数部が15ビットであることに比べて格段に劣る(double-doubleの1.8\times10^{308}に対し、binary128は1.2\times10^{4932})。
具体的には、double-double/四倍精度の値**q**をdouble-doubleの技術で表す場合、2つの倍精度数**x**と**y**の対を用いて**q**=**x**+**y**という和の形で表現する。この各々の仮数部は**q**の仮数部の上下半分ずつを持つ。つまり、**q**の代わりに(**x**,**y**)という対で格納されており、**q**に対する演算（加減乗除等）は**x**と**y**に対する等価な（ただし複雑な）演算に変換される。これにより、四倍精度の演算は（複数の）倍精度の演算の組み合わせに還元される。倍精度演算は多くの場合ハードウェアとして実装されているため、double-double演算は通常は、一般の任意精度演算よりは十分に高速である。
なお、double-double演算には次のような特徴がある。

* 値の絶対値が減少すると、追加の精度も減少する。そのため、正規化数の範囲での最小の数は倍精度よりも狭い。完全な精度を持つ最小の数は1000...02 (ゼロが106個) × 2−1074、あるいは 1.000...02 (ゼロが106個) × 2−968である。

* 実際の精度は変化する。一般に、対の下位パートの絶対値は上位パートのULPの半分を超えない。下位パートが上位パートのULPの半分より小さいならば、上位と下位の仮数の間に黙示的な全0あるいは全1のビットが存在する。仮数のビット数が固定であることに依存するアルゴリズムは128ビット長の倍精度数を使うとき失敗する可能性がある。

* 上の理由により、1 + 2−1074のような値を表現する事ができる。これは1よりも大きい最小の表現可能な値である。

更に高い精度が要求されるならば、triple-double（159または161ビット精度）やquad-double（212または215ビット精度）演算を考えることもできる。

同様の手法で2つの四倍精度数からなるdouble-quad演算を考えることもできる。これは少なくとも226または227ビットの精度を持つ。

## 八倍精度

四倍精度より精度が高く、任意精度演算より高性能な計算のためのフォーマットとして、6倍から8倍までの精度が検討されており、いくつかの方式によるソフトウェアによる実装の報告、GRAPE-MP 4/6/8の試作や、FPGAへの実装といったハードウェアの報告などがある。

## 関連項目
* コンピュータの数値表現
* 浮動小数点数
* 
* IEEE 754

## 注釈

## 出典

## 外部リンク
* High-Precision Software Directory
* QPFloat, a free software (GPL) software library for quadruple-precision arithmetic
* HPAlib, a free software (LGPL) software library for quad-precision arithmetic
* libquadmath, the GCC quad-precision math library
* IEEE-754 Analysis, Interactive web page for examining Binary32, Binary64, and Binary128 floating-point values
* IBM Advances Against x86 with Power9 (2016年8月30日 HPCwire記事)