1000億以内素数カウントアルゴリズム
転載先:
とても良い文章で、効率はかなり高くて、残念ながら注釈が少なくて、少し怒っているように見えます
1000億以内素数カウントアルゴリズム
/******************************************************************
copyright (C) 2007 Huang Yuanbing
version 1.1, 2007 PrimeNumber
mailto: [email protected] (remove the last digit 1 for "laji" mail)
free use for non-commercial purposes
******************************************************************
main ideal come from a paper by M.DELEGLISE ADN J.LAGARIAS
"COMPUTING PI(x): THE MEISSEL, LEHMER, LAGARIAS, MILLER, ODLYZKO METHOD"
a = PI(y); (y >>= x^(1/3) and y <= x^(1/2))
PI(x) = phi(x, a) + a - 1 - P2Xa(x,a);
phi(x, a) = S0 + S
= S0 + S1 + S3 + S2
= S0 + S1 + S3 + U + V
= S0 + S1 + S3 + U + V1 + V2
= S0 + S1 + S3 + U + V1 + W1 + W2 + W3 + W4 + W5
it need O(n^2/3) space and MAXN > 10 and MAXN <= 1e11
**********************************************************************/
＃i nclude <time.h>
＃i nclude <stdio.h>
＃i nclude <math.h>
＃i nclude <memory.h>
＃i nclude <assert.h>
#define COMP 4
#define MASKN(n) (1 << ((n >> 1) & 7))
#define min(a, b) (a) < (b) ? (a) : (b)
#define MULTI_THREAD
#define THRED_NUMS 4
//#define PRINT_DEBUG
typedef unsigned int uint;
#ifdef _WIN32
typedef __int64 int64;
#else
typedef long long int64;
#endif
int *Prime;
int *PI;
int64 MAXN = (int64)1e11;
int pt[130000]; // MAXN < 1e11, size = np
int *phiTable;
int *minfactor;
int x23, x12, x13, x14, y;
int64 x;
int ST = 7, Q = 1, phiQ = 1;
//Q = 2*3*5*7*11*13 = 30030
//phiQ = 1*2*4*6*10*12 = 5760
int phi(int x, int a)
{
if (a == ST)
return (x / Q) * phiQ + phiTable[x % Q];
else if (a == 0)
return x;
else if (x < Prime[a])
return 1;
else
return phi(x, a - 1) - phi(x / Prime[a], a - 1);
}
#ifdef MULTI_THREAD
struct ThreadInfo{
int pindex;
int theadnums;
int64 pnums;
}Threadparam[2 * THRED_NUMS + 2];
#ifdef _WIN32
＃i nclude <windows.h>
DWORD WINAPI WIN32ThreadFun(LPVOID pinfo)
#else
＃i nclude <pthread.h>
void* POSIXThreadFun(void *pinfo)
#endif
{
ThreadInfo *pThreadInfo = (ThreadInfo *) (pinfo);
int addstep = pThreadInfo->theadnums * 2;
for (int i = pThreadInfo->pindex; pt[i] != 0; i += addstep){
if (pt[i] > 0)
pThreadInfo->pnums -= phi(pt[i], pt[i + 1]);
else
pThreadInfo->pnums += phi(-pt[i], pt[i + 1]);
}
// printf("ThreadID = %4d, phi(%d) = %d
", GetCurrentThreadId(), pThreadInfo->pindex, pThreadInfo->pnums);
// printf("ThreadID = %4d, phi(%d) = %d
", pthread_self(), pThreadInfo->pindex, pThreadInfo->pnums);
return 0;
}
int64 multiThread(int theadnums)
{
int i;
int64 pnums = 0;
assert(theadnums < 2 * THRED_NUMS);
for (i = 0; i < theadnums; i++){
Threadparam[i].pindex = 2 * i + 1;
Threadparam[i].theadnums = theadnums;
Threadparam[i].pnums = 0;
}
#ifdef _WIN32
HANDLE tHand[THRED_NUMS * 2];
DWORD threadID[THRED_NUMS * 2];
for (i = 0; i < theadnums; i++){
tHand[i] = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE)WIN32ThreadFun, (LPVOID)(&Threadparam[i]), 0, &threadID[i]);
if (tHand[i] == NULL)
printf("create Win32 thread error
");
}
WaitForMultipleObjects(theadnums, tHand, true, INFINITE);
for (i = 0; i < theadnums; i++){
pnums += Threadparam[i].pnums;
CloseHandle(tHand[i]);
}
#else
pthread_t tid[THRED_NUMS];
for (i = 0; i < theadnums; i++){
int error = pthread_create(&tid[i], NULL, POSIXThreadFun, &Threadparam[i]);
if ( error != 0 )
printf("Create pthread error: %d
", error);
}
for (i = 0; i < theadnums; i++){
pthread_join(tid[i], NULL);
pnums += Threadparam[i].pnums;
}
#endif
return pnums;
}
#endif
int freememory(int alloc)
{
int memsize = (x / y) + 100;
if (alloc == 0){
int psize = (int) (memsize / log(memsize));
//printf("psize = %d memeszie = %d
", psize, (int)(psize * 1.2) );
PI = new int [memsize + 100];
Prime = new int[(int)(psize * 1.2) + 100];
assert(PI && Prime);
}else{
delete phiTable;
delete minfactor;
delete Prime;
delete PI;
}
return alloc;
}
void init_phiTable( )
{
clock_t start = clock();
int p, i, j;
if (x < 1e10)
ST = 6;
if (ST > PI[y])
ST = PI[y];
for (i = 1; i <= ST; ++i){
Q *= Prime[i];
phiQ *= Prime[i] - 1;
}
phiTable = new int[Q + 10];
for (i = 0; i < Q; ++i)
phiTable[i] = i;
for (i = 1; i <= ST; ++i){
p = Prime[i];
for (j = Q - 1; j >= 0; --j)
phiTable[j] -= phiTable[j / p];
}
printf(" Q = %d, PhiQ = %d
", Q, phiQ);
printf(" init_phiTable time = %ld ms
", clock() - start);
}
void init_minFactor( )
{
clock_t start = clock();
int i, j, maxy = y + 10;
int sqrty = (int)sqrt(maxy) + 1;
minfactor = new int[maxy + 10];
for (i = 0; i < maxy; i++)
minfactor[i] = i;
for (i = 1; Prime[i] <= maxy; i++){
for (j = Prime[i]; j <= maxy; j += Prime[i]){
if (minfactor[j] == -j || minfactor[j] == j)
minfactor[j] = -Prime[i];
else
minfactor[j] = -minfactor[j];
}
}
for (i = 1; Prime[i] <= sqrty; i++){
int powp = Prime[i] * Prime[i];
for (j = powp; j <= maxy; j += powp)
minfactor[j] = 0;
}
printf(" init_minFactor time = %ld ms
", clock() - start);
}
int sieve( )
{
clock_t start = clock();
int primes = 1;
int maxp = (x / y) + 10;
Prime[1] = 2;
unsigned char *bitMask = (unsigned char *) PI;
memset(bitMask, 0, (maxp + 64) >> COMP);
for (int p = 3; p < maxp; p += 2){
if ( !(bitMask[p >> COMP] & MASKN(p)) ){
Prime[++primes] = p;
if (p > maxp / p)
continue;
for (int j = p * p; j < maxp; j += p << 1)
bitMask[j >> COMP] |= MASKN(j);
}
}
Prime[0] = primes;
Prime[primes] = 0;
printf("pi(%d) = %d
", maxp, primes);
printf(" sieve time = %ld ms
", clock() - start);
return primes;
}
void init_x( )
{
x = (int64)MAXN;
x23 = (int)(pow(x, 2.0 / 3) + 0.01);
x12 = (int)(pow(x, 1.0 / 2) + 0.01);
x13 = (int)(pow(x, 1.0 / 3) + 0.01);
x14 = (int)(pow(x, 1.0 / 4) + 0.01);
assert((int64)x12 * x12 <= x);
assert((int64)x13 * x13 * x13 <= x);
assert((int64)x14 * x14 * x14 * x14 <= x);
y = x13;
printf(" x14 = %d, x13 = %d, x12 = %d x23 = %d, y = %d
", x14, x13, x12, x23, y);
freememory(0);
}
int64 cal_S0( )
{
int64 S0 = x;
for (int j = 2; j <= y; j++){
if (minfactor[j] > 0)
S0 += x / j;
else if (minfactor[j] < 0)
S0 -= x / j;
}
printf("S0 = %I64d
", S0);
return S0;
}
// so bad performance in this function
int64 cal_S3( )
{
clock_t start = clock();
int i, p, a;
int np = 1;
int64 S3 = 0;
for (i = 1; i <= PI[x14]; i++){
p = Prime[i];
a = PI[p] - 1;
#ifdef PRINT_DEBUG
assert(p <= x14 && a >= 0);
#endif
for (int m = y / p + 1; m <= y; m++){
int xx = x / (int64)(m * p);
#ifndef MULTI_THREAD
if (minfactor[m] > p)
S3 -= phi(xx, a);
else if (minfactor[m] < -p)
S3 += phi(xx, a);
#else
if (minfactor[m] > p){
pt[np++] = xx;
pt[np++] = a;
}
else if (minfactor[m] < -p){
pt[np++] = -xx;
pt[np++] = a;
}
#endif
}
}
#ifdef MULTI_THREAD
printf("np = %d
", np);
if (np < 100)
S3 = multiThread(1);
else
S3 = multiThread(THRED_NUMS);
#endif
printf("S3 = %I64d
", S3);
printf(" cal S3 time = %ld ms
", clock() - start);
return S3;
}
void init_PI( )
{
clock_t start = clock();
int Px = 1;
PI[0] = PI[1] = 0;
for (int i = 1; Prime[i]; i++, Px++){
for (int j = Prime[i]; j <= Prime[i + 1]; j++)
PI[j] = Px;
}
//printf(" Px = %d, primes = %d
", Px, Prime[0]);
printf(" PI[%d] = %d
", Px, PI[Px - 1]);
printf(" init_PI time = %ld ms
", clock() - start);
}
int64 cal_P2xa( )
{
int64 phi2 = (int64)PI[y] * (PI[y] - 1) / 2 - (int64)PI[x12] * (PI[x12] - 1) / 2;
for (int i = PI[y] + 1; i <= PI[x12] + 0; i++){
int p = Prime[i];
#ifdef PRINT_DEBUG
assert(p > y && p <= x12);
#endif
phi2 += PI[x / p];
}
printf("P2xa(%I64d, %d) = %I64d
", x, PI[y], phi2);
return phi2;
}
int64 cal_S1( )
{
int64 temp = PI[y] - PI[x13];
int64 S1 = temp * (temp - 1) / 2;
printf("S1 = %I64d
", S1);
return S1;
}
int64 cal_U( )
{
int64 p, U = 0;
int sqrt_xdivy = (int)sqrt(x / y);
for (int i = PI[sqrt_xdivy] + 1; i <= PI[x13]; i++){
p = Prime[i];
U += PI[y] - PI[x / (p * p)];
}
printf("U = %I64d
", U);
return U;
}
int64 cal_V1( )
{
int64 V1 = 0;
int MINq, p;
for (int i = PI[x14] + 1; i <= PI[x13]; i++){
p = Prime[i];
MINq = min((uint)y, x / ((int64)p * p));
#ifdef PRINT_DEBUG
assert(p > x14 && p <= x13 && MINq >= p);
#endif
V1 += (int64)(2 - PI[p]) * (PI[MINq] - PI[p]); //!!!!!!
}
printf("V1 = %I64d
", V1);
return V1;
}
int64 cal_V2( )
{
int64 V2 = 0;
int i, sqrt_xdivy = (int)sqrt(x / y);
// int sqrt_xdivy2 = (int)sqrt(x / (y * y));
// int sqrt_xdivp = (int)sqrt(x / 1);
int xdivy2 = x / (y * y);
#if 0
uint p, q;
uint W[6] = {0};
//W1
for (p = x14 + 1; p <= xdivy2; p++)
for (q = p + 1; q <= y; q++)
W[1] += PI[x / (p * q)];
//W2
for (p = xdivy2 + 1; p <= sqrt_xdivy; p++){
sqrt_xdivp = (uint)sqrt(x / p);
for (q = p + 1; q <= sqrt_xdivp; q++)
W[2] += PI[x / (p * q)];
}
//W3
for (p = xdivy2 + 1; p <= sqrt_xdivy; p++)
for (q = sqrt(x / p) + 1; q <= y; q++)
W[3] += PI[x / (p * q)];
//W4
for (p = sqrt_xdivy + 1; p <= x13; p++){
sqrt_xdivp = (uint)sqrt(x / p);
for (q = p; q <= sqrt_xdivp; q++)
W[4] += PI[x / (p * q)];
}
//W5
for (p = sqrt_xdivy + 1; p <= x13; p++)
for (q = sqrt(x / p) + 1; q <= xdivy2; q++)
W[5] += PI[x / (p * q)];
V2 = W[1] + W[2] + W[3] + W[4] + W[5];
#else
for (i = PI[x14] + 1; i <= PI[sqrt_xdivy]; i++){
int64 p = Prime[i];
#ifdef PRINT_DEBUG
assert(p > x14 && p <= sqrt_xdivy);
#endif
for (int j = PI[p] + 1; j <= PI[y]; j++){
int q = Prime[j];
#ifdef PRINT_DEBUG
assert(q > p && q <= y);
#endif
V2 += PI[x / (p * q)];
}
}
for (i = PI[sqrt_xdivy] + 1; i <= PI[x13]; i++){
int64 p = Prime[i];
#ifdef PRINT_DEBUG
assert(p > sqrt_xdivy && p <= x13);
#endif
xdivy2 = x / (p * p);
for (int j = PI[p] + 1; j <= PI[xdivy2]; j++){
int q = Prime[j];
V2 += PI[x / (p * q)];
}
}
#endif
printf("V2 = %I64d
", V2);
return V2;
}
int main(int argc, char* argv[])
{
clock_t start = clock();
if (argc > 1)
MAXN = atoll(argv[1]); //mingw
init_x( );
sieve( );
init_PI( );
init_minFactor( );
init_phiTable( );
int64 pix = PI[y] - 1;
pix += cal_S3( ); //mul thread
pix -= cal_P2xa( );
pix += cal_S0( );
if (y != x13){
pix += cal_S1( );
pix += cal_U( );
}
pix += cal_V2( );
pix += cal_V1( );
#if 0
printf("phi(%u, 209) = %d, time = %ld ms
", x, phi(x, 209), clock() - start);
delete phiTable;
return 0;
#endif
puts("
Pi(x) = phi(x, a) + a - 1 - phi2(x,a):");
printf("pi(%I64d) = %I64d, time use %ld ms
", x, pix, clock() - start);
freememory(1);
return 0;
}
//benchmark, Windows Xp SP2, 1G, Mingw G++ 3.4.5
// pi( 2147483647) = 105097565, time use 108 ms (machine AMD 3600+ 2.0G)
// pi(100000000000) = 4118054813, time use 2171 ms (machine intel PD 940 3.2G)最後に論文を添付します.....
次はCSDNの上にもう1篇の比較的に細かい文章が素数を分析します:
http://blog.csdn.net/hexiios/article/details/4400068
下にも貼って、みんなで見学
素数を探して自然数Nより小さいすべての素数を求めます.
ソリューション
プログラム1-1クラシックアルゴリズム
古典的な素数判定アルゴリズムは、正の整数nを与え、2〜sqrt(n)の間のすべての整数でnを除去し、除去可能であればnは素数ではなく、除去できなければnは素数である.したがって、Nより小さいすべての素数を求めるプログラムは以下の通りである.
#include
#include
#define N 1000000
int main(int argc, char *argv[])
{
int m, n;
for (n =2 ; n < N; n++)
{
for (m = 2; m * m <= n; m++)
if (n % m == 0) break;
if (m * m > n) printf("%6d",n);
}
system("PAUSE");
return 0;
アルゴリズムの時間的複雑さはO(N*sqrt(N))であり,解の規模の小さい入力は問題ない.しかし規模の大きいNに対しては,アルゴリズムは力不足である.エラドセふるい(sieve of Eratosthenes)というアルゴリズムがあり,解く際にかなりの効率を有するが,より多くの空間を犠牲にする.
プログラム1-2アラマルチプラグふるいアルゴリズム
このプログラムの目標は,iが素数であればa[i]=1を設定することである.素数でない場合は0に設定します.まず、すべての配列要素を1に設定し、既知の非素数がないことを示します.次に、非素数(既知の素数の倍数)として知られているインデックスに対応する配列要素を0に設定し、すべての小さい素数の倍数を0に設定した後もa[i]が1のままであれば、探した素数であると判断することができる.
#include
#include
#define N 1000000
int a[N];
int main(int argc, char *argv[])
{
int i, j;
for (i = 2; i < N; i++) a[i] = 1;
for (i = 2; i < N; i++)
{
if (a[i])
for (j = i + i; j < N; j += i)
a[j] = 0;
}
for (i =2; i < N; i++)
if (a[i]) printf("%6d ",i);
system("PAUSE");
return 0;
}
たとえば、32未満の素数を計算すると、まずすべての配列項目を1(次の2列目)に初期化し、非素数が発見されるごとの数を表します.次に、2、3、5倍のインデックスを持つ配列項目を0に設定します.2、3、5倍の数が非素数であるためです.1に保持されている配列項目は、素数(次の右列)に対応するインデックスを付けます.
i 2 3 5 a[i]
2 1 1
3 1 1
4 1 0
5 1 1
6 1 0
7 1 1
8 1 0
9 1 0
10 1 0
11 1 1
12 1 0
13 1 1
14 1 0
15 1 0
16 1 0
17 1 1
18 1 0
19 1 1
20 1 0
21 1 0
22 1 0
23 1 1
24 1 0
25 1 0
26 1 0
27 1 0
28 1 0
29 1 1
30 1 0
31 1 1
エラドセふるいアルゴリズムをどのように理解しますか?
小学校の時に学んだ素数と合数の性質を覚えているに違いない.いずれの合数もいくつかの素数の積として表すことができる.たとえば、4=2*2、6=2*3、12=2*2*3です.すなわち、合数Nは必ずNより小さい(またはいくつかの)素数の整数倍であり、逆に、Nがそれより小さい素数の倍数でなければ、Nは必然的に素数である.
プログラム1-3クラシックアルゴリズムの改良バージョン
古典的な素数判定アルゴリズムでは,2からsqart(n)までのすべての整数でnを除去する必要はなく,その間の素数だけで十分であり,合数が必ずいくつかの素数の積として表されるためである.アルゴリズムの改良バージョンは次のとおりです.
#include
#include
#define N 1000000
int prime[N];
int main(int argc, char *argv[])
{
int i, n, flag;
prime[0] = 0;//素数の合計を保存する
for (n =2 ; n < N; n++)
{
flag = 1;
for (i = 1; i <= prime[0] && prime[i] * prime[i] <= n; i++)
if (n % prime[i] == 0)
{
flag = 0;
break;
}
if (flag)
{
prime[0]++;
prime[prime[0]] = n;
printf("%6d ",n);
}
}
system("PAUSE");
return 0;
}
アルゴリズムは効率的に以前のバージョンより向上したが,決定された素数を記憶するために空間を犠牲にする必要がある.
プログラム1-4アラドプラグふるいアルゴリズムの改良版
プログラム1-2は、最も簡単な要素タイプ、0と1の2つの値を含む配列を使用し、整数の配列ではなくビットの配列を使用すると、より高い空間的有効性を得ることができます.
#include
#include
#define N 1000000
#define BITSPERWORD 32
#define SHIFT 5
#define MASK 0x1F
#define COUNT 1+N/BITSPERWORD//Bit i対写数i
void set(int a[],int i);
void clr(int a[],int i);
int test(int a[],int i);
int a[COUNT];
int main(int argc, char *argv[])
{
int i, j;
for (i = 0; i < COUNT; i++) a[i] = -1;
for (i = 2; i < N; i++) {
if (test(a,i))
for (j = i + i; j < N; j += i)
clr(a,j);
}
for (i =2; i < N; i++)
if (test(a,i)) printf("%4d ",i);
system("PAUSE");
return 0;
}
void set(int a[],int i) {
a[i >> SHIFT] |= (1<<(i & MASK));
}
void clr(int a[],int i) {
a[i >> SHIFT] &= ~(1<<(i & MASK));
}
int test(int a[],int i) {
return a[i >> SHIFT] & (1<<(i & MASK));
}