FMA (fused multiply-add) の話

この記事では**融合積和** (fused multiply-add; 以下もっぱら **FMA**) について扱います。

この記事は元々「その辺で提供されている `fma` 関数の実装が正しいかチェックしてみた」記事だったのがだんだん守備範囲を広げていったものなので、FMAを単に使うだけではなく、**実装する**側の視点が多く含まれています。

筆者が実装したHaskell向けのFMA関数はこちらにあります：

* [haskell-floating-point/FMA.hs at master · minoki/haskell-floating-point](https://github.com/minoki/haskell-floating-point/blob/master/fp-ieee/src/Numeric/Floating/IEEE/Internal/FMA.hs)

# 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つの関数が導入されました。

```c
#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`（丸めモードが「負の無限大方向」以外の場合）

という規則によって 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命令の実行に置き換えられ、関数呼び出しは発生しません。

```c
#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命令に変換されました。

```c
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` に対応していないのは興味深いです（セマンティクスが違うとかそういう話なのでしょうが）。

* [c++ - Automatically generate FMA instructions in MSVC - Stack Overflow](https://stackoverflow.com/questions/34265982/automatically-generate-fma-instructions-in-msvc)

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` です。

```c
#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` を指定するべきなのです。

* 参考：[Other Builtins (Using the GNU Compiler Collection (GCC))](https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html#Other-Builtins)

# FMAの実装の落とし穴

Cやその他言語での `fma` 関数を実装する側のことを考えてみましょう。

まず、ターゲットとなるアーキテクチャが必ずFMA命令を持っている場合（例：AArch64）はその命令を直接呼び出せばOKです。ハードウェア実装されたFMA命令はまず間違いなく正しく実装されています（さもなくば「CPUのバグ」として報告されます）。

一方、ターゲットがFMA命令を持っているとは限らない場合（例：x86）は、何らかの方法でFMAをソフトウェア実装する必要があります。標準ライブラリーを実装する人はFMAのことを隅々まで把握しているとは限らないため、そういう人が書いた `fma` 関数の実装にはバグが潜んでいる可能性があります。

標準ライブラリーの実装者が正しく `fma` 関数を実装できるように、あるいは、利用者が自分の使う `fma` 関数のバグを早期発見できるように、ここではFMAの典型的な「間違い方」をいくつか挙げてみます。

## オーバーフロー

浮動小数点数はかなり（指数部が）広い範囲の数を表せますが、それでも表せる範囲には限りがあります。
例えば、よく使われる倍精度（`double` 型）では絶対値が $2^{-1074}$ 未満の数を表現することはできませんし、絶対値が $2^{1024}$ 以上[^double-max]の数は表現できません。

[^double-max]: 正確には、絶対値が $2^{1024}-2^{971}$ より大きい数。

さて、FMAも浮動小数点演算ですから、結果（の絶対値）があまりにも大きいとオーバーフローを起こして無限大を返します。

しかし、浮動小数点数での乗算 `x * y` がオーバーフローしても FMA `fma(x, y, z)` がオーバーフローするとは限りません。

具体的には以下の場合です。

1. `z` が無限大の場合
2. `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` を使って

```c
(float)((double)x * (double)y + (double)z)
```

と計算した方が精度が上がりそうな気がします。何だったら、FMAの実装だってこれでいいんじゃないのという気がしてしまいます。

ですが、これはダメです。通常の積和ではなく（`double` に対する）FMAを使ってもダメです。

```c
// ダメな実装の例
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メーカーがこういう間違いをやらかしたら[最悪リコール](https://en.wikipedia.org/wiki/Pentium_FDIV_bug)ですよ（最近はマイクロコードの修正で済むのかもしれませんが）。

ちなみに、 `float` に対するFMAを実装する際に `double` が全く活用できないかというとそういうわけでもなくて、積の部分を `double` で計算することによってアルゴリズムを簡単にすることができます。積を計算した後の和の部分はどうにか工夫する必要があります。

## 高精度な型を使った手抜き

`float` の場合と同じ理屈で、多少高精度な `long double` 型があったからといって `double` に対するFMAを

```c
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](https://github.com/minoki/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](https://opensource.apple.com/source/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上で数値計算するプログラムを書く人は注意してください<small>（そんな人いるわけないか）</small>。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](https://developercommunity.visualstudio.com/content/problem/89929/c-stdfma-should-use-mm-fmadd-ss-when-avx2-is-enabl.html)
* [fmaf function may return wrong results on machines without hardware FMA](https://social.msdn.microsoft.com/Forums/en-US/cef029f6-315c-4c97-af2a-3f255b28eac9/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を呼び出すこともなく常にソフトウェア実装を使いますが、その実装が間違っているようです。

一応筆者の方でバグ報告しておきましたが、こういうのは修正パッチを添付して初めて対応される可能性がある系のアレなので、放っておいても修正される可能性は低いでしょう。

* [MinGW-w64 - for 32 and 64 bit Windows / Bugs / #848 fma returns an incorrect result](https://sourceforge.net/p/mingw-w64/bugs/848/)

ライブラリーの `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` が使えない環境への配慮が欠けていました。

* [emscripten/fmaf.c at master · emscripten-core/emscripten](https://github.com/emscripten-core/emscripten/blob/master/system/lib/libc/musl/src/math/fmaf.c)

# 他の言語

最近の言語のいくつかはFMAを計算する関数を持っています。それらの実装が正しいのか検査してみました。

検証に使ったコード類はGitHubの [test-fma](https://github.com/minoki/test-fma) に置いておきます。

## Java

OracleのJRE (14.0.2)で確かめてみたところ、実行時にFMA命令が利用できる場合はそれを利用するようです。

FMA命令が利用できない場合、ソフトウェア実装を利用するようです。このソフトウェア実装は `double` に関する方は良さそうですが、 `float` の方はbuggyです。実行結果を見た感じ、

```java
static float fma(float x, float y, float z) {
    return (float)((double)x * (double)y + (double)z);
}
```

という感じで実装されている可能性が高いです。この実装がダメな理由はすでに説明しました。

OpenJDKのソースも見てみたところ、`double` に対しては `BigDecimal` を使うようになっていました。コストはかかりますが、リファレンス実装としてはそれが一番わかりやすいですからね。

* [java/lang/Math.java -- double fma (double, double, double)](https://hg.openjdk.java.net/jdk/jdk/file/d72a0f0dc013/src/java.base/share/classes/java/lang/Math.java#l1757)

一方、`float` に対する `fma` は `double` で積和を計算するようになっていました。これはダメですね。

* [java/lang/Math.java -- float fma(float, float, float)](https://hg.openjdk.java.net/jdk/jdk/file/d72a0f0dc013/src/java.base/share/classes/java/lang/Math.java#l1858)

（筆者はJavaには詳しくないのですが、Oracleが配布してるJavaランタイムってOpenJDKをビルドしたものと同一なんでしたっけ？）

## C# / .NET Core

.NET Coreでは `Math.FusedMultiplyAdd` / `MathF.FusedMultiplyAdd` としてFMA演算が利用できるようです。

これらは、FMA命令が利用できる場合はそれを利用するようです。

FMA命令が利用できない場合は、libcの `fma` 関数を呼び出すようです。つまり、**Windows上ではbuggyです**。

[先の記事](https://qiita.com/mod_poppo/items/3fa4cdc35f9bfb352ad5)で書きましたが、C\#は16進浮動小数点数リテラルに対応していないので、筆者のテストプログラムを移植するのに手間がかかりました。<small>ファッキン</small>

## Rust

`f32`, `f64` 型のそれぞれに `mul_add` という名前でFMA演算が用意されています。

* [f32::mul_add](https://doc.rust-lang.org/std/primitive.f32.html#method.mul_add)
* [f64::mul_add](https://doc.rust-lang.org/std/primitive.f64.html#method.mul_add)

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にはなりません）。

* [func FMA - math - The Go Programming Language](https://golang.org/pkg/math/#FMA)

FMA命令が利用できない場合のFMA関数は独自に実装しているようです。ソースはおそらく [src/math/fma.go](https://github.com/golang/go/blob/master/src/math/fma.go) です。

筆者はGoに詳しくないので、FMA命令を有効化してビルドする方法はよくわかりませんでした。あと、Windows上でGo製のバイナリをIntel SDE上で動かそうとしたらクラッシュしました。

## D

std.math にそれっぽい関数がありますが、未実装のようです。未実装なものにどうこう言っても仕方がないのですが、 `fma` 関数は `real` だけではなくそれぞれの型について用意してあげないとダメです（2段階で丸めが行われることを回避するため）。

* [std.math - D Programming Language](https://dlang.org/phobos/std_math.html#.fma)

## Julia

Juliaは数学関数の実装を自前で持っているらしいです（openlibm）。おそらく正しく実装されていることが期待できます。

後で確認します。

# まとめ

この記事では

* FMAとは何か
* FMAを使うための書き方（主に標準Cに則った方法で）
* 一部の環境でFMA関数が間違った値を返す話

を書きました。それでは、皆さんたのしいFMAライフを！