通过试验划分进行素数条件测试

Question

我的问题是关于试用部门的条件测试。关于采用什么条件测试似乎存在一些争论。让我们看看 RosettaCode 中的代码.

int is_prime(unsigned int n)
{
    unsigned int p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    for (p = 3; p <= n/p; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

轮分解或使用预先确定的素数列表不会改变我问题的本质。

条件测试我能想到的三种情况:

1. p<=n/p
2. p*p<=n
3. int cut = sqrt(n);对于(p = 3；p <= 剪切；p += 2)

案例 1:适用于所有 n 但它必须在每次迭代时进行额外的除法(编辑:实际上它不需要额外的除法但它仍然较慢。我不确定为什么。请参阅汇编输出下)。我发现对于较大的 n 值(在我的 Sandy Bridge 系统上)，它比情况 2 慢两倍。

案例 2:明显比案例 1 快，但它有一个问题，即它会溢出大 n 并进入无限循环。它可以处理的最大值是

(sqrt(n) + c)^2 = INT_MAX  //solve
n = INT_MAX -2*c*sqrt(INT_MAX) + c^2
//INT_MAX = 2^32 -> n = 2^32 - c*s^17 + c^2; in our case c = 2

例如，对于 uint64_t，情况 2 将进入无限循环，因为 x =-1L-58 (x^64-59) 是素数。

情况 3:每次迭代都不需要进行除法或乘法运算，并且不会像情况 2 那样溢出。它也比情况 2 稍快一些。唯一的问题是如果sqrt(n) is accurate enough .

有人可以向我解释为什么情况 2 比情况 1 快得多吗？情况 1 没有像我一样使用额外的除法，但尽管如此它仍然慢很多。

这是素数 2^56-5 的次数；

case 1 9.0s
case 2 4.6s
case 3 4.5s

这是我用来测试 http://coliru.stacked-crooked.com/a/69497863a97d8953 的代码.我还在这个问题的末尾添加了函数。

这里是 GCC 4.8 的汇编输出形式，其中 -O3 用于情况 1 和情况 2。它们都只有一个除法。案例 2 也有一个乘法，所以我的第一个猜测是案例 2 会更慢，但它在 GCC 和 MSVC 上大约快两倍。我不知道为什么。

案例一:

.L5:
  testl %edx, %edx
  je  .L8
.L4:
  addl  $2, %ecx
  xorl  %edx, %edx
  movl  %edi, %eax
  divl  %ecx
  cmpl  %ecx, %eax
  jae .L5

案例二:

.L20:
  xorl  %edx, %edx
  movl  %edi, %eax
  divl  %ecx
  testl %edx, %edx
  je  .L23
.L19:
  addl  $2, %ecx
  movl  %ecx, %eax
  imull %ecx, %eax
  cmpl  %eax, %edi
  jae .L20

以下是我用来测试时间的函数:

int is_prime(uint64_t n)
{
    uint64_t p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    for (p = 3; p <= n/p; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

int is_prime2(uint64_t n)
{
    uint64_t p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    for (p = 3; p*p <= n; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

int is_prime3(uint64_t n)
{
    uint64_t p;
    if (!(n & 1) || n < 2 ) return n == 2;

    /* comparing p*p <= n can overflow */
    uint32_t cut = sqrt(n);
    for (p = 3; p <= cut; p += 2)
        if (!(n % p)) return 0;
    return 1;
}

赏金之后添加的内容。

Aean 发现案例 1 中保存商和余数的速度与案例 2 一样快(或稍快)。我们称此案例为 4。以下代码的速度是案例 1 的两倍。

int is_prime4(uint64_t n)
{
    uint64_t p, q, r;
    if (!(n & 1) || n < 2 ) return n == 2;

    for (p = 3, q=n/p, r=n%p; p <= q; p += 2, q = n/p, r=n%p)
        if (!r) return 0;
    return 1;
}

我不确定为什么这有帮助。在任何情况下都不再需要使用案例 2。对于案例 3，大多数版本的硬件或软件中的 sqrt 函数都能正确生成完美的正方形，因此通常可以安全使用。案例 3 是唯一适用于 OpenMP 的案例。

Answer 1

UPD:这显然是一个编译器优化问题。虽然 MinGW 在循环体中只使用了一条 div 指令，但 Linux 上的 GCC 和 MSVC 都未能重用前一次迭代中的商。

我认为我们能做的最好的事情就是显式定义 quo 和 rem 并在相同的基本指令 block 中计算它们，以向编译器显示我们需要商和余数.

int is_prime(uint64_t n)
{
    uint64_t p = 3, quo, rem;
    if (!(n & 1) || n < 2) return n == 2;

    quo = n / p;
    for (; p <= quo; p += 2){
        quo = n / p; rem = n % p;
        if (!(rem)) return 0;
    }
    return 1;
}

---

我试过你来自 http://coliru.stacked-crooked.com/a/69497863a97d8953 的代码在 MinGW-w64 编译器上，case 1 比 case 2 更快。​​

所以我猜测您正在针对 32 位架构进行编译并使用 uint64_t 类型。您的程序集显示它不使用任何 64 位寄存器。

如果我做对了，那就是原因。

在 32 位架构上，64 位数字在两个 32 位寄存器中表示，您的编译器将完成所有连接工作。做 64 位加法、减法和乘法很简单。但是模和除法是通过一个小函数调用完成的，在 GCC 中名为 ___umoddi3 和 ___udivdi3，aullrem 和 aulldiv 在 MSVC 中。

所以实际上你需要一个 ___umoddi3 和一个 ___udivdi3 用于case 1 中的每个迭代，一个 ___udivdi3 和案例 2 中 64 位乘法的一个串联。这就是为什么在您的测试中 case 1 似乎比 case 2 慢两倍。

在案例 1 中您真正得到的是:

L5:
    addl    $2, %esi
    adcl    $0, %edi
    movl    %esi, 8(%esp)
    movl    %edi, 12(%esp)
    movl    %ebx, (%esp)
    movl    %ebp, 4(%esp)
    call    ___udivdi3         // A call for div
    cmpl    %edi, %edx
    ja  L6
    jae L21
L6:
    movl    %esi, 8(%esp)
    movl    %edi, 12(%esp)
    movl    %ebx, (%esp)
    movl    %ebp, 4(%esp)
    call    ___umoddi3        // A call for modulo.
    orl %eax, %edx
    jne L5

在案例 2 中您真正得到的是:

L26:
    addl    $2, %esi
    adcl    $0, %edi
    movl    %esi, %eax
    movl    %edi, %ecx
    imull   %esi, %ecx
    mull    %esi
    addl    %ecx, %ecx
    addl    %ecx, %edx
    cmpl    %edx, %ebx
    ja  L27
    jae L41
L27:
    movl    %esi, 8(%esp)
    movl    %edi, 12(%esp)
    movl    %ebp, (%esp)
    movl    %ebx, 4(%esp)
    call    ___umoddi3         // Just one call for modulo
    orl %eax, %edx
    jne L26

MSVC 未能重用 div 的结果。 return 破坏了优化。试试这些代码:

__declspec(noinline) int is_prime_A(unsigned int n)
{
    unsigned int p;
    int ret = -1;
    if (!(n & 1) || n < 2) return n == 2;

    /* comparing p*p <= n can overflow */
    p = 1;
    do {
        p += 2;
        if (p >= n / p) ret = 1; /* Let's return latter outside the loop. */
        if (!(n % p)) ret = 0;
    } while (ret < 0);
    return ret;
}

__declspec(noinline) int is_prime_B(unsigned int n)
{
    unsigned int p;
    if (!(n & 1) || n < 2) return n == 2;

    /* comparing p*p <= n can overflow */
    p = 1;
    do {
        p += 2;
        if (p > n / p) return 1; /* The common routine. */
        if (!(n % p)) return 0;
    } while (1);
}

在 Windows 的 MSVC/ICC 上，is_prime_B 将比 is_prime_A 慢两倍。