快速整除测试（2、3、4、5、..、16）？

Question

最快的整除测试是什么？比如说，给定一个 little-endian 架构和一个 32 位有符号整数：如何快速计算出一个数字可以被 2、3、4、5、... 最多 16 整除？

警告：给定的代码仅为示例。每条线都是独立的！在许多没有 DIV 硬件（如许多 ARM）的处理器上，使用模运算的明显解决方案很慢。一些编译器也无法进行此类优化（例如，如果除数是函数的参数或依赖于某些东西）。

Divisible_by_1 = do();
Divisible_by_2 = if (!(number & 1)) do();
Divisible_by_3 = ?
Divisible_by_4 = ?
Divisible_by_5 = ?
Divisible_by_6 = ?
Divisible_by_7 = ?
Divisible_by_8 = ?
Divisible_by_9 = ?
Divisible_by_10 = ?
Divisible_by_11 = ?
Divisible_by_12 = ?
Divisible_by_13 = ?
Divisible_by_14 = ?
Divisible_by_15 = ?
Divisible_by_16 = if(!number & 0x0000000F) do();

和特殊情况：

Divisible_by_2k = if(number & (tk-1)) do();  //tk=2**k=(2*2*2*...) k times

Answer 1

可分性的快速测试很大程度上取决于表示数字的基数。如果基数为 2，我认为您只能对可被 2 的幂整除进行“快速测试”。如果该数字的最后 n 个二进制数字为 0，则二进制数可被 2ⁿ 整除. 对于其他测试，我认为您通常找不到比% 更快的东西。

Answer 2

在所有情况下（包括可被 2 整除）：

if (number % n == 0) do();

使用低位掩码进行与运算只是混淆，使用现代编译器不会比以可读方式编写代码快。

如果您必须测试所有案例，您可以通过将一些案例放在if 中来提高性能：例如，如果被 2 整除失败，那么测试被 4 整除是没有意义的.

Answer 3

您应该只使用 (i % N) == 0 作为您的测试。

我的编译器（一个相当旧的 gcc 版本）为我尝试的所有情况生成了很好的代码。 在适合位测试的地方，它做到了。当 N 是一个常数时，它在任何情况下都不会产生明显的“除法”，它总是使用一些“技巧”。

只要让编译器为你生成代码，它几乎肯定会比你更了解机器的架构 :) 这些都是简单的优化，你不可能想出比编译器更好的东西。

不过，这是一个有趣的问题。我无法列出编译器对每个常量使用的技巧，因为我必须在不同的计算机上编译.. 但是如果没有人能打败我，我稍后会更新这个回复:)

Answer 4

有点开玩笑，但假设你得到了其余的答案：

Divisible_by_6  = Divisible_by_3 && Divisible_by_2;
Divisible_by_10 = Divisible_by_5 && Divisible_by_2;
Divisible_by_12 = Divisible_by_4 && Divisible_by_3;
Divisible_by_14 = Divisible_by_7 && Divisible_by_2;
Divisible_by_15 = Divisible_by_5 && Divisible_by_3;

Answer 5

正如@James 提到的，让编译器为您简化它。如果n 是一个常数，那么任何下降编译器都能够识别该模式并将其更改为更有效的等价物。

例如代码

#include <stdio.h>

int main() {
    size_t x;
    scanf("%u\n", &x);
    __asm__ volatile ("nop;nop;nop;nop;nop;");
    const char* volatile foo = (x%3 == 0) ? "yes" : "no";
    __asm__ volatile ("nop;nop;nop;nop;nop;");
    printf("%s\n", foo);
    return 0;
}

用g++-4.5 -O3编译，x%3 == 0的相关部分会变成

mov    rcx,QWORD PTR [rbp-0x8]   # rbp-0x8 = &x
mov    rdx,0xaaaaaaaaaaaaaaab
mov    rax,rcx
mul    rdx
lea    rax,"yes"
shr    rdx,1
lea    rdx,[rdx+rdx*2]
cmp    rcx,rdx
lea    rdx,"no"
cmovne rax,rdx
mov    QWORD PTR [rbp-0x10],rax

翻译回C代码的意思是

(hi64bit(x * 0xaaaaaaaaaaaaaaab) / 2) * 3 == x ? "yes" : "no"
// equivalatent to:                 x % 3 == 0 ? "yes" : "no"

这里不涉及分工。 （注意0xaaaaaaaaaaaaaaab == 0x20000000000000001L/3）

---

编辑：

