当前位置： 文江博客话题详情

使用埃拉托斯特尼筛的 Spoj PRIME1（C 中）

发布于 2025-01-07 11:14:26 字数 7101 浏览 0 评论 0原文

我正在尝试使用埃拉托斯特尼分段筛来解决问题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;
}

更新：感谢丹尼尔的回复。我实现了您的一些建议，我的代码现在产生了正确的输出：-

/*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*/

                        }
                }
        }
}
原文

I'm trying to solve the problem PRIME1 using segmented sieve of Eratosthenes. My program works correctly with the normal sieve that is up to NEW_MAX. But there is some problem with cases n > NEW_MAX, where segmented sieving comes into play. In such cases it merely prints all the numbers. Here is the link to the code with relevant test cases: 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;
}

Update: Thanks fr your response Daniel. I implemented some of ur suggestions, my code now produces correct output :-

/*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;
}

but my sieve function below further taking into account ur other suggestions produces incorrect output:-

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*/

                        }
                }
        }
}

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（2）

离笑几人歌 2025-01-14 11:14:26

你的问题是，

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; } 在循环体中（更好的是，直接计算 p 不小于 m 的第一个倍数，完全避免这种情况）。

并且不使用 '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) 的素数倍数”并不意味着你不需要消除较小素数的倍数，这意味着你的初始版本没有t，这导致该范围内的几乎所有数字都没有被消除。您应该以 p = 2 开始外层循环 - 或者对 2 的倍数进行单独的传递，并以 3 开始该循环，将 p 增加 2，然后开始内部循环位于p*p中的较大者且p的最小倍数不小于m。后者的代码有效，因此将其包装在 an 中

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[] 您可以获得更好的缓存局部性。通过每个标志使用一位，您可以获得更好的缓存局部性，但这会使读取/设置/清除标志仍然更慢。什么能产生最佳性能取决于筛子的结构和尺寸。

Your problem is

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

You only ever eliminate the multiples of 2 as such if the lower end of the range is 4 or less, and you only eliminate multiples of primes larger than sqrt(m). Remove the i >= m part from the loop condition and replace it with an if (i < m) { continue; } in the loop body (better, calculate the first multiple of p not less than m directly and avoid that condition completely).

And instead of using 't' and 'f' as flags, use 1 and 0 as DMR intended, that will be understood better.

Re update: This

/*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)

will hurt you if m == 2. If m == 2, you must start with i = 4.

Concerning

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

it seems you misunderstood me, "and you only eliminate multiples of primes larger than sqrt(m)" doesn't mean you needn't eliminate multiples of smaller primes, it means that your initial version didn't, which resulted in almost all numbers in the range not eliminated. You should start the outer loop with p = 2 -or have a separate pass for the multiples of 2 and start that loop with 3, incrementing p by 2, and start the inner loop at the larger of p*p and the smallest multiple of p not less than m. Your code for the latter works, so wrapping it in an

if ((i = p*p) < m) {
    /* set i to smallest multiple of p which is >= m */
}

would work (you can make it a bit faster avoiding a branch and using only one division, but the difference will be minuscule).

Finally, your choice of what you represent by 0 and 1 is uncanonical, this is C, so 0 evaluates to false in conditions and everything else to true, so the canonical replacement would have been 't' -> 1 and 'f' -> 0 and in contexts like this, where the array entries are flags, one would check

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

also there's no need to change the array types from char[] to int[], char is an integer type too, so 0 and 1 are perfectly fine char values. The choice of type for the arrays has performance implications. On the one hand, int-sized loads and stores are typically faster than byte-sized, so that would favour int flags[] or even long int flags[] for word-sized reads and writes. On the other hand, with the smaller char flags[] you get better cache locality. You would get even better cache locality with using a single bit per flag, but that would make reading/setting/clearing flags still slower. What yields the best performance depends on the architecture and size of the sieves.

回复 原文
江南烟雨〆相思醉 2025-01-14 11:14:26

丹尼尔·费舍尔似乎已经澄清了一些主要的麻烦点。

如果您有兴趣查看此素数问题的一些更简洁的代码/解释，请查看：

Daniel Fischer seems to have clarified some of the main trouble points.

If you're interested in seeing some more concise code/explanation for this primes problem, check out:

http://www.swageroo.com/wordpress/spoj-problem-2-prime-generator-prime1/

回复 原文
~没有更多了~

关于作者

邮友

暂无简介

文章
评论
25 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

西西弗的石头怪

文章 0 评论 0

5397313

文章 0 评论 0

烟沫凡尘

文章 0 评论 0

一个破名字

文章 0 评论 0

萌︼了一个春

文章 0 评论 0

当爱已成负担

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文