はじめに
お久しぶりです。Sugawara Yuutaです。前の記事を書いてからもう5、6ヶ月も経っているのですね。
大変なことばかりでしたが、僕は先日無事に高校を卒業することができました。これからどうするかはぼやっとしか決まっていませんが、今からも少しずつ技術を磨いていこうと思います。
さて、今回はすでに投稿した記事とは一風変わって、文字列を浮動小数点数（英: floating-point number)へ変換するアルゴリズムの高速化に取り組んできました。
高校時代全くと言ってよいほどできなかった数学の学び直しも兼ねて行ったので、最後まで読んでいただけると嬉しいです。
なぜ重要なのか...
floatは世界中で使われています。ただこの仕組みの基本的なアイディアはシンプルだというのに、floatから文字列に変換したり、その逆は任意精度の演算が必要になり、どうしても複雑になってしまいます。なによりこの動作は時間のかかるものであり、この数十年の間で開発されたアルゴリズムも数個しかありません。
これを少しでも高速化し、シンプルにすることが今回の目標です。そんなことできるのでしょうか？と思った方もいるかもしれません。結論から言うと可能です。記事の終わりまでにその疑問を解消できれば幸いです。読んでくださっている方にはこの記事を、僕が今回開発した手法を通して、floatやそれに関するアルゴリズムを再発見する材料にしていただければ良いと思っています。
floatの仕組み
ここではfloatが一体何かをご紹介します。もしfloatが大まかにどのように動作するかを知っていたり、IEEE-754のフォーマットを知っている方は次の部分から読むことをおすすめします。
コンピューターでは複数の方法で数字を表すことができます。一番よく使われるのは、おそらく整数でしょう。どこかで聞いたことがあるでしょうが、2進数が使われ、下記のように数えます。これは4ビットの場合で、15まで表すことができ、0を含めて$2^4$パターンを表すことができます。
┌─┬─┬─┬─┐
│ │ │ │ │ 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桁目までは逆向きに${1 \over 2}$, ${1 \over 4}$, ${1 \over 8}$...というふうに数えます。
ただしちょっと待ってください。これでは表せられる最大の数が大きく下がってしまいます。ということでプログラムを書いて最大の数がいくつか見てみます。
// 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）って知っていますか？アボガドロ定数を$6.02 * 10^{23}$のように表す書き方のことです。この考え方をもとにして作られたスタンダードがIEEE-754です。
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）を使っているため同様のことは起きません。
次に真ん中は指数であり、先程出したアボガドロ定数の場合は$10^{23}$の部分に当たります。ベースは10進数では10ですが、もちろん2進数では2固定になります。小さい数も同様に表すために、指数はマイナスになることもあります。
最後に仮数部、英語ではmantissaと呼ばれます（この用法が好まれているかはおいておいて）。整数部分は正規化するときに0以外の一桁にします。ところで、2進数では0以外の1桁になりえる数が1しかないことに注目し、このビットはフォーマット上には存在しないものの、あるとして扱われます。
問題提起
ここまで二進数を見てきましたが、僕たちはもちろん10進数を使うのでこれらを行き来する方法が必要となります。具体的に言うと
- $10^x$を$2^e$に正確に変換する
必要があります。式にすると
$$ 2^e = 10^x $$
両辺から二進数の対数をとると
$$ \log_2(2^e) = \log_2(10^x) $$
$\log_{base}({base}^e) = e$より、
$$ e = \log_2(10^x) $$
にできます。ただ、これだとまだ計算が容易ではない$10^x$を評価する必要があります。
ここで$\log_{base_0}({base_1}^{e}) = e * \log_{base_0}({base_1})$を使って
$$ e = x * \log_2(10)$$
ただ、floatの指数は整数のみなので、仮数に何かをかけて調整する必要があります。式にすると
$$ m * 2^{\lfloor{e}\rfloor} = 2^e $$
これは難しくなく、並べ替えをして${{base}^{e_0} \over {base}^{e_1}} = {base}^{e_0 - e_1}$を使えば
$$ m = 2^{e - \lfloor{e}\rfloor} $$
これで問題を狭めることができます。ここで問題は
- $2^x (0 \le x < 1)$ を求める
になりました。
2^xを求める: 1. "exponentiation by squaring"
$2^x$を計算する上で一番最初に思いついたのは"exponentiation by squaring"アルゴリズムでした。これは一般的に$x$が整数であるときの指数関数の実装ですが、全て逆向きにすれば同じことが可能です。仕組み自体は難しくないので聞いたことない方は一旦記事を読む手を止めてWikipedia等読んでみてください。
ただ、"全て逆にする"せいで2乗の代わりに平方根を求めることが必要になってきます。おそらく平方根を求めるアルゴリズムでシンプルながら高速なのはニュートン法でしょう。
ニュートン法はあたえられた$f(x)$のルート(0に評価される場所)を探します。すごいところは、収束が2次収束であり、毎ループごとに正解の桁が大まかに2倍になっていくことです。実装は簡単で、以下のようにします。
$$ x_{n+1} = x_n - {f(x_n) \over f'(x_n)} $$
ちなみに、複雑にはなりますがHalley's methodを使えば3時収束も可能だそうです...
この場合、$x = \sqrt{a}$の場合、$f(x) = x^2 - a$が0になることがわかると思います。この場合の微分係数は$f'(x) = 2x$であり、これを使ってルートを計算できます。鮮明には覚えていませんが、これを試しに作った時思ったより遅くて結構落胆しました...原因はもちろん平方根で、ニュートン法を使うと複数のループが必要になること、精度を高くすると（特殊な計算方法を除いて）計算効率が$n^2$にスケールすることなどが挙げられます。必要な正確性がもっと少なくて良いならハードウェア実装も使えるのですが、そこまではやいわけではありません。
2^xを求める: 2. "テイラー展開"
テイラー展開・マクローリン展開は多項式で完全に表すことのできない関数をある点を中心に多項式で近似するという考えに基づいて存在しています。多項式で表されたほうが、たとえそれが近似であるとしても計算が楽というアドバンテージがついてくることになります。
テイラー展開は以下のように表されます。
$$ \sum_{n=1}^\infty {f^{(n)}(a) \over n!} * (x-a)^n $$
$f^{(n)}(x)$は$f(x)$のn次微分係数です。そして2^xを微分すると$2^x * \log(2)^n$となるので、
$$ 2^x \approx 1 + x\log(2) + {x^2\log(2)^2 \over 2!} + {x^3\log(2)^3 \over 3!} $$
ここで、$x^n$の計算コストを減らすために、ホーナー法を利用することができます。これを使えば
$$ p(x)=a_{n}x^{n}+a_{n-1}x^{n-1}+\cdots +a_{1}x+a_{0}$$
のn次多項式を以下のようにn回の乗算で済ませることができます。
$$ p(x)=(\cdots (a_{n}x+a_{n-1})x+\cdots +a_{1})x+a_{0} $$
さらに、テイラーの定理を使うと誤差の上限を求めることができるため、大体何項必要か知りたいときには便利かもしれません。
ちなみに、入力を意図的に減らすことによって収束を早くできます。今の入力の最大値は1($r = 1$)なのでラグランジュの剰余高を使って先ほど言ったテイラーの定理を表現すると、
$$ R(r) = {\log(2)^k * r^k \over (k+1)!} * r^{k+1} $$
となり、64bitフルで精度を要すると18項(掛け算18回)も使ってしまうということがわかります。しかし、入力を上位数ビットと残りのビットに分けて、上位ビットをループしながら上からnつ目のビットが存在したら $\sqrt[2^n]{2}$をかけるということをした後に、残りのビットでテイラー展開を行うと先にする処理も合わせて掛け算が13になるよう抑えることができます。
2^xを求める: 3. "Remezのアルゴリズム"
テイラー展開よりも掛け算の数を減らすことはできるのでしょうか？できます。ミニマックスという考え方を利用した多項式近似を生成することで達成できます。
ミニマックスとは最大の誤差の絶対値が最小になるような係数の組み合わせのことです。実際にこの分野はすでによく研究されているため、今回はremezコマンドをsollyaで利用しました。
一様ノルム（つまり最大の誤差の絶対値）がだいたい$2^{-64}$に近くなるように計算します。今回は64bitに10項、32bitに5項を使っています。
ちなみにこの項数は仮数に掛けるmを求めるときに切り捨ての代わりに四捨五入をしているので、そのまま行った時よりもさらに項数は少なくなっているはずです。(こうすると入力が $0 \le x < 1$から $0 \le x \le {1 \over 2}$になります)
ここにテイラー展開を試した時に使用したホーナー法を組み合わせると現在使っているものになります。
完成したもの
この手法は現在のGoの標準ライブラリ、strconvは使っているルックアップテーブルを使わずに、さらに標準ライブラリよりほとんどのケースでより高速に計算することができます。
今回は関数の動作も標準ライブラリに添わせたので、ParseFloatのみならインポート文を置き換えるだけで利用できるようになっていて、標準ライブラリのテストも通過しています。
完成したものはGitHubにて三条項BSDライセンスで公開しています。
https://github.com/sugawarayuuta/refloat
小噺
この後数億の様々な有効/無効なデータを入力して予期しない動作を早期発見するファジングを試したのですが、約4時間後、標準ライブラリとのミスマッチが発生し、よくみたところ標準ライブラリ側の予期しない動作のようでした。入力は800+桁の仮数があるものだったので、しょうがないかもしれないですが（標準ライブラリは大きい入力に対して800桁の10進数の配列を使って計算します...）。
おわりに
今回はfloatパーサーを紹介しました。何かわかりにくいところがあったり質問・意見、気軽にコメントなどどうぞ。今回作ったものの改善などの協力もお待ちしています。読んでくれてありがとうございました。それでは。