#include <stdio.h>
#include <math.h>


long double psi(long double);
long double rs3(long double);
long double z(long double);
long double R3(long double);
long double d1psi(long double);
long double d2psi(long double);
long double d3psi(long double);
long double d6psi(long double);
long double d9psi(long double);
long double theta(long double);
long double u(long double);
long double p(long double);
int N(long double);

int blocksBegin[205000000];
int blocksEnd  [205000000];
int blocksZeros[205000000];

long double ep = 0.00064L ;

long double pi = 3.14159265358979323846L ;



int rs3_uses = 0;



int main(void)
{

long double t;
long double last;
long double new;
int count;
long double delta = 0.005L;
long double deriv;
long double guess;
int k;
int j;
int jmax;
int l;
int iter;
int found;
int index1;
int index2;
long double g1;
long double g2;
int zeroscount;
int multiplier;
int quit;
int ngriddivisor;
long double del;
long double rs3guess;
long double rs3g1;
long double rs3g2;
int L;
int blockcount[100];
int gramitv; 
int newmultiplier;
int nmulti;
int blocknumber = 0;
int rosser_rule_exceptions = 0;
int rreindex;
int rosser_array[1000];
int kk;
long double y;
int success;
int missing;
int answer;
long double tmin;
long double tmax;
long double tmid;
long double left;
long double right;







for(j=0; j<= 100; j++)
{
  blockcount[j] = 0;
}



// 40870156


//  43628107


// 46082042


// 49,624,541

// 50,799,238


//  55,221,454

// 56,948,780

//  60,515,663

//  61,331,766

// 69,784,844


//  jmax  = 75000000;

//  jmax  = 70000000;

  jmax = 75000000 ;




zeroscount = 1;

found = 0;


tmin=10.0L;


//    tmax = 35000000.0L;

tmax = 90000000.0L;



// 37592215


for(j = 0; j< jmax ; j++)
{ 


left = tmin;
right = tmax;



for(iter=1; iter <= 80; iter++)
{

  tmid = (left+right)/2.0L;

  y = theta(tmid) - ((long double) j)*pi;

  if( y > 0.0L)
  {
     right = tmid;
  }
  else
  {
     left = tmid;
  } 

 //   printf("theta(%.16Lf) = %.16Lf\n", tmid, theta(tmid) );

  

//   deriv = (theta(guess+delta) - theta(guess-delta))/(2.0L*delta);

//  guess =guess - ( theta(guess) - ((long double)j)*pi )/deriv;

}



guess = tmid;




rs3guess = rs3(guess);








  if( rs3guess > 0.0L && (1 == (j%2))  && j==(jmax-1)  )
  {
    printf("bad gram point at %d\n", j);
  }

  if( rs3guess > 0.0L && (1 == (j%2))   )
  {
    L = L+1;
  }




  if( rs3guess > 0.0L && (0 == (j%2)) )
  {

/************************************************
      if( 0 == (j%100000))
      {
         printf("good gram point at %d\n", j);


         printf("The cumulative zeros count is: %d\n\n\n", zeroscount);
 
      }
********************************************/



    if(found == 1)
    { 
       found = 2;
       L=L+1;
       index2 = j;
       g2 = guess;
       rs3g2=rs3guess;
    }



    if(found == 0)
    { 
       found = 1;
       L = 0;
       index1 = j;
       g1 = guess;
       rs3g1 = rs3guess;
    }


      
  }



  if( rs3guess < 0.0L && (0 == (j%2))   && j==(jmax-1)    )
  {
    printf("bad gram point at %d\n", j);
  }

  if( rs3guess < 0.0L && (0 == (j%2))      )
  {
   // printf("bad gram point at %d\n", j);
   L = L+1;
  }





 if( rs3guess < 0.0L && (1 == (j%2)) )
  {
   // printf("good gram point at %d\n", j);

    if(found == 1)
    { 
       found = 2;
       L = L+1;
       index2 = j;
       g2 = guess;
       rs3g2=rs3guess;
    }


    if(found == 0)
    { 
       found = 1;
       L = 0;
       index1 = j;
       g1 = guess;
       rs3g1 = rs3guess;
    }



  }



if( found == 2)
{
  /*** test the Gram block ***/



  blockcount[L] ++ ;

  blocknumber++;



  if(0 == (j%1000000)     )
  {
    gramitv = 0;
    printf("Gram point %d:  %d calls to rs3\n", j,  rs3_uses);
    printf("Distribution of Gram blocks with lengths 1 to 8:\n");
    for(k=1; k<= 8; k++)
    {
      printf("length %d:  %d blocks\n", k, blockcount[k] );
      gramitv = gramitv + k*blockcount[k];
    }
    printf("Total of %d Gram intervals\n", gramitv);
    printf("\n\n\n");
  }



if( index2-index1 == 1)
{
   zeroscount = zeroscount+1;





/********** paste below *******/


      if( 0 == (j%100000))
      {
         printf("good gram point at %d\n", j);
         printf("The cumulative zeros count is: %d\n\n", zeroscount); 
      }


blocksBegin[blocknumber-1] = index1;

blocksEnd[blocknumber-1] = index2;

blocksZeros[blocknumber-1] = 1;

}



if( index2-index1 > 1)
{
   

multiplier = 1;
newmultiplier = 2;


quit = 0;

while( quit == 0)
{




ngriddivisor = multiplier*(1+index2-index1);

del = (g2-g1)/ ((long double) ngriddivisor    ) ;


last = rs3g1;

t = g1;

count = 0;

for( l = 1; l <= ngriddivisor; l++)
{
  t = t+del;

  new = rs3(t);
  if(new*last < 0.0L)
  {
    count = count + 1;
  }
  last = new;
}

  if( count < index2-index1 )
  {
     nmulti = multiplier+newmultiplier;
     multiplier = newmultiplier; 
     newmultiplier = nmulti;
  }

  if( count > index2-index1 )
  {
    printf("We have done an overcount...\n\n");

    scanf("%d", &answer);
  }




  if( count == index2-index1 )
  {
     quit = 1;

     zeroscount = zeroscount + count;   

// printf("Gram block from %d to %d, blocknum=%d zeros = %d, in total %d\n", index1, index2, blocknumber, index2-index1, zerosfound );


/****** paste below ***********/


      if( 0 == (j%100000))
      {
         printf("good gram point at %d\n", j);
         printf("The cumulative zeros count is: %d\n\n", zeroscount); 
      }






blocksBegin[blocknumber-1] = index1;

blocksEnd[blocknumber-1] = index2;

blocksZeros[blocknumber-1] =  index2- index1;


  }

  if( multiplier >  1024 )
  {
    quit = 1;
    printf("Rosser rule violation... \n");
    printf("index1 = %d\n", index1);
    printf("index2 = %d\n", index2);

    blocksBegin[blocknumber-1] = index1;
    blocksEnd[blocknumber-1] = index2;
    blocksZeros[blocknumber-1] = count ;

printf("%d zeros are missing from this Gram block\n", index2-index1- count);

//     printf("ready to tackle the missing zeros...\n\n");

   /*****  Rosser rule exception handling  below ***/


/**************************************************** 
printf("First, look in next Gram block ... \n");

printf("Second, look in previous Gram block...\n");

printf("Make a note of the Block number, which is in index 0 notation: %d\n", blocknumber-1);
*******************************************************/


rosser_rule_exceptions ++;

printf("rosser_rule_exceptions = %d\n\n", rosser_rule_exceptions);

rreindex = blocknumber-1;


rosser_array[rosser_rule_exceptions-1] = rreindex;

// printf("rosser_array[%d] = %d \n ", rosser_rule_exceptions-1, rreindex);


// printf("revisit after first pass completes \n\n");


  }

}



}



found = 1;

L = 0;


g1 = g2;

index1 = index2;

rs3g1 = rs3g2;


}


}


printf("\n\ncount is %d  t is %.16Lf\n\n",zeroscount, g2 );



printf("We now revisit Rosser rule exceptions... \n\n");

printf("Rosser's rule says that, in a Gram block of length L >= 1, there will be at least L zeros of the Z function\n");



for(kk = 1; kk <= rosser_rule_exceptions; kk++)
{


   rreindex = rosser_array[kk-1]; 

   printf("exception number %d is for Gram block %d \n", kk, rreindex);

   index1 = blocksBegin[rreindex];

   index2 = blocksEnd[rreindex];

   count = blocksZeros[rreindex];

   missing = index2-index1 - count;


   printf("We found %d missing zeros in the block from gram point %d to gram point %d\n", missing, index1, index2);




   printf("We now look at the following Gram Block\n");


   index1 = blocksBegin[rreindex+1];

   index2 = blocksEnd[rreindex + 1];



   y = pi*((long double) index1);

   guess = y/(logl(y)/2.0L);




for(iter=1; iter <= 15; iter++)
{




deriv = (theta(guess+delta) - theta(guess-delta))/(2.0L*delta);

guess =guess - ( theta(guess) - ((long double)index1)*pi )/deriv;

if(iter == 15)
{
   printf("guess = %.16Lf\n", guess);
}





}

printf("\n");





rs3guess = rs3(guess);


g1 = guess;
rs3g1 = rs3guess;





   y = pi*((long double) index2);

   guess = y/(logl(y)/2.0L);




for(iter=1; iter <= 15; iter++)
{




deriv = (theta(guess+delta) - theta(guess-delta))/(2.0L*delta);

guess =guess - ( theta(guess) - ((long double)index2)*pi )/deriv;

if( iter == 15)
{

   printf("guess = %.16Lf\n", guess); 
}




}   

printf("\n");




rs3guess = rs3(guess);


g2 = guess;
rs3g2 = rs3guess;





if( 1 ==  1)
{
   

multiplier = 1;



quit = 0;

while( quit == 0)
{




ngriddivisor = multiplier*(1+index2-index1);

del = (g2-g1)/ ((long double) ngriddivisor    ) ;


last = rs3g1;

t = g1;

count = 0;

for( l = 1; l <= ngriddivisor; l++)
{
  t = t+del;

  new = rs3(t);
  if(new*last < 0.0L)
  {
    count = count + 1;
  }
  last = new;
}

  if( count < missing + index2-index1 )
  {

     multiplier = 2*multiplier; 

  }

  if( count > missing + index2-index1 )
  {

    printf("We have an overcount...\n");
    scanf("%d", &answer);


  }





  if( count == missing +  index2-index1 )
  {
     quit = 1;
     printf("We found the %d missing zeros in the following Gram block\n", missing);

    printf("index1 = %d\n", index1);
    printf("index2 = %d\n", index2);


     success = 1;
     zeroscount = zeroscount+missing;
     printf("The updated zeros count is: %d\n", zeroscount);

  }



  if( multiplier >  1024 )
  {
    quit = 1;
    success = 0;
    printf("We didn't find all the missing zeros in the next Gram block \n");
    printf("index1 = %d\n", index1);
    printf("index2 = %d\n", index2);
  }

}



}






  if(success == 0 )
  {
    printf("success = %d\n", success);

    printf("We now look in the preceding Gram Block\n");


   index1 = blocksBegin[rreindex-1];

   index2 = blocksEnd[rreindex - 1];



   y = pi*((long double) index1);

   guess = y/(logl(y)/2.0L);


for(iter=1; iter <= 15; iter++)
{




deriv = (theta(guess+delta) - theta(guess-delta))/(2.0L*delta);

guess =guess - ( theta(guess) - ((long double)index1)*pi )/deriv;


if(iter == 15)
{
   printf("Gram point guess = %.16Lf\n", guess);
}



}

printf("\n");

rs3guess = rs3(guess);


g1 = guess;
rs3g1 = rs3guess;





   y = pi*((long double) index2);

   guess = y/(logl(y)/2.0L);




for(iter=1; iter <= 15; iter++)
{




deriv = (theta(guess+delta) - theta(guess-delta))/(2.0L*delta);

guess =guess - ( theta(guess) - ((long double)index2)*pi )/deriv;


if(iter == 15)
{
   printf("Gram point guess = %.16Lf\n", guess);
}



}   

printf("\n");




rs3guess = rs3(guess);


g2 = guess;
rs3g2 = rs3guess;



if( 1 ==  1)
{
   

multiplier = 1;



quit = 0;

while( quit == 0)
{




ngriddivisor = multiplier*(1+index2-index1);

del = (g2-g1)/ ((long double) ngriddivisor    ) ;


last = rs3g1;

t = g1;

count = 0;

for( l = 1; l <= ngriddivisor; l++)
{
  t = t+del;

  new = rs3(t);
  if(new*last < 0.0L)
  {
    count = count + 1;
  }
  last = new;
}

  if( count < missing + index2-index1 )
  {

     multiplier = 2*multiplier; 

  }

  if( count > missing + index2-index1 )
  {

    printf("We have an overcount...\n");
    scanf("%d", &answer);


  }






  if( count == missing +  index2-index1 )
  {
     quit = 1;
     printf("We found the %d missing zeros in the preceding Gram block\n", missing );

    printf("index1 = %d\n", index1);
    printf("index2 = %d\n", index2);

     success = 1;
     zeroscount = zeroscount+ missing;
     printf("The updated zeros count is: %d\n", zeroscount);
  }



  if( multiplier >  1024 )
  {
    quit = 1;
    success = 0;
    printf("We didn't find all the  missing zeros in the preceding Gram block \n");
    printf("index1 = %d\n", index1);
    printf("index2 = %d\n", index2);
  }

}





if(success == 0)
{

   printf("we need to look further for the missing zeros\n\n");
}






}





}



}



/********************************************

    blocksBegin[blocknumber-1] = index1;
    blocksEnd[blocknumber-1] = index2;
    blocksZeros[blocknumber-1] = count ;


*********************************************/
   


/************************************
rosser_array[rosser_rule_exceptions-1] = rreindex;
************************************************/




/*******************************

count = 0;


last = rs3(t);

  for(k=1; k <= 97945588 ; k++)
  {
    t=t+delta; 
    new=rs3(t);
    
    if(last*new < 0.0L ) 
    {  
      count=count+1;
    }

    last = new;
  }

printf("count is %d  t is %Lf\n", count, t);

****************************/


return 0;

}

