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

long double rs3(long double);
long double Z(long double);
long double R3taylor(long double);
long double theta(long double);
long double u(long double);
long double p(long double);
int N(long double);
long double grampoint(long);
long double C0(long double);
long double C1(long double);
long double C2(long double);
long double C3(long double);
long double C4(long double);
long double rs3taylor(long double);

long blocksBegin[100000000];
long blocksEnd  [100000000];
long blocksZeros[100000000];

long INITIAL_GRAM_POINT = 290000000;
long LAST_GRAM_POINT   =  300000000;


long double ep = 0.00064L ;

long double pi = 3.14159265358979323846L ;
long double tmin = 7.0L;
long double tmax = 10000000000000.0L;


int rs3_uses = 0;



int main(void)
{

long double t;
long double last;
long double new;
long count;
long double delta = 0.005L;
int k;
long j;
long l;
int found;
long index1;
long index2;
long double g1;
long double g2;
long zeroscount;
long multiplier;
int quit;
long ngriddivisor;
long double del;
long double rs3g1;
long double rs3g2;
int L;
int blockcount[100];
int gramitv; 
long newmultiplier;
long nmulti;
int blocknumber = 0;
int rosser_rule_exceptions = 0;
int rreindex;
int rosser_array[1000];
int kk;
int success;
long missing;
int answer;
long double gp;
long double rs3g; 
long double diff;



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



zeroscount = 0;

found = 0;

printf("Initial Gram point = %ld\n", INITIAL_GRAM_POINT);
printf("Last Gram point = %ld\n\n", LAST_GRAM_POINT);


printf("sizeof int = %d\n", sizeof(int));
printf("sizeof long = %d\n", sizeof(long));


j = INITIAL_GRAM_POINT;
gp = grampoint(j);
rs3g = rs3(gp);

if( rs3g > 0.0L && (1 == ((j+2)%2)) )
{
  printf("bad gram point at %ld\n", j);
  scanf("%d", &answer);
}

  if( rs3g < 0.0L && (0 == ((j+2)%2)) )
  {
     printf("bad gram point at %ld\n", j);
     scanf("%d", &answer);
  }





j = LAST_GRAM_POINT;
gp = grampoint(j);
rs3g = rs3(gp);

if( rs3g > 0.0L && (1 == ((j+2)%2)) )
{
  printf("bad gram point at %ld\n", j);
  scanf("%d", &answer);
}

  if( rs3g < 0.0L && (0 == ((j+2)%2)) )
  {
     printf("bad gram point at %ld\n", j);
     scanf("%d", &answer);
  }








for(j = INITIAL_GRAM_POINT ; j<= LAST_GRAM_POINT ; j++)
{ 



gp = grampoint(j);



rs3g = rs3(gp);


  if( rs3g > 0.0L && (1 == ((j+2)%2))  && ((j==INITIAL_GRAM_POINT) || (j==LAST_GRAM_POINT)) )
  {
    printf("bad gram point at %ld\n", j);
  }




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




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





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



    if(found == 0)
    { 
       found = 1;
       L = 0;
       index1 = j;
       g1 = gp;
       rs3g1 = rs3g;
    }


      
  }



  if( rs3g < 0.0L && (0 == ((j+2)%2))&&((j==INITIAL_GRAM_POINT) || (j==LAST_GRAM_POINT))   )
  {
    printf("bad gram point at %ld\n", j);
  }

  if( rs3g < 0.0L && (0 == ((j+2)%2))      )
  {
   L = L+1;
  }



 if( rs3g < 0.0L && (1 == ((j+2)%2)) )
  {
   // good gram point at j

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


    if(found == 0)
    { 
       found = 1;
       L = 0;
       index1 = j;
       g1 = gp;
       rs3g1 = rs3g;
    }



  }



if( found == 2)
{
  /*** scan the Gram block for zeta zeros  ***/



  blockcount[L] ++ ;

  blocknumber++;





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





      if( 0 == (j%100000))
      {
         printf("good gram point at %ld\n", j);
         printf("The cumulative zeros count is: %ld\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");
    printf("index1 = %ld\n", index1);
    printf("index2 = %ld\n", index2);
    printf("adding %ld zeros to the count\n\n", index2-index1);
    quit = 1;


    zeroscount = zeroscount +   (index2-index1);

    blocksBegin[blocknumber-1] = index1;

    blocksEnd[blocknumber-1] = index2;

    blocksZeros[blocknumber-1] =  index2- index1;

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




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

     zeroscount = zeroscount + count;   






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






blocksBegin[blocknumber-1] = index1;

blocksEnd[blocknumber-1] = index2;

blocksZeros[blocknumber-1] =  index2- index1;


  }

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

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

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

     zeroscount = zeroscount + count;
// new  


rosser_rule_exceptions ++;

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

rreindex = blocknumber-1;


rosser_array[rosser_rule_exceptions-1] = rreindex;




  }

}



}



