C++ 优化该算法

Question

在观看了一些 Terence Tao 的视频后，我想尝试将算法实现到 C++ 代码中，以找到一个数 n 内的所有质数。在我的第一个版本中，我简单地测试了从 2 到 n 的每个整数是否可以被从 2 到 sqrt(n) 中的任何数字整除，并成功让程序在约 52 秒内找到了 1-10,000,000 之间的质数。

为了优化程序并实现我现在所知道的欧拉筛法，我认为任务完成的速度会比 51 秒快得多，但不幸的是，情况并非如此。即使是到 1,000,000，也需要相当长的时间（没有计时）。

#include <iostream>
#include <vector>
using namespace std;

void main()
{
    vector<int> tosieve = {};        
    for (int i = 2; i < 1000001; i++) 
    {                                       
        tosieve.push_back(i);               
    }                                       
        for (int j = 0; j < tosieve.size(); j++)
        {
            for (int k = j + 1; k < tosieve.size(); k++)
            {
                if (tosieve[k] % tosieve[j] == 0)
                {
                    tosieve.erase(tosieve.begin() + k);
                }
            }
        }
    //for (int f = 0; f < tosieve.size(); f++)
    //{
    //  cout << (tosieve[f]) << endl;
    //}
    cout << (tosieve.size()) << endl;
    system("pause");
}

这个程序的运行速度为什么如此缓慢？是因为向量的重复引用吗？即使我可能完全忽略了什么（我是完全的初学者），我认为使用这种可怕而低效的方法在2和1,000,000之间找到质数，应该比我最初的在2到10,000,000之间找到它们的方式要快得多。希望有人能给出清晰的答案 - 希望将来在优化使用大量递归的程序时可以使用所获得的任何知识。

- vinny

---

请修正您的缩进。 - cdonat

从向量的末尾以外的任何位置删除元素都会变慢。 - emlai

埃拉托色尼筛法不测试可除性，它只使用加法。（您无需了解任何算术即可手动执行筛法。） - molbdnilo

优化：不要构建向量并删除非质数；从空向量开始，添加质数。 - molbdnilo

Answer 1

问题在于“erase”会将向量中的每个元素向下移动一位，这意味着它是一个O(n)的操作。

有三个替代选择：

1）只需将删除的元素标记为“空白”（例如将它们变成0）。这意味着未来的传递将不得不经过那些空位置，但这并不太昂贵。

2）创建一个新的向量，并将新值添加到其中。

3）使用std::remove_if：这将使元素下移，但会在单个传递中完成，因此效率更高。如果使用std::remove_if，则必须记住它不会重新调整向量大小。

- Chris Jefferson

Answer 2

大多数 vector 操作，包括 erase()，都具有 O(n) 的线性时间复杂度。

由于您有两个大小为 10^6 的循环和一个大小为 10^6 的 vector，因此您的算法执行高达 10^18 次操作。

对于如此大的 N，立方算法将需要大量时间。
N = 10^6 对于二次算法来说已经足够大了。

请仔细阅读关于埃拉托色尼筛法的内容。两种算法花费相同的时间意味着您做错了第二种算法。

- Yeldar Kurmangaliyev

Answer 3

我看到这里有两个性能问题:

首先，push_back() 会不时重新分配动态内存块。使用 reserve()：

vector<int> tosieve = {};
tosieve.resreve(1000001);       
for (int i = 2; i < 1000001; i++) 
{                                       
    tosieve.push_back(i);               
}

第二个erase()必须移动您尝试删除的元素后面的所有元素。相反，您可以将元素设置为0，并在最后对向量进行操作（未经测试的代码）：

for (auto& x : tosieve) {
    for (auto y = tosieve.begin(); *y < x; ++y) // this check works only in
                                                // the case of an ordered vector
        if (y != 0 && x % y == 0) x = 0;
}
{ // this block will make sure, that sieved will be released afterwards
    auto sieved = vector<int>{};
    for(auto x : tosieve)
        sieved.push_back(x);
    swap(tosieve, sieved);
} // the large memory block is released now, just keep the sieved elements.

考虑使用标准算法而不是手写循环，它们可以帮助您表达意图。在这种情况下，我建议使用 std::transform() 作为筛选器的外部循环，std::any_of() 作为内部循环，std::generate_n() 用于在开始时填充 tosieve，std::copy_if() 用于填充 sieved（未经测试的代码）：

vector<int> tosieve = {};
tosieve.resreve(1000001);
generate_n(back_inserter(tosieve), 1000001, []() -> int {
    static int i = 2; return i++;
});

transform(begin(tosieve), end(tosieve), begin(tosieve), [](int i) -> int {
    return any_of(begin(tosieve), begin(tosieve) + i - 2,
                  [&i](int j) -> bool {
                      return j != 0 && i % j == 0;
                  }) ? 0 : i;
});
swap(tosieve, [&tosieve]() -> vector<int> {
    auto sieved = vector<int>{};
    copy_if(begin(tosieve), end(tosieve), back_inserter(sieved),
            [](int i) -> bool { return i != 0; });
    return sieved;
});

编辑：

还有另一种方法可以完成这个任务：

vector<int> tosieve = {};
tosieve.resreve(1000001);
generate_n(back_inserter(tosieve), 1000001, []() -> int {
    static int i = 2; return i++;
});
swap(tosieve, [&tosieve]() -> vector<int> {
    auto sieved = vector<int>{};
    copy_if(begin(tosieve), end(tosieve), back_inserter(sieved),
            [](int i) -> bool {
                return !any_of(begin(tosieve), begin(tosieve) + i - 2,
                               [&i](int j) -> bool {
                                   return i % j == 0;
                               });
            });
    return sieved;
});

