はじめに
こんにちは．たいみょーじんです．
この記事は，自作OS Advent Calendar 2022の10日目の記事です．
対象読者:OS自作をやっていたり，興味がある人．特に30日でできる! OS自作入門(以下「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でFPUを用いた小数演算の実装を紹介します．
FPUの内部構造
Intel Software Developer's ManualのVolume 1，Chapter 8を読めばだいたい分かります．
以下のレジスタから構成されています．
R0からR7は，8個の浮動小数を保持できるスタックで，小数演算を行うスタックマシンとして機能します．
スタックは以下のようにR7が底で，プッシュすると数字の小さい方に伸びていきます．
スタック内の各レジスタの指定は絶対アドレス指定であるR0からR7の他に，スタックのトップからの相対アドレス指定ST0からST7もあります．
Control WordはFPUの設定に関するレジスタで，整数への丸めの方法(天井関数，床関数，銀行丸め，0への丸め)に関する設定(Rounding Control)や，ゼロ除算などの各種例外を発生させるかどうか(Exception Masks)などが設定できます．
Status WordはEFLAGSのFPU版のようなもので，各種例外が発生したかどうかを表すフラグや，スタックポインタなどを含みます．
FPU命令
以下FPUの制御方法を紹介しますが，基本となるFPU命令は全てIntel Software Developer's ManualのVolume2命令一覧の中に載っています．
また，全てのFPU命令の集合と，全ての'F'で始まるニーモニックの集合が等しいという，x86にしては珍しく整理された設計になっています．
hariboslinuxでは，以下のように使いたいFPU命令について，同名の関数をアセンブリで作成することで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);
(アセンブリ)
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のコンテキストを切り替えます．(ソースコード)
例外ハンドラ内で，CR0レジスタのTSフラグが立っていたら，CLTS命令でTSフラグをクリアし，FPUのコンテキストを切り替えるtake_fpu()を実行します．
また，FNSTSW命令でStatus Wordの値を取得し，例外フラグを確認します．
例外が起きていれば，FNCLEXでStatus Wordの例外フラグをクリアし，各種例外処理に進みます．
以下のソースコードではFPUでゼロ除算例外が発生したときにアプリケーションを強制終了しています．
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のコンテキストの切り替え処理です．(ソースコード)
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を更新します．
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(¤t_task->fpu_registers);
// 自分のタスクが今回初めてFPUを使用する場合，FPUを初期化する．
else init_fpu();
// FPUは俺がもらった！
fpu_user_task = current_task;
}
}
FPUの初期化はタスクごとに行われます．(ソースコード)
FNINIT命令を実行するとFPUの設定値が初期化され，FPUスタックは空になります．
次に，FNSTCW命令でControl Wordの値をメモリに書き写し，ゼロ除算例外のマスクを外してからFLDCW命令でControl Wordを書き戻します．
こうすることで，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命令と組み合わせれば，以下のように円周率を返す関数ができます．(ソースコード)
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では床関数と銀行丸めを実装しています．(ソースコード)
#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}$をプッシュします．
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$をプッシュします．
これらの命令を使用して以下のように三角関数を実装します．
(ソースコード)
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$がプッシュされます．
以下のようにして任意の底の対数関数を実装できます．(ソースコード)
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$は以下のように実装できます．
(ソースコード)
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$を求める関数を作ることができます．(ソースコード)
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が負の場合虚数などを考慮する必要があるのでとりあえずゼロ除算させる．
}
これを使って，ネイピア数を返す関数を作ることもできます．
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を用いた浮動小数演算のイメージが何となくでも分かっていただけたら嬉しいです．