给定素数N,计算下一个素数?
一位同事告诉我,C#字典集合由于与哈希有关的神秘原因而被素数调整。 而我当前的问题是,“它怎么知道下一个主要是什么?他们是一个巨大的桌面还是dynamic的计算?这是一个可怕的非确定性的运行时插入导致resize”
所以我的问题是,如果N是素数,那么计算下一个素数的最有效方法是什么?
连续素数之间的差距已知是相当小的,第一个差距超过100发生在素数370261。这意味着,即使是一个简单的蛮力在大多数情况下将足够快,以O(ln(p)* sqrt(p))。
对于p = 10,000这是O(921)操作。 考虑到我们将每O(ln(p))插入执行一次(粗略地说),这在大多数现代硬件上的大部分问题的约束下处于毫秒数量级的限制之内。
大约一年前,我在为C ++ 11实现无序(散列)容器的同时为libc ++工作。 我以为我会在这里分享我的经验。 这个经验支持marcog为“暴力”的合理定义所接受的答案 。
这意味着在大多数情况下,即使是一个简单的蛮力也足够快,平均取O(ln(p)* sqrt(p))。
我开发了几个size_t next_prime(size_t n)
的实现,这个函数的规范是:
返回:大于或等于
n
的最小素n
。
next_prime
每个实现都伴随着一个辅助函数is_prime
。 is_prime
应该被认为是一个私有的实现细节; 不意味着被客户直接调用。 当然每一个实现都是正确的,但是也要testing下面的性能testing:
int main() { typedef std::chrono::high_resolution_clock Clock; typedef std::chrono::duration<double, std::milli> ms; Clock::time_point t0 = Clock::now(); std::size_t n = 100000000; std::size_t e = 100000; for (std::size_t i = 0; i < e; ++i) n = next_prime(n+1); Clock::time_point t1 = Clock::now(); std::cout << e/ms(t1-t0).count() << " primes/millisecond\n"; return n; }
我应该强调,这是一个性能testing,并不反映典型的使用情况,看起来更像是:
// Overflow checking not shown for clarity purposes n = next_prime(2*n + 1);
所有的性能testing都是用
clang++ -stdlib=libc++ -O3 main.cpp
实施1
有七个实现。 显示第一个实现的目的是为了certificate如果你没有停止testingsqrt(x)
的候选素数x
,那么你甚至未能实现一个可以被分类为powershell的实现。 这个实现是非常慢的 。
bool is_prime(std::size_t x) { if (x < 2) return false; for (std::size_t i = 2; i < x; ++i) { if (x % i == 0) return false; } return true; } std::size_t next_prime(std::size_t x) { for (; !is_prime(x); ++x) ; return x; }
对于这个实现,我只需要将e
设置为100而不是100000,只是为了获得合理的运行时间:
0.0015282 primes/millisecond
实施2
这个实现是蛮力实现中最慢的一个,与实现1唯一的不同之处在于,当因子超过sqrt(x)
时,它停止testing素数。
bool is_prime(std::size_t x) { if (x < 2) return false; for (std::size_t i = 2; true; ++i) { std::size_t q = x / i; if (q < i) return true; if (x % i == 0) return false; } return true; } std::size_t next_prime(std::size_t x) { for (; !is_prime(x); ++x) ; return x; }
请注意, sqrt(x)
不是直接计算的,而是由q < i
推断的。 这加速了几千倍:
5.98576 primes/millisecond
并validation了马cocoa的预测:
在大多数现代硬件上,大部分问题的时间约为1毫秒。
实施3
通过避免使用%
运算符,可以将速度提高一倍(至less在我使用的硬件上):
bool is_prime(std::size_t x) { if (x < 2) return false; for (std::size_t i = 2; true; ++i) { std::size_t q = x / i; if (q < i) return true; if (x == q * i) return false; } return true; } std::size_t next_prime(std::size_t x) { for (; !is_prime(x); ++x) ; return x; } 11.0512 primes/millisecond
实施4
到目前为止,我甚至没有使用2是唯一的偶数的常识。 这个实施结合了这种知识,再次提高了将近一倍的速度:
bool is_prime(std::size_t x) { for (std::size_t i = 3; true; i += 2) { std::size_t q = x / i; if (q < i) return true; if (x == q * i) return false; } return true; } std::size_t next_prime(std::size_t x) { if (x <= 2) return 2; if (!(x & 1)) ++x; for (; !is_prime(x); x += 2) ; return x; } 21.9846 primes/millisecond
实施4可能是大多数人在思考“蛮力”时想到的。
实施5
使用下面的公式可以很容易地select所有可以被2和3整除的数字:
6 * k + {1, 5}
其中k> = 1。下面的实现使用这个公式,但用一个可爱的xor技巧来实现:
bool is_prime(std::size_t x) { std::size_t o = 4; for (std::size_t i = 5; true; i += o) { std::size_t q = x / i; if (q < i) return true; if (x == q * i) return false; o ^= 6; } return true; } std::size_t next_prime(std::size_t x) { switch (x) { case 0: case 1: case 2: return 2; case 3: return 3; case 4: case 5: return 5; } std::size_t k = x / 6; std::size_t i = x - 6 * k; std::size_t o = i < 2 ? 1 : 5; x = 6 * k + o; for (i = (3 + o) / 2; !is_prime(x); x += i) i ^= 6; return x; }
这实际上意味着algorithm只需检查1/3的整数而不是1/2的数字,性能testing显示预期的加速将近50%:
32.6167 primes/millisecond
实施6
这个实现是实现5的逻辑扩展:它使用下面的公式来计算不能被2,3和5整除的所有数字:
30 * k + {1, 7, 11, 13, 17, 19, 23, 29}
它还展开了is_prime中的内部循环,并创build了一个“小素数”列表,该列表对于处理小于30的数字非常有用。
static const std::size_t small_primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 }; static const std::size_t indices[] = { 1, 7, 11, 13, 17, 19, 23, 29 }; bool is_prime(std::size_t x) { const size_t N = sizeof(small_primes) / sizeof(small_primes[0]); for (std::size_t i = 3; i < N; ++i) { const std::size_t p = small_primes[i]; const std::size_t q = x / p; if (q < p) return true; if (x == q * p) return false; } for (std::size_t i = 31; true;) { std::size_t q = x / i; if (q < i) return true; if (x == q * i) return false; i += 6; q = x / i; if (q < i) return true; if (x == q * i) return false; i += 4; q = x / i; if (q < i) return true; if (x == q * i) return false; i += 2; q = x / i; if (q < i) return true; if (x == q * i) return false; i += 4; q = x / i; if (q < i) return true; if (x == q * i) return false; i += 2; q = x / i; if (q < i) return true; if (x == q * i) return false; i += 4; q = x / i; if (q < i) return true; if (x == q * i) return false; i += 6; q = x / i; if (q < i) return true; if (x == q * i) return false; i += 2; } return true; } std::size_t next_prime(std::size_t n) { const size_t L = 30; const size_t N = sizeof(small_primes) / sizeof(small_primes[0]); // If n is small enough, search in small_primes if (n <= small_primes[N-1]) return *std::lower_bound(small_primes, small_primes + N, n); // Else n > largest small_primes // Start searching list of potential primes: L * k0 + indices[in] const size_t M = sizeof(indices) / sizeof(indices[0]); // Select first potential prime >= n // Known a-priori n >= L size_t k0 = n / L; size_t in = std::lower_bound(indices, indices + M, n - k0 * L) - indices; n = L * k0 + indices[in]; while (!is_prime(n)) { if (++in == M) { ++k0; in = 0; } n = L * k0 + indices[in]; } return n; }
这可以说是超越了“蛮力”,有利于提高27.5%的速度:
41.6026 primes/millisecond
实施7
把上面的游戏重复一次,开发一个不能被2,3,5和7整除的数字公式是很实际的:
210 * k + {1, 11, ...},
这里没有显示源代码,但与实现6非常相似。这是我select实际用于libc ++的无序容器的实现,并且源代码是开放源代码(可在链接中find)。
这个最后的迭代对于另一个14.6%的速度提升是有好处的:
47.685 primes/millisecond
使用这个algorithm可以确保libc ++的哈希表的客户端可以select他们认为对他们的情况最有利的任何主数据,而且这个应用程序的性能是完全可以接受的。
以防万一有人好奇:
使用reflection器我确定.Net使用一个静态类,其中包含一个72个素数的硬编码列表,其中至less有7199369个扫描最小的素数,至less是当前大小的两倍,而对于大于它计算的大小通过对所有奇数进行试划分,直到潜在数字的sqrt。 这个类是不可变的,线程安全的(即更大的素数不会被存储以供将来使用)。
只是几个实验与素间距离。
这是可视化其他答案的补充。
我从第100.000(= 1,299,709)到第200,000(= 2,750,159)
一些数据:
Maximum interprime distance = 148 Mean interprime distance = 15
相距距离频率图:
相距距离与素数
只是看它是“随机的”。 不过…
一个很好的窍门是使用部分筛子。 例如,N = 2534536543556之后的下一个素数是多less?
检查相对于小素数列表的N的模数。 从而…
mod(2534536543556,[3 5 7 11 13 17 19 23 29 31 37]) ans = 2 1 3 6 4 1 3 4 22 16 25
我们知道N之后的下一个素数必须是一个奇数,我们可以立即抛弃这个小素数列表的所有奇数倍。 这些模数允许我们筛出这些小素数的倍数。 如果我们使用200以下的小质数,我们可以使用这种scheme立即丢弃大于N的大多数潜在素数,除了一个小列表。
更明确地说,如果我们正在寻找一个2534536543556以上的素数,它不能被2整除,所以我们只需要考虑超出这个值的奇数。 上面的模数表明2534536543556与2模3是一致的,因此2534536543556 + 1与0模3一致,必须是2534536543556 + 7,2534536543556 + 13等。实际上,我们可以筛选出许多数字而不需要任何需要要testing他们的素质和没有任何审判分裂。
同样的事实
mod(2534536543556,7) = 3
告诉我们2534536543556 + 4与0 mod 7一致。当然,这个数字是偶数,所以我们可以忽略它。 但是2534536543556 + 11是一个可以被7整除的奇数,2534536543556 + 25等也是可以被整除的。再次,我们可以把这些数字清除为复合的(因为它们可以被7整除),所以不是整数。
只使用最多37个素数的小列表,我们可以排除紧随2534536543556之后的大部分数字,除了less数:
{2534536543573 , 2534536543579 , 2534536543597}
在这些数字中,它们是否是最好的?
2534536543573 = 1430239 * 1772107 2534536543579 = 99833 * 25387763
我已经做了提供列表中前两个数字的素数分解的努力。 看到它们是复合的,但主要的因素很大。 当然,这是有道理的,因为我们已经确保没有剩下的数字可以有小的主要因素。 在我们的短名单(2534536543597)中的第三个实际上是超过N的第一个素数。我所描述的筛选scheme往往会导致数字是否为素数,或者是由大量素数因子组成。 所以我们需要在find下一个素数之前实际应用一个明确的素数testing,只有几个数字。
类似的scheme很快就会产生超过N = 100000000000000000000000的下一个素数,如100000000000000000000000010。
没有函数f(n)来计算下一个素数。 相反,一个数字必须经过素数testing。
当find第n个素数时,从第1个到第(n-1)个已知所有素数也是非常有用的,因为那些是唯一需要testing的数字。
由于这些原因,如果有一组预先计算的大素数,我不会感到惊讶。 如果需要重复计算某些素数,对我来说没有意义。
正如其他人已经指出的那样,目前尚未findfind下一个素数的方法。 因此,大多数algorithm更注重于使用快速检查素数的方法,因为您必须检查已知素数和下一个素数之间的n / 2。
根据应用程序的不同,如Paul Wheeler所指出的那样,您也可以通过对查找表进行硬编码来避开。
纯粹的新奇,总是有这种方法:
#!/usr/bin/perl for $p ( 2 .. 200 ) { next if (1x$p) =~ /^(11+)\1+$/; for ($n=1x(1+$p); $n =~ /^(11+)\1+$/; $n.=1) { } printf "next prime after %d is %d\n", $p, length($n); }
这当然会产生
next prime after 2 is 3 next prime after 3 is 5 next prime after 5 is 7 next prime after 7 is 11 next prime after 11 is 13 next prime after 13 is 17 next prime after 17 is 19 next prime after 19 is 23 next prime after 23 is 29 next prime after 29 is 31 next prime after 31 is 37 next prime after 37 is 41 next prime after 41 is 43 next prime after 43 is 47 next prime after 47 is 53 next prime after 53 is 59 next prime after 59 is 61 next prime after 61 is 67 next prime after 67 is 71 next prime after 71 is 73 next prime after 73 is 79 next prime after 79 is 83 next prime after 83 is 89 next prime after 89 is 97 next prime after 97 is 101 next prime after 101 is 103 next prime after 103 is 107 next prime after 107 is 109 next prime after 109 is 113 next prime after 113 is 127 next prime after 127 is 131 next prime after 131 is 137 next prime after 137 is 139 next prime after 139 is 149 next prime after 149 is 151 next prime after 151 is 157 next prime after 157 is 163 next prime after 163 is 167 next prime after 167 is 173 next prime after 173 is 179 next prime after 179 is 181 next prime after 181 is 191 next prime after 191 is 193 next prime after 193 is 197 next prime after 197 is 199 next prime after 199 is 211
除了所有的趣味和游戏之外,众所周知,最佳哈希表大小严格certificate是4N−1
forms的素数。 所以find下一个素数是不够的。 你也必须做其他的检查。
据我记得,它使用当前大小的两倍旁边的素数。 它不会计算素数 – 表中预加载的数字达到一个很大的值(不完全是大约10,000,000左右)。 当达到这个数字时,它使用一些天真的algorithm来获得下一个数字(如curNum = curNum + 1),并使用一些validation它,如果这些方法: http : //en.wikipedia.org/wiki/Prime_number#Verifying_primality