【发布时间】：2021-02-10 05:49:01
【问题描述】：
有一个existing question“3 个长整数的平均值”专门关注三个有符号整数的平均值的有效计算。
但是，使用无符号整数可以进行其他优化，但不适用于上一个问题中涉及的场景。这个问题是关于三个 unsigned 整数平均值的有效计算，其中平均值向零舍入，即在数学术语中我想计算 ⌊ (a + b + c) / 3 ⌋。
计算这个平均值的简单方法是
avg = a / 3 + b / 3 + c / 3 + (a % 3 + b % 3 + c % 3) / 3;
首先，现代优化编译器会将除法转换为具有倒数加移位的乘法，并将模运算转换为反向乘法和减法，其中反向乘法可能使用 scale_add em> 成语可用于许多架构，例如lea 在 x86_64 上，add 和 lsl #n 在 ARM 上，iscadd 在 NVIDIA GPU 上。
在尝试以适合许多常见平台的通用方式优化上述内容时，我观察到整数运算的成本通常在关系 logical ≤ (add | sub) ≤ shift ≤ scale_add ≤ mul。这里的成本是指所有延迟、吞吐量限制和功耗。当处理的整数类型比本机寄存器宽度宽时，任何此类差异都会变得更加明显，例如在 32 位处理器上处理 uint64_t 数据时。
因此，我的优化策略是尽量减少指令数，并在可能的情况下用“廉价”操作替换“昂贵”操作，同时不增加寄存器压力并为广泛的乱序处理器保留可利用的并行性。
第一个观察结果是，我们可以通过首先应用产生一个和值和一个进位值的 CSA（进位保存加法器）将三个操作数的总和减少为两个操作数的总和，其中进位值的权重是两倍的总和值。在大多数处理器上，基于软件的 CSA 的成本是 5 个逻辑。一些处理器，比如 NVIDIA GPU，有一个 LOP3 指令，可以一举计算三个操作数的任意逻辑表达式，在这种情况下，CSA 会压缩为两个 LOP3s（注意：我还没有说服 CUDA 编译器发出这两个LOP3s；它目前产生四个LOP3s！）。
第二个观察结果是，因为我们正在计算除以 3 的模数，所以我们不需要反向乘法来计算它。我们可以改用 dividend % 3 = ((dividend / 3) + dividend) & 3，将模数减少为 add 加上 logical，因为我们已经有了除法结果。这是通用算法的一个实例：股息 % (2n-1) = ((dividend / (2n-1) + 股息) & (2n-1).
最后对于修正项(a % 3 + b % 3 + c % 3) / 3中的除以3，我们不需要通用除以3的代码。由于除数很小，在[0, 6]中，我们可以将x / 3简化为(3 * x) / 8 只需要一个 scale_add 加上一个 shift。
下面的代码显示了我当前正在进行的工作。使用 Compiler Explorer 检查为各种平台生成的代码显示了我所期望的紧凑代码（使用 -O3 编译时）。
但是，在使用 Intel 13.x 编译器对我的 Ivy Bridge x86_64 机器上的代码进行计时时，一个缺陷变得明显：与简单的代码相比，我的代码提高了延迟（uint64_t 数据从 18 个周期到 15 个周期）版本，吞吐量恶化（uint64_t 数据从每 6.8 个周期一个结果到每 8.5 个周期一个结果）。更仔细地查看汇编代码，很明显为什么会这样：我基本上设法将代码从大致三向并行度降低到大致双向并行度。
是否有一种普遍适用的优化技术，对常见处理器特别是所有类型的 x86 和 ARM 以及 GPU 都有益，并且可以保留更多的并行性？或者，是否有一种优化技术可以进一步减少总体操作数以弥补并行度的降低？校正项的计算（下面代码中的tail）似乎是一个很好的目标。简化 (carry_mod_3 + sum_mod_3) / 2 看起来很诱人，但对于九种可能的组合之一给出的结果不正确。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#define BENCHMARK (1)
#define SIMPLE_COMPUTATION (0)
#if BENCHMARK
#define T uint64_t
#else // !BENCHMARK
#define T uint8_t
#endif // BENCHMARK
T average_of_3 (T a, T b, T c)
{
T avg;
#if SIMPLE_COMPUTATION
avg = a / 3 + b / 3 + c / 3 + (a % 3 + b % 3 + c % 3) / 3;
#else // !SIMPLE_COMPUTATION
/* carry save adder */
T a_xor_b = a ^ b;
T sum = a_xor_b ^ c;
T carry = (a_xor_b & c) | (a & b);
/* here 2 * carry + sum = a + b + c */
T sum_div_3 = (sum / 3); // {MUL|MULHI}, SHR
T sum_mod_3 = (sum + sum_div_3) & 3; // ADD, AND
if (sizeof (size_t) == sizeof (T)) { // "native precision" (well, not always)
T two_carry_div_3 = (carry / 3) * 2; // MULHI, ANDN
T two_carry_mod_3 = (2 * carry + two_carry_div_3) & 6; // SCALE_ADD, AND
T head = two_carry_div_3 + sum_div_3; // ADD
T tail = (3 * (two_carry_mod_3 + sum_mod_3)) / 8; // ADD, SCALE_ADD, SHR
avg = head + tail; // ADD
} else {
T carry_div_3 = (carry / 3); // MUL, SHR
T carry_mod_3 = (carry + carry_div_3) & 3; // ADD, AND
T head = (2 * carry_div_3 + sum_div_3); // SCALE_ADD
T tail = (3 * (2 * carry_mod_3 + sum_mod_3)) / 8; // SCALE_ADD, SCALE_ADD, SHR
avg = head + tail; // ADD
}
#endif // SIMPLE_COMPUTATION
return avg;
}
#if !BENCHMARK
/* Test correctness on 8-bit data exhaustively. Should catch most errors */
int main (void)
{
T a, b, c, res, ref;
a = 0;
do {
b = 0;
do {
c = 0;
do {
res = average_of_3 (a, b, c);
ref = ((uint64_t)a + (uint64_t)b + (uint64_t)c) / 3;
if (res != ref) {
printf ("a=%08x b=%08x c=%08x res=%08x ref=%08x\n",
a, b, c, res, ref);
return EXIT_FAILURE;
}
c++;
} while (c);
b++;
} while (b);
a++;
} while (a);
return EXIT_SUCCESS;
}
#else // BENCHMARK
#include <math.h>
// A routine to give access to a high precision timer on most systems.
#if defined(_WIN32)
#if !defined(WIN32_LEAN_AND_MEAN)
#define WIN32_LEAN_AND_MEAN
#endif
#include <windows.h>
double second (void)
{
LARGE_INTEGER t;
static double oofreq;
static int checkedForHighResTimer;
static BOOL hasHighResTimer;
if (!checkedForHighResTimer) {
hasHighResTimer = QueryPerformanceFrequency (&t);
oofreq = 1.0 / (double)t.QuadPart;
checkedForHighResTimer = 1;
}
if (hasHighResTimer) {
QueryPerformanceCounter (&t);
return (double)t.QuadPart * oofreq;
} else {
return (double)GetTickCount() * 1.0e-3;
}
}
#elif defined(__linux__) || defined(__APPLE__)
#include <stddef.h>
#include <sys/time.h>
double second (void)
{
struct timeval tv;
gettimeofday(&tv, NULL);
return (double)tv.tv_sec + (double)tv.tv_usec * 1.0e-6;
}
#else
#error unsupported platform
#endif
#define N (3000000)
int main (void)
{
double start, stop, elapsed = INFINITY;
int i, k;
T a, b;
T avg0 = 0xffffffff, avg1 = 0xfffffffe;
T avg2 = 0xfffffffd, avg3 = 0xfffffffc;
T avg4 = 0xfffffffb, avg5 = 0xfffffffa;
T avg6 = 0xfffffff9, avg7 = 0xfffffff8;
T avg8 = 0xfffffff7, avg9 = 0xfffffff6;
T avg10 = 0xfffffff5, avg11 = 0xfffffff4;
T avg12 = 0xfffffff2, avg13 = 0xfffffff2;
T avg14 = 0xfffffff1, avg15 = 0xfffffff0;
a = 0x31415926;
b = 0x27182818;
avg0 = average_of_3 (a, b, avg0);
for (k = 0; k < 5; k++) {
start = second();
for (i = 0; i < N; i++) {
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
avg0 = average_of_3 (a, b, avg0);
b = (b + avg0) ^ a;
a = (a ^ b) + avg0;
}
stop = second();
elapsed = fmin (stop - start, elapsed);
}
printf ("a=%016llx b=%016llx avg=%016llx",
(uint64_t)a, (uint64_t)b, (uint64_t)avg0);
printf ("\rlatency: each average_of_3() took %.6e seconds\n",
elapsed / 16 / N);
a = 0x31415926;
b = 0x27182818;
avg0 = average_of_3 (a, b, avg0);
for (k = 0; k < 5; k++) {
start = second();
for (i = 0; i < N; i++) {
avg0 = average_of_3 (a, b, avg0);
avg1 = average_of_3 (a, b, avg1);
avg2 = average_of_3 (a, b, avg2);
avg3 = average_of_3 (a, b, avg3);
avg4 = average_of_3 (a, b, avg4);
avg5 = average_of_3 (a, b, avg5);
avg6 = average_of_3 (a, b, avg6);
avg7 = average_of_3 (a, b, avg7);
avg8 = average_of_3 (a, b, avg8);
avg9 = average_of_3 (a, b, avg9);
avg10 = average_of_3 (a, b, avg10);
avg11 = average_of_3 (a, b, avg11);
avg12 = average_of_3 (a, b, avg12);
avg13 = average_of_3 (a, b, avg13);
avg14 = average_of_3 (a, b, avg14);
avg15 = average_of_3 (a, b, avg15);
b = (b + avg0) ^ a;
a = (a ^ b) + avg0;
}
stop = second();
elapsed = fmin (stop - start, elapsed);
}
printf ("a=%016llx b=%016llx avg=%016llx", (uint64_t)a, (uint64_t)b,
(uint64_t)(avg0 + avg1 + avg2 + avg3 + avg4 + avg5 + avg6 + avg7 +
avg8 + avg9 +avg10 +avg11 +avg12 +avg13 +avg14 +avg15));
printf ("\rthroughput: each average_of_3() took %.6e seconds\n",
elapsed / 16 / N);
return EXIT_SUCCESS;
}
#endif // BENCHMARK
【问题讨论】：
-
我不认为你的
carry = (a_xor_b & c) | (a & b)赶上a&b&c。为什么 (4/3sum)%4 等于sum%3？ -
最好将这些测试程序编辑到问题中；它们足够长，以至于在 SO cmets 中变得一团糟。也许链接一个简短的版本来测试每个输入的低 3 位及其 MSB 的所有组合，可能在 tio.run 或 godbolt.org 上，因为现在 Godbolt 有一个编译和运行选项。
-
输入的分布是什么？
-
@Sopel：无法对此做出任何假设，因此假设均匀随机。因此非常需要无分支解决方案。
-
@njuffa IMO 值得针对特定架构和数据大小对其进行优化。如果我们考虑 64 位系统和 32 位数据（在大多数 64 位系统上整数是 32 位），“天真”方法是最有效的。 godbolt.org/z/dh1jrPIMO 如果我们正在考虑微优化和 uOptimizations，它不能从硬件和实现中抽象出来。
标签： c algorithm bit-manipulation micro-optimization extended-precision