质数算法

Question

有人能告诉我如何在C语言中实现埃拉托色尼筛法算法吗？我需要生成质数，但我的算法速度很慢。

我的代码：

#include <stdio.h>

int prime(long int i)
{
    long int j;
    int state = 1;
    for(j=2;j<i;j++)
    {
        if((i%j)==0){state=0;break;}
    }
    return state;
}

int main()
{
    int t;
    long int m,n,i;
    scanf("%d", &t);
    while(t--) {
        scanf("%d %d", &m,&n);
        for(i=m;i<=n;i++)
        {
            if(i==1){
                //do nothing for 1
            } else{
                if(prime(i))printf("%d\n",i);
            }
        }
    }
    return 0;
}

t是测试用例的数量，m和n是打印质数的范围。

- jaykumarark

---

是的，直接的筛法是慢的。如果你把你目前的代码发出来，也许有人可以帮助你改进它。 - aschepler

5秒谷歌搜索：http://www.dreamincode.net/code/snippet3315.htm这段话是一个网址，指向DreamInCode网站上的一个代码片段。 - Marc B

你想要作为结果的最大整数是多少？ - Benoit

看一下你在帖子中提到的维基页面底部，你会找到一些链接，这些链接是Sieve of Eratosthenes算法在C语言中的实现。 - DReJ

尤其是在无限排除每个偶数后，再无限排除每个3的倍数等情况下，程序可能会非常缓慢，虽然它仍在运行，但我不明白为什么。 - Benoit

你的代码根本不是筛法，而只是一个子优化的试除法循环，其复杂度为n^2。将测试限制在被测试数的平方根范围内将使复杂度降至n^1.5。仅通过质数进行测试而不像你所做的那样测试所有数字，将进一步降低复杂度约一个对数因子。但是埃拉托色尼筛法的复杂度为n log log n。 - Will Ness

Answer 1

您需要创建一个布尔数组，其大小与您要查找的最大质数一样大。在开始时，它完全被初始化为真。

这个数组的第i个单元格如果是素数，则为true；否则为false。

从 i=2 开始迭代：它是素数，则将任何索引为2的倍数的单元格设置为false。继续下一个素数(i=3)，然后再次进行相同的操作。继续前进到下一个素数（现在是i=5:i=4不是素数：在处理 i=2 时已将 array[4] 设置为false），并重复该过程。

- peoro

---

当你测试第i个数字时，为什么不按照i步长前进呢？（counter += i） - BlackBear

@BlackBear：你不行。如果你在 i=2，然后跳到 4，那么你就跳过了 3... 无论如何，你可以找到一些类似于你提出的优化方法，以快速“移动到下一个质数”，但我认为它们不能改善你的算法复杂度（甚至你的算法也不能）。 - peoro

@BlackBear 是的，你说得对！通过迭代增加索引（通过 counter += whatever）正是埃拉托斯特尼筛法的本质。这就是如何“将索引倍数为2的任何单元格设置为false”必须完成的方式。如果每个索引都被测试以发现它是否是倍数，则它将成为一个试除法筛法，其时间复杂度要差得多。 - Will Ness

所以这个答案对于算法的关键点保持沉默，作者在评论中的建议将使其变得完全错误。 - Will Ness

Answer 2

在我看来，你的算法缓慢是因为你计算了不必要的数字。 尝试使用以下代码：

int isPrime(int number){

    if(number < 2) return 0;
    if(number == 2) return 1;
    if(number % 2 == 0) return 0;
    for(int i=3; (i*i)<=number; i+=2){
        if(number % i == 0 ) return 0;
    }
    return 1;

}

- Toomtarm Kung

Answer 3

这里是使用埃拉托斯特尼筛法的非常简单的代码，适用于所有正整数int。

int is_prime(int n){
  int p;
  for(p = 2; p < n; p++){
    if(n % p ==0 && p != n)
      return 0;    
  }
  return 1;
}

- user7658994

---

不错，我给你的答案点了赞。修改建议：你对质数的定义是错误的，而且 #include <stdio.h> 是多余的。 - MarianD

不，这不是埃拉托斯特尼筛法。这是一种非常低效的试除法素性测试实现。 - Will Ness

Answer 4

马克·B的链接展示了一个简单正确的算法，由NSLogan编写。我对其进行了轻微修改以显示一些尺寸/速度权衡。出于兴趣，我想分享一下。