• 魔术常数 0xaaaaaaaaaaaaaaab 可以在http://www.hackersdelight.org/magic.htm 中计算
• 对于 2ⁿ - 1 形式的除数，请检查 http://graphics.stanford.edu/~seander/bithacks.html#ModulusDivision

Answer 6

您可以用乘法代替非二次幂常数的除法，本质上是乘以除数的倒数。用这种方法得到准确结果的细节很复杂。

Hacker's Delight 在第 10 章中详细讨论了这一点（不幸的是，网上没有提供）。

从商可以通过另一个乘法和减法得到模数。

Answer 7

以下是一些我还没有看到其他人建议的提示：

一个想法是使用switch 语句，或预先计算一些数组。然后，任何体面的优化器都可以简单地直接索引每个案例。例如：

// tests for (2,3,4,5,6,7)
switch (n % 8)
{
case 0: break;
case 1: break;
case 2: do(2); break;
case 3: do(3); break;
case 4: do(2); do(4) break;
case 5: do(5); break;
case 6: do(2); do(3); do(4); break;
case 7: do(7); break;
} 

您的应用程序有点模棱两可，但您可能只需要检查小于 n=16 的素数。这是因为所有数字都是当前或先前素数的因数。因此，对于 n=16，您可能只需以某种方式检查 2, 3, 5, 7, 11, 13 就可以逃脱。只是一个想法。

Answer 8

这些数字的 LCM 似乎是 720720。它非常小，因此您可以执行单模运算并将余数用作预计算 LUT 中的索引。

Answer 9

在previous question 中，我展示了一种快速算法，用于在基数 N 中检查作为 N-1 因数的除数。 2 的不同幂之间的基本变换是微不足道的；这只是一些分组。

因此，在基数 4 中检查 3 很容易；在 base 16 中检查 5 很容易，在 base 64 中检查 7（和 9）很容易。

非素因数是微不足道的，所以只有 11 和 13 是困难的情况。对于 11，您可以使用基数 1024，但此时它对于小整数并不是很有效。

Answer 10

找出除法指令的替代方案（包括 x86/x64 上的模数）并不是一个坏主意，因为它们非常慢。比大多数人意识到的要慢（甚至慢得多）。那些建议 n 是变量的“% n”的人给出了愚蠢的建议，因为它总是会导致使用除法指令。另一方面，“% c”（其中 c 是一个常数）将允许编译器确定其曲目中可用的最佳算法。有时是除法指令，但很多时候不是。

在this documentTorbjörn Granlund 中显示，无符号 32 位 mults:divs 所需的时钟周期比在 Sandybridge 上为 4:26 (6.5x)，在 K10 上为 3:45 (15x)。对于 64 位，相应的比率为 4:92 (23x) 和 5:77 (14.4x)。

“L”列表示延迟。 “T”列表示吞吐量。这与处理器并行处理多条指令的能力有关。 Sandybridge 可以每隔一个周期发出一个 32 位乘法或每个周期发出一个 64 位乘法。对于 K10，相应的吞吐量是相反的。对于分区，K10 需要完成整个序列才能开始另一个序列。我怀疑 Sandybridge 也是如此。

以 K10 为例，这意味着在 32 位除法 (45) 所需的周期内，可以发出相同数量 (45) 的乘法，并且其中的倒数第二个和最后一个将完成除法完成后的一个和两个时钟周期。 45次乘法可以完成很多工作。

值得注意的是，随着从 K8-K9 到 K10 的演进，div 的效率变得越来越低：32 位和 64 位从 39 到 45 和 71 到 77 个时钟周期。

Granlund 在 gmplib.org 的 page 和斯德哥尔摩的 Royal Institute of Technology 包含更多好东西，其中一些已被纳入 gcc 编译器。

Answer 11

这可能对你的代码没有帮助，但在某些情况下，有一个巧妙的技巧可以帮助你做到这一点：

对于除以 3：对于一个十进制表示的数字，您可以将所有数字相加，并检查和是否可以被 3 整除。

示例：12345 => 1+2+3+4+5 = 15 => 1+5 = 6，可被 3 整除 (3 x 4115 = 12345)。

更有趣的是，相同的技术适用于 X-1 的所有因子，其中 X 是表示数字的基数。因此，对于十进制数，您可以检查除以 3 或 9。对于十六进制，您可以检查除以 3,5 或 15。对于八进制数，您可以检查除以 7。

Answer 12

一点邪恶的、令人困惑的位旋转可以让你除以 15。

对于 32 位无符号数：

