/*
*  Utility & Plot for Riemann's Zeta
*
*  Authors: Rosario Turco (IT)
*  Version: 2.5  (C) 2010
*
*  Product: Tested with PARI/GP  version 2.4.2
*
*  Blog:	http://MATHBuildingBlock.blogspot.com
*
*  Site:	www.gruppoeratostene.com
*  Papers: 	Database CNR Solar
*
*
*   
*  SOFTWARE WITHOUT ANY WARRANTY WHATSOEVER.
* 
*  Attention: for Gram points you must also load lambert.txt
*  
*/

debug=0;


/* Complex Utility */

ModComplex(z)=local();{
 return(sqrt((real(z))^2 + (imag(z))^2));
}

TetaComplex(z)=local();{
 return(atan(imag(z)/real(z)));
}

/* r*exp(I*teta) */
ModTetaComplex(z)=local();{
    r = ModComplex(z);
	teta = TetaComplex(z);
    print("z= ",z, "-> ", r, "*exp(j", teta,")");
}

/* 
   Cerca gli zeri della zeta di Riemann sulla retta critica sigma=1/2.
*/

GetCriticalZero(n,pz=14)=local(k=0,i=1);{
  pz2=pz+1;
  b = n+1;
  v=vector(n);
	  
  while(i<b,
        z=-1;  
        while(z == -1,	
		        pz2 = pz2 + 1;
				z = GetZeroZetaHardy(pz,pz2);
				if(z != -1,
                   v[i]=z; 
				   i = i + 1;
				   z = 0;
 	            );	
                pz = pz2;
            ); 
	);		
	return(v);  
}

GetZeroZetaHardy(b1,b2)=local(k=0,a=-1);{
  trap( , return(-1),
          a=solve(k=b1,b2,real(RiemannSiegelZ(k)));
		  return(a);
  );

}

/*
** for s=1/2+I*t
*/
RiemannSiegelZ(t)=local(Teta=0,Zeta=0);{
  Teta = RiemannSiegelTeta(t);
  Zeta = exp(I*Teta) * zeta(1/2+I*t);
  return(Zeta);
}

RiemannSiegelTeta(t)=local(Teta=0);{
  Teta = imag(lngamma(1/4+1/2*I*t)) - ((t/2) * log(Pi));
  return(Teta);
}

/*
** for s=si+I*t
*/

RiemannSiegelZSigma(t,si)=local(Teta=0,Zeta=0);{
  Teta = RiemannSiegelTetaSigma(t,si);
  Zeta = exp(I*Teta) * zeta(si+I*t);
  return(Zeta);
}

RiemannSiegelTetaSigma(t,si)=local(Teta=0);{
  Teta = imag(lngamma(1/4+si*I*t)) - ((t/2) * log(Pi));
  return(Teta);
}

DerivRiemannSiegelSigma(s)=local(DRS=0);{
  DRS = lngamma(s) - log(2*Pi)/2;
  return(DRS);
}

Argzeta(t)=local(Arg=0);{
  z = zeta(1/2+I*t);
  Arg = atan(imag(z)/real(z));
  return(Arg);
}

S(t)=local(Si=0);{  
  Si = 1/Pi * (Argzeta(t)); 
  return(Si);
}

N(T)=local(Num=0);{  
 
  Num = (T/(2*Pi)) * log(T/(2*Pi)) - (T/(2*Pi)) + 1;
    
  return(ceil(Num));
}


/*
** The first Derivative zeta with
** s=si+I*t
**
*/
TheFirstDerivzeta(s, delta=1.0 E-20)=local(DerivZeta=0);{
  DerivZeta = (zeta(s+delta/2) - zeta(s-delta/2))/delta;
  return(DerivZeta);
}

/*
** The second Derivative zeta with
** s=si+I*t
**
*/
TheSecondDerivzeta(s, delta=1.0 E-20)=local(SecDerivZeta=0);{
  SecDerivZeta = (TheFirstDerivzeta(s+delta/2) - TheFirstDerivzeta(s-delta/2))/delta;
  return(SecDerivZeta);
}

/*
** Newton-Raphson Method for zeta 
** s=si+I*t
**
*
*/

SecantMzeta(s0,delta=1.0E-20, c=200) = local(i=0); { 
 s = s0;
 if( c > 1, 
     for(i=2,c,
	     trap( , return(s),
               a = zeta(s)/TheFirstDerivzeta(s,delta);
               s = s - a;
         );			   
     );
 );
 return(s);
}

/*
** Newton-Raphson Method for derivative of zeta
** s=si+I*t
**
*
*/

SecantMTSzeta(s0,delta=1.0E-20, c=200) = local(i=0); { 
 s = s0;
 if( c > 1, 
     for(i=2,c,
	     trap( , return(s),
               a = TheFirstDerivzeta(s,delta)/TheSecondDerivzeta(s,delta);
               s = s - a;
         );			   
     );
 );
 return(s);
}

/*
** It returns a vector of n Gram points
*/

GetGramPoints(n)=local(i=0);{
 g=vector(n);
 for(i=1,n,
    g[i]=GramPoint(i-1);
 );
 return(g);
}

/*
** It returns the gn Gram point
** 
*/

GramPoint(n)=local(gn=0);{
 gn=2*Pi*exp(1+lambertw((8*n+1)/(8*exp(1))));  
 return(gn);
}




/* 
**
** Plotting Area 
**
*/

/* Imag of zeta */
PlotImZeta(b1=0.0,b2=55.0)=local();{
  ploth(X=b1,b2,imag(zeta(1/2+X*I)),"Recursive");
 }

/* Modulus of zeta */
PlotModZeta(b1=0.0,b2=55.0)=local();{
    ploth(X=b1,b2,sqrt(real(zeta(1/2+X*I))^2 + imag(zeta(1/2+X*I))^2),"Recursive");
}

PlotModSquareZeta(b1=0.0,b2=55.0)=local();{
    ploth(X=b1,b2,real(zeta(1/2+X*I))^2 + imag(zeta(1/2+X*I))^2,"Recursive");
}
/* Real and Imag of zeta */
PlotZeta(b1=0.0,b2=55.0)=local();{
    ploth(X=b1,b2,[real(zeta(1/2+X*I)),imag(zeta(1/2+X*I))]);
}

/* Plot Real Zeta 
*
*  example: 
*  PlotRealZeta(-5,5) or 
*  PlotRealZeta(-18.5,10)  or 
*  PlotRealZeta(-15.5,0) 
*/

PlotRealZeta(b1,b2)=local();{
   ploth(t=b1,b2,zeta(t), ,100);
} 

PlotHardyZeta(b1,b2)=local(t=0);{
   ploth(t=b1,b2,real(RiemannSiegelZ(t)));
}