found = 1;

L = 0;


g1 = g2;

index1 = index2;

rs3g1 = rs3g2;


}
  /*** done scanning  the Gram block for zeta zeros  ***/

}
/***** done big loop for all Gram points in range ***/



printf("\n\ncount is %ld zeros\n",zeroscount);
printf("in [%.4Lf, %.4Lf]\n", grampoint(INITIAL_GRAM_POINT), grampoint(LAST_GRAM_POINT) );




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++)
{

printf("\n");



   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 have now %ld missing zeros in the block from gram point %ld to gram point %ld\n", missing, index1, index2);




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


   index1 = blocksBegin[rreindex+1];

   index2 = blocksEnd[rreindex + 1];

   g1 = grampoint(index1);


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


printf("\n");

rs3g1 = rs3(g1);


g2 = grampoint(index2);


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

printf("\n");

rs3g2 = rs3(g2);





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 %ld missing zeros in the following Gram block\n", missing);

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


     success = 1;
     zeroscount = zeroscount+missing;
     printf("The updated zeros count is: %ld\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 = %ld\n", index1);
    printf("index2 = %ld\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];
   g1= grampoint(index1);
   printf("Gram point g1 = %.16Lf\n", g1);
   printf("\n");
   rs3g1 = rs3(g1);
   g2= grampoint(index2);
   printf("g2 = %.16Lf\n", g2);
   printf("\n");
   rs3g2 = rs3(g2);
   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 %ld missing zeros in the preceding Gram block\n", missing );

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

     success = 1;
     zeroscount = zeroscount+ missing;
     printf("The updated zeros count is: %ld\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 = %ld\n", index1);
    printf("index2 = %ld\n", index2);
  }

}





if(success == 0)
{

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



}

}
/*********** end of loop: 
for(kk = 1; kk <= rosser_rule_exceptions; kk++) **********/

return 0;

}

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

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

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 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 R3taylor(long double t)
{
  long double r3;
  int m;
  long double sign;
  long double p2;
  long double u2;
  long double z;


  m = N(t)-1;

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

  p2 = p(t); 
  u2 = u(t);
  z = 2.0L*p2-1.0L;

  r3 = C0(z)/u2;

  r3 = r3 + C1(z)/powl(u2, 3.0L);

  r3 = r3 + C2(z)/powl(u2, 5.0L);

  r3 = r3 + C3(z)/powl(u2, 7.0L);

  r3 = r3 + C4(z)/powl(u2, 9.0L);


  r3 = sign*r3; 

  return r3;




}





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

  y = Z(t) + R3taylor(t);

  rs3_uses = rs3_uses+1;

  return y;


}


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

  y = Z(t) + R3taylor(t);


  return y;


}





long double C0(long double z)
{
long double w;

w = 0.3826834323650897717285L;

w = w + 0.4372404680775204493603L*powl(z, 2.0L);

w = w+ 0.1323765754803435233240L*powl(z, 4.0L);

w = w -0.01360502604767418865498L*powl(z, 6.0L);

w = w -0.01356762197010358088792L*powl(z, 8.0L);

w = w -0.001623725323144465282855L*powl(z, 10.0L);

w = w +   0.0002970535373337969078313L*powl(z, 12.0L);

w = w+   0.00007943300879521469588016L*powl(z, 14.0L);

w = w +   0.0000004655612461450450503706L*powl(z, 16.0L);

w = w -0.000001432725163095510575408L*powl(z, 18.0L);

w = w -0.0000001035484711231294607501L*powl(z, 20.0L);

w = w +   0.00000001235792708386173805613L*powl(z, 22.0L);

w = w + 0.000000001788108385795490498567L*powl(z, 24.0L);

w = w -0.00000000003391414389927035906941L*powl(z, 26.0L);

w = w - 0.00000000001632663390256590510137L*powl(z, 28.0L);

w = w - 0.0000000000003785109318541220382855L*powl(z, 30.0L);

w = w +  0.00000000000009327423259201724845662L*powl(z, 32.0L);

w = w +   0.000000000000005221843015978136855314L*powl(z, 34.0L);

w = w -0.0000000000000003350673072744263789515L*powl(z, 36.0L);

w = w -0.00000000000000003412426522811726494081L*powl(z, 38.0L);

w = w + 0.0000000000000000005751203341432399160340L*powl(z, 40.0L);

w = w + 0.0000000000000000001489530136321150545476L*powl(z, 42.0L);

w = w+ 0.000000000000000000001256537271702141685330L*powl(z, 44.0L);

w = w -0.0000000000000000000004721295250143425668954L*powl(z, 46.0L);

w = w -0.00000000000000000000001326906936303961999277L*powl(z, 48.0L);

return w;

}


