各种高性能算法 求素数 太牛 挑战算法

发个某牛根据intfree改写的程序，计算10000000000（100亿）以素数个数，只要18秒多，我的CPU是AMD1.8G，要在INTEL的3.0G机子上跑肯定在10秒内。。。

C/C++ code
#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <time.h>

__int64 maxn = 10000000L;
const long maxlen = 100000;
long m = 5;
const long dots = 80;

BYTE *mask;
long primes[maxlen+2], mprimes[maxlen+2], ans[maxlen+2];
long bits[256];
long len, r;

void solve(long a, long b, long c, long &x, long &y)
{
if(a == 0)
{
x = 0;
y = c/b;
}
else
{
solve(b%a, a, c, y, x);
x -= b/a*y;
}
}

// primes数组从下标1开始存放素数，例如primes[1]为2
void init()
{
long p, k;

k = 2;
len = 0;

// 求待求范围的平凡根范围内的素数
do
{
++len;
primes[len] = k;

do
{
++k;
p = 1;
while((primes[p]*primes[p] < k) && (k%primes[p] != 0))
{
++p;
}
}while(k%primes[p] == 0);
}while(1.0*k*k <= maxn);

bits[0] = 0;
for(k=1; k<=255; ++k)
{
bits[k] = bits[k>>1] + (k & 1);
}
}

void makeprimes()
{
long i, k, n, size, kmax, min;
long last, p, x, dot;

// 以2*3*5*7*11*13为基，循环完后n为30030
n = 1;
for(i=1; i<=m; ++i)
n *= primes[i];

for(i=m+1; i<=len; ++i)
{
p = primes[i];
solve(p, -n, 1, x, mprimes[i]);
mprimes[i] = (mprimes[i]%p+p)%p;
}

r = len;
size = long(maxn/n+7)>>3;
last = 0;
dot = 0;

memset(ans, 0, sizeof(ans));
mask = new BYTE[size];
assert(mask);

for(i=1; i<n; ++i)
{
k = 1;
while(k < m && (i%primes[k] != 0))
++k;
if(i%primes[k] == 0) // 最小值即为2，3，5，7，11，13的倍数的段可以直接排除
continue;

p = i-last;
last = i;
kmax = (long)((maxn-i)/n);
for(k=m+1; k<=len; ++k)
ans[k] = (ans[k] + mprimes[k]*p)%primes[k];

// 重置mask数组为0
memset(mask, 0, sizeof(BYTE)*size);
for(k=m+1; k<=len; ++k)
{
p = primes[k];
if(ans[k] == 0)
min = p-1;
else
min = ans[k]-1;

// if((i>(maxn-maxn/n*n) && (min > maxn/n-2)) || (i<=(maxn-maxn/n*n) && (min > maxn/n-1)))
// {
// printf("min>size*8: i=%d, k = %d, min=%d, size=%d,prime=%d\n", i, k, min,size, primes[k]);
// continue;
// }

// 下面的汇编代码用来给所有的合数设置标志位
__asm
{
mov eax, mask
mov ebx, kmax
mov ecx, min
mov edx, p
$001:
bts [eax], ecx // 置mask数组的第min位为1
add ecx, edx // min位置指向下一个合数
cmp ecx, ebx // 是否到了kmax位置
jl $001 // 不到则转$001位置
}
}

// 去掉本段内所有的合数的个数
for(k=0; k<size; ++k)
{
assert(k <= size);
kmax -= bits[mask[k]];
}
// 加上本段内的素数个数
r += kmax;
if(i/n > dot/dots)
{
++dot;
putchar('.');
}
}

delete []mask;

putchar('\n');
}

int main()
{
long start, end, t;

printf("Input a range: ");
scanf("%I64d", &maxn);

start = clock();
init();
end = clock();
printf("Init time=%lf\n", (double)(end-start)/1000);

start = clock();
makeprimes();
end = clock();

t = end-start;
printf("sum=%ld, cacl time=%lf\n", r, (double)t/1000);

return 0;
}

结果：
Input a range: 10000000000
Init time=0.032000

sum=455052511, cacl time=18.578000
Press any key to continue

#include "iostream"
#include "windows.h"
using namespace std;

