【问题标题】：Fast Euclidean division in CC中的快速欧几里得除法
【发布时间】：2010-11-11 05:54:47
【问题描述】：

我有兴趣得到Euclidean除法的余数，即对于一对整数（i，n），求r如：

i = k * n + r, 0 <= r < |k|

简单的解决方案是：

int euc(int i, int n)
{
    int r;

    r = i % n;
    if ( r < 0) {
        r += n;
    }
    return r;
}

但是由于我需要执行数千万次（它在多维数组的迭代器中使用），所以我想尽可能避免分支。要求：

分支但更快也是可取的。
只适用于正 n 的解决方案是可以接受的（但它必须适用于负 i）。
n 是事先不知道的，可以是任何值 > 0 和

编辑

实际上很容易得到错误的结果，所以这里是一个预期结果的例子：

euc(0, 3) = 0
euc(1, 3) = 1
euc(2, 3) = 2
euc(3, 3) = 0
euc(-1, 3) = 2
euc(-2, 3) = 1
euc(-3, 3) = 0

有些人还担心优化这个没有意义。我需要一个多维迭代器，其中超出范围的项目被重复原始数组的“虚拟数组”中的项目替换。所以如果我的数组 x 是 [1, 2, 3, 4]，虚拟数组是 [...., 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4]，比如x[-2]就是x1等...

对于维度为 d 的 nd 数组，我需要对每个点进行 d 欧几里得除法。如果我需要在 n^d 数组与 m^d 内核之间进行关联，我需要 n^d * m^d * d 欧几里得除法。对于 100x100x100 点的 3d 图像和 5*5*5 点的内核，这已经是约 4 亿欧几里得分割。

【问题讨论】：

N 是否总是正数 (>0)，如您的示例所示？或者我们可能会遇到负 N 值？
对我来说，n 总是正数，是的。
带环绕的多维数组？嗯，甜甜圈:-)
所有的部门都使用相同的n吗？
标准 % 运算符有什么问题？真的太慢了吗？

标签： c bit-manipulation micro-optimization

【解决方案1】：

编辑：没有乘法或分支。

int euc(int i, int n)
{
    int r;

    r = i % n;
    r += n & (-(r < 0));

    return r;
}

这是生成的代码。根据 MSVC++ 检测分析器（我的测试）和 OP 的测试，它们的性能几乎相同。

; Original post
00401000  cdq              
00401001  idiv        eax,ecx 
00401003  mov         eax,edx 
00401005  test        eax,eax 
00401007  jge         euc+0Bh (40100Bh) 
00401009  add         eax,ecx 
0040100B  ret              

; Mine
00401020  cdq              
00401021  idiv        eax,ecx 
00401023  xor         eax,eax 
00401025  test        edx,edx 
00401027  setl        al   
0040102A  neg         eax  
0040102C  and         eax,ecx 
0040102E  add         eax,edx 
00401030  ret

【讨论】：

很有可能编译器实际上会使用分支来计算 (r
不幸的是，这种方法也不起作用。例如，euc(-3, 3) 返回 3（如果 n == 3，则返回值应在 [0, 2] 中，且 -3 = -1 * 3 + 0）。
糟糕，这是因为我将 r 更改为 i。固定。
在 core 2 duo 上，两者基本相同，使用带有 rdtsc 的简单循环进行测试。看起来要显着加快速度是相当困难的。
我在 MSVC++ 检测分析器中得到了类似的结果。

【解决方案2】：

int euc(int i, int n)
{
    return (i % n) + (((i % n) < 0) * n);
}

【讨论】：

【解决方案3】：

如果您还可以保证 i 永远不会小于 -n，您可以简单地将可选加法放在模数之前。这样一来，您就不需要分支了，如果您不需要，模会删除您添加的内容。

int euc(int i, int n)
{
    return (i + n) % n;
}

如果i小于-n，你仍然可以使用这个方法。在这种情况下，您可能确切地知道您的值将在什么范围内。因此，您可以将 x*n 添加到 i，而不是添加 n 到 i，其中 x 是给您足够范围的任何整数。为了提高速度（在没有单周期乘法的处理器上），您可以左移而不是乘法。

【讨论】：

这是一个简单而好的解决方案。在 core duo 上它也不是真的更快，但它比 pentium 4 上的其他方法快得多（快 15 %）。
@arke：我相信我的代码对于任何 i 都是正确的：您的代码为 f(-5, 3) 返回 -2，但我的代码返回 1，正如预期的那样。
哎呀，你是对的。不过，您始终可以在添加 n 之前对其进行缩放 :)。现在我考虑一下，我会左移而不是乘法。在 P4 上应该仍然更快。

