有效的模255计算

Question

我试图找到最有效的方法来计算32位无符号整数的模255。我的主要关注点是找到一种在x86和ARM平台上运行良好的算法，并着眼于除此之外的适用性。首先，我试图避免内存操作（这可能很昂贵），因此我在避免表的同时寻找位细微的方法。我还试图避免潜在的昂贵操作，例如分支和乘法，并尽量减少使用的操作和寄存器的数量。

下面的ISO-C99代码捕获了我迄今为止尝试的八种变体。它包括一个用于穷举测试的框架。我附加了一些粗略的执行时间测量，它似乎工作得很好，足以获得第一个性能印象。在我尝试的几个平台上（都是快速整数乘法），变体WARREN\u MUL\u Shur\u 2、WARREN\u MUL\u Shur\u 1和DIGIT\u SUM\u CARRY\u OUT\u 1似乎是最有效的。我的实验表明，我在编译器资源管理器上尝试的x86、ARM、PowerPC和MIPS编译器都很好地利用了平台特定的功能，例如三输入LEA、字节扩展指令、乘法累加和指令预测。

变体NAIVE_USING_DIV使用整数除法，用除数反乘，后跟减法。这是基线情况。现代编译器知道如何有效地实现无符号整数除法255（通过乘法），并将在适当的情况下使用离散替换反乘。要计算模base-1，可以对base数字求和，然后折叠结果。例如3334 mod 9：sum3 3 3 4 = 13，折叠1 3=4。如果折叠后的结果是base-1，我们需要生成0。DIGIT_SUM_THEN_FOLD使用此方法。

A、 Cockburn，“使用8/16位算术有效实现OSI传输协议校验和算法”，ACM SIGCOM计算机通信评论，第17卷，第3期，7月/8月。1987年，第13-20页

显示了在校验和计算模255的上下文中有效添加数字模base-1的不同方式。计算数字的字节总和，并在每次加法后，也添加加法中的任何执行。因此，这将是一个ADD a， b，ADC a，0序列。使用base 256数字写出加法链，很明显，计算基本上是与0x0101 ... 0101相乘。结果将位于最高有效数字位置，只是需要单独捕获该位置的加法执行。此方法仅在base数字包含2^k位时有效。这里我们有k=3。我尝试了三种不同的方法将base-1的结果重新映射为0，从而产生变体DIGIT_SUM_CARRY_OUT_1、DIGIT_SUM_CARRY_OUT_2、DIGIT_SUM_CARRY_OUT_3。

Joe Keane在1995/07/09的新闻组comp.lang.c中展示了一种有趣的高效计算模63的方法。虽然线程参与者Peter L. Montgomery证明了算法是正确的，但不幸的是Keane先生没有回应解释其推导的请求。该算法也在H. Warren的Hacker's Delight第2版中复制。我能够以纯粹的机械方式将其扩展到模127和模255。这是（适当命名的）KEANE_MAGIC变体。更新：由于我最初发布了这个问题，我已经发现Keane的方法基本上是以下内容的巧妙定点实现：返回（uint32_t）（fmod（x*256.0 / 255.0 0.5, 256.0) * (255.0 / 256.0))；。这使它成为下一个变体的近亲。

Henry S. Warren，Hacker's Delight第2版，第272页显示了一种“乘移右”算法，大概是由作者自己设计的，它基于n mod 2^k-1=楼层（2^k/2^k-1*n）mod 2^k。固定点计算用于与因子2^k/2^k-1相乘。我构造了两个不同的变体，它们如何处理base-1到0的初步结果的映射。它们是变体WARREN_MUL_SHR_1和WARREN_MUL_SHR_2。

有没有比我迄今为止确定的三个顶尖竞争者更高效的模255计算算法，特别是对于整数乘法较慢的平台？在这种情况下，对基恩的无乘法算法进行有效修改，对四个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(__html" target="_blank">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;
}

Answer 1

我想你可能不是在寻找需要快速64位乘法的解决方案，而是为了记录在案：

return (x * 0x101010101010102ULL) >> 56;

Answer 2

这是我对最快答案的感觉。我还不知道基恩是否可以改进或容易推广。

给定整数x≥ 0，设q=⌊x/255⌋ （在C中，q=x/255；）和r=x− 255 q（在C中，r=x%255；），因此q≥ 0和0≤ r

此方法计算（x x/255）mod 2⁸（在C中，（x x/255）

请注意，x x/255=x x/255=（28/255）x，其中第一步从x是整数开始。此方法使用乘数（202−82−162−242−32）而不是28/255，这是无限级数202−82−162−242−32......由于近似值略低，因此此方法必须检测残差28−1=255。

该方法的直观性是计算y=（28/255）x mod 28，它等于（28/255）（255 q r）mod 28=（28=（28/255）mod 28=（28/255）r，并返回y− y/28，等于r。

由于这些公式没有使用（28/255）r=r这一事实，因此Keane可以为两个保护位从28切换到210。理想情况下，这些始终为零，但由于定点截断和210/255的近似值，它们不是。Keane添加2以从截断切换到舍入，这也避免了Warren中的特殊情况。

这种方法有点使用乘数22（202−82−162−242−322−40）=22（202−162−32）（202−8）。C语句x=（（（x

我现在没有时间正式进行错误分析。非正式地说，第一次计算的误差区间具有宽度

关于改进，2−40近似项似乎没有必要（？），但是我们也可以拥有它，除非我们可以去掉2−32术语。下降2−32将近似质量推到了规格之外。

Answer 3

对于任意无符号整数x和n，计算模表达式x%n涉及（至少在概念上）三个操作：除法、乘法和减法：

quotient = x / n;
product = quotient * n;
modulus = x - product;

然而，当n是2的幂（n=2^{p）时，可以更快速地确定模，只需屏蔽除较低p位之外的所有位。}

在大多数CPU上，加法、减法和位屏蔽是非常“便宜”（快速）的操作，乘法更“昂贵”，除法非常昂贵——但请注意，大多数优化编译器会将编译时常数除法转换为乘法（通过不同的常数）和位移位（见下文）。

因此，如果我们可以将模255转换为模256，而不需要太多开销，我们可能会加快这个过程。我们可以注意到，x%n相当于（x x/n）%（n 1）†。因此，我们的概念操作现在是：除法、加法和掩蔽。

在屏蔽较低8位的特定情况下，基于x86/x64的CPU（和其他？）可能会执行进一步的优化，因为它们可以访问（大多数）寄存器的8位版本。

下面是clang cl编译器为原始模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

下面是使用上述“技巧”生成的代码（明显更快）：

unsigned Trick255(unsigned x)
{
    return (x + x / 255) & 0xFF;
}

    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位）平台（Intel®Core）上测试此代码™ i7-8550U CPU）表明，它显著（但不是很大）优于问题中提出的其他算法。

David Eisenstat给出的答案解释了这种等价是如何/为什么有效的。