この記事では融合積和 (fused multiply-add; 以下もっぱら FMA) について扱います。
この記事は元々「その辺で提供されている fma 関数の実装が正しいかチェックしてみた」記事だったのがだんだん守備範囲を広げていったものなので、FMAを単に使うだけではなく、実装する側の視点が多く含まれています。
筆者が実装したHaskell向けのFMA関数はこちらにあります：
FMAとは、FMAを含む規格
融合積和（fused multiply-add; FMA）とは、積と和が合体した演算です。数式で書けば $xy+z$ です。この際、丸め（浮動小数点数で正確に表せない数をそれに近い浮動小数点数で代用すること）を1回しか行わないため、「積→和」の2回に分けて演算するよりも精度が向上します。
FMAは精度が向上するだけでなく、対応したCPUやその他演算器であれば積和を普通に（乗算→加算の2命令で）計算するよりも高速に計算できます。
FMAは行列乗算、ベクトルの内積や、多項式の評価などで活用できます。
C言語（C99）
高級言語の中でFMAを規格に取り込んだのは、筆者の知る限りC言語が初です。C言語へは1999年の改訂版（C99）で、 fma, fmaf, fmal の3つの関数が導入されました。
#include <math.h>
double fma(double x, double y, double z);
float fmaf(float x, float y, float z);
long double fmal(long double x, long double y, long double z);
この他、関係するマクロとして FP_FAST_FMA, FP_FAST_FMAF, FP_FAST_FMAL が規定されています。これらは、コンパイル先の環境でFMAを通常の積和と同等かより高速に計算できる場合（i.e. CPUがFMA命令を持っている場合）に <math.h> で定義されます。
後述しますが、C言語のソースコード中で a * b + c と書いてもFMA命令が利用されるとは限りません。
一般には a * b + c とFMAは意味が違う（結果が異なる）ので、「値が変わる可能性のある最適化」が有効でないとFMAを利用するわけにはいかないのです。C標準では、この辺に関係する #pragma として #pragma STDC FP_CONTRACT が規定されています。
IEEE 754
FMA演算はIEEE 754-2008以降で fusedMultiplyAdd として規定されています。
FMAのコーナーケース
浮動小数点数の中には普通の実数として解釈できないやつがちょいちょいあります。それらに対するFMAの挙動がどう規定されているか確認しておきます。
0の符号
IEEE 754の浮動小数点数には0の符号の区別があるので、各演算は「結果が（数学的には）0の時にどちらの符号を選ぶか」を決める必要があります。
fma の結果が浮動小数点数として ±0 の場合、
- 数学的な（丸めを行わない、正確な）値として $xy+z\ne 0$ ならば
fma(x, y, z)はその符号を保つ。 - 数学的な（丸めを行わない、正確な）値として $xy+z=0$ ならば、
x * yとzの浮動小数点数としての足し算の規則に従う。つまり、- 数学的な値として $(xy,z)=(0,0)$ の場合
-
x * y = +0,z = +0の場合、fma(x, y, z) = +0 + +0 = +0 -
x * y = -0,z = +0の場合、fma(x, y, z) = -0 + +0 = +0（丸めモードが「負の無限大方向」以外の場合） -
x * y = +0,z = -0の場合、fma(x, y, z) = +0 + -0 = +0（丸めモードが「負の無限大方向」以外の場合） -
x * y = -0,z = -0の場合、fma(x, y, z) = -0 + -0 = -0 -
x * yの符号はxの符号とyの符号のxorです。
-
- 数学的な値として $(xy,z)\ne(0,0)$ の場合、
fma(x, y, z) = +0（丸めモードが「負の無限大方向」以外の場合）
- 数学的な値として $(xy,z)=(0,0)$ の場合
という規則によって 0 の符号が定められます。入力に無限大やNaNが含まれる場合は浮動小数点数としての演算結果は 0 にはならないのでここでは考える必要はありません。
無限大
z = ±∞ の時、 fma(x, y, z) は、
-
xとyのいずれも有限であればz（無限大） -
xとyのいずれかが無限大、NaNであればx * y + zに従う
となります。
後でまた触れますが、たとえ浮動小数点数としての積 x * y がオーバーフローにより無限大となっても、数学的な積が有限であればFMAの計算においては有限として扱われることに注意してください。
NaN
入力のいずれかがNaNの場合は結果はquiet NaNになります。
浮動小数点例外については「入力にquiet NaNが含まれる場合は浮動小数点例外は起こらない」という原則に反して、 fma(0, ±∞, NaN) が 0 * (±∞) に相当する例外（Invalid operation）を投げる可能性があります（投げられない可能性もある）。
FMAを使うコードの書き方
コード中の x * y + z は最適化でFMA演算になるか
x * y + z を fma(x, y, z) にするのは「値を変えうる変更」であり、単にコンパイラーの最適化レベルを上げるだけではそういう最適化が行われるとは限りません。なので、そういう最適化が行われなくてもコンパイラーを責めないでください。
確実にFMAで計算させたい場合は、 fma 関数を明示的に呼び出すべきです。まともなコンパイラーなら、ターゲットの命令セットがFMA命令を持っていれば fma への直接呼び出しは最適化によってFMA命令の実行に置き換えられ、関数呼び出しは発生しません。
#include <math.h>
// coeff[0] + x * (coeff[1] + x * (... * (coeff[n-2] + x * coeff[n-1]) ...))
double horner(size_t n, double coeff[], double x)
{
double acc = 0;
for (size_t m = n; m > 0; --m) {
// たとえ最適化が有効でもFMA命令が使われるとは限らないし、
// 処理系次第ではFMA命令が使われるかもしれない
acc = coeff[m - 1] + x * acc;
}
return acc;
}
double horner_fma(size_t n, double coeff[], double x)
{
double acc = 0;
for (size_t m = n; m > 0; --m) {
// 必ずFMA演算が使われる
acc = fma(x, acc, coeff[m - 1]);
}
return acc;
}
コンパイラーに「多少値が変わってもいいからFMA命令が使えるなら x * y + z をFMAで計算してくれ」と伝える場合は #pragma STDC FP_CONTRACT を指定します。
筆者が試したところ、Clangではデフォルトで x * y + z にはFMAを使わないが、 #pragma STDC FP_CONTRACT ON を指定し、かつFMA命令が利用できる場合は x * y + z がFMA命令に変換されました。
double horner_contract(size_t n, double coeff[], double x)
{
#pragma STDC FP_CONTRACT ON
double acc = 0;
for (size_t m = n; m > 0; --m) {
// FMA命令が利用可能で、処理系が #pragma STDC FP_CONTRACT に対応していればFMA命令が利用される可能性が高い
acc = coeff[m - 1] + x * acc;
}
return acc;
}
double horner_no_contract(size_t n, double coeff[], double x)
{
#pragma STDC FP_CONTRACT OFF
double acc = 0;
for (size_t m = n; m > 0; --m) {
// 処理系が #pragma STDC に対応していれば、たとえFMA命令が利用可能であってもFMA命令は使われない
acc = coeff[m - 1] + x * acc;
}
return acc;
}
一方、GCCでは、デフォルトで積和がFMAに変換されました。そして、バージョン10.1の時点では #pragma STDC FP_CONTRACT はどうやら未実装のようで、単に無視されました。とはいえ、そういう最適化を行うかどうかはコンパイルオプションで指定できるので、コンパイルオプション -ffp-contract={fast,on,off} を指定するか、対応する #pragma GCC あるいは __attribute__ で制御することは可能です。
MSVCの場合は、/O2 /arch:AVX2 でFMA命令を有効化した上で /fp:fast もしくはそれに相当する #pragma を指定する必要があるようです。MSVCには #pragma fp_contract という #pragma があるのに標準の #pragma STDC FP_CONTRACT に対応していないのは興味深いです（セマンティクスが違うとかそういう話なのでしょうが）。
Intel C++ Compilerの対応状況も気になりますが、筆者はICCを使える環境にないので確かめられません。
FMA命令が利用可能かコンパイル時に判断する
x86系がターゲットの場合、FMA命令が登場したのは比較的最近なので、コンパイラーのデフォルトではFMA命令は利用されません。FMA命令を利用するには、 -mfma 等のオプションを指定する必要があります。
「コンパイラーの設定でFMA命令を有効化できているか不安だ」「FMA命令の代わりに遅いソフトウェア実装が使われるくらいならコンパイルに失敗した方がマシだ」という場合は、C言語の場合は FP_FAST_FMA マクロ（double の場合。float の場合は FP_FAST_FMAF）を利用することができます。すでに書いたように、 FP_FAST_FMA マクロは fma の呼び出しが x * y + z と同等かそれ以上に速い場合に定義されます。つまり、ターゲットがFMA命令を持っているかの判定に使えます。
x86系に限れば __FMA__ 等のマクロも定義されていたりしますが、C標準で定義されているのは FP_FAST_FMA です。
#include <math.h>
// 「遅いFMAが使われるくらいならコンパイルエラーになった方がマシだ」という人のためのコード
#if !defined(FP_FAST_FMA)
#error FMA instruction must be available
#endif
double horner_fast_fma(size_t n, double coeff[], double x)
{
double acc = 0;
for (size_t m = n; m > 0; --m) {
acc = fma(x, acc, coeff[m - 1]); // 必ずFMA命令が使われる
}
return acc;
}
ただし、コンパイラーとlibcの組み合わせによってはFMA命令が使える場合でも FP_FAST_FMA が定義されない場合があるので注意してください。
fused multiply-sub
ハードウェアによっては、$xy + z$ を計算するfused multiply-addの他にfused multiply-subtract $x y - z$ や、積の符号を反転させた $- x y - z$, $-x y + z$ を計算する命令を持っていることがあります。というかここで紹介するx86とAArch64の両方がそうです。
一方、C言語やIEEE 754で規定されているのはfused multiply-addのみです。fused multiply-subtract等の亜種を利用するには、環境依存な方法を使うしかないのでしょうか？
実は、FMSUB(x, y, z) = FMADD(x, y, -z) という関係が常に成り立つので、C言語なら fma(x, y, -z) や fma(-x, y, z) とかけば符号の反転とFMA命令が融合されてFMSUB等の命令が使用されることを期待できます。
$\pm xy\pm z$ のそれぞれの融合演算を持っている環境なら、 ±fma(±x, ±y, ±z) は一命令にコンパイルできます。具体的には、
<fused x*y+z> = fma(x, y, z) = fma(-x, -y, z) = -fma(-x, y, -z) = -fma(x, -y, -z): fused multiply-add
<fused x*y-z> = fma(x, y, -z) = fma(-x, -y, -z) = -fma(-x, y, z) = -fma(x, -y, z): fused multiply-sub
<fused -x*y+z> = -fma(x, y, -z) = -fma(-x, -y, -z) = fma(-x, y, z) = fma(x, -y, z): negated FMSUB, fused negated multiply-add
<fused -x*y-z> = -fma(x, y, z) = -fma(-x, -y, z) = fma(-x, y, -z) = fma(x, -y, -z): negated FMA, fused negated multiply-sub
という感じです。
逆に言えば、FMA命令と符号反転があれば他の変種は表せるので、命令セットとしてFMAの変種を用意しておくのは冗長なのですが、あえてそれらの変種が用意されているのはコードサイズの削減という目的が考えられます。あるいは、NaNの符号ビットやペイロードを気にしているのかもしれません。
__builtin_fma
GCCやClangは __builtin_fma という組み込み関数を提供しています。これは通常の fma 関数とどう違うのでしょうか？
ぶっちゃけると、これらは同じです。つまり、__builtin_fma 関数の呼び出しは fma 関数の呼び出しに置き換えられますし、 __builtin_fma に対する最適化は fma にも適用されます。
唯一異なるのは、 -fno-builtin オプションを指定した場合に fma の方には最適化が効かなくなることです。この最適化には
- 定数畳み込み
- ハードウェアFMAが利用可能なアーキテクチャー（x86の場合は要
-mfma）で関数呼び出しが当該命令列に置き換えられる - （試した感じではClangのみ）
fma(x, 1.0, y)がx + yに置き換えられる
が含まれます。「ハードウェアFMAのために __builtin_fma を使う」というのは意味のない行為で、ハードウェアFMAが使いたかったら（x86系の場合） -mfma を指定するべきなのです。
FMAの実装の落とし穴
Cやその他言語での fma 関数を実装する側のことを考えてみましょう。
まず、ターゲットとなるアーキテクチャが必ずFMA命令を持っている場合（例：AArch64）はその命令を直接呼び出せばOKです。ハードウェア実装されたFMA命令はまず間違いなく正しく実装されています（さもなくば「CPUのバグ」として報告されます）。
一方、ターゲットがFMA命令を持っているとは限らない場合（例：x86）は、何らかの方法でFMAをソフトウェア実装する必要があります。標準ライブラリーを実装する人はFMAのことを隅々まで把握しているとは限らないため、そういう人が書いた fma 関数の実装にはバグが潜んでいる可能性があります。
標準ライブラリーの実装者が正しく fma 関数を実装できるように、あるいは、利用者が自分の使う fma 関数のバグを早期発見できるように、ここではFMAの典型的な「間違い方」をいくつか挙げてみます。
オーバーフロー
浮動小数点数はかなり（指数部が）広い範囲の数を表せますが、それでも表せる範囲には限りがあります。
例えば、よく使われる倍精度（double 型）では絶対値が $2^{-1074}$ 未満の数を表現することはできませんし、絶対値が $2^{1024}$ 以上1の数は表現できません。
さて、FMAも浮動小数点演算ですから、結果（の絶対値）があまりにも大きいとオーバーフローを起こして無限大を返します。
しかし、浮動小数点数での乗算 x * y がオーバーフローしても FMA fma(x, y, z) がオーバーフローするとは限りません。
具体的には以下の場合です。
-
zが無限大の場合 -
x * yがギリギリでオーバーフローするが、正確な値 $xy + z$ の丸めの際にオーバーフローが起こらない場合
1.は fma(0x1p1000, 0x1p1000, INFINITY) みたいなやつです。一応注釈しておくと、 0x1p1000 というのは $2^{1000}$ のことで、これを二乗すると $2^{2000}$ となって倍精度で表せなくなります。この場合は結局無限大を返しますが、浮動小数点例外のことを考えると違いがあります。
また、例を少し変えて fma(-0x1p1000, 0x1p1000, INFINITY) とするとどうでしょうか。x * y + z の場合に同じような計算をすると x * y がオーバーフローにより負の無限大となり、正の無限大と足した時にNaNが発生します。しかし、FMAの計算では途中計算は無限の精度を持つものとして行われるため、積 $xy=-2^{2000}$ はオーバーフローしません。この場合のFMAは、 $-2^{2000}$ という有限の数と正の無限大を足すので、NaNではなく正の無限大が返るべきです。
2.の例は fma(0x1p512, 0x1p512, -0x1p971) です。x = 0x1p512 と y = 0x1p512 の積は 0x1p1024 つまり $2^{1024}$ となって倍精度では表せませんが、その後 z を足してやると 0x1p1024 - 0x1p971 = 0x1.ffff ffff ffff fp1023 という有限の倍精度浮動小数点数で表せる数が得られます。
float に対するFMA
float に対する積和 $xy+z$ は、double を使って
(float)((double)x * (double)y + (double)z)
と計算した方が精度が上がりそうな気がします。何だったら、FMAの実装だってこれでいいんじゃないのという気がしてしまいます。
ですが、これはダメです。通常の積和ではなく（double に対する）FMAを使ってもダメです。
// ダメな実装の例
float fmaf(float x, float y, float z)
{
return (float)((double)x * (double)y + (double)z); // ダメ
// return (float)fma((double)x, (double)y, (double)z); // これもダメ
}
この場合、 float から double へのキャストは正確に行われます。また、単精度（float）の精度は24ビット、倍精度（double）は精度53ビットなので、単精度どうしの積（これはせいぜい48ビットあれば表せる）は double を使えば正確に計算できます。なので、積 (double)x * (double)y も正確に計算されます。
丸め（による誤差）が発生する可能性があるのは、 (double)z と足す部分です。さらにその後、 double として得られた値を float にキャストする場面でも丸めが行われます。
つまり、上記のコードでは
<xy+z の正確な値> → double → float
という風に「2段階で」丸めが行われるわけです。一方、FMAで要請されているのは「ただ1回の」丸めです。
<xy+z の正確な値> → float
2段階で丸めるのと、ただ1回丸めるのでは結果は一致するのでしょうか？残念ながら、一般には否です。
例として、最近接偶数丸めで、「正確な値」が t = 0x1.0000 02ff ffff fcp52 だった場合を考えてみましょう。16進ネイティブではない人向けに $t$ の2進表記も書いておきます。
t = 0x1. 0000 02ff ffff fc p52
= 0b1.0000 0000 0000 0000 0000 0010 1111 1111 1111 1111 1111 1111 1111 1100 p52
^ ^-- 53ビット目
24ビット目
この実数 $t$ を単精度（精度24ビット）へ丸めることを考えましょう。25ビット目以降を切り捨てるか、25ビット目以降を切り上げるかを決める必要があります。25ビット目以降を見ると、切り捨てた方が元の値に近いことがわかります。よって、 $t$ を単精度浮動小数点数へ丸めた結果は 0x1.0000 02p52 となります。
一方、$t$ を一旦倍精度（精度53ビット）へ丸めることにすると、54ビット目以降は切り上げるのが適切です。この際、盛大に繰り上がりが発生して、倍精度で表した結果は t' = 0x1.0000 03p52 となります。この $t'$ をさらに単精度に丸めると、25ビット目が1でそれ以降は0なので、24ビット目が0になるように丸めが行われます（最近接「偶数」丸め）。よって、$t$ を double 経由で float へ丸めた結果は 0x1.0000 04p52 となります。
ということで、最近接丸めを2段階でやるとマズいということがお分かり頂けたでしょうか？
小難しく書きましたが、一言で言えば「四捨五入を2段階でやるとマズいぞ」ということです。1.49 を整数へ丸めるのに2段階で四捨五入（1.49 → 1.5 → 2）したらおかしいですよね？10進じゃなくて2進浮動小数点数でも同じことなのです。
ちなみに、この問題が発生するのは「最近接丸め」の場合で、「方向丸め (directed rounding)」の場合は丸めが2段階でも関係ありません。
この「2段階で丸めるとヤバイ」のは一般論ですが、それがFMAの場合に実際に問題となることを確認してみましょう。具体例です。
x = 0x1.fffffep23, y = 0x1.000004p28, z = 0x1.fep5
の時、 $xy+z$ はさっきの $t$ となり、正しい fmaf(x, y, z) の値 0x1.000002p52 と double で計算した積和の値 0x1.000004p52 に違いが出ます。
些末な違いだと思いましたか？あなたにとって些末な違いでも、別の誰かにとっては重大な違いかもしれないのです。CPUメーカーがこういう間違いをやらかしたら最悪リコールですよ（最近はマイクロコードの修正で済むのかもしれませんが）。
ちなみに、 float に対するFMAを実装する際に double が全く活用できないかというとそういうわけでもなくて、積の部分を double で計算することによってアルゴリズムを簡単にすることができます。積を計算した後の和の部分はどうにか工夫する必要があります。
高精度な型を使った手抜き
float の場合と同じ理屈で、多少高精度な long double 型があったからといって double に対するFMAを
double fma(double x, double y, double z)
{
return (double)((long double)x * (long double)y + (long double)z); // ダメ
}
と実装するのはダメです。long double の仮数部が159ビット以上ある環境なら問題ないのかもしれませんが……（検証求む）。
ハードウェアFMA
x86系
x86系では、比較的最近（2010年代の途中）までFMA命令がありませんでした。Intelの場合は実装されたのはHaswellですが、AMDの方は一旦Intelとは互換性のない形で実装されたそうです。AMD系のCPUでIntelと互換性のあるFMA命令が実装されたのはRyzen以降のようです（AMDに関しては筆者は実際のマシンを持っていないので、Wikipedia等の情報です）。
この記事ではx86系のFMA命令としては、Intelが最初に実装したFMA命令（3オペランドの方）のみを紹介します。
IntelのFMA命令セットには、fused multiply-addの他にfused multiply-subや符号を反転させた変種、それからSIMDでレーンごとに異なる演算をするやつが含まれています。
| mneomnic | 表す演算 | 備考 |
|---|---|---|
| VFMADD{132,213,231}[PS][DS] | $xy+z$ | |
| VFMSUB{132,213,231}[PS][DS] | $xy-z$ | |
| VFNMADD{132,213,231}[PS][DS] | $- xy + z$ | fused negative multiply-add |
| VFNMSUB{132,213,231}[PS][DS] | $- xy - z$ | fused negative multiply-subtract |
| VFMADDSUB{132,213,231}P[DS] | $x_iy_i-(-1)^iz_i$ ($i=0,1,\ldots$) | |
| VFMSUBADD{132,213,231}P[DS] | $x_iy_i+(-1)^iz_i$ ($i=0,1,\ldots$) |
これらの命令はAVXで導入されたVEX prefixでエンコードされます。これらの命令は3つのオペランドを取り、結果を入力の1つに代入しますが、それらの組み合わせによって 132, 213, 231 の3つのバリエーションがあります。
浮動小数点数の積はほぼ可換なので、オペランドの順番については
x ← FMA(x, y, z)
z ← FMA(x, y, z)
の2種類あれば良さそうな気がしますが、3種類あるのはNaNのペイロード等を考慮したのかもしれません（要検証）。
ARM (AArch64)
AArch64には最初から（ARMv8-Aから）FMA系の命令があります。
| mneomnic | 表す演算 | 備考 |
|---|---|---|
| FMADD | $xy + z$ | |
| FMSUB | $xy - z$ | |
| FNMADD | $- (xy + z) = - x y - z$ | negated fused multiply-add |
| FNMSUB | $- (xy - z) = - x y + z$ | negated fused multiply-sub |
AArch64のFMA系の命令は4オペランドです。
x86とAArch64では FNMADD, FNMSUB っぽい名前の命令の意味が異なるので要注意です。
libcの実装の検査
各種libcに含まれる fma 関数の実装が正しいかどうかチェックしてみました。いずれもx86_64で動作確認しています。
検証に使ったコード類はGitHubの test-fma に置いておきます。
筆者はIntel C++ Compilerを使える環境を持っていないので、誰かICCを触れる人がいたらそれの動作確認をお願いしたいです（特にWindows上で）。
Linux / glibc
glibcの実装を筆者のテストプログラムで検証したところ、問題は見つかりませんでした。
glibcは実行時（ロード時）にCPUIDで各種命令の有無を検出して実装を切り替える、という機構を持っています。なので、コンパイル時に -mfma を指定し忘れても、実行時のCPUにFMA命令があれば fma 関数の呼び出しの際にはFMA命令が使われます。
glibc的にGCC以外で利用する場合のスタンスがどうなっているのか知りませんが、clangからglibcを使った場合は -mfma を指定しても FP_FAST_FMA が定義されませんでした。それでも fma 関数の呼び出し自体はFMA命令になっています。
Linux / musl
muslの実装を筆者のテストプログラムで検証したところ、問題は見つかりませんでした。
muslは実行時にCPUIDを使ってFMA命令の有無を検出するようなことはなさそうです。常にソフトウェア実装が使われます。
動作確認したわけではありませんが、muslもclangで利用した場合は FP_FAST_FMA が定義されない気がします。
FreeBSD
FreeBSDの実装を筆者のテストプログラムで検証したところ、問題は見つかりませんでした。
FreeBSDのlibmは実行時にCPUIDを使ってFMA命令の有無を検出するようなことはしなさそうです。
FreeBSDの math.h ではコンパイラーによらず常に FP_FAST_FMAF が定義され、 FP_FAST_FMA は定義されません。謎です。
macOS
macOSに付属するlibm（libSystemの一部）の fma の実装を筆者のテストプログラムで検証したところ、問題は見つかりませんでした。
昔のmacOSに付属する（libSystem内の）libmはソースが公開されていましたが、最近のものはソース非公開です。昔のものは以下で参照できます。
公開されている中で最も新しいのはLibm-2026で、これはMac OS X 10.7に付属したもののようです。
筆者の持っているMac (MacBook Pro Late 2013, Catalina, Haswell搭載) では、常にFMA命令が使われるようです。「常に」というのは、古いCPUでの動作を確認しようとしてIntel SDEでIvy Bridgeあたりを指定するとHaswell以降の命令を使おうとして死にます。
Catalina自身は（FMA未搭載の）Ivy Bridge搭載マシンにも対応しているはずなので、FMA命令がない環境向けの実装もどこかに存在するはずです。ということは、CPUID以外の方法でCPUの種類を判別しているか、あるいはシステム起動時にCPUの種類を判別していると考えられます。
（ところでRosetta 2はAVX命令をサポートしないようですが、そうするとFMA命令も使えないのでしょうか）
Windows / UCRT
Universal CRTの fma 関数は、実行時のCPUID等によりFMA命令が利用できるか判断し、利用できる場合はそれを利用するようです。
FMA命令が利用される場合はもちろん問題はありませんが、FMA命令が利用できない場合のソフトウェア実装はbuggyです。Windows上で数値計算するプログラムを書く人は注意してください（そんな人いるわけないか）。libcはC/C++以外の言語からも利用されることがあるので、C/C++以外の人も要注意です。
GCCやClangの場合は -O2 -mfma で fma 関数の呼び出しがFMA命令に置き換えられましたが、Visual C++ではどうでしょうか。VS2019以降のMSVCは /O2 /arch:AVX2 以上を指定すると、FMA命令が使えると判断して fma 関数の呼び出しをFMA命令に置き換えるようです。この際に FP_FAST_FMA は定義されません。
- [C++] std::fma should use _mm_fmadd_ss when AVX2 is enabled - Developer Community
- fmaf function may return wrong results on machines without hardware FMA
Windows / mingw-w64
結論から言うと、mingw-w64の fma 関数もbuggyです。
MinGW / mingw-w64はCランタイムとしてMSVCRTに依存しますが、一部の関数（フォーマット関連、C99で追加された関数など）はMSVCRTに頼らず自前で実装しています。これらの、MinGW / mingw-w64独自の部分のソースは公開されています。
mingw-w64の fma 関数はCPUIDを呼び出すこともなく常にソフトウェア実装を使いますが、その実装が間違っているようです。
一応筆者の方でバグ報告しておきましたが、こういうのは修正パッチを添付して初めて対応される可能性がある系のアレなので、放っておいても修正される可能性は低いでしょう。
ライブラリーの fma 関数がbuggyとは言っても、コンパイラーがFMA命令に置き換える分には関係ないので、GCCやClangで -mfma を指定した場合に fma 関数の直接呼び出しは「正しい」FMAになるはずです。
また、GCC / Clangいずれの場合も FP_FAST_FMA は定義されないようでした。
Emscripten
Emscriptenの場合、 double に対する fma は正しく実装されているようでしたが、 float に対する fmaf は微妙に間違っていました。
数学関数の実装では fesetround 等の浮動小数点環境をいじる関数が使われている場合がありますが、JavaScript / WebAssemblyターゲットの場合はそういう関数を呼び出すことはできません（呼び出したとしても失敗します）。
ソースを見てみたところ、 float に対する fmaf は fesetround(FE_TOWARDZERO) を使っており、 fesetround が使えない環境への配慮が欠けていました。
他の言語
最近の言語のいくつかはFMAを計算する関数を持っています。それらの実装が正しいのか検査してみました。
検証に使ったコード類はGitHubの test-fma に置いておきます。
Java
OracleのJRE (14.0.2)で確かめてみたところ、実行時にFMA命令が利用できる場合はそれを利用するようです。
FMA命令が利用できない場合、ソフトウェア実装を利用するようです。このソフトウェア実装は double に関する方は良さそうですが、 float の方はbuggyです。実行結果を見た感じ、
static float fma(float x, float y, float z) {
return (float)((double)x * (double)y + (double)z);
}
という感じで実装されている可能性が高いです。この実装がダメな理由はすでに説明しました。
OpenJDKのソースも見てみたところ、double に対しては BigDecimal を使うようになっていました。コストはかかりますが、リファレンス実装としてはそれが一番わかりやすいですからね。
一方、float に対する fma は double で積和を計算するようになっていました。これはダメですね。
（筆者はJavaには詳しくないのですが、Oracleが配布してるJavaランタイムってOpenJDKをビルドしたものと同一なんでしたっけ？）
C# / .NET Core
.NET Coreでは Math.FusedMultiplyAdd / MathF.FusedMultiplyAdd としてFMA演算が利用できるようです。
これらは、FMA命令が利用できる場合はそれを利用するようです。
FMA命令が利用できない場合は、libcの fma 関数を呼び出すようです。つまり、Windows上ではbuggyです。
先の記事で書きましたが、C#は16進浮動小数点数リテラルに対応していないので、筆者のテストプログラムを移植するのに手間がかかりました。ファッキン
Rust
f32, f64 型のそれぞれに mul_add という名前でFMA演算が用意されています。
Rustでx86のFMA命令を有効にするには target-feature=+fma という感じの指定をするらしいです。
FMA命令が使えない場合はその環境のlibcの fma 関数を使うようで、Windows上ではMSVC版にせよGNU版にせよbuggyということになります。
ちなみに、Rustも16進浮動小数点数リテラルに対応していないので、ry
Go
math パッケージに float64 用の FMA 関数が用意されています。float32 用のものはなさそうです（すでに説明しましたが、「float64 用のFMAの結果をキャストする」のでは float32 用のFMAにはなりません）。
FMA命令が利用できない場合のFMA関数は独自に実装しているようです。ソースはおそらく src/math/fma.go です。
筆者はGoに詳しくないので、FMA命令を有効化してビルドする方法はよくわかりませんでした。あと、Windows上でGo製のバイナリをIntel SDE上で動かそうとしたらクラッシュしました。
D
std.math にそれっぽい関数がありますが、未実装のようです。未実装なものにどうこう言っても仕方がないのですが、 fma 関数は real だけではなくそれぞれの型について用意してあげないとダメです（2段階で丸めが行われることを回避するため）。
Julia
Juliaは数学関数の実装を自前で持っているらしいです（openlibm）。おそらく正しく実装されていることが期待できます。
後で確認します。
まとめ
この記事では
- FMAとは何か
- FMAを使うための書き方（主に標準Cに則った方法で）
- 一部の環境でFMA関数が間違った値を返す話
を書きました。それでは、皆さんたのしいFMAライフを！
-
正確には、絶対値が $2^{1024}-2^{971}$ より大きい数。 ↩