首先，是代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <time.h>

#define USE_BITS

#ifdef USE_BITS
#define alloc_prime char *prime = calloc(i/8+1,sizeof(*prime));
#define set_not_prime(x) prime[x/8]|= (1<<(x%8))
#define is_prime(x) (!(prime[x/8]&(1<<(x%8))))
#endif

#ifdef USE_CHAR
#define alloc_prime char *prime = calloc(i+1,sizeof(*prime));
#define set_not_prime(x) prime[x] = 1
#define is_prime(x) (prime[x] == 0)
#endif

#ifdef USE_SIZE_TYPE
#define alloc_prime size_t *prime = calloc(i+1,sizeof(*prime));
#define set_not_prime(x) prime[x] = 1
#define is_prime(x) (prime[x] == 0)
#endif

 
int main(){
    int i;
    printf("Find primes up to: ");
    scanf("%i",&i);
 
    clock_t start, stop;
    double t = 0.0;
       
    assert((start = clock())!=-1);
       
    //create prime list
    alloc_prime;
    int c1, c2, c3;
    if(!prime){
    printf("Can't allocate %zu bytes.\n",i*sizeof(*prime));
    exit(1);
    }
       
    //set 0 and 1 as not prime
    set_not_prime(0);
    set_not_prime(1);

    //find primes then eliminate their multiples (0 = prime, 1 = composite)
    for(c2 = 2;c2 <= (int)sqrt(i)+1;c2++){
      if(is_prime(c2)){
        c1=c2;
        for(c3 = 2*c1;c3 <= i+1; c3 += c1){
          set_not_prime(c3);
        }
      }
    }
       
    stop = clock();
    t = (double) (stop-start)/CLOCKS_PER_SEC;
       
    //print primes
    for(c1 = 0; c1 < i+1; c1++){
        if(is_prime(c1))printf("%i\n",c1);
    }
    printf("Run time: %f\n", t); //print time to find primes

    return 0;
}

由于使用了sqrt，编译时请使用：gcc prime.c -lm

我比较了三种不同的存储布尔变量的方法，这些方法是peoro建议的。其中一种方法实际上使用了位，第二种方法占用整个字节，最后一种方法占用整个机器字。一个关于哪种方法最快的天真猜测可能是机器字方法，因为每个标志/布尔值都更“自然”地被处理。在我的机器上计算到100,000,000的质数所需的时间如下：

位：3.8秒 字符：5.8秒 M-字：10.8秒

有趣的是，即使所有丑陋的位移也只能表示一个比特的布尔值，总体上仍然更快。我的猜测是缓存和引用局部性似乎超过了额外的约3条指令。此外，您最终会使用n分之1的n位机器字方法的内存！

- Rooke

Answer 5

尽管这是一篇非常古老的帖子，但下面是我使用"埃氏筛法"算法生成素数的尝试。

#include <stdio.h>

#define NUM 8000        /* Prime Numbers in the Range.  'Range + 2' actually. */

int main()
{
  int a[NUM] = {0};         /* Array which is going to hold all the numbers */
  int i , j;

  /* initializing array from 2 to given number + 2, assuming the first prime number is 2 */
  for(i = 2,j=0; i < NUM+2, j<NUM; i++, j++) 
  {
    a[j] =i;
  }


  for(i = 0; i < NUM; i++ ) 
  {
    int num = a[i];

    /* If number is not 0 then only update the later array index. */
    if(num != 0) 
    {
      for (j = i+1; j < NUM; j++) 
      {
        if( (a[j]%num == 0) ) 
        {
            a[j]=0;
        }
      }
    }
  }


  for(i = 0; i < NUM; i++) 
  {
    /* Print all the non Zero *Prime numbers* */
    if(a[i] != 0) 
    {
      printf("%d \n", a[i]);
    }
  }

}

希望这能帮助到某个人。

- pankaj

---

不，这不是埃拉托斯特尼筛法。它是试除法筛法。其复杂度比埃氏筛法更差（随着NUM的增大，运行速度变得越来越慢）。 - Will Ness

Answer 6

