Lua で遊ぶ浮動小数点数

(倍精度)浮動小数点数でいろいろ遊ぶ際に、 Lua が便利なのではないかと思った。以下、特に断らない限り「浮動小数点数」と言ったら倍精度のものを指す。

なぜ Lua を使うのか?

Lua はC言語で実装されていて、数値の扱いについてはC言語と近い挙動を示すと考えられる。C言語と違って累乗の演算子 (x^y) があるのが地味に便利である。

Lua 5.2 以降では、浮動小数点数の16進表記をサポートするようになった。ソースコード中にリテラルで 0x1.fp2 と書けるし、文字列から数値に変換するときに tonumber("0x1.fp2") と書けるし、数値から文字列に変換するときは %a または %A を使って string.format("%a",7.75) と書ける。

あとは、標準ライブラリに足りない機能があったときに簡単にC言語で拡張ライブラリを書ける。

Lua の数値型について。Lua 5.2 までは、数値は全て倍精度浮動小数点数(C言語の double)を使って実装されていた。JavaScript と同じような感じである。

Lua 5.3 では、内部的に整数型(64ビット符号付き)と浮動小数点数型が分かれた。ただし、整数同士の割り算の際には自動的に浮動小数点数への変換が行われるので、どこぞの言語と違って 1/2 が 0 になるようなトラップはない。(ただし、整数同士の掛け算等では変換は行われないので、 0x100000000*0x100000000 が 0 になるというトラップがある。まあ、この記事では基本的に浮動小数点数しか扱わないので関係ないはず。)

準備

Lua の標準ライブラリは貧弱なので、 Lua の作者の一人が公開している mathx ライブラリを使う。

tarball を落としてきて、Makefile をいじってコンパイルすると、 mathx.so ができる。

OS X の場合は、 Makefile 中の MAKESO 変数の定義を -bundle -undefined dynamic_lookup を含むものに(コメントアウトを外して)書き換える。

遊ぶ

mathx.so と同じディレクトリで、ターミナルで

$ lua -lmathx

という感じにすると、 mathx ライブラリをロードした状態で対話的に実行できる。

プロンプトでは、入力の先頭に = を書くと式の値を表示できる。(Lua 5.3では、 = を書かなくても式の値を表示してくれるようになった)

16進表記

仮数部は53ビット、先頭の1を除けば52ビットである。16進表記では1桁が4ビットに対応するので、小数点以下の部分は16進13桁で表せる。

指数部の範囲は、正規化数は −1022 以上 1023 以下である。指数部が 1023 を超えると無限大となる。

> = 0x1p1023
8.9884656743116e+307
> = 0x1p1024
inf

一方、指数部が −1022 を下回ると、直ちに 0 とはならずに、非正規化数になる。

> = 0x1p-1022
2.2250738585072e-308
> = mathx.isnormal(0x1p-1022)
true
> = 0x1p-1023
1.1125369292536e-308
> = mathx.isnormal(0x1p-1023)
false

非正規化数として表せる正の数で絶対値が一番小さいのは、\(2^{-1022-52}\) である。

> = 2^(-1022-52)
4.9406564584125e-324
> = 2^(-1022-53)
0.0

マイナス0

ゼロには符号がある。0.0-0.0 を文字列化してみればそれが反映される。ただし、比較演算ではこれらは等しいものとしてみなされる。

> = 0.0
0.0
> = -0.0
-0.0
> = 0.0==-0.0
true

ちなみに、Lua 5.2 まではマイナス0を得るのに -0 と書けば良かったが、Lua 5.3では小数点なしの 0 は整数扱いになったので、小数点をつけて -0.0 と書かなければならなくなった。

ゼロの逆数をとると、無限大に符号が反映される。

> = 1/0.0
inf
> = 1/-0.0
-inf

copysign 関数を使うと、ある浮動小数点数の符号を別の浮動小数点数にコピーすることができる。当然、 Lua の標準ライブラリにはないので、 mathx ライブラリに入っているのを使う。

