64ビット整数のパック8ビット整数を1ずつ並列に減算する、ハードウェアSIMDなしのSWAR

質問

64ビット整数があり、それを8つの要素を持つパックされた8ビット整数の配列として解釈している場合。定数 1 を引き、ある要素の結果が他の要素の結果に影響を与えることなくオーバーフローを処理する必要があります。

私は現在このコードを持っていて、それは動作しますが、並列で各パックされた8ビット整数の減算を行い、メモリアクセスを行わないソリューションが必要です。 x86 では、次のような SIMD 命令を使用できました。 psubb のようなSIMD命令が使えるのですが、私がコーディングしているプラットフォームはSIMD命令をサポートしていません。 (この場合は RISC-V) です。

そこで、私は SWAR (レジスタ内SIMD) のバイト間のキャリー伝搬を手動で打ち消すことを試みています。 uint64_t のバイト間のキャリー伝搬を手動でキャンセルすることができます。

uint64_t sub(uint64_t arg) {
    uint8_t* packed = (uint8_t*) &arg;

    for (size_t i = 0; i < sizeof(uint64_t); ++i) {
        packed[i] -= 1;
    }

    return arg;
}

ビット演算子でできそうな気がするけど、よくわからない。SIMD命令を使用しない解決策を探しています。C または C++ でのソリューションで、かなり移植性の高いもの、あるいは私自身のソリューションを実装できるように、その背後にある理論だけを探しています。

どのように解決するのですか?

効率的なSIMD命令、SSE/MMXを搭載したCPUをお持ちの場合 paddb ( _mm_add_epi8 ) も実行可能です。 Peter Cordesの回答 は、GNU C (gcc/clang) のベクトル構文、およびストリクト・エイリアシング UB に対する安全性についても説明しています。 私はその回答もレビューすることを強くお勧めします。

を使って自分でやる uint64_t にアクセスする際には、アラインメントの問題やストリクト・エイリアシングのUBを避けるために注意が必要です。 uint8_t 配列に uint64_t* . その部分を省いたのは、データを最初から uint64_t で始めることでその部分を省きましたが、GNU Cでは may_alias 型定義は問題を解決します (それについてはピーターの答えか memcpy ).

そうでなければ、データの割り当てや宣言は uint64_t として割り当て、それに uint8_t* でアクセスします。 unsigned char* は何でもエイリアスにすることができるので、 8ビット要素特有の問題を回避することができます。 (もし uint8_t が存在する場合、それはおそらく unsigned char .)

これは以前の不正なアルゴリズムからの変更であることに注意してください（改訂履歴を参照）。

これは、任意の引き算のためにループすることなく可能であり、以下のような既知の定数に対してより効率的になります。 1 のような既知の定数に対してより効率的になります。 主なトリックは、ハイビットを設定することで各バイトからのキャリーアウトを防ぎ、その後、減算結果を修正することです。

与えられた引き算のテクニックを少し最適化することにします。 ここで . 定義しているのです。

SWAR sub z = x - y
    z = ((x | H) - (y &~H)) ^ ((x ^~y) & H)

と H と定義されています。 0x8080808080808080U (すなわち、各パックされた整数のMSB)です。デクリメントの場合は y は 0x0101010101010101U .

私たちは、以下のことを知っています。 y はそのMSBがすべてクリアされているので、マスクのステップをひとつスキップすることができます（つまり y & ~H と同じです。 y と同じです）。計算は次のように進みます。

1. の各コンポーネントの MSB を設定します。 x の各コンポーネントの MSB を 1 に設定し、借用が MSB を越えて次のコンポーネントに伝搬しないようにします。これを調整済み入力と呼びます。
2. 各コンポーネントから 1 を減算します。 0x01010101010101 を補正された入力から引きます。これは、ステップ1のおかげでコンポーネント間の借用が発生しません。これを調整済み出力と呼びます。
3. 次に結果の MSB を修正する必要があります。結果を修正するために、元の入力のMSBを反転させたものと調整された出力をxorします。

この演算は次のように書くことができます。

#define U64MASK 0x0101010101010101U
#define MSBON 0x8080808080808080U
uint64_t decEach(uint64_t i){
      return ((i | MSBON) - U64MASK) ^ ((i ^ MSBON) & MSBON);
}

好ましくは、コンパイラによってインライン化されることです( コンパイラディレクティブ を使って強制する)、または式が他の関数の一部としてインラインで書かれていることが望ましいです。

テストケースです。

in:  0000000000000000
out: ffffffffffffffff

in:  f200000015000013
out: f1ffffff14ffff12

in:  0000000000000100
out: ffffffffffff00ff

in:  808080807f7f7f7f
out: 7f7f7f7f7e7e7e7e

in:  0101010101010101
out: 0000000000000000

パフォーマンス詳細

この関数を1回だけ呼び出した場合のx86_64アセンブリです。パフォーマンスを向上させるために、定数をできるだけ長くレジスタに保存できるようにインライン化する必要があります。定数がレジスタにあるようなタイトなループでは、実際のデクリメントには、最適化後に or+not+and+add+xor の5つの命令が必要です。コンパイラの最適化を打ち負かすような代替案は見当たりません。

uint64t[rax] decEach(rcx):
    movabs  rcx, -9187201950435737472
    mov     rdx, rdi
    or      rdx, rcx
    movabs  rax, -72340172838076673
    add     rax, rdx
    and     rdi, rcx
    xor     rdi, rcx
    xor     rax, rdi
    ret

以下のスニペットのいくつかのIACAテストと。

// Repeat the SWAR dec in a loop as a microbenchmark
uint64_t perftest(uint64_t dummyArg){
    uint64_t dummyCounter = 0;
    uint64_t i = 0x74656a6d27080100U; // another dummy value.
    while(i ^ dummyArg) {
        IACA_START
        uint64_t naive = i - U64MASK;
        i = naive + ((i ^ naive ^ U64MASK) & U64MASK);
        dummyCounter++;
    }
    IACA_END
    return dummyCounter;
}

Skylake マシンで、デクリメント、xor、および compare+jump を実行すると、1 反復あたりわずか 5 サイクル未満で実行できることを示すことができます。

Throughput Analysis Report
--------------------------
Block Throughput: 4.96 Cycles       Throughput Bottleneck: Backend
Loop Count:  26
Port Binding In Cycles Per Iteration:
--------------------------------------------------------------------------------------------------
|  Port  |   0   -  DV   |   1   |   2   -  D    |   3   -  D    |   4   |   5   |   6   |   7   |
--------------------------------------------------------------------------------------------------
| Cycles |  1.5     0.0  |  1.5  |  0.0     0.0  |  0.0     0.0  |  0.0  |  1.5  |  1.5  |  0.0  |
--------------------------------------------------------------------------------------------------

(もちろん、x86-64 では、単にロードする、あるいは movq をXMMのREGにロードして paddb というように、RISC-VのようなISAでどのようにコンパイルされるかを見る方が興味深いかもしれません)。