/*************************************************
? rstheta(32585736.40025939423)/Pi
%105 = 74999999.00000000000

*********************************************/






/*********************************************

[david@localhost ~]$ gcc -O2 -lm riemannsiegel3gp14b.c
[david@localhost ~]$ ./a.out 
count is 400000  t is 260877.210000

**********************************************/

/***************************************************
$ gcc -O2 -lm riemannsiegel3gp14b.c
[david@localhost ~]$ time ./a.out 
count is 799998  t is 489737.940000

real	34m49.345s
user	34m49.029s
sys	0m0.099s
*************************************/


/**************************************************************
[david@localhost ~]$ gcc -O2 -lm gramriemannsiegel3gp14b.c 
[david@localhost ~]$  ./a.out   
count is 1000000  t is 600269.9119103966204989
[david@localhost ~]$ 
**************************************************************/




/****************************************************************
[david@localhost ~]$ gcc -O2 -lm gramriemannsiegel3gp14b.c 
[david@localhost ~]$ time  ./a.out
count is 2000000  t is 1131944.7903375072875178

real	3m45.673s
user	3m45.666s
sys	0m0.002s
************************************************************/



/****************************************************************
[david@localhost ~]$ gcc -O2 -lm gramriemannsiegel3gp14b.c 
[david@localhost ~]$ 
[david@localhost ~]$ time  ./a.out
count is 12193874  t is 6000000.4859992099723058

real	51m33.528s
user	51m33.112s
sys	0m0.101s

**********************************************************/







