一位同事告诉我,C# Dictionary集合因与哈希相关的神秘原因而通过质数调整大小。我的第一个问题是,“它如何知道下一个质数是什么?他们是否存储了一个巨大的表格或即时计算?这是插入导致调整大小的可怕的非确定性运行时。”
因此,我的问题是,鉴于N是一个质数,最有效的方法是计算下一个质数是什么?
一位同事告诉我,C# Dictionary集合因与哈希相关的神秘原因而通过质数调整大小。我的第一个问题是,“它如何知道下一个质数是什么?他们是否存储了一个巨大的表格或即时计算?这是插入导致调整大小的可怕的非确定性运行时。”
因此,我的问题是,鉴于N是一个质数,最有效的方法是计算下一个质数是什么?
我开发了几个这意味着即使是简单的暴力,在大多数情况下也足够快,平均需要O(ln(p)*sqrt(p))。
size_t next_prime(size_t n)
的实现,其中该函数的规范为:
每个返回: 大于或等于
n
的最小质数。
next_prime
的实现都伴随着一个辅助函数is_prime
。应将is_prime
视为私有实现细节;不应直接由客户端调用。当然,每个实现都经过了正确性测试,但也进行了以下性能测试: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;
}
我应该强调这是一项性能测试,并不代表典型用法,后者可能更像:
// Overflow checking not shown for clarity purposes
n = next_prime(2*n + 1);
所有性能测试都使用以下编译:
clang++ -stdlib=libc++ -O3 main.cpp
实现方法1
共有七种实现方法。展示第一种实现方法的目的是为了证明,如果你在判断候选质数x
是否有因子时没有在sqrt(x)
处停止测试,那么你甚至都没有达到可以被归类为暴力破解的实现方法。这种实现方法极其缓慢。
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
这种实现方式是最慢的"暴力"实现之一,唯一的区别在于当因子超过sqrt(x)
时停止测试素性。
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
并验证了marcog的预测:
...这远远在大多数问题的约一毫秒级别的硬件上操作限制范围内。
实现3
可以通过避免使用%
运算符来将速度几乎提高一倍(至少在我所使用的硬件上):
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
使用以下公式,您可以轻松选择所有既不被2整除也不被3整除的数字:
6 * k + {1, 5}
当 k >= 1 时。以下实现使用该公式,但是采用了一种可爱的异或技巧来实现:
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;
}
这实际上意味着算法只需要检查1/3的整数是否为质数,而不是1/2的数字,性能测试显示预期加速近50%:
32.6167 primes/millisecond
实现6
这个实现是实现5的一个逻辑延伸:它使用以下公式来计算所有不能被2、3和5整除的数字:
30 * k + {1, 7, 11, 13, 17, 19, 23, 29}
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;
}
41.6026 primes/millisecond
实现 7
对于不被2、3、5和7整除的数字,继续进行上述游戏一个迭代是很实际的,可以发展出一个公式:
210 * k + {1, 11, ...},
这里没有展示源代码,但它与第6种实现非常相似。 这是我选择在libc ++的无序容器中实际使用的实现, 该源代码是开源的(在链接中找到)。
这个最终版本可以再获得14.6%的速度提升:
47.685 primes/millisecond
使用该算法可以确保libc++哈希表的客户端可以选择任何他们认为最有利的质数,而且在这种情况下,性能是相当可接受的。
万一有人好奇:
使用反射器,我确定.NET使用一个静态类,其中包含一个硬编码列表,其中包含了范围在7199369以内的大约72个素数,该列表扫描寻找最小的至少是当前大小两倍的质数,并且对于大于此大小的大小,它通过试除法计算所有奇数直到潜在数字的平方根来计算下一个质数。这个类是不可变的和线程安全的(即不会为将来使用存储更大的质数)。
连续质数之间的差距通常很小,第一个超过100的差距出现在质数370261。这意味着即使是简单粗暴的方法,在大多数情况下也足够快,平均需要O(ln(p)*sqrt(p))的时间。
对于p=10,000,这是O(921)的操作。请注意,考虑到我们将每隔O(ln(p))插入一次(粗略地说),这远远在大多数现代硬件上花费约1毫秒的问题限制范围内。
FILL_FACTOR * n
个条目复制并重新哈希到新数组中,实际上速度非常快。 - maaartinus一个好的技巧是使用部分筛法。例如,接下来的质数是多少,它跟数字N = 2534536543556有关呢?
通过与一些小质数的模数检查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
我们知道接下来的素数必须是奇数,因此我们可以立即排除这些小素数列表中的所有奇数倍数。这些模数允许我们筛选出那些小质数的倍数。如果我们使用小于200的小质数,我们可以使用这个方案立即丢弃大于N的大多数潜在素数,除了一个小列表。
更明确地说,如果我们要寻找2534536543556之后的质数,它不能被2整除,因此我们只需要考虑该值之后的奇数。上面的模数表明2534536543556同余于2模3,因此2534536543556+1同余于0模3,而2534536543556+7、2534536543556+13等也都同余于0模3。实际上,我们可以筛选出许多数字,而无需测试它们是否为素数或进行任何试除法。
同样,事实上
mod(2534536543556,7) = 3
告诉我们,2534536543556 + 4 模 7 同余于 0。当然,那个数是偶数,所以我们可以忽略它。但是,2534536543556 + 11 是一个奇数,可以被 7 整除,同样的,2534536543556 + 25 也可以被整除,等等。同样地,我们可以排除这些数字,因为它们明显是合数(因为它们可以被 7 整除),而不是质数。
仅使用小于 37 的质数列表,我们可以排除大部分紧随我们起始点 2534536543556 的数字,只有少数例外:
{2534536543573 , 2534536543579 , 2534536543597}
这些数字中,有质数吗?
2534536543573 = 1430239 * 1772107
2534536543579 = 99833 * 25387763
我已经努力提供了列表中前两个数字的质因数分解。请注意它们是合数,但质因数很大。当然,这是有道理的,因为我们已经确保剩下的数字没有小的质因数。在我们短列表中的第三个数字(2534536543597)实际上是超过N的第一个质数。我描述的筛法倾向于产生要么是质数,要么是由通常很大的质因数组成的数字。因此,在找到下一个质数之前,我们需要对只有少数数字进行显式的素性测试。
类似的方案迅速得出了超过N = 1000000000000000000000000000的下一个质数,即1000000000000000000000000103。
这是为了可视化其他答案而进行的补充。
我从第100,000个(=1,299,709)到第200,000个(=2,750,159)获取了质数。
一些数据:
Maximum interprime distance = 148
Mean interprime distance = 15
素数间距的频率图:
素数间距与素数之间的关系:
仅仅为了看它是否“随机”。不过......
没有函数 f(n) 用来计算下一个质数,必须测试一个数字是否为质数。
当查找第 n 个质数时,已知从第 1 个到第 (n-1) 个所有质数非常有用,因为只需要测试这些数字作为因子。
由于这些原因,如果存在预先计算的大质数集合,我不会感到惊讶。对我而言,某些质数需要反复重新计算似乎没有太多意义。
如果要追求新奇,总有这种方法:
#!/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
4N-1
的质数。 因此,仅找到下一个质数是不够的。 您还需要进行其他检查。next_prime = prime + 2
。你可能是对的,而且没有人能证明一旦你足够高,你总是会错。所以去试试吧。 :) - tchrist