unsigned int prime(unsigned int lmax)
{
bool* prime_num = new bool[lmax];
memset(prime_num, true, sizeof(bool) * lmax);

for (unsigned int i = 2; i < lmax; ++i)
{
if (prime_num[i])
for (unsigned int j = 2; i * j < lmax; ++j)
prime_num[i * j] = false;
}

unsigned int pnum = 0;
for (unsigned int i = 2; i < lmax; ++i)
{
if (prime_num[i])
++pnum;
}

delete[] prime_num;
return pnum;
}

int main()
{
int a = 1000000000;
DWORD s = GetTickCount();
unsigned int pnum = prime(1000000000);
float use_time = ((float)(GetTickCount() - s)) / 1000.f;

cout<<"1000000000 以内有素数个数 : "<<pnum<<endl;
cout<<"用时 : "<<use_time<<"秒"<<endl;

return 0;
}

1000000000 以内有素数个数 : 50847534
用时 : 353.641秒
请按任意键继续. . .

前段日子我也在做素数，下面是我的一个总结，看到那个2秒搞定的，我还真想去看看

当数字小于1000000时,可以用简单的判断

C/C++ code
int isprime(int n)
{
int i;
for(i=2;i<=sqrt(n);i++)
if (n%i==0)
return 0;
return 1;
}

但当数值在1000000到100000000时
介绍一种方法
prime={2,3,5,7,11,13,17,……}为事先做好的素数表
如果需要判断的数最大为100000000，则prime的最大元素为大于10000的最小素数即可

C/C++ code

//素数表生成法：
int t[10010]={0},prime[3000];//3000可能有点大，具体多少运行了，就知道了
int len;
void getprime(void)
{
int i,j;
t[0]=t[1]=1;          //t[i]=1,表示该数不是素数，被筛除
for(i=2;i  <=sqrt(10010);i++)
{
if(!t[i])
{
for(j=i+i;j <10010;j+=i)
t[j]=1;
}
}
len=0;
for(i=2;i  <10010;i++)
{
if(!t[i])
prime[len++]=i;
}
}

int isprime(long n)
{
int i;
while(prime[i]*prime[i]  <=n)
{
if(n%prime[i]==0)
return 0;
i++;
}
return 1;
}

如果数值大于100000000时
可以用Miller-Rabbin素数测试法，判断是否为素数

C/C++ code
int Miller_Rabbin(long long n)
{
long long i,s,a;
s=10; //s的值可以根据需要变大
// randomize();
for(i=0;i <s;i++)
{
a=long long(rand()%(n-1)+2); //自动生成受限
if(modular_exp(a,n-1,n)>1)
return 0;
}
return 1;
}
long long modular_exp(long long a,long long b,long long c)//求a^b%c该函数受限
{
if(a==0)
return 0;
if(b==0)
return 1;
if(b==1)
return a%c;
return (a*modular_exp(a,b-1,c))%c;
}

刷选法：

#include <stdio.h>
#include <dos.h>

typedef struct node
{
unsigned long i;
struct node *next;
} TPm;

void main()
{
struct time tm;
unsigned long i,n;
bool q;
TPm *p,*t,*t1,*t2;
printf( "Input N: ");
scanf( "%lu ",&n);
gettime(&tm);
printf( "The program started at %02d:%02d:%02d.%02d\n\n ",tm.ti_hour, tm.ti_min, tm.ti_sec, tm.ti_hund);
p=new TPm;
t=new TPm;
p-> i=2;
p-> next=t;
t-> i=3;
t-> next=NULL;
for(i=5;i <=n;i+=2)
{
t=p;
q=true;
while(t)
{
if(!(i%t-> i))
{
q=false;
break;
}
t1=t;
t=t-> next;
}
if(q)
{
t2=new TPm;
t2-> i=i;
t2-> next=NULL;
t1-> next=t2;
printf( "%lu\t ",i);
}
}
t=p;
i=0;
while(t)
{
i++;
p=t;
// printf( "%lu\t ",t-> i);
t=t-> next;
delete p;
}
printf( "\n\nTotal: %lu\n ",i);
gettime(&tm);
printf( "The program stopped at %02d:%02d:%02d.%02d\n ",tm.ti_hour, tm.ti_min, tm.ti_sec, tm.ti_hund);

}