./index.html ../index.html

アセンブラで浮動小数点演算

気まぐれにアセンブラの勉強～。

Delphiのインラインアセンブラなのは単にデバッグが楽だからなので気にしない。
呼びだし規約はstdcallに準ずる事にします。

初歩…定数を返す

function A_e: Extended; stdcall;
const
  Value: Extended = 0.12345678901234567890123;
asm
  fld Value
end;

function A_d: Double; stdcall;
const
  Value: Double = 0.12345678901234567890123;
asm
  fld Value
end;

function A_s: Single; stdcall;
const
  Value: Single = 0.12345678901234567890123;
asm
  fld Value
end;

begin
  WriteLn(FloatToStrF(A_e, ffFixed, 18, 18));
  WriteLn(FloatToStrF(A_d, ffFixed, 18, 18));
  WriteLn(FloatToStrF(A_s, ffFixed, 18, 18));
end.

実行結果

0.123456789012345679
0.123456789012345677
0.123456791043281555

戻り値は、FPUのスタックというものに積んで返します。 FPUのスタックは、コントロールワードで決定された精度を持ちます。 Delphiの場合、何もしなくてもスタートアップルーチンがExtended精度に設定してくれます。 つまり、実は戻り値の型の精度はアセンブラ上では意味がありません。

メモリから値を取り出してスタックへ積むには、fldを使います。

結果の精度も見て下さい。 FloatToStrFは18桁まで表示可能なのですが、一番近い値を示しているのはExtendedです。

関係ないと書いたばっかですが、Delphiのインラインアセンブラ賢いですね～。
定数の型をちゃんとわかっていて、それぞれ tword ptr, qword ptr, dword ptr 扱いしてくれます。

足し算

function A_e(const A, B: Extended): Extended; stdcall;
asm
  fld A
  fld B
  fadd
end;

function A_d(const A, B: Double): Double; stdcall;
asm
  fld A
  fld B
  fadd
end;

function A_s(const A, B: Single): Single; stdcall;
asm
  fld A
  fld B
  fadd
end;

begin
  WriteLn(FloatToStrF(A_e(PI, Exp(1)), ffFixed, 18, 18));
  WriteLn(FloatToStrF(A_d(PI, Exp(1)), ffFixed, 18, 18));
  WriteLn(FloatToStrF(A_s(PI, Exp(1)), ffFixed, 18, 18));
end.

実行結果

5.859874482048838470
5.859874482048838210
5.859874486923217770

FPUのスタックは、スタックですから、値のロードを繰り返すと、順番に積まれて、最後に積んだのが一番上に来ます。

faddは、スタックの上二つの値を取り出して、足した結果を積みます。

引き算

function A_e(const A, B: Extended): Extended; stdcall;
asm
  fld A
  fld B
  fsub
end;

begin
  WriteLn(FloatToStrF(A_e(PI, Exp(1)), ffFixed, 18, 18));
end.

実行結果

0.423310825130748003

浮動小数点型(Extended, Double, Single)というのは、メモリとFPUのスタックの間でのやりとりでのみ意味を持ちますが、Delphiのインラインアセンブラは前述のように自動判別してくれます。 コードが記述上変わらないのに三通り書くのが馬鹿らしくなってきたので、以後理由が無ければExtendedのみに統一します。

fsubは、最後に積んだ値を、その前に積んだ値から引きます。

負符号

function A_e(const X: Extended): Extended; stdcall;
asm
  fld X
  fchs
end;

begin
  WriteLn(FloatToStrF(A_e(1), ffFixed, 18, 18));
end.

実行結果

-1.000000000000000000

fchsはスタックトップの符号を反転させます。

剰余算

function A_e(const A, B: Extended): Extended; stdcall;
asm
  fld A
  fld B
  fmul
end;

function B_e(const A, B: Extended): Extended; stdcall;
asm
  fld A
  fld B
  fdiv
end;

begin
  WriteLn(FloatToStrF(A_e(2, Sqrt(2)), ffFixed, 18, 18));
  WriteLn(FloatToStrF(B_e(2, Sqrt(2)), ffFixed, 18, 18));
end.

実行結果

2.828427124746190100
1.414213562373095050

fmulがかけ算、fdivが割り算。

定数

function Zero: Extended; stdcall;
asm
  fldz
end;

function One: Extended; stdcall;
asm
  fld1
end;

function Pi: Extended; stdcall;
asm
  fldpi
end;

function Log_10_2: Extended; stdcall;
asm
  fldlg2
end;

function Log_2_10: Extended; stdcall;
asm
  fldl2t
end;

function Log_e_2: Extended; stdcall;
asm
  fldln2
end;

function Log_2_e: Extended; stdcall;
asm
  fldl2e
end;

begin
  WriteLn(FloatToStrF(Zero, ffFixed, 18, 18));
  WriteLn(FloatToStrF(One, ffFixed, 18, 18));
  WriteLn(FloatToStrF(Pi, ffFixed, 18, 18));
  WriteLn(FloatToStrF(Log_10_2, ffFixed, 18, 18));
  WriteLn(FloatToStrF(Log_2_10, ffFixed, 18, 18));
  WriteLn(FloatToStrF(Log_e_2, ffFixed, 18, 18));
  WriteLn(FloatToStrF(Log_2_e, ffFixed, 18, 18));
