【发布时间】：2010-11-04 23:43:54
【问题描述】：
地板除法是当结果总是向下（朝向-∞）而不是朝向0时：
是否可以在 C/C++ 中有效地实现整数除法或欧式整数除法？
（显而易见的解决方案是检查被除数的符号）
【问题讨论】：
标签： c++ c math division integer-division
【发布时间】：2010-11-04 23:43:54
【问题描述】：
我编写了一个测试程序来对这里提出的想法进行基准测试：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#define N 10000000
#define M 100
int dividends[N], divisors[N], results[N];
__forceinline int floordiv_signcheck(int a, int b)
{
return (a<0 ? a-(b-1) : a) / b;
}
__forceinline int floordiv_signcheck2(int a, int b)
{
return (a - (a<0 ? b-1 : 0)) / b;
}
__forceinline int floordiv_signmultiply(int a, int b)
{
return (a + (a>>(sizeof(a)*8-1))*(b-1)) / b;
}
__forceinline int floordiv_floatingpoint(int a, int b)
{
// I imagine that the call to floor can be replaced to a cast
// if you can get FPU rounding control to work (I couldn't).
return floor((double)a / b);
}
void main()
{
for (int i=0; i<N; i++)
{
dividends[i] = rand();
do
divisors[i] = rand();
while (divisors[i]==0);
}
LARGE_INTEGER t0, t1;
QueryPerformanceCounter(&t0);
for (int j=0; j<M; j++)
for (int i=0; i<N; i++)
results[i] = floordiv_signcheck(dividends[i], divisors[i]);
QueryPerformanceCounter(&t1);
printf("signcheck : %9llu\n", t1.QuadPart-t0.QuadPart);
QueryPerformanceCounter(&t0);
for (int j=0; j<M; j++)
for (int i=0; i<N; i++)
results[i] = floordiv_signcheck2(dividends[i], divisors[i]);
QueryPerformanceCounter(&t1);
printf("signcheck2 : %9llu\n", t1.QuadPart-t0.QuadPart);
QueryPerformanceCounter(&t0);
for (int j=0; j<M; j++)
for (int i=0; i<N; i++)
results[i] = floordiv_signmultiply(dividends[i], divisors[i]);
QueryPerformanceCounter(&t1);
printf("signmultiply : %9llu\n", t1.QuadPart-t0.QuadPart);
QueryPerformanceCounter(&t0);
for (int j=0; j<M; j++)
for (int i=0; i<N; i++)
results[i] = floordiv_floatingpoint(dividends[i], divisors[i]);
QueryPerformanceCounter(&t1);
printf("floatingpoint: %9llu\n", t1.QuadPart-t0.QuadPart);
}
结果：
signcheck : 61458768
signcheck2 : 61284370
signmultiply : 61625076
floatingpoint: 287315364
所以，根据我的结果，检查标志是最快的：
(a - (a<0 ? b-1 : 0)) / b
【讨论】：
-
除第一个以外的所有时间都包括
printf用于上一个答案。printf有点慢，不知道是不是慢到影响你的结果。 -
您可能还会受到编译器做出的内联决策的影响，值得查看生成的程序集。尝试
double而不是float，因为那里可能也会发生转化。 -
感谢您的建议，我已经更新了程序。使用
double会使浮点运行得更快一些，但不会快很多。否则，结果并没有太大变化。 -
这个程序不是只测试正股利吗？
rand()被指定为只返回正整数。 -
嘿，那是真的。也许这需要重新审视。
五年后我将重新审视这个问题，因为这对我也很重要。我对 x86-64 的两个纯 C 版本和两个内联汇编版本进行了一些性能测量，结果可能很有趣。
已测试的地板除法变体是：
- 我已经使用了一段时间的实现；
- 上面介绍的轻微变体，仅使用一种划分；
- 前一个，但在 inline-assembly 中手动实现；和
- 在程序集中实现的
CMOV版本。
以下是我的基准程序：
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#ifndef VARIANT
#define VARIANT 3
#endif
#if VARIANT == 0
#define floordiv(a, b) (((a) < 0)?((((a) + 1) / (b)) - 1):((a) / (b)))
#elif VARIANT == 1
#define floordiv(a, b) ((((a) < 0)?((a) - ((b) - 1)):(a)) / (b))
#elif VARIANT == 2
#define floordiv(a, b) ({ \
int result; \
asm("test %%eax, %%eax; jns 1f; sub %1, %%eax;" \
"add $1, %%eax; 1: cltd; idivl %1;" \
: "=a" (result) \
: "r" (b), \
"0" (a) \
: "rdx"); \
result;})
#elif VARIANT == 3
#define floordiv(a, b) ({ \
int result; \
asm("mov %%eax, %%edx; sub %1, %%edx; add $1, %%edx;" \
"test %%eax, %%eax; cmovs %%edx, %%eax; cltd;" \
"idivl %1;" \
: "=a" (result) \
: "r" (b), \
"0" (a) \
: "rdx"); \
result;})
#endif
double ntime(void)
{
struct timeval tv;
gettimeofday(&tv, NULL);
return(tv.tv_sec + (((double)tv.tv_usec) / 1000000.0));
}
void timediv(int n, int *p, int *q, int *r)
{
int i;
for(i = 0; i < n; i++)
r[i] = floordiv(p[i], q[i]);
}
int main(int argc, char **argv)
{
int n, i, *q, *p, *r;
double st;
n = 10000000;
p = malloc(sizeof(*p) * n);
q = malloc(sizeof(*q) * n);
r = malloc(sizeof(*r) * n);
for(i = 0; i < n; i++) {
p[i] = (rand() % 1000000) - 500000;
q[i] = (rand() % 1000000) + 1;
}
st = ntime();
for(i = 0; i < 100; i++)
timediv(n, p, q, r);
printf("%g\n", ntime() - st);
return(0);
}
我使用 GCC 4.9.2 使用 gcc -march=native -Ofast 编译了这个，结果在我的 Core i5-2400 上如下。每次运行的结果都具有相当的可重复性——至少它们总是以相同的顺序降落。
- 变体 0：7.21 秒
- 变体 1：7.26 秒
- 变体 2：6.73 秒
- 变体 3：4.32 秒
所以CMOV 的实现至少将其他人从水中吹了出来。令我惊讶的是，变体 2 以相当大的优势超越了它的纯 C 版本（变体 1）。我原以为编译器应该能够发出至少和我一样高效的代码。
这里有一些其他平台，供比较：
AMD 速龙 64 X2 4200+，GCC 4.7.2：
- 变体 0：26.33 秒
- 变体 1：25.38 秒
- 变体 2：25.19 秒
- 变体 3：22.39 秒
至强 E3-1271 v3，GCC 4.9.2：
- 变体 0：5.95 秒
- 变体 1：5.62 秒
- 变体 2：5.40 秒
- 变体 3：3.44 秒
作为最后一点，我或许应该警告不要过于认真地对待CMOV 版本的明显性能优势，因为在现实世界中，其他版本中的分支可能不会像在这个基准测试中那样完全随机，并且如果分支预测器可以做一个合理的工作，分支版本可能会变得更好。但是，实际情况在很大程度上取决于实际使用的数据，因此尝试进行任何通用基准测试可能毫无意义。
【讨论】：
-
请注意，当除数为正时，变体 0 和 1（未检查其他变体）仅给出正确的地板除法结果。 （在这种情况下，结果也与欧几里得除法一致。）例如，如果 a,b = 5,-3，这些公式表示商为 -1，但地板除法的结果应为 -2。 （5,-3 的欧几里得除法是 -1，但是当两个操作数都是负数时，这些公式给出的答案与下限和欧几里得除法不同。）
-
您是否对除数为正常数的版本进行了基准测试（但除数是在运行时计算的）？对于二的幂的红利，右移运算符几乎肯定会和任何替代方法一样快，并且可能比任何有效的方法都更具可读性。对于其他红利，用于地板/欧几里得除法的最佳代码应该比用于截断的代码更快，但我不知道有什么可靠的方法可以说服编译器生成好的代码。
-
@supercat：对于编译时常量除数，编译器将能够使用this trick 生成更好的代码，无论它是否是二次幂。所以是的，具有编译时常数除数的除法肯定比任何处理动态除数的算法都要快。
-
@Dolda2000：我的问题是，在正常数除数的情况下，必须做些什么才能为地板/欧几里得除法生成有效代码？使用带符号的
/运算符通常会导致编译器生成额外的代码来创建接近零的不均匀性，这将使程序员有必要编写更多的额外代码。我认为最好的方法是转换为无符号，加一个偏移量，做一个除法，然后减去一个不同的偏移量，但这似乎有点恶心。 -
@supercat：如果你想研究这项技术，可以在 publicly available 上找到一篇关于该主题的广泛论文。它包括截断除法和取整除法的情况。
提出一些免费的分支来根据符号更正结果可能会更有效，因为分支很昂贵。
请参阅Hacker's Delight 中Chapter 2 的第 20 页，了解如何访问该标志。
【讨论】：
-
我唯一能想到的就是符号位与除数的乘法。
-
您可以用 0 或 1 表示某些条件并将它们添加到结果中。
请注意：x86 sar 指令在涉及 2 的幂时执行地板除法。
【讨论】：
由于 IEEE-754 将向 -inf 的舍入指定为所需的舍入模式之一，我想您的问题的答案是肯定的。但是也许你可以解释一下，如果你想知道如果一个人正在编写编译器，你将如何实现这个过程，或者想知道如何使用一个特定的编译器来执行这个操作？
【讨论】：
-
哎呀，我的意思是整数除法！我会修改问题。我不确定编写编译器与这个问题有什么关系，因为我已经检查过并且 C/C++ 会截断除法。
-
@CyberShadow：所以将除数和除数转换为某种浮点数，使用 floor 并将结果转换回整数。还是您在寻找其他类型的答案？
-
我不认为使用浮点数会比检查符号更快...
-
hardwaresecrets.com/fullimage.php?image=6761 似乎表明双精度除法可以与整数除法竞争（仍然存在转换成本，并且可能每次都切换舍入模式）。但是，这种方法仍然不适用于全部 64 位整数，并且如果目标平台没有（完全）IEEE-754 浮点计算（包括更改舍入的函数），它（完全）不起作用模式）。
-
我测试了它，性能很糟糕，但我可能做错了。看看我在自己的答案中发布的程序 - 弄乱舍入模式会使事情变得更快吗？
是否可以在 C/C++ 中有效地实现下限或欧几里得整数除法？
是的。
（显而易见的解决方案是检查被除数的符号）
我完全同意，并且很难相信存在明显更快的替代方案。
【讨论】：