logo标签列表

使用埃拉托斯特尼筛法解决Spoj PRIME1问题（C语言版）。

cprimessieve-of-eratosthenes
3
我正在尝试使用埃拉托色尼筛法解决问题PRIME1。我的程序在使用普通筛法时，即在NEW_MAX范围内，能够正确运行。但是在一些n > NEW_MAX的情况下，需要使用分段筛法。在这种情况下，它只会打印出所有数字。以下是带有相关测试用例的代码链接：http://ideone.com/8H5lK#view_edit_box。
/* segmented sieve */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define MAX_LIMIT 1000000000  //10^9
#define NEW_MAX  31623 /*SQUARE ROOT OF 1000000000*/
#define MAX_WIDTH 100000   //10^5
char flags[NEW_MAX+100];  /*TO PREVENT SEGMENTATION FAULT*/ 

void initialise(char flagarr[], long int n) //initialise all elements to true from 1 to n
{
    long int i;
    for (i = 1; i <= n; i++)
        flagarr[i] = 't';
}

void sieve(unsigned long long m, unsigned long long n, char seg_flags[])
{
    unsigned long long p, i, limit;
    if (m == 1)
        seg_flags[1] = 'f';
    /*Seperate inner loop for p=2 so that evens are not iterated*/
    for (i = 4; i >= m && i <= n; i += 2)
    {
        seg_flags[i-m+1] = 'f';
    }
    if (seg_flags == flags)
        limit = NEW_MAX;
    else
        limit = sqrt(n);
    for (p = 3; p <= limit+1; p += 2)  //initial p+=2 bcoz we need not check even
    {
        if (flags[p] == 't')
        {
            for (i = p*p; i >= m && i <= n; i += p)  //start from p square since under it would have been cut
            seg_flags[i-m+1] = 'f';          /*p*p,  p*p+p,  p*p + 2p   are not primes*/
        }
    }
}

void print_sieve(unsigned long long m,unsigned long long n,char flagarr[])
{
    unsigned long long i;
    if (flags == flagarr)    //print non-segented sieve 
    {
        for (i = m; i <= n; i++)
            if (flagarr[i] == 't')
                printf("%llu\n", i);
    }
    else
    {
        //print segmented
        for (i = m; i <= n; i++)
            if (flagarr[i-m+1] == 't')
                printf("%llu\n", i);
    }
}

int main()
{
    unsigned long long m, n;
    int t;
    char seg_flags[MAX_WIDTH+100];
    /*setting of flags for prime nos. by sieve of erasthromas upto NEW_MAX*/
    initialise(flags, NEW_MAX);
    flags[1] = 'f'; /*1 is not prime*/
    sieve(1, NEW_MAX, flags);
    /*end of initial sieving*/
    scanf("%d", &t);
    while (t--)
    {
        scanf("%llu %llu", &m, &n);
        if (n <= NEW_MAX)
            print_sieve(m, n, flags); //NO SEGMENTED SIEVING NECESSARY
        else if (m > NEW_MAX)
        {
            initialise(seg_flags, n-m+1);  //segmented sieving necessary
            sieve(m, n, seg_flags);
            print_sieve(m, n, seg_flags);
        }
        else if (m <= NEW_MAX && n > NEW_MAX)  //PARTIAL SEGMENTED SIEVING NECESSARY
        {
            print_sieve(m, NEW_MAX, flags);
            /*now lower bound for seg sieving is new_max+1*/
            initialise(seg_flags, n-NEW_MAX);
            sieve(NEW_MAX+1, n, seg_flags);
            print_sieve(NEW_MAX+1, n, seg_flags);
        }
        putchar('\n');
    }
    system("pause");
    return 0;
}

更新：感谢您的回复，Daniel。我实施了您的一些建议，现在我的代码可以产生正确的输出 :-

/*segmented sieve*/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#define MAX_LIMIT 1000000000  /*10^9*/
#define NEW_MAX  31623 /*SQUARE ROOT OF 1000000000*/
#define MAX_WIDTH 100000   /*10^5*/
int flags[NEW_MAX+1];  /*TO PREVENT SEGMENTATION FAULT goblal so initialised to 0,true*/    

void initialise(int flagarr[],long int n) 
/*initialise all elements to true from 1 to n*/
{
    long int i;
    for(i=3;i<=n;i+=2)
        flagarr[i]=0;
}

void sieve(unsigned long m,unsigned long n,int seg_flags[])
{

    unsigned long p,i,limit;  

    /*Seperate inner loop for p=2 so that evens are not iterated*/
    if(m%2==0)
        i=m;
    else
        i=m+1;
    /*i is now next even > m*/
    for(;i<=n;i+=2)
    {

        seg_flags[i-m+1]=1;

    }
    if(seg_flags==flags)
        limit=NEW_MAX;
    else
        limit=sqrt(n);
    for(p=3;p<=limit+1;p+=2)  /*initial p+=2 bcoz we need not check even*/
    {
        if(flags[p]==0)
        {


            for(i=p*p; i<=n ;i+=p)  
            /*start from p square since under it would have been cut*/
            {
                if(i<m)
                    continue;
                seg_flags[i-m+1]=1; 
                     /*p*p,  p*p+p,  p*p + 2p   are not primes*/

            }
        }
    }
}

