求素数

正在看《algorithm in C》，看到了Sieve of Eratosthenes的素数筛法，于是在网上查了下，然后自己写了几个。

这里的问题是：给定N，求出小于N的所有素数，N > 0。至于求出前N个素数，利用素数定理，也可以转换为该问题

关于相应的方法，最近在编程随想的博客里看到一篇文章，写的不错。下面是链接：http://blog.csdn.net/program_think/article/details/7032600，自行前往观赏。我这里讲的会比较简略。

方法1： 除法

用一个数组来保存以求得的素数，然后判断下一个数是不是素数，就用数组里的素数去除。

例如刚开始数组里有[2, 3, 5, 7]，要判断9是不是素数，就先用9出2，然后再除3，发现能除尽，因此9不是素数。相应的能发现11是素数，该方法比较简单。

方法2： 筛法

大致过程是这样的，对于数2,3,...N-1. 对每个数给个标记，标记为1表示素数。初始时所有数的标记都为1.

首先知道2是素数，然后将里面2*2, 2*3, 2*4.... 这些数的标记都设为0，也就是把2的倍数都给筛掉了。

然后知道3的标记为1，3是素数，将3的倍数的标记都设为0。

依次下去，最后保留的标记为1的数，就是素数。

筛法的实现的时候，可以用字符数组来保存标记，也可以用比特串来保存标记。

下面是具体的代码。

先看除法的。 前面的代码里面CLOCK相关的宏是用来计时的，宏M(op, arrp, us)是用来测试函数的。这里的代码和下面给出的代码都放到一个文件里就能编译运行了，所以这里看起来很复杂的样子，其实只要看里面的div_prime函数的实现就行了。

#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <sys/time.h>         /* linux里的头文件 */


#define CLOCK_BEGIN(start) (gettimeofday(&(start), NULL))
#define CLOCK_END(end) (gettimeofday(&(end), NULL))
#define CLOCK_TIME(start, end) ((1e6*((end).tv_sec-(start).tv_sec)+(end).tv_usec-(start).tv_usec))

#define M(op, arrp, us) \
    CLOCK_BEGIN(start);\
    (arrp) = (op);\
    CLOCK_END(end);\
    free((arrp));\
    (us) = CLOCK_TIME(start, end);

unsigned *div_prime(unsigned N, unsigned *count)
{
    int i, j, g;
    int isprime = 1;
    double t;
    unsigned total;     /* count of prime number, init is 3 */
    unsigned *prime;

    /* estimate the number of primes */
    /* 分配内存，分配的大小用N*1.17/log(N)来计算，这个公式是根据素数定理来的，具体参见编程随想的那篇文章，这里要注意的是1.17这个数字，如果使用这个数字，则N最好要大于1000，不然会出问题 */
    prime = malloc(sizeof(unsigned) * (int)(N*1.17/log(N)));
    if(prime == NULL)
    {
        fprintf(stderr, "div_prime: Insufficent memory!\n");
        exit(0);
    }
    prime[0] = 2; prime[1] = 3, prime[2] = 5;
    total = 3;          /* found prime's number */
    /* 明显偶数不是素数（2除外），3的倍数也不是素数，2*3的倍数也不是素数，即我们不用判断6n形式的数，只需判断6n+1,6n+2,6n+3,6n+4,6n+5，但这里面6n+2和6n+4是2的倍数，不需判断，6n+3是3的倍数不需判断，因此只需要判断具有6n+1和6n+5的数是不是素数，即只判断数7, 11, 13, 17, 19, 23....，这里两个数的间隔是2和4的交替。 */

     g = 2;     

    for(i = 7; i < N; i+=g)
    {
        g = 6 - g;     /* only check 6i+1 and 6i+5 */
        for(j = 0; prime[j] * prime[j] <= i; j++)
            if(i % prime[j] == 0)
            {
                isprime = 0; /* non-prime */
                break;
            }
        if(isprime == 0)
            isprime = 1;
        else
            prime[total++] = i;  /* isprime is 1, i is prime */
    }
    *count = total;       /* count 保存找到的素数的个数，这个通过参数返回 */
    return prime;
}

这里关于内存分配和只判断6n+1和6n+5形式的数需要注意一下，其他的就比较容易懂了。

下面是筛法，这里采用unsigned char型数组来保存标记

unsigned char *sieve_prime(unsigned N)
{
    unsigned int i, j;
    unsigned char *prime;
    unsigned int of;   /* test overflow */

    prime = malloc(sizeof(unsigned char) * N);
    if(prime == NULL)
    {
        fprintf(stderr, "sieve_prime: Insufficient memory.\n");
        exit(0);
    }
    for(i = 0; i < N; i++)
        prime[i] = 1;
    prime[0] = 0; prime[1] = 0;
    for(i = 2; i < (sqrt(N)); i++) /* 原来这里使用(unsigned)(sqrt(N)), 但后来发现不正确 */
    {
        if(prime[i] == 1)
        {
            /* be caution about i*j overflow */
            for(j = i,of=i*j; i==of/j && of < N; j++,of+=i)
                prime[of] = 0;
        }
    }
    return prime;
}