long double theta(long double t)
{ 
   long double sum;
   
   sum = (t/2.0L)*logl(t/(2.0*pi)) - t/2.0L - pi/8.0L ;
   sum = sum + 1.0L/(48.0L*t) + 7.0L/(5760.0L*powl(t, 3.0L)) ;
   sum = sum + 31.0L/(80640.0L*powl(t, 5.0L)) ;
   sum = sum + 381.0L/(1290240.0L*powl(t, 7.0L)) ;

   return sum;

/***********************************************
rstheta =
  (X)->(X/2.0)*log(X/(2.0*Pi))-X/2.0-Pi/8.0+1/(48.0*X)+7.0/(5760.0*X^3)+31.0/(80640*X^5)+381.0/(1290240*X^7)
**************************************************/
}



long double z(long double t)
{
  long double sum;
  long double th;
  int k;
  int kmax;

  th = theta(t);


  kmax =  (int)  sqrtl(t/(2.0L*pi)) ;
  sum = (long double) 0;

  for(k=1; k<= kmax; k++)
  {
    sum = sum+ cosl( th - t*logl((long double) k) )/sqrtl((long double) k);
  }

  sum = 2.0L*sum;

  return sum;

/*************************************
Z =
  (X)->2*sum(Y=1,floor(sqrt(X/(2*Pi))),cos(rstheta(X)-X*log(Y))/sqrt(Y))
*************************************/
}