//  ? for(X=0, 48, print(X,"   ",polcoeff(c1, X)))

long double C1(long double z)
{
long double w;

w = -0.02682510262837534702999L*z;

w = w + 0.01378477342635185304987L*powl(z,3.0L);

w = w +  0.03849125048223508222874L*powl(z, 5.0L);

w = w +  0.009871066299062076472012L*powl(z, 7.0L);

w = w -0.003310759760858404332909L*powl(z, 9.0L);

w = w  -0.001464780857795415082498L*powl(z,11.0L);

w = w  -0.00001320794062487696367516L*powl(z,13.0L);

w = w + 0.00005922748701847141323223L*powl(z,15.0L);

w = w+ 0.000005980242585373448587711L*powl(z,17.0L);

w = w -0.0000009641322456169826352673L*powl(z,19.0L);

w = w -0.0000001833473372271441176002L*powl(z, 21.0L);

w = w + 0.000000004467087562717833599561L*powl(z, 23.0L);

w = w + 0.000000002709635082177274321693L*powl(z, 25.0L);

w = w + 0.00000000007785288654315851046295L*powl(z,27.0L);

w = w -0.00000000002343762601089368853248L*powl(z,29.0L);

w = w  -0.000000000001583017278998752164216L*powl(z,31.0L);

w = w + 0.0000000000001211994157372379124665L*powl(z, 33.0L);

w = w +  0.00000000000001458378116110830701758L*powl(z,35.0L);

w = w  -0.0000000000000002878630525813191750456L*powl(z,37.0L);

w = w -0.00000000000000008662862902123724122528L*powl(z,39.0L);

w = w -0.0000000000000000008430722727137041271560L*powl(z,41.0L);

w = w+  0.0000000000000000003630807223097346200173L*powl(z,43.0L);

w = w +  0.00000000000000000001162669821283829671945L*powl(z,45.0L);

w = w -0.000000000000000000001097548671152753181440L*powl(z,47.0L);

return w;
}



// ? for(X=0, 48, print(X,"   ",polcoeff(c2, X)))

long double C2(long double z)
{
long double w;

w =  0.005188542830293168493785L;

w = w+ 0.0003094658388063474603346L*powl(z,2.0L);

w = w -0.01133594107822937338218L*powl(z,4.0L);

w = w + 0.002233045741958144772057L*powl(z,6.0L);

w = w + 0.005196637408862330205117L*powl(z,8.0L);

w = w+  0.0003439914407620833669466L*powl(z,10.0L);

w = w  -0.0005910648427470582821732L*powl(z,12.0L);

w = w  -0.0001022997254793585745443L*powl(z,14.0L);

w = w +  0.00002088839221699275540807L*powl(z,16.0L);

w = w+ 0.000005927665493096535957892L*powl(z,18.0L);

w = w -0.0000001642383836243627597769L*powl(z,20.0L);

w = w -0.0000001516119970094068286173L*powl(z,22.0L);

w = w -0.000000005907803698206667962923L*powl(z,24.0L);

w = w + 0.000000002091151485947818897775L*powl(z,26.0L);

w = w + 0.0000000001781564958329235105380L*powl(z,28.0L);

w = w -0.00000000001616407245535383075286L*powl(z,30.0L);

w = w -0.000000000002380696249666761570721L*powl(z,32.0L);

w = w+ 0.00000000000005398265295542594918182L*powl(z,34.0L);

w = w+ 0.00000000000001975014219696951527331L*powl(z,36.0L);

w = w + 0.0000000000000002333286873288263483105L*powl(z,38.0L);

w = w -0.0000000000000001118751761004808020820L*powl(z,40.0L);

w = w -0.000000000000000004164009488883767188513L*powl(z,42.0L);

w = w + 0.0000000000000000004446081109291883028286L*powl(z,44.0L);

w = w +  0.00000000000000000002854611478363714422952L*powl(z,46.0L);

w = w -0.000000000000000000001191323143003791035106L*powl(z,48.0L);

return w;

}




