# This file contains some procedures to check the calculations of the paper
# "Estimates of li(theta(x))-pi(x) and the Riemann hypothesis"

Digits:=20;

lambda:=evalf(gamma)+2-log(4*evalf(Pi));

L:=proc(N,tt) local t,k,L;
t:=evalf(tt);
Li(t)-sum((k-1)!*t/log(t)^k,k=1..N);
end;

theta:=proc(x) local p,pf,logp,S; # Chebichev function
S:=0.; p:=2; pf:=evalf(p); logp:=log(pf);
while p <= x do
   S:=S+logp;
   p:=nextprime(p); pf:=evalf(p); logp:=log(pf);
od;
S;
end:

B:=proc(yy) local y; # calculation of  B(y)=pi(y)-theta(y)/log(y), cf. § 3.2
y:=evalf(yy);
numtheory[pi](floor(y))-theta(y)/log(y);
end;

#######################################################
# VALUES OF A(p)
########################################################

A:=proc(x) Li(theta(x))-numtheory[pi](x) end; # calculates A(x)

valA:=proc(pimax) local p,i;
# print the values of A(p) for the first pimax primes 
# and for p= 401,409, 3643,33647.
print(); print("VALUES of A(p)"); print();
for i from 1 to pimax do p[i]:= ithprime(i); od;
for i from 1 by 2 to pimax do
   print("p=",p[i],evalf(A(p[i]),12),"     p=",p[i+1],evalf(A(p[i+1]),12)); 
od;
print();
print("p=",401,evalf(A(401),12), "     p=",409,evalf(A(409),12));
print("p=",3643,evalf(A(3643),12), "     p=",33647,evalf(A(33647),12)); 
end;

#######################################################
# LEMMA 3.3
########################################################


lem33:=proc() local y0,y1,y2,y1f,y2f,pif;
y0:=8.3; y1:=599; y2:=2903; y1f:=evalf(y1); y2f:=evalf(y2);pif:=evalf(Pi);
print("NUMERICAL VALUES in LEMMA 3.3"); print();
print("y0=",y0, "B(y0)-L_1(y0)=", evalf(B(y0)-L(1,y0),12));
print("y1=",y1, "B(y1)-L_1(y1)-sqrt(y1)/(4*pi)=",
       evalf(B(y1f)-L(1,y1f)-sqrt(y1f)/4/pif,12));
print("y2=",y2, "B(y2)-L_1(y2)+sqrt(y2)/(4*pi)=",
       evalf(B(y2f)-L(1,y2f)+sqrt(y2f)/4/pif,12));
print();
end;

###########################################################
# COROLLARY 3.1
############################################################

coro31:=proc() local x0,U,pif;
print("NUMERICAL VALUE in COROLLARY 3.1"); print();
x0:=10.^8; pif:=evalf(Pi);
U:=7.993-log(x0)^3/8/pif/x0^(1/4)-18/10000*log(x0)^5/sqrt(x0);
print("7.993-log(10^8)^3/8/pi/(10^8)^(1/4)-18/10000*log(10^8)^5/sqrt(10^8)=",
          evalf(U,12)); print();
end;


###########################################################
# PROPOSITION 3.4
############################################################

F1:=proc(tt) local t; t:=evalf(tt);L(1,t)/(t/log(t)^2); end;

F2:=proc(tt) local t;t:=evalf(tt); L(2,t)/(t/log(t)^3); end;

F1til:=proc(t); # cf. (3.15)
if t > 95 then F1(t) else 1.785; fi;
end;

F2til:=proc(t); # cf. (3.15)
if t > 381 then F2(t) else 4.05; fi;
end;

Q:=proc(ka1,xx) local x,k,sqx,logx,T1,T2,T3,T4,T5,T6,T7,lisT4,T4som,T;
# computes Q defined by (3.24)
print("CALCULATION of Q(kappa_1,x) by FORMULA (3.21)");print();
x:=evalf(xx); logx:=log(x); sqx:=sqrt(x);
print("ka1=",ka1,"x=",xx,"log x=",evalf(logx,12),"rac x=",sqx);
T1:=4*F2til(sqx);
print("T1=",evalf(T1,12));
T2:=2/300;
print("T2=",T2);
T3:=3.05*logx^3/sqx;
print("T3=",evalf(T3,12));
for k from 3 to ka1 do
   T4[k]:=k*F1til(x^(1/k))*logx/x^(1/2-1/k);
   print("k=",k,"T4[k]=",evalf(T4[k],12));
od;
lisT4:=[T4[kk] $ kk=3..ka1]; T4som:=add(i,i=lisT4);
print("T4=",evalf(T4som,12));
T5:=7.23*ka1^3/x^(1/2-1/ka1);
print("T5=",evalf(T5,12));
T6:=0.94/logx^2;
print("T6=",evalf(T6,12));
T7:=9*logx^5/10000/sqx;
print("T7=",evalf(T7,12));
T:=T1+T2+T3+T4som+T5+T6+T7;
print("T=Q(kappa_1,x)=",evalf(T,12));
T;
end;