def mod_15ish(unsigned int x) {
  // returns a number between 0 and 21 that is either x % 15
  // or 15 + (x % 15), and returns 0 only for x == 0
  x = (x & 0xF0F0F0F) + ((x >> 4) & 0xF0F0F0F);
  x = (x & 0xFF00FF) + ((x >> 8) & 0xFF00FF);  
  x = (x & 0xFFFF) + ((x >> 16) & 0xFFFF);
  // *1
  x = (x & 0xF) + ((x >> 4) & 0xF);
  return x;
}

def Divisible_by_15(unsigned int x) {
  return ((x == 0) || (mod_15ish(x) == 15));
}

您可以基于mod_15ish 为3 和5 构建类似的除法例程。

如果你有 64 位无符号整数要处理，以显而易见的方式扩展 *1 行上方的每个常量，并在 *1 行上方添加一行以使用掩码右移 32 位的0xFFFFFFFF。 （最后两行可以保持不变）mod_15ish 然后遵循相同的基本合同，但返回值现在在0 和31 之间。 （所以维护的是x % 15 == mod_15ish(x) % 15）

Answer 13

需要考虑的一件事：由于您只关心最多 16 的整除性，因此您实际上只需要检查最多 16 的素数的整除性。它们是 2、3、5、7、11 和 13。

将您的数字除以每个素数，并使用布尔值进行跟踪（例如 div2 = true）。数字二和三是特殊情况。如果 div3 为真，请尝试再次除以 3，设置 div9。二和它的力量非常简单（注意：'&'是处理器可以做的最快的事情之一）：

if n & 1 == 0:
    div2 = true
    if n & 3 == 0: 
        div4 = true
        if n & 7 == 0: 
            div8 = true
            if n & 15 == 0:
                div16 = true

您现在拥有布尔值 div2、div3、div4、div5、div7、div8、div9、div11、div13 和 div16。全部 其他数字是组合；例如 div6 与 (div2 && div3) 相同

因此，您只需要进行 5 或 6 次实际除法（仅当您的数字可以被 3 整除时才需要 6 次）。

就我自己而言，我可能会使用单个寄存器中的位作为布尔值；例如 bit_0 表示 div2。然后我可以使用面具：

if (flags & (div2+div3)) == (div2 + div3): do_6()

注意 div2+div3 可以是预先计算的常数。如果 div2 是 bit0，div3 是 bit1， 然后 div2+div3 == 3。这使得上面的 'if' 优化为：

if (flags & 3) == 3: do_6()

所以现在...没有分隔线的模组：

def mod(n,m):
    i = 0
        while m < n:
            m <<= 1
            i += 1
        while i > 0:
            m >>= 1
            if m <= n: n -= m
            i -= 1
     return n

div3 = mod(n,3) == 0
...

顺便说一句：上述代码的最坏情况是 32 位数字通过任一循环 31 次

仅供参考：刚刚看了上面的 Msalter 的帖子。对于某些素数，可以使用他的技术代替 mod(...)。

Answer 14

首先，我提醒您，二进制中 bn...b2b1b0 形式的数字具有值：

number = bn*2^n+...+b2*4+b1*2+b0

现在，当你说数字 %3 时，你有：

number%3 =3= bn*(2^n % 3)+...+b2*1+b1*2+b0

（我使用 =3= 表示模 3 同余）。还要注意b1*2 =3= -b1*1

现在，我将使用 + 和 - 以及可能的乘法来编写所有 16 个除法（请注意，乘法可以写为移动到不同位置的相同值的移位或总和。例如，5*x 表示您计算的 x+(x<<2) x 一次）

我们将号码称为n，假设Divisible_by_i 是一个布尔值。作为一个中间值，假设Congruence_by_i 是一个与n 模i 一致的值。

另外，假设 n0 表示 n 的第 0 位，n1 表示第 1 位等，即

ni = (n >> i) & 1;