关于筛法需要注意的是，需要判断溢出，原来没有判断溢出，并且用int型来保存i和j时，当i=46349就溢出了，因为46349*46349 > 2147483647 = INT_MAX.   

另外这里采用of = i*j, 然后根据i == of/j来判断是否溢出。  也许有更好的溢出判断方法，但我不知道。

老实说，我还真不能保证这代码完全正确，虽然我进行了一些测试，比如比较它与div_prime的结果是否相同等，但我还是怕这里会出问题，主要是溢出这个东西有点麻烦。

下面是利用比特串的方法实现的筛法

#define SHIFT 5
#define MASK 0x1F
#define BITPERWORD 32  /* 8*sizeof(int) */

/* set the ith bit to 1 */
void seti(unsigned *a, unsigned i)
{
    /* i>>SHIFT means i/32
     * i & MASK means i%32 */
    a[i>>SHIFT] |= (1 << (i & MASK));
}

/* clear the ith bit to 0 */
void clri(unsigned *a, unsigned i)
{
    a[i>>SHIFT] &= ~(1 << (i & MASK));
}

/* if the i-th bit is 0, return 0, otherwise 1 */
int testi(unsigned *a, unsigned i)
{
    unsigned t;
    t = a[i>>SHIFT] & (1 << (i & MASK));
    return ((t > 0) ?  1 : 0);
}

unsigned *bitsieve_prime(unsigned N)
{
    unsigned int i, j;
    unsigned *prime;
    unsigned int of; /* test overflow */

    /* store N bits need at least 1+N/BITPERWORD int */
    prime = (unsigned *) malloc(sizeof(unsigned) * (1+N/BITPERWORD));
    if(prime == NULL)
    {
        fprintf(stderr, "bitsieve_prime: Insufficient memory.\n");
        exit(0);
    }
    for(i = 2; i < N; i++)
        seti(prime, i);  /* set all to 1 */
    clri(prime, 0);clri(prime, 1); /* set prime[0] and prime[1] to 0 */

    for(i = 2; i < (sqrt(N)); i++)
    {
        if(testi(prime, i))
        {
            /* be caution about i*j overflow */
            for(j = i,of=i*j; i==of/j && of < N; j++,of+=i)
                clri(prime, of);
        }
    }
    return prime;
}

采用bit串后，除了置0和置1的操作要另外实现，其他和unsigned char型的筛法一样。

这里采用的bit串的操作的代码，是在《programming pearls》的第一章中看到的。比较简洁。

我这里采用的函数实现的方法，应该也可以采用宏来实现。

下面是对这三个函数的时间测试函数

void time_test(void)
{
    struct timeval start, end;
    unsigned *div_p, *bit_p;
    unsigned char *sieve_p;
    unsigned div_us, bit_us, sieve_us;
    unsigned count;
    int i;
/* 仅测试1000到3百万之间的数，并间隔为9999，每个数都测就太花时间了 */
    for(i = 1000; i < 3000000; i += 9999)
    {
        M(div_prime(i, &count), div_p, div_us);  /* 这里的M宏就是div_prime代码里的宏了 */
        M(sieve_prime(i), sieve_p, sieve_us);
        M(bitsieve_prime(i), bit_p, bit_us);
        fprintf(stdout, "%d\t%lf\t%lf\t%lf\n", i, div_us/1e3, sieve_us/1e3, bit_us/1e3);
    }
}

M宏并不复杂，就不解释了。求出的div_us,sieve_us, bit_us都是微秒，最后输出的时候，就改成毫秒了。

要注意，这里的i是从1000开始的，，这个数字与div_prime里的数字1.17有关。如果从小于1000的数开始，内存分配上可能会有问题。

为完整性，下面是main函数，就是调用time_test。

int main(int argc, char *argv[])
{
    time_test();

    return 0;
}

下面是用time_test测试时间后，绘出的图，采用gnuplot工具，学一个绘图工具，对于编程性能分析有很大帮助，:-)，我是初学gnuplot工具。

明显还是筛法更快一些，其中采用unsigned char保存素数的筛法要快一些，不过当N变大之后，就和采用bit串的越来越接近了。有可能bit串的筛法慢的原因是我采用函数的方式实现的seti,clri。如果用宏来实现可能会快一些，但我没进行测试。 不过bit串的空间消耗要比unsigned char的小8倍了，所以综合来看选择比特串的筛法要更实惠一些。

额， 图中的div表示除法的div_prime， sieve表示unsigned char的sieve_prime，bit表示bit串的bitsieve_prime.

别看这是个简单的问题，但实际写起代码来还是有困难的，而且很容易犯错。我看《algorithm in C》一书的官网上给出的筛法的代码就是错误的。这里的代码如有错误还望指出，非常感谢。