Qminka1:=proc(xx)
local x,logx,sqx,ka1,ka1m,k,minQ,T,T1,T2,T3,T4,T5,T6,T7,lisT4,T4som;
# This procedure calculates the value of kappa_1 that
# minimises Q(kappa_1,x). Choose "Digits" large if xx is large.
print("DETERMINATION of kappa MINIMIZING  Q(kappa_1,x) for x=",xx);print();
x:=evalf(xx); logx:=log(x); sqx:=sqrt(x); minQ:=infinity;
for ka1 from 3 to log(x)/log(10.4) do # ka1 is < log(x)/log(10.4)
   T1:=4*F2til(sqx);
   T2:=2/300;
   T3:=3.05*logx^3/sqx;
   for k from 3 to ka1 do
      T4[k]:=k*F1til(x^(1/k))*logx/x^(1/2-1/k);
   od;
   lisT4:=[T4[kk] $ kk=3..ka1]; T4som:=add(i,i=lisT4);
   T5:=7.23*ka1^3/x^(1/2-1/ka1);
   T6:=0.94/logx^2;
   T7:=9*logx^5/10000/sqx;
   T:=T1+T2+T3+T4som+T5+T6+T7;
   if T < minQ then minQ:=T; ka1m:=ka1; fi;
   print("kappa_1=",ka1,"T=Q(kappa_1,x)=",evalf(T,12));
od;
print();
print("x=",xx,"log10(x)=",evalf(log10(x),3),
      "minimum of Q=",evalf(minQ,12),"for ka1=",ka1m); print();
end;


####################################################
# PROPOSITION 3.6 : COMPUTATION of A(x) FOR x SMALL
####################################################
#####################################################
#  PROPOSITION 3.6 : SOLUTION of (i)
#####################################################

prop36i:=proc()
print("SOLUTION of PROPOSITION 3.6 (i)"); priny();
print("A(11)=",evalf(A(11),12),"A(7)=",evalf(A(7),12)); print();
end;

##########################################################
#  PROPOSITION 3.6 : SOLUTION of (ii), (iii) and (v) FOR 409 <= x <=10^8
##########################################################

prop36ii:=proc(pmax,pas,const) 
local A,p,pp,pf,ppf,logp,logpp,losqp,losqpp,pi,theta,ti,
ppas, # intermediary results are printed 
# when p exceeds ppas that takes the values multiple of pas
CM, # defined by A(p)=sqrt(p)/log^2(p)*(2+lambda+CM/log(p))
maxCM,maxCMp,ACM,
# maxCM=max(C(p), 409<= p <= pmax, maxCM=C(maxCMp,ACM=A(maxCMp)
p0,A0, # p0 is the largest prime s.t. C(p0)>=const, 
# A0=A(p0), pp0 is the prime following p0
cm, # is c tilda(p) defined by A(p)=sqrt(pp)/log^2(pp)*(2-lambda+cm/log(pp))
mincm, mincmp, Acm, 
# mincm is the min of cm, mincm is attained in mincmp,Acm=A(mincmp)
pp0,pp0f,logpp0,C0,Cpp0,t0;
global lambda;
# calculation of A(p) and C(p) and ctilda(p) for p <= 10^8
ti:=time();
print("SOLUTION of PROPOSITION 3.6 (ii), (iii) and (v)");
print("pmax=",pmax,"log10(pmax)=",evalf(log10(pmax),2),
      "pas=",pas,"const=",const);print();
ppas:=pas;
maxCM:=-infinity; p0:=-2; mincm:=infinity; #initialisation
p:=2; pi:=0; theta:=0.;
pf:=2.;logp:=log(pf); losqp:=logp^2/sqrt(2.);
while p <=pmax do
   pp:=nextprime(p); ppf:=evalf(pp); logpp:=log(ppf); losqpp:=logpp^2/sqrt(ppf);
   pi:=pi+1; theta:=theta+logp;
   A:=Li(theta)-pi;
   if p >= 409 then
      CM:=(A*losqp-2-lambda)*logp;
      if CM > maxCM then maxCM:=CM; maxCMp:=p; ACM:=A; fi;
      if CM > const then p0:=p; A0:=A; C0:=CM; fi;
      cm:=(A*losqpp-2+lambda)*logpp;
      if cm < mincm then mincm:=cm; mincmp:=p; Acm:=A; fi;
   fi;
   if pp > ppas then
      ppas:=ppas+pas;
      print("p=",p,"pp=",pp,"A(p)=",evalf(A,12));
      print("C(p)=",evalf(CM,12),"ctilda(p)=",evalf(cm,12),"time=",time()-ti);
      print();
   fi;
   p:=pp; pf:=ppf; logp:=logpp; losqp:=losqpp;