long double psi(long double w)
{
  long double y;

  y = cosl(2.0L*pi*(w*w-w-1.0L/16.0L))/cosl(2.0L*pi*w);

  return y;

/*****************************************
Psi =
  (Z)->cos(2*Pi*(Z^2-Z-1.0/16))/cos(2*Pi*Z)
*******************************************/

}


long double u(long double t)
{
  long double y;
  long double w;

  w = t/(2.0L*pi);
  y = pow(w, 0.25L);

  return y;

/*****************************
u =
  (X)->(X/(2*Pi))^(0.25)
***************************/

}


int N(long double t)
{
  int m;

  m = (int) floorl(u(t)*u(t));

  return m;

/************************
N =
  (X)->floor(u(X)^2)  
****************************/

}



long double p(long double t)
{
  long double u1;
  int m;
  long double p2;


  u1 = u(t);

  m = N(t);

  p2 = u1*u1 - ((long double) m);

  return p2;

/************************************
p =
  (X)->u(X)^2-N(X)
************************************/

}


long double d3psi(long double t)
{
  long double d3;

  d3 = -0.5L*psi(t - 2.0L*ep) + psi(t-ep) - psi(t+ep) + 0.5L*psi(t+2.0L*ep);
  d3=d3/powl(ep, 3.0L);

  return d3;

/****************************************
d3psi =
  (Z)->(-0.5*Psi(Z-epsilon*2)+Psi(Z-epsilon)-Psi(Z+epsilon)+0.5*Psi(Z+epsilon*2))/(epsilon^3)
****************************************/

}




