如果你的目标系统支持SSE指令集，那么使用它可以获得良好的性能提升。我已经测试过其他发布的方法并发现这是最快的方法。
大量数值比较的示例结果如下：
diff0: 188.020679 ms // branching
diff1: 118.934970 ms // max min
diff2: 97.087710 ms // branchless mul add
diff3: 54.495269 ms // branchless signed
diff4: 31.159628 ms // sse
diff5: 30.855885 ms // sse v2
以下是我的完整测试代码。我使用了SSE2指令，这些指令现在在x86ish CPU中广泛使用，通过SSE内置函数，应该相当可移植（MSVC、GCC、Clang、Intel编译器等）。
注意事项：
- 实际上，这个计算方法是先计算最大值，然后计算最小值，最后减去两者之差，但每次执行16个值的计算，每个指令都会执行。
- 在
diff5中展开循环貌似没有太大效果，但可能可以进行调整。
- 对于最后15个或更少的值，当前的回退方法是在一个循环中使用带符号的技巧方法，但可能可以通过展开和/或SSE进一步加速。
- 这些函数本身非常简单，因此它们应该很容易移植到任何具有SSE内置函数或asm的东西上。
- 我使用了Windows特定的计时函数，因为
std::chrono::high_resolution_clock在MSVC实现中精度较低，对于C/C++测试代码的混合使用感到抱歉。
- 在测量性能之后，将根据参考分支实现测试结果，因此它们应该是正确的。
如果对代码或这种方法有任何问题/建议，请留下评论。
#include <cstdlib>
#include <cstdint>
#include <cstdio>
#include <cmath>
#include <random>
#include <algorithm>
#define WIN32_LEAN_AND_MEAN
#define NOMINMAX
#include <Windows.h>
#include <emmintrin.h>
void diff0(const std::uint8_t* a, const std::uint8_t* b, std::uint8_t* res,
std::size_t n)
{
for (std::size_t i = 0; i < n; i++) {
res[i] = a[i] > b[i] ? a[i] - b[i] : b[i] - a[i];
}
}
void diff1(const std::uint8_t* a, const std::uint8_t* b, std::uint8_t* res,
std::size_t n)
{
for (std::size_t i = 0; i < n; i++) {
res[i] = std::max(a[i], b[i]) - std::min(a[i], b[i]);
}
}
void diff2(const std::uint8_t* a, const std::uint8_t* b, std::uint8_t* res,
std::size_t n)
{
for (std::size_t i = 0; i < n; i++) {
res[i] = (a[i] > b[i]) * (a[i] - b[i]) + (a[i] < b[i]) * (b[i] - a[i]);
}
}
void diff3(const std::uint8_t* a, const std::uint8_t* b, std::uint8_t* res,
std::size_t n)
{
for (std::size_t i = 0; i < n; i++) {
std::int16_t diff = a[i] - b[i];
std::uint16_t mask = diff >> 15;
res[i] = (diff + mask) ^ mask;
}
}
void diff4(const std::uint8_t* a, const std::uint8_t* b, std::uint8_t* res,
std::size_t n)
{
auto pA = reinterpret_cast<const __m128i*>(a);
auto pB = reinterpret_cast<const __m128i*>(b);
auto pRes = reinterpret_cast<__m128i*>(res);
std::size_t i = 0;
for (std::size_t j = n / 16; j--; i++) {
__m128i max = _mm_max_epu8(_mm_load_si128(pA + i), _mm_load_si128(pB + i));
__m128i min = _mm_min_epu8(_mm_load_si128(pA + i), _mm_load_si128(pB + i));
_mm_store_si128(pRes + i, _mm_sub_epi8(max, min));
}
for (i *= 16; i < n; i++) {
std::int16_t diff = a[i] - b[i];
std::uint16_t mask = diff >> 15;
res[i] = (diff + mask) ^ mask;
}
}
void diff5(const std::uint8_t* a, const std::uint8_t* b, std::uint8_t* res,
std::size_t n)
{
auto pA = reinterpret_cast<const __m128i*>(a);
auto pB = reinterpret_cast<const __m128i*>(b);
auto pRes = reinterpret_cast<__m128i*>(res);
std::size_t i = 0;
const std::size_t UNROLL = 2;
for (std::size_t j = n / (16 * UNROLL); j--; i += UNROLL) {
__m128i max0 = _mm_max_epu8(_mm_load_si128(pA + i + 0), _mm_load_si128(pB + i + 0));
__m128i min0 = _mm_min_epu8(_mm_load_si128(pA + i + 0), _mm_load_si128(pB + i + 0));
__m128i max1 = _mm_max_epu8(_mm_load_si128(pA + i + 1), _mm_load_si128(pB + i + 1));
__m128i min1 = _mm_min_epu8(_mm_load_si128(pA + i + 1), _mm_load_si128(pB + i + 1));
_mm_store_si128(pRes + i + 0, _mm_sub_epi8(max0, min0));
_mm_store_si128(pRes + i + 1, _mm_sub_epi8(max1, min1));
}
for (std::size_t j = n % (16 * UNROLL) / 16; j--; i++) {
__m128i max = _mm_max_epu8(_mm_load_si128(pA + i), _mm_load_si128(pB + i));
__m128i min = _mm_min_epu8(_mm_load_si128(pA + i), _mm_load_si128(pB + i));
_mm_store_si128(pRes + i, _mm_sub_epi8(max, min));
}
for (i *= 16; i < n; i++) {
std::int16_t diff = a[i] - b[i];
std::uint16_t mask = diff >> 15;
res[i] = (diff + mask) ^ mask;
}
}
int main() {
const std::size_t ALIGN = 16;
const std::size_t N = 10 * 1024 * 1024 * 3;
auto a = static_cast<uint8_t*>(_mm_malloc(N, ALIGN));
auto b = static_cast<uint8_t*>(_mm_malloc(N, ALIGN));
{
std::mt19937 engine(std::random_device{}());
std::uniform_int<std::uint8_t> distribution(0, 255);
for (std::size_t i = 0; i < N; i++) {
a[i] = distribution(engine);
b[i] = distribution(engine);
}
}
auto res0 = static_cast<uint8_t*>(_mm_malloc(N, ALIGN));
auto resX = static_cast<uint8_t*>(_mm_malloc(N, ALIGN));
LARGE_INTEGER f, t0, t1;
QueryPerformanceFrequency(&f);
QueryPerformanceCounter(&t0);
diff0(a, b, res0, N);
QueryPerformanceCounter(&t1);
printf("diff0: %.6f ms\n",
static_cast<double>(t1.QuadPart - t0.QuadPart) / f.QuadPart * 1000);
#define TEST(diffX)\
QueryPerformanceCounter(&t0);\
diffX(a, b, resX, N);\
QueryPerformanceCounter(&t1);\
printf("%s: %.6f ms\n", #diffX,\
static_cast<double>(t1.QuadPart - t0.QuadPart) / f.QuadPart * 1000);\
for (std::size_t i = 0; i < N; i++) {\
if (resX[i] != res0[i]) {\
printf("error: %s(%03u, %03u) == %03u != %03u\n", #diffX,\
a[i], b[i], resX[i], res0[i]);\
break;\
}\
}
TEST(diff1);
TEST(diff2);
TEST(diff3);
TEST(diff4);
TEST(diff5);
_mm_free(a);
_mm_free(b);
_mm_free(res0);
_mm_free(resX);
getc(stdin);
return 0;
}
(a - b)*(a - b)- chux - Reinstate Monica(a-b)*(a-b) = (b-a)*(b-a)。 - Niklas B.abs函数？你只有在处理最大规范之类的问题时才需要使用abs。 - nwellnhofpsadbw指令，可为8个Word8SAD操作提供总和。 因此，如果您将2个输入字节零扩展到XMM寄存器中，则psadbw会执行所需的操作。 它专门为视频编解码器执行多个像素块运动搜索而设计，但对于您的用例，您可以使用SSE4.1的pmovzxbq将2个字节加载到2个qword中以并行检查2个像素组件。 同时使用pmuludq来平方这2个结果。 我不知道如何让Haskell编译器发出它，我完全不懂Haskell。 - Peter Cordes