seqfun Yahoo Group Count primes in a range up to 10^18 =============================================== Cino Hilliard Message 1 of 1 Sep 23, 2007 ----------------------------------------------- //gcc 4.1.2 //Countprimes uses Achim Flammenkamp's algorithm below. //Count the primes from start and stop up to 2^64. //Comments below with Cino show some of my additions. // Cino Hilliard // Sept 2007 // http://groups.yahoo.com/group/seqfun/files/ #include <stdio.h> /* Achim Flammenkamp, Bielefeld, D, achim@... */ #include <malloc.h> /* ifdefs: LONG=uns64 _ANSI_EXTENSION true_64bit_words */ #include <stdlib.h> /* SUM_{p<x} 1/p= ln(ln(x))+0.5772156649015-0.315718451893 */ #include <math.h> /* extern double sqrt(double), log(double), floor(double); */ #include <string.h> #include <windows.h> #define true_64bit_words 0 // Cino This is necessary #define _ANSI_EXTENSION 0 // So is this #ifndef uns32 /* Version 1.0a 1998-05-19 source availible */ #define uns32 unsigned /* Version 2.0a 1998-05-23 5 years in the WWW */ #endif /* Version 2.0b 2003-06-03 0-1-typo fixed Argh! */ #ifndef uns64 /* Version 2.0c 2003-06-04 works up to 2^64-1 */ # ifdef _ANSI_EXTENSION /* uses < 1938 MB for minimal sieve size to 2^64 */ # define uns64 unsigned long long #define LONG uns64 # else # define uns64 unsigned long # endif #endif #ifdef LONG # ifdef _ANSI_EXTENSION # define UL "%llu" # else # define UL "%lu" # endif # define HALF uns32 #else # define LONG uns32 # define UL "%u" # define HALF unsigned short #endif #define BITS 8 #define LCM (2*3*5) #define SIEVE_SIZE (16<<17) #define chr64(a) _ui64toa(a,str,10) #define ld long double #define digits 40 //Cino add the timer #define timer GetTickCount()/1000.0 long u,v,ai,j,aa[50],s[50],ln,lna,lnb; unsigned long long p[50],cr,z,x,y,xx; static char a1[50],a2[50],str[50]; static FILE* fp1; static FILE* fp2; LONG anz,pi2; static long a[]={1,2,4,8,16,32,64,128}; static long b1[]={1,7,11,13,17,19,23,29}; static long r; long digits2; static char *pEnd; char stop1[30],start1[30]; int paricalc(void); long double sum,est,xld; unsigned long long p[50],s1[20],p1[9],x0; LONG ii2; LONG stop, start, stlcm, ii, jj, hh; int main(unsigned argc, char *argv[]) { float t; LONG size, hi, h, i, ji, js; uns32 k, hj, c, cj; LONG j,*index; unsigned char *sieve, *prime_diff; unsigned char bi, b, m; HALF psize; int lim; printf(" Count primes in a range from m to n up to 10^18 \n"); while(1) { sum=0; anz=0; pi2=0; digits2=0; for(u=0;u<digits;u++)s[u]=0; printf("start => "); gets(start1); start = strtoull(start1,&pEnd,10); printf("stop => "); gets(stop1); stop = strtoull(stop1,&pEnd,10); //for(ai=ln;ai>=0;ai--) //printf("%i",a1[ai]); size = SIEVE_SIZE; if (stop < start) goto finish; t = timer; system("date/t"); system("time/t"); j =(LONG)sqrt((long double)stop+0.1); if (!(j&1)) j--; if ((j-1)/2 >= (5-1)/2*(LCM/5)*size) /* to prevent overflow if j == 2^32-1 */ size = (j/48+1)*2; if (!(sieve = (unsigned char*)calloc(size, sizeof(unsigned char)))) return fprintf(stderr,"error: can't get %u kB storage for sieve\n", ((unsigned)sizeof(unsigned char)*size+(1<<10)-1)>>10), 2; psize = (int)((1.015*j)/(log((double)j)-1)); /* upper bound for Pi(stop^1/2) */ if (!(index = (LONG*) malloc((sizeof(LONG)+1)*psize))) return fprintf(stderr,"error: can't get storage for %u primes\n",psize), 4; //printf("Memory Used %d\n",index); prime_diff = (unsigned char*)(index+psize); if (start <= 2){anz++;} //{} are necessary! if (start <= 3){anz++;} if (start%2 == 0)start += 1; if (start%3 == 0)start += 2; stlcm = start/LCM; hh = size*(stlcm/size); /*** sifting the sieve primes ***/ k=0; m=0; h= i= cj= 0; js=jj=0; b=1; c=1; for(;;) { switch (c) { do { case 1: i++; if ((m+=1) >= LCM/3) m -= LCM/3; jj += h; h++; jj += h; if (!(sieve[i]&b)) { c= 2; break; } case 2: i++; if (i == size) { i=0; cj += size; b += b; if (!b) break; } if ((m+=1) >= LCM/3) m -= LCM/3; jj += h; if (!(sieve[i]&b)) { c= 1; break; } } while (1); } if (8*jj > stop/3) break; j = 3*(cj+i)+c; bi= m - !!m - m/BITS; prime_diff[k]= BITS*2*(j/LCM-js); /* difference of successive primes < 480 */ prime_diff[k] |= bi; js= j/LCM; index[k]= ((bi&1) != (bi>>2&1) ? 5 : 0); ii= (8*jj)/(LCM/3); if (ii < stlcm) ii += j*((stlcm-ii)/j); if (ii < stlcm) { hi = (index[k] ? 19 : 1); /* inverse zu index[k] --- 0 <--> 1 Argh!! */ hj= js+js; ji= 2*(3*m+c); if (ji >= LCM) ji-=LCM, hj++; switch (c) { do { case 1: ii += 2*hj; hi += 2*ji; while (hi >= LCM) hi -= LCM, ii++; if (ii >= stlcm && hi!=5 && hi!=25) break; case 2: ii += hj; hi += ji; while (hi >= LCM) hi -= LCM, ii++; if (ii >= stlcm && hi!=5 && hi!=25) break; } while(1); } index[k] = (BITS*hi)/LCM; } index[k] |= BITS*(ii-hh); k++; if (jj >= size) continue; /*** sift with prime pre-sieve ***/ #ifdef true_64bit_words ii = 8*jj; #define hi ii #else hi = 8*jj; #endif ji= 2*(cj+i)+1; hj= 2*(ji+c)-3; bi= 1; switch(m) { do { case 0: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += 2*j; case 2: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += hj; case 3: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += ji; case 4: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += hj; case 5: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += ji; case 6: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += hj; case 7: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += 2*j; case 9: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += ji; lim= size - LCM/3*j; while ((int)hi < lim) { sieve[hi] |= bi; hi += 2*j; sieve[hi] |= bi; hi += hj; sieve[hi] |= bi; hi += ji; sieve[hi] |= bi; hi += hj; sieve[hi] |= bi; hi += ji; sieve[hi] |= bi; hi += hj; sieve[hi] |= bi; hi += 2*j; sieve[hi] |= bi; hi += ji; } } while (1); do { case 1: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += ji; case 8: while(hi >= size) { hi -= size; bi+=bi; if (!bi) goto go_on; } sieve[hi] |= bi; hi += hj; lim= size - LCM/3* 5; while ((int)hi < lim) { sieve[hi] |= bi; hi += ji; sieve[hi] |= bi; hi += hj; sieve[hi] |= bi; hi += ji; sieve[hi] |= bi; hi += hj; sieve[hi] |= bi; hi += ji; sieve[hi] |= bi; hi += hj; sieve[hi] |= bi; hi += ji; sieve[hi] |= bi; hi += hj; sieve[hi] |= bi; hi += ji; sieve[hi] |= bi; hi += hj; } } while (1); } go_on: ; #ifdef true_64bit_words #undef hi #endif } /****** main sifting starts *****/ for (i=size;i--;) sieve[i]=0; sieve[0] |= !hh; /* 1 is not prime */ if (start <= 5){anz++;} if (start%5 == 0) start += 2*(3-start%3); hh *= LCM; while (1) { j= prime_diff[0]/BITS; for (h=1;h<k;h++) /* sieve with next sieve prime (h=0 is prime 5) */ { j += prime_diff[h]/BITS; ii = index[h]/BITS; if (ii >= size) { index[h] -= size*BITS; continue; } hj= (size <= LCM/2*j ? 0 : size - LCM/2*j); i=ji=js= j; ji +=js; i += ji; #ifdef true_64bit_words #define hi ii #else hi = ii; #endif //Cino - equated index[h] in place and added break in place of goto switch( prime_diff[h]%BITS ) { case 0: switch( index[h]%BITS ) { do { case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += i; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += ji; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += js; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += ji; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += js; case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += ji; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += i; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += js+1; lim= hj-1; while ((int)hi < lim) { sieve[hi] |= 1; hi += i; sieve[hi] |= 2; hi += ji; sieve[hi] |= 4; hi += js; sieve[hi] |= 8; hi += ji; sieve[hi] |= 16; hi += js; sieve[hi] |= 32; hi += ji; sieve[hi] |= 64; hi += i; sieve[hi] |=128; hi += js+1; } } while (1); } case 1: js+=1; i+=1; ji+=1; switch( index[h]%BITS ) { do { case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += ji; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += js; case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += ji-1; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += js; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += ji; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += i; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += js; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += i; lim= hj-7; while ((int)hi < lim) { sieve[hi] |= 32; hi += ji; sieve[hi] |= 16; hi += js; sieve[hi] |= 1; hi += ji-1; sieve[hi] |=128; hi += js; sieve[hi] |= 8; hi += ji; sieve[hi] |= 4; hi += i; sieve[hi] |= 64; hi += js; sieve[hi] |= 2; hi += i; } } while (1); } case 2: i+=2; ji+=2; switch( index[h]%BITS ) { do { case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += js; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += ji; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += js; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += ji; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += i; case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += js+1; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += i; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += ji; lim= hj-11; while ((int)hi < lim) { sieve[hi] |= 1; hi += js; sieve[hi] |= 64; hi += ji; sieve[hi] |= 2; hi += js; sieve[hi] |=128; hi += ji; sieve[hi] |= 8; hi += i; sieve[hi] |= 32; hi += js+1; sieve[hi] |= 4; hi += i; sieve[hi] |= 16; hi += ji; } } while (1); } case 3: js+=1; i+=3; ji+=1; switch( index[h]%BITS ) { do { case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += ji+1; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += js; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += ji; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += i; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += js; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += i; case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += ji; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += js; lim= hj-13; while ((int)hi < lim) { sieve[hi] |= 32; hi += ji+1; sieve[hi] |= 4; hi += js; sieve[hi] |= 2; hi += ji; sieve[hi] |=128; hi += i; sieve[hi] |= 16; hi += js; sieve[hi] |= 8; hi += i; sieve[hi] |= 1; hi += ji; sieve[hi] |= 64; hi += js; } } while(1); } case 4: js+=1; i+=3; ji+=3; switch( index[h]%BITS ) { do { case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += js; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += ji; case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += i; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += js; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += i; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += ji; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += js; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += ji-1; lim= hj-17; while ((int)hi < lim) { sieve[hi] |= 32; hi += js; sieve[hi] |= 64; hi += ji; sieve[hi] |= 1; hi += i; sieve[hi] |= 8; hi += js; sieve[hi] |= 16; hi += i; sieve[hi] |=128; hi += ji; sieve[hi] |= 2; hi += js; sieve[hi] |= 4; hi += ji-1; } } while (1); } case 5: js+=2; i+=4; ji+=2; switch( index[h]%BITS ) { do { case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += ji; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += i; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += js-1; case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += i; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += ji; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += js; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += ji; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += js; lim= hj-19; while ((int)hi < lim) { sieve[hi] |= 1; hi += ji; sieve[hi] |= 16; hi += i; sieve[hi] |= 4; hi += js-1; sieve[hi] |= 32; hi += i; sieve[hi] |= 8; hi += ji; sieve[hi] |=128; hi += js; sieve[hi] |= 2; hi += ji; sieve[hi] |= 64; hi += js; } } while (1); } case 6: js+=1; i+=5; ji+=3; switch( index[h]%BITS ) { do { case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += i; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += js; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += i; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += ji; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += js; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += ji+1; case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += js; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += ji; lim= hj-23; while ((int)hi < lim) { sieve[hi] |= 32; hi += i; sieve[hi] |= 2; hi += js; sieve[hi] |= 64; hi += i; sieve[hi] |= 4; hi += ji; sieve[hi] |= 8; hi += js; sieve[hi] |=128; hi += ji+1; sieve[hi] |= 1; hi += js; sieve[hi] |= 16; hi += ji; } } while (1); } case 7: js+=2; i+=6; ji+=4; switch( index[h]%BITS ) { do { case 0: if (hi >= size) {index[h] = 0; break;} sieve[hi] |= 1; hi += js-1; case 7: if (hi >= size) {index[h] = 7; break;} sieve[hi] |=128; hi += i; case 6: if (hi >= size) {index[h] = 6; break;} sieve[hi] |= 64; hi += ji; case 5: if (hi >= size) {index[h] = 5; break;} sieve[hi] |= 32; hi += js; case 4: if (hi >= size) {index[h] = 4; break;} sieve[hi] |= 16; hi += ji; case 3: if (hi >= size) {index[h] = 3; break;} sieve[hi] |= 8; hi += js; case 2: if (hi >= size) {index[h] = 2; break;} sieve[hi] |= 4; hi += ji; case 1: if (hi >= size) {index[h] = 1; break;} sieve[hi] |= 2; hi += i; lim= hj-29; while ((int)hi < lim) { sieve[hi] |= 1; hi += js-1; sieve[hi] |=128; hi += i; sieve[hi] |= 64; hi += ji; sieve[hi] |= 32; hi += js; sieve[hi] |= 16; hi += ji; sieve[hi] |= 8; hi += js; sieve[hi] |= 4; hi += ji; sieve[hi] |= 2; hi += i; } } while (1); } hi -= size; index[h] |= BITS*hi; } } /* output of remaining (prime) numbers cino version with no switch error*/ char str[24],str1[24]; i=(start<=hh ? 0 : (start-hh)/LCM); hh += LCM*i; bi= sieve[i]; for (;i<size;i++) { bi= sieve[i]; //This gets the next 8 numbers 30k+r, r=1,7,11,13,17,19,23,29 for(r=0;r<8;r++) { if(!(bi&a[r])) //test for 0 bit { ii=hh+b1[r]; //track the numbers in the wheel if (ii > stop)goto end; if(ii >= start&&ii<=stop) //absolutely necessary! { anz++; } } } hh += LCM; sieve[i]=0; } } end: free(index); free(sieve); finish: t=timer-t; digits2=digits; printf("# {%s <= Primes <= %s} = %s\n",start1,stop1,chr64(anz)); printf("Sec = %.06f\n",t); } return anz; } =============================================== Cached by Georg Fischer at Nov 14 2019 12:47 with clean_yahoo.pl V1.4