コンパイラを浮動小数点数に対応させる（それと数学関数を実装）

Cコンパイラの実装を始めた際の興味が制御構造や関数定義・呼び出しをどう実装するかということに向いていたため int などの固定長の整数型（ポインタ含む）のみの実装にとどまっていて、 浮動小数点数は煩雑になるんじゃないかという懸念から未実装だった。

セルフホスティングやレジスタ割付などある程度の機能が実装できたので、ついに浮動小数点数の実装に手を出してみた。

ターゲットアーキテクチャ

x86で浮動小数点数を扱う方法は歴史的経緯でいくつか方法があるみたいだが、近頃のマシンではまず搭載しているだろうということでSSE2をターゲットにした

整数型の汎用レジスタとは別に、SSE2で使えるレジスタは xmm0からxmm15までの16個となる(*)。

型の表現

コンパイラ内部に型の種類として整数型というものがあり、その中で int や char, long といった型を持っている。 浮動小数点数の型をそれらとまとめて数値型の一種として扱うか、それとも整数型とは別に分けて扱うのがいいか？という方針の選択に少し悩んだ。

結論としては浮動小数点数型を新たに作成するのがいいと思う。 理由としては整数型と区別して扱うことがが大半なのと、キャストの際に別の型になっていてもコード上まったく問題ないというところだった。

値の表現

コード生成で中間表現に変換した際には式の値を仮想レジスタというもので扱っているが、それを浮動小数点数型用の仮想レジスタを別途作るかどうかでも迷った。 実際のアーキテクチャでは整数を扱う汎用レジスタと浮動小数点数のレジスタが分かれているので、同様に分けたほうがいいのかなと悩んだりした。

こっちは型を分けずに同じ仮想レジスタを使用して、浮動小数点数型かどうかのフラグを持つようにした。 理由はコード生成時に仮想レジスタが同じ型じゃないと処理する関数を別々に用意する必要があるのを避けるため。

使用する命令

生成するコードで必要となったx86の命令は限られていて、

• movsd: xmm レジスタ同士またはメモリとのやりとり
• addsd, subsd, mulsd, divsd: 四則演算
• cvtsi2sd: 整数型からの変換（汎用レジスタからxmmへ）
• cvttsd2si: 整数型への変換（xmmから汎用レジスタへ）
• ucomisd: 比較（uはUnorderedとのことだけど、付かない場合との違いがわかってない…）

だけで事が足りる。 float の場合には d の代わりに s となる。

式の取り扱いで整数型と浮動小数点数型が混在した場合などを考慮すると取り扱いが煩雑になるんじゃないかと懸念していたが、そんなことはなかった。 パースしてASTを構築した段階で式の各項の型が判明するので、 + 演算子などで左辺と右辺の方が違う場合には大きい方に合わせるよう暗黙的にキャストされて揃った状態に変換されている。

なので演算自体は整数/浮動小数点数の同じ型同士だけとなり混在した場合を考える必要がなく、生成するためのコードはシンプルで済む。

関数呼出規約(ABI)

x86-64では関数への引数をレジスタ経由で受け渡す。 浮動小数点数型の場合には整数型で使用する汎用レジスタ6個とは別に、xmm0〜xmm7 で受け渡すという規約になっている。 別々に考慮してやる必要があるので地味に面倒で、場合漏れでバグを仕込みがちだった。

戻り値は %rax の代わりに %xmm0 となる。

Callee Save、呼び出された関数側ではxmmレジスタを保存する必要がないので、破壊されたくないレジスタは呼び出し元で保存する必要がある。

関数が可変長引数の場合ははまだ対応してない。 va_arg で取り出す際に浮動小数点数型かどうかによって処理を分ける必要があり、単純なマクロでは済まないのでどうしたものか…。 （gccでは __builtin_va_XXXX と置換されて組込みの動作で処理している模様）

レジスタ割付

整数型の仮想レジスタを物理レジスタに割り付けるのと同じように、浮動小数点数用の仮想レジスタも物理レジスタに割り付ける必要がある。 が処理的にはまったく同じ操作で割付できるので、流用で解決できる。

数学関数ライブラリ

自作Cコンパイラではリンクを行わずに直接実行ファイルを生成してしまい、gccなどのホスト開発環境のライブラリを利用できないので、数学関数も自作する必要がある。 これが結構大変だし地味だしうまい方法がわからずに苦労する。

とりあえず適当に実装してみたが、本格的にするには先人の知恵を借りてnewlibだかmuslなど既存のソースを参考にするのがいいかもしれない。

当面必要になった関数をば：

floor, ceil

整数に変換して戻す。 問題点：範囲外の場合壊れる。

使用頻度高そうなので、専用命令で用意してくれてもよさそうなんだが…。

fmod

floor を利用： x - floor(x / m) * m。 問題点：除算、乗算を避けられないか？

fabs

if で符号を判定して条件分岐。 問題点：ビット操作の方がよい？

sqrt

