無を犠牲にするfloatパーサーを作る
はじめに
お久しぶりです。Sugawara Yuutaです。前の記事を書いてからもう5、6ヶ月も経っているのですね。
大変なことばかりでしたが、僕は先日無事に高校を卒業することができました。これからどうするかはぼやっとしか決まっていませんが、今からも少しずつ技術を磨いていこうと思います。
さて、今回はすでに投稿した記事とは一風変わって、文字列を浮動小数点数(英: floating-point number)へ変換するアルゴリズムの高速化に取り組んできました。
高校時代全くと言ってよいほどできなかった数学の学び直しも兼ねて行ったので、最後まで読んでいただけると嬉しいです。
なぜ重要なのか...
floatは世界中で使われています。ただこの仕組みの基本的なアイディアはシンプルだというのに、floatから文字列に変換したり、その逆は任意精度の演算が必要になり、どうしても複雑になってしまいます。なによりこの動作は時間のかかるものであり、この数十年の間で開発されたアルゴリズムも数個しかありません。
これを少しでも高速化し、シンプルにすることが今回の目標です。そんなことできるのでしょうか?と思った方もいるかもしれません。結論から言うと可能です。記事の終わりまでにその疑問を解消できれば幸いです。読んでくださっている方にはこの記事を、僕が今回開発した手法を通して、floatやそれに関するアルゴリズムを再発見する材料にしていただければ良いと思っています。
floatの仕組み
コンピューターでは複数の方法で数字を表すことができます。一番よく使われるのは、おそらく整数でしょう。どこかで聞いたことがあるでしょうが、2進数が使われ、下記のように数えます。これは4ビットの場合で、15まで表すことができ、0を含めて
┌─┬─┬─┬─┐
│ │ │ │ │ 14=8+4+2
│ │ │ │ │ 5 =4+1
└─┴─┴─┴─┘ 7 =4+2+1
▲ ▲ ▲ ▲ ...
│ │ │ └─────1
│ │ └─────2
│ └─────4
└─────8
ここで問題になってくるのが少数です。日常生活でも0.25
や0.1
などの数が登場します。例えば8ビット与えられたとして、あなたなら、これをどのように実装するでしょうか?
いちばんシンプルな方法としては、真ん中に小数点をおいてしまうことでしょう。やっていることは変わりません。小数点の後、5桁目からは1, 2, 4, 8...というように桁を使い、4桁目までは逆向きに
ただしちょっと待ってください。これでは表せられる最大の数が大きく下がってしまいます。ということでプログラムを書いて最大の数がいくつか見てみます。
// https://go.dev/play/p/12BNqmE2Wi4
package main
import "fmt"
func main() {
printFrac(^uint8(0))
}
func printFrac(u8 uint8) {
fmt.Print(u8>>4, ".")
u8 &= 0xf
for u8 != 0 {
u8 *= 10
fmt.Print(u8 >> 4)
u8 &= 0xf
}
fmt.Println()
}
^
はこの場合NOT
と呼ばれ、桁を逆にする(0 -> 1, 1 -> 0)役割があります。
つまり^uint8(0)
は8ビット整数の最大値を表しているだけです。
8ビット整数の最大値は255なわけで、表示された15.9375では大きく減少していることが見て取れると思います。もっと広い範囲を表しながら、また小さい範囲も同時に表すことができるフォーマットはないんでしょうか?それがfloat、IEEE-754です。
IEEE-754
科学的表記(英: scientific notation)って知っていますか?アボガドロ定数を
Fresheneesz at the English Wikipedia project, CC 表示-継承 3.0, https://commons.wikimedia.org/w/index.php?curid=3357169による
まず、一番左がサインビット、符号でありこれが0であれば数字は正の値、1であれば負の値です。ここで気をつけないといけないのが、数字が0であっても符号には2つパターンがあることです(0, -0)。
ただ、一般的な符号付き整数は2の補数表現(英: two's complement)を使っているため同様のことは起きません。
次に真ん中は指数であり、先程出したアボガドロ定数の場合は
最後に仮数部、英語ではmantissaと呼ばれます(この用法が好まれているかはおいておいて)。整数部分は正規化するときに0以外の一桁にします。ところで、2進数では0以外の1桁になりえる数が1しかないことに注目し、このビットはフォーマット上には存在しないものの、あるとして扱われます。
問題提起
ここまで二進数を見てきましたが、僕たちはもちろん10進数を使うのでこれらを行き来する方法が必要となります。具体的に言うと
-
を10^x に正確に変換する2^e
必要があります。式にすると
両辺から二進数の対数をとると
にできます。ただ、これだとまだ計算が容易ではない
ここで
ただ、floatの指数は整数のみなので、仮数に何かをかけて調整する必要があります。式にすると
これは難しくなく、並べ替えをして
これで問題を狭めることができます。ここで問題は
-
を求める2^x (0 \le x < 1)
になりました。
2^x を求める: 1. "exponentiation by squaring"
ただ、"全て逆にする"せいで2乗の代わりに平方根を求めることが必要になってきます。おそらく平方根を求めるアルゴリズムでシンプルながら高速なのはニュートン法でしょう。
ニュートン法はあたえられた
ちなみに、複雑にはなりますがHalley's methodを使えば3時収束も可能だそうです...
この場合、
2^x を求める: 2. "テイラー展開"
テイラー展開・マクローリン展開は多項式で完全に表すことのできない関数をある点を中心に多項式で近似するという考えに基づいて存在しています。多項式で表されたほうが、たとえそれが近似であるとしても計算が楽というアドバンテージがついてくることになります。
テイラー展開は以下のように表されます。
ここで、
のn次多項式を以下のようにn回の乗算で済ませることができます。
さらに、テイラーの定理を使うと誤差の上限を求めることができるため、大体何項必要か知りたいときには便利かもしれません。
ちなみに、入力を意図的に減らすことによって収束を早くできます。今の入力の最大値は1(
となり、64bitフルで精度を要すると18項(掛け算18回)も使ってしまうということがわかります。しかし、入力を上位数ビットと残りのビットに分けて、上位ビットをループしながら上からnつ目のビットが存在したら
2^x を求める: 3. "Remezのアルゴリズム"
テイラー展開よりも掛け算の数を減らすことはできるのでしょうか?できます。ミニマックスという考え方を利用した多項式近似を生成することで達成できます。
ミニマックスとは最大の誤差の絶対値が最小になるような係数の組み合わせのことです。実際にこの分野はすでによく研究されているため、今回はremez
コマンドをsollyaで利用しました。
一様ノルム(つまり最大の誤差の絶対値)がだいたい
ちなみにこの項数は仮数に掛けるm
を求めるときに切り捨ての代わりに四捨五入をしているので、そのまま行った時よりもさらに項数は少なくなっているはずです。(こうすると入力が
ここにテイラー展開を試した時に使用したホーナー法を組み合わせると現在使っているものになります。
完成したもの
この手法は現在のGoの標準ライブラリ、strconv
は使っているルックアップテーブルを使わずに、さらに標準ライブラリよりほとんどのケースでより高速に計算することができます。
今回は関数の動作も標準ライブラリに添わせたので、ParseFloat
のみならインポート文を置き換えるだけで利用できるようになっていて、標準ライブラリのテストも通過しています。
完成したものはGitHubにて三条項BSDライセンスで公開しています。
小噺
この後数億の様々な有効/無効なデータを入力して予期しない動作を早期発見するファジングを試したのですが、約4時間後、標準ライブラリとのミスマッチが発生し、よくみたところ標準ライブラリ側の予期しない動作のようでした。入力は800+桁の仮数があるものだったので、しょうがないかもしれないですが(標準ライブラリは大きい入力に対して800桁の10進数の配列を使って計算します...)。
おわりに
今回はfloatパーサーを紹介しました。何かわかりにくいところがあったり質問・意見、気軽にコメントなどどうぞ。今回作ったものの改善などの協力もお待ちしています。読んでくれてありがとうございました。それでは。
Discussion