星海's Blog

老头初学编程
入门排序算法(逐步更新)
数据结构与算法分析数学和文本部份

C素数筛法(更新中)

星海 posted @ 2011年2月11日 09:20 in 数据结构与算法分析 , 3683 阅读
#include 	<stdio.h>
#include	<stdlib.h>
#include	<math.h>
#include	<stdbool.h>
#include	<time.h>

unsigned int prime(unsigned int max)
{
	unsigned int count = 0;	/* 返回count: 素数个数 */
	unsigned int i, j, k;

// 申请内存
	bool *sieve = (bool *) malloc(sizeof(bool) * (max));
	if (sieve == NULL) {
		fprintf(stderr, "\ndynamic memory allocation failed\n");
		exit(EXIT_FAILURE);
	}

	sieve[0] = false;
	sieve[1] = false;
	sieve[2] = true;

//除2之外,所有偶数都是合数
	for (i = 3; i < max; i++)
		if (i % 2)
			sieve[i] = true;
		else
			sieve[i] = false;

	k = (int)sqrt(max);

	for (i = 2; i < k; i++) {
		if (sieve[i])
			for (j = i + i; j < max; j += i) {
				sieve[j] = false;
			}
	}

	for (i = 0; i < max; i++)
		if (sieve[i])
			count++;

	free(sieve);
	sieve = NULL;

	return count;
}

int main(void)
{
	unsigned int count;
	unsigned int N = 100000000;
	clock_t start, end;

	start = clock();
	count = prime(N);

	end = clock();
	printf("[%u]以内素数个数%u 计算用时:%g 秒\n", N, count,
	       (double)(end - start) / (double)CLOCKS_PER_SEC);
}

[100000000]以内素数个数5761455 计算用时:5.86 秒

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

/*
 * ===  FUNCTION  ======================================================================
 *         Name:  prime
 *  Description:  只存奇数的bool型数组,素数筛法。[0]表示自然数3
 * =====================================================================================
 */
unsigned int prime(unsigned int max)
{
    unsigned int count = 1;    /* 返回count: 素数个数,已经包括2,所以count为1 */
    unsigned int limit = max / 2 - 1;    /* max以内的素数用到的数组需要空间 max/2-1 */
    unsigned int i, k;
    unsigned int multiprime;

// 申请内存,       
    bool *sieve = (bool *) malloc(sizeof(bool) * (limit));
    if (sieve == NULL) {
        fprintf(stderr, "\ndynamic memory allocation failed\n");
        exit(EXIT_FAILURE);
    }

    for (i = 0; i < limit; i++)
        sieve[i] = true;

// 2*i+3 为数组下标内容所代表的数值

    k = (unsigned int)sqrt(max);
    for (i = 0; 2 * i + 3 < k; i++) {
        if (sieve[i]) {
            int temp = 2 * i + 3;    /* temp代表数组下标所代表的真实数值 */
            /* 索引i所代表数值的奇数倍索引,3倍时为3*temp */
            for (multiprime = temp * 3; multiprime < max;
                 multiprime += 2 * temp)
                sieve[(multiprime - 3) / 2] = false;
        }

    }

    for (i = 0; i < limit; i++)
        if (sieve[i])
            count++;

    free(sieve);
    sieve = NULL;

    return count;
}

int main(void)
{
    unsigned int count;
    unsigned int N = 100000000;
    clock_t start, end;

    start = clock();
    count = prime(N);

    end = clock();
    printf("[%u]以内素数个数%u 计算用时:%g 秒\n", N, count,
           (double)(end - start) / (double)CLOCKS_PER_SEC);
}

[100000000]以内素数个数5761455 计算用时:2.85 秒

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

unsigned int prime(unsigned int max)
{
	unsigned int count = 0;	/* 返回count: 素数个数 */
	unsigned int i, j, k;

// 申请内存
	bool *sieve = (bool *) malloc(sizeof(bool) * (max));
	if (sieve == NULL) {
		fprintf(stderr, "\ndynamic memory allocation failed\n");
		exit(EXIT_FAILURE);
	}

	sieve[0] = false;
	sieve[1] = false;
	sieve[2] = true;

//除2之外,所有偶数都是合数
	for (i = 3; i < max; i++)
		if (i % 2)
			sieve[i] = true;
		else
			sieve[i] = false;

	k = (int)sqrt(max);

	/*------------------------------------------------------------------
	 *  台湾ACMTino同学的算法:i是素数,则下一个起点是i*i,把后面所有的i*i+2*x*i筛掉
         *  我的实现还不完备
	 *----------------------------------------------------------------*/
	for (i = 2; i < k; i++) {
		if (sieve[i])
			for (j = i * i; j < max; j += 2 * i) {
				sieve[j] = false;
			}
	}

	for (i = 0; i < max; i++)
		if (sieve[i])
			count++;

	free(sieve);
	sieve = NULL;

	return count;
}

int main(void)
{
	unsigned int count;
	unsigned int N = 100000000;
	clock_t start, end;

	start = clock();
	count = prime(N);

	end = clock();
	printf("[%u]以内素数个数%u 计算用时:%g 秒\n", N, count,
	       (double)(end - start) / (double)CLOCKS_PER_SEC);
}

[100000000]以内素数个数5761455 计算用时:3.53 秒

#include     <stdio.h>
#include    <stdlib.h>
#include    <math.h>
#include    <string.h>
#include    <time.h>
#include    <stdbool.h>

