给定素数N,计算下一个素数?

72

一位同事告诉我,C# Dictionary集合因与哈希相关的神秘原因而通过质数调整大小。我的第一个问题是,“它如何知道下一个质数是什么?他们是否存储了一个巨大的表格或即时计算?这是插入导致调整大小的可怕的非确定性运行时。”

因此,我的问题是,鉴于N是一个质数,最有效的方法是计算下一个质数是什么?


3
这个问题应该放在数学溢出网站上讨论。 - Kirk Broadhurst
2
也许你的同事是错误的,或者它使用了一些预先计算的质数而不是找到下一个质数。 - President James K. Polk
20
@Kirk:我不同意——这是一个算法问题,而不是数学问题。 - Jim Lewis
4
这段话的意思是:“这都属于理论计算机科学范畴,它处于编程和数学的交叉点。因此,我认为在任何一个网站上发布这个问题都没有问题。” - moinudin
11
@Kirk:这绝对不适合在MathOverflow上发布,该网站仅适用于“研究级问题”。我也不同意它需要放在http://math.stackexchange.com上,但那里至少也是适当的。 - Antal Spector-Zabusky
显示剩余4条评论
9个回答

79
大约一年前,我在实现C++11的无序(哈希)容器时,为libc ++工作。在这里,我想分享我的经验。这个经验支持marcog的采纳答案对“暴力”的“合理”定义。

这意味着即使是简单的暴力,在大多数情况下也足够快,平均需要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}

它还展开了is_prime内部循环,并创建了一个“小质数”列表,用于处理小于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种实现非常相似。 这是我选择在libc ++的无序容器中实际使用的实现, 该源代码是开源的(在链接中找到)。

这个最终版本可以再获得14.6%的速度提升:

47.685 primes/millisecond

使用该算法可以确保libc++哈希表的客户端可以选择任何他们认为最有利的质数,而且在这种情况下,性能是相当可接受的。


5
由于几乎所有CPU架构上的除法指令都会产生余数和商,因此实现3比实现2的速度快一倍表明你的编译器没有很好地进行优化,而是在实现2中使用了两个除法指令。在GCC 4.7.1上,实现2确实比实现3快4%,因为不需要额外的乘法。如果您交换两个if语句,则您的编译器可能会做得更好。实际上,“返回false”情况更有可能发生,因此仅出于这个原因就值得切换(提高1%的速度)。 - han
1
O(ln(p)sqrt(p))是完全偏离的。O(sqrt(p))是试除法测试质数的复杂度,它们在O(ln(p))步骤中间有间隔,但是所有这些数字都是合成数*,而不是质数。如果我正确地阅读了这个,那么平均复杂度为O(ln(ln(p))*sqrt(p)/ln(p))。 - Will Ness
2
@WillNess:我只是引用了被接受的答案。确切的复杂度对于我的回答来说是一个次要问题,我的回答存在的目的是展示各种实现策略对性能的影响。 - Howard Hinnant

44

万一有人好奇:

使用反射器,我确定.NET使用一个静态类,其中包含一个硬编码列表,其中包含了范围在7199369以内的大约72个素数,该列表扫描寻找最小的至少是当前大小两倍的质数,并且对于大于此大小的大小,它通过试除法计算所有奇数直到潜在数字的平方根来计算下一个质数。这个类是不可变的和线程安全的(即不会为将来使用存储更大的质数)。


好的回答。检查每个奇数并不是确定质数最有效的方法,但大多数字典可能包含不到350万个键。 - Kirk Broadhurst
2
STLPort 也使用了一个存储列表:really-ugly-long-gitweb-url - Marcus Borkenhagen
3
我忘了提到它首先测试是否可被2整除,尽管它似乎只从奇数开始尝试,因为在查找下一个质数时会增加2。然而,还有一种奇怪的情况,就是如果你调用HashHelpers.GetPrime(7199370),它将循环遍历从7199370到2147483646的所有偶数而找不到一个质数,然后返回7199370。这很有趣,但当然这是内部使用的,所以可能永远不会以这种方式调用。 - Paul Wheeler
3
抱歉,我必须在此澄清一个错误。原来输入值为偶数的情况下,有一种巧妙的按位或运算(OR)与数字1相结合的方法,可以将其转换成比它大的下一个奇数。 - Paul Wheeler
1
不过,这并不是我问题的答案。 - John Shedletsky

35

连续质数之间的差距通常很小,第一个超过100的差距出现在质数370261。这意味着即使是简单粗暴的方法,在大多数情况下也足够快,平均需要O(ln(p)*sqrt(p))的时间。

