试验分裂对素数的条件测试
我的问题是关于试验部门的条件测试。 关于采用什么条件测试似乎存在争议。 我们来看看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; } 车轮分解或使用预定的素数列表不会改变我的问题的本质。
我可以想到三种情况来进行条件测试：
- P <= N / P
- P * P <= N
- 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的代码。 我还在这个问题的最后添加了这些function。
下面是assembly输出formsGCC 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，大多数版本的sqrtfunction在硬件或软件中都可以获得完美的正方形，因此一般情况下使用它是安全的。 案例3是唯一适用于OpenMP的案例。
UPD：显然这是一个编译器优化问题。 虽然MinGW在循环体中只使用了一个div指令，但Linux和MSVC上的GCC都无法重用上一次迭代中的商。
我认为我们能做的最好是明确定义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中命名为___umoddi3和___udivdi3 ，在MSVC中命名为___umoddi3和___udivdi3 。
实际上，在case 1 ，每次迭代需要一个___umoddi3和一个___udivdi3 ，在case 1需要一个___umoddi3和一个64位乘法串联。 这就是为什么case 1在测试中看起来比case 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， is_prime_A MSVC / ICC上的is_prime_B慢两倍。
sqrt(n)只要你的sqrt单调递增就足够准确，它可以得到完美的正方形，并且每个unsigned int都可以完全表示为double 。 在我所知道的每个平台上都是这三种情况。
您可以通过实现函数unsigned int sqrti(unsigned int n)来解决这些问题（如果您认为它们是问题），该函数使用Newton的方法返回unsigned int平方根的底限。 （如果你以前从未这样做过，这是一个有趣的练习！）
回答这篇文章的一小部分。
案例2修复以处理溢出。
#include int is_prime(unsigned n) { unsigned p; if (!(n & 1) || n < 2) return n == 2; #define UINT_MAX_SQRT (UINT_MAX >> (sizeof(unsigned)*CHAR_BIT/2)) unsigned limit = n; if (n >= UINT_MAX_SQRT * UINT_MAX_SQRT) limit = UINT_MAX_SQRT * UINT_MAX_SQRT - 1; for (p = 3; p * p < limit; p += 2) if (!(n % p)) return 0; if (n != limit) if (!(n % p)) return 0; return 1; }
如果sizeof(unsigned)和CHAR_BIT都是奇数，则极限计算失败 - 这是一种罕见的情况。
关于你的第一个问题：为什么（2）比（1）更快？
嗯，这取决于编译器，也许。
然而，通常可以预期分割是比乘法更昂贵的操作。
关于你的第二个问题：sqrt（）是一个准确的函数吗？
一般来说，这是准确的。
可能给你带来问题的唯一情况是sqrt(n)是一个整数。
例如，如果你的系统中n == 9和sqrt(n) == 2.9999999999999 ，那么你就麻烦了，因为整数部分是2，但确切的值是3。
然而，这种罕见的情况很容易通过增加一个不那么小的双常数来处理，比如说。
因此，您可以写：
double stop = sqrt(n) + 0.1; for (unsigned int d = 2; d <= stop; d += 2) if (n % d == 0) break; /* not prime!! */
添加的术语0.1可以为您的算法添加一次迭代，这根本不是一个大问题。
最后，你的算法的明显选择是（3），即sqrt()方法，因为没有任何计算（乘法或除法），并且值stop只计算一次。
您可以获得的另一项改进如下：
- 注意，每个素数
p >= 5的forms为6n - 1或6n + 1。
因此，您可以将变量d的增量替换为2,4,2,4等。