C99のfma関数の挙動を調査

2014/12/16(火)C99のfma関数の挙動を調査

以前、

HaswellとFMA

の記事で、Haswellに新しく実装されたFMA命令について書きました。某Mさんがたまたまfmaという命令がmath.hにあるのを見つけてくれてその存在を知ったので、どういうものなのか調査してみました。かなりややこしいのですがせっかくなので記憶が薄れないうちに記録を残しておきます。

まず、fmaという命令はC99で追加されていて、使うときにはmath.hが必要らしいです。そこで、

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

int main()
{
        double u, x, y, z;
        int i;

        u = 1;
        for (i=0; i<52; i++) u /= 2;

        x = 1+u;
        y = 1-u;
        z = -1;

        printf("%.15g\n", fma(x, y, z));
}

というプログラムを作って、その挙動を観察してみました。

まずコンパイルの可否。次の表のようになりました。

	オプション無し	-lm	-mfma
gcc4.1	x	x	x
gcc4.4	x	o	x
gcc4.8	x	o	o

この結果を受けて、コンパイル出来た3通りについて、-staticで静的リンクした上でIvyBridge(FMA命令を持たない)とHaswell(FMA命令を持つ)で実行してみました。

	IvyBridge	Haswell
gcc4.4 -lm	0	-4.93038065763132e-32
gcc4.8 -lm	-4.93038065763132e-32	-4.93038065763132e-32
gcc4.8 -mfma	実行できず	-4.93038065763132e-32

u=2^-52として(1+u)(1-u)+(-1)を計算しているので、

FMA命令がちゃんと機能していれば、真の値-2^-104が、
FMAでなくただの乗算+加算だと、1-2^-104が1に丸められ、それに(-1)を足すので0が、

結果となるはずです。

これらの実験結果を見ると、次のような推定が出来ます。

gcc4.1では、C99のfma関数のサポートそのものが無い。
gcc4.4では、C99のfma関数が実装され、libm内にfma関数が存在する。-mfmaによるHaswellのFMA命令の有効化は実装されていない。libm内のfma関数は、CPUがFMA命令を持っていればそれを使い、持っていなければ単にfma(a,b,c)=a*b+cを実行する。
gcc4.8では、-mfmaによるfma命令の有効化が実装された。これを付けてコンパイルすると、CPUがFMA命令を持つことが前提のコードが生成され、main関数中のfmaはlibm呼び出しではなく直接HaswellのFMA命令にコンパイルされる。
gcc4.8で更に注目するべきは、libm呼び出しの場合IvyBridgeでも何故か正しい値が計算されたことである。これは普通ではなかなか出来ないはずで、CPUがFMA命令を持たなかった場合は(例えばmpfrを用いた)正確なFMA命令のエミュレーションを行うような実装が行われたと思われる。

この最後の推測を裏付けるため、10億回実行させた場合の実行時間(秒)を計測してみました。

	IvyBridge	Haswell
gcc4.4 -lm	2.46	2.236
gcc4.8 -lm	127.42	2.828
gcc4.8 -mfma	-	2.214

50倍ほどの速度の低下が認められ、何らかのエミュレーションが行われていると推定出来ます。