11. 浮動小数点演算 その１
今回は Linux に限定した話題ではありませんが，浮動小数点数の取り扱いに 関しての解説です．486DX より上位の CPU はすべて CPU に内蔵された浮動小数点 ユニット (以下 FPU) で有効数字 18 桁で 10の±4931乗の範囲の数値の演算が可能と なっています．せっかく用意されているのでアセンブラから FPU を直接使って 80 bit の浮動小数点演算をしてみます. (C の double は 64 bit, 勝った! :-)
x86系の CPU では，浮動小数点演算の実行はレジスタ中心の整数演算部分 とは違ってスタック中心のアーキテクチャを採用しています．浮動小数点 演算専用のスタックに数値を設定した後に演算命令を実行する方法を取ります． forth 言語を知っていると楽ですね ;-)
FPUのレジスタ構造
486DX より上位の CPU は浮動小数点ユニットとして以下のレジスタ群を持っています.
いくつものレジスタがありますが, 重要なレジスタは最初にある浮動小数点 データレジスタスタックです. 浮動小数点演算を行うために 80bit のレジスタを 8 本持っています. データレジスタは st0, st1,.. st7 のように 指定しますが, 整数型の命令と異なって, これらのレジスタはスタックを構成し, スタックポインタに相当する Status Word レジスタの bit 11-13 が指すレジスタ をスタックトップとして, そこから相対的に指定することになります. NASM ではスタックトップのレジスタを st0, １つ下(古い)レジスタを st1 として指定します.
Status Word, Control Word, Tag Word レジスタの各ビットの意味を簡単に示しますが, あまり気にする必要はないでしょう.
- Status Word レジスタ
B ビジーフラグ C3 条件コード3 (ZF) ST スタックトップポインタ C2 条件コード2 (PF) C1 条件コード1 C0 条件コード0 (CF) ES エラーステータス SF スタックフォルト PE 精度例外フラグ UE アンダーフロー例外フラグ OE オーバーフロー例外フラグ ZE ゼロ除算例外フラグ DE デノーマル演算例外フラグ IE 無効演算例外フラグ
- Control Word レジスタ
RC 丸めモード 00:至近値(偶数) 01:近-∞ 10:近+∞ 11:切捨て PC 仮数精度制御 00:24bit 01:予約 10:53bit 11:64bit PM 精度例外マスク UM アンダーフロー例外マスク OM オーバーフロー例外マスク ZM ゼロ除算例外マスク DM デノーマル演算例外マスク IM 無効演算例外マスク
- Tag Word レジスタ
- データレジスタの状態を保持しています. 物理レジスタ番号に対応しています.
値 意味 00 有効値 01 ゼロ 10 特殊数（無限大, デノーマル, 非数） 11 空
詳細は Intel のオンラインマニュアル(日本語)を参照して下さい．
演算とデータレジスタスタック
FPU はデータレジスタスタックに積んだ数値に対して演算を行います． 次の図の (a) (b) は 16 - 4 = 12 の計算をする場合のデータレジスタスタック の様子を示しています．(c) は 16 と 4 を比較する場合の例です.
実行しているのは次のコードです．fld でメモリにある80ビットの浮動小数を スタックにつみます. (a) の fsub 命令でスタックの上部の2つの数値を取り出して 減算し, 結果をスタックに戻します. (b) の fsubp 命令では fsub の後にスタック をポップ(スタックトップを空にした後,ポインタをインクリメント)します. (c) の例では, fcompp で大きさを比較した結果を C0, C2, C3 の条件コードに設定 してスタックを2回ポップしています. 比較に使った数値が不用の場合です.
Num16 dt 16.0
Num4 dt 4.0
fld tword[Num16]
fld tword[Num4]
fsub st1, st0 ; (a)
; または fsubp st1,st0 (b)
; または fcompp st1,st0 (c)
使いやすいのは (b) の方法です. 直前の演算結果を次の処理で利用できるためです. FPU で「スタックのポップ」は「数値の廃棄」に相当することに注意して下さい.
fsub 命令の代わりに fadd なら減算, fmul なら乗算, fdiv なら除算ができます. fsin のような三角関数, 平方根を求める fsqrt などの関数はスタックトップに設定した引数が関数値に置き換わります.
浮動小数の入出力
浮動小数点演算自体は簡単ですが，難しい部分は浮動小数の入出力の部分です．
浮動小数点専用の命令を持っているため計算自体は非常に簡単です. しかし， 文字として入力された浮動小数点数を演算に使用できる形式に変換したり， 演算結果を (文字として) 表示するのは，計算に比べて非常に面倒です. 演算は内部的に 80 ビットの2進数値で行われるため，数値は80 ビットの ２進数値で与え，結果も 80 ビットの ２進数値で返ります.
class="verbatim-code"実数を表す文字列と 80 ビットの浮動小数点数の２進表現との相互変換をする ためのサブルーチンを用意します.
FPU の 80bit の 2進浮動小数点数で表現された内部データを表示のために 10進数に変換する部分と文字列で表現された浮動小数点数を FPU の 80bit ２進浮動小数点数に変換するためのプログラムです．
;---------------------------------------------------------------------
; Floating Point
; file : float.inc
; 2000/11/04 (C) Jun Mizutani
;---------------------------------------------------------------------
%ifndef __FLOAT_INC
%define __FLOAT_INC
section .bss
Sign resb 1 ; 符号
BCD rest 1 ; 仮数値 BCD 表現(10バイト)
ExpSign resb 1 ; 指数部符号
Exponent resd 1 ; 指数値(4バイト)
DigitStr resb 20 ; BCD変換用数字列 buffer
section .text
;---------------------------------------------------------------------
; FPU のスタックトップのレジスタ中の浮動小数点数をバッファに文字列
; として格納する. データレジスタスタックは pop されるため, 必要なら
; call 前に fld st0 で dup しておくこと.
; BCDはリトルエンディアンで格納されていることに注意
; レジスタは保存される
; edi に文字列用のバッファ先頭アドレスを指定して呼ぶ (29byte以上)
;---------------------------------------------------------------------
Float2Ascii: pusha
call Float2BCD ; 浮動小数点数
ffree st0 ; スタックトップを空に
fincstp ; pop
cmp byte[Sign], 0 ; 符号をチェック
je .positive
mov al, "-"
mov [edi], al
inc edi
.positive
mov esi, BCD
add esi, 8 ; BCD先頭無効バイトをスキップ
mov bh, [esi]
mov bl, bh
shr bh, 4
or bh, 0x30
mov [edi], bh ; 小数点以下第1桁
inc edi
mov al, "." ; 小数点
mov [edi], al
inc edi
and bl, 0x0F
or bl, 0x30
mov [edi], bl ; 第1桁
inc edi
.next
dec esi ; 次のBCD1バイトを取得
mov bh, [esi]
mov bl, bh
shr bh, 4
or bh, 0x30
mov [edi], bh
inc edi
and bl, 0x0F
or bl, 0x30
mov [edi], bl
inc edi
cmp esi, BCD
jne .next
; 指数部の表示
mov ebx, [Exponent]
or ebx, ebx ; 指数値をチェック
je .exit
mov al, "e"
mov [edi], al
inc edi
or ebx, ebx ; 指数値をチェック
jg .e_positive
neg ebx
mov al, "-"
mov [edi], al
inc edi
.e_positive
; 指数値を数字に変換
mov eax, ebx
mov ebx, 100000000
mov ecx, 9
mov esi, 0
.exp_loop:
xor edx, edx ; edx = 0
div ebx ; eax = edx:eax / ebx
push edx ; modulo
test eax, eax
jz .ps1
add al, '0'
mov [edi], al
inc edi
inc esi ; 0 以外の数値出現
jmp short .ps2
.ps1: test esi, esi ; 前の 0 は表示しない
jz .ps2
add al, '0' ; 10/29
mov [edi], al
inc edi
.ps2:
pop eax ; modulo
push edi
xor edx, edx ; edx = 0
xchg eax, ebx ; : save eax
mov edi, 10 ; ebx = ebx / 10
div edi ; :
xchg eax, ebx ; : restore eax
pop edi
loop .exp_loop
.exit
xor eax, eax
mov [edi], al
popa
ret
;---------------------------------------------------------------------
; FPU のスタックトップのレジスタ中の浮動小数点数を
; Sign, BCD, Exponent に分解して格納
; 2進表現の 80bit 浮動小数点数を10進表現に変換
; Float2Ascii から呼び出される
;---------------------------------------------------------------------
Float2BCD: push eax
mov byte [Sign], 0 ; 正の数値を仮定
fldz ; 0 をプッシュ
fcomp st1 ; st0 と st1 を比較して pop
fstsw ax ; フラグをaxに転送
sahf ; ah をフラグレジスタに転送
jz near .zero
push ebx
push edx
.notzero ; ゼロでなければ, 符号を設定
jb .Positive
fabs ; 符号を記録し，絶対値を求める
mov byte [Sign], 1
.Positive ; 符号が正
mov dword [Exponent], 0 ; 指数を 0 で初期化
; 指数の正負を判定
fld1 ; 1.0 をプッシュ
fcomp st1 ; (1.0 - n) 判定が逆になる
fstsw ax
sahf
jbe .PosExp ; １以下なら.NegExp
;----------------------------------------------------------
; 負の指数の場合
; ここで 0 から 1 の間の数値が得られている．
; 適当な10のべき乗を乗算して数値fが 1 < f < 10 となるようにする．
mov ebx, TblSize*10 ; 表の直後に設定
mov edx, 8192
.NegExp
sub ebx, 10 ; 表の１つ前に移動
shr edx, 1
test edx, edx ; べき数が0なら次へ
jz .mul10
fld tword [NegTbl + ebx]; 10の負のべき乗
fcomp st1 ; st0(10^n) - st1
fstsw ax ; フラグをaxに転送
sahf ; ah をフラグレジスタに転送
jbe .NegExp ; 小さい間繰り返し 11/1
sub [Exponent], edx ; 指数 -= ExpTbl [di]
fld tword [PosTbl + ebx]; 10の正のべき乗を乗算
fmulp st1, st0
jmp .NegExp ; さらに繰り返し
.mul10
fld tword [PosTbl] ; 10を乗算
fmulp st1, st0
dec dword [Exponent] ; 指数 -= 1
jmp .range1_10
;----------------------------------------------------------
; 正の指数の数値を扱う
; 適当な10のべき乗を除算して数値fが 1< f <10 となるようにする．
.PosExp
mov ebx, TblSize*10 ; 表を後ろから参照
mov edx, 8192
.PosExp2
sub ebx, 10 ; 表の１つ前に移動
shr edx, 1
test edx, edx
jz .range1_10
fld tword [PosTbl + ebx]; 10の正のべき乗
fcomp st1 ; st0(10^n) - st1
fstsw ax ; フラグをaxに転送
sahf ; ah をフラグレジスタに転送
ja .PosExp2 ; 大きい間繰り返し
add [Exponent], edx
fld tword [PosTbl + ebx]; 1 x 2^n
fdivp st1, st0 ; 除算 st1/st0
jmp .PosExp2 ; さらに繰り返し
.range1_10
; この時点で数値は 1 <= x < 10 の範囲になっている,
; 10^17 を乗じて(残りの17桁分)整数化し，BCDに変換する
fld tword [MaxDigits]
fmul st1
; 10^18で fbstp に失敗する場合に対処
fld tword [MaxBCD]
fcomp st1
fstsw ax ; フラグをaxに転送
sahf ; ah をフラグレジスタに転送
jae .bcd
fld tword [PosTbl]
fdivp st1, st0
inc dword [Exponent] ; 指数 += 1
.bcd
fbstp tword [BCD] ; stack topをBCD18桁でメモリ10byteに格納
; 変換前にRCで指定される丸めモードで整数化
; BCD には 18桁のBCDとして返る．
pop edx
pop ebx
pop eax
ret
; ゼロの場合は特別，定数としてのゼロを設定
.zero
xor eax, eax
mov [BCD], eax ; BCDの10バイトを0に設定
mov [BCD + 4], eax
mov [BCD + 8], ax
mov dword[Exponent], 0
mov byte[Sign], ' '
pop eax
ret
;---------------------------------------------------------------------
; edi に設定された浮動小数点数の文字列を 80bit で FPU のスタックトップ
; データレジスタに設定する．先頭部分の空白やタブは除いておく.
; 数値が "yyy.xxxx" (y <> 0) の形式のとき
; "yyy" の桁が18桁以上の場合は単に指数を増加させて数字は無視
; このコードは文字を DigitStr 配列に逆順に格納していることに注意.
; リトルエンディアンの場合この方が有利
; 10進表現の浮動小数文字列を 2進表現 80bit 浮動小数に変換
; edi は読み込み済みのバッファの次の位置が返る. 他のレジスタは保存.
;---------------------------------------------------------------------
Ascii2Float: push eax
push ebx
push edx
xor eax, eax
mov ebx, DigitStr
.Clear mov [ebx], eax
add ebx, 4
cmp ebx, DigitStr + 20
jne .Clear
mov edx, -18 ; BCD変換のための指数の初期値
mov ebx, 18 ; 処理すべき残っている有効桁数
mov byte[Sign], 0
call GetChar
jae .SkipZero
je .FractionPart
cmp al, "+" ; 符号付?
je .SkipSign
cmp al, "-"
jne near .ConvErr
mov byte[Sign], 1
.SkipSign: call GetChar
jae .SkipZero
cmp al, "." ; 小数点から始まる?
je .FractionPart
jmp near .ConvErr
.SkipZero
cmp al, 0 ; 先頭部分の 0 を除く
jne .Digits
call GetChar
jb .DecimalPoint
jmp .SkipZero
.Digits:
inc edx ; 小数点以上の数値桁数を調整
test ebx, ebx ; 有効桁数を超えていないか
jz .SkipDigit ; 18桁を越えたら数字はスキップ
mov [DigitStr + ebx], al; 後ろから詰める
dec ebx
.SkipDigit
call GetChar
jb .DecimalPoint
jmp .Digits
.DecimalPoint:
cmp al, "."
jne .BCDConv
;----------------------------------------------------------
; 小数点以下の桁の処理
; 小数点以下の各桁では１桁ごとに指数を減少させる必要がある．
; FPU の BCD変換では小数点以下の値は許されないため．
; 例えば 0.12345 --> 12345e-5 とする．
.FractionPart:
call GetChar
jb .BCDConv
cmp ebx, 18 ; すでに0以外の数値 10/29
jne .Fp2 ; 10/29
cmp al, 0 ; .0000xx をチェック
jne .Fp2
dec edx
jmp .FractionPart
.FractionPart2:
call GetChar
jb .BCDConv
.Fp2: test ebx, ebx ; 有効桁数を超えたか?
jz .FractionPart2 ; 超えていたら無視する
mov [DigitStr + ebx], al; 超えていなければ格納
dec ebx
jmp .FractionPart2
;----------------------------------------------------------
; DigitStr 変数を BCD に変換
; 数字の文字列の数字２つの下位4ビットからBCDの1バイトに変換
.BCDConv:
push eax
push edi
push esi
mov edi, BCD
add edi, 8
mov esi, DigitStr
add esi, 18
.BCDLoop: mov al, [esi]
dec esi
shl al,4
or al, [esi]
dec esi
mov [edi], al
dec edi
cmp esi, DigitStr
ja .BCDLoop
fbld [BCD] ; BCD の値を FPU にロード
pop esi
pop edi
pop eax
;----------------------------------------------------------
; 指数部の処理
.ExponentPart
cmp al, "e" ; e または E のチェック
je .Exponent
cmp al, "E"
jne .SetExponent
;----------------------------------------------------------
; 指数の値を読んで BCD用に調整された指数値に加える
.Exponent
mov byte[ExpSign], 0 ; 正の指数を仮定
xor ebx, ebx
xor eax, eax
call GetChar ; E の後ろの文字取得
jae .ExpNum ; 正の指数
cmp al, "+"
je .EatExpSign
cmp al, "-"
jne near .ConvErr ; +,-,Digit でない
mov byte[ExpSign], 1 ; 指数は負
.EatExpSign:
call GetChar ; 符号の後ろの文字取得
jb near .ConvErr ; 符号の後ろに数字必要
.ExpNum:
imul ebx, 10 ; 指数部の数字を数値化
add ebx, eax
call GetChar
jae .ExpNum
cmp byte[ExpSign], 0 ; 指数は負か？
je .PosExp
neg ebx ; 負にする
.PosExp:
add edx, ebx ; BCD化用の指数に加える
.SetExponent
cmp edx, 4930 ; 指数が範囲外かチェック
jg .ConvErr
cmp edx, -4930
jl .ConvErr
mov ebx, TblSize*10 ; 表の直後に設定
test edx, edx ; 指数が負かどうか
js .NegativeExp
;----------------------------------------------------------
; 指数部は正
shl edx, 19 ; キャリーを使う
.ExpLoop1:
test edx, edx ; 指数が0なら終了
jz .Done
sub ebx, 10
shl edx, 1
jnc .ExpLoop1
fld tword[PosTbl + ebx]
fmulp st1, st0
jmp .ExpLoop1
;----------------------------------------------------------
; 指数部は負
.NegativeExp:
neg edx ; 絶対値を求める
shl edx, 19 ; キャリーを使う
.ExpLoop2: test edx, edx ; 0 なら終了
jz .Done
sub ebx, 10
shl edx, 1
jnc .ExpLoop2
fld tword[PosTbl + ebx]
fdivp st1, st0
jmp .ExpLoop2
.Done:
cmp byte[Sign], 0
je .Exit
fchs
.Exit
pop edx
pop ebx
pop eax
clc
ret
.ConvErr: pop edx
pop ebx
pop eax
stc
ret
;---------------------------------------------------------------------
; 2進10進変換用の定数
;---------------------------------------------------------------------
MaxBCD: dt 9.99999999999999999e17
MaxDigits dt 1.0e+17
NegTbl dt 1.0e-1
dt 1.0e-2
dt 1.0e-4
dt 1.0e-8
dt 1.0e-16
dt 1.0e-32
dt 1.0e-64
dt 1.0e-128
dt 1.0e-256
dt 1.0e-512
dt 1.0e-1024
dt 1.0e-2048
dt 1.0e-4096
PosTbl: dt 10.0
dt 1.0e+2
dt 1.0e+4
dt 1.0e+8
dt 1.0e+16
dt 1.0e+32
dt 1.0e+64
dt 1.0e+128
dt 1.0e+256
dt 1.0e+512
dt 1.0e+1024
dt 1.0e+2048
dt 1.0e+4096
TblSize equ ($-PosTbl)/10
;---------------------------------------------------------------------
; AL の文字が数字かどうかのチェック
; 数字なら整数に変換して AL 返す. 非数字ならキャリーセット
;---------------------------------------------------------------------
IsDigit:
cmp al, "0" ; 数字か？
jb NotDigit
cmp al, "9"
ja NotDigit
sub al, "0" ; 整数に変換
ctc
ret
NotDigit:
stc
ret
;---------------------------------------------------------------------
; バッファの次の数字を AL に取得
; 数字なら整数に変換して AL 返す. 非数字ならキャリーセット
;---------------------------------------------------------------------
GetChar: mov al, [edi]
inc edi ; 次の文字を AL に取得
call IsDigit
ret
%endif
結構複雑ですが, 興味のある方はコメントをたよりにコードを追って見て下さい．
使うだけなら覚えることは:
- edi に浮動小数点数の文字列の先頭アドレスを設定して call Ascii2Float
- edi に 29 バイト以上のバッファの先頭アドレスを設定して call Float2Ascii
です.
使用例
float.inc を使った簡単な例です．100πを計算し，文字列から数値への変換をテスト しています．
;---------------------------------------------------------------------
; testfpu.asm
; 2000/10/03 Jun Mizutani
;---------------------------------------------------------------------
%include "stdio.inc"
%include "float.inc"
section .bss
fstring resb 32
section .data
Num100 dt 100.0
TestNum0 db "000012345678901234567890123.456",0
TestNum1 db "0000.0000000000012345678901234567890123",0
TestNum2 db ".123456789012345678E-7",0
cset db ' Conversion error!', 0x0A, 0
cclear db ' OK!', 0x0A, 0
section .text
global _start
;------------------------------------
; start here
;------------------------------------
_start:
finit ; FPU を初期化
fldpi ; πをpush
fld tword [Num100] ; 100をpush
fmulp st1, st0 ; 100π
call print_float ; 表示
mov edi, TestNum0
call print_atof ; 文字列を浮動小数に変換表示
mov edi, TestNum1
call print_atof ; 文字列を浮動小数に変換表示
mov edi, TestNum2
call Exit
;------------------------------------
; ediを先頭とする文字列を浮動小数に変換して表示
print_atof
mov eax, edi
call OutAsciiZ
call Ascii2Float ; 浮動小数文字列をFPUにセット
jb .error
mov eax, cclear ; 変換失敗
jmp .cend
.error
mov eax, cset ; 変換エラー
.cend
call OutAsciiZ
; FPU スタックトップを表示
print_float
mov edi, fstring
call Float2Ascii
mov eax, fstring
call OutAsciiZ
call NewLine
call NewLine
ret
実行結果
jm:~/ex_asm$ ./asm testfpu jm:~/ex_asm$ ./testfpu 3.14159265358979324e2 000012345678901234567890123.456 OK! 1.23456789012345678e22 0000.0000000000012345678901234567890123 OK! 1.23456789012345678e-12
浮動小数点の計算は最近の CPU にとっては簡単な仕事です.