はじめに

この記事はひとりNEONアドベントカレンダー2020 7日目の記事です
昨日は比較命令系を紹介したが、本日は割り算命令を紹介
と言いつつ、このアドベントカレンダー始まって以来の一番少ない命令系

$ grep ^vdiv /usr/lib/gcc/aarch64-linux-gnu/7.5.0/include/arm_neon.h
vdiv_f32 (float32x2_t __a, float32x2_t __b)
vdiv_f64 (float64x1_t __a, float64x1_t __b)
vdivq_f32 (float32x4_t __a, float32x4_t __b)
vdivq_f64 (float64x2_t __a, float64x2_t __b)
vdiv_f16 (float16x4_t __a, float16x4_t __b)
vdivq_f16 (float16x8_t __a, float16x8_t __b)

div命令は、浮動小数点数型にしか提供されておらず、整数型の除算命令はNEONには無い。
浮動小数点演算ではなく、固定小数点演算しろ、というメッセージを感じる（感想）

`div`命令

何回も書いているが、qが命令のあとにつくのは、128bit幅の演算。
vdiv_で始まる命令は64bit幅のレジスタを2つ引数にとり、 floatを2つか、doubleを1つか、float16_t(いわゆるfp16)を4つまとめたベクトルを引数に取る¹
vdivq_で始まる命令は128bit幅のレジスタを2つ引数に取り、 floatを4つか、doubleを2つか、float16_t(いわゆるfp16)を8つまとめたベクトルを引数に取る
各要素において、1つ目の引数/2つ目の引数 を計算する
こちらも何度も書いているが、f16命令は Arm v8.2 の拡張命令

OpenCVでの参考

OpenCV では、universal intrinsicという名前で各アーキテクチャのSIMD命令のラッパーが提供されている
NEON実装もあるので、参考に見てみよう

intrin_neon.hpp

#define OPENCV_HAL_IMPL_NEON_BIN_OP(bin_op, _Tpvec, intrin) \
inline _Tpvec operator bin_op (const _Tpvec& a, const _Tpvec& b) \
{ \
    return _Tpvec(intrin(a.val, b.val)); \
} \
inline _Tpvec& operator bin_op##= (_Tpvec& a, const _Tpvec& b) \
{ \
    a.val = intrin(a.val, b.val); \
    return a; \
}


#if CV_SIMD128_64F
OPENCV_HAL_IMPL_NEON_BIN_OP(/, v_float32x4, vdivq_f32)
OPENCV_HAL_IMPL_NEON_BIN_OP(+, v_float64x2, vaddq_f64)
OPENCV_HAL_IMPL_NEON_BIN_OP(-, v_float64x2, vsubq_f64)
OPENCV_HAL_IMPL_NEON_BIN_OP(*, v_float64x2, vmulq_f64)
OPENCV_HAL_IMPL_NEON_BIN_OP(/, v_float64x2, vdivq_f64)
#else
inline v_float32x4 operator / (const v_float32x4& a, const v_float32x4& b)
{
    float32x4_t reciprocal = vrecpeq_f32(b.val);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    return v_float32x4(vmulq_f32(a.val, reciprocal));
}
inline v_float32x4& operator /= (v_float32x4& a, const v_float32x4& b)
{
    float32x4_t reciprocal = vrecpeq_f32(b.val);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    a.val = vmulq_f32(a.val, reciprocal);
    return a;
}
#endif

最初はマクロの定義である。二項演算子に対しNEONの命令を1つ対応させる。


#define OPENCV_HAL_IMPL_NEON_BIN_OP(bin_op, _Tpvec, intrin) \
inline _Tpvec operator bin_op (const _Tpvec& a, const _Tpvec& b) \
{ \
    return _Tpvec(intrin(a.val, b.val)); \
} \
inline _Tpvec& operator bin_op##= (_Tpvec& a, const _Tpvec& b) \
{ \
    a.val = intrin(a.val, b.val); \
    return a; \
}

ここで、CV_SIMD128_64Fというマクロが登場する

#if CV_SIMD128_64F

これは当該SIMD命令がdoubleの演算をサポートするか否かを表すdefine
NEONにおいては、Arm v8 か否か、を表すフラグでもある

#if CV_SIMD128_64F
OPENCV_HAL_IMPL_NEON_BIN_OP(/, v_float32x4, vdivq_f32)
OPENCV_HAL_IMPL_NEON_BIN_OP(+, v_float64x2, vaddq_f64)
OPENCV_HAL_IMPL_NEON_BIN_OP(-, v_float64x2, vsubq_f64)
OPENCV_HAL_IMPL_NEON_BIN_OP(*, v_float64x2, vmulq_f64)
OPENCV_HAL_IMPL_NEON_BIN_OP(/, v_float64x2, vdivq_f64)
#else
inline v_float32x4 operator / (const v_float32x4& a, const v_float32x4& b)
{
    float32x4_t reciprocal = vrecpeq_f32(b.val);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    return v_float32x4(vmulq_f32(a.val, reciprocal));
}
inline v_float32x4& operator /= (v_float32x4& a, const v_float32x4& b)
{
    float32x4_t reciprocal = vrecpeq_f32(b.val);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    a.val = vmulq_f32(a.val, reciprocal);
    return a;
}
#endif

で本題のifdefの中身だが、Arm v8の場合はストレートにvdivq_f32およびvdivq_f64命令を使っている
問題はArm v7 以前の場合。
何やら複雑に複数の命令を呼んでいる

inline v_float32x4 operator / (const v_float32x4& a, const v_float32x4& b)
{
    float32x4_t reciprocal = vrecpeq_f32(b.val);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    return v_float32x4(vmulq_f32(a.val, reciprocal));
}
inline v_float32x4& operator /= (v_float32x4& a, const v_float32x4& b)
{
    float32x4_t reciprocal = vrecpeq_f32(b.val);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    reciprocal = vmulq_f32(vrecpsq_f32(b.val, reciprocal), reciprocal);
    a.val = vmulq_f32(a.val, reciprocal);
    return a;
}

これは後日ネタにする予定だが、Arm v7以前のNEONには、浮動小数点演算の除算命令が一切無い。
doubleのサポート云々関係なく、そもそもfloatでも割り算命令自体がArm v7のNEONには存在しない。
予想だが、Armはもともと組み込み向けのプロセッサであり、消費電力や回路規模には厳しい顧客向けに作られていた。
で、浮動小数点回路はただでさえ大きくなりがちなのだが、除算回路はその中でも大きい。なのでそこを回路規模小さくするために浮動小数点の除算命令を除外したのだと思う。
それでも計算したい場合向けにrecpsq命令があるのだが、それはまた別の日のお話

おわりに

今日は数少ないdiv命令を紹介しました
明日も手島の予定で、subかloadを紹介します

double1つを引数にとって割り算って、もはやSI M D命令では無い疑惑。 ↩