> = mathx.copysign(1.0, 0.0)
1.0
> = mathx.copysign(1.0, -0.0)
-1.0
> = mathx.copysign(1.0, mathx.nan)
1.0
> = mathx.copysign(1.0, -mathx.nan)
-1.0

ちなみに、C言語には signbit マクロなるものがあるようだが、 mathx ライブラリには対応する関数はないようだ。

nextafter

浮動小数点数というのは(マイナス0や無限大や NaN を除けば)所詮実数の有限な離散部分集合でしかないので、ある浮動小数点数に対して「隣の」浮動小数点数が存在するはずである。例えば、 1 の次に大きい浮動小数点数は \(1 + \frac{1}{2^{52}}\) で、1の次に小さい浮動小数点数は \(1-\frac{1}{2^{53}}\) である。

nextafter(x, y) を使うと、x の「隣の」浮動小数点数を得られる。方向は y で指定する。y には -inf か 0 か inf を指定しておけば良さそう。

> = string.format("%a", mathx.nextafter(1.0, mathx.inf))
0x1.0000000000001p+0
> = string.format("%a", mathx.nextafter(1.0, -mathx.inf))
0x1.fffffffffffffp-1
> = string.format("%a", mathx.nextafter(0.0, mathx.inf))
0x1p-1074
> = string.format("%a", mathx.nextafter(0.0, -mathx.inf))
-0x1p-1074
> = mathx.nextafter(-0.0, 0.0)
0.0
> = mathx.nextafter(-0.0, -0.0)
-0.0
> = mathx.nextafter(0.0, -0.0)
-0.0

0.0 の次に小さい浮動小数点数は -0.0 ではなく -0x1p-1074 となっている。

0x1p-n (\(=2^{-n}\)) の「次の」浮動小数点数は、通常は 0x1.0000000000001p-n となる。

> = string.format("%a", mathx.nextafter(0x1p-1019, mathx.inf))
0x1.0000000000001p-1019
> = string.format("%a", mathx.nextafter(0x1p-1020, mathx.inf))
0x1.0000000000001p-1020
> = string.format("%a", mathx.nextafter(0x1p-1021, mathx.inf))
0x1.0000000000001p-1021
> = string.format("%a", mathx.nextafter(0x1p-1022, mathx.inf))
0x1.0000000000001p-1022

しかし、 n が 1023 を超えると、様子が違ってくる。

> = string.format("%a", mathx.nextafter(0x1p-1023, mathx.inf))
0x1.0000000000002p-1023
> = string.format("%a", mathx.nextafter(0x1p-1024, mathx.inf))
0x1.0000000000004p-1024
> = string.format("%a", mathx.nextafter(0x1p-1025, mathx.inf))
0x1.0000000000008p-1025
> = string.format("%a", mathx.nextafter(0x1p-1070, mathx.inf))
0x1.1p-1070
> = string.format("%a", mathx.nextafter(0x1p-1072, mathx.inf))
0x1.4p-1072
> = string.format("%a", mathx.nextafter(0x1p-1074, mathx.inf))
0x1p-1073

非正規化数になると、仮数部の桁数が少なくなっていくのが見て取れる。


おまけ:他の言語での事情

Lua を使う理由は先に書いたが、他のスクリプト言語ではダメなのか。ダメな理由を以下に書く。

C言語: C99 では、浮動小数点数の扱いがいろいろ強化された。しかし、C言語は残念ながらスクリプト言語ではないため、対話的に遊べない。

Python: 浮動小数点数の16進表記をサポートしていない。1.0/0.0 が例外を投げる。

Ruby: 浮動小数点数の16進表記をサポートしていない。…と思いきや、%a フォーマット指定子には対応している。しかし、リテラルには対応していない。1/2 が 0 になるという罠がある。

JavaScript: 浮動小数点数の16進表記に対応していない。C言語にある諸々の関数にアクセスする手段がない。


コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です