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

Posted

技术标签:

【中文标题】给定素数 N,计算下一个素数?【英文标题】:Given Prime Number N, Compute the Next Prime? 【发布时间】:2011-05-27 10:26:27 【问题描述】:

一位同事刚刚告诉我,由于与散列相关的神秘原因,C# 字典集合按质数调整大小。我的直接问题是,“它怎么知道下一个素数是什么?他们是在运行一个巨大的表还是计算?这是一个可怕的非确定性运行时插入导致调整大小”

所以我的问题是,给定 N,它是一个素数,计算下一个素数的最有效方法是什么?

【问题讨论】:

这确实属于 mathoverflow。 也许你的同事不正确,或者它使用了一些预先计算的素数而不是找到下一个素数。 @Kirk:我不同意——这是一个算法问题,而不是数学问题。 @Kirk 这一切都属于理论计算机科学,它非常处于编程和数学的中间。所以老实说,我认为在这两个网站上发布这个问题都没有问题。 @Kirk:这绝对不属于 MathOverflow,它仅适用于研究级问题。我也不同意它需要在math.stackexchange.com 上,但它至少也适合在那里。 【参考方案1】:

大约一年前,我在为libc++ 工作,同时实施 C++11 的无序(散列)容器。我以为我会分享 我在这里的经历。此体验支持marcog's accepted answer 合理“蛮力”的定义。

这意味着即使是简单的蛮力在大多数情况下也足够快 情况下,平均取 O(ln(p)*sqrt(p))。

我开发了size_t next_prime(size_t n) 的几个实现,其中 此功能的规格是:

返回:大于或等于n的最小素数。

next_prime 的每个实现都伴随着一个辅助函数 is_primeis_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

有七种实现。显示第一个的目的 实施是为了证明如果你未能停止测试候选人 对于sqrt(x) 的因素,素数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

此实现是 brute force 实现中最慢的,并且 与实现 1 的唯一区别是它停止测试素数 当因子超过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 &lt; 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 中的内部循环,并创建一个“小 primes",对于处理小于 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++ 的哈希表的客户端可以选择 他们决定的任何素数对他们的情况最有利,并且表现 对于这个应用程序是完全可以接受的。

【讨论】:

由于几乎任何 CPU 架构上的除法指令都会产生余数和商,因此实现 3 的速度比实现 2 快一倍这一事实表明您的编译器没有很好地优化并使用了两个实现 2 中的除法指令。在 GCC 4.7.1 上,实现 2 确实比实现 3 快 4%,因为不需要额外的乘法。如果您切换两个 if 子句,您的编译器可能会做得更好。事实上,return false 的情况更有可能发生,因此仅出于这个原因就值得切换(1 % 加速)。 O(ln(p)*sqrt(p)) 很遥远。 O(sqrt(p)) 是 prime 的试除法测试的复杂度,它们 间隔为 O(ln(p)) 步,但所有这些数字两者之间是复合的,而不是素数。如果我正确阅读this,则平均复杂度为 O( ln(ln(p)) * sqrt(p) / ln(p) )。 @WillNess: 我只是在引用已接受的答案。确切的复杂性是我的答案的一个附带问题,它的存在是为了证明各种实施策略对性能的影响。【参考方案2】:

以防万一有人好奇:

使用反射器,我确定 .Net 使用一个静态类,该类包含一个硬编码列表,其中包含 ~72 个素数,范围高达 7199369,它扫描至少是当前大小两倍的最小素数,以及大于该大小的大小它通过将所有奇数除以潜在数的 sqrt 来计算下一个素数。这个类是不可变的和线程安全的(即不存储更大的素数以供将来使用)。

【讨论】:

好答案。检查每个奇数并不是确定素数的最有效方法,但可能绝大多数字典包含不到 350 万个键。 STLPort 也使用存储列表:really-ugly-long-gitweb-url 我没有提到它首先测试可被 2 整除,尽管它似乎只尝试从奇数开始,因为它在寻找下一个素数时会增加 2。然而,也有一个奇怪的地方,那就是如果你调用 HashHelpers.GetPrime(7199370) 它将遍历从 7199370 到 2147483646 的所有偶数而不找到素数,然后只返回 7199370。有趣,但当然是内部,所以它可能永远不会以这种方式调用。 呸,我错了,有一个偷偷摸摸的按位或 1 将偶数输入值变成下一个更大的奇数。 不过,实际上并不是我的问题的答案。【参考方案3】:

众所周知,gaps between consecutive prime numbers 非常小,第一个超过 100 的间隙出现在素数 370261。这意味着在大多数情况下,即使是简单的蛮力也足够快,需要 O(ln(p) *sqrt(p)) 平均。

对于 p=10,000,这是 O(921) 次操作。请记住,我们将在每次 O(ln(p)) 插入(粗略地说)时执行一次,这完全在大多数现代硬件上以毫秒为单位的大多数问题的约束范围内。

【讨论】:

在增长字典的情况下,我不会称其为“快速”。 同意复杂性并不过分,但每个操作都是相对繁重的余数检查; & 随着 p 的增加,检查本身的复杂性也会增加。 @GregS 查看我的编辑。 @Kirk 当然,意识到这些费用是成为一名经验丰富的程序员的事情之一。 @marcog 除非我还在睡觉,对于 p = 10000,ln(p) = 9.2 和 sqrt(p) = 100,=> O(920)。 @Kirk 不,睡着的是我。修复!【参考方案4】:

一个不错的技巧是使用部分筛子。例如,数字 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

我们知道 N 之后的下一个素数必须是奇数,我们可以立即丢弃这个小素数列表的所有奇数倍。这些模数允许我们筛选出这些小素数的倍数。如果我们使用不超过 200 的小素数,我们可以使用这个方案立即丢弃大多数大于 N 的潜在素数,除了一个小列表。

更明确地说,如果我们正在寻找一个超过 2534536543556 的素数,它不能被 2 整除,所以我们只需要考虑超过该值的奇数。上面的模数表明 2534536543556 与 2 mod 3 一致,因此 2534536543556+1 与 0 mod 3 一致,2534536543556+7、2534536543556+13 等也必须如此。有效地,我们可以筛选出许多数字而无需任何需要测试它们的素数,并且不进行任何试验。

同样的,

mod(2534536543556,7) = 3

告诉我们 2534536543556+4 等于 0 mod 7。当然,这个数字是偶数,所以我们可以忽略它。但是 2534536543556+11 是一个可以被 7 整除的奇数,2534536543556+25 等也是如此。同样,我们可以将这些数字排除在明显的合数之外(因为它们可以被 7 整除),因此不是质数。

仅使用最多 37 个素数的小列表,我们可以排除紧跟在我们的起点 2534536543556 之后的大部分数字,仅排除少数几个:

2534536543573 , 2534536543579 , 2534536543597

在这些数字中,它们是质数吗?

2534536543573 = 1430239 * 1772107
2534536543579 = 99833 * 25387763

我已经努力提供列表中前两个数字的素数分解。看到它们是复合的,但主要因素很大。当然,这是有道理的,因为我们已经确保剩下的数字不会有小的质因数。我们短名单中的第三个(2534536543597)实际上是 N 之外的第一个素数。我描述的筛选方案将倾向于产生素数,或者由通常较大的素数因子组成。因此,在找到下一个素数之前,我们实际上只需要对少数几个数字应用一个显式的素数测试。

类似的方案很快产生超过 N = 1000000000000000000000000000 的下一个素数,即 1000000000000000000000000103。

【讨论】:

“我们知道 N 之后的下一个素数必须是奇数,我们可以立即丢弃这个小素数列表的所有奇数倍数。这些模数允许我们筛选出那些小素数的倍数。” ? @zander - 我添加了更多解释。 这有点道理!谢谢你【参考方案5】:

只是几个素数间距离的实验。


这是对其他答案的可视化的补充。

我得到了从第 100.000 (=1,299,709) 到第 200.000 (=2,750,159) 的素数

一些数据:

Maximum interprime distance = 148

Mean interprime distance = 15  

质距频率图:

质数距离与质数

只是为了看看它是“随机的”。 However ...

【讨论】:

【参考方案6】:

没有函数 f(n) 来计算下一个素数。相反,必须测试一个数字的素性。

在找到第 n 个素数时,已经知道从第 1 个到第 (n-1) 个的所有素数也非常有用,因为这些是唯一需要作为因子进行测试的数字。

由于这些原因,如果有一组预先计算好的大素数,我不会感到惊讶。如果某些素数需要反复重新计算,这对我来说真的没有意义。

【讨论】:

你不需要知道从 sqrt(p(n)) 到 p(n-1) 的素数。 @Sparr 同意,但您需要计算 p(m),其中 p(m) >= p(n)。在编写素数算法时,您会保留所有遇到的素数以供“稍后”使用。 有没有provably没有这样的功能?还是缺乏想象力的证明? 没有已知功能,因此没有可用的功能,因此无论出于何种意图和目的,都没有功能。如果有这样的功能,那么就不需要研究或研究素数了,对吧?【参考方案7】:

正如其他人已经指出的那样,尚未找到在给定当前素数的情况下找到下一个素数的方法。因此,大多数算法更侧重于使用快速检查 primality 的方法,因为您必须检查已知素数和下一个素数之间的 n/2 个数字。

根据应用程序,您也可以像Paul Wheeler 所指出的那样,只对查找表进行硬编码。

【讨论】:

【参考方案8】:

为了纯粹的新颖性,总是有这种方法:

#!/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 形式的素数。所以仅仅找到下一个素数是不够的。您还必须进行另一项检查。

【讨论】:

谁知道你可以使用正则表达式来搜索素数:/ 有趣的是,Paul Wheeler 的回答表明 MS 至少使用了一个 4N+1 素数。【参考方案9】:

据我所知,它使用当前大小的两倍旁边的素数。它不计算那个素数 - 有预加载数字的表达到某个大值(不完全是,大约 10,000,000 左右)。当达到该数字时,它会使用一些简单的算法来获取下一个数字(如 curNum=curNum+1)并使用以下方法对其进行验证:http://en.wikipedia.org/wiki/Prime_number#Verifying_primality

【讨论】:

只有一对素数使得 curNum 和 curNum+1 都是素数。加倍努力。 试试next_prime = prime + 2。你可能是对的,没有人能证明一旦你足够高,你就永远是错的。所以去吧。 :)

以上是关于给定素数 N,计算下一个素数?的主要内容,如果未能解决你的问题,请参考以下文章

素数和

素数和

第四周:循环控制

素数和

中国大学MOOC-C程序设计(浙大翁恺)—— 素数和

0073-简单的素数问题