c - Spoj PRIME1 使用埃拉托色尼筛(在 C 中)

Question

我正在尝试解决问题 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*/

                        }
                }
        }
}

Answer 1

你的问题是

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'作为标志，使用 1和 0正如 DMR 的意图，这将得到更好的理解。

重新更新:这个

/*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 中较大的位置开始内循环和 p 的最小倍数不少于 m .后者的代码有效，因此将其包装在

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