AVX/AVX2によるRCPとRSQRT

# はじめに
低速な演算であるdivとsqrtを高速化するrcpとrsqrtについて説明します．

## 逆数（RCP）
![v1_32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/67908/5d0e8b7d-609e-9235-09f3-016bc75b9d2f.png)

逆数(reciprocal)である 1/x を高速に計算します．除算をするよりも，rcp結果を乗算したほうが高速です．
ただし，精度は11から12ビット程度です（正確には，1.5x2^-12）．
double精度のrcp命令はありません．

AVX-512ではrcp14_psと14ビット精度のrcpが使え，通常のrcpとレイテンシも同じです．
なお，11ビットrcp命令はAVX-512にはありません．

```cpp
_mm_rcp14_ps(__m128 src) //AVX-512
_mm256_rcp14_ps(__m256 src) //AVX-512
_mm512_rcp14_ps(__m512 src) //AVX-512
```

28ビット精度のrcpもありますが，AVX512-ERで，これはXeon Phi専用です．
現在の通常のCPUではサポートしていません．

また，数値計算精度は，ニュートンラフソン（Newton-Raphson）法により高精度化できます．
これは，4回の乗算，加減算か，2回のFMA演算の追加演算で実現できます．

RCPの話
https://math-koshimizu.hatenablog.jp/entry/2017/07/28/083000

#### _mm256_rcp_ps (AVX)
```cpp
__m256 _mm256_rcp_ps (__m256 a)
asm: vrcpps ymm, ymm
```
**動作**  
_mm256_rcp_ps
![v1_32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/67908/27405d70-39b4-8ff7-9c54-a6d2faa51e1e.png)

**CPI/Uops**  

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

**説明**  
11~12ビット精度で逆数を高速に求めます．
除算結果はa/b = a*(1/b)で求まります．

精度を高速化するときは，下記，ニュートンラフソン法を使います．
1度のイタレーションで22ビット精度と，23ビットの浮動小数点演算結果とほぼ同等精度になります．

```math
rcp_{t+1}(x) = rcp_{t}(x)*(2-rcp_{t}(x)*x)
```

実装を以下に示します．
式を展開すれば，2倍は加算で表現できるため，最初の実装ができます．
元の式を使えばFMAで計算できるためより高速ですが，定数値２をロードする必要があり，2の定数を保持したレジスタを使いまわすことができない場合は，ロードのレイテンシがボトルネックとなります．
2回以上繰り返せばわずかに高精度になりますが，浮動小数点の桁を踏まえるとそれほど高精度化しません．

```cpp
//rcp with newton-raphson 1-iteration
__m256 _mm256_rcpnr_ps(__m256 x)
{
	__m256 res = _mm256_rcp_ps(x);
	//rcp*(2-rcp*x)->(rcp+rcp)-rcp*rcp*x
	return res = _mm256_sub_ps(_mm256_add_ps(res, res), _mm256_mul_ps(x, _mm256_mul_ps(res, res)));
}

//rcp with newton-raphson 1-iteration (FMA ver) requided set2
__m256 _mm256_rcpnr_fma_ps(__m256 x)
{
	__m256 rcp = _mm256_rcp_ps(x);
	//rcp*(2-rcp*x)
	return _mm256_mul_ps(rcp, _mm256_fnmadd_ps(x, rcp, _mm256_set1_ps(2.f)));
}

//rcp with newton-raphson 2-iteration
__m256 _mm256_rcpnr2_ps(__m256 x)
{
	__m256 rcp = _mm256_rcp_ps(x);
	rcp = _mm256_sub_ps(_mm256_add_ps(rcp, rcp), _mm256_mul_ps(x, _mm256_mul_ps(rcp, rcp)));
	return rcp = _mm256_sub_ps(_mm256_add_ps(rcp, rcp), _mm256_mul_ps(x, _mm256_mul_ps(rcp, rcp)));
}
```

### 逆数平行根（RSQRT）
除算と同様に，１１ビット精度で平方根の逆数を高速に計算します．
doubleはありません．
こちらも，ニュートンラフソン法で高精度化可能です．

AVX-512では，14ビット精度のrsqrtが使え，通常のrsqrtとレイテンシも同じです．
なお，11ビットrsqrt命令はAVX-512にはありません．

```cpp
_mm_rsqrt14_ps(__m128 src) //AVX-512
_mm256_rsqrt14_ps(__m256 src) //AVX-512
_mm512_rsqrt14_ps(__m512 src) //AVX-512
```

#### _mm256_rsqrt_ps (AVX)
```cpp
__m256 _mm256_rsqrt_ps (__m256 a)
asm: vrsqrtps ymm, ymm
```

**動作**  
_mm256_rsqrt_ps
![v1_32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/67908/27405d70-39b4-8ff7-9c54-a6d2faa51e1e.png)

**CPI/Uops**  

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

**説明**  
float型の逆数平![v1_32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/67908/799ed943-2821-674c-b7bc-d919707838b0.png)
方根を計算します．パフォーマンスは，rcpと同じです．
double型の命令はありません．

平方根が必要な場合は，逆数を取るのではなく，乗算をしてください．
平方根の逆数は，自信との乗算で戻ります（$\sqrt{a} = a\frac{1}{\sqrt{a}}$）．
11ビット精度しかいらない場合は，rcp(rsqrt)がもっとも速く動きます．

22ビット精度で必要なら，ニュートンラフソン法で高精度化します．
下記の漸化式で更新します．

```math
N^0.5 = X_1*N = 1/2*N*X_0*(3-N*X_0^2)
```

実装は下記を使用してください．

```cpp
inline __m128 _mm_rsqrtnr_ps(__m128 x) 
{
	__m128 three = _mm_set1_ps(3.0f), half = _mm_set1_ps(0.5f);
	__m128 rsqrt = _mm_rsqrt_ps(x);
	__m128 sqrt = _mm_mul_ps(x, rsqrt);
	return rsqrt = _mm_mul_ps(_mm_mul_ps(half, rsqrt), _mm_fnmadd_ps(sqrt, rsqrt, three));
}

inline __m256 _mm256_rsqrtnr_ps(__m256 x)
{
	__m256 three = _mm256_set1_ps(3.0f), half = _mm256_set1_ps(0.5f);
	__m256 rsqrt = _mm256_rsqrt_ps(x);
	__m256 sqrt = _mm256_mul_ps(x, rsqrt);
	return rsqrt = _mm256_mul_ps(_mm256_mul_ps(half, rsqrt), _mm256_fnmadd_ps(sqrt, rsqrt, three));
}
```

また，SVMLの`_mm256_invsqrt_ps`は，ニュートンラフソン法で高精度化した逆数平方根です．
SVMLについては後に説明します．