埃拉托斯特尼筛法算法

Question

我正在阅读《使用C++编程：原理和实践》一书，第四章有一个练习题：

我需要编写一个程序，使用埃拉托色尼筛法计算1到100之间的质数。

这是我想出来的程序：

#include <vector>
#include <iostream>

using namespace std;

//finds prime numbers using Sieve of Eratosthenes algorithm
vector<int> calc_primes(const int max);

int main()
{
    const int max = 100;

    vector<int> primes = calc_primes(max);

    for(int i = 0; i < primes.size(); i++)
    {
        if(primes[i] != 0)
            cout<<primes[i]<<endl;
    }

    return 0;
}

vector<int> calc_primes(const int max)
{
    vector<int> primes;

    for(int i = 2; i < max; i++)
    {
        primes.push_back(i);
    }

    for(int i = 0; i < primes.size(); i++)
    {
        if(!(primes[i] % 2) && primes[i] != 2)
             primes[i] = 0;
        else if(!(primes[i] % 3) && primes[i] != 3)
             primes[i]= 0;
        else if(!(primes[i] % 5) && primes[i] != 5)
             primes[i]= 0;
        else if(!(primes[i] % 7) && primes[i] != 7)
             primes[i]= 0;
    }   

    return primes;
}

虽然我在这本书中还不是很了解C++，但我会尽力帮忙。

现在的问题是，只有当max大于500时，才能将所有值打印到控制台上，如果max > 500，则不会打印所有内容。

我做错了什么吗？

P.S.：任何建设性的批评都将不胜感激。

- RaouL

---

这不是埃拉托色尼筛法。你似乎只筛选可被2、3、5和7整除的值，但其他质数的倍数呢？提示：你根本不需要用到取模运算。参考维基百科上的“埃拉托色尼筛法”页面。 - UncleBens

请问您能否提供一个开箱即用的编译示例？ - anon

@UncleBens-我这样做是因为维基百科上的例子从表格中删除了所有2、3、5、7的倍数，但我想我理解错了。 - RaouL

参见Melissa O'Neill关于筛法算法的真正含义的论述，以及Lambda the Ultimate讨论(http://lambda-the-ultimate.org/node/3127)。 - Charles Stewart

@AndreyT：这篇文章的标题是“埃拉托斯特尼筛法算法”。上面的代码并没有实现筛法。 - jon hanson

Answer 1

把筛子看作一个集合。
按顺序遍历集合。对于筛子中的每个值，删除所有可被它整除的数字。

#include <set>
#include <algorithm>
#include <iterator>
#include <iostream>


typedef std::set<int>   Sieve;

int main()
{
    static int const max = 100;

    Sieve   sieve;

    for(int loop=2;loop < max;++loop)
    {
        sieve.insert(loop);
    }


    // A set is ordered.
    // So going from beginning to end will give all the values in order.
    for(Sieve::iterator loop = sieve.begin();loop != sieve.end();++loop)
    {
        // prime is the next item in the set
        // It has not been deleted so it must be prime.
        int             prime   = *loop;

        // deleter will iterate over all the items from
        // here to the end of the sieve and remove any
        // that are divisable be this prime.
        Sieve::iterator deleter = loop;
        ++deleter;

        while(deleter != sieve.end())
        {
            if (((*deleter) % prime) == 0)
            {
                // If it is exactly divasable then it is not a prime
                // So delete it from the sieve. Note the use of post
                // increment here. This increments deleter but returns
                // the old value to be used in the erase method.
                sieve.erase(deleter++);
            }
            else
            {
                // Otherwise just increment the deleter.
                ++deleter;
            }
        }
    }

    // This copies all the values left in the sieve to the output.
    // i.e. It prints all the primes.
    std::copy(sieve.begin(),sieve.end(),std::ostream_iterator<int>(std::cout,"\n"));

}

- Martin York

Answer 2

我不知道为什么你没有得到所有的输出，因为看起来你应该得到所有的内容。你缺少哪些输出？

筛法实现有误。类似这样的代码：

vector<int> sieve;
vector<int> primes;

for (int i = 1; i < max + 1; ++i)
   sieve.push_back(i);   // you'll learn more efficient ways to handle this later
sieve[0]=0;
for (int i = 2; i < max + 1; ++i) {   // there are lots of brace styles, this is mine
   if (sieve[i-1] != 0) {
      primes.push_back(sieve[i-1]);
      for (int j = 2 * sieve[i-1]; j < max + 1; j += sieve[i-1]) {
          sieve[j-1] = 0;
      }
   }
}

我会实现筛法。(上面的代码是我凭空写的，不保证能够工作或编译。我认为它没有超出第四章的范围。)

像往常一样返回primes，并打印出全部内容。

- David Thornley

---

假设我将max初始化为10000，那么它只会打印8000到9000之间的整数。 - RaouL