end.

実行結果

0.000000000000000000
1.000000000000000000
3.141592653589793240
0.301029995663981195
3.321928094887362350
0.693147180559945309
1.442695040888963410

よく使うような定数は、コード上に埋めこまなくても、使用可能です。 fldXXは、定数をスタックに積みます。 (fldpiとか見ると、C言語に限った話ではありませんが、#define M_PIの類が馬鹿らしく思えてきません？…機械語って偉大だ)

式

ここらで、具体的な式、例えば…適当ですけど、 1 / X + 2 * Y を計算してみます。

function A_e(const X, Y: Extended): Extended; stdcall;
const
  _2: Extended = 2;
asm
  fld1
  fld X
  fdiv
  fld _2
  fld Y
  fmul
  fadd
end;

begin
  WriteLn(FloatToStrF(A_e(2, Pi), ffFixed, 18, 18));
end.

実行結果

6.783185307179586480

FPUはスタックを用いた計算を行ないますので、まず 1 / X + 2 * Y を逆ポーランド記法に直します。 すると 1 X / 2 Y * + 。後はこの通りに命令を並べるだけです。

逆ポーランド記法がわからない人は、拙作「いい電卓」で普通の式から変換できます。(宣伝)

代入

var
  R: Extended;

procedure A_e; stdcall;
asm
  fld1
  fstp R
end;

begin
  A_e;
  WriteLn(FloatToStrF(R, ffFixed, 18, 18));
end.

実行結果

1.000000000000000000

fstpでスタックから値をポップしてメモリ上に代入できます。

ここではDelphiのインラインアセンブラだから型を省略できていますが、普通はtword ptr, qword ptr, dword ptrを書かなければなりません。

整数→実数

function IntToFloat(X: Integer): Extended; stdcall;
asm
  fild X
end;

function UnsignedToFloat(X: Cardinal): Extended; stdcall;
var
  Work: Int64;
asm
  mov Work.Int64Rec.Hi, 0
  mov eax, X
  mov Work.Int64Rec.Lo, eax
  fild Work
end;

begin
  WriteLn(FloatToStrF(IntToFloat(-1), ffFixed, 18, 18));
  WriteLn(FloatToStrF(UnsignedToFloat($ffffffff), ffFixed, 18, 18));
end.

実行結果

-1.000000000000000000
4294967295.000000000000000000

fildで整数をスタックに積めますが、この命令は、常に符号付きの扱いをしてくれます。 なので、符号無し整数の場合は、一旦、一回り大きな整数型を経由する必要があります。

実数→整数

function Round(const X: Extended): Integer; stdcall;
var
  R: Integer;
asm
  fld X
  fistp R
  mov eax, R
end;

function Trunc(const X: Extended): Integer; stdcall;
var
  SaveCW, WorkCW: Word;
  R: Integer;
asm
  fld X
  fnstcw SaveCW
  fnstcw WorkCW
  or WorkCW, $0c00
  fldcw WorkCW
  fistp R
  fldcw SaveCW
  mov eax, R
end;

function Ceil(const X: Extended): Integer; stdcall;
var
  SaveCW, WorkCW: Word;
  R: Integer;
asm
  fld X
  fnstcw SaveCW
  fnstcw WorkCW
  and WorkCW, $f3ff
  or WorkCW, $0800
  fldcw WorkCW
  fistp R
  fldcw SaveCW
  mov eax, R
end;

function Floor(const X: Extended): Integer; stdcall;
var
  SaveCW, WorkCW: Word;
  R: Integer;
asm
  fld X
  fnstcw SaveCW
  fnstcw WorkCW
  and WorkCW, $f3ff
  or WorkCW, $0400
  fldcw WorkCW
  fistp R
  fldcw SaveCW
  mov eax, R
end;

begin
  WriteLn(Round(1.7), '/', Round(-1.7), '/', Round(1.3), '/', Round(-1.3));
  WriteLn(Trunc(1.7), '/', Trunc(-1.7), '/', Trunc(1.3), '/', Trunc(-1.3));
  WriteLn(Ceil(1.7), '/', Ceil(-1.7), '/', Ceil(1.3), '/', Ceil(-1.3));
  WriteLn(Floor(1.7), '/', Floor(-1.7), '/', Floor(1.3), '/', Floor(-1.3));
end.

実行結果

2/-2/1/-1
1/-1/1/-1
2/-1/2/-1
1/-2/1/-2

fistpでスタックから値を取り出し整数としてメモリへ書き出せますが、小数点以下の扱い方が、銀行家の丸め、ゼロ方向へ切り捨て、+∞方向へ切り上げ、-∞方向へ切り下げ、の4種類あります。

デフォルトが銀行家の丸めなので、それ以外は、FPUのコントロールワードを一時的に書き換えてやります。

銀行家の丸めとは、原則四捨五入ですが、丁度0.5の時は、偶数側を選ぶ…らしいです。 Delphiのヘルプには「JIS ハンドブック標準化「数値の丸め方」（Z8401-1961）を参照してください」なんて書いてありますが…。
普通の四捨五入が欲しい時は、当たり前ですが、0.5を足して切り捨てでおっけーです。

コントロールワードについては、僕もよくわかって無いので詳しい説明は省かせて下さい。 fnstcwで取得、fldcwで設定です。

比較

function G(const X, Y: Extended): Boolean; stdcall;
asm
  fld Y
  fld X
  fcompp
  fnstsw ax
  sahf
  seta al
end;

begin
  WriteLn(G(1, 2), '/', G(2, 1));
end;

実行結果

FALSE/TRUE

fcomppで比較を行ないますが、比較結果がFPUのステータスワードという特別な場所に格納されるため、fstsw axでステータスワードをaxに取得、sahfでax中の比較結果をフラグレジスタに映します。 後は、通常のja等で分岐したり、seta等で結果を得たりできるという寸法です。
a, b, ae, be, e, ne系は動作を確認しましたが、g, ge, l, leは上手く行きませんでした。 ですので、前者を使っておけば問題無い様です。

後から積んだ方が、比較の基準になるみたいなので、注意。ここだけ今までの感覚と逆なんですね。

各種関数

function Sin(const X: Extended): Extended; stdcall;
asm
  fld X
  fsin
end;

function Cos(const X: Extended): Extended; stdcall;
asm
  fld X
  fcos
end;

function Tan(const X: Extended): Extended; stdcall;
asm
  fld X
  fptan
  fstp st(0)
end;

function ArcTan(const Y: Extended): Extended;
asm
  fld Y
  fld1
  fpatan
end;

function ArcTan2(const Y, X: Extended): Extended;
asm
  fld Y
  fld X
  fpatan
end;

function Abs(const X: Extended): Extended; stdcall;
asm
  fld X
  fabs
end;

function Sqrt(const X: Extended): Extended; stdcall;
asm
  fld X
  fsqrt
end;

function Log(const Base, X: Extended): Extended; stdcall;
asm
  fld1
  fld X
  fyl2x
  fld1
  fld Base
  fyl2x
  fdiv
end;

function E: Extended; stdcall;
asm
  fld1 { 2 ^ (log 2 e) }
  fldl2e
  fld1
  fsub
  f2xm1
  fld1
  fadd
  fscale
  fstp st(1)
end;

function Log2(const X: Extended): Extended; stdcall;
asm
  fld1
  fld X
  fyl2x
end;

function Pow2(const X: Extended): Extended; stdcall;
asm
  fld X
  fld1
  fscale
  fstp st(1)
end;

begin
  WriteLn(FloatToStrF(Sin(Pi / 6), ffFixed, 18, 18));
  WriteLn(FloatToStrF(Cos(Pi / 6), ffFixed, 18, 18));
  WriteLn(FloatToStrF(Tan(Pi / 6), ffFixed, 18, 18));
  WriteLn(FloatToStrF(ArcTan(1), ffFixed, 18, 18));
  WriteLn(FloatToStrF(ArcTan2(1, Sqrt(3)), ffFixed, 18, 18));
  WriteLn(FloatToStrF(Abs(-2), ffFixed, 18, 18));
  WriteLn(FloatToStrF(Sqrt(2), ffFixed, 18, 18));
  WriteLn(FloatToStrF(Log(2, 8), ffFixed, 18, 18));
  WriteLn(FloatToStrF(E, ffFixed, 18, 18));
  WriteLn(FloatToStrF(Pow2(Log2(2) * 3), ffFixed, 18, 18));
end.

実行結果

0.500000000000000000
0.866025403784438647
0.577350269189625764
0.785398163397448310
0.523598775598298873
2.000000000000000000
1.414213562373095050
3.000000000000000000
2.718281828459045240
8.000000000000000000

eを得るやつだけ、やけに苦労していますが、2 ^ X - 1を計算するf2xm1が、Xの小数部分しか見てくれない(X > 1の時何もしない)のが原因です。 fscaleは、スタックトップに、2 ^ (スタックのひとつ前に積まれた数) を乗算します。単に2をかけ算しても同じですが、2を用意する方が面倒(定数にしてメモリから読むか、fld1, fld1, faddか)なので、こんな感じです。 fscaleはスタックのひとつ前に積まれた数を捨ててくれないので、fstp st(1)で捨ててやります。
要するに、実際には (2 ^ ((log 2 e) - 1) - 1 + 1) * 2 ^ 1 を計算しています。

Tanでも使っていますが、fstpは、fstp st(n)と書くと、スタックのn番目の値を捨てる事ができます。

これぐらいの基本的な関数があれば、後は組み合わせです。 例えば、最後の例では、Log 2 X を求める方法と、2 ^ Xを求める方法を組み合わせ、2 ^ 3 を計算しています。