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

Question

一位同事刚刚告诉我,C#字典集合按照与散列有关的神秘理由按素数调整大小.而我当前的问题是,"它如何知道下一个素数是什么？他们是故事还是一个巨大的表格或动态计算？这是一个可怕的非确定性运行时插入导致调整大小"

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

Answer 1

大约一年前,我正在为libc ++工作,同时为C++ 11实现无序(哈希)容器.我想我会在这里分享我的经历.这种经验支持了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

这种实现是蛮力实现中最慢的,与实现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 < 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.以下实现使用此公式,但使用可爱的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;
}

这实际上意味着算法必须仅检查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 ++的哈希表的客户端可以选择他们认为对他们的情况最有利的任何素数,并且此应用程序的性能是完全可以接受的.

Answer 2

以防万一有人好奇:

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

Answer 3

已知连续素数之间的间隙非常小,第一个间隙超过100时出现素数370261.这意味着即使是一个简单的蛮力在大多数情况下也会足够快,取O(ln(p)*sqrt(p))平均而言.

对于p = 10,000,即O(921)运算.请记住,每次插入O(ln(p))时我们都会执行此操作(粗略地说),这完全符合大多数现代硬件上大多数问题的约束.

Answer 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.

Answer 5

只是几个实验与素数间距离.

这是对其他答案进行可视化的补充.

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

一些数据:

Maximum interprime distance = 148

Mean interprime distance = 15  

Interprime距离频率图:

Interprime距离与素数

只是为了看它"随机".但是......

Answer 6

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

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

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