【解决方案4】：

整数乘法比除法快得多。对于已知 N 的大量调用，您可以通过乘以 N 的伪逆来代替除以 N。

我将举例说明这一点。取 N=29。然后计算一次伪逆 2^16/N：K=2259（从 2259.86 截断...）。我假设 I 是正数并且 I*K 适合 32 位。

Quo = (I*K)>>16;   // replaces the division, Quo <= I/N
Mod = I - Quo*N;   // Mod >= I%N
while (Mod >= N) Mod -= N;  // compensate for the approximation

在我的示例中，假设 I=753，我们得到 Quo=25 和 Mod=28。（无需补偿）

编辑。

在您的 3D 卷积示例中，对 i%n 的大多数调用将在 0..n-1 中使用 i，因此在大多数情况下，第一行就像

if (i>=0 && i<n) return i;

将绕过昂贵且在这里无用的 idiv。

此外，如果您有足够的 RAM，只需将所有维度对齐到 2 的幂并使用位操作（移位和）而不是除法。

编辑 2.

我实际上在 10^9 次通话中尝试过。我％n：2.93s，我的代码：1.38s。请记住，这意味着对 I 的限制（I*K 必须适合 32 位）。

另一个想法：如果你的值是 x+dx，x 在 0..n-1 和 dx 小，那么以下将涵盖所有情况：

if (i<0) return i+n; else if (i>=n) return i-n;
return i;

【讨论】：

你说得对，我的分析非常粗略，事实上，在迭代器中，我在处理越界情况之前就做了这个测试。我不能假设任何关于内存的事情（如果内存不是问题，那么在大多数情况下，基于 FFT 的卷积/相关无论如何都会快得多）。

【解决方案5】：

我在 gcc -O3 中使用 TSC 对每个人的提案进行计时（除了常数 N 的提案），它们都花费了相同的时间（在 1% 以内）。

我的想法是 ((i%n)+n)%n（无分支）或 (i+(n

【讨论】：

我要补充一点，到目前为止，这仍然是正确的——没有人有更快的实现（如果你在没有优化的情况下编译会有很大的不同；没有）。

【解决方案6】：

我真的很喜欢这个表达：

r = ((i%n)+n)%n;

反汇编很短：

r = ((i%n)+n)%n;

004135AC  mov         eax,dword ptr [i] 
004135AF  cdq              
004135B0  idiv        eax,dword ptr [n] 
004135B3  add         edx,dword ptr [n] 
004135B6  mov         eax,edx 
004135B8  cdq              
004135B9  idiv        eax,dword ptr [n] 
004135BC  mov         dword ptr [r],edx

它没有跳转（2 个 idiv，这可能很昂贵），并且可以完全内联，避免函数调用的开销。

你怎么看？

【讨论】：

比我的还长，也没有跳转，可以内联？这是另一个聪明的。 :)
顺便说一句：那是未优化的程序集。通过优化，编译器应该避免 [n] 的冗余负载，而是使用寄存器。
这使用了两个idiv 指令，因此比所有其他建议的成本都要高！

【解决方案7】：

如果您有足够低的范围，请创建一个查找表 - 两个暗淡的数组。您也可以将函数设为内联，并通过查看生成的代码来确保它是内联的。

【讨论】：

【解决方案8】：

我认为 280Z28 和 Christopher 比我更好地涵盖了组装高尔夫，并且处理随机访问。

不过，您实际上在做的似乎是处理整个数组。显然，出于内存缓存的原因，您已经希望尽可能按顺序执行此操作，因为避免缓存未命中比避免小分支要好很多很多倍。

在这种情况下，首先通过适当的边界检查，您可以执行内部循环，我将称之为“破折号”。检查下一个 k 增量不会导致任一数组的最小维度溢出，然后使用新的更内部循环“破折号”k 步，该循环每次只将“物理”索引增加 1做另一个idiv。您或编译器可以展开此循环，使用 Duff 的设备等。

如果内核很小，特别是如果它是固定大小的，那么那个（或它的倍数，适当展开以偶尔减去而不是加法），可能是用于“破折号”长度的值.编译时常量破折号长度可能是最好的，因为那时您（或编译器）可以完全展开破折号循环并省略延续条件。只要这不会使代码变得太大而无法快速运行，它实质上就是将整个正模运算替换为整数增量。

如果内核不是固定大小，但其最后一维通常非常小，请考虑为最常见的大小设置不同版本的比较函数，并在每个版本中完全展开破折号循环。

另一种可能性是计算将发生溢出的下一个点（在任一数组中），然后冲到该值。您在破折号循环中仍然有一个延续条件，但它只使用增量来尽可能长。

或者，如果您正在执行的操作是数字相等或其他一些简单的操作（我不知道“相关性”是什么），您可以查看 SIMD 指令或其他什么，在这种情况下，破折号长度应该是架构上最广泛的单指令比较（或适当的 SIMD 操作）的倍数。不过，这不是我有过的经验。

【讨论】：

我认为问题的原始版本没有提到应用程序，但这是迄今为止唯一有用的答案。没有人有更快的 euc 实现，但它绝对可以加快更高级别的计算。
原来的问题只是暗示，说它被使用了数百万次。我想我是在编辑说我们有（至少有时）顺序访问之后才发布的。公平地说，有用！= 回答问题：OP 要求无分支实现，而不是快速实现，并说“分支但更快也是可取的”;-) 即使使用 -O3，OP 的代码在 GCC 3 上也有一个分支（我现在这台机器上的所有东西）。

【解决方案9】：

没有分支，但有点摆弄：

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (((unsigned int)r) >> 31) * n;
    return r;
}

没有乘法：

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (r >> 31) & n;
    return r;
}

这给出了：

; _i$ = eax
; _n$ = ecx

cdq
idiv   ecx
mov eax, edx
sar eax, 31
and eax, ecx
add eax, edx

【讨论】：

第二个版本只有在有符号右移是算术的情况下才有效
它还假定 32 位整数，因此源必须附带移植指南（或使用 stdint.h）。您可以通过强制转换为无符号来强制进行逻辑转换，然后操作：~((((unsigned)r)>>31)-1)。什么的。
您可以轻松替换 31。例如(sizeof( int ) * 8 - 1) 当不存在算术移位时，你必须放弃它，并使用更合适的东西。 ~((((unsigned)r)>>31)-1) 乍一看还不错。
其实，如果你想真正便携，并不能保证一个char是8位的。 CHAR_BIT 解决了这个问题。但是，也不能保证 int 存储表示的所有位都参与该值。 IIRC 只有 char 有这个保证。一个符合要求的实现可以有一个 9 位字节和 32 位整数（在这种情况下 sizeof(int)*8 有效，但 sizeof(int)*CHAR_BIT 无效），或一个 9 位字节和 36 位整数（在这种情况下，sizeof(int)*CHAR_BIT 有效，但 sizeof(int)*8 无效't. 我认为这就是 C++ 有 numeric_limits::digits 的原因。
... 总而言之，我认为您只需要搬运工做一些特定于平台的事情，或者如果他们不想或无法击败它，就选择一些不太好的便携式版本.在某些平台上，搬运工无论如何都希望在 asm 中执行此操作。毕竟，这几乎就是静态内联函数的用途:-)

【解决方案10】：

您在回答 Eric Bainville 时说，大多数时候0 <= i < n 并且您有

if (i>=0 && i<n) return i;

无论如何作为euc() 的第一行。

无论如何，当您进行比较时，您不妨使用它们：

int euc(int i, int n)
{
    if (n <= i)            return i % n;
    else if (i < 0)        return ((i + 1) % n) + n - 1;
    else /* 0 <= i < n */  return i;  // fastest possible response for common case
}

【讨论】：

不适用于i = -n，因为它将返回n 而不是0

【解决方案11】：

如果您可以保证数组的维度始终是 2 的幂，那么您可以这样做：

r = (i & (n - 1));

如果您可以进一步保证您的维度来自给定的子集，您可以这样做：

template<int n>
int euc(int i) {
    return (i & (n - 1));
}

int euc(int i, int n) {
    switch (n) {
        case 2: return euc<2>(i);
        case 4: return euc<4>(i);
    }
}

【讨论】：

【解决方案12】：

这里是Christopher's version，如果右移不是算术，则回退到Jason's。

#include <limits.h>
static inline int euc(int i, int n)
{
    // check for arithmetic shift
    #if (-1 >> 1) == -1
        #define OFFSET ((i % n >> (sizeof(int) * CHAR_BIT - 1)) & n)
    #else
        #define OFFSET ((i % n < 0) * n)
    #endif

    return i % n + OFFSET;
}

回退版本应该更慢，因为它使用imul 而不是and。

【讨论】：

正如我前面提到的，这些都是与 g++ -O3 相同的速度（甚至是原始代码）。