你是怎么进行检查的？你是打印到窗口吗？有些输出可能会滚动到顶部。 - David Thornley

我正在使用cout向控制台打印内容，奇怪的是它从8000开始打印。 - RaouL

使用新的筛法实现后，一切都正常打印 :=) - RaouL

Answer 3

来自算法和数据结构：

void runEratosthenesSieve(int upperBound) {
      int upperBoundSquareRoot = (int)sqrt((double)upperBound);
      bool *isComposite = new bool[upperBound + 1];
      memset(isComposite, 0, sizeof(bool) * (upperBound + 1));
      for (int m = 2; m <= upperBoundSquareRoot; m++) {
            if (!isComposite[m]) {
                  cout << m << " ";
                  for (int k = m * m; k <= upperBound; k += m)
                        isComposite[k] = true;
            }
      }
      for (int m = upperBoundSquareRoot; m <= upperBound; m++)
            if (!isComposite[m])
                  cout << m << " ";
      delete [] isComposite;
}

- Eric J.

---

该问题特别要求使用 vector<>。 - David Thornley

哪里？他的实现恰巧使用了它，但我没有看到他的问题中有关于vector<>的特定问题。 - Eric J.

Answer 4

有趣的是，似乎没有人回答你关于输出问题的问题。我在代码中没有看到任何应该根据max值影响输出的东西。

说实话，在我的Mac上，我得到了所有的输出。当然是错误的，因为算法不正确，但我确实得到了所有的输出。如果你继续有输出问题，你没有提到你运行的平台，这可能会有用。

以下是您的代码版本，最小程度地修改以遵循实际的Sieve算法。

#include <vector>
#include <iostream>

using namespace std;

//finds prime numbers using Sieve of Eratosthenes algorithm
vector<int> calc_primes(const int max);

int main()
{
    const int max = 100;

    vector<int> primes = calc_primes(max);

    for(int i = 0; i < primes.size(); i++)
    {
        if(primes[i] != 0)
            cout<<primes[i]<<endl;
    }

    return 0;
}

vector<int> calc_primes(const int max)
{
    vector<int> primes;

    // fill vector with candidates
    for(int i = 2; i < max; i++)
    {
        primes.push_back(i);
    }

    // for each value in the vector...
    for(int i = 0; i < primes.size(); i++)
    {
        //get the value
        int v = primes[i];

        if (v!=0) {
            //remove all multiples of the value
            int x = i+v;
            while(x < primes.size()) {
                primes[x]=0;
                x = x+v;
            }
        }
    }
    return primes;
}

- Mark Bessey

Answer 5

在下面的代码片段中，数字在插入vector之前进行了过滤。除数来自于向量。

我还通过引用传递了向量。这意味着巨大的向量不会从函数复制到调用者中。(大块内存复制需要很长时间)

vector<unsigned int> primes;

void calc_primes(vector<unsigned int>& primes, const unsigned int MAX)
{
    // If MAX is less than 2, return an empty vector
    // because 2 is the first prime and can't be placed in the vector.
    if (MAX < 2)
    {
         return;
    }

    // 2 is the initial and unusual prime, so enter it without calculations.
    primes.push_back(2);
    for (unsigned int number = 3; number < MAX; number += 2)
    {
        bool is_prime = true;
        for (unsigned int index = 0; index < primes.size(); ++index)
        {
            if ((number % primes[k]) == 0)
            {
                is_prime = false;
                break;
            }
        }

        if (is_prime)
        {
            primes.push_back(number);
        }
    }    
}

这不是最高效的算法，但它遵循了筛法算法。

- Thomas Matthews

---

第一条评论说你返回2作为第一个质数，但实际上并没有。为什么你要在MAX之前停止？ - Eric J.

实际上，这个注释说我正在返回，因为#2是第一个质数，而MAX小于2。 - Thomas Matthews

Answer 6

这里是一个简洁明了的实现，使用bool类型：

#include <iostream>
#include <cmath>

void find_primes(bool[], unsigned int);
void print_primes(bool [], unsigned int);

//=========================================================================
int main() 
{
    const unsigned int max = 100;
    bool sieve[max];

    find_primes(sieve, max);

    print_primes(sieve, max);
}
//=========================================================================

/*
    Function: find_primes()
    Use: find_primes(bool_array, size_of_array);

    It marks all the prime numbers till the
    number: size_of_array, in the form of the
    indexes of the array with value: true.

    It implemenets the Sieve of Eratosthenes,
    consisted of:

    a loop through the first "sqrt(size_of_array)"
    numbers starting from the first prime (2).

    a loop through all the indexes < size_of_array,
    marking the ones satisfying the relation i^2 + n * i
    as false, i.e. composite numbers, where i - known prime 
    number starting from 2.
*/
void find_primes(bool sieve[], unsigned int size)
{
    // by definition 0 and 1 are not prime numbers
    sieve[0] = false;
    sieve[1] = false;

    // all numbers <= max are potential candidates for primes
    for (unsigned int i = 2; i <= size; ++i)
    {
        sieve[i] = true;
    }

    // loop through the first prime numbers < sqrt(max) (suggested by the algorithm)
    unsigned int first_prime = 2;
    for (unsigned int i = first_prime; i <= std::sqrt(double(size)); ++i)
    {
        // find multiples of primes till < max
        if (sieve[i] = true)
        {
            // mark as composite: i^2 + n * i 
            for (unsigned int j = i * i; j <= size; j += i)
            {
                sieve[j] = false;
            }
        }
    }
}