现在我们不再标记那些后续不需要复制的元素，而是直接复制我们想要的元素。这种方法不仅比上面的建议更快，而且表达意图更清晰。

- cdonat

Answer 4

您有一个非常有趣的任务。谢谢！

我很乐意从头开始实现自己版本的解决方案。

我创建了三个独立的函数，都基于埃拉托色尼筛法。这三个版本在复杂度和速度上都不同。

简单注意一下，我最简单（最慢）的版本在短短的0.025秒内找到了所有小于您所需极限的质数10'000'000。

我还测试了所有3个版本，以找到小于2^32 （4'294'967'296）的质数，" simple "版本用47秒解决，" intermediate "版本用30秒解决，" advanced "版本用12秒解决。因此，在仅仅12秒内，它就可以找到小于40亿的所有质数（在2^32下有203'280'221个这样的质数，请参见OEIS sequence）！

为简单起见，我只会详细描述其中的简单版本，这是代码：

template <typename T>
std::vector<T> GenPrimes_SieveOfEratosthenes(size_t end) {
    // https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
    if (end <= 2)
        return {};
    size_t const cnt = end >> 1;
    std::vector<u8> composites((cnt + 7) / 8);
    auto Get = [&](size_t i){ return bool((composites[i / 8] >> (i % 8)) & 1); };
    auto Set = [&](size_t i){ composites[i / 8] |= u8(1) << (i % 8); };
    std::vector<T> primes = {2};
    size_t i = 0;
    for (i = 1; i < cnt; ++i) {
        if (Get(i))
            continue;
        size_t const p = 2 * i + 1, start = (p * p) >> 1;
        primes.push_back(p);
        if (start >= cnt)
            break;
        for (size_t j = start; j < cnt; j += p)
            Set(j);
    }
    for (i = i + 1; i < cnt; ++i)
        if (!Get(i))
            primes.push_back(2 * i + 1);
    return primes;
}

这段代码实现了寻找素数的最简单但速度最快的算法，称为埃拉托斯特尼筛法。为了优化速度和内存，我只搜索奇数。这种奇数优化使我能够存储2倍少的内存并且执行2倍少的步骤，因此精确地提高了速度和内存消耗2倍。

算法很简单，我们分配一个位数组，该数组在位置K处具有位1，如果K是复合数，则具有0位，如果K可能是质数，则具有0位。最后，数组中所有0位表示明确的质数（对于确定的质数）。由于奇数优化，这个位数组仅存储奇数，因此第K位实际上是数字2 * K + 1。

然后从左到右遍历这个位数组，如果我们在位置K处遇到0位，那么它意味着我们找到了一个素数P = 2 * K + 1，现在从位置（P * P）/2开始，我们将每个P位标记为1。这意味着我们标记大于P * P的所有组合数，因为它们可以被P整除。

我们只执行这个程序，直到P * P变得大于或等于我们的限制End（我们正在找到所有小于End的质数）。这个限制保证在达到它之后，数组中的所有零位都表示质数。

代码的第二个版本只对这个简单版本进行了一项优化，即将其全部多核（多线程）。但是，这种优化使代码变得更加庞大和复杂。基本上，它将整个位范围切片到所有核心中，以便它们并行地将位写入内存。

我将仅解释我的第三个高级版本，它是3个版本中最复杂的。它不仅具有多线程优化，还具有所谓的Primorial优化。

什么是Primorial？它是前几个最小质数的乘积，例如我选择primorial 2 * 3 * 5 * 7 = 210。

我们可以看到，任何primorial都可以通过模除该primorial的整数范围将无限范围的整数分成轮子。例如，primorial 210分成范围[0; 210)，[210; 2210)，[2210; 3*210)等。

现在，我们可以通过数学方法证明，在所有的Primorial范围内，我们可以标记相同的数字位置为复合数，确切地说，我们可以将所有2或3或5或7的倍数标记为复合数。

我们可以看到，在210个余数中，有162个余数肯定是复合数，只有48个余数可能是质数。

因此，我们只需要检查整个搜索空间中的48/210=22.8%的质性。这种减少搜索空间的方法使任务速度提高了4倍以上，并且消耗的内存也减少了4倍以上。

我们可以看到，我的第一个简化版本实际上由于仅使用奇数优化而实际上使用了Primorial等于2的优化。是的，如果我们选择Primorial 2而不是Primorial 210，则会得到第一个版本（Simple）算法。

我所有的3个版本都经过了正确性和速度测试。尽管仍然可能存在一些微小的错误。请注意，除非经过彻底测试，否则不建议直接在生产中使用我的代码。

所有3个版本都通过重复使用彼此的答案进行正确性测试。我通过将所有限制（end值）从0到2^18进行全面测试来彻底测试正确性。这需要一些时间。

请参见main()函数以了解如何使用我的函数。 在线试用！ 源代码在此处。由于StackOverflow每个帖子的限制为30K个符号，因此我无法在此处内联源代码，因为它的大小几乎为30K，并且与上面的英文帖子一起超过了30K。因此，我在单独的Github Gist服务器上提供源代码链接。请注意，上面的在线试用！链接也包含完整的源代码，但由于GodBolt运行时间的限制为3秒，因此我将搜索限制减小到了2 ^ 32以下。 Github Gist代码 输出：

10M time 'Simple' 0.024 sec
Time 2^32 'Simple' 46.924 sec, number of primes 203280221
Time 2^32 'Intermediate' 30.999 sec
Time 2^32 'Advanced' 11.359 sec
All checked till 0
All checked till 5000
All checked till 10000
All checked till 15000
All checked till 20000
All checked till 25000

- Arty