AVX/AVX2によるFMA

# はじめに
Fused-multily-add (FMA)の説明をします．

## 融合積和（FMA）
Fused-multily-add (FMA)は，積和演算を1命令で行います．
乗算と加算を同時に実行するため，演算量を半分にするだけでなく，浮動小数点の演算誤差も抑制できます．
FMAは，IntelのCPUだとHaswellから使えます．AMDのZenシリーズは，すべてで使えます．

FMAは，MAD, MADD (multiply-add) やMAC, MACC (multiply-accumulate)とも呼ばれます．
なお，FMAを漢字の融合積和と呼ぶことはまれです．

FMA命令には，下記**6種類**ありますが，演算コストは全て同じです．
AVXで使えるFMAは，すべて浮動小数点命令です．
整数FMAは，IcelakeからAVX512IFMAでサポートされ，52ビット精度の整数FMAが計算可能です．

表：Intel/AMDのFMA

|Instruction|Expression   |
|:---------:|:-----------:|
|fmadd      | a*b+c       |
|fmsub      | a*b-c       |
|fnmadd     | -a*b+c      |
|fnmsub     | -a*b-c      |
|fmaddsub   | a*b+c,a*b-c |
|fmsubadd   | a*b-c,a*b+c |

FMAは通常d=ab+cを表現するために，4つの変数が必要となります．
例えば，下記のような形を取らなければなりません．

```cpp
//d=ab+c
fma ymm ymm ymm ymm
```

この形を持っていたのがRyzenより前のAMDのCPUで，2012年のBulldozerから実装されています．
しかしながら，この形式だと命令が長くなるため，Intelは3オペランドのFMAを2013年のHaswellで設計し，どれか一つの引数を破壊して出力とすることで3オペランドのFMAを実現しました．
そのため，同じFMA命令であっても下記のように，どこを破壊するかによって3通りの命令があります．
つまり，以下のように演算します．

```cpp
a = a*b+c;
b = a*b+c;
c = a*b+c;
```

対応するアセンブリ命令は以下のようになります．

```cpp
vfmadd132ps ymm, ymm, ymm
vfmadd213ps ymm, ymm, ymm
vfmadd231ps ymm, ymm, ymm
```

FMAには，計算式で6種類あり，破壊するレジスタを選択するのに3種類あるため，アセンブリとしては18種類あります．
intrinsicsで書く限りレジスタ管理をする必要がないため，どの命令が呼ばれたか意識する必要はありませんが，どれかのレジスタが破壊されるということは覚えておいたほうが良いでしょう．

なお，AMDもPiledriver，Kabini以降から3オペランドFMAに対応しています．
また，現在のAMDのCPUでは，4オペランドのFMAは利用が非推奨であったり，Zen3ではすでに使用することすらできません．
使用する場合も，intrinsicsから呼び出すことはできず，アセンブリから呼び出す形でしか使用できません．

参考までに，ARMのFMA命令も示します．
Intel/AMDとARMのfnmadd/subの命令は，マイナスの符号のかかる場所が違うので注意してください．
* [fnmadd](https://developer.arm.com/documentation/100076/0100/a64-instruction-set-reference/a64-floating-point-instructions/fnmadd?lang=en)
* [fnmsub](https://developer.arm.com/documentation/100076/0100/a64-instruction-set-reference/a64-floating-point-instructions/fnmsub?lang=en)

表：ARMのFMA

|Instruction|Expression   |
|:---------:|:-----------:|
|fmadd      | a*b+c       |
|fmsub      | a*b-c       |
|fnmadd     | -(a*b+c)    |
|fnmsub     | -(a*b-c)    |

## _mm256_fmadd_ps|pd (FMA)
```cpp
__m256  _mm256_fmadd_ps (__m256 a, __m256 b, __m256 c)
__m256d _mm256_fmadd_pd (__m256d a, __m256d b, __m256d c)
asm: vfmadd132ps ymm, ymm, ymm
     vfmadd213ps ymm, ymm, ymm
     vfmadd231ps ymm, ymm, ymm
asm: vfmadd132pd ymm, ymm, ymm
     vfmadd213pd ymm, ymm, ymm
     vfmadd231pd ymm, ymm, ymm
```
**動作**  
作成中

**CPI/Uops**  

| Architecture | Latency | Throughput |Uops|
|--------------|:-------:|:----------:|:--:|
| Alderlake    |    4    |    0.5     |  - |
| Icelake      |    4    |    0.5     |  1 |
| Skylake      |    4    |    0.5     |  1 |
| Broadwell    |    5    |    0.5     |  1 |
| Haswell      |    5    |    0.5     |  1 |
| Zen3         |    4    |    0.5     |  1 |
| Zen2         |    5    |    0.5     |  1 |
| Zen          |    5    |     1      |  2 |

* _mm_fmadd_psなどの128ビットに対する命令も，呼んでいるアセンブリは同じのためコストは同一です．
* FMAの全18命令，コストは同じです．

**説明**  


参考までに，加算，乗算，FMAのパフォーマンスを並べます．
Alderlake以外，どの命令であってもFMA，乗算と加算はFMAで混ぜたほうが高速です．

また，AVX-512のFMAは，サーバー系とクライアント向けで性能が違います．
Icelake, Tigerlake, Rocketlakeはクライアント向けですがAVX-512が使えます．しかし，FMAユニットの数からパフォーマンスがサーバー向けの半分です．

ADD 

| Architecture | Latency | Throughput |Uops|
|--------------|:-------:|:----------:|:--:|
| Alderlake    |    2    |    0.5     |  - |
| Icelake      |    4    |    0.5     |  1 |
| Skylake      |    4    |    0.5     |  1 |
| Broadwell    |    3    |     1      |  1 |
| Haswell      |    3    |     1      |  1 |
| Ivy Bridge   |    3    |     1      |  1 |
| Sandy Bridge |    3    |     1      |  1 |
| Zen3         |    3    |    0.5     |  1 |
| Zen2         |    3    |    0.5     |  1 |
| Zen          |    3    |     1      |  2 |

MUL

| Architecture | Latency | Throughput |Uops|
|--------------|:-------:|:----------:|:--:|
| Alderlake    |    4    |    0.5     |  - |
| Icelake      |    4    |    0.5     |  1 |
| Skylake      |    4    |    0.5     |  1 |
| Broadwell    |    3    |    0.5     |  1 |
| Haswell      |    5    |    0.5     |  1 |
| Ivy Bridge   |    5    |     1      |  1 |
| Sandy Bridge |    5    |     1      |  1 |
| Zen3         |    3    |    0.5     |  1 |
| Zen2         |    3    |    0.5     |  1 |
| Zen          |   3(4)  |     1      |  2 |

AVX512 ADD/MUL/FMA

| Architecture | Latency | Throughput |Uops|
|--------------|:-------:|:----------:|:--:|
| Icelake      |    4    |     1      |  1 |
| Icelake-SP   |    4    |    0.5     |  1 |
| Cascade Lake |    4    |    0.5     |  1 |
| Skylake-X    |    4    |    0.5     |  1 |