x87命令を使ってみる

はじめに

少し前の記事ですが、x87命令であるfsinfcosが遅いとか精度がどうとか話題になっているようです。

参考→サイン、コサインをインテルの CPU で計算すると少しバグっているらしい

x87の歴史については別記事(例えばx86の浮動小数計算とSIMD命令の変遷)に譲りますが、ようするにもともと外付けとして実装されたFPUが、そのうち内部に統合されたのがx87です。そんなわけで、もともとはコプロセッサのことを指しましたが、今は命令セットとして扱われています。内部はスタックマシンになっており、しかもFPUデータレジスタはMMXレジスタと共有なので、一緒に使うとややこしいことになります。

前置きはともかく、さっそくx87命令を使ってみようと思います。

そもそも最近のコンパイラはx87命令吐くの?

吐きます。先の記事のコメント欄にもあるように、long doubleを使うと、デフォルトでx87命令を吐きます。

long double func(long double a, long double b){
  return a + b;
}
$ g++ -O3 -S ld.cpp
func(long double, long double):
LFB0:
        fldt    24(%rsp)
        fldt    8(%rsp)
        faddp   %st, %st(1)
        ret

ここで出てるfldtとかfaddpはx87命令です。

なんでlong double使うとx87命令吐くの?

この辺詳しくないのですが、x86上のコンパイラの多くがlong doubleの精度が80bitであり、x87に合わせてあるからだと思います。sizeof(long double)は16を返しますが、内部表現としては80bit精度の拡張倍精度浮動小数点数になっています(多分)1。なので、long doubleの計算をする際には、そのままの精度で80bit精度の計算ができるx87命令が呼ばれます。

ここでfldtは、メモリの値をレジスタスタックに積む命令ですが、最後の”t”はsource formatが「Extended」、すなわち80bitであることを表しています。これは、メモリ上の拡張倍精度浮動小数点数をそのままレジスタに積みます、という意味です。

x87命令を使ってみる

というわけでx87命令を使ってみましょう。先述したように、x87命令の内部精度は80bitであり、通常扱う倍精度浮動小数点数は64bitですから、変換が必要になります。これはx87のスタックに積む命令が勝手にやってくれます。

例えば、doubleを受け取ってsinの値を返す関数はこんな感じになるでしょう。

double
 __attribute__ ((noinline))
fsin(double a){
  double r;
  asm volatile (
      "movsd %%xmm0, %0\n\t" // 1.
      "fldl %0\n\t"          // 2.
      "fsin\n\t"             // 3.
      "fstpl %1\n\t"         // 4.
      :"=m"(a)
      :"m"(r)
    );
  return r;
}

やってることはこんな感じです。

  1. xmm0に入ってきた引数を、一度メモリにコピーします 2
  2. それをfldl命令でx87スタックレジスタに積みます。ここでfldlは「64bit倍精度を積む」命令で、このとき内部で80bit表現に変換されます。
  3. fsinで、スタックの一番上にあるデータのsin値を計算します。
  4. 計算された値がスタックの一番上に積んであるので、それをメモリに書き戻します。ここで80bit表現が64bit表現にコンバートされます。

まぁその、これみて遅そうだな、と思っていただければ。いや、実際に速度を測定していないのでアレですが。

ちゃんとsinを計算できてるか確認しましょう。

#include <stdio.h>
#include <math.h>

int
main(void){
  const int N = 10000;
  const double s = 2.0*M_PI/static_cast<double>(N);
  double x = 0.0;
  for(int i=0;i<N;i++){
    printf("%f %f\n",x, fsin(x));
    x += s;
  }
}
$ g++ test.cpp
$ ./a.out > test.dat

image0.png

計算できているみたいですね。

まとめ

x87命令を使ってsin関数を計算してみました。x87は内部表現が80bitで通常の倍精度実数と異なるし、外付けFPUだったころの仕様を引きずっていたりしますが、まだちゃんと使えます(ダイの上でどう実装されているかは知りません)。おそらく過去のプログラムの互換性のためだけに残されているのだと思いますが、たまにこうやって使ってあげるとx87命令も喜ぶんじゃないでしょうか。

参考文献

  1. 筆者の手元の記録では、昔のx86系の石+gcc/iccsizeof(long double)は12を返していたようですが、DECのAlpha21264とDECのコンパイラは8を返していたようです。 

  2. ちょっとアホなことをやっていますが、コンパイラの最適化で変なことが起きないようにするためです。もっと賢い方法があるかもしれません。