SSE整数除法？
Posted
技术标签:
【中文标题】SSE整数除法？【英文标题】：SSE integer division? 【发布时间】：2013-05-29 19:58:05 【问题描述】：_mm_div_ps 用于浮点值除法，_mm_mullo_epi16 用于整数乘法。但是有整数除法（16位值）吗？我该如何进行这样的划分？
【问题讨论】：
没有足够的人需要它，所以... 不，它不存在。这可能是需要它的人数不足以及整数除法单元的难度（和模具空间）的组合。 要除以常数还是除以变量？ 一种快速检查方法：使用 gcc 编译typedef short vec __attribute__((vector_size(16))); vec f(vec x) return x/7; vec g(vec x,vec y) return x/y; 并查看它如何模拟常数除法，但将全除法降低为标量。如果您能找到更好的代码，请通过向 gcc 提交错误报告来提供帮助。 （clang 应该是一样的，但我没有检查）
另见：libdivide
【参考方案1】：
数学说确实可以走得更快
Agner Fog 的 (http://www.agner.org/optimize/#vectorclass) 方法在使用单个除数进行除法时效果很好。此外，如果除数在编译时已知，或者在运行时不经常更改，则此方法还有更多好处。
但是，当对__m128i 元素执行 SIMD 除法时，编译时除数和被除数都没有可用信息，我们别无选择，只能转换为浮点数并执行计算。另一方面，使用_mm_div_ps 不会给我们带来惊人的速度提升，因为这条指令在大多数微架构上具有 11 到 14 个周期的可变延迟，如果我们考虑 Knights Landing，有时可以达到 38 个周期。最重要的是，该指令没有完全流水线化，并且根据微架构具有 3-6 个周期的倒数吞吐量。
但是我们可以避免使用_mm_div_ps，而使用_mm_rcp_ss。
不幸的是，__m128 _mm_rcp_ss (__m128 a) 之所以快只是因为它提供了近似值。即（取自Intel Intrnisics Guide）：
计算 a 中低位单精度（32 位）浮点元素的近似倒数，将结果存储到 dst 的低位元素中，并将高位 3 个压缩元素从 a 复制到 dst 的高位元素.此近似值的最大相对误差小于 1.5*2^-12。
因此，要从_mm_rcp_ss 中受益，我们需要补偿近似造成的损失。在这个方向上做了很棒的工作，在Improved division by invariant integers by Niels Möller and Torbjörn Granlund：
由于当前处理器中缺乏有效的除法指令，除法是使用预先计算的除数倒数的单字近似值进行乘法运算，然后执行几个调整步骤。
要计算 16 位有符号整数除法，我们只需要一个调整步骤，并据此确定我们的解决方案。
SSE2
static inline __m128i _mm_div_epi16(const __m128i &a_epi16, const __m128i &b_epi16)
//
// Setup the constants.
//
const __m128 two = _mm_set1_ps(2.00000051757f);
const __m128i lo_mask = _mm_set1_epi32(0xFFFF);
//
// Convert to two 32-bit integers
//
const __m128i a_hi_epi32 = _mm_srai_epi32(a_epi16, 16);
const __m128i a_lo_epi32_shift = _mm_slli_epi32(a_epi16, 16);
const __m128i a_lo_epi32 = _mm_srai_epi32(a_lo_epi32_shift, 16);
const __m128i b_hi_epi32 = _mm_srai_epi32(b_epi16, 16);
const __m128i b_lo_epi32_shift = _mm_slli_epi32(b_epi16, 16);
const __m128i b_lo_epi32 = _mm_srai_epi32(b_lo_epi32_shift, 16);
//
// Convert to 32-bit floats
//
const __m128 a_hi = _mm_cvtepi32_ps(a_hi_epi32);
const __m128 a_lo = _mm_cvtepi32_ps(a_lo_epi32);
const __m128 b_hi = _mm_cvtepi32_ps(b_hi_epi32);
const __m128 b_lo = _mm_cvtepi32_ps(b_lo_epi32);
//
// Calculate the reciprocal
//
const __m128 b_hi_rcp = _mm_rcp_ps(b_hi);
const __m128 b_lo_rcp = _mm_rcp_ps(b_lo);
//
// Calculate the inverse
//
#ifdef __FMA__
const __m128 b_hi_inv_1 = _mm_fnmadd_ps(b_hi_rcp, b_hi, two);
const __m128 b_lo_inv_1 = _mm_fnmadd_ps(b_lo_rcp, b_lo, two);
#else
const __m128 b_mul_hi = _mm_mul_ps(b_hi_rcp, b_hi);
const __m128 b_mul_lo = _mm_mul_ps(b_lo_rcp, b_lo);
const __m128 b_hi_inv_1 = _mm_sub_ps(two, b_mul_hi);
const __m128 b_lo_inv_1 = _mm_sub_ps(two, b_mul_lo);
#endif
//
// Compensate for the loss
//
const __m128 b_hi_rcp_1 = _mm_mul_ps(b_hi_rcp, b_hi_inv_1);
const __m128 b_lo_rcp_1 = _mm_mul_ps(b_lo_rcp, b_lo_inv_1);
//
// Perform the division by multiplication
//
const __m128 hi = _mm_mul_ps(a_hi, b_hi_rcp_1);
const __m128 lo = _mm_mul_ps(a_lo, b_lo_rcp_1);
//
// Convert back to integers
//
const __m128i hi_epi32 = _mm_cvttps_epi32(hi);
const __m128i lo_epi32 = _mm_cvttps_epi32(lo);
//
// Zero-out the unnecessary parts
//
const __m128i hi_epi32_shift = _mm_slli_epi32(hi_epi32, 16);
#ifdef __SSE4_1__
//
// Blend the bits, and return
//
return _mm_blend_epi16(lo_epi32, hi_epi32_shift, 0xAA);
#else
//
// Blend the bits, and return
//
const __m128i lo_epi32_mask = _mm_and_si128(lo_epi32, const_mm_div_epi16_lo_mask);
return _mm_or_si128(hi_epi32_shift, lo_epi32_mask);
#endif
此解决方案只能使用 SSE2，如果可用，将使用 FMA。但是，使用普通除法可能与使用近似值一样快（甚至更快）。
在存在AVX 的情况下，可以改进此解决方案，因为可以使用一个AVX 寄存器同时处理高和低部分。
验证
由于我们只处理 16 位，因此我们可以使用蛮力测试在几秒钟内轻松验证解决方案的正确性：
void print_epi16(__m128i a)
int i; int16_t tmp[8];
_mm_storeu_si128( (__m128i*) tmp, a);
for (i = 0; i < 8; i += 1)
printf("%8d ", (int) tmp[i]);
printf("\n");
bool run_mm_div_epi16(const int16_t *a, const int16_t *b)
const size_t n = 8;
int16_t result_expected[n];
int16_t result_obtained[n];
//
// Derive the expected result
//
for (size_t i = 0; i < n; i += 1)
result_expected[i] = a[i] / b[i];
//
// Now perform the computation
//
const __m128i va = _mm_loadu_si128((__m128i *) a);
const __m128i vb = _mm_loadu_si128((__m128i *) b);
const __m128i vr = _mm_div_epi16(va, vb);
_mm_storeu_si128((__m128i *) result_obtained, vr);
//
// Check for array equality
//
bool eq = std::equal(result_obtained, result_obtained + n, result_expected);
if (!eq)
cout << "Testing of _mm_div_epi16 failed" << endl << endl;
cout << "a: ";
print_epi16(va);
cout << "b: ";
print_epi16(vb);
cout << endl;
cout << "results_obtained: ";
print_epi16(vr);
cout << "results_expected: ";
print_epi16(_mm_loadu_si128((__m128i *) result_expected));
cout << endl;
return eq;
void test_mm_div_epi16()
const int n = 8;
bool correct = true;
//
// Brute-force testing
//
int16_t a[n];
int16_t b[n];
for (int32_t i = INT16_MIN; correct && i <= INT16_MAX; i += n)
for (int32_t j = 0; j < n; j += 1)
a[j] = (int16_t) (i + j);
for (int32_t j = INT16_MIN; correct && j < 0; j += 1)
const __m128i jv = _mm_set1_epi16((int16_t) j);
_mm_storeu_si128((__m128i *) b, jv);
correct = correct && run_mm_div_epi16(a, b);
for (int32_t j = 1; correct && j <= INT16_MAX; j += 1)
const __m128i jv = _mm_set1_epi16((int16_t) j);
_mm_storeu_si128((__m128i *) b, jv);
correct = correct && run_mm_div_epi16(a, b);
if (correct)
cout << "Done!" << endl;
else
cout << "_mm_div_epi16 can not be validated" << endl;
AVX2
有了上面的解决方案，AVX2 的实现就很简单了：
static inline __m256i _mm256_div_epi16(const __m256i &a_epi16, const __m256i &b_epi16)
//
// Setup the constants.
//
const __m256 two = _mm256_set1_ps(2.00000051757f);
//
// Convert to two 32-bit integers
//
const __m256i a_hi_epi32 = _mm256_srai_epi32(a_epi16, 16);
const __m256i a_lo_epi32_shift = _mm256_slli_epi32(a_epi16, 16);
const __m256i a_lo_epi32 = _mm256_srai_epi32(a_lo_epi32_shift, 16);
const __m256i b_hi_epi32 = _mm256_srai_epi32(b_epi16, 16);
const __m256i b_lo_epi32_shift = _mm256_slli_epi32(b_epi16, 16);
const __m256i b_lo_epi32 = _mm256_srai_epi32(b_lo_epi32_shift, 16);
//
// Convert to 32-bit floats
//
const __m256 a_hi = _mm256_cvtepi32_ps(a_hi_epi32);
const __m256 a_lo = _mm256_cvtepi32_ps(a_lo_epi32);
const __m256 b_hi = _mm256_cvtepi32_ps(b_hi_epi32);
const __m256 b_lo = _mm256_cvtepi32_ps(b_lo_epi32);
//
// Calculate the reciprocal
//
const __m256 b_hi_rcp = _mm256_rcp_ps(b_hi);
const __m256 b_lo_rcp = _mm256_rcp_ps(b_lo);
//
// Calculate the inverse
//
const __m256 b_hi_inv_1 = _mm256_fnmadd_ps(b_hi_rcp, b_hi, two);
const __m256 b_lo_inv_1 = _mm256_fnmadd_ps(b_lo_rcp, b_lo, two);
//
// Compensate for the loss
//
const __m256 b_hi_rcp_1 = _mm256_mul_ps(b_hi_rcp, b_hi_inv_1);
const __m256 b_lo_rcp_1 = _mm256_mul_ps(b_lo_rcp, b_lo_inv_1);
//
// Perform the division by multiplication
//
const __m256 hi = _mm256_mul_ps(a_hi, b_hi_rcp_1);
const __m256 lo = _mm256_mul_ps(a_lo, b_lo_rcp_1);
//
// Convert back to integers
//
const __m256i hi_epi32 = _mm256_cvttps_epi32(hi);
const __m256i lo_epi32 = _mm256_cvttps_epi32(lo);
//
// Blend the low and the high-parts
//
const __m256i hi_epi32_shift = _mm256_slli_epi32(hi_epi32, 16);
return _mm256_blend_epi16(lo_epi32, hi_epi32_shift, 0xAA);
我们可以使用上述相同的方法来执行代码验证。
性能
我们可以使用每个周期的测量触发器 (F/C) 来评估性能。在这种情况下，我们希望查看每个周期可以执行多少个除法。为此，我们定义了两个向量 a 和 b 并执行逐点除法。 a 和 b 都使用 xorshift32 填充随机数据，初始化为 uint32_t state = 3853970173;
我使用RDTSC 测量周期，使用暖缓存执行 15 次重复，并使用中值作为结果。为了避免频率缩放和资源共享对测量的影响，Turbo Boost 和 Hyper-Threading 被禁用。为了运行代码，我使用Intel Xeon CPU E3-1285L v3 3.10GHz Haswell，32GB RAM 和 25.6GB/s 带宽到主内存，运行 Debian GNU/Linux 8 (jessie)，kernel 3.16.43-2+deb8u3。 gcc 使用的是 4.9.2-10。结果如下：
纯SSE2实现
我们将plain division 与上面提出的算法进行比较：
===============================================================
= Compiler & System info
===============================================================
Current CPU : Intel(R) Xeon(R) CPU E3-1285L v3 @ 3.10GHz
CXX Compiler ID : GNU
CXX Compiler Path : /usr/bin/c++
CXX Compiler Version : 4.9.2
CXX Compiler Flags : -O3 -std=c++11 -msse2 -mno-fma
--------------------------------------------------------------------------------
| Size | Division F/C | Division B/W | Approx. F/C | Approximation B/W |
--------------------------------------------------------------------------------
| 128 | 0.5714286 | 26911.45 MB/s | 0.5019608 | 23634.21 MB/s |
| 256 | 0.5714286 | 26909.17 MB/s | 0.5039370 | 23745.44 MB/s |
| 512 | 0.5707915 | 26928.14 MB/s | 0.5039370 | 23763.79 MB/s |
| 1024 | 0.5707915 | 26936.33 MB/s | 0.5039370 | 23776.85 MB/s |
| 2048 | 0.5709507 | 26938.51 MB/s | 0.5039370 | 23780.25 MB/s |
| 4096 | 0.5708711 | 26940.56 MB/s | 0.5039990 | 23782.65 MB/s |
| 8192 | 0.5708711 | 26940.16 MB/s | 0.5039370 | 23781.85 MB/s |
| 16384 | 0.5704735 | 26921.76 MB/s | 0.4954040 | 23379.24 MB/s |
| 32768 | 0.5704537 | 26921.26 MB/s | 0.4954639 | 23382.13 MB/s |
| 65536 | 0.5703147 | 26914.53 MB/s | 0.4943539 | 23330.13 MB/s |
| 131072 | 0.5691680 | 26860.21 MB/s | 0.4929539 | 23264.40 MB/s |
| 262144 | 0.5690618 | 26855.60 MB/s | 0.4929187 | 23262.22 MB/s |
| 524288 | 0.5691378 | 26858.75 MB/s | 0.4929488 | 23263.56 MB/s |
| 1048576 | 0.5677474 | 26794.14 MB/s | 0.4918968 | 23214.34 MB/s |
| 2097152 | 0.5371243 | 25348.39 MB/s | 0.4700511 | 22183.07 MB/s |
| 4194304 | 0.5128146 | 24200.82 MB/s | 0.4529809 | 21377.28 MB/s |
| 8388608 | 0.5036971 | 23770.36 MB/s | 0.4438345 | 20945.84 MB/s |
| 16777216 | 0.5005390 | 23621.14 MB/s | 0.4409909 | 20811.32 MB/s |
| 33554432 | 0.4992792 | 23561.90 MB/s | 0.4399777 | 20763.49 MB/s |
--------------------------------------------------------------------------------
我们可以观察到普通除法将如何比建议的近似步骤略快。在这种情况下，我们可以得出结论，在 Haswell 微架构上使用SSE2 近似值将不是最佳的。
但是，如果我们在旧的 Sandy Bridge 机器上运行相同的结果，例如 Intel(R) Xeon(R) CPU X5680 @ 3.33GHz，我们已经可以看到近似值的好处：
===============================================================
= Compiler & System info
===============================================================
Current CPU : Intel(R) Xeon(R) CPU X5680 @ 3.33GHz
CXX Compiler ID : GNU
CXX Compiler Path : /usr/bin/c++
CXX Compiler Version : 4.8.5
CXX Compiler Flags : -O3 -std=c++11 -msse2 -mno-fma
--------------------------------------------------------------------------------
| Size | Division F/C | Division B/W | Approx. F/C | Approximation B/W |
--------------------------------------------------------------------------------
| 128 | 0.2857143 | 14511.41 MB/s | 0.3720930 | 18899.89 MB/s |
| 256 | 0.2853958 | 14512.51 MB/s | 0.3715530 | 18898.91 MB/s |
| 512 | 0.2853958 | 14510.53 MB/s | 0.3715530 | 18896.44 MB/s |
| 1024 | 0.2853162 | 14511.81 MB/s | 0.3700759 | 18824.00 MB/s |
| 2048 | 0.2853162 | 14511.04 MB/s | 0.3708130 | 18860.31 MB/s |
| 4096 | 0.2852964 | 14511.16 MB/s | 0.3711826 | 18879.27 MB/s |
| 8192 | 0.2852666 | 14510.23 MB/s | 0.3713172 | 18886.39 MB/s |
| 16384 | 0.2852616 | 14509.86 MB/s | 0.3712920 | 18885.60 MB/s |
| 32768 | 0.2852244 | 14507.93 MB/s | 0.3712709 | 18884.86 MB/s |
| 65536 | 0.2851003 | 14501.41 MB/s | 0.3701114 | 18826.14 MB/s |
| 131072 | 0.2850711 | 14499.95 MB/s | 0.3685017 | 18743.58 MB/s |
| 262144 | 0.2850745 | 14500.47 MB/s | 0.3684799 | 18742.78 MB/s |
| 524288 | 0.2848062 | 14486.66 MB/s | 0.3681040 | 18723.63 MB/s |
| 1048576 | 0.2846679 | 14479.64 MB/s | 0.3671284 | 18674.02 MB/s |
| 2097152 | 0.2840133 | 14446.52 MB/s | 0.3664623 | 18640.01 MB/s |
| 4194304 | 0.2745241 | 13963.13 MB/s | 0.3488823 | 17745.24 MB/s |
| 8388608 | 0.2741900 | 13946.39 MB/s | 0.3476036 | 17680.37 MB/s |
| 16777216 | 0.2740689 | 13940.32 MB/s | 0.3477076 | 17685.97 MB/s |
| 33554432 | 0.2746752 | 13970.75 MB/s | 0.3482017 | 17711.36 MB/s |
--------------------------------------------------------------------------------
如果 Nehalem 说的是，看看这在更老的机器上会如何表现会更好（假设它有 RCP 支持）。
SSE41 + FMA 实施
我们将the plain division 与上面提出的算法进行比较，启用FMA 和SSE41
===============================================================
= Compiler & System info
===============================================================
Current CPU : Intel(R) Xeon(R) CPU E3-1285L v3 @ 3.10GHz
CXX Compiler ID : GNU
CXX Compiler Path : /usr/bin/c++
CXX Compiler Version : 4.9.2
CXX Compiler Flags : -O3 -std=c++11 -msse4.1 -mfma
--------------------------------------------------------------------------------
| Size | Division F/C | Division B/W | Approx. F/C | Approximation B/W |
--------------------------------------------------------------------------------
| 128 | 0.5714286 | 26884.20 MB/s | 0.5423729 | 25506.41 MB/s |
| 256 | 0.5701559 | 26879.92 MB/s | 0.5412262 | 25503.95 MB/s |
| 512 | 0.5701559 | 26904.68 MB/s | 0.5423729 | 25584.65 MB/s |
| 1024 | 0.5704735 | 26911.46 MB/s | 0.5429480 | 25622.57 MB/s |
| 2048 | 0.5704735 | 26915.03 MB/s | 0.5433802 | 25640.09 MB/s |
| 4096 | 0.5703941 | 26917.72 MB/s | 0.5435965 | 25651.63 MB/s |
| 8192 | 0.5703544 | 26915.85 MB/s | 0.5436687 | 25656.76 MB/s |
| 16384 | 0.5699972 | 26898.44 MB/s | 0.5262583 | 24834.54 MB/s |
| 32768 | 0.5699873 | 26898.93 MB/s | 0.5262076 | 24833.21 MB/s |
| 65536 | 0.5698882 | 26894.48 MB/s | 0.5250567 | 24778.35 MB/s |
| 131072 | 0.5697024 | 26885.50 MB/s | 0.5224302 | 24654.59 MB/s |
| 262144 | 0.5696950 | 26884.72 MB/s | 0.5223095 | 24649.49 MB/s |
| 524288 | 0.5696937 | 26885.37 MB/s | 0.5223308 | 24650.21 MB/s |
| 1048576 | 0.5690340 | 26854.14 MB/s | 0.5220133 | 24634.71 MB/s |
| 2097152 | 0.5455717 | 25746.56 MB/s | 0.5041949 | 23794.65 MB/s |
| 4194304 | 0.5125461 | 24188.11 MB/s | 0.4756604 | 22447.05 MB/s |
| 8388608 | 0.5043430 | 23800.67 MB/s | 0.4659974 | 21991.51 MB/s |
| 16777216 | 0.5017375 | 23677.94 MB/s | 0.4614457 | 21776.58 MB/s |
| 33554432 | 0.5005865 | 23623.50 MB/s | 0.4596277 | 21690.63 MB/s |
--------------------------------------------------------------------------------
FMA + SSE4.1 确实给了我们一些改进，但这还不够好。
AVX2 + FMA 实现
最后，我们可以看到将 AVX2 plain division 与近似方法进行比较的真正好处：
===============================================================
= Compiler & System info
===============================================================
Current CPU : Intel(R) Xeon(R) CPU E3-1285L v3 @ 3.10GHz
CXX Compiler ID : GNU
CXX Compiler Path : /usr/bin/c++
CXX Compiler Version : 4.9.2
CXX Compiler Flags : -O3 -std=c++11 -march=haswell
--------------------------------------------------------------------------------
| Size | Division F/C | Division B/W | Approx. F/C | Approximation B/W |
--------------------------------------------------------------------------------
| 128 | 0.5663717 | 26672.73 MB/s | 0.9481481 | 44627.89 MB/s |
| 256 | 0.5651214 | 26653.72 MB/s | 0.9481481 | 44651.56 MB/s |
| 512 | 0.5644983 | 26640.36 MB/s | 0.9463956 | 44660.99 MB/s |
| 1024 | 0.5657459 | 26689.41 MB/s | 0.9552239 | 45044.21 MB/s |
| 2048 | 0.5662151 | 26715.40 MB/s | 0.9624060 | 45405.33 MB/s |
| 4096 | 0.5663717 | 26726.27 MB/s | 0.9671783 | 45633.64 MB/s |
| 8192 | 0.5664500 | 26732.42 MB/s | 0.9688941 | 45724.83 MB/s |
| 16384 | 0.5699377 | 26896.04 MB/s | 0.9092624 | 42909.11 MB/s |
| 32768 | 0.5699675 | 26897.85 MB/s | 0.9087077 | 42883.21 MB/s |
| 65536 | 0.5699625 | 26898.59 MB/s | 0.9001456 | 42480.91 MB/s |
| 131072 | 0.5699253 | 26896.38 MB/s | 0.8926057 | 42124.09 MB/s |
| 262144 | 0.5699117 | 26895.58 MB/s | 0.8928610 | 42137.13 MB/s |
| 524288 | 0.5698622 | 26892.87 MB/s | 0.8928002 | 42133.63 MB/s |
| 1048576 | 0.5685829 | 26833.13 MB/s | 0.8894302 | 41974.25 MB/s |
| 2097152 | 0.5558453 | 26231.90 MB/s | 0.8371921 | 39508.55 MB/s |
| 4194304 | 0.5224387 | 24654.67 MB/s | 0.7436747 | 35094.81 MB/s |
| 8388608 | 0.5143588 | 24273.46 MB/s | 0.7185252 | 33909.08 MB/s |
| 16777216 | 0.5107452 | 24103.19 MB/s | 0.7133449 | 33664.28 MB/s |
| 33554432 | 0.5101245 | 24074.10 MB/s | 0.7125114 | 33625.03 MB/s |
--------------------------------------------------------------------------------
结论
这种方法绝对可以提供对普通除法的加速。实际上可以达到多少加速，这实际上取决于底层架构，以及该部门如何与应用程序逻辑的其余部分交互。
【讨论】：
如果您正在检查 FMA 支持，您可能还应该检查 SSE4.1 并使用pmovzx 进行符号扩展，至少对于 128 位向量的低半部分。顺便说一句，SSE 移位计数饱和，所以我认为您可以使用 srai_epi16(a_epi16, 16) 然后 punpckl/h，因此将符号扩展到 32 位整数，并设置为使用 packs/usdw 重新打包。 （如果您想要截断而不是饱和，请先屏蔽。shift/and/OR 或 shift/pblendw 比 and/and/pack 更有效，但使用 shuffle 解包可能仍然是一个胜利。）
同意。我尝试做一个纯粹的SSE2 版本，而FMA 检查是以前基准测试的遗留问题。显然，在 ISA 方面往上走，可以比最初的代码做得更好。我可能应该编辑代码并引入带有基准的AVX2 版本，以获得完整的答案。
之前忘了说：如果除法是更大算法的一部分，divps is often better than rcpps + Newton iteration for overall throughput on modern CPUs。 （除了 KNL，您打算在其中使用 AVX512ER。）如果您在非完全流水线分隔器上没有瓶颈，则 128 位 divps 只是一个 uop，而不是几个。 rcpps+Newton 与 divps 具有相似的延迟。如果没有循环携带，OoO exec 无论如何都会隐藏延迟。
@PeterCordes - 感谢您的建议。我已经更新了我的答案，添加了将普通除法与近似步骤进行比较的性能配置文件。【参考方案2】：
请参阅 Agner Fog 的矢量类，他实现了一个快速算法，用 SSE/AVX 对 8 位、16 位和 32 位字（但不是 64 位）进行整数除法http://www.agner.org/optimize/#vectorclass
在文件 vectori128.h 中查找代码和算法的描述，作为他写得很好的手册 VectorClass.pdf
这是他手册中描述算法的片段。
"整数除法 x86 指令集中没有指令及其扩展 对整数向量除法很有用，如果它们这样的指令会很慢 存在。因此，向量类库使用的是快速整数算法 分配。该算法的基本原理可以用这个公式表示： a / b ≈ a * (2n / b) >> n 此计算经过以下步骤： 1. 为 n 找到一个合适的值 2.计算2n/b 3. 计算舍入误差的必要修正 4. 进行乘法和右移并应用修正以进行舍入 错误
如果多个数字除以相同的除数，这个公式是有利的 湾。步骤 1、2 和 3 只需执行一次，而对每个步骤重复步骤 4 股息的价值 a.数学细节在文件中描述 矢量i128.h。 （另见 T. Granlund 和 P. L. Montgomery：不变量除法 使用乘法的整数，SIGPLAN 的程序。”...
编辑：vectori128.h 文件末尾附近显示了如何使用标量变量进行短除法 “计算用于快速除法的参数比计算 做除法。因此，使用相同的除数对象是有利的 多次。例如，将 80 个无符号短整数除以 10：
short x = 10;
uint16_t dividends[80], quotients[80]; // numbers to work with
Divisor_us div10(x); // make divisor object for dividing by 10
Vec8us temp; // temporary vector
for (int i = 0; i < 80; i += 8) // loop for 4 elements per iteration
temp.load(dividends+i); // load 4 elements
temp /= div10; // divide each element by 10
temp.store(quotients+i); // store 4 elements
"
编辑：整数除以短裤向量
#include <stdio.h>
#include "vectorclass.h"
int main()
short numa[] = 10, 20, 30, 40, 50, 60, 70, 80;
short dena[] = 10, 20, 30, 40, 50, 60, 70, 80;
Vec8s num = Vec8s().load(numa);
Vec8s den = Vec8s().load(dena);
Vec4f num_low = to_float(extend_low(num));
Vec4f num_high = to_float(extend_high(num));
Vec4f den_low = to_float(extend_low(den));
Vec4f den_high = to_float(extend_high(den));
Vec4f qf_low = num_low/den_low;
Vec4f qf_high = num_high/den_high;
Vec4i q_low = truncate_to_int(qf_low);
Vec4i q_high = truncate_to_int(qf_high);
Vec8s q = compress(q_low, q_high);
for(int i=0; i<8; i++)
printf("%d ", q[i]);
printf("\n");
【讨论】：
这仅对除以常量非常有用-OP 说他想除以变量，在这种情况下，我认为唯一的选择是解包，转换为浮点数并执行浮点除法。 @PaulR，我添加了一个编辑以回应您的评论。在文件 vectori128.h 的末尾，他举例说明了如何使用变量对短裤进行除法。我还没有详细介绍，因为整数除法还不是我的代码中的瓶颈（我不得不使用它的时候我通过右移位来做到这一点，这非常快）。 @PaulR，所以我更仔细地研究了这一点。 VectorClass 中的快速方法仅适用于标量。对于向量/向量整数除法，我认为你是正确的。如果你必须做整数向量/向量，唯一的解决方案是转换为浮点数，进行除法，然后转换回短裤。我添加了一个编辑这样做。我想您也可以将短裤保存到数组中进行标量除法，然后也将它们加载回来。 是的，总结得差不多了——你可以在不使用浮点的情况下除以一个（单个）常量，但其他任何事情都需要转换为浮点然后再返回。 @PaulR，我想我们现在同意了。【参考方案3】：对于8位除法，可以通过创建幻数表来实现。
见"Hacker's Delight", page 238
签名：
__m128i _mm_div_epi8(__m128i a, __m128i b)
__m128i abs_b = _mm_abs_epi8(b);
static const uint16_t magic_number_table[129] =
0x0000, 0x0000, 0x8080, 0x5580, 0x4040, 0x3380, 0x2ac0, 0x24c0, 0x2020, 0x1c80, 0x19c0, 0x1760, 0x1560, 0x13c0, 0x1260, 0x1120,
0x1010, 0x0f20, 0x0e40, 0x0d80, 0x0ce0, 0x0c40, 0x0bb0, 0x0b30, 0x0ab0, 0x0a40, 0x09e0, 0x0980, 0x0930, 0x08e0, 0x0890, 0x0850,
0x0808, 0x07d0, 0x0790, 0x0758, 0x0720, 0x06f0, 0x06c0, 0x0698, 0x0670, 0x0640, 0x0620, 0x05f8, 0x05d8, 0x05b8, 0x0598, 0x0578,
0x0558, 0x0540, 0x0520, 0x0508, 0x04f0, 0x04d8, 0x04c0, 0x04b0, 0x0498, 0x0480, 0x0470, 0x0458, 0x0448, 0x0438, 0x0428, 0x0418,
0x0404, 0x03f8, 0x03e8, 0x03d8, 0x03c8, 0x03b8, 0x03ac, 0x03a0, 0x0390, 0x0388, 0x0378, 0x0370, 0x0360, 0x0358, 0x034c, 0x0340,
0x0338, 0x032c, 0x0320, 0x0318, 0x0310, 0x0308, 0x02fc, 0x02f4, 0x02ec, 0x02e4, 0x02dc, 0x02d4, 0x02cc, 0x02c4, 0x02bc, 0x02b4,
0x02ac, 0x02a8, 0x02a0, 0x0298, 0x0290, 0x028c, 0x0284, 0x0280, 0x0278, 0x0274, 0x026c, 0x0268, 0x0260, 0x025c, 0x0258, 0x0250,
0x024c, 0x0248, 0x0240, 0x023c, 0x0238, 0x0234, 0x022c, 0x0228, 0x0224, 0x0220, 0x021c, 0x0218, 0x0214, 0x0210, 0x020c, 0x0208,
0x0202
;
Uint8 load_den[16];
_mm_storeu_si128((__m128i*)load_den, abs_b);
uint16_t mul[16];
for (size_t i = 0; i < 16; i++)
uint16_t cur_den = load_den[i];
mul[i] = magic_number_table[cur_den];
// for denominator 1, magic number is 0x10080 that 16bit-overflow occurs.
__m128i one = _mm_set1_epi8(1);
__m128i is_one = _mm_cmpeq_epi8(abs_b, one);
// -128/-128 is a special case where magic number does not work.
__m128i v80 = _mm_set1_epi8(0x80);
__m128i is_80_a = _mm_cmpeq_epi8(a, v80);
__m128i is_80_b = _mm_cmpeq_epi8(b, v80);
__m128i is_80 = _mm_and_si128(is_80_a, is_80_b);
// __m128i zero = _mm_setzero_si128();
// __m128i less_a = _mm_cmpgt_epi8(zero, a);
// __m128i less_b = _mm_cmpgt_epi8(zero, b);
// __m128i sign = _mm_xor_si128(less_a, less_b);
__m128i abs_a = _mm_abs_epi8(a);
#if 0
__m128i p = _mm_unpacklo_epi8(abs_a, zero);
__m128i q = _mm_unpackhi_epi8(abs_a, zero);
__m256i c = _mm256_castsi128_si256(p);
c = _mm256_insertf128_si256(c, q, 1);
#else
// Thanks to Peter Cordes
__m256i c = _mm256_cvtepu8_epi16(abs_a);
#endif
__m256i magic = _mm256_loadu_si256((const __m256i*)mul);
__m256i high = _mm256_mulhi_epu16(magic, c);
__m128i v0h = _mm256_extractf128_si256(high, 0);
__m128i v0l = _mm256_extractf128_si256(high, 1);
__m128i res = _mm_packus_epi16(v0h, v0l);
__m128i div = _mm_blendv_epi8(res, abs_a, is_one);
// __m128i neg = _mm_sub_epi8(zero, div);
// __m128i select = _mm_blendv_epi8(div, neg, sign);
__m128i select = _mm_sign_epi8(div, _mm_or_si128(_mm_xor_si128(a, b), one));
return _mm_blendv_epi8(select, one, is_80);
无符号：
__m128i _mm_div_epu8(__m128i n, __m128i den)
static const uint16_t magic_number_table[256] =
0x0001, 0x0000, 0x8000, 0x5580, 0x4000, 0x3340, 0x2ac0, 0x04a0, 0x2000, 0x1c80, 0x19a0, 0x0750, 0x1560, 0x13c0, 0x0250, 0x1120,
0x1000, 0x0f10, 0x0e40, 0x0d80, 0x0cd0, 0x0438, 0x03a8, 0x0328, 0x0ab0, 0x0a40, 0x09e0, 0x0980, 0x0128, 0x00d8, 0x0890, 0x0048,
0x0800, 0x07c8, 0x0788, 0x0758, 0x0720, 0x06f0, 0x06c0, 0x0294, 0x0668, 0x0640, 0x021c, 0x05f8, 0x05d8, 0x01b4, 0x0194, 0x0578,
0x0558, 0x013c, 0x0520, 0x0508, 0x04f0, 0x04d8, 0x04c0, 0x04a8, 0x0094, 0x0480, 0x006c, 0x0458, 0x0448, 0x0034, 0x0024, 0x0014,
0x0400, 0x03f4, 0x03e4, 0x03d4, 0x03c8, 0x03b8, 0x03ac, 0x039c, 0x0390, 0x0384, 0x0378, 0x036c, 0x0360, 0x0354, 0x014a, 0x0340,
0x0334, 0x032c, 0x0320, 0x0318, 0x010e, 0x0304, 0x02fc, 0x02f4, 0x02ec, 0x02e4, 0x02dc, 0x02d4, 0x02cc, 0x02c4, 0x02bc, 0x02b4,
0x02ac, 0x02a4, 0x02a0, 0x0298, 0x0290, 0x028c, 0x0284, 0x007e, 0x0278, 0x0072, 0x026c, 0x0066, 0x0260, 0x025c, 0x0254, 0x0250,
0x004a, 0x0244, 0x0240, 0x023c, 0x0036, 0x0032, 0x022c, 0x0228, 0x0224, 0x001e, 0x001a, 0x0016, 0x0012, 0x000e, 0x000a, 0x0006,
0x0200, 0x00fd, 0x01fc, 0x01f8, 0x01f4, 0x01f0, 0x01ec, 0x01e8, 0x01e4, 0x01e0, 0x01dc, 0x01d8, 0x01d6, 0x01d4, 0x01d0, 0x01cc,
0x01c8, 0x01c4, 0x01c2, 0x01c0, 0x01bc, 0x01b8, 0x01b6, 0x01b4, 0x01b0, 0x01ae, 0x01ac, 0x01a8, 0x01a6, 0x01a4, 0x01a0, 0x019e,
0x019c, 0x0198, 0x0196, 0x0194, 0x0190, 0x018e, 0x018c, 0x018a, 0x0188, 0x0184, 0x0182, 0x0180, 0x017e, 0x017c, 0x017a, 0x0178,
0x0176, 0x0174, 0x0172, 0x0170, 0x016e, 0x016c, 0x016a, 0x0168, 0x0166, 0x0164, 0x0162, 0x0160, 0x015e, 0x015c, 0x015a, 0x0158,
0x0156, 0x0154, 0x0152, 0x0051, 0x0150, 0x014e, 0x014c, 0x014a, 0x0148, 0x0047, 0x0146, 0x0144, 0x0142, 0x0140, 0x003f, 0x013e,
0x013c, 0x013a, 0x0039, 0x0138, 0x0136, 0x0134, 0x0033, 0x0132, 0x0130, 0x002f, 0x012e, 0x012c, 0x012a, 0x0029, 0x0128, 0x0126,
0x0025, 0x0124, 0x0122, 0x0021, 0x0120, 0x001f, 0x011e, 0x011c, 0x001b, 0x011a, 0x0019, 0x0118, 0x0116, 0x0015, 0x0114, 0x0013,
0x0112, 0x0110, 0x000f, 0x010e, 0x000d, 0x010c, 0x000b, 0x010a, 0x0009, 0x0108, 0x0007, 0x0106, 0x0005, 0x0104, 0x0003, 0x0102
;
static const uint16_t shift_table[256] =
0x0001, 0x0100, 0x0100, 0x0080, 0x0100, 0x0040, 0x0040, 0x0020, 0x0100, 0x0080, 0x0020, 0x0010, 0x0020, 0x0040, 0x0010, 0x0020,
0x0100, 0x0010, 0x0040, 0x0080, 0x0010, 0x0008, 0x0008, 0x0008, 0x0010, 0x0040, 0x0020, 0x0080, 0x0008, 0x0008, 0x0010, 0x0008,
0x0100, 0x0008, 0x0008, 0x0008, 0x0020, 0x0010, 0x0040, 0x0004, 0x0008, 0x0040, 0x0004, 0x0008, 0x0008, 0x0004, 0x0004, 0x0008,
0x0008, 0x0004, 0x0020, 0x0008, 0x0010, 0x0008, 0x0040, 0x0008, 0x0004, 0x0080, 0x0004, 0x0008, 0x0008, 0x0004, 0x0004, 0x0004,
0x0100, 0x0004, 0x0004, 0x0004, 0x0008, 0x0008, 0x0004, 0x0004, 0x0010, 0x0004, 0x0008, 0x0004, 0x0020, 0x0004, 0x0002, 0x0040,
0x0004, 0x0004, 0x0020, 0x0008, 0x0002, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004,
0x0004, 0x0004, 0x0020, 0x0008, 0x0010, 0x0004, 0x0004, 0x0002, 0x0008, 0x0002, 0x0004, 0x0002, 0x0020, 0x0004, 0x0004, 0x0010,
0x0002, 0x0004, 0x0040, 0x0004, 0x0002, 0x0002, 0x0004, 0x0008, 0x0004, 0x0002, 0x0002, 0x0002, 0x0002, 0x0002, 0x0002, 0x0002,
0x0100, 0x0001, 0x0004, 0x0008, 0x0004, 0x0010, 0x0004, 0x0008, 0x0004, 0x0020, 0x0004, 0x0008, 0x0002, 0x0004, 0x0010, 0x0004,
0x0008, 0x0004, 0x0002, 0x0040, 0x0004, 0x0008, 0x0002, 0x0004, 0x0010, 0x0002, 0x0004, 0x0008, 0x0002, 0x0004, 0x0020, 0x0002,
0x0004, 0x0008, 0x0002, 0x0004, 0x0010, 0x0002, 0x0004, 0x0002, 0x0008, 0x0004, 0x0002, 0x0080, 0x0002, 0x0004, 0x0002, 0x0008,
0x0002, 0x0004, 0x0002, 0x0010, 0x0002, 0x0004, 0x0002, 0x0008, 0x0002, 0x0004, 0x0002, 0x0020, 0x0002, 0x0004, 0x0002, 0x0008,
0x0002, 0x0004, 0x0002, 0x0001, 0x0010, 0x0002, 0x0004, 0x0002, 0x0008, 0x0001, 0x0002, 0x0004, 0x0002, 0x0040, 0x0001, 0x0002,
0x0004, 0x0002, 0x0001, 0x0008, 0x0002, 0x0004, 0x0001, 0x0002, 0x0010, 0x0001, 0x0002, 0x0004, 0x0002, 0x0001, 0x0008, 0x0002,
0x0001, 0x0004, 0x0002, 0x0001, 0x0020, 0x0001, 0x0002, 0x0004, 0x0001, 0x0002, 0x0001, 0x0008, 0x0002, 0x0001, 0x0004, 0x0001,
0x0002, 0x0010, 0x0001, 0x0002, 0x0001, 0x0004, 0x0001, 0x0002, 0x0001, 0x0008, 0x0001, 0x0002, 0x0001, 0x0004, 0x0001, 0x0002
;
static const uint16_t mask_table[256] =
0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0x0000, 0xffff,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0xffff, 0x0000,
0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0xffff, 0xffff,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000,
0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff,
0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000,
0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000,
0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff,
0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000
;
uint8_t load_den[16];
_mm_storeu_si128((__m128i*)load_den, den);
uint16_t mul[16];
uint16_t mask[16];
uint16_t shift[16];
for (size_t i = 0; i < 16; i++)
const uint16_t cur_den = load_den[i];
mul[i] = magic_number_table[cur_den];
mask[i] = mask_table[cur_den];
shift[i] = shift_table[cur_den];
#if 0
__m128i a = _mm_unpacklo_epi8(n, _mm_setzero_si128());
__m128i b = _mm_unpackhi_epi8(n, _mm_setzero_si128());
__m256i c = _mm256_castsi128_si256(a);
c = _mm256_insertf128_si256(c, b, 1);
#else
// Thanks to Peter Cordes
__m256i c = _mm256_cvtepu8_epi16(n);
#endif
__m256i magic = _mm256_loadu_si256((const __m256i*)mul);
__m256i high = _mm256_mulhi_epu16(magic, c);
__m256i low = _mm256_mullo_epi16(magic, c);
__m256i low_down = _mm256_srli_epi16(low, 8);
__m256i high_up = _mm256_slli_epi16(high, 8);
__m256i low_high = _mm256_or_si256(low_down, high_up);
__m256i target_up = _mm256_mullo_epi16(c, _mm256_loadu_si256((const __m256i*)shift));
__m256i cal1 = _mm256_sub_epi16(target_up, low_high);
__m256i cal2 = _mm256_srli_epi16(cal1, 1);
__m256i cal3 = _mm256_add_epi16(cal2, low_high);
__m256i cal4 = _mm256_srli_epi16(cal3, 7);
__m256i res = _mm256_blendv_epi8(high, cal4, _mm256_loadu_si256((const __m256i*)mask));
__m128i v0h = _mm256_extractf128_si256(res, 0);
__m128i v0l = _mm256_extractf128_si256(res, 1);
return _mm_packus_epi16(v0h, v0l);
【讨论】：
如果你打算使用 AVX2，你可以并且应该使用 VPMOVZXBW 用一条指令创建c，而不是单独解包 lo/hi 并将它们混在一起。 _mm256_cvtepu8_epi16
还可以优化符号计算：由于blendv只需要高位正确，你可以__m128i sign = _mm_xor_si128(a,b);。或者更好的是，使用vpsignb 有条件地否定。除了a^b 对a==b 来说为零，结果应该是+1 而不是0。所以(a^b) | 1 会起作用，因为你已经有一个set1(1) 向量，它会使其非零而不影响符号少量。所以__m128i select = _mm_sign_epi8(div, a^b | one);。 （vblendv 是 2 uops，只在 Haswell 的 shuffle 端口上运行）
@Peter 谢谢你的建议。【参考方案4】：
感谢原始解决方案，这适用于新手： 这个Agner Fog's subroutine library 在优化方面对我施展了魔法
这是你多次除以同一个变量值的情况（比如在大循环中）
#include <asmlib.h>
unsigned int a, b, d;
unsigned int divisor = any_random_value;
div_u32 OptimumDivision(divisor);
a/OptimumDivision;
b/OptimumDivision;
这是 unsigned int - 如果您需要负值，请使用div_i32，而不是在我的测试中更快，即使手册说相反
我的性能大约是 3 倍或更多
【讨论】：
以上是关于SSE整数除法？的主要内容，如果未能解决你的问题，请参考以下文章