sqrtsd 命令を使用。 問題点：精度は？

sin, cos

-π〜π の範囲でマクローリン展開。 問題点：fmod 、 forループ使用、速度、精度。

$$\begin{align} \sin x &=& x - \frac{x^3}{3!} + \frac{x^5}{5!} + \cdots &=& \sum_{n=0}^{\infty} -1^n \frac{x^{2n+1}}{(2n+1)!} \\ \cos x &=& 1 - \frac{x^2}{2!} + \frac{x^4}{4!} + \cdots &=& \sum_{n=0}^{\infty} -1^n \frac{x^{2n}}{2n!} \end{align}$$

10項程度。

tan

sin / cos。 問題点：専用に用意する。

atan

範囲を22.5度以内に絞ってマクローリン展開：

$$\begin{align} \tan^{-1} x &= x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \cdots \\ &= \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{2n+1} \end{align}$$

入力が $x$ 、求める角度を $\theta$ とすると、

$$\begin{align} \tan^{-1} x &= \theta \\ \tan \theta &= x \end{align}$$

角度を絞るためにタンジェントの加法定理：

$$\begin{align} \tan(\alpha - \beta) &= \frac{\tan \alpha - \tan \beta}{1 + \tan \alpha \tan \beta} \\ &= \frac{1 - \frac{\tan \beta}{\tan \alpha}}{\frac{1}{\tan \alpha} + \tan \beta} \end{align}$$

を利用する。

• $1 \lt x$ の場合（45度以上）
  • $\alpha=\frac{\pi}{2}$ とすると $\tan \frac{\pi}{2} = \infty$ から $\therefore \tan (\frac{\pi}{2} - \theta) = \frac{1}{x}$
• $\tan \frac{\pi}{8} \lt x \le 1$ の場合（22.5度以上45度以下）
  • $\alpha=\frac{\pi}{4}$ とすると $\tan \frac{\pi}{4} = 1$ から $\therefore \tan (\frac{\pi}{4} - \theta) = \frac{1 - x}{1 + x}$

と変換して値の範囲を狭める。

exp

$$\begin{align} e^x &= 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \frac{x^4}{4!} + \cdots \\ &= \sum_{n=0}^{\infty} \frac{x^n}{n!} \end{align}$$

で求める。

入力 $x$ を整数部 $i$ と小数部 $d$ に分けて、 $e^x = e^{i+d} = e^i e^d$ で $d$ を0.0〜1.0の範囲に絞って $e^d$ をマクローリン展開で求める。 整数部 $e^i$ は整数の指数計算のように2進数表現して求める。 負の場合は逆数にする。

ちょっと値が大きかったらすぐに発散するので、整数部がintに収まらなかったらという心配はしなくてよいかも。

log

入力$x$ を2のべき乗との積で0.5〜1.0の範囲に絞る：

$$\log x = \log (2^n \cdot d) = \log 2^n + \log d = n \log 2 + \log d$$

このままマクローリン展開：

$$\log (1 + x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \cdots$$

しても収束が遅くて打ち切り誤差がひどいので、 $x$ の符号を反転した

$$-\log (1 - x) = x + \frac{x^2}{2} + \frac{x^3}{3} + \frac{x^4}{4} + \cdots$$

との和を取ると偶数の項が消えて

$$\log (1 + x) - \log (1 - x) = \log \frac{1 + x}{1 - x} = 2 \left( x + \frac{x^3}{3} + \frac{x^5}{5} + \cdots \right) \\ $$ $y = \frac{1 + x}{1 - x}$ と置いて $x$ について解いて、 $$\begin{align} 1 + x &= y \cdot (1 - x) \\ x \cdot (y + 1) &= y - 1 \\ \therefore x &= \frac{y - 1}{y + 1} \end{align}$$

で求めた $x$ を使って求めると2乗で減っていくのでマシになる。

2じゃなくて $e$ のべき乗にすれば $\log 2$ との乗算が必要なくなるが、 IEEE-754表現が2のべき乗の形で持っていて簡単に得られるのでこちらのほうがいいかも。

pow

exp と同様に、小数部のみをマクローリン展開する。 指数関数のマクローリン展開は $\sum_{n} \frac{\log ^n a}{n !} x^n$ 。

問題点： log も使用するので処理が倍かかる。

感想

• コンピュータだから計算が得意、なのかと思っていたが数学関数は人間様の知恵で作り上げられたもので、それらを知らずにナイーブに実装すると精度も速度も実用となる数値計算はできたものではないというのが新たな発見だった。
• 現代ではハードウェアで浮動小数点数の演算をサポートしてくれているから必要はないんだけどFPUがどのように処理しているのかというのにはちょっと興味あるので、そこもソフトウェア実装してみたいと思わなくもなかったり。 てか昔のBASICやコンパイラはよく用意してくれたね…。

リンク

• リポジトリ：https://github.com/tyfkda/xcc
• SSE2 - Wikipedia
• The Newlib Homepage
• musl libc