Congruence_by_1 = 0
Congruence_by_2 = n&0x1
Congruence_by_3 = n0-n1+n2-n3+n4-n5+n6-n7+n8-n9+n10-n11+n12-n13+n14-n15+n16-n17+n18-n19+n20-n21+n22-n23+n24-n25+n26-n27+n28-n29+n30-n31
Congruence_by_4 = n&0x3
Congruence_by_5 = n0+2*n1-n2-2*n3+n4+2*n5-n6-2*n7+n8+2*n9-n10-2*n11+n12+2*n13-n14-2*n15+n16+2*n17-n18-2*n19+n20+2*n21-n22-2*n23+n24+2*n25-n26-2*n27+n28+2*n29-n30-2*n31
Congruence_by_7 = n0+2*n1+4*n2+n3+2*n4+4*n5+n6+2*n7+4*n8+n9+2*n10+4*n11+n12+2*n13+4*n14+n15+2*n16+4*n17+n18+2*n19+4*n20+n21+2*n22+4*n23+n24+2*n25+4*n26+n27+2*n28+4*n29+n30+2*n31
Congruence_by_8 = n&0x7
Congruence_by_9 = n0+2*n1+4*n2-n3-2*n4-4*n5+n6+2*n7+4*n8-n9-2*n10-4*n11+n12+2*n13+4*n14-n15-2*n16-4*n17+n18+2*n19+4*n20-n21-2*n22-4*n23+n24+2*n25+4*n26-n27-2*n28-4*n29+n30+2*n31
Congruence_by_11 = n0+2*n1+4*n2+8*n3+5*n4-n5-2*n6-4*n7-8*n8-5*n9+n10+2*n11+4*n12+8*n13+5*n14-n15-2*n16-4*n17-8*n18-5*n19+n20+2*n21+4*n22+8*n23+5*n24-n25-2*n26-4*n27-8*n28-5*n29+n30+2*n31
Congruence_by_13 = n0+2*n1+4*n2+8*n3+3*n4+6*n5-n6-2*n7-4*n8-8*n9-3*n10-6*n11+n12+2*n13+4*n14+8*n15+3*n16+6*n17-n18-2*n19-4*n20-8*n21-3*n22-6*n3+n24+2*n25+4*n26+8*n27+3*n28+6*n29-n30-2*n31
Congruence_by_16 = n&0xF

或者当因式分解时：

Congruence_by_1 = 0
Congruence_by_2 = n&0x1
Congruence_by_3 = (n0+n2+n4+n6+n8+n10+n12+n14+n16+n18+n20+n22+n24+n26+n28+n30)-(n1+n3+n5+n7+n9+n11+n13+n15+n17+n19+n21+n23+n25+n27+n29+n31)
Congruence_by_4 = n&0x3
Congruence_by_5 = n0+n4+n8+n12+n16+n20+n24+n28-(n2+n6+n10+n14+n18+n22+n26+n30)+2*(n1+n5+n9+n13+n17+n21+n25+n29-(n3+n7+n11+n15+n19+n23+n27+n31))
Congruence_by_7 = n0+n3+n6+n9+n12+n15+n18+n21+n24+n27+n30+2*(n1+n4+n7+n10+n13+n16+n19+n22+n25+n28+n31)+4*(n2+n5+n8+n11+n14+n17+n20+n23+n26+n29)
Congruence_by_8 = n&0x7
Congruence_by_9 = n0+n6+n12+n18+n24+n30-(n3+n9+n15+n21+n27)+2*(n1+n7+n13+n19+n25+n31-(n4+n10+n16+n22+n28))+4*(n2+n8+n14+n20+n26-(n5+n11+n17+n23+n29))
// and so on

如果这些值最终为负数，请将其与 i 添加直到它们变为正数。

现在你应该做的是通过我们刚才所做的相同过程递归地输入这些值，直到Congruence_by_i 变得小于i（很明显>= 0）。这类似于我们想要找到一个数字的余数 3 或 9 时所做的事情，记得吗？将数字相加，如果它有多个数字，则将结果的数字再次加起来，直到您只得到一个数字。

现在i = 1, 2, 3, 4, 5, 7, 8, 9, 11, 13, 16：

Divisible_by_i = (Congruence_by_i == 0);

剩下的：

Divisible_by_6 = Divisible_by_3 && Divisible_by_2;
Divisible_by_10 = Divisible_by_5 && Divisible_by_2;
Divisible_by_12 = Divisible_by_4 && Divisible_by_3;
Divisible_by_14 = Divisible_by_7 && Divisible_by_2;
Divisible_by_15 = Divisible_by_5 && Divisible_by_3;

编辑：请注意，有些添加可以从一开始就避免。例如n0+2*n1+4*n2 与n&0x7 相同，n3+2*n4+4*n5 与(n>>3)&0x7 相同，因此对于每个公式，您不必单独获取每个位，为了清晰和相似，我这样写在操作中。要优化每个公式，您应该自己处理它；分组操作数和分解操作。

Answer 15

一种可以帮助对所有整数值进行模减少的方法使用位切片和弹出计数。

mod3 = pop(x & 0x55555555) + pop(x & 0xaaaaaaaa) << 1;  // <- one term is shared!
mod5 = pop(x & 0x99999999) + pop(x & 0xaaaaaaaa) << 1 + pop(x & 0x44444444) << 2;
mod7 = pop(x & 0x49249249) + pop(x & 0x92492492) << 1 + pop(x & 0x24924924) << 2;
modB = pop(x & 0x5d1745d1) + pop(x & 0xba2e8ba2) << 1 + 
       pop(x & 0x294a5294) << 2 + pop(x & 0x0681a068) << 3;