// ? for(X=0, 48, print(X,"   ",polcoeff(c3, X)))

long double C3(long double z)
{
long double w;

w = -0.001339716090719456904270L*z;

w = w+ 0.003744215136379393704664L*powl(z,3.0L);

w = w -0.001330317891932146812032L*powl(z,5.0L);

w = w -0.002265466076547178711476L*powl(z,7.0L);

w = w+ 0.0009548499998506730415112L*powl(z,9.0L);

w = w + 0.0006010038458963603912076L*powl(z,11.0L);

w = w -0.0001012885828677662195334L*powl(z,13.0L);

w = w -0.00006865733449299825642457L*powl(z,15.0L);

w = w+ 0.0000005985366791538598159306L*powl(z,17.0L);

w = w + 0.000003331659851239947129044L*powl(z,19.0L);

w = w+ 0.0000002191928910243508105718L*powl(z,21.0L);

w = w -0.00000007890884245681494410555L*powl(z,23.0L);

w = w -0.000000009414685081295262151652L*powl(z,25.0L);

w = w+ 0.0000000009570116210883480301881L*powl(z,27.0L);

w = w+ 0.0000000001876313745347066279681L*powl(z,29.0L);

w = w -0.000000000004437837679323399327465L*powl(z,31.0L);

w = w -0.000000000002242673850561735324841L*powl(z,33.0L);

w = w -0.00000000000003627686865735243689408L*powl(z,35.0L);

w = w+  0.00000000000001763980955082158160783L*powl(z,37.0L);

w = w + 0.0000000000000007960765246786777757313L*powl(z,39.0L);

w = w -0.00000000000000009419651490589690762535L*powl(z,41.0L);

w = w -0.000000000000000007133103854569657743497L*powl(z,43.0L);

w = w+ 0.0000000000000000003289910584554629007871L*powl(z,45.0L);

w = w + 0.00000000000000000004180730374898726358829L*powl(z,47.0L);

return w;

}



//  ? for(X=0, 48, print(X,"   ",polcoeff(c4, X)))

long double C4(long double z)
{
long double w;


w =  0.0004648338936176338185363L;

w = w -0.001005660736534047075978L*powl(z,2.0L);

w = w +  0.0002404485657372579302245L*powl(z,4.0L);

w = w + 0.001028308614970232187826L*powl(z,6.0L);

w = w -0.0007657861071755644186600L*powl(z,8.0L);

w = w -0.0002036528680308481762148L*powl(z,10.0L);

w = w+ 0.0002321229049106872789514L*powl(z,12.0L);

w = w + 0.00003260214424386519760774L*powl(z,14.0L);

w = w -0.00002557906251794952514025L*powl(z,16.0L);

w = w -0.000004107464438915744753982L*powl(z,18.0L);

w = w + 0.000001178111364037129388130L*powl(z,20.0L);

w = w+  0.0000002445656142248457854232L*powl(z,22.0L);

w = w -0.00000002391582476734432243033L*powl(z, 24.0L);

w = w -0.000000007505214207035755288539L*powl(z,26.0L);

w = w+ 0.0000000001331227941625842819291L*powl(z, 28.0L);

w = w+ 0.0000000001344062675422561971870L*powl(z, 30.0L);

w = w + 0.000000000003513770042430485928693L*powl(z,32.0L);

w = w -0.000000000001519154453370391933574L*powl(z, 34.0L);

w = w -0.00000000000008915417681447087305522L*powl(z,36.0L);

w = w+ 0.00000000000001119589116522853577137L*powl(z,38.0L);

w = w+ 0.000000000000001051601332991481483697L*powl(z,40.0L);

w = w -0.00000000000000005178655273646692084491L*powl(z,42.0L);

w = w -0.000000000000000008065874861917114521441L*powl(z,44.0L);

w = w +  0.0000000000000000001060820453021358190156L*powl(z,46.0L);

w = w+ 0.00000000000000000004433680672100193810607L*powl(z,48.0L);

return w;
}



long double grampoint(long k)
{
int iter;
long double left;
long double right;
long double tmid;
long double y;


left = tmin;
right = tmax;



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

  tmid = (left+right)/2.0L;

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

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

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

} 

return tmid;

}