読者です 読者をやめる 読者になる 読者になる

おぺんcv

画像処理エンジニアのブログ

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 から減算して、最終的な結果をデスティネーションベクタの要素に返します。

最初この説明を見たときはなんのこっちゃと思いましたが
ニュートン法で逆数の近似を求める際に使うようです

ある数 d とその逆数の現在の近似値 x_nがあるとき、ニュートン法によって
 x_{n+1} = x_n(2 - dx_n)
という反復を繰り返すことで近似値が 1/dに収束します
この (2 - dx_n)を計算してくれるのが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のみ)

vdiv[q]_f32(va, vb)

32bit浮動小数ベクタvaとvbの除算を行います

サンプル

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の使い方 初期化編