Alpha Max Plus Beta Min Algorithm

画像処理フィルタを書いていると頻出する計算に次のようなものがあります。

X = sqrt(A * A + B * B)

要するに二つの数の二乗の和の平方根を求めるわけです。
普通にC等で書くならば特に問題ないわけですが、いざSIMD化しようとなるととたんに面倒になります。

まず整数演算の場合、掛け算は16bit同士からしかできません。
よって8bit同士の掛け算をしたい場合は、両者をまず16bitに変換しなければなりません。

__m128i xmm0 = _mm_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
__m128i zero = _mm_setzero_si128()
__m128i xmm1 = _mm_unpackhi_epi16(xmm0, zero);
xmm0 = _mm_unpacklo_epi16(xmm0, zero);

更に、16bit同士の掛け算は最大で32bitになる可能性があるため(65,535 x 65,535 = 4,294,836,225)、出力は32bitになります。
しかも素直に32bitで出力してくれるような命令は存在せず、上位16bitと下位16bitがわかれた状態で出力されるか、積和計算を利用するかしなければなりません。

__m128i xmm2 = _mm_mullo_epi16(xmm0, xmm0);
__m128i xmm3 = _mm_mulhi_epi16(xmm0, xmm0);
xmm0 = _mm_unpacklo_epi16(xmm2, xmm3);
xmm2 = _mm_unpackhi_epi16(xmm2m xmm3);

xmm3 = _mm_unpackhi_epi16(xmm1, zero);
xmm1 = _mm_unpacklo_epi16(xmm1, zero);
xmm1 = _mm_madd_epi16(xmm1, xmm1);
xmm3 = _mm_madd_epi16(xmm3, xmm3);

ただ二乗するだけで、最初は128bitのXMMレジスタ一本で済んでいたのがいつのまにやら4本まで膨れあがり、コードもすでに12行です。
この操作をもう一つの数に対しても行い全部で8本になったレジスタを足しあわせて再び4本のレジスタになったあと、今度はこれらの平方根を求めなければなりません。
しかし、平方根を求める命令は整数演算にはないので、32bit整数を単精度浮動小数点に変換してから行うことになります。

__m128 flt0 = _mm_cvtepi32_ps(xmm0)
flt0 = _mm_sqrt_ps(flt0)
xmm0 = _mm_cvtps_epi32(flt0)

これを4回行った後、最後に再び8bit整数に戻します。

xmm0 = _mm_packs_epi32(xmm0, xmm1);
xmm1 = _mm_packs_epi32(xmm2, xmm3);
xmm0 = _mm_packus_epi16(xmm0, xmm1);

なんでたかがsqrt(A*A+B*B)を行うためにこんな辛い目にあわねばならないのでしょう。
しかもこれで計算できるのは一度にたったの16個です(SSE/SSE2の場合)。
普通にCで書いて16回計算するのとスピードは変わらないか下手すると遅くなってしまいます。
例えばTEdgeMaskはtriticalプラグインでは珍しくクリティカルパートがSIMD化されていないのですが、おそらくはこの計算が出てくるためにやらなかったのでしょう。

というわけで、もしどうしてもスピードがほしいなら普通はこんなに馬鹿正直な計算はしません。

sqrt(A*A+B*B) -> abs(A)+abs(B)

正直、手抜きにも程があるとは思います。
しかし、実用上これで問題が出るケースは意外と少ないので、これで済ませてしまったりするのですね。

とまあ、コストとパフォーマンスの兼ね合いという人類普遍のテーマについてつらつらと考えていたところ、prunedtree氏が"alpha max plus beta min algorithm"を教えてくれました。

sqrt(A*A+B*B) -> max(abs(A),abs(B))*alpha+min(abs(A),abs(B))*beta

alphaとbetaの選び方によっては最大誤差3.96%、平均誤差1.3%……これはイイ！
流石はgradfun2dbやFRfunシリーズを書いた人です。こういうテクニックをよく知ってらっしゃる。
抜け道というものは探せばあるものなんだなぁと思いました。


ちなみにprunedtree氏は、昔はmarcFDという名前でした。
なんでも若気の至りによって犯した過ち(DeenやaWarpSharpのことだと思われる)を悔いて改名したそうです。