SIMD命令によるプログラム高速化

　こんにちは。ディマージシェアの技術担当です。[以前の記事](https://qiita.com/ob_nullpo/items/f6037f8622c848d8cc3b)で小数点演算を整数演算で近似してパフォーマンス改善を図る方法を紹介しました。今回は、更にSIMD命令を用いてパフォーマンスを改善してみようと思います。

# SIMD命令
　SIMDはsingle instruction, multiple dataの略で、その名の通り、1つの命令で複数のデータを処理します。最近のIntel系CPUではAVX2と呼ばれる命令セットが使用でき、1命令で16個の16bit整数（256bitレジスタ）を処理できます。詳細はgoogleで調べると資料はたくさんあるので、ここでは割愛します。
　SIMD命令は1回の演算で複数のデータを処理できるので高速化に利用できます。理想はコンパイラが自動で良いコードを出力してくれれば良いのですが、現状は難しいようで、ソースコード内にSIMD命令を使用するように明示します。C言語では、SIMD命令と関数がほぼ1:1対応の組み込み関数が用意されています。「avx2 intrinsics」のようなキーワードで調べると出てきます。

# SIMDを用いたコーディング
```cpp
       for (i = 0; i < dot_ct; i++) {
           int16_t y = ((ss * dots_x[i]) >> 16) + ((cc * dots_y[i]) >> 16) + ooy;
           horizontal_sum[x][y] += 1;
       }
```

このfor文をSIMD化してみます。
```cpp
       __m256i ymm0, ymm1, ymm_ss, ymm_cc, ymm_ooy;

       ymm_ss = _mm256_set1_epi16(ss);
       ymm_cc = _mm256_set1_epi16(cc);
       ymm_ooy = _mm256_set1_epi16(ooy);

       for (i = 0; i < (dot_ct & 0xFFFFFFF0); i += 16) {
           ymm0 = _mm256_load_si256((__m256i*)&dots_x[i]);
           ymm1 = _mm256_load_si256((__m256i*)&dots_y[i]);

           ymm0 = _mm256_mulhi_epi16(ymm0, ymm_ss);
           ymm1 = _mm256_mulhi_epi16(ymm1, ymm_cc);
           ymm0 = _mm256_add_epi16(ymm0, ymm1);
           ymm0 = _mm256_add_epi16(ymm0, ymm_ooy);

           ymm1 = _mm256_permute2f128_si256(ymm0, ymm0, 1);

           uint64_t tmp = _mm256_cvtsi256_si64(ymm0);
           int16_t y;
           
           y = tmp & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 16) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 32) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 48) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           
           ymm0 = _mm256_srli_si256(ymm0, 8);
           tmp = _mm256_cvtsi256_si64(ymm0);
           y = tmp & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 16) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 32) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 48) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           
           tmp = _mm256_cvtsi256_si64(ymm1);
           y = tmp & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 16) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 32) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 48) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;

           ymm1 = _mm256_srli_si256(ymm1, 8);
           tmp = _mm256_cvtsi256_si64(ymm1);
           y = tmp & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 16) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 32) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
           y = (tmp >> 48) & 0x0000FFFF;
           horizontal_sum[x][y] += 1;
       }
       for (; i < dot_ct; i++) {
           int16_t y = ((ss * dots_x[i]) >> 16) + ((cc * dots_y[i]) >> 16) + ooy;
           horizontal_sum[x][y] += 1;
       }
```
SIMD化するとこのようになります。もともとはシンプルだったfor文ですが、ちょっと複雑になります。

# ベンチマーク
SIMD化前
```
time = 0.864623 sec
```
SIMD化後
```
time = 0.500570 sec
```
1.73倍の性能になりました。

# まとめ
　プログラムの高速化は、様々な手法があり、今回はSIMD命令を採用してみました。SIMD命令は1命令でたくさんのデータを処理できますが、計算結果を取り出すのに少し手間がかかります。今回は、乗算命令は加算やビット演算と比べるとコストが高いことに着目しています。乗算がSIMD命令で1/16回で済む代わりに、結果の取り出しにたくさんのビット演算を行っています。それでも結果としてパフォーマンスが改善しています。CPUの命令レベルでは同じ1命令でも、命令ごとにコストが違います。命令ごとのコスト差を意識することで、より高速なコードを書くことができます。

# 付録
[ソースコードへのリンク](https://github.com/HirokiHayakawa/inclination_detection)