nju*_*ffa 21 c algorithm assembly bit-manipulation micro-optimization
我试图找到最有效的方法来计算 32 位无符号整数的模 255。我的主要重点是找到一种可以在 x86 和 ARM 平台上运行良好的算法，并着眼于除此之外的适用性。首先，我试图避免内存操作（这可能很昂贵），所以我正在寻找有点复杂的方法，同时避免使用表格。我还试图避免可能昂贵的操作，例如分支和乘法，并尽量减少使用的操作和寄存器的数量。
下面的 ISO-C99 代码捕获了我迄今为止尝试过的八个变体。它包括一个用于详尽测试的框架。我对这个粗略的执行时间测量进行了猛烈抨击，这似乎工作得很好，可以获得第一次性能印象。在一些平台上我试过（全部具有快速整数倍）的变种WARREN_MUL_SHR_2，WARREN_MUL_SHR_1和DIGIT_SUM_CARRY_OUT_1似乎是最高效的。我的实验表明，我在Compiler Explorer 中尝试的 x86、ARM、PowerPC 和 MIPS 编译器都很好地利用了特定于平台的功能，例如三输入LEA、字节扩展指令、乘法累加和指令预测。
该变体NAIVE_USING_DIV使用整数除法，与除数反乘，然后减法。这是基线情况。现代编译器知道如何有效地实现 255 的无符号整数除法（通过乘法），并将在适当的情况下使用离散替换反乘。要计算模数，base-1可以对base数字求和，然后折叠结果。例如3334 mod 9: sum 3+3+3+4 = 13, fold 1+3 = 4. 如果折叠后的结果是base-1，我们需要生成0来代替。DIGIT_SUM_THEN_FOLD使用这种方法。
A. Cockburn，“使用 8/16 位算法有效实现 OSI 传输协议校验和算法”，ACM SIGCOMM 计算机通信评论，卷。17, No. 3, 七月/八月 1987 年，第 13-20 页
展示了base-1在校验和计算模 255 的上下文中有效地添加数字模数的不同方法。计算数字的逐字节总和，并且在每次添加之后，也添加来自加法的任何进位。所以这将是一个ADD a, b,ADC a, 0序列。使用base 256数字为此写出加法链，很明显，计算基本上是乘以0x0101 ... 0101。结果将在最高有效数字位置，除了需要单独捕获该位置的加法的进位。该方法仅在一个base数字包含2k位时有效。在这里，我们有k=3. 我尝试了三种不同的方法将 的结果重新映射base-1到 0，从而产生了变体DIGIT_SUM_CARRY_OUT_1, DIGIT_SUM_CARRY_OUT_2,DIGIT_SUM_CARRY_OUT_3.
1995 年 7 月 9 日，Joe Keane 在新闻组 comp.lang.c 中演示了一种有效计算模 63 的有趣方法。虽然线程参与者 Peter L. Montgomery 证明了该算法是正确的，但不幸的是 Keane 先生没有回应解释其推导的请求。该算法也在 H. Warren 的Hacker's Delight 2nd ed 中重现。我能够以纯粹的机械方式将它扩展到模 127 和模 255。这是（适当命名的）KEANE_MAGIC 变体。更新：自从我最初发布这个问题以来，我发现基恩的方法基本上是以下内容的巧妙定点实现：return (uint32_t)(fmod (x * 256.0 / 255.0 + 0.5, 256.0) * (255.0 / 256.0));. 这使它成为下一个变体的近亲。
亨利 S. 沃伦，《黑客的喜悦》第 2 版。，第 272 显示了一种“乘法右移”算法，大概是由作者自己设计的，它基于数学属性 n mod 2 k-1 = floor (2 k / 2 k-1 * n) mod 2 k。定点计算用于乘以因子 2 k / 2 k-1。我构建了两个变体，它们的不同之处在于它们如何处理从base-10的初步结果的映射。它们是变体WARREN_MUL_SHR_1和WARREN_MUL_SHR_2。
是否存在比我目前确定的三个顶级竞争者更有效的模 255 计算算法，尤其是对于整数乘法较慢的平台？base 256在这种情况下，对基恩的四位数字求和的无乘法算法的有效修改似乎特别有趣。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#define NAIVE_USING_DIV (1)
#define DIGIT_SUM_THEN_FOLD (2)
#define DIGIT_SUM_CARRY_OUT_1 (3)
#define DIGIT_SUM_CARRY_OUT_2 (4)
#define DIGIT_SUM_CARRY_OUT_3 (5)
#define KEANE_MAGIC (6) // Joe Keane, comp.lang.c, 1995/07/09
#define WARREN_MUL_SHR_1 (7) // Hacker's Delight, 2nd ed., p. 272
#define WARREN_MUL_SHR_2 (8) // Hacker's Delight, 2nd ed., p. 272
#define VARIANT (WARREN_MUL_SHR_2)
uint32_t mod255 (uint32_t x)
{
#if VARIANT == NAIVE_USING_DIV
return x - 255 * (x / 255);
#elif VARIANT == DIGIT_SUM_THEN_FOLD
x = (x & 0xffff) + (x >> 16);
x = (x & 0xff) + (x >> 8);
x = (x & 0xff) + (x >> 8) + 1;
x = (x & 0xff) + (x >> 8) - 1;
return x;
#elif VARIANT == DIGIT_SUM_CARRY_OUT_1
uint32_t t;
t = 0x01010101 * x;
t = (t >> 24) + (t < x);
if (t == 255) t = 0;
return t;
#elif VARIANT == DIGIT_SUM_CARRY_OUT_2
uint32_t t;
t = 0x01010101 * x;
t = (t >> 24) + (t < x) + 1;
t = (t & 0xff) + (t >> 8) - 1;
return t;
#elif VARIANT == DIGIT_SUM_CARRY_OUT_3
uint32_t t;
t = 0x01010101 * x;
t = (t >> 24) + (t < x);
t = t & ((t - 255) >> 8);
return t;
#elif VARIANT == KEANE_MAGIC
x = (((x >> 16) + x) >> 14) + (x << 2);
x = ((x >> 8) + x + 2) & 0x3ff;
x = (x - (x >> 8)) >> 2;
return x;
#elif VARIANT == WARREN_MUL_SHR_1
x = (0x01010101 * x + (x >> 8)) >> 24;
x = x & ((x - 255) >> 8);
return x;
#elif VARIANT == WARREN_MUL_SHR_2
x = (0x01010101 * x + (x >> 8)) >> 24;
if (x == 255) x = 0;
return x;
#else
#error unknown VARIANT
#endif
}
uint32_t ref_mod255 (uint32_t x)
{
volatile uint32_t t = x;
t = t % 255;
return t;
}
// timing with microsecond resolution
#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
int main (void)
{
double start, stop;
uint32_t res, ref, x = 0;
printf ("Testing VARIANT = %d\n", VARIANT);
start = second();
do {
res = mod255 (x);
ref = ref_mod255 (x);
if (res != ref) {
printf ("error @ %08x: res=%08x ref=%08x\n", x, res, ref);
return EXIT_FAILURE;
}
x++;
} while (x);
stop = second();
printf ("test passed\n");
printf ("elapsed = %.6f seconds\n", stop - start);
return EXIT_SUCCESS;
}
Adr*_*ica 14
对于任意无符号整数x和n，求模表达式x % n涉及（至少在概念上）三个运算：除法、乘法和减法：
quotient = x / n;
product = quotient * n;
modulus = x - product;
然而，当Ñ为2（的功率N = 2 p），模可确定多更加迅速简单地通过屏蔽掉所有，但越低，p位。
在大多数 CPU 上，加法、减法和位屏蔽是非常“便宜”（快速）的运算，乘法更“昂贵”而除法非常昂贵——但请注意，大多数优化编译器会将除法转换为编译时常量乘法（通过不同的常数）和位移（见下文）。
因此，如果我们可以将模 255 转换为模 256，而无需太多开销，我们可能会加快该过程。我们可以通过注意到它x % n等价于(x + x / n) % (n + 1)†来做到这一点。因此，我们现在的概念操作是：除法、加法和屏蔽。
在屏蔽低 8 位的特定情况下，基于 x86/x64 的 CPU（和其他？）可能能够执行进一步优化，因为它们可以访问（大多数）寄存器的 8 位版本。
以下是 clang-cl 编译器为一个 naïve modulo 255 函数（传入ecx并返回的参数eax）生成的内容：
unsigned Naive255(unsigned x)
{
return x % 255;
}
mov edx, ecx
mov eax, 2155905153 ;
imul rax, rdx ; Replacing the IDIV with IMUL and SHR
shr rax, 39 ;
mov edx, eax
shl edx, 8
sub eax, edx
add eax, ecx
这是使用上述“技巧”生成的（显然更快）代码：
mov edx, ecx
mov eax, 2155905153 ;
imul rax, rdx ; Replacing the IDIV with IMUL and SHR
shr rax, 39 ;
mov edx, eax
shl edx, 8
sub eax, edx
add eax, ecx
mov eax, ecx
mov edx, 2155905153
imul rdx, rax
shr rdx, 39
add edx, ecx
movzx eax, dl ; Faster than an explicit AND mask?
在 Windows-10（64 位）平台（英特尔® 酷睿™ i7-8550U CPU）上测试此代码表明，它显着（但不是很大）优于问题中提出的其他算法。
† David Eisenstat 给出的答案解释了这种等价性如何/为什么有效。
这是我对最快答案如何运作的看法。我还不知道基恩是否可以改进或易于推广。
给定一个整数 x ? 0, 让 q = ?x/255? (在 C, q = x / 255;) 和 r = x ? 255 q (in C, r = x % 255;) 所以 q ? 0 和 0 ? r < 255 是整数，x = 255 q + r。
此方法计算 (x + ?x/255?) mod 2 8（在 C 中，(x + x / 255) & 0xff），它等于 (255 q + r + q) mod 2 8 = (2 8 q + r) mod 2 8 = r。
注意 x + ?x/255? = ?x + x/255? = ?(2 8 /255) x?，其中第一步是从 x 是一个整数。此方法使用乘数 (2 0 + 2 ?8 + 2 ?16 + 2 ?24 + 2 ?32 ) 而不是 2 8 /255，后者是无限级数 2 0 + 2 ?8 + 2 ?16的总和+ 2 ?24 + 2 ?32 + .... 由于近似值稍低，该方法必须检测残差 2 8 ? 1 = 255。
这种方法的直觉是计算 y = (2 8 /255) x mod 2 8，它等于 (2 8 /255) (255 q + r) mod 2 8 = (2 8 q + (2 8 /255) r) mod 2 8 = (2 8 /255) r，然后返回 y ? y/2 8，等于 r。
由于这些公式不使用 ?(2 8 /255) r? = r，对于两个保护位，基恩可以从 2 8切换到 2 10。理想情况下，这些将始终为零，但由于定点截断和 2 10 /255的近似值，它们不是。基恩加 2 以从截断切换到舍入，这也避免了沃伦的特殊情况。
这种方法使用乘数 2 2 (2 0 + 2 ?8 + 2 ?16 + 2 ?24 + 2 ?32 + 2 ?40 ) = 2 2 (2 0 + 2 ?16 + 2 ?32 ) (2 0 + 2 ?8 )。C 语句x = (((x >> 16) + x) >> 14) + (x << 2);计算 x? = ?2 2 (2 0 + 2 ?16 + 2 ?32 ) x? 模式 2 32。然后((x >> 8) + x) & 0x3ff是x？？= ?(2 0 + 2 ?8 ) x?? 模式 2 10 .
我现在没有时间正式进行错误分析。非正式地，第一次计算的错误区间宽度<1；第二个，宽度 < 2 + 2 ?8；第三个，width < ((2 ? 2 ?8 ) + 1)/2 2 < 1，它允许正确的舍入。
关于改进，近似值的 2 ?40项似乎没有必要 (?)，但我们不妨拥有它，除非我们可以删除 2 ?32项。降低 2 ?32会使近似质量超出规范。
猜猜您可能不是在寻找需要快速 64 位乘法的解决方案，而是为了记录：
return (x * 0x101010101010102ULL) >> 56;
这种方法（自上次编辑后略有改进）融合了沃伦和基恩。在我的笔记本电脑上，它比 Keane 快，但不如 64 位乘法和移位快。它避免了乘法，但受益于单个旋转指令。与原始版本不同，它在 RISC-V 上可能没问题。
像沃伦一样，这种方法近似于 ?(256/255) x mod 256? 在 8.24 定点。Mod 256，每个字节 b 贡献一个项 (256/255) b，大约是 b.bbb 基数 256。此方法的原始版本只是对所有四个字节旋转求和。（我稍后会讲到修订版。）这个总和总是低估实际价值，但最后低估了不到 4 个单位。通过在截断之前添加 4/2 ?24，我们保证了正确的答案，就像在基恩中一样。
修订版通过放宽近似质量来节省工作。我们写 (256/255) x = (257/256) (65536/65535) x，在 16.16 定点计算 (65536/65535) x（即，将 x 添加到其 16 位旋转），然后乘以 257 /256 和 mod 由 256 变成 8.24 定点。第一个乘法在 16.16 的最后一位的误差小于 2 个单位，第二个是精确的 (!)。总和低于 (2/2 16 ) (257/256)，因此常数项 514/2 24足以修复截断。如果不同的立即数操作数更有效，也可以使用更大的值。
uint32_t mod255(uint32_t x) {
x += (x << 16) | (x >> 16);
return ((x << 8) + x + 514) >> 24;
}