void print_sieve(unsigned long m,unsigned long n,int flagarr[])
{
    unsigned long i;
    if(m<3)
    {printf("2\n");m=3;}
    if(m%2==0)
        i=m+1;
    else
        i=m;
if(flags==flagarr)    /*print non-segented sieve*/  
{
    for(;i<=n;i+=2)
        if(flagarr[i]==0)
                printf("%lu\n",i);
}
else {
 //print segmented

    for(;i<=n;i+=2)
        if(flagarr[i-m+1]==0)
                printf("%lu\n",i);
}}

int main()
{
    unsigned long m,n;
    int t;
    int seg_flags[MAX_WIDTH+100];
    /*setting of flags for prime nos. by sieve of erasthromas upto NEW_MAX*/
    sieve(1,NEW_MAX,flags);
    /*end of initial sieving*/
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lu %lu",&m,&n);
        if(n<=NEW_MAX)
            print_sieve(m,n,flags); 
            /*NO SEGMENTED SIEVING NECESSARY*/
        else if(m>NEW_MAX)
        {
            initialise(seg_flags,n-m+1);  
            /*segmented sieving necessary*/
            sieve(m,n,seg_flags);
            print_sieve(m,n,seg_flags);
        }
        else if(m<=NEW_MAX && n>NEW_MAX) 
         /*PARTIAL SEGMENTED SIEVING NECESSARY*/
        {
            print_sieve(m,NEW_MAX,flags);
            /*now lower bound for seg sieving is new_max+1*/
            initialise(seg_flags,n-NEW_MAX);
            sieve(NEW_MAX+1,n,seg_flags);
            print_sieve(NEW_MAX+1,n,seg_flags);
        }
        putchar('\n');
    }
    system("pause");
    return 0;
}

但是，我下面的筛选函数在考虑了您的其他建议后仍会产生错误的输出：-
void sieve(unsigned long m,unsigned long n,int seg_flags[])
{

        unsigned long p,i,limit;  
        p=sqrt(m);
        while(flags[p]!=0)
                p++;
        /*we thus get the starting prime, the first prime>sqrt(m)*/

        if(seg_flags==flags)
                limit=NEW_MAX;
        else
                limit=sqrt(n);
        for(;p<=limit+1;p++)  /*initial p+=2 bcoz we need not check even*/
        {
                if(flags[p]==0)
                {
                        if(m%p==0) /*to find first multiple of p>=m*/
                                i=m;
                        else
                                i=(m/p+1)*p;

                        for(; i<=n ;i+=p)  
                        /*start from p square since under it would have been cut*/
                        {

                                seg_flags[i-m+1]=1;     
                                         /*p*p,  p*p+p,  p*p + 2p   are not primes*/

                        }
                }
        }
}
2个回答
2

你的问题是

for (i = 4; i >= m && i <= n; i += 2)
for (i = p*p; i >= m && i <= n; i += p)

如果范围的下限为4或更低，则只消除2的倍数，如果要消除的素数大于sqrt(m)。从循环条件中删除i >= m部分，并在循环体中替换为if (i < m) { continue; }（更好的方法是直接计算第一个不小于m的p的倍数，并完全避免该条件）。

而且，不要使用't'和'f'作为标志，而是使用DMR预期的1和0，这样更容易理解。

关于更新：

/*Seperate inner loop for p=2 so that evens are not iterated*/
if(m%2==0)
    i=m;
else
    i=m+1;
/*i is now next even > m*/
for(;i<=n;i+=2)

如果 m == 2，会对您造成影响。如果 m == 2，您必须从 i = 4 开始。

关于

unsigned long p,i,limit;  
p=sqrt(m);
while(flags[p]!=0)
    p++;
/* snip */
for(;p<=limit+1;p++)

看起来你误解了我的意思，“而且你只消除大于sqrt(m)的质数的倍数”并不意味着你不需要消除更小的质数的倍数，它意味着你的初始版本没有这样做，这导致几乎所有范围内的数字都没有被消除。你应该从p = 2开始外循环，或者有一个单独的循环用于2的倍数，并从3开始该循环，每次将p增加2，并在p*p和不小于m的最小p倍数中选择较大的值作为内循环的起点。你对后者的代码是正确的，所以将其包装在一个

标记中即可。
if ((i = p*p) < m) {
    /* set i to smallest multiple of p which is >= m */
}

这段代码可以运行（你可以通过避免使用分支和仅使用一个除法来加快速度，但差别微不足道）。

最后，你选择用0和1表示的方式并不是规范的。在C语言中，0在条件语句中被视为false，其他所有值都被视为true。因此，规范的替换应该是't' -> 1和'f' -> 0。在像这样的上下文中，数组条目是标志位时，应该检查

if (flags[p])   // instead of: if (flags[p] != 0)
if (!flags[p])  // instead of: if (flags[p] == 0)

此外，没有必要将数组类型从char[]更改为int[]，因为char也是一种整数类型，因此0和1是完全可以作为char值的。选择数组类型会影响性能。一方面，int大小的加载和存储通常比字节大小的快，所以会倾向于使用int flags[]甚至long int flags[]进行字大小的读写。另一方面，较小的char flags[]可以获得更好的缓存局部性。使用每个标志位一个比特位可以获得更好的缓存局部性，但这会使读取/设置/清除标志仍然变慢。最佳性能取决于架构和筛子的大小。
1
非常感谢您的巨大帮助，我终于在 SPOJ 上获得了 AC。 - ksb
0
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接