浮動小数点演算

定義と性質

浮動小数点演算 (floating-point arithmetic) は実数をある範囲の精度で近似した数式表現で行う算術演算のこと。この近似表現を浮動小数点数という。科学技術計算では非常に大きな値や小さな値を扱う代わりに、観測精度の限界から必要な有効数字が保証されていれば十分であることが多い。

浮動小数点数は有効数字の固定された仮数 (significand) を基数 (base) の指数表記で表現される。基数は一般的に人間を対象とするなら 10、コンピュータを対象とするなら 2 または 16 である。\[ 123.45 = 1.2345 \times 10^2 \]

浮動小数点表記は指数の取り方で同一の値に対して複数の表現が可能である。例えば \(12.345\times 10^1\) は上記の値の別の表現である。このため表記上は有効数字とスケールが分かるように \(1.2345 \times 10^2\) のように先頭桁以外を小数点以下とする位置で表記するのが通例である。

IEEE 754 での浮動小数点表現

コンピュータ上での浮動小数点数は 2 進数を用いて符号 \(s\) (正負)、指数 \(e\) (整数)、仮数 \(M\) (正の整数) で構成される。この表す数値は \((-1)^s \times M \times 2^{e-E}\) である。\(E\) は指数の下駄 (bias) と呼ばれ表現形式に固有の定数である。

同じ値に対して指数の取り方によって異なる表現が発生するのを避けるため、仮数部の最も左のビットが 1 となるように正規化 (normalization) を行う。さらに、最上位ビットが 1 である前提で省略し有効桁数を 1 ビット余分に利用することができる (ケチ表現)。

123.45 を IEEE 754 binary 32 で表してみよう。IEEE 754 binary 32 は符号 1bit、指数部 8bit、仮数部 23bit で構成され指数の下駄は \(E=127\) で定義されている。まず正の値であるため符号ビットは \(s=0\) である。\[ 123.45 = (-1)^0 \times 123.45 \times 2^{127-127} \] 正規化のため仮数部が基数の 2 より小さくなるまで 2 で割り、その回数を指数部に加算してゆく。\[ \begin{align*} (-1)^0 \times 123.45 \times 2^{127-127} & = (-1)^0 \times 61.725 \times 2^{128-127} \\ & = \cdots \\ & = (-1)^0 \times 1.92890625 \times 2^{133-127} \\ & = (-1)^0 \times (1 + 0.92890625) \times 2^{133-127} \end{align*} \] 仮数部から 1 を引いた値が仮数ビットとなる。従って仮数 \(M=\)0.92890625 → 0.11101101110011001100110, 指数 \(e=\)133 → 10000101 となる。

以上より 123.45 の IEEE 754 binary 32 表現は 0 10000101 11101101110011001100110 のビットの並びとなる。

julia> bits(Float32(123.45))
"01000010111101101110011001100110"

julia> significand(123.45)
1.92890625

julia> exponent(123.45)
6

scala> import java.lang.Float
import java.lang.Float

scala> Integer.toString(Float.floatToIntBits(123.45f), 2)
res19: String = 1000010111101101110011001100110

丸め誤差

julia> a = 0.3 - 0.2
0.09999999999999998

julia> b = 0.2 - 0.1
0.1

julia> a == b
false

プログラミング

倍精度浮動小数点

倍精度浮動小数点は符号部 1bit、指数部 11bit、仮数部 52bit で構成される IEEE 754 に準拠した 64bit 浮動小数点型データ。Scala では Double、Julia では Float64、Java/C/C++ では double として使用できる。

Scala において Double の最小値、最大値、0より大きく最も小さい数は以下の定数で得られる。

scala> Double.MinValue
res112: Double = -1.7976931348623157E308

scala> Double.MaxValue
res113: Double = 1.7976931348623157E308

scala> Double.MinPositiveValue
res114: Double = 4.9E-324

仮数部で表現できるのが 53bit までの整数であることから、有効数字はそのビット数で表現が保証できる15桁 (10進数) である。またその最大値 0x1fffffffffffff を超える Long 値から Double への変換は桁落ちが発生する。

scala> val fb53 = (1L << 53) - 1
fb53: Long = 9007199254740991

scala> (fb53 + 0).toDouble
res99: Double = 9.007199254740991E15

scala> (fb53 + 1).toDouble
res100: Double = 9.007199254740992E15

scala> (fb53 + 2).toDouble
res101: Double = 9.007199254740992E15

scala> (fb53 + 1).toDouble == (fb53 + 2).toDouble
res102: Boolean = true

単精度浮動小数点

単精度浮動小数点は IEEE 754 で符号 1bit、指数 8bit、仮数 23bit で構成される 32bit 浮動小数点型データ。概ねどのような言語でもサポートしており、Scala では Float、Julia では Float32、Java, C/C++ では float として使用できる。

大量の演算に CPU の SIMD 命令や GPU などのベクトル演算ユニットを利用できる場合、倍精度よりも単精度浮動小数点を使用した方が速度やメモリ使用量の面で有利となるケースが多い。

ただし、並列演算用にベクトル化されていない実装 (演習や汎用ライブラリ目的) であれば取り立てて優位性はないため倍精度浮動小数点を使用すればよい (数年あれば状況は変わり得るためその時々の CPU アーキテクチャで判断してもらいたい)。

Scala において Float の最小値、最大値、0より大きく最も小さい数は以下の定数で得られる。

scala> (Float.MinValue, Float.MaxValue, Float.MinPositiveValue)
res119: (Float, Float, Float) = (-3.4028235E38, 3.4028235E38, 1.4E-45)

Float は仮数部で表現できるのが 24bit 整数であるため 0xffffff を超える整数からの変換は誤差を含む。

scala> val fb24 = (1L << 24) - 1
fb24: Long = 16777215

scala> (fb24 + 0).toFloat
res106: Float = 1.6777215E7

scala> (fb24 + 1).toFloat
res107: Float = 1.6777216E7

scala> (fb24 + 2).toFloat
res108: Float = 1.6777216E7

scala> (fb24 + 1).toFloat == (fb24 + 2).toFloat
res109: Boolean = true

半精度浮動小数点

半精度浮動小数点は IEEE 754 で符号 1bit、指数 5bit、仮数 10bit で構成される 16bit 浮動小数点型データ。サポートしている言語は少なく Julia での Float16 程度だが、その小ささから GPU や CPU 命令が対応している場合に単精度浮動小数点よりさらに高速である。