/*
 * ===  FUNCTION  ======================================================================
 *  Description:  只计算奇数数组里的奇数是否为素数,[0]为3
 *          prime(n)中的n为n以内的素数,不包括n本身。
 * =====================================================================================
 */
int prime(unsigned int n)
{
    unsigned int i, j, k;
    unsigned int limit = n / 2 - 1;

    unsigned int count = 1;    /* 数组中不包含的2也为素数 */
    bool *sieve = (bool *) malloc(sizeof(bool) * limit);
    if (sieve == NULL) {
        fprintf(stderr, "\ndynamic memory allocation failed\n");
        exit(EXIT_FAILURE);
    }
    for (i = 0; i < limit; i++)
        sieve[i] = true;

    j = (unsigned int)sqrt(n);
    for (i = 0; 2 * i + 3 <= j; i++)
        if (sieve[i]) {
            int temp = 2 * i + 3;    /* 2*i+3为数组下标实际代表的数值 */
            for (k = temp * temp; k <= n; k += 2 * temp)    /* x是素数,则下一个起点是x*x,并把后面所有的x*x+2*x*i筛掉 */
                sieve[(k - 3) / 2] = false;
        }
    for (i = 0; i < limit; i++)
        if (sieve[i])
            count++;

    free(sieve);
    sieve = NULL;

    return count;

}

int main(int argc, char *argv[])
{
    unsigned int count;
    unsigned int N = 100000001;
    clock_t start, end;

    start = clock();
    count = prime(N);

    end = clock();
    printf("[%u]以内素数个数%u 计算用时:%g 秒\n", N, count,
           (double)(end - start) / (double)CLOCKS_PER_SEC);
    return 0;
}

[100000000]以内素数个数5761455 计算用时:2.56 秒

/*-----------------------------------------------------------------------------
 *  可参考论文COMPUTING π(x): THE MEISSEL-LEHMER METHOD
 *  以下函数为百度C语言贴吧 5230所做 (我完全不懂撒 -__-)
 *-----------------------------------------------------------------------------*/
#include	<stdio.h>
#include	<string.h>
#include	<stdlib.h>
#include	<time.h>
#include	<math.h>

int *primarr, *v;
int q = 1, p = 1;

//π(n)
int pi(int n, int primarr[], int len)
{
	int i = 0, mark = 0;
	for (i = len - 1; i > 0; i--) {
		if (primarr[i] < n) {
			mark = 1;
			break;
		}
	}
	if (mark)
		return i + 1;
	return 0;
}

//Φ(x,a)
int phi(int x, int a, int m)
{
	if (a == m)
		return (x / q) * p + v[x % q];
	if (x < primarr[a - 1])
		return 1;
	return phi(x, a - 1, m) - phi(x / primarr[a - 1], a - 1, m);
}

int prime(int n)
{
	char *mark;
	int mark_len;
	int count = 0;
	int i, j, m = 7;
	int sum = 0, s = 0;
	int len, len2, len3;

	mark_len = (n < 10000) ? 10002 : ((int)exp(2.0 / 3 * log(n)) + 1);

	//筛选n^(2/3)或n内的素数
	mark = (char *)malloc(sizeof(char) * mark_len);
	memset(mark, 0, sizeof(char) * mark_len);
	for (i = 2; i < (int)sqrt(mark_len); i++) {
		if (mark[i])
			continue;
		for (j = i + i; j < mark_len; j += i)
			mark[j] = 1;
	}
	mark[0] = mark[1] = 1;

	//统计素数数目
	for (i = 0; i < mark_len; i++)
		if (!mark[i])
			count++;

	//保存素数
	primarr = (int *)malloc(sizeof(int) * count);
	j = 0;
	for (i = 0; i < mark_len; i++)
		if (!mark[i])
			primarr[j++] = i;

	if (n < 10000)
		return pi(n, primarr, count);

	//n^(1/3)内的素数数目
	len = pi((int)exp(1.0 / 3 * log(n)), primarr, count);
	//n^(1/2)内的素数数目
	len2 = pi((int)sqrt(n), primarr, count);
	//n^(2/3)内的素数数目
	len3 = pi(mark_len - 1, primarr, count);

	//乘积个数
	j = mark_len - 2;
	for (i = (int)exp(1.0 / 3 * log(n)); i <= (int)sqrt(n); i++) {
		if (!mark[i]) {
			while (i * j > n) {
				if (!mark[j])
					s++;
				j--;
			}
			sum += s;
		}
	}
	free(mark);
	sum = (len2 - len) * len3 - sum;
	sum += (len * (len - 1) - len2 * (len2 - 1)) / 2;

	//欧拉函数
	if (m > len)
		m = len;
	for (i = 0; i < m; i++) {
		q *= primarr[i];
		p *= primarr[i] - 1;
	}
	v = (int *)malloc(sizeof(int) * q);
	for (i = 0; i < q; i++)
		v[i] = i;
	for (i = 0; i < m; i++)
		for (j = q - 1; j >= 0; j--)
			v[j] -= v[j / primarr[i]];

	sum = phi(n, len, m) - sum + len - 1;
	free(primarr);
	free(v);
	return sum;
}

int main()
{
	int n;
	int count;
	clock_t start, end;
	scanf("%d", &n);

	start = clock();
	count = prime(n);
	end = clock() - start;
	printf
	    ("[%d]内的素数个数为%d,计算用时:%d毫秒\n",
	     n, count, (int)end / 1000);
	getchar();
	return 0;
}

 


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter