8

我有兴趣获得欧几里得除法的余数，即对于一对整数 (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 且 < MAX_INT 的任意值

编辑

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

  • 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] 是 x 1等...

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

4

12 回答 12

7

编辑：没有乘法或分支woot。

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
于 2009-07-16T04:36:10.263 回答
5

我认为 280Z28 和 Christopher 比我更好地覆盖了组装高尔夫，这涉及随机访问。

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

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

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

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

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

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

于 2009-07-16T09:26:09.257 回答
4

没有分支，但有点摆弄：

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
于 2009-07-16T10:57:18.823 回答
3

整数乘法比除法快得多。对于已知 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;
于 2009-07-16T06:13:27.383 回答
1

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

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

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

于 2009-07-16T05:32:28.183 回答
1 
int euc(int i, int n)
{
    return (i % n) + (((i % n) < 0) * n);
}
于 2009-07-16T04:55:56.290 回答
1

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

我的想法是 ((i%n)+n)%n （无分支）或 (i+(n<<16))%n （对于大 n 或极负 i 显然失败）会更快，但他们都花了同样的时间。

于 2009-07-16T06:27:53.343 回答
1

我真的很喜欢这个表达：

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，这可能很昂贵），并且可以完全内联，避免函数调用的开销。

你怎么看？

于 2009-07-16T06:38:39.393 回答
1

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

于 2009-07-16T06:40:56.813 回答
0

这是克里斯托弗的版本，如果右移不是算术，则可以回退到杰森的版本。

#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.

于 2009-07-16T12:18:44.667 回答
0

你在给 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
}
于 2009-07-16T11:09:16.380 回答
0

如果您可以保证数组的维度始终是 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);
    }
}
于 2009-07-16T11:42:57.797 回答