/*
    Function: print_primes()
    Use: print_primes(bool_array, size_of_array);

    It prints all the prime numbers, 
    i.e. the indexes with value: true.
*/
void print_primes(bool sieve[], unsigned int size)
{
    // all the indexes of the array marked as true are primes
    for (unsigned int i = 0; i <= size; ++i)
    {
        if (sieve[i] == true) 
        {
            std::cout << i <<" ";
        }
    }
}

涵盖数组情况。一个std::vector实现将包括一些小的更改，例如将函数减少到一个参数，通过引用传递向量，并且循环将使用向量的size()成员函数而不是减少的参数。

- Ziezi

Answer 7

以下是我的版本，基本上使用了一个布尔位向量，然后通过奇数和快速加法找到要设置为 false 的倍数。最后构造一个向量并将素数值返回给客户端。

std::vector<int>  getSieveOfEratosthenes ( int max )
{
  std::vector<bool> primes(max, true);
  int sz = primes.size();

  for ( int i = 3; i < sz ; i+=2 )
    if ( primes[i] ) 
      for ( int j = i * i; j < sz; j+=i)
        primes[j] = false;

  std::vector<int> ret;
  ret.reserve(primes.size());
  ret.push_back(2);

  for ( int i = 3; i < sz; i+=2 )
    if ( primes[i] )
      ret.push_back(i);

  return ret;
}

- bjackfly

---

我在max=1000000处遇到了一个错误，它出现在第一个循环中。我用以下代码替换了原来的代码：int m = (int)sqrt((double)sz); for (int i = 3; i < m; i += 2)，现在它可以快速运行了。 - MrHIDEn

Answer 8

我现在正在跟随同一本书。我已经想出了算法的以下实现。

#include<iostream>
#include<string>
#include<vector>
#include<algorithm>
#include<cmath>
using namespace std;
inline void keep_window_open() { char ch; cin>>ch; }

int main ()
{
    int max_no = 100;

    vector <int> numbers (max_no - 1);
    iota(numbers.begin(), numbers.end(), 2);

    for (unsigned int ind = 0; ind < numbers.size(); ++ind)
    {
        for (unsigned int index = ind+1; index < numbers.size(); ++index)
        {
            if (numbers[index] % numbers[ind] == 0)
            {
                numbers.erase(numbers.begin() + index);
            }
        }
    }
    cout << "The primes are\n";
    for (int primes: numbers)
    {
        cout << primes << '\n';
    }
}

- scientific_explorer

Answer 9

这里有一个经典的方法来实现这个功能，

int main()
{
    int max = 500;
    vector<int> array(max); // vector of max numbers, initialized to default value 0

    for (int i = 2; i < array.size(); ++ i) // loop for rang of numbers from 2 to max
    {
        // initialize j as a composite number; increment in consecutive composite numbers
        for (int j = i * i; j < array.size(); j +=i)
            array[j] = 1;  // assign j to array[index] with value 1
    }

    for (int i = 2; i < array.size(); ++ i) // loop for rang of numbers from 2 to max
        if (array[i] == 0)  // array[index] with value 0 is a prime number
        cout << i << '\n';  // get array[index] with value 0

    return 0;
}

- Khaled Eid

Answer 10

这是我实现的更高效的埃拉托斯特尼筛法算法版本。

#include <iostream>
#include <cmath>
#include <set>

using namespace std;

void sieve(int n){
    set<int> primes;
    primes.insert(2);
    for(int i=3; i<=n ; i+=2){
        primes.insert(i);
    }       

    int p=*primes.begin();
    cout<<p<<"\n";
    primes.erase(p);

    int maxRoot = sqrt(*(primes.rbegin()));

    while(primes.size()>0){
        if(p>maxRoot){
            while(primes.size()>0){
                p=*primes.begin();
                cout<<p<<"\n";
                primes.erase(p);        
            }
            break;
        }

        int i=p*p;  
        int temp = (*(primes.rbegin()));
        while(i<=temp){
            primes.erase(i);
            i+=p;
            i+=p;
        }
        p=*primes.begin();
        cout<<p<<"\n";
        primes.erase(p);
    }
}

int main(){
    int n;
    n = 1000000;
    sieve(n);                
    return 0;
}

- agasheesh