はじめに
低速な演算であるdivとsqrtを高速化するrcpとrsqrtについて説明します．
逆数（RCP）
逆数(reciprocal)である 1/x を高速に計算します．除算をするよりも，rcp結果を乗算したほうが高速です．
ただし，精度は11から12ビット程度です（正確には，1.5x2^-12）．
double精度のrcp命令はありません．
AVX-512ではrcp14_psと14ビット精度のrcpが使え，通常のrcpとレイテンシも同じです．
なお，11ビットrcp命令はAVX-512にはありません．
_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)
__m256 _mm256_rcp_ps (__m256 a)
asm: vrcpps ymm, ymm
動作
_mm256_rcp_ps
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ビットの浮動小数点演算結果とほぼ同等精度になります．
rcp_{t+1}(x) = rcp_{t}(x)*(2-rcp_{t}(x)*x)
実装を以下に示します．
式を展開すれば，2倍は加算で表現できるため，最初の実装ができます．
元の式を使えばFMAで計算できるためより高速ですが，定数値２をロードする必要があり，2の定数を保持したレジスタを使いまわすことができない場合は，ロードのレイテンシがボトルネックとなります．
2回以上繰り返せばわずかに高精度になりますが，浮動小数点の桁を踏まえるとそれほど高精度化しません．
//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にはありません．
_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)
__m256 _mm256_rsqrt_ps (__m256 a)
asm: vrsqrtps ymm, ymm
動作
_mm256_rsqrt_ps
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型の逆数平
方根を計算します．パフォーマンスは，rcpと同じです．
double型の命令はありません．
平方根が必要な場合は，逆数を取るのではなく，乗算をしてください．
平方根の逆数は，自信との乗算で戻ります（$\sqrt{a} = a\frac{1}{\sqrt{a}}$）．
11ビット精度しかいらない場合は，rcp(rsqrt)がもっとも速く動きます．
22ビット精度で必要なら，ニュートンラフソン法で高精度化します．
下記の漸化式で更新します．
N^0.5 = X_1*N = 1/2*N*X_0*(3-N*X_0^2)
実装は下記を使用してください．
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については後に説明します．