この記事はOpenCV Advent Calendar 2023の7日目の記事です。
他の記事は目次にまとめられています。
TL;DR
- Round処理で「0.5」ぴったりのとき、切り上げるか？切り下げるか？
- IEE754では「推奨」はあっても、「厳密にこれを使え」が決まっていない！
- 偶数側に寄せるというのが一般的
- Rounding{9.5, 10.5, 11.5, 12.5} => {10, 10, 12, 12}
- OpenCV ARM向け実装では実装不足だったので、パッチ書いたけどマージされない
この記事を書いた人
- JTCで、元組み込みエンジニア。今は(準)社内ニート
■ (1+0)/2 = 1 から始まる物語
その日、非常に面白いissueが発行された。 (1+0)/2 の結果が、ARM MacとIntel Linuxで異なる、というのだ、な、なんだってーー！！
cv::Mat A = cv::Mat::ones(3, 3, CV_8UC1);
cv::Mat B = cv::Mat::zeros(3, 3, CV_8UC1);
cv::Mat mean = (A + B) / 2;
std::cout << "A = \n" << A << std::endl;
std::cout << "B = \n" << B << std::endl;
std::cout << "(A+B)/2 = \n" << mean << std::endl;
(A+B)/2 =
[ 1, 1, 1;
1, 1, 1;
1, 1, 0]
(A+B)/2 =
[ 0, 0, 0;
0, 0, 0;
0, 0, 0]
◯全部1ではない！
ここで、ARM Macの結果をよく見ると、なんと「全部1」ではない。「先頭8つが1で、最後が0」になっている。ということは……これはOptimize処理に関する話だと推定できる。更に、最適化処理を外したら、ちゃんと動くようになった。 勝ったな、ガハハハハ!! (負けフラグ「やぁ」)
■何がこの問題を引き起こしたのか
さて、そこで原因調査をしてみたのだが、なかなかNEON実装の闇が深かった。ARMのHAL実装(Hardware Abstraction Layer) は、3rdpart/catotene/ 以下に実装されている。その最後の加算部分。
◯実装を確認してみる
(1+0)/2 は、(1+0) * 0.5 -> 1 * 0.5 + 0 * 0.5 -> 0.5 + 0.0 で処理される。
template <> struct wAdd<u32>
{
typedef u32 type;
f32 alpha, beta, gamma;
float32x4_t valpha, vbeta, vgamma;
wAdd(f32 _alpha, f32 _beta, f32 _gamma):
alpha(_alpha), beta(_beta), gamma(_gamma)
{
valpha = vdupq_n_f32(_alpha);
vbeta = vdupq_n_f32(_beta);
vgamma = vdupq_n_f32(_gamma + 0.5);
}
void operator() (const VecTraits<u32>::vec128 & v_src0,
const VecTraits<u32>::vec128 & v_src1,
VecTraits<u32>::vec128 & v_dst) const
{
float32x4_t vs1 = vcvtq_f32_u32(v_src0);
float32x4_t vs2 = vcvtq_f32_u32(v_src1);
vs1 = vmlaq_f32(vgamma, vs1, valpha);
vs1 = vmlaq_f32(vs1, vs2, vbeta);
v_dst = vcvtq_u32_f32(vs1);
}
順番を追っていくと、そこまで難しくない。
- valpha = vbata = 0.5 ※ 1/2 なので0.5
- vgamma = 0.5 ※後述
- vs1 に、v_src0のアドレスに入っているデータ、f32として読み出す
- vs2 に、v_src1のアドレスに入っているデータ、f32として読み出す
- vs1 に、( vs1 * valpha ) + vgamma を代入
- vs1 に、( vs2 * vbeta ) + vs1 を代入
- v_dstに、vs1をf32からu32に変換(!)して代入
◯バグの根幹
このロジックのバグは、vcvtq_u32_f32()にある。ARM社のマニュアルを確認する。
(邦訳) 浮動小数点を符号なし固定小数点に変換し、ゼロ (ベクトル) に向かって丸めます。この命令は、ゼロ方向への丸めモードを使用して、スカラーまたはベクトル内の各要素を浮動小数点から固定小数点の符号なし整数に変換し、その結果を汎用デスティネーション レジスタに書き込みます。
英語原文
Floating-point Convert to Unsigned fixed-point, rounding toward Zero (vector). This instruction converts a scalar or each element in a vector from floating-point to fixed-point unsigned integer using the Round towards Zero rounding mode, and writes the result to the general-purpose destination register.
平坦な言い方を「（正数でいえば）切り捨て処理(truncation)」なのである。そして、事前にgamma=0.5のバイアスを加える事で、「四捨五入」になっている。
| (1) src | 0.00 | 0.25 | 0.50 | 0.75 | 1.00 |
|---|---|---|---|---|---|
| (2) src+0.5 | 0.50 | 0.75 | 1.00 | 1.25 | 1.50 |
| (3) trunc | 0 | 0 | 1 | 1 | 2 |
この時点で、２つのバグがある、お気づきでしょうか……
◯負の数
「負の数」の場合、gamma= -0.5のバイアスで考えてから、0方向に丸める必要がある。
＜現状実装＞ +0.5のバイアスだと、-1.0も0になる！！
| (1) src | 0.00 | -0.25 | -0.50 | -0.75 | -1.00 |
|---|---|---|---|---|---|
| (2) src+0.5 | 0.50 | 0.25 | 0.00 | -0.25 | -0.50 |
| (3) trunc | 0 | 0 | 0 | 0 | 0 |
＜理想＞ 符号に応じて+0.5/-0.5を使い分ければ、ただしい丸めになる！
| (1) src | 0.00 | -0.25 | -0.50 | -0.75 | -1.00 |
|---|---|---|---|---|---|
| (2) src-0.5 | 0.50 | 0.75 | -1.00 | -1.25 | -1.50 |
| (3) trunc | 0 | 0 | -1 | -1 | -1 |
◯最近接丸め（偶数）
biasの符号だけでは、元々のバグは説明できない。ではなぜか？ IEE754に関するwikipediaページを確認してみよう。
浮動小数点数の丸め
IEEE 754-2008標準では5種類の丸めアルゴリズムが定義されている。うち2種類は最近接な値に丸める方法であり、残り3種類は方向丸めと呼ぶ。最近接丸め
- 最近接丸め（偶数）
- 最も近くの表現できる値へ丸める。表現可能な2つの値の中間の値であったら、一番低い仮数ビット（桁）が0になるほうを採用する。これは二進での標準動作かつ十進でも推奨となっている。
つまり、以下のようになる。
- 余りが無い、偶数・奇数は、そのままの値
- 余りが0.5だと、偶数側に近づく
- 10.5だと 10 or 11 が丸め結果の候補。偶数である10を選択
- 11.5だと 11 or 12 が丸め結果の候補。偶数である12を選択
| src | 9.5 | 10.0 | 10.5 | 11.0 | 11.5 | 12.0 | 12.5 | 13.0 | 13.5 | 14.0 |
|---|---|---|---|---|---|---|---|---|---|---|
| dst | 10 | 10 | 10 | 11 | 12 | 12 | 12 | 13 | 14 | 14 |
この 「二進での標準動作かつ十進でも推奨となっている」 がポイント！なんと、実行環境によっては別の振る舞いをする場合があるのだ。そして、悲劇的な事に、ARMv7には最近接丸め（偶数）のNEON命令は存在しない！
■それでは、どうやって問題を解決するのか？
ここでは、ARMv7とARMv8それぞれで、異なるアプローチで問題解決を図る。
◯ARMv7
「符号を外して+0.5 biasで計算、偶数に丸め、最後に符号を戻す！」と、まあ、なかなかメンドクサイ。とはいえ、符号に応じてbiasをプラス井マイナス切り替えて作るよりかは楽（のような気がする）。
static const float32x4_t v0_0_f32 = vdupq_n_f32(0.0);
static const float32x4_t v0_5_f32 = vdupq_n_f32(0.5);
static const int32x4_t v1_0_s32 = vdupq_n_s32(1);
const int32x4_t val_positive = vreinterpretq_s32_u32( vcgtq_f32( val, v0_0_f32 ) );
const int32x4_t ret_signs = vsubq_s32(
vandq_s32( v1_0_s32, val_positive ),
vbicq_s32( v1_0_s32, val_positive ) );
const float32x4_t val_abs = vabsq_f32( val );
const int32x4_t round = vcvtq_s32_f32( vaddq_f32( val_abs, v0_5_f32 ) );
const int32x4_t odd = vandq_s32( round, v1_0_s32 );
const float32x4_t diff = vsubq_f32( vcvtq_f32_s32(round), val_abs);
const int32x4_t round_down = vandq_s32( odd, vreinterpretq_s32_u32( vceqq_f32( diff,v0_5_f32 ) ) );
const int32x4_t ret_abs = vsubq_s32( round, round_down );
const int32x4_t ret = vmulq_s32( ret_abs, ret_signs );
return ret;
- 変換したい値を、valとする
- val_positive に、valが0.0より大きいかどうかの判定を代入
- ret_signs に、（最終結果の符号）0.0より大きければ+1.0、そうでなければ-1.0を代入
- val_abs に。元の値の絶対値を代入
- round に、元の値に0.5を足してtruncateした値を代入（仮）
- odd に、奇数か判定結果を代入
- diff に、roundとval_absの差分を代入。つまり、少数点部の取得
- round_down に、少数点部が0.5であるかを判定
- ret_abs に、結果が偶数になるように丸めを実行
- ret に、ret_absに、ret_signの符号をつけて結果を求める。
◯ARMv8
ARMv8以後のNEONを使うと、これが一行。
return vcvtnq_s32_f32(val);
（邦訳）浮動小数点を符号付き整数に変換し、偶数 (ベクトル) に近い値に丸めます。この命令は、Round to Nearest 丸めモードを使用して、スカラーまたはベクトル内の各要素を浮動小数点値から符号付き整数値に変換し、結果を SIMD&FP デスティネーション レジスタに書き込みます。
英語原文
Floating-point Convert to Signed integer, rounding to nearest with ties to even (vector). This instruction converts a scalar or each element in a vector from a floating-point value to a signed integer value using the Round to Nearest rounding mode, and writes the result to the SIMD&FP destination register.
■パッチの状況は？
まーじされませんなー(12/07)
■まとめ
- Round処理で「0.5」ぴったりのとき、切り上げるか？切り下げるか？
- IEE754では「推奨」はあっても、「厳密にこれを使え」が決まっていない！
- 偶数側に寄せるというのが一般的
- Rounding{9.5, 10.5, 11.5, 12.5} => {10, 10, 12, 12}
- OpenCV ARM向け実装では実装不足だったので、パッチ書いたけどマージされない
以上です、ご精読ありがとうございました！