//Sum both twins < n using a twin index file to read records from file
//h:\sievedata\prime3-1trill.bin which was created by the sieve 30k+r
//Achim wheel program modified in f:\eratosieve\primesieve8byte.c
//This replaces the basic program referred to in the link which gave the
//wrong result for the 10th term which should be
//Sum of 524733511 both twin primes < 252097800623 = 126670338180923086268
//
#include <windows.h>
#include <stdio.h>
#define uns64 unsigned long long
#define chr64(a) ulltoa(a,str,10)
#define timer GetTickCount()/1000.0

static char    fil[50],fil2[50];
static HANDLE  H;
static DWORD   Read;
static int     recl;
long  high;
float t;
static FILE*  fp1;
static FILE*  fp2;
unsigned long long offset,b;
unsigned long long p1,p2;
unsigned long long c,x,end;
unsigned long long n;
char start1[30],end1[30],str[30];
long start,j;

char sump[50],gp[50];
uns64 twin(long);
int paricalc(void);
unsigned long long s1,s2,s3,s4;

int main(int argc, char *argv[])
{
strcpy((char*)fil,"h:\\sievedata\\twin3-1trill.bin");  //Lower twin primes < 1 trillion
 printf("Sum both twin primes < n. Read a file of lower twin primes < 10^12\n");
  printf("\n");
  H=CreateFile(fil,GENERIC_READ,FILE_SHARE_READ,NULL,OPEN_EXISTING,0,0);
  b=536870913;
  recl=8;
long double sum;
while(1)
{
sum=0;s1=0,s2=0;s3=0;s4=0;
c=0;
//  printf("Number => ");
//  gets(start1);
//  start=atoll(start1);
start = 1870585220;
  printf("Max   => ");
  gets(end1);
  end=atoll(end1);
t=timer;
for(j=1;j<=start;j++)
  {
  n = twin(j);
if(n > end)break;
x=n+n+2;
//  printf("%s ",chr64(n));
    s1+= x/10000;
     x = x%10000;
    s2+= x/1000;
     x = x%1000;
    s3+= x/100;
    s4+= x%100;
  sum+=x;
  c++;
  }
  printf("\n");
paricalc();
  printf("Time = %.06f\n",timer-t);
}//wend
  CloseHandle(H);
  return 0;
}

uns64 twin(long m)
  {
  high=m/b;
  offset=recl*(m-1);
  SetFilePointer(H,offset,&high,FILE_BEGIN);
  ReadFile(H,&p1,recl,&Read,NULL);
  return p1;
  }
int paricalc(void)
     {
      sprintf(sump,"%.18G0000+%.18G000+%.18G00+%.18G",(double)s1,(double)s2,(double)s3,(double)s4);
      fp1=fopen("fgp.gp","w");
      fprintf(fp1,"%s\n",sump);
      if(fp1)
       {
         fflush(fp1);
         fclose(fp1);
       }
      strcpy((char*)gp,"gp.exe -f -q < fgp.gp > pdata.txt");
//      printf("%s %s\n","Pari string ",gp);
      system(gp);
      fp2=fopen("pdata.txt","r");
      sump[0]=0;
      fgets(sump,50,fp2);
       printf("Sum of %s both twin primes < %s = %s\n",chr64(c),end1,sump);
      if(fp2)
       {
         fflush(fp2);
         fclose(fp2);
       }
return 0;
    }