ARM NEONの使い方 除算編
四則演算ラスト!
今回は除算編です
アーキテクチャによる違い
32bit ARM
32bit ARMには除算を直接行うSIMD命令がないため
a/bという除算を行う場合、まず逆数を計算するSIMD命令によってbの逆数を計算し
これをaと掛け算することによって除算を表現します
64bit ARM
64bit ARMには除算を行うSIMD命令が追加されています
逆数の計算
vrecpe[q]_f32(va)
32bit浮動小数ベクタvaの各要素の逆数の近似値を計算します
近似値なのでCPU(IEEE754)よりも精度が低いです
サンプル
32bit浮動小数ベクタva, vbの除算をva*(vbの逆数)という形でやってみます
vaとvbの各レーンを同じ値に設定しているので、除算後の値が1になるのが望ましい結果です
#include <stdio.h> #include <stdint.h> #include <arm_neon.h> static void print_vector(float32x4_t value) { float lane[4]; vst1q_f32(lane, value); for (int i = 0; i < 4; i++) printf("lane[%d]: %f\n", i, lane[i]); } int main() { float a[4] = { 1, 2, 3, 4 }; float b[4] = { 1, 2, 3, 4 }; float32x4_t va = vld1q_f32(a); float32x4_t vb = vld1q_f32(b); float32x4_t vrecipb = vrecpeq_f32(vb); // approximate reciprocal of vb float32x4_t vc = vmulq_f32(va, vrecipb); printf("va:\n"); print_vector(va); printf("\nvb:\n"); print_vector(vb); printf("\n1/vb:\n"); print_vector(vrecipb); printf("\nva * (1/vb):\n"); print_vector(vc); return 0; }
実行結果
va: lane[0]: 1.000000 lane[1]: 2.000000 lane[2]: 3.000000 lane[3]: 4.000000 vb: lane[0]: 1.000000 lane[1]: 2.000000 lane[2]: 3.000000 lane[3]: 4.000000 1/vb: lane[0]: 0.998047 lane[1]: 0.499023 lane[2]: 0.333008 lane[3]: 0.249512 va * (1/vb): lane[0]: 0.998047 lane[1]: 0.998047 lane[2]: 0.999023 lane[3]: 0.998047
近似値なので結構精度低いなーという印象
まあ用途次第でしょうけど
ベクタ逆数のステップ
よくわからないタイトルですが
先ほどの逆数の精度を良くする方法です
vrecps[q]_f32(va)
ARM Information Centerの説明が以下です
VRECPS (ベクタ逆数のステップ)は、あるベクタの要素を対応する別のベクタの要素で乗算し、その結果を 2 から減算して、最終的な結果をデスティネーションベクタの要素に返します。
最初この説明を見たときはなんのこっちゃと思いましたが
ニュートン法で逆数の近似を求める際に使うようです
ある数とその逆数の現在の近似値があるとき、ニュートン法によって
という反復を繰り返すことで近似値がに収束します
このを計算してくれるのがvrecpsというわけ
サンプル
先ほどの逆数のサンプルにニュートン法の反復計算を追加してみました
#include <stdio.h> #include <stdint.h> #include <arm_neon.h> static void print_vector(float32x4_t value) { float lane[4]; vst1q_f32(lane, value); for (int i = 0; i < 4; i++) printf("lane[%d]: %f\n", i, lane[i]); } int main() { float a[4] = { 1, 2, 3, 4 }; float b[4] = { 1, 2, 3, 4 }; float32x4_t va = vld1q_f32(a); float32x4_t vb = vld1q_f32(b); float32x4_t vrecipb = vrecpeq_f32(vb); // approximate reciprocal of vb float32x4_t vc = vmulq_f32(va, vrecipb); printf("va:\n"); print_vector(va); printf("\nvb:\n"); print_vector(vb); printf("\n1/vb:\n"); print_vector(vrecipb); printf("\nva * (1/vb):\n"); print_vector(vc); // Newton-Raphson iteration int steps = 3; for (int i = 0; i < steps; i++) { printf("\nNewton step: %d\n", i + 1); float32x4_t vtmp = vrecpsq_f32(vb, vrecipb); vrecipb = vmulq_f32(vrecipb, vtmp); vc = vmulq_f32(va, vrecipb); printf("\n1/vb:\n"); print_vector(vrecipb); printf("\nva * (1/vb):\n"); print_vector(vc); } return 0; }
実行結果
va: lane[0]: 1.000000 lane[1]: 2.000000 lane[2]: 3.000000 lane[3]: 4.000000 vb: lane[0]: 1.000000 lane[1]: 2.000000 lane[2]: 3.000000 lane[3]: 4.000000 1/vb: lane[0]: 0.998047 lane[1]: 0.499023 lane[2]: 0.333008 lane[3]: 0.249512 va * (1/vb): lane[0]: 0.998047 lane[1]: 0.998047 lane[2]: 0.999023 lane[3]: 0.998047 Newton step: 1 1/vb: lane[0]: 0.999996 lane[1]: 0.499998 lane[2]: 0.333333 lane[3]: 0.249999 va * (1/vb): lane[0]: 0.999996 lane[1]: 0.999996 lane[2]: 0.999999 lane[3]: 0.999996 Newton step: 2 1/vb: lane[0]: 1.000000 lane[1]: 0.500000 lane[2]: 0.333333 lane[3]: 0.250000 va * (1/vb): lane[0]: 1.000000 lane[1]: 1.000000 lane[2]: 1.000000 lane[3]: 1.000000 Newton step: 3 1/vb: lane[0]: 1.000000 lane[1]: 0.500000 lane[2]: 0.333333 lane[3]: 0.250000 va * (1/vb): lane[0]: 1.000000 lane[1]: 1.000000 lane[2]: 1.000000 lane[3]: 1.000000
(定量的に評価したわけではありませんが)
ステップを進めると精度が良くなっていくのがわかります
除算 (64bit ARMのみ)
サンプル
vdivq_f32を使用して32bit浮動小数ベクタva, vbの除算をやってみます
コンパイルする際は64bit用のコンパイラ(aarch64-linux-gnu-g++など)を使ってね
#include <stdio.h> #include <stdint.h> #include <arm_neon.h> static void print_vector(float32x4_t value) { float lane[4]; vst1q_f32(lane, value); for (int i = 0; i < 4; i++) printf("lane[%d]: %f\n", i, lane[i]); } int main() { float a[4] = { 1, 2, 3, 4 }; float b[4] = { 1, 2, 3, 4 }; float32x4_t va = vld1q_f32(a); float32x4_t vb = vld1q_f32(b); float32x4_t vc = vdivq_f32(va, vb); printf("va:\n"); print_vector(va); printf("\nvb:\n"); print_vector(vb); printf("\nva / vb:\n"); print_vector(vc); return 0; }
実行結果
va: lane[0]: 1.000000 lane[1]: 2.000000 lane[2]: 3.000000 lane[3]: 4.000000 vb: lane[0]: 1.000000 lane[1]: 2.000000 lane[2]: 3.000000 lane[3]: 4.000000 va / vb: lane[0]: 1.000000 lane[1]: 1.000000 lane[2]: 1.000000 lane[3]: 1.000000
そのまま除算ができるとありがたいですね
また64bit ARM(ARMv8)では浮動小数のSIMD演算がIEEE754準拠になったため
結果もCPUと変わらないはずです
参考
NEONの32bitと64bitの違いが解説されてます
今回の記事を書くにあたりとても参考になりました
community.arm.com
次回
次回は初期化編です
ARM NEONの使い方 初期化編