各种高性能算法 求素数 太牛 挑战算法
发个某牛根据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);
}
浙公网安备 33010602011771号