自作OSにおけるFPUの制御

# はじめに

こんにちは．[たいみょーじん](https://twitter.com/Egh2Deywos)です．
この記事は，[自作OS Advent Calendar 2022](https://adventar.org/calendars/8392)の10日目の記事です．

対象読者:OS自作をやっていたり，興味がある人．特に[30日でできる! OS自作入門](https://book.mynavi.jp/ec/products/detail/id=22078)(以下「30日本」)を読んだことがある人．

# FPUとは？

みなさんの自作OSは，整数演算だけでなく，小数演算に対応しているでしょうか？
最近のCPUの場合SSE命令によって小数演算を行うことが多いですが，ひと昔前はCPUとは別にFPU(Floating Point Unit，浮動小数点演算装置)という小数演算専用のチップが用いられていました．
x86系列の場合，最初に8086と繋げて小数演算を行う8087というFPUが開発され，その後も87で終わる名前のFPUが続いたため，これらはx87 FPUと呼ばれます．
さらに時代が進むとx87はx86の内臓ユニットとして取り込まれ，その後SIMDの小数演算が主流になったことでFPUはその役割を終えました．
今回は，30日本で作成されるharibote OSベースの自作OS，[hariboslinux](https://github.com/TaiseiIto/hariboslinux)でFPUを用いた小数演算の実装を紹介します．

# FPUの内部構造

[Intel Software Developer's Manual](https://www.intel.com/content/www/us/en/developer/articles/technical/intel-sdm.html)のVolume 1，Chapter 8を読めばだいたい分かります．
以下のレジスタから構成されています．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2354177/fb68132a-e185-e1ea-5e49-c071d039f7f8.png)

R0からR7は，8個の浮動小数を保持できるスタックで，小数演算を行うスタックマシンとして機能します．
スタックは以下のようにR7が底で，プッシュすると数字の小さい方に伸びていきます．
スタック内の各レジスタの指定は絶対アドレス指定であるR0からR7の他に，スタックのトップからの相対アドレス指定ST0からST7もあります．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2354177/85fbb6bd-df60-c431-657d-3a1812d20848.png)

Control WordはFPUの設定に関するレジスタで，整数への丸めの方法(天井関数，床関数，銀行丸め，0への丸め)に関する設定(Rounding Control)や，ゼロ除算などの各種例外を発生させるかどうか(Exception Masks)などが設定できます．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2354177/c599e6a2-1fe4-3815-f0f9-6ba0c18a25ca.png)

Status WordはEFLAGSのFPU版のようなもので，各種例外が発生したかどうかを表すフラグや，スタックポインタなどを含みます．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2354177/4736f1bd-4640-c75f-2471-5d64bd78bbc0.png)

# FPU命令

以下FPUの制御方法を紹介しますが，基本となるFPU命令は全て[Intel Software Developer's Manual](https://www.intel.com/content/www/us/en/developer/articles/technical/intel-sdm.html)のVolume2命令一覧の中に載っています．
また，全てのFPU命令の集合と，全ての'F'で始まるニーモニックの集合が等しいという，~~x86にしては珍しく~~整理された設計になっています．
hariboslinuxでは，以下のように使いたいFPU命令について，同名の関数をアセンブリで作成することでC言語からFPU命令を呼び出せるようにしています．
([ヘッダ](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/io.h))

```c:FPU命令を呼び出す関数の宣言
// FPU instructions
void f2xm1(void);
void fcos(void);
void fincstp(void);
void fld1(void);
void fldcw(unsigned short *control);
void fldl(double *x);
void fldl2e(void);
void fldpi(void);
void fnstcw(unsigned short *control);
void fpatan(void);
void fptan(void);
void frndint(void);
void fsin(void);
void fsqrt(void);
void fstpl(double *x);
void fyl2x(void);
```

([アセンブリ](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/io.s))

```:FPU命令を呼び出す関数の実装(一部)
f2xm1:				# void f2xm1(void);
0:
	pushl	%ebp
	movl	%esp,	%ebp
	f2xm1
	leave
	ret

fcos:				# void fcos(void);
0:
	pushl	%ebp
	movl	%esp,	%ebp
	fcos
	leave
	ret

fincstp:			# void fincstp(void);
0:
	pushl	%ebp
	movl	%esp,	%ebp
	fincstp
	leave
	ret
```


# CPUとの同期

haribote OSのようなマルチタスクOSの場合，タスクスイッチによりCPUのコンテキストが切り替わりますが，それに合わせてFPUのコンテキストも入れ替える必要があります．
CPUがタスクスイッチを行うとCR0レジスタのTS(Task Switched)フラグがセットされます．
CPUは次のタスクに切り替わりましたが，FPUは直前のタスクのコンテキストを保持したままです．
別の言い方をすると，切替後のタスクから見てFPUのレジスタが他のタスクによって破壊されたように見えるのです．
この状態でFPU命令を実行しようとすると，例外`INT 0x07`が発生するので，これを拾ってFPUのコンテキストを切り替えます．([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/kernel/exception.c))
例外ハンドラ内で，CR0レジスタのTSフラグが立っていたら，`CLTS`命令でTSフラグをクリアし，FPUのコンテキストを切り替える`take_fpu()`を実行します．
また，`FNSTSW`命令でStatus Wordの値を取得し，例外フラグを確認します．
例外が起きていれば，`FNCLEX`でStatus Wordの例外フラグをクリアし，各種例外処理に進みます．
以下のソースコードではFPUでゼロ除算例外が発生したときにアプリケーションを強制終了しています．

```c:INT 0x07の例外ハンドラ
void device_not_available_exception_handler(void)
{
	unsigned short fpu_status_word;
	switch_polling_serial_mode();
	if(get_cr0() & CR0_TASK_SWITCHED)
	{               // タスクが切り替わっている場合
    	clts();     // TSフラグをクリアし，
		take_fpu(); // FPUのコンテキストを切り替える．
	}
	fpu_status_word = fnstsw(); // Status Wordレジスタの値を取得
	if(fpu_status_word & (FPU_STATUS_EXCEPTION_INVALID_OPERATION | FPU_STATUS_EXCEPTION_DENORMALIZED_OPERAND | FPU_STATUS_EXCEPTION_ZERO_DIVIDE | FPU_STATUS_EXCEPTION_OVERFLOW | FPU_STATUS_EXCEPTION_UNDERFLOW | FPU_STATUS_EXCEPTION_PRECISION))
	{ // 浮動小数演算の例外が発生したならば
		fnclex(); // 例外フラグをクリア
		if(fpu_status_word & FPU_STATUS_EXCEPTION_ZERO_DIVIDE)
		{ // ゼロ除算例外だった場合
			printf_shell(get_current_shell(), "FPU ZERO DIVIDE!!!\n");
			release_fpu();
            // アプリケーションを強制終了
			exit_application(-1, get_current_task()->task_status_segment.esp0);
		}
	}
	switch_interrupt_serial_mode();
}
```

以下がFPUのコンテキストの切り替え処理です．([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/kernel/fpu.c))
`fpu_user_task`は直前までFPUを使用していたタスクを指し示すポインタで，これが自分自身であった場合，FPUは他のタスクに破壊されていないため，FPUのコンテキストの切り替えは不要になります．
FPUが他のタスクに使用されていた場合，`FNSAVE`命令で`fpu_user_task`のために現在のFPU内容をメモリに保存してから，FPUのコンテキストを切り替えます．
切替後のタスクが過去にFPUを使用していた場合，`FRSTOR`命令でFPUのコンテキストをメモリからFPUに読み込むことで，FPUのコンテキストを復元します．
切替後のタスクが今回初めてFPU命令を実行しようとした場合，`init_fpu`関数でFPUのコンテキストの初期化を行います．
最後に，`fpu_user_task`を更新します．

```c:FPUのコンテキストの切り替え処理
void take_fpu(void)
{
	Task *current_task = get_current_task();
	if(fpu_user_task != current_task) // 最後にFPUを使用したのが自分のタスクではない場合
	{
        // 他のタスクがFPUを使用していた場合，そのタスクにおけるFPUのコンテキストをメモリに保存する．
		if(fpu_user_task)fnsave(&fpu_user_task->fpu_registers);
        // 自分のタスクが以前からFPUを使用していた場合，自分のタスクにおけるFPUのコンテキストをメモリから復元する．
		if(current_task->flags & FPU_INITIALIZED)frstor(&current_task->fpu_registers);
        // 自分のタスクが今回初めてFPUを使用する場合，FPUを初期化する．
		else init_fpu();
        // FPUは俺がもらった！
		fpu_user_task = current_task;
	}
}
```

FPUの初期化はタスクごとに行われます．([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/kernel/fpu.c))
`FNINIT`命令を実行するとFPUの設定値が初期化され，FPUスタックは空になります．
次に，`FNSTCW`命令でControl Wordの値をメモリに書き写し，ゼロ除算例外のマスクを外してから`FLDCW`命令でControl Wordを書き戻します．
こうすることで，FPUでゼロ除算が起きた時にその例外を拾うことができるようになります．

```c:FPUの初期化
void init_fpu(void)
{
	unsigned short control;
	prohibit_switch_task();
	fninit(); // FPUを初期化
	fnstcw(&control);
	control &= ~FPU_STATUS_EXCEPTION_ZERO_DIVIDE; // ゼロ除算例外だけは拾うようにする．
	fldcw(&control);
	fpu_user_task = get_current_task();
	fpu_user_task->flags |= FPU_INITIALIZED;
	allow_switch_task();
}
```

# 各種演算

さて，これでマルチタスク環境でFPUを運用する準備が整いました．
次はFPUを用いて浮動小数の各種演算を実装していきます．
しかし，C言語でfloat型，double型の演算を行う場合，gccは標準でFPU命令ではなくSSE命令を生成します．
gccにコマンドライン引数`-mno-sse -mno-sse2`を与えることでSSE命令の生成が抑制され，代わりにFPU命令が生成されます．
これで四則演算や整数(からの/への)キャストくらいはgccが勝手にコード生成してくれます．

## 定数のプッシュとポップ

FPUには以下の定数をスタックにプッシュする命令が用意されています．

|命令|プッシュされる定数|
|-|-|
|FLD1|$1$|
|FLDL2T|$\log_2 10$|
|FLDL2E|$\log_2 e$|
|FLDPI|$\pi$|
|FLDLG2|$\log_{10} 2$|
|FLDLN2|$\log_e 2$|
|FLDZ|$0$|

スタックからメモリに値をポップする`FSTPL`命令と組み合わせれば，以下のように円周率を返す関数ができます．([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/fpu.c))

```c:円周率を返す関数
double fpu_pi(void)
{
	double pi;
	fldpi();    // FPUスタックに円周率をプッシュ
	fstpl(&pi); // FPUスタックから変数piに円周率をポップ
	return pi;
}
```

## 整数への丸めの実装

x87 FPUにおける整数への丸め関数には以下の4つがあります．

* 銀行丸め(四捨五入．ただし小数部がちょうど0.5の場合最も近い偶数に丸める)
* 床関数(その数を超えない最大の整数)
* 天井関数(その数を下回らない最小の整数)
* 0への丸め(0とその数の間(境界を含む)にあるもっともその数に近い整数)

`FRNDINT`命令はスタックから値をポップし，その値を整数に丸めてスタックにプッシュしますが，どの関数で丸めるかはControl WordのRounding Controlの設定によって決まります．

hariboslinuxでは床関数と銀行丸めを実装しています．([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/fpu.c))

```c:整数への丸め関数
#define FPU_CONTROL_ROUND_TOWARD_NEAREST	0x0000
#define FPU_CONTROL_ROUND_DOWN				0x0400
#define FPU_CONTROL_ROUND_UP				0x0800
#define FPU_CONTROL_ROUND_TOWARD_ZERO		0x0c00

double fpu_floor(double x) // 床関数
{
	double result;
	unsigned short original_control;
	unsigned short changed_control;
	fnstcw(&original_control); // 現在のControl Wordの値をoriginal_controlに保存する
    // Cotrol Wordの丸め方法を床関数にする
	changed_control = (original_control & ~(FPU_CONTROL_ROUND_TOWARD_NEAREST | FPU_CONTROL_ROUND_DOWN | FPU_CONTROL_ROUND_UP | FPU_CONTROL_ROUND_TOWARD_ZERO)) | FPU_CONTROL_ROUND_DOWN;
	fldcw(&changed_control); // Control Wordに書き込む
	fldl(&x); // xをプッシュする
	frndint(); // xをポップし，整数に丸め，プッシュする．
	fstpl(&result); // 丸められたxをresultにポップする．
	fldcw(&original_control); // Control Wordを元に戻す．
	return result;
}

double fpu_trunc(double x) // 銀行丸め
{
	double result;
	unsigned short original_control;
	unsigned short changed_control;
	fnstcw(&original_control); // 現在のControl Wordの値をoriginal_controlに保存する
    // Cotrol Wordの丸め方法を銀行丸めにする
	changed_control = (original_control & ~(FPU_CONTROL_ROUND_TOWARD_NEAREST | FPU_CONTROL_ROUND_DOWN | FPU_CONTROL_ROUND_UP | FPU_CONTROL_ROUND_TOWARD_ZERO)) | FPU_CONTROL_ROUND_TOWARD_ZERO;
	fldcw(&changed_control); // Control Wordに書き込む
	fldl(&x); // xをプッシュする
	frndint(); // xをポップし，整数に丸め，プッシュする．
	fstpl(&result); // 丸められたxをresultにポップする．
	fldcw(&original_control); // Control Wordを元に戻す．
	return result;
}
```

## ルートの実装

`FSQRT`命令は，FPUスタックから$x$をポップし，$\sqrt{x}$をプッシュします．

```c:ルートを返す関数
double fpu_sqrt(double x)
{
	double result;
	fldl(&x);       // xをプッシュ
	fsqrt();        // xをポップし，sqrt(x)をプッシュ
	fstpl(&result); // sqrt(x)をresultにポップ
	return result;
}
```

## 三角関数の実装

FPUでは以下の三角関数が使えます．
角度の単位はいずれもラジアンです．

|関数|命令|
|-|-|
|正弦|FSIN|
|余弦|FCOS|
|正接|FPTAN|
|逆正接|FPATAN|

ただし，逆正接関数`FPATAN`はスタックから値$x$をポップし，さらに値$y$をポップし，$xy$平面において，原点から座標$\(x,y\)$へのベクトルの角度をプッシュします．
それに合わせて正接関数`FPTAN`はちょうど`FPATAN`の逆の操作になるようになっています．
つまり`FPTAN`命令はスタックから角度$\theta$をポップし，$\tan{\theta}$をプッシュし，さらに$1$をプッシュします．
これらの命令を使用して以下のように三角関数を実装します．
([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/fpu.c))

```c:三角関数の実装
double fpu_sin(double x) // 正弦関数
{
	double result;
	fldl(&x);       // xをプッシュ
	fsin();         // xをポップし，sin(x)をプッシュ
	fstpl(&result); // sin(x)をresultにポップ
	return result;
}

double fpu_cos(double x) // 余弦関数
{
	double result;
	fldl(&x);       // xをプッシュ
	fcos();         // xをポップし，cos(x)をプッシュ
	fstpl(&result); // cos(x)をresultにポップ
	return result;
}

double fpu_tan(double x) // 正接関数
{
	double result;
	fldl(&x);       // xをプッシュ
	fptan();        // xをポップし，tan(x)をプッシュし，1をプッシュ
	fstpl(&result); // 1をポップ
	fstpl(&result); // tan(x)をresultにポップ
	return result;
}

double fpu_asin(double x) // 逆正弦関数
{
	if(x == 1.0)return fpu_pi() / 2.0;   // asin(1) = pi / 2
	if(x == -1.0)return -fpu_pi() / 2.0; // asin(-1) = -pi / 2
	else if(0.0 <= x && x < 1.0)return fpu_atan(fpu_sqrt(x * x / (1.0 - x * x))); // 0 <= x < 1 ならば asin(x) = atan(sqrt(x^2/(1-x^2)))
	else if(-1.0 < x && x < 0.0)return -fpu_atan(fpu_sqrt(x * x / (1.0 - x * x))); // 0 <= x < 1 ならば asin(x) = -atan(sqrt(x^2/(1-x^2)))
	else return 0.0 / 0.0; // 定義域外
}

double fpu_acos(double x) // 逆余弦関数
{
	if(x == 0.0)return fpu_pi() / 2.0; // acos(0) = pi
	else if(0.0 < x && x <= 1.0)return fpu_atan(fpu_sqrt((1.0 - x * x) / (x * x))); // 0 < x <= 1 ならば acos(x) = atan(sqrt((1-x^2)/x^2))
	else if(-1.0 <= x && x < 0.0)return fpu_pi() - fpu_atan(fpu_sqrt((1.0 - x * x) / (x * x))); // -1 <= x < 0 ならば acos(x) = -atan(sqrt((1-x^2)/x^2))
	else return 0.0 / 0.0; // 定義域外
}

double fpu_atan(double x) // 逆正接関数
{
	return fpu_atan2(x, 1.0);
}
```

## 対数関数の実装

`FYL2X`命令は，FPUスタックからポップした値を$x$とし，さらにもう一度ポップした値を$y$とし，$y\log_2x$がプッシュされます．
以下のようにして任意の底の対数関数を実装できます．([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/fpu.c))

```c:対数関数
double fpu_log(double base, double power) // 底の変換の公式を使い，任意の底の対数関数を計算
{
	return fpu_log2(power) / fpu_log2(base);
}

double fpu_log2(double power)
{
	double result;
	fld1();         // 1をプッシュ
	fldl(&power);   // powerをプッシュ
	fyl2x();        // powerと1をポップし，1 * log_2 powerをプッシュ
	fstpl(&result); // log_2 powerをresultにポップ
	return result;
}
```

## 指数関数の実装

`F2XM1`命令は，FPUスタックからポップした値を$x$とし，$2^x-1$をプッシュします．
ただし，$-1 \le x \le 1$である必要があります．
出力がなぜ$2^x$ではなく$2^x-1$なのかはインテルのおじさんのみぞ知る所ですが，おそらく結果の絶対値を小さくすることで丸め誤差を小さくするためでしょう．
定義域が$-1 \le x \le 1$なので値域は$0.5 \le 2^x \le 2$となりますが，$-0.5 \le 2^x-1 \le 1$の方が絶対値が小さくなります．
浮動小数は絶対値が小さいほど高精度になるので，丸め誤差を小さくできます．
さて，この命令を使うことで，$2^x$は以下のように実装できます．
([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/fpu.c))

```c:2のx乗を求める関数
double fpu_2_to_the(double exponent)
{
	double result;
	if(-1.0 < exponent && exponent < 1.0) // 小数乗の部分だけF2XM1で計算する．
	{
		fldl(&exponent);     // 指数をプッシュする．
		f2xm1();             // 指数をポップし，2^exponent - 1をプッシュする．
		fstpl(&result);      // 2^exponent - 1をresultにポップする．
		return result + 1.0; // 2^exponentを返す．
	}
	else if(exponent == fpu_trunc(exponent)) // 指数が整数の場合
	{
		if(exponent == 1.0)return 2.0;       // 2の1乗は2
		else if(exponent == -1.0)return 0.5; // 2の-1乗は0.5
		else
		{
			int exponent_integer = (int)exponent;            // 指数を整数にキャスト
			int exponent_half = exponent_integer / 2;        // 指数を半分にする
			int exponent_remainder = exponent_integer % 2;   // 指数を半分にしたときの余り
			double power_half = fpu_2_to_the(exponent_half); // 指数を半分にして再帰呼び出し
			return power_half * power_half * fpu_2_to_the(exponent_remainder); // 全てかけ合わせる．
		}
	}
	else // 指数が整数部分と小数部分両方持っている場合
	{
		double exponent_integer = fpu_trunc(exponent);         // 指数の整数部分
		double exponent_decimal = exponent - exponent_integer; // 指数の小数部分
		return fpu_2_to_the(exponent_integer) * fpu_2_to_the(exponent_decimal); // 整数場と小数場をそれぞれ求めて掛け合わせる．
	}
}

```

これを利用して，$x^y$を求める関数を作ることができます．([ソースコード](https://github.com/TaiseiIto/hariboslinux/blob/main/src/libs/fpu.c))

```c:xのy乗を求める関数
double fpu_pow(double base, double exponent)
{
	if(base == 0.0)return 0.0; // 0のexponent乗は0(嫌な予感がするのであまり考えない)
	else if(0.0 < base) // baseが正の場合
	{
		if(-1.0 < exponent && exponent < 1.0) // 小数乗の部分は基数を2に変換してfpu_2_to_the関数に計算させる．
		{
			fldl(&exponent); // exponentをプッシュする．
			fldl(&base); // baseをプッシュする．
			fyl2x(); // baseとexponentをポップし，exponent * log_2 baseをプッシュする．
			fstpl(&exponent); // exponent * log_2 baseをexponentにポップする．exponentは基数2の場合の指数に書き換わる．
			return fpu_2_to_the(exponent); // 2のexponent乗を計算する．
		}
		else if(exponent == fpu_trunc(exponent)) // 指数が整数の場合
		{
			if(exponent == 1.0)return base; // baseの1乗はbase
			else if(exponent == -1.0)return 1.0 / base; // baseの-1乗は1/base
			else
			{
				int exponent_integer = (int)exponent;             // 指数を整数にキャスト
				int exponent_half = exponent_integer / 2;         // 指数を半分にする
				int exponent_remainder = exponent_integer % 2;    // 指数を半分にしたときの余り
				double power_half = fpu_pow(base, exponent_half); // 指数を半分にして再帰呼び出し
				return power_half * power_half * fpu_pow(base, exponent_remainder); // 整数場と小数場をそれぞれ求めて掛け合わせる．
			}
		}
		else // 指数が整数部分と小数部分両方持っている場合
		{
			double exponent_integer = fpu_trunc(exponent);         // 指数の整数部分
			double exponent_decimal = exponent - exponent_integer; // 指数の小数部分
			return fpu_pow(base, exponent_integer) * fpu_pow(base, exponent_decimal); // 整数場と小数場をそれぞれ求めて掛け合わせる．
		}
	}
	else return 0.0 / 0.0; // baseが負の場合虚数などを考慮する必要があるのでとりあえずゼロ除算させる．
}
```

これを使って，ネイピア数を返す関数を作ることもできます．

```c:ネイピア数を返す関数
double fpu_e(void)
{
	return fpu_pow(2.0, fpu_log2e()); // 2^(log_2 e)
}

double fpu_log2e(void)
{
	double result;
	fldl2e();       // log_2 eをプッシュする
	fstpl(&result); // log_2 eをresultにポップする
	return result;
}
```

以上です．FPUを用いた浮動小数演算のイメージが何となくでも分かっていただけたら嬉しいです．

# 参考文献

* [Intel Software Developer's Manual](https://www.intel.com/content/www/us/en/developer/articles/technical/intel-sdm.html)
* [Intel 8087](https://ja.wikipedia.org/wiki/Intel_8087)
* [ストリーミングSIMD拡張命令](https://ja.wikipedia.org/wiki/%E3%82%B9%E3%83%88%E3%83%AA%E3%83%BC%E3%83%9F%E3%83%B3%E3%82%B0SIMD%E6%8B%A1%E5%BC%B5%E5%91%BD%E4%BB%A4)
* [Control register](https://en.wikipedia.org/wiki/Control_register)
* [OSDev Exceptions](https://wiki.osdev.org/Exceptions)