对于p=10,000,这是O(921)的操作。请注意,考虑到我们将每隔O(ln(p))插入一次(粗略地说),这远远在大多数现代硬件上花费约1毫秒的问题限制范围内。


同意复杂度并不过高,但每个操作都是相对较重的余数检查;并且检查本身的复杂度随 p 的增加而增加。 - Kirk Broadhurst
除非我还在睡觉,否则对于 p=10000,ln(p)=9.2且sqrt(p)=100,则 O(920)。 - Kirk Broadhurst
@GregS 相比将所有的 FILL_FACTOR * n 个条目复制并重新哈希到新数组中,实际上速度非常快。 - maaartinus
哈希表的工作方式是每插入O(p)个元素就增长一次。当哈希表变得太满时,只需将其大小加倍即可。选择新的质数作为哈希表大小的实际操作非常简单(甚至可以说是微不足道),因为我们并不太关心它是否是一个好的质数,甚至几乎不关心它是否是质数。使用质数作为哈希表大小只是最小化冲突的机会。实际上,调整哈希表本身的大小是O(p),但由于每插入O(p)个元素才执行一次,因此平均时间复杂度为O(1)。 - Dorje
https://dev59.com/cm855IYBdhLWcg3wNhZ1#rs7nnYgBc1ULPQZF1GHn - Will Ness
显示剩余6条评论

12

一个好的技巧是使用部分筛法。例如,接下来的质数是多少,它跟数字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。


“我们知道,在N之后的下一个质数必须是奇数,因此我们可以立即丢弃这些小质数列表中的所有奇数倍数。这些模数允许我们筛选出那些小质数的倍数。” 你确切指的是什么? - venkatKA
@zander - 我增加了更多的解释。 - user85109
有点道理!谢谢 - venkatKA

12

对互质数之间距离做一些实验。


这是为了可视化其他答案而进行的补充。

我从第100,000个(=1,299,709)到第200,000个(=2,750,159)获取了质数。

一些数据:

Maximum interprime distance = 148

Mean interprime distance = 15  

素数间距的频率图:

alt text

素数间距与素数之间的关系:

alt text

仅仅为了看它是否“随机”。不过......


5

没有函数 f(n) 用来计算下一个质数,必须测试一个数字是否为质数。

当查找第 n 个质数时,已知从第 1 个到第 (n-1) 个所有质数非常有用,因为只需要测试这些数字作为因子。

由于这些原因,如果存在预先计算的大质数集合,我不会感到惊讶。对我而言,某些质数需要反复重新计算似乎没有太多意义。


你不需要知道从sqrt(p(n))到p(n-1)的质数。 - Sparr
@Sparr 同意,但您需要这些来计算 p(m),其中 p(m) >= p(n)。编写质数算法时,您会保留所有遇到的质数以供“以后使用”。 - Kirk Broadhurst
有没有可以证明不存在这样的函数?还是这是因为缺乏想象力而得出的结论? - John Shedletsky
没有已知的函数,因此没有可用的函数,因此在所有意义上都没有函数。如果有这样的函数,那么就不需要研究质数了,对吧? - Kirk Broadhurst

3

如果要追求新奇,总有这种方法:

#!/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 的质数。 因此,仅找到下一个质数是不够的。 您还需要进行其他检查。

2
谁会知道,您可以使用正则表达式来搜索质数:/ - Nikita Rybak
1
有趣的是,保罗·惠勒的答案表明微软至少使用了一个4N+1的质数。 - Steve Jessop

3
正如其他人已经指出的那样,目前还没有找到一种方法可以在给定当前质数的情况下找到下一个质数。因此,大多数算法更关注使用快速检查素性的方法,因为您必须检查已知质数和下一个质数之间的n/2个数字。
根据应用程序的不同,您也可以仅使用硬编码的查找表,就像Paul Wheeler所指出的那样。

0
据我所记,它使用当前大小的两倍旁边的质数。它不计算该质数 - 有一个预加载数字表,直到某个大值(不确定,大约在1000万左右)。当达到该数字时,它使用一些天真的算法来获取下一个数字(例如curNum = curNum + 1),并使用以下方法之一进行验证:http://en.wikipedia.org/wiki/Prime_number#Verifying_primality

只存在一对质数,使得curNum和curNum+1都是质数。再努力尝试。 - tchrist
尝试使用next_prime = prime + 2。你可能是对的,而且没有人能证明一旦你足够高,你总是会错。所以去试试吧。 :) - tchrist

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接