od;
print(); print("SOLUTION of PROPOSITION 3.6 (ii)"); print();
print("max C(p)=",evalf(maxCM,12),"for p=",maxCMp, "A(p)=",evalf(ACM,12)); 
print(); print("SOLUTION of PROPOSITION 3.6 (iii)"); print();
pp0:=nextprime(p0); pp0f:=evalf(pp0); logpp0:=log(pp0f);
Cpp0:=logpp0*(A0*logpp0^2/sqrt(pp0f)-2-lambda);
if Cpp0 < const then
   t0:=fsolve(log(t)*(A0*log(t)^2/sqrt(t)-2-lambda)=const,t=p0..pp0);
fi;
print("p0=",p0,"p0+=",pp0,"A(p0)=",evalf(A0,12),"C(p0)=",evalf(C0,12));
print("lim_{x->p0+} C(x)=",evalf(Cpp0,12),"t=",evalf(t0,12)); print();
print(); print("SOLUTION of PROPOSITION 3.6 (v) for 409 <= p < 10^8"); print();
print("Minimum of ctilda(p) for 409 <= p < 10^8=", 
       evalf(mincm,12),"for p=",mincmp);
print("Total time=",time()-ti); print();
end;

##########################################################
#  PROPOSITION 3.6 : SOLUTION of (v) FOR 83 <= x < 409
##########################################################

prop36v:=proc(pmax) 
local p,pp,pf,ppf,logp,logpp,pi,theta,f,maxf,t0,t,mindif,dif,mindifp;
global A,lambda;
# Choose pmax >= 401.
f:=proc(tt) local t; t:=evalf(tt); sqrt(t)/log(t)^2*(2-lambda+5.12/log(t)); end;
mindif:=infinity;
p:=2; pi:=0; theta:=0.;
pf:=2.;logp:=log(pf); 
print("SOLUTION of PROPOSITION 3.6 (v) for 83 <= p < 409"); print();
while p <=pmax do
   pi:=pi+1; theta:=theta+logp;
   A[p]:=Li(theta)-pi;
   pp:=nextprime(p); ppf:=evalf(pp); logpp:=log(ppf);
   if p = 83 then
      t0:=fsolve(f(t)=A[p],t=83..89);
      print("A(83)=",evalf(A[p],12),"f(83)=", evalf(f(83),12),
            "f(89)=",evalf(f(89),12));
      print("f(t0)=A(83) for t0=",evalf(t0,12)); print();
   fi;
   if 89 <= p and p < 409 then
      dif:=A[p]-max(f(p),f(pp));
      if dif < mindif then mindif:=dif; mindifp:=p; fi
 # print("p=",p,"Ap-fp=",evalf(A[p]-f(p),12),"Ap-fp+=",evalf(A[p]-f(pp),12));
   fi;
   p:=pp; pf:=ppf; logp:=logpp;
od;
print("Minimum for 89 <= p < 409 of A(p)-max(f(p),f(p+))=",evalf(mindif,12),
        "for p=",mindifp); print();
end;


##############################################################
#  PROPOSITION 3.6 : SOLUTION of (iv) FOR x <=10^4 and (vi) for 37 <= x <89
##############################################################

prop36iv:=proc(pmax) local p,logp,pp,pi,theta,phi,
M,MaxM,Mp; # MaxM is the max of M=A(p)*log^2(p)/sqrt(p). It is attained in Mp
global A,lambda;
phi:=proc(tt) local t; t:=evalf(tt); log(t)^2/sqrt(t); end;
MaxM:=0;
p:=2; pi:=0; theta:=0.; logp:=log(2.);
while p <=pmax do
   pi:=pi+1; theta:=theta+logp;
   A[p]:=Li(theta)-pi;
   M[p]:=A[p]*phi(p);
   if  M[p] > MaxM then MaxM:=M[p]; Mp:=p; fi;
   p:=nextprime(p); logp:=log(evalf(p));
od;
print("SOLUTION of PROPOSITION 3.6 (iv)"); print();
print("max for 59 <= p < 10000 of A(p)*log^2(p)/sqrt(p)=",evalf(MaxM,12),
        "for p =",Mp); print();
print("SOLUTION of PROPOSITION 3.6 (vi)"); print();
print("2-lambda=",evalf(2-lambda,12));
p:=11; pp:=13;
while p <= 83 do
   print("p=",p,"pp=",pp,"A(p)phi(p)=",evalf(A[p]*phi(p),12),
         "A(p)phi(p+)=",evalf(A[p]*phi(pp),12));
   p:=pp; pp:=nextprime(pp);
od;
print();
end;