seqfun Yahoo Group Sum of Isolated primes =============================================== Cino Hilliard Message 1 of 1 Jun 5, 2007 ----------------------------------------------- Hi, Here is a Gcc program that computes the sum of isolated primes The listing and the windows xp pro p4 exe file is in the files section. Yahoo makes you join to read that section. Cino // sumIsoAchim uses Achim Flammenkamp's algorithm below. Sum the Isolated non // twin primes from start and stop up to 2^64. The routine splits the sums into // smaller parts for a Pari call to handle big numbers whose sum > 10^16. In // addition the estimated sum = stop^2/(2*log(stop)-1)- start^2/(2*log (start)-1) // is computed for comparison. Comments below with Cino show some of my additions. // Cino Hilliard // Oct 2006,June 2007 // Feb 2007 compile with c:\cygwin\bin\g++ but must use sprintf. ulltoa not allowed // http:// groups.yahoo.com/group/seqfun/files/ #include <stdio.h> /* FAchimlammenkamp, Bielefeld, D, achim@uni- bielefeld.de */ #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 1 // 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 fact .2565 #define chr64(a) ulltoa(a,str,10) // Cino save partial sums to be added by Pari #define use(x) do {\ anz++;\ x1=x;\ s4+=x1/1000000;\ x1= x1%1000000;\ s3+=x1/100000;\ x1= x1%100000;\ s2+=x1/10000;\ x1= x1%10000;\ s1+=x1/1000;\ s0+=x1%1000;\ } while(0); // ulltoa(x,buffer,10);\ // printf("%s ",buffer);\ // Cino add the timer #define timer GetTickCount()/1000.0 unsigned long long j,s0=0,s1=0,s2=0,s3=0,s4=0,s5=0,s6=0,s7=0,s8=0,s9=0,s10=0,x1,ct; static FILE* fp1; static FILE* fp2; static FILE* fp3; char buffer[200]; 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; char str[24],str1[24]; char *binary(long); int paricalc(); char start1[20],stop1[20]; LONG anz,pi2,start,stop,len; uns64 lastprime; int main(unsigned argc, char *argv[]) { float t; uns32 k, hj, c, cj; unsigned char *sieve, *prime_diff; unsigned char bi, b, m; LONG stlcm, *index,ii,ii2,ii3,jj, hh; LONG size, hi, h, i, ji, js; HALF psize; int lim,ln; printf(" Sieve, Compute the sum of twin primes in a range primes between 1 and 10^18\n"); printf(" Cino Hilliard\n"); printf(" June, 2007\n"); printf(" Input start and stop as prompted below\n"); printf("\n"); while(1) { j=0;s0=0;s1=0;s2=0;s3=0;s4=0;s5=0;s6=0;s7=0;s8=0;s9=0;s10=0;x1; anz=0;pi2=0;*start1=0;*stop1=0;start=0;stop=0;size=0; ii2=5;ii3=3; printf("start => "); gets(start1); start = atoll(start1); printf("stop => "); gets(stop1); stop = atoll(stop1); t = timer; j =(LONG)sqrt((long double)stop+0.1)+10; printf("Sqrt stop = %llu\n",j); if (!(j&1)) j--; size = (j/48+1)*2; printf("Extending sieve size to %d\n",size); 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 %llu\n",index); prime_diff = (unsigned char*)(index+psize); if (start <= 3) use(3); 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)use(5); 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 out0-out7. 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 ***/ 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! { if(ii2+2!=ii&&ii3+2!=ii2) // If it is not a twin part, add it. { use(ii2); // printf("%s ",chr64(ii2)); } ii3=ii2; ii2=ii; } } } hh += LCM; sieve[i]=0; } } end: printf("%s ",chr64(ii3)); printf("%s ",chr64(ii2)); printf("%s\n",chr64(ii)); if(ii2+2!=ii&&ii3+2!=ii2)use(ii2); s0-=6; free(index); free(sieve); finish: t=timer-t; // Cino - Estimate of sum of primes in an interval start stop printf("# {%s <= Iso primes <= %s} = %.18G\n",start1,stop1,(double) anz); // printf("%s\n","Sum estimate = stop^2/(2*log(stop)-1) - start^2/ (2*log(start)-1)"); paricalc(); printf("Sec = %.6f\n",t); }// wend(1) return anz; } int paricalc() { char sump[200],gp[100],logdata[100]; long ln; // Cino Here we combine the partial sums and call Pari to get the total. // This assumes gp.exe is in the path. sprintf(sump,"%.18G000000+%.18G00000+%.18G0000+%.18G000+%.18G", (double)s4,(double)s3,(double)s2,(double)s1,(double) s0); printf("%s\n %s\n" ,"String to fgp.gp for Pari to evaluate. ",sump); if((fp1=fopen("fgp.gp","w"))==0) { fprintf(stderr,"Can't open file %s\n","fgp.gp");exit(1); } fprintf(fp1,"%s\n",sump); if(fp1) { fflush(fp1); fclose(fp1); } strcpy((char*)gp,"gp.exe -q -f < fgp.gp > pdata.txt"); system(gp); if((fp2=fopen("pdata.txt","r"))==0) { fprintf(stderr,"Can't open file %s\n","pdata.txt");exit(1); } sump[0]=0; fgets(sump,50,fp2); ln = strlen(sump); fseek(fp2,0,0); fgets(sump,ln,fp2); if(fp2) { fflush(fp2); fclose(fp2); } printf("%s%s\n","Sum of Iso or non twin primes = ",sump); return 0; } char *binary(long n) // return the binary representation of n { char *BCX_RetStr={0}; char str[30]; char *zeros; zeros = (char*) calloc (10,sizeof(char)); // This is a must. char zeros[10] will not work; char z[10][10] = {"","0","00","000","0000","00000","000000","0000000"}; long ln; { itoa(n,str,2); ln = 8-strlen(str); strcpy(zeros,z[ln]); // set preceeding zeros of the binary strcat(zeros,str); } return zeros; // Return the number in base r1 to base r2 } =============================================== Cached by Georg Fischer at Nov 14 2019 12:47 with clean_yahoo.pl V1.4