可能な代替案は、代わりに計算するsqrt((x*x+y*y)/2ことです。これは、すべての可能なベクトルの大きさを0..255の範囲にスケーリングします。
2つの（高速）アルゴリズムはほぼ完璧な結果をもたらすようです。1つはCordicを使用し、もう1つは最大の内積を使用します。
void cordic_it(uint16 &x, uint16 &y, int n) {
auto X = x + y >> n; // vsraq_n_u16(x, y, n) in arm neon
y = abs(y - x >> n); // vabdq_u16(y, x >> n) in arm neon
}
uint16_t scaled_magnitude_cordic(uint8_t x, uint8_t y) {
const int kRound = 1;
if (x < y) std::swap(x,y);
// multiply by factor of 256/sqrt(2) == 181.02
// then reduce by the gain of the cordic iterations of 1.16
// - with prescaling we also ensure, that the cordic iterations
// do not lose too much significant bits when shifting right
uint16_t X = x * 156, Y = y * 156;
// exactly 4 iterations. 3 is too little, 5 causes too much noise
for (int j = 1; j <= 4; j++) cordic_it(X,Y,j);
return (X+kRound) >> 8;
}
kRoundを変更することで、結果を調整できます。
Histogram of real - approx: -1 0 1
kRound == 0 -> smaller code 1 46617 18918
kRound == 1 -> approx >= real 0 46378 19158
kRound == -73 -> balanced error 3695 58301 3540
を選択するkRound == 1と、次の方法ですべての結果を修正できます。
uint8_t fix_if_larger_by_one(uint8_t sqrt, uint8_t x, uint8_t y) {
auto P = (x*x + y*y) / 2;
auto Q = sqrt*sqrt;
return sqrt - (P < Q);
}
x a + y bの内積をいくつかの角度で近似することにより、平方根を計算することもできます。従来のアプローチでは、単一の角度を使用しますa = 1, b = 1/2。
5つの固有の角度で、およそ[0 10 20 30 40]またはの角度に対して[5 15 25 35 45]、いずれかの係数のセットが考えられます。どちらも、最大で1単位ずれたほぼ完全な結果を生成します。
1) [181 0], [178 31], [170 62], [157 91], [139 116]
2) [180 18], [175 46], [164 76], [148 104], [128 128]
オプション1には9つの自明でない係数があります（ただし、62 == 31 * 2）。オプション2には8つの自明でない係数があり、次の実装に役立ちます。
int approx(uint8_t x, uint8_t y) {
if (x < y) std::swap(x,y); // sort so that x >= y
auto a4 = (x + y) / 2; // vhaddq_u8(x,y) on Arm Neon
auto a0 = (x * 180 + y * 18) >> 8;
auto a1 = (x * 175 + y * 46) >> 8;
auto a2 = (x * 164 + y * 76) >> 8;
auto a3 = (x * 148 + y * 104) >> 8;
return max_of_five_elements(a0,a1,a2,a3,a4);
}
_mm_maddubs_epi16このほぼ偶数の係数のセットは、インストリンシクスを使用したSSSE3命令セットに非常にうまく変換され_mm_max_epu16ます。各ドット積ですがa1、インターリーブされたx、yおよびインターリーブされた係数から1つの命令で簡単に計算できます。_mm_packus_epi16当然、レイテンシーに対抗し、uint8_t入力からの計算、並べ替え、または平均化を無駄にしないために、16個の隣接する近似を同時に計算する方が理にかなっています。
auto a0 = _mm_maddubs_epi16(xy, coeffs0); // coeffs0 = 90 9 90 9 ...
auto a1 = _mm_maddubs_epi16(xy, coeffs1); // coeffs1 = 87 23 87 23 ...
auto a2 = _mm_maddubs_epi16(xy, coeffs2); // coeffs2 = 82 38 82 38 ...
auto a3 = _mm_maddubs_epi16(xy, coeffs3); // coeffs3 = 74 52 74 52 ...
auto a4 = _mm_maddubs_epi16(xy, coeffs4); // coeffs4 = 64 64 64 64 ...
a1 = _mm_add_epi16(a1, x_per_2); // LSB of the coefficient 87.5
// take the maximum, shift right by 7 and pack to uint8_t
a0 = _mm_max_epu16(a0, a1);
a0 = _mm_max_epu16(a0, a2);
a0 = _mm_max_epu16(a0, a3);
a0 = _mm_max_epu16(a0, a4);
a0 = _mm_srli_epi16(a0, 7);
a0 = _mm_packus_epi16(a0, a0);
わずか8つの係数を使用することは、ARM Neonの実装にも適しています。これにより、16ビット×16ビットのスカラー乗算を使用して、すべての係数を単一の全幅レジスタに格納できます。
完璧な結果を得るには、ドット積アルゴリズムを他の方向に補正する必要があります。これは、次のリファレンス実装の1つ下の要素である値を与える可能性があるためですfloor(sqrt((x*x+y*y)/2)。
uint8_t fix_if_smaller_by_one(uint8_t sqrt, uint8_t x, uint8_t y) {
auto P = (x*x + y*y) / 2;
auto Q = (sqrt+1)*(sqrt+1);
return sqrt + (Q <= P);
}
他の近似アルゴリズムは通常、除算またはスケーリングのいずれかを使用します。これは、レーンごとの可変シフトがないため、AVX2より前のIntelではベクトル化が困難です。