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