四倍精度浮動小数点数のソースを表示

'''四倍精度浮動小数点形式'''（よんばいせいどふどうしょうすうてん、{{lang-en|quadruple-precision floating-point format}}）は、[[浮動小数点数]]の形式の1つで、よく使われている通常の倍精度形式と比して、仮数部の長さが約2倍である。
<!--
この[[128ビット]]四倍精度浮動小数点数は結果に倍精度以上の精度を必要とするアプリケーションのためだけに設計されたのではなく<ref>{{cite web|url=http://crd.lbl.gov/~dhbailey/dhbpapers/dhb-jmb-acat08.pdf|title=High-Precision Computation and Mathematical Physics|author=David H. Bailey and	Jonathan M. Borwein|date=July 6, 2009|accessdate=2012-11-13}}</ref>、途中計算における一時変数のオーバーフローと累積誤差の抑制により、結果が倍精度として出力される数値計算処理に対して安定性と精度を与えるために設計された。
--><!--
用途はいくらでもあるはずなので、「何々のために設計された」というのはおかしいはず。というか設計としては単に指数部を仮数部を伸ばしただけ。
-->

[[ニコラス・ハイアム|Nicholas J. Higham]]『Accuracy and Stability of Numerical Algorithms』（2nd.e）が、[[ウィリアム・カハン|カハン]]の言として引用している文によれば、
<blockquote>
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.<ref>{{cite book|first=Nicholas | last=Higham |title=Accuracy and Stability of Numerical Algorithms (2 ed)| publisher=SIAM|year=2002 | pages=43 }}</ref>

（訳）これまでのところは、10バイトの[[拡張倍精度]]は、高精度演算の価値と高速に動作させるための実装コストとの間の妥協点としては許容できるものです。そうしてすぐにさらに2バイト分長い精度が妥当になるでしょう。そうして究極的には16バイトフォーマットもが…。この種の、より広い精度への段階的な進化は[[IEEE 754]]の規格が形成された時には既に想定されていました。
</blockquote>
とのことである。（編注: 英語版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)<ref group="注釈">ケチ表現</ref>

これは十進換算で33桁〜36桁に相当する(十進小数→binary128→十進小数のラウンドトリップが可能な最大桁数が33桁、binary128→十進小数→binary128のラウンドトリップが可能な最大桁数が36桁<ref name=whyieee>{{cite web|url=http://www.cs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF|title=Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic| author=William Kahan |date=1 October 1987|accessdate=2012-11-13}}</ref>)。

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

[[Image:IEEE_754_Quadruple_Floating_Point_Format.svg]]

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

* E<sub>min</sub> = 0001<sub>16</sub>−3FFF<sub>16</sub> = −16382
* E<sub>max</sub> = 7FFE<sub>16</sub>−3FFF<sub>16</sub> = 16383
* バイアス = 3FFF<sub>16</sub> = 16383

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

指数部が0000<sub>16</sub>と7FFF<sub>16</sub>の場合は特別な場合で、それぞれ、ゼロ・非正規数、無限大・NaNの表現である。

{|class="wikitable" style="text-align:center"
! 指数 !! 仮数ゼロ !! 仮数非ゼロ !! 式
|-
| 0000<sub>16</sub> || [[0]], [[−0]] || [[非正規化数]] || <math>(-1)^{\text{signbit}} \times 2^{-16382} \times 0.\text{significandbits}_2</math>
|-
| 0001<sub>16</sub>, ..., 7FFE<sub>16</sub> ||colspan=2| 正規化数 || <math>(-1)^{\text{signbit}} \times 2^{{\text{exponentbits}_2} - 16383} \times 1.\text{significandbits}_2</math>
|-
| 7FFF<sub>16</sub> || ±[[無限|∞]] || [[NaN]] (quiet, signalling)
|}

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

=== 仮数部 ===
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 × 10<sup>4932</sup> (最大の四倍精度浮動小数点数)

 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は[[倍精度]]と同様、切り捨てられる。これは仮数のビット数が奇数であるため、丸め位置以降のビット列が<code>0101...</code>となり、1/2ulpより小さいためである。

== サポート ==

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

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

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

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

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

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

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

== double-double演算 ==
通常の倍精度浮動小数点形式のペアを用いて、四倍精度による処理を擬似的に実現する手法がある<ref name=Hida>Yozo Hida, X. Li, and D. H. Bailey, [http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.4.5769 Quad-Double Arithmetic: Algorithms, Implementation, and Application], Lawrence Berkeley National Laboratory Technical Report LBNL-46996 (2000). Also Y. Hida et al., [http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf Library for double-double and quad-double arithmetic] (2007).</ref><ref name=Shewchuk>J. R. Shewchuk, [https://www.cs.cmu.edu/~quake/robust.html Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates], Discrete & Computational Geometry 18:305-363, 1997.</ref><ref name="Knuth-4.2.3-pr9">{{cite book |last=Knuth |first=D. E. |title=The Art of Computer Programming |edition=2nd |at=chapter 4.2.3. problem 9. }} </ref>。その手法を「double-double演算」などといい、またそれによって擬似的に実現されている四倍精度が「擬似四倍精度」などと言われることがある。「double-double演算」は固定精度の計算となるため、[[任意精度演算|任意精度]]の計算と比較して高速に計算できる<ref name="Hida" />。

真の四倍精度計算とは異なり、仮数53ビットのIEEE倍精度数の1対を用いることで、double-double演算は少なくとも<ref name=Hida/>2×53=106ビットの仮数(もしくは可能性としては符号ビットの巧妙な扱いにより107ビット<ref>Robert Munafo[http://mrob.com/pub/math/f161.html F107 and F161 High-Precision Floating-Point Data Types] (2011).</ref>)を実現する計算手法である。これはIEEEのbinary128の113ビットの仮数よりも7ビット短いだけである。一方、指数部は11ビットのままであるので、表現できる値の大きさの範囲は基本的に倍精度と同じである<ref name=Hida />。これは四倍精度の指数部が15ビットであることに比べて格段に劣る(double-doubleの<math>1.8\times10^{308}</math>に対し、binary128は<math>1.2\times10^{4932}</math>)。
具体的には、double-double/四倍精度の値'''q'''をdouble-doubleの技術で表す場合、2つの倍精度数'''x'''と'''y'''の対を用いて'''q'''='''x'''+'''y'''という和の形で表現する。この各々の仮数部は'''q'''の仮数部の上下半分ずつを持つ<ref name=Shewchuk/>。つまり、'''q'''の代わりに('''x''','''y''')という対で格納されており、'''q'''に対する演算（加減乗除等）は'''x'''と'''y'''に対する等価な（ただし複雑な）演算に変換される。これにより、四倍精度の演算は（複数の）倍精度の演算の組み合わせに還元される。倍精度演算は多くの場合ハードウェアとして実装されているため、double-double演算は通常は、一般の[[任意精度演算]]よりは十分に高速である<ref name=Shewchuk/><ref name=Hida/>。
なお、double-double演算には次のような特徴がある<ref>[http://pic.dhe.ibm.com/infocenter/aix/v7r1/index.jsp?topic=%2Fcom.ibm.aix.genprogc%2Fdoc%2Fgenprogc%2F128bit_long_double_floating-point_datatype.htm 128-Bit Long Double Floating-Point Data Type]</ref>。

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

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

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

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

同様の手法で2つの四倍精度数からなるdouble-quad演算を考えることもできる。これは少なくとも226または227ビットの精度を持つ。<ref>sourceware.org [http://sourceware.org/ml/libc-alpha/2012-03/msg01024.html Re: The state of glibc libm]</ref>

== 八倍精度 ==
{{see also|[[:en:Octuple-precision floating-point format]]}}
四倍精度より精度が高く、[[任意精度演算]]より高性能な計算のためのフォーマットとして、6倍から8倍までの精度が検討されており<ref>[http://www.kek.jp/ja/NewsRoom/Release/20120123163000/ KEK 「多倍長精度計算フォーラム」第2回研究会を開催]</ref><ref>[http://suchix.kek.jp/mpcomp/20111210-forum/mpcomp-forum.pdf http://suchix.kek.jp/mpcomp/20111210-forum/mpcomp-forum.pdf 「多倍長精度演算の性能評価」(PDF)]</ref>、いくつかの方式によるソフトウェアによる実装の報告<ref>[http://www.apple.com/acg/pdf/oct3a.pdf http://www.apple.com/acg/pdf/oct3a.pdf "Octuple-precision floating point on Apple G4" (PDF)]</ref><ref name="#1">[https://cir.nii.ac.jp/crid/1572543026893254272 https://cir.nii.ac.jp/crid/1572543026893254272 「高性能な8倍精度浮動小数点演算機構の実現」]</ref><ref>[http://id.nii.ac.jp/1001/00071791/ http://id.nii.ac.jp/1001/00071791/ 情報処理学会電子図書館「GPUによる4倍・8倍精度BLASの実装と評価」]</ref><ref>[http://kaken.nii.ac.jp/d/p/21540123 http://kaken.nii.ac.jp/d/p/21540123 「Web上での簡便な8倍精度数値計算システムの構築」]</ref>、[[GRAPE]]-MP 4/6/8の試作や、[[FPGA]]への実装<ref name="#1"/>といったハードウェアの報告などがある。

== 関連項目 ==
* [[コンピュータの数値表現]]
* [[浮動小数点数]]
* {{コンピュータの数値精度}}
* [[IEEE 754]]

== 注釈 ==
{{脚注ヘルプ}}
{{Notelist}}
== 出典 ==
{{reflist|colwidth=30em}}

== 外部リンク ==
* [http://crd.lbl.gov/~dhbailey/mpdist/ High-Precision Software Directory]
* [http://sourceforge.net/p/qpfloat/home/Home/ QPFloat], a [[free software]] ([[GPL]]) software library for quadruple-precision arithmetic
* [http://www.nongnu.org/hpalib/ HPAlib], a free software ([[LGPL]]) software library for quad-precision arithmetic
* [http://gcc.gnu.org/onlinedocs/libquadmath libquadmath], the [[GNU Compiler Collection|GCC]] quad-precision math library
* [http://babbage.cs.qc.cuny.edu/IEEE-754 IEEE-754 Analysis], Interactive web page for examining Binary32, Binary64, and Binary128 floating-point values
* [https://www.hpcwire.com/2016/08/30/ibm-unveils-power9-details/ IBM Advances Against x86 with Power9 (2016年8月30日 HPCwire記事)]

{{データ型}}
{{DEFAULTSORT:よんはいせいとふとうしようすうてんすう}}
[[Category:コンピュータの算術]]
[[Category:プログラミング]]