long double d2psi(long double t)
{
  long double d2;

  d2 = ( psi(t-ep)-2.0L*psi(t) + psi(t+ep)  )/(ep*ep) ;
  
  return d2;

/*******************************
d2psi =
  (Z)->(Psi(Z-epsilon)-2*Psi(Z)+Psi(Z+epsilon))/(epsilon^2)
********************************/

}


long double d6psi(long double t)
{
  long double d6;

  d6 = psi(t-3.0L*ep) - 6.0L*psi(t-2.0L*ep) + 15.0L*psi(t-ep);
  d6 = d6 - 20.0L*psi(t) + 15.0L*psi(t+ep) - 6.0*psi(t+2.0L*ep);
  d6 = d6 + psi(t+3.0L*ep);
  d6= d6/powl(ep, 6.0L) ;

  return d6;

  
    
/*********************************************  
d6psi =
  (Z)->(Psi(Z-3*epsilon)-6*Psi(Z-2*epsilon)+15*Psi(Z-epsilon)-20*Psi(Z)+15*Psi(Z+epsilon)-6*Psi(Z+2*epsilon)+Psi(Z+3*epsilon))/(epsilon^6)
********************************************/

}


long double R3(long double t)
{
  long double r3;
  int m;
  long double sign;
  long double p2;
  long double u2;


  m = N(t)-1;

  if( (m%2) == 0 )
  {
    sign = 1.0L;
  }
  else
  {
    sign = -1.0L;
  }

  p2 = p(t); 
  u2 = u(t);


  r3 = psi(p2)/u2 - (1.0L/(96.0L*pi*pi))*d3psi(p2)/(powl(u2, 3.0L)) ;
r3 =r3+(d2psi(p2)/(64.0L*pi*pi)+d6psi(p2)/(18432.0L*pi*pi*pi*pi))/(powl(u2,5.0L));
  r3 = sign*r3; 

  return r3;



/*************************************
R3 =
  (X)->(-1)^(N(X)-1)*(Psi(p(X))/u(X)-(1/(96*Pi^2))*d3psi(p(X))/(u(X)^3)+(d2psi(p(X))/(64*(Pi^2))+d6psi(p(X))/(18432*(Pi^4)))/(u(X)^5))
*****************************************/

}


long double rs3(long double t)
{
  long double y;

  y = z(t) + R3(t);

  rs3_uses = rs3_uses+1;

  return y;

/*************************
rs3 =
  (X)->Z(X)+R3(X)
******************************/

}