第一步是认识到将其分割成任意大小的块是微不足道的；而且您不需要为每个数字使用数组（或位域）。例如，如果您只关心100000到110000之间的数字，则可以将所有可被3整除的数字标记为“非质数”，方法是"for( index = 3 * (100000+3-1)/3; index < 110000; index += 3) {"。以下是更灵活的示例：

    for( value = 2; value < sqrt( block_end_value-1); value++ ) {
        for( index = value  * (block_state_value+value -1)/value; index < block_end_value; index += value ) {
            mark_not_prime(index - block_state_value);
        }
    }

第二步是认识到你不需要关心每个数字（上面的for( value = 2; value < sqrt( block_end_value-1); value++)是低效的）。例如，如果你已经标记了可被2整除的数字为“非质数”，那么就没有必要关心能否被4、6或8整除；如果你已经标记了可被3整除的数字为“非质数”，那么就没有必要关心能否被6、9或12整除；......基本上，你只想测试一个数字是否能被另一个质数整除。这意味着要在100000到110000范围内找到质数，你首先要在0到sqrt(110000)范围内找到质数；如果你想在0到sqrt(110000)范围内找到质数，你需要在0到sqrt(sqrt(110000))范围内找到质数；......

第三步是意识到可以通过复制重复模式来加速。您可以创建一个2位模式（表示“可被2整除”），并将这两位复制到每个位置。同样，您可以创建一个6位模式（表示“可被2或3整除”），并将这6位复制到每个位置。同样，您可以创建一个39916800位模式（表示“可被2、3、4、5、6、7、8、9、10和11整除”），并将该39916800位模式复制到每个位置。当然，您也可以预先生成模式并将其存储在程序代码中。

第四步是意识到2的倍数太过简单，不用存储它们，这样可以减半所有表格和模式（以及任何存储/预生成的模式）的内存消耗。

通过结合以上第三步和第四步，表示“可被2、3、4、5、6、7、8、9、10和11整除”的预生成模式将花费19958400比特，约为2.38 MiB。独立使用此模式的第一部分也可用于查找从1到第一个大于11的质数（在本例中为1至13的数字）。

第五步是认识到如果你已经有了一个模式，你可以使用它来找到“N =下一个“未标记”的质数”，将现有的模式复制N次，然后将N的倍数标记为“非质数”；最终得到一个更大的模式。例如；如果你有一个表示“可被2、3、4、5、6、7、8、9、10和11整除”的模式，你会跳过12（因为根据现有模式它不是质数）；将该模式复制13次，然后将可被13整除的数字标记为“非质数”，最终得到一个表示“可被2、3、4、5、6、7、8、9、10、11、12和13整除”的模式。

第六步是意识到一旦您拥有足够大的模式覆盖您所需的范围，您可以填补缺失的除数而无需复制。例如，如果您只想在1到6227020800的范围内得到质数，则可以采用代表“可被2、3、4、5、6、7、8、9、10、11、12和13整除”的模式，然后将在14到6227020800范围内可被质数整除的数字标记为“非质数”。

通过结合上述所有内容，如果你想在1000000000到11000000000的范围内找到质数，则需要首先在1到sqrt(11000000000)的范围内查找质数；为此，您需要复制预生成的模式并扩展它，直到获得足够大的表格来表示1到sqrt(11000000000)范围内的质数；然后复制该扩展模式并填充缺失的数字以生成表示1到sqrt(11000000000)范围内的质数的表格；接着创建1000000000到11000000000范围内的质数表，并通过从“扩展预生成的模式”中复制数据来初始化内存，然后使用1到sqrt(11000000000)范围内的质数表来查找尚未合并到“扩展预生成的模式”中的质数，以查找仍需标记为“非质数”的数字，这些数字是生成1000000000到11000000000范围内数字表所需的。

- Brendan

Answer 7

Toomtarm Kung的回答很好，非常感谢。

然而，我仍然遇到了一个问题：在32位上，当i>46340时，(i*i)可能会溢出；在64位上，当i>3037000499时，也会产生错误的结果，即使用32位整数时，2147483647将无法被识别为质数。

为了避免整数溢出，可以将(i * i)<=number替换为i <= number / i

现在，为了避免双重除法，代码可以写成如下形式：

int isPrime(int number){
    if(number < 2) return 0;
    if(number == 2) return 1;
    if(number % 2 == 0) return 0;
    int j;
    for(int i=3; ((j=number/i) >= i); i+=2){
        if(number - (j * i) == 0 ) return 0;
    }
    return 1;
}

- joheirba