modD = pop(x & 0x91b91b91) + pop(x & 0xb2cb2cb2) << 1 +
       pop(x & 0x64a64a64) << 2 + pop(x & 0xc85c85c8) << 3;

这些变量的最大值是 48、80、73、168 和 203，它们都适合 8 位变量。第二轮可以并行进行（也可以应用一些LUT方法）

      mod3 mod3 mod5 mod5 mod5 mod7 mod7 mod7 modB modB modB modB modD modD modD modD
mask  0x55 0xaa 0x99 0xaa 0x44 0x49 0x92 0x24 0xd1 0xa2 0x94 0x68 0x91 0xb2 0x64 0xc8
shift  *1   *2   *1   *2   *4   *1   *2   *4   *1   *2   *4   *8   *1   *2   *4   *8
sum   <-------> <------------> <----------->  <-----------------> <----------------->

Answer 16

假设number 是unsigned（32 位）。那么以下是计算高达 16 的可除性的非常快速的方法。（我没有测量，但汇编代码表明是这样。）

bool divisible_by_2 = number % 2 == 0;
bool divisible_by_3 = number * 2863311531u <= 1431655765u;
bool divisible_by_4 = number % 4 == 0;
bool divisible_by_5 = number * 3435973837u <= 858993459u;
bool divisible_by_6 = divisible_by_2 && divisible_by_3;
bool divisible_by_7 = number * 3067833783u <= 613566756u;
bool divisible_by_8 = number % 8 == 0;
bool divisible_by_9 = number * 954437177u <= 477218588u;
bool divisible_by_10 = divisible_by_2 && divisible_by_5;
bool divisible_by_11 = number * 3123612579u <= 390451572u;
bool divisible_by_12 = divisible_by_3 && divisible_by_4;
bool divisible_by_13 = number * 3303820997u <= 330382099u;
bool divisible_by_14 = divisible_by_2 && divisible_by_7;
bool divisible_by_15 = number * 4008636143u <= 286331153u;
bool divisible_by_16 = number % 16 == 0;

关于被d 整除，以下规则成立：

• 当d 是 2 的幂时：

作为pointed out by James Kanze，您可以使用is_divisible_by_d = (number % d == 0)。编译器足够聪明，可以将其实现为 (number & (d - 1)) == 0，这非常有效但被混淆了。

但是，当d 不是 2 的幂时，看起来上面显示的混淆比当前编译器所做的更有效。 （稍后会详细介绍）。

• d 为奇数时：

该技术采用is_divisible_by_d = number * a <= b 的形式，其中a 和b 是cleverly obtained constants。请注意，我们只需要 1 次乘法和 1 次比较：

• 当d 是偶数但不是 2 的幂时：

然后，写d = p * q，其中p是2的幂，q是奇数，并使用unpythonic建议的"tongue in cheek"，即is_divisible_by_d = is_divisible_by_p && is_divisible_by_q。同样，仅执行 1 次乘法（在计算 is_divisible_by_q 时）。

许多编译器（我使用godbolt 测试了clang 5.0.0、gcc 7.3、icc 18 和msvc 19）将number % d == 0 替换为(number / d) * d == number。他们使用一种巧妙的技术（参见Olof Forshell 的answer 中的参考资料）用乘法和位移来代替除法。他们最终做了2次乘法。相比之下，上述技术只执行 1 次乘法。

2018 年 10 月 1 日更新

看起来上面的算法很快就会出现在 GCC 中（已经在主干中）：

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=82853

GCC 的实现似乎更加高效。实际上，上面的实现有三个部分：1）被除数的偶数部分整除； 2) 被除数的奇数整除； 3）&&连接前两步的结果。通过使用在标准 C++ 中无法有效使用的汇编指令 (ror)，GCC 将这三个部分包装成一个单独的部分，这与奇数部分的可分性非常相似。好东西！有了这个实现，最好（为了清晰和性能）一直回退到%。

2020 年 5 月 5 日更新

我关于这个主题的文章已经发表：

Quick Modular Calculations (Part 1)，Overload Journal 154，2019 年 12 月，第 11-15 页。

Quick Modular Calculations (Part 2)，Overload Journal 155，2020 年 2 月，第 14-17 页。

Quick Modular Calculations (Part 3)，Overload Journal 156，2020 年 4 月，第 10-13 页。