kubo39's blog

ただの雑記です。

浮動小数点数の10進表記についてまとめてみた

浮動小数点数

リソースをまとめたのみで、これといってアルゴリズムを解説しているわけではない。 後々時間をみつけて気になったものは記事化するかもしれない(しないかもしれない)。

はじめに

以下の前提をおく。

  • IEEE754準拠

IEEE 754-2008

IEEE754とかなにかというのは、IEEE 754(wikipedia)をみるとだいたいわかる。

符号化形式

D言語でdouble型の値から指数部・仮数部・符号を抜き出す疑似コードは以下のようになる。

// 浮動小数点数を整数表現で持つ.
ulong ival = *cast(ulong*) &val;

// 指数部は11ビットで表される.
enum log2_max_exp = cast(int) log2(1024);

// 下位52ビットは仮数部を表す.
ulong mantissa = ival & ((1L << 52) - 1);

// 真ん中11ビットは指数部:
//   52ビット右にシフトして仮数部を除いたものと,
//   (1 << 11) - 1
//      = 1 << 11 - 1
//      = 1023
//   でorをとって指数部を取り出す.
int exponent = (ival >> 52) & ((1L << 11) - 1);

// 先頭の1ビットは符号を表す.
bool negative = (ival >> (53 + 10)) & 1;

// 指数部が0のとき(0 or 非正規化数),もしくは2047のとき(無限大 or NaN)は
// 仮数部はケチ表現であり(正規化数), 先頭の1を省略している.
// そのため先頭の1を復活させる.
if (exponent != 0 && exponent != 2 * 1024 - 1)
    mantissa |= 1L << 52;
// 非正規化数の場合は指数を1にする.
else if (exponent == 0)
    exponent = 1;
// バイアスを引く. doubleの場合は1023固定.
exponent -= 1023;

double: 1.23の場合を考える。

まず単純にdoubleの整数から1ビット・11ビット・52ビットを取り出すと以下のような値となる。

s0: 0
e0: 1023
m0: 10011101011100001010001111010111000010100011110101110

s0は符号ビットで計算に関係ないので今後無視する。

次にバイアスや正規化を行う。

上記で示したようにexponentはバイアス(doubleの場合は-1023だが、exponentが0の場合は非正規化数となるので-1022になる)がかかり、 正規化数であるので仮数部も1.10011101011100001010001111010111000010100011110101110となる。

e1: 0
m1: 1.10011101011100001010001111010111000010100011110101110

m * 2 ^^ (e-bias) は 1.10011101011100001010001111010111000010100011110101110 * 2 ^^ 0 となる。

これを整数表現のみにすると

11001110101110000101000111101011100001010001111010111 * (2 ^^ -51)

のようになる。

アルゴリズム

浮動小数点数の10進表記変換を行うアルゴリズムはかなり昔から存在している。

Dragon4 (1990)

実際にはこのアルゴリズム(の元になったもの)はずいぶん以前から存在していたようだが、論文として体系化されたものは1990年の上述した論文が最初のようだ。

Grisu(Grisu1/Grisu2/Grisu3) (2010)

長らくDragon4(もしくは近しいアルゴリズム)が用いられてきた中で、2010年のGrisuの登場は当時としては画期的だったようだ。 正確な情報はわからないが、内部的にfloatを多用するJavaScriptの普及がこういった最適化研究の後押しとなったようだ。

Grisu2/Grisu3は高速だが、正しくない結果となるようなケースもあるためフォールバックとしてDragon4を用いる実装が一般的なようだ。

Errol (2016)

こちらはDragon4よりははやいがGrisu3よりは遅いらしい。

Ryu (2018)

基本的な手法はGrisu3に近い(10進リテラルへの置き換えを行う際に上限と下限から最も近似している最短の10進表示ものを探す)が、よりシンプルな整数の演算への置き換えでパフォーマンスが上がり正確な結果を得ることができるらしい。

Grisu-Exact (2020)

Grisu2 + Ryu のアプローチのようだ。作者は後述するDragonboxと同一の人。

Schubfach (2017~2018?)

Dragonboxのpaperで絶賛されている。

The beauty of Schubfach is that it does not perform any iterative search for the shortest representation. Rather, it just magically finds the "optimal number" at only one trial.

という点が画期的であるとのこと。

元々Javaのアルゴリズム(不正確な結果になってしまう問題が以前から知られていた)の改善案として提案されていたもので、2022年にようやく本家に取り込まれた。

Dragonbox (2020)

おそらく現在最も新規の手法。Schubfachの弱点にGrisuの手法を組み合わせることで改良したアルゴリズム。

  • Schubfachの弱点として以下をあげている。

there are three 128bit x 64bit -> high 64bit multiplications which are expensive. (This is for the case of double; for float, we instead have 64bit x 32bit -> high 32bit multiplications, but this is still expensive.)

  • Grisuの手法を使うことでより計算効率のよい形に変換可能

by applying some Grisu-like method, it is possible to reduce the amortized number of 128bit x 64bit -> high 64bit multiplications into one, at the cost of having some branches and divisions-by-constants.

言語・ライブラリ実装

実際にどういった言語で標準ライブラリであったり3rd patry製ライブラリであったりで用いられているかを調べてみた。

C++

Java

Go

Rust

Swift

Zig

Nim

Haskell

OCaml

Julia

Python

Hare

その他参考