试验部门对素性的条件测试

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); for(p = 3; p <= cut; p + = 2)

案例1:适用于所有n但它必须在每次迭代时进行额外的除法(编辑:实际上它不需要额外的除法但它仍然较慢.我不知道为什么.请参阅下面的汇编输出).我发现它的速度是情况2的两倍,因为大的n值是素数(在我的Sandy Bridge系统上).

案例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)是否足够准确.

有人可以向我解释为什么案例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上的速度大约是两倍.我不知道为什么.

情况1:

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

案例2:

.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计算它们，以向编译器显示我们想要商和余数。

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;
}

我在 MinGW-w64 编译器上尝试了http://coliru.stacked-crooked.com/a/69497863a97d8953中的代码，case 1速度比case 2.

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

如果我猜对了，那就是有原因的。

在 32 位架构上，64 位数字在两个 32 位寄存器中表示，您的编译器将完成所有串联工作。进行 64 位加法、减法和乘法很简单。但求模和除法是通过一个小函数调用来完成的，在 GCC和MSVC 中，该函数调用名为___umoddi3and 。___udivdi3aullremaulldiv

因此，实际上，中的每次迭代都需要一___umoddi3加一，中的 64 位乘法的一加一串联。这就是为什么看起来比你的测试慢两倍的原因。___udivdi3case 1___udivdi3case 2case 1case 2

你真正得到的case 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

你真正得到的case 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