#COMPUTATION FOR THE ARTICLE: hdenHR.tex

############################################
# constants
#############################################

Digits:=20;

const:=proc() local x1,x0, n0,L0,la0,nu0,c;
x0:=nextprime(10^10);
n0:=2220822442581729257;
L0:=log(evalf(n0));
la0:=log(L0);
nu0:=la0/L0;
c:=0.046117644421509023827;
print("x0=",x0);
print("n0=",n0);
print("L0=",L0);
print("lambda_0=",la0);
print("nu0=",nu0);
print("c=",c);
end;

############################################
# Lemma 2.1
#############################################

lemma2_1:=proc() local lis,noli,p,q,b,c,i,lisf,maxi;
lis:=[]; p:=2;
while p <= 11 do
   q:=p^2;
   while q <=121 do
      lis:=[op(lis),[q,p]];
      q:=q*p;
   od;
   p:=nextprime(p);
od;
lis:=sort(lis,proc(u,v) u[1] < v[1] end);
noli:=nops(lis);
# here "lis" contains the ordered list of the powers of primes with exponent > 1
print("lis=",lis);
b[0]:=1;lisf:=[2];
for i from 1 to noli-1 do
   b[i]:=b[i-1]*lis[i][2];#log(b[i]) is equal to psi(x)-theta(x) for x=lis[i][1]
   c[i]:=sqrt(lis[i+1][1])-log(b[i]);
   lisf:=[op(lisf),c[i]];
od;
# For x < 4, psi(x) is equal to theta(x) and the maximum of
# sqrt(x)+theta(x)-psi(x) is 2.
#  On the interval [lis[i][1],lis[i+1][1]), psi(x)-theta(x)
# is constant and equal to log(b[i] and sqrt(x) is maximal
# for x=lis[i+1]-epsilon. Therefore, the maximum of
# sqrt(x)+theta(x)-psi(x) is c[i].
print("lisf=",lisf); print("lisf=",evalf(lisf,5));
maxi:=max(op(lisf));
print("maxi=",maxi,"=",evalf(maxi,5));
end;

########################################################
# § 2.2 Logarithmic integral
########################################################

liVal:=proc() local lis,lis1;
lis:=[2.,exp(2.)];
print("Values of li = ",lis,map(u->Li(u),lis));
lis1:=[0.,1.];
print("values of li^(-1) = ",lis1,map(u->Lim1Newtonx0(u,u*log(u),0),lis1));
end;

formula2_7:=proc(Nmax) local N,f,k,f1; global s,t;
for N from 1 to Nmax do
   f:=1/(N-1)!*(s^(N-1)*Li(t^s)-sum((k-1)!*s^(N-1-k)*t^s/(log(t))^k,k=1..N-1));
   f1:=simplify(diff(f,t),assume=positive);
   print("  N=",N," f= ",f," f ' = ",f1);
od;
end;

lemma2_2:=proc() local f,f1,t0,mini,maxi,f0,t1; global s,t;
f:=Li(t)-t/log(t); f1:=diff(f,t);
t0:=fsolve(f=0,t=2..5);
print("f=",f,"f ' =",f1,"t0=",t0);
#######
f:=t-Li(t); f1:=diff(f,t); mini:=evalf(subs(t=exp(1.),f));
print("f=",f,"f ' =",f1,"min=",mini);
#######
f:=Li(t)-1.49*t/log(t); f1:=normal(diff(f,t));
maxi:=evalf(subs(t=exp(1.49/0.49),f));
print("f=",f,"f ' =",f1,"max=",maxi);
#######
f:=Li(t)-t/log(t)-1.101*t/log(t)^2; f1:=normal(diff(f,t));
t1:=fsolve(f1,t=10^9..10^10); f0:=evalf(subs(t=10.^10,f));
print("f=",f,"f ' =",f1); print("t1=",t1,"f(10^10)=",f0);
end;

lemma2_3:=proc() local f,f1,t0,u; global s,t;
f:=Li(t)-t/log(t)-t/log(t)^2-2*t/log(t)^3-u*t/log(t)^4;
f1:=normal(diff(f,t));
print("f=",f,"f ' =",f1);
t0:=fsolve(subs(u=6,f)=0,t=10..100);
print("u=6","t0=",t0);
print("u=7","exp(28)=",exp(28.),"f(4.96 10^12)=",
         evalf(subs(u=7,t=4.96*10^12,f)));
print("u=40/3","exp(80/11)=",exp(80./11),"f(exp(80/11))=",
         evalf(subs(u=40/3,t=exp(80./11),f)));
end;

lemma2_4:=proc() local f,f328,num,y,t0; global t,L;
f:=Li(t*(log(t)+log(log(t))))-t;
f328:=evalf(subs(t=3.28,f));
print("f=",f,"f(3.28)=",f328);
f:=Li(t*log(t))-t;
print("f=",f,"f(41)=",evalf(subs(t=41.,f)));
num:=(L-1)*(log(L)-2)-1;
print("num=",num,"subs(L=log(4678),num)=",evalf(subs(L=log(4678.),num)));
f:=t*(log(t)+log(log(t))-1);
t0:=fsolve(f=1,t=3..4);
print("f=",f,"t0=",t0);
y:=t-Li(f);
print("y(12218)=",evalf(subs(t=12218.,y)));
print();
end;

lemma2_5:=proc() local f11,f12,f2,f21,f22,B,B103,Li103; global t,u,L;
print("The derivative of li^(-1) (t) is log u with u=li^(-1) (t)"); print();
f11:=log(u)/2/sqrt(u);
f12:=normal(diff(f11,u))*log(u);
print("df1/dt=",f11," d2f1/dt2=",f12);
f2:=(t*log(t))^(1/4);
f21:=normal(diff(f2,t));
f22:=normal(diff(f21,t));
B:=1/(4*u^(1/4))*(1-log(L)/L)^(-7/4)*(3+8/L+19/L/(L-2));
print("d2f2/dt2=",f22," B=",B);
B103:=evalf(subs(u=103.,L=log(103.),B));
print("Li(103)=",Li(103.),"subs(u=103,L=log(103),B)=",B103);
end;

Lim1Newtonx0:=proc(xf,x0,imp) local x,yanc,ynou,digits,dif,eps;
# computes the inverse function of the logarithmic integral
# by the Newton method.
x:=evalf(xf);
if imp > 5 then print("calcul de Li(x), x=",x); print(); fi;
eps:=10.^(-Digits-2);
digits:=Digits; Digits:=digits+5;
if x <= 0. then
   yanc:=1.+0.25*exp(x);
   dif:=infinity;
   while abs(dif) > eps do
      dif:=(x-Li(yanc))*log(yanc);
      ynou:=yanc+dif;
      if imp > 5 then print("x=",x,"yanc=",yanc,"ynou=",ynou,"dif=",dif); fi;
      yanc:=ynou;
   od;
else
   if x < 41. then yanc:=1.4; else yanc:=evalf(x0); fi;
   dif:=infinity;
   while abs(dif)/yanc > eps do
      dif:=(x-Li(yanc))*log(yanc);
      ynou:=yanc+dif;
      if imp > 5 then print("x=",x,"yanc=",yanc,"ynou=",ynou,"dif=",dif); fi;
      yanc:=ynou;
   od;
fi;
dif:=(x-Li(yanc))*log(yanc);
ynou:=yanc+dif;
Digits:=digits;
evalf(ynou);
end;

########################################################
# pi_r(x)
########################################################

pirdex:=proc(x) local i,f,ff,fp,fh,ffh,A,s,k,S,A4,C0,x1,Ah1,A4h1,
Ah2,A4h2,C0h,pix1,thx1,x1n,x1f,thx1n,pix1n,x2n,x2f,thx2n,pix2n,
x3n,x3f,thx3n,pix3n,alpha,ff220,ff221,ff222,ff226,ff227,phi,r0;
global t,r,L,L1;
# computes the expansion of pi_r(x) of Proposition 2.6
for i from 1 to 5 do
   S:=sum((k-1)!*s^(i-1-k)*t^s/(log(t))^k,k=1..i-1);
   f[i]:=1/(i-1)!*(s^(i-1)*Li(t^s)-S);
   fp[i]:=simplify(diff(f[i],t),assume=positive);
   print("i=",i,"f[i]=",f[i],"f'[i]=",fp[i]);
od;
x1n:=89967803; x1f:=evalf(x1n);
thx1n:=89953175.41601372639; pix1n:=227325003281036;
x2n:=767135587;   x2f:=evalf(x2n);
thx2n:= 767090801.98491111998  ; pix2n:=14752057842920905;
x3n:=19035709163; x3f:=evalf(x3n);
thx3n:= 19035493858.482419137;  pix3n:=7823414443039054263;
# these values have been computed by another program in C
print("################  upper bound #################");
ff:=-(s-1)*f[1]+f[2]+(s-1)*alpha*f[4]-alpha*f[5];
ff:=collect(subs(log(t)=L,ff),[L,s]);
ff226[0]:=24*coeff(ff,Li(t^s),1);
for i from 1 to 4 do
   ff226[i]:=24*normal(coeff(-ff,L,-i)/t^s);
od;
ff:=ff226[0]*Li(t^s)/24-sum(ff226[ii]*t^s/(24*L^ii),ii=1..4);
print("eq. (2.26) f=",ff);
A:=collect(x^s/L+alpha*x^s/L^4+subs(t=x,Li(x^s)=
      x^s/s/L+x^s/s^2/L^2+2*x^s/s^3/L^3+7*x^s/s^4/L^4,ff),[L,s]);
for i from 1 to 4 do ff227[i]:=24*normal(coeff(s^i*A,L,-i)); od;
A:=sum(ff227[ii]/(24*s^ii*L^ii),ii=1..4);
print("eq. (2.27)", A);
print("eq. (2.19)", collect(subs(s=r+1,A),[L,r]));
C0:=collect(subs(s=r+1,pix1-x1^(s-1)*thx1/log(x1)-subs(t=x1,L=L1,ff)),[L1,r]);
ff220[0]:=24*coeff(C0,Li(x1^(r+1)),1);
for i from 1 to 4 do
   ff220[i]:=24*normal(coeff(C0/x1^(r+1),L1,-i));
od;
C0:=pix1-x1^r*thx1/log(x1)+ff220[0]*Li(x1^(r+1))/24+sum(ff220[ii]*x1^(r+1)/(24*L1^ii),ii=1..4);
print("eq. (2.20) C0=",C0);
print("##############  lower bound #######################");
fh:=-(s-1)*f[1]+f[2]-alpha*(s-1)*f[4]+alpha*f[5];
ffh:=subs(alpha=-alpha,ff);
print("eq. (2.28) fh=",ffh);
C0h:=subs(alpha=-alpha,C0);
print("eq. (2.23) C0h=",C0h);
print(" ################ lower bound for r smaller than r0 #############");
Ah1:=collect(x^s/L-alpha*x^s/L^4+subs(t=x,Li(x^s)=
      x^s/s/L+x^s/s^2/L^2+2*x^s/s^3/L^3+6*x^s/s^4/L^4,ffh),[L,s]);
for i from 1 to 4 do ff221[i]:=24*normal(coeff(s^i*Ah1,L,-i)); od;
Ah1:=sum(ff221[ii]/(24*s^ii*L^ii),ii=1..4);
print("eq. (2.21) in s=", Ah1);
print("eq. (2.21) in r", collect(subs(s=r+1,Ah1),[L,r]));
print(" ############ lower bound for r greater than r0 ##################");
Ah2:=collect(x^s/L-alpha*x^s/L^4+subs(t=x,Li(x^s)=
      x^s/s/L+x^s/s^2/L^2+2*x^s/s^3/L^3+7*x^s/s^4/L^4,ffh),[L,s]);
for i from 1 to 4 do ff222[i]:=24*normal(coeff(s^i*Ah2,L,-i)); od;
Ah2:=sum(ff222[ii]/(24*s^ii*L^ii),ii=1..4);
print("eq. (2.22) in s", Ah2);
print("eq. (2.22) in r", collect(subs(s=r+1,Ah2),[L,r]));
print(" ############### Corollary 2.7: r=1 ####################");
print("eq. (2.30) with alpha=1", subs(s=2,alpha=1,A));
print("eq. (2.30) with alpha=0.5", subs(s=2,alpha=1/2,A));
print("eq. (2.30) with alpha=0.15", subs(s=2,alpha=3/20,A));
print("eq. (2.31) with alpha=1", subs(s=2,alpha=1,Ah1));
print("eq. (2.31) with alpha=0.5", subs(s=2,alpha=1/2,Ah1));
print("eq. (2.31) with alpha=0.15", subs(s=2,alpha=3/20,Ah1));
print("f(x1)=with alpha=1", evalf(subs(L=log(x1f),t=x1f,s=2,alpha=1,ff)));
print("f(x2)= with alpha=0.5", evalf(subs(L=log(x2f),t=x2f,s=2,alpha=1/2,ff)));
print("f(x3)= with alpha=0.15", evalf(subs(L=log(x3f),t=x3f,s=2,alpha=3/20,ff)));
print("C0= with alpha=1",evalf(subs(L1=log(x1f),pix1=pix1n,thx1=thx1n,x1=x1f,alpha=1,r=1,C0)));
print("C0= with alpha=0.5",evalf(subs(L1=log(x2f),pix1=pix2n,thx1=thx2n,x1=x2f,alpha=1/2,r=1,C0)));
print("C0= with alpha=0.15",evalf(subs(L1=log(x3f),pix1=pix3n,thx1=thx3n,x1=x3f,alpha=3/20,r=1,C0)));
print("C0h= with alpha=1",evalf(subs(L1=log(x1f),pix1=pix1n,thx1=thx1n,x1=x1f,alpha=1,r=1,C0h)));
print("C0h= with alpha=0.5",evalf(subs(L1=log(x2f),pix1=pix2n,thx1=thx2n,x1=x2f,alpha=1/2,r=1,C0h)));
print("C0h= with alpha=0.15",evalf(subs(L1=log(x3f),pix1=pix3n,thx1=thx3n,x1=x3f,alpha=3/20,r=1,C0h)));
phi:=3*t^4+8*t^3+6*t^2-1;
r0[1]:=fsolve(phi=24.,t=1..3);
r0[2]:=fsolve(phi=24/0.5,t=1..3);
r0[3]:=fsolve(phi=24/0.15,t=1..3);
print("r0(1)=",r0[1],"r0(1/2)=",r0[2],"r0(0.15)=",r0[3]);
end;

###############################################
# LEMMA 2.8
###############################################

lemma2_8:=proc()
local A,A1,A2,B,B1,B2,B3,B4,n0,L0,la0,x0,x0f,logx0,a,anum,b,bnum,c0,c0num,d,dnum,
c1,c2,c3,c1num,c2num1,c2num2,c3num1,c3num2,fdown,B3num,B4num;
global L,lambda;
print("##############  upper bound #######################");
A:=L*(1+lambda/(2*L))^2*(L+lambda+lambda/L+1)-(L+lambda+lambda/L)^2;
print("A=",A);
A:=collect(expand(A),lambda);
print("A ordered in lambda =",A);
print("##############  lower bound #######################");
n0:=2220822442581729257;
L0:=log(evalf(n0)); la0:=log(L0);
x0:=10^10+19; x0f:=evalf(x0); logx0:=log(x0f);
fdown:=sqrt(n*L)*(1+b*lambda/L);
print("f(n)=",fdown);
print("1+1/logx0+107/40/logx0^2=",1+1/logx0+107/40/logx0^2);
anum:=1.049; bnum:=0.365;c0num:=0.7;dnum:=0.88;
c1num:=0.44; c2num1:=0.26; c3num1:=0.; c2num2:=0.18; c3num2:=0.08;
print("anum=",anum,"bnum=",bnum);
print("2b/(1+b*la0/L0)=",2*bnum/(1+bnum*la0/L0),"c0num=",c0num);
A1:=L*(1+b*lambda/L)^2*(L+lambda+c0*lambda/L+a)-(L+lambda+c0*lambda/L)^2;
print("A1=",A1);
A2:=collect(expand(A1/(lambda*L)),[c0,L,lambda]);
print("A2=expand(A1/(lambda*L)",A2);
print("c0(2-2*b-b^2*lambda0/L0)=",c0num*(2-2*bnum-bnum^2*la0/L0),"dnum=",dnum);
B:=collect(coeff(A2,c0,0)-c0/L-d*lambda/L^2,[L,lambda]);
print("B=",B);
B1:=coeff(B,L,0)-c1/L; B2:=(coeff(B,L,-1)+c0-c2)/L;B3:=(coeff(B,L,-2)-c3*L)/L^2;
print("B=",B1+B2+B3);
print("a/lambda0-c1/L0=",anum/la0-c1num/L0);
print("********** case 1, lambda0 <= lambda <=4.3, c2=0.26, c3=0");
print("(b^2+2b-1)*lambda0+2ab-c2=",(bnum^2+2*bnum-1)*la0+2*anum*bnum-c2num1);
print("4.3*b^2*+(ab^2-d)=",4.3*bnum^2+anum*bnum^2-dnum);
print("********** case 2, lambda > 4.3, c2=0.18, c3=0.08");
print("4.3*(b^2+2b-1)+2ab-c2=",4.3*(bnum^2+2*bnum-1)+2*anum*bnum-c2num2);
print("4b^2/e^2-c3=",4*bnum^2*exp(-2.)-c3num2);
end;


########################################################
# § 2.4 The Riemann Zeta function
########################################################

deltak:=proc(mmax,dig,imp) local m; global del,fac,gaf,bb;
# calculates the sequence delta[m] for m <= mmax
# See Cohen, p. 207
Digits:=dig;
gaf[0]:=evalf(gamma); del[1]:=gaf[0];
fac[0]:=1; # factorial
for m from 1 to mmax do
   fac[m]:=m*fac[m-1]; gaf[m]:=evalf(gamma(m));
   del[m+1]:=(m+1)*gaf[m]/fac[m]+add(gaf[kk]*del[m-kk]/fac[kk],kk=0..m-1);
   bb[m+1]:=evalf(Zeta(m+1))*(1-1/2^(m+1))-1-del[m+1];
   if imp > 3 and m>= 2 then
      print("m=",m,"gamma=",gaf[m], "delta=",del[m],"sum 1/rho^m =",bb[m]);;
   fi;
od;
print();
end;

calculcn:=proc(nmax,imp) local n,f,tay,pol; global cn,t;
# computes the sequence c_n defined in Lemma 2.9
f:=((1-t^2)*(1-2*t))^(-1/2);
tay:=taylor(f,t=0,nmax+1);
pol:=convert(tay,polynom);
for n from 0 to nmax do
   cn[n]:=coeff(pol,t,n);
   if imp > 3 then
print("n=",n,"c_n=",cn[n],"=",evalf(cn[n],7),"c_n/2^n=",evalf(cn[n]/2^n,7));
   fi;
od;
end;

somrho:=proc(nmax,dig,imp) local S,n,a;
global bb,cn,t;
# This procedure computes the sum 1/|rho*(1+rho)| over the non trivial
# zeroes of the Riemann Zeta function by using the method of
# Henri Cohen, t. II, p. 272.
Digits:=dig;
if imp < 0 then
   calculcn(nmax,imp);
   deltak(nmax+2,dig,imp);
fi;
S[-1]:=0.;
for n from 0 to nmax do
   a:=cn[n]*bb[n+2];
   S[n]:=S[n-1]+a;
   if imp=-2 then
      print("n=",n,S[n],"a=", evalf(a,7));
   fi;
od;
S[nmax];
end;

##################################################
# Lemma 5.5, 5.6, 5.7, Proposition 5.9
###################################################

lemma5_5:=proc() local x0,S,k0,k,logx0s2;
x0:=evalf(10^10+19); logx0s2:=log(x0)/log(2.);
print("(1+7.5*10^(-7))*1+1.101/log(x0))=",  (1+7.5*10^(-7))*(1+1.101/log(x0)) );
for k0 from 3 to 33 do
   S:=1.05*add(x0^(1/k-1/3)/(1+1/k),k=3..k0-1)+logx0s2/x0^(1/3-1/k0);
   print("k0=",k0,"S=",S);
od;
end;

lemma5_6:=proc() local a,af;
for a in {0,1/3,1/2,1} do
   af:=evalf(a);
   print("a=",a,"delta_a=",evalf(af/2*Li(2^(af+2))-2^(af+2)/(2*log(2.))));
od;
end;

lemma5_7:=proc() local a,b;
a:=2.^(3/2)*(2-1/log(2.))/log(2.);
b:=evalf(Li(2.^(5/2)));
print("Estimate of the constant in formula (5.13)");
print("2^(3/2)*(2-1/log(2))/log(2)-li(2^(5/2)=",a-b);
end;

prop5_9:=proc() local P,lis;global nu;
P:=expand((1+0.365*nu)^3*(1+nu)*(1-0.37*nu)^2-(1+1.018*nu)^2*(1-0.3405*nu)^2);
print("(1+0.365 nu)^3(1+nu)(1-0.37 nu)^2-(1+1.018 nu)^2(1-0.3405)^2=",P);
lis:=fsolve(P=0,nu);
print("roots of P=",lis);
plot(P,nu=-1.2..3);
end;

###############################################
# Computation for small n
###############################################

hden:=proc(nmax,h,flot,imp) local n,p,pmax,j,k,t,a;
# This is the naive algorithm.
# compute h(n) for n<= nmax. Make pinit=2, pmax=1.4 sqrt(n log n)
# flot=true (to get a real h(n)), =false (for an integer h(n))
t:=time();
if flot then
   for n from 0 to nmax do h[n]:=1. od;
        else
   for n from 0 to nmax do h[n]:=1 od;
fi;
p:=2; pmax:=1.4*sqrt(nmax*log(evalf(nmax)));
while p <= pmax do
    for n from nmax by -1 to p do
      a:=p*h[n-p];
      if a > h[n] then h[n]:=a ; fi;
    od;
    if imp > 9 then print(`p=`,p,time()-t); fi;
p:=nextprime(p);
od;
if imp >= 2 then
   for n from 2 to nmax do
      if h[n]>h[n-1] then
         print("n=",n,"h(n)=",h[n]);
      fi;
   od;
fi;
end;

calculbn:=proc(n,imp) local x0,ti,nf,limm,rac,lis,liso;
global h,b; # assumes that h(n) has been computed
nf:=evalf(n); x0:=nf*log(nf); # ti:=time();
limm:=Lim1Newtonx0(n,x0,0); # calcul de li^(-1) (n)
rac:=sqrt(limm);
b[n]:=(rac-log(evalf(h[n])))/(nf*log(nf))^(1/4);
if imp > 5 then print("n=",n,"bn=",b[n]); fi;
b[n];
end;

bninterval:=proc(nmin,nmax) local n,nf,logn,mu,t; global h,b,c;
print("b_n is defined in (1.12) and mu_n by b_n=2/3-c+mu_n*(log log n)/log n");
hden(nmax,h,false,0); # computes h(n) for n <= nmax.
c:=0.046117644421509023827;
for n from nmin to nmax do
   nf:=evalf(n); logn:=log(nf);
   b[n]:=calculbn(n,0);
   mu:=fsolve(b[n]=2/3-c+t*log(logn)/logn,t);
   print("n=",n,"h=",h[n],"b=",evalf(b[n],10),"mu=",evalf(mu,10));
od;
end;

#########################################################
#  OK h(n)
##########################################################

ok:=proc(n) local bn,nf,logn,loglogn,loghn,rac; global c,compOK;
compOK:=compOK+1;
print("Procedure ok(n), n=",n);
nf:=evalf(n); logn:=log(nf); loglogn:=log(logn);
rac:=sqrt(Lim1Newtonx0(nf,nf*logn,0));
loghn:=loghden(n,0);
bn:=(rac-loghn)/(nf*logn)^(1/4);
if bn >= 2/3+c+0.77*loglogn/logn then
   print("n=",n,"does not satisfy inequality (iv) of Thm 1.1");
   return false;
else
   return true;
fi;
end;

good_interval:=proc(n1,n2) local n1f,n2f,logn1,logn2,loglogn2,rac,loghn1;
global superhden,c,comp;
comp:=comp+1;
if irem(comp,10000)=0 then
   print("comp=",comp,[n1,n2],"temps=",time()-ti);
fi;
n1f:=evalf(n1); n2f:=evalf(n2);
logn1:=log(n1f); logn2:=log(n2f); loglogn2:=log(logn2);
rac:=sqrt(Lim1Newtonx0(n2f,n2f*logn2,0));
loghn1:=loghden(n1,0);
evalb((rac-loghn1)/(n1f*logn1)^(1/4) <= 2/3+c+0.77*loglogn2/logn2);
end;

ok_rec:=proc(n1,n2) local nmed;
if n2-n1 >= 2 then
   if good_interval(n1,n2) then
      return true;
   else
      nmed:=iquo(n1+n2,2);
      if not ok_rec(nmed,n2) then
         return false;
      else
         return ok_rec(n1,nmed);
      fi;
   fi;
elif n1=n2 then
   return ok(n1);
else
   if ok(n2) then
      return ok(n1);
   else
      return false;
    fi;
fi;
end;

testOK:=proc(n1,n2) global comp,compOK,ti,c,superhden;
c:=0.046117644421509023827;
comp:=0; compOK:=0; ti:=time();
ok_rec(n1,n2);
print("comp=",comp,"compOK=",compOK);
end;

##########################################################
# Calcul de log h(n)
##########################################################

loghden:=proc(n,imp)  local indicemax,nmax,i,som,logh,p,pk,m;
global superhden,comp;  # Lire superhden.m
# This procedure computes log h(n).
indicemax:=nops(superhden);  nmax:=superhden[indicemax][2];
if n > nmax then error "n is too large"; fi;
if n < 5 then
   if n=1 then return 0.; fi;
   if n=2 then return log(2.) fi;
   if n= 3 or n=4 then return log(3.); fi;
fi;
i:=cherchez(superhden,n,1,indicemax,2);
if n-superhden[i][2] < superhden[i+1][2]-n then
   som:=superhden[i][2]; p:=superhden[i][1]; logh:=superhden[i][3];
   p:=nextprime(p); som:=som+p;
   while som <= n do
      logh:=logh+ log(evalf(p));
      p:=nextprime(p); som:=som + p;
   od;
   pk:=prevprime(p); m:=n-(som - p);
else
   som:=superhden[i+1][2]; p:=superhden[i+1][1]; logh:=superhden[i+1][3];
   som:=som-p;
   while som > n do
      logh:=logh - log(evalf(p));
      p:=prevprime(p); som:=som - p:
   od;
   logh:=logh-log(evalf(p)); pk:=prevprime(p); m:=n-som ;
fi;
logh:=logh+log(evalf(Gpkm(pk,m,0)[1]));
if imp > 4 then
   print("n=",n,"pk=",pk,"m=",m,"logh:=",logh);
fi;
logh;
end;

############################################################
# Fonction G(pk,m)
##############################################################

cherche:=proc(liste,m,ini,ifi) local a,b,c;
# "cherche" means look for; "liste" is a list of numbers in ascending order
# If  m >= liste[ifi], returns  ifi; if m < liste[ini] returns "error";
# if not,  returns c such that ini <= c < ifi and liste[c] <=m < liste[ifi].
# If the list is not-decreasing, returns the maximal such c
if nops(liste)=0 or liste[ini] > m or ini <= 0 or ifi > nops(liste) then error
"bad choice of m=%1, ini=%2, ifi=%3 or nopsliste=%4",m,ini,ifi,nops(liste);
fi;
if m >= liste[ifi] then return ifi; fi;
a:=ini; b:=ifi; # here, liste[ini] <= m and liste[ifi] < m
while b-a > 1 do
   c:=iquo(a+b,2);
   if liste[c] <= m then a:=c; else b:=c; fi; # liste[b] > m holds
od;
a;
end;

cherchez:=proc(liste,m,ini,ifi,j) local a,b,c;
# it is the same procedure than the preceding one;
# only, "liste" is a list of lists liste[i]=[a_i,b_i,c_i,...]
# For j given, the numbers liste[i][j] are increasing or not-decreasing;
# returns c such that ini<= c < ifi and liste[c][j] <= m < liste[c+1][j]
# in the same conditions that in the preceding procedure
if nops(liste)=0 or liste[ini][j] > m or ini <= 0 or ifi > nops(liste) then
   error "bad choice of m=%1, ini=%2, ifi=%3 or nopsliste=%4",
          m,ini,ifi,nops(liste);
fi;
if m >= liste[ifi][j] then return ifi; fi;
a:=ini; b:=ifi;
while b-a > 1 do
   c:=iquo(a+b,2);
   if liste[c][j] <= m then a:=c; else b:=c; fi;
od;
a;
end;


mj:=proc(j,K,r) local a; global S;
# This short procedure calculates m_j(P_r) defined in the article [9]
if j <= r then
   a:=S[K+j]-S[K]-S[r]+S[r-j];
else
   a:=infinity;
fi;
a;
end;

Gcomb:=proc(pk,Mmax,imp)
# imp is a printing parameter. If imp < 5 nothing is printed.
# This procedure computes G(pk,m) for m even, 0 <= m <= Mmax <= pkp1-3
# where pk is a prime and pkp1 is the prime following pk.
# The result is returned in a table GG[m]. Mmax is supposed not too large.
# First, we compute the primes P[1],...,P[K]=pk,...,P[Rhat]
# This procedure uses the combinatorial method explained in the article [9], § 8.1
# By induction on r and j, H(j,P[r],m) is saved in the memory HH[j,m].
local pkp1,P,r,Rhat,K,i,j,j1,mj1,m,bou,maxP,maxPgen,pkjr,
GG,HH,pro,rmes,mmes,nu,
tem,ti;   # variables used for the time
global S;
tem:=time(); ti:=tem; # initialization of the time
if not isprime(pk) then error "pk=%1 is not prime",pk;fi;
pkp1:=nextprime(pk);
if Mmax > pkp1-3 then error "Mmax=%1 is too big",Mmax; fi;
if Mmax < pkp1-pk then
   for m from 0 by 2 to Mmax do GG[m]:=1; od;
   return GG;
fi;
# Determination of P[1] (cf. (94))
P[1]:=nextprime(pkp1-Mmax-1); r:=1;
# Note that Mmax >= pkp1-pk implies P[1] <= pk.
# Determination of the primes P[r], r=1..K
while P[r] < pk do r:=r+1; P[r]:=nextprime(P[r-1]); od;
K:=r; # here, we have P[K]=pk
# Determination of large primes.
# These primes are limited to P[Rhat] (cf. § 8.4)
Rhat:=K+10; # 10 seems a reasonnable value.
for r from K+1 to Rhat+1 do P[r]:=nextprime(P[r-1]); od;
# P[Rhat+1] will only be used to compute maxPgen
S[0]:=0; #S is the sommatory function of P[r].
for r from 1 to Rhat do S[r]:=S[r-1]+P[r]; od;
if imp >= 6 then
   print(`Function Gcomb, pk=`,pk,`Mmax=`,Mmax,`Rhat=`,Rhat,
         `K=`,K,`P[K]=`,P[K],`small primes=`,  [P[r1] $ r1=1..K-1],
         `large primes=`, [P[r1] $ r1=K+1..Rhat],
         `time of computation of the primes=`,time()-ti); ti:=time();
fi;
pro:=1; # computation of the product pro:=P[K+1]...P[Rhat]
for r from K+1 to Rhat do pro:=pro*P[r]; od;
# The boolean variable "bou" will be used to eventually start again the
# "while" loop, if maxPgen (= rhs of (109)) is not large enough.
bou:=true;
while bou do
   for m from 0 by 2 to Mmax do HH[0,m]:=1; od;  # formula (108)
   for r from 1 to Rhat do
   # computation by induction of H(j,p[r],m)=HH[j,m]
      for j from min(Rhat-K,r) by -1 to max(1,r-K) do
         j1:=j-1; mj1:=mj(j,K,r-1); pkjr:=P[K+j]-P[r];
         for m from mj(j,K,r) by 2 to min(Mmax,mj1-2) do
            HH[j,m]:= P[r]*HH[j1,m-pkjr]; # formula (7.21)
         od;
         for m from mj1 by 2 to Mmax do
            HH[j,m]:= min(HH[j,m],P[r]*HH[j1,m-pkjr]); # formula (106)
         od;
      od; # (j)
      if imp >= 8 then
         print(`r=`,r,`time=`,time()-ti,time()-tem); ti:=time();
      fi;
   od;  # (r)
   # Determination of maxPgen. This parameter is an upper bound
   # for the large primes (cf. formula (109)).
   maxPgen:=0;
   for m from 0 by 2 to Mmax do
      GG[m]:=pro/HH[Rhat-K,m]; # computation of G(pk,m) by formula (102)
      if GG[m] > 1 then
         maxP:=min(m+pk,floor(GG[m]*m/(GG[m]-1))); # maximum for m
      else
         maxP:=m+pk; # maximum for m
      fi;
      maxPgen:=max(maxPgen,maxP); # maximum for all m's
   od; # (m)
   # If maxPgen >= P[Rhat+1], P[Rhat+1] can occur, we have to continue.
   if maxPgen < P[Rhat+1] then
      bou:=false; # the computation is over
      if imp >= 5 then
         print(`P[Rhat+1]=`,P[Rhat+1],`is > maxPgen=`,maxPgen,
               `Rhat-K=`,Rhat-K,`real increment=`,
                cherche([P[r1] $ r1=K..Rhat+1],maxPgen,1,Rhat+2-K)-1,
                `the computation is over`,`time=`,time()-ti);  ti:=time();
      fi;
   else
      if imp >= 5 then
         print(`P[Rhat+1]=`,P[Rhat+1],`is <= maxPgen=`,maxPgen,
             `we start again with more primes`,`time=`,time()-ti); ti:=time();
      fi;
      r:=Rhat+1; # We have to add primes
      while P[r] <= maxPgen do
         pro:=pro*P[r]; S[r]:=S[r-1]+P[r];
         r:=r+1; P[r]:=nextprime(P[r-1]);
      od;  # here, P[r] > maxPgen holds
      Rhat:=r-1;
   fi;
   # If maxPgen was too small, we start again the loop on "bou"
od;  # (bou)
if imp >= 5  then # we measure the largest prime which has been used
   rmes:=K+1;
   for m from 0 by 2 to Mmax do
      nu:=numer(GG[m]);
      for r from rmes+1 to Rhat do
         if irem(nu,P[r])=0 then
            rmes:=r; mmes:=m;
         fi;
      od; # (r)
   od; # (m)
   print(`pmaxreel=`,P[rmes],`increm=`,rmes-K,`mmes=`,mmes,`total time=`,
          time()-tem);
fi;
GG;
end;

Gpkm:=proc(pk,m,imp) local mm,G1,G2,delta1max,maxd,maxdelta,pkp1,pkp2,q,d,oldd,
Gprov,Gpr,Gout,mot,bou,qhat,delta,pkp1mmm,dmax;
# This procedure computes G(pk,m) as explained in the article, § 9.
if not isprime(pk) then error "pk=%1 is not prime",pk;fi;
pkp1:=nextprime(pk);
if m < 0 or m >= pkp1 or not type(m,integer) then error "m=%1 is wrong",m; fi;
if pkp1-m in {1,2}  then
   return [pkp1/2,["m is equal to 1 or 2"]];
fi;
mm:=m-irem(m,2); # mm is the greatest even number <= m
pkp1mmm:=pkp1-mm;
if isprime(pkp1mmm) then
   return [pkp1/(pkp1mmm),[`p_{k+1}-m is prime, delta=0`]];
fi;
pkp2:=nextprime(pkp1);
delta1max:=2*trunc(1.275*log(evalf(pk))^2);
# delta1max is an upper bound for delta_1(pk) defined in the paper, cf. (117)
if irem(m,9,'qq')=0 then maxdelta:=2*qq-2 else maxdelta:=2*qq; fi;
# here, maxdelta is the greatest even integer < 2m/9, cf. (112)
maxdelta:=min(maxdelta,delta1max);
maxd:=pkp2-pkp1+3*maxdelta/2; maxd:=maxd+irem(maxd,2); # maxd is even
# maxd is an upper bound for the d's for which we need to calculate G(pkp1,d)
G1:=Gcomb(pkp1,maxd,imp);
q:=nextprime(pkp1mmm); # here, we are sure that pkp1-mm is not prime
# q is the smallest prime which may occur in the denominator of G(pk,m)
d:=q-pkp1mmm; bou:=false; Gprov:=1;
while d <= maxdelta and not bou do
   Gpr:=pkp1/q*G1[d]; # formula (113)
   if Gpr > Gprov then Gprov:=Gpr; dmax:=d; fi;
   bou:=evalb(G1[d] >= 1+d/pkp1); # formula (111)
   q:=nextprime(q); oldd:=d; d:=q-pkp1mmm;
od;
if bou then # here, we have found delta
   delta:=oldd;
   qhat:=iquo(pkp1*pkp2*(pkp1mmm+delta),(pkp1+delta)*(pkp1-3*delta/2));
   while q <= qhat do
      Gpr:=pkp1/q*G1[d]; # formula (113)
      if Gpr > Gprov then Gprov:=Gpr; dmax:=d; fi;
      q:=nextprime(q); d:=q-pkp1mmm;
   od;
   # we expect that dmax=delta often holds; below, we check that
   if dmax < delta then
      mot:=`exceptionnal case one`;
   elif dmax > delta then
      mot:=`exceptionnal case two`;
   else
      mot:=``;
   fi;
   Gout:=[Gprov,[`maxd=`,maxd,`delta=`,delta,mot]];
else # here, we have not found delta
   if mm > 9*delta1max/2 then
      mot:=`We should have found delta`;
   else
      mot:=``;
   fi;
   G2:=Gcomb(pk,mm,imp);
   Gprov:=G2[mm];
   Gout:=[Gprov,[`maxm=`,mm,`delta not found`,mot]];
fi;
Gout;
end;

###############################################
# Computation of Prop. 5.5
###############################################

testppi:=proc(i) local lis,cp23,cm23,pk,sigk,
thk,pkp1,sigkp1,logsigkp1,thkp1,raclik,raclikp1,bsigk,tauk,muk,
sigklog14,logsigk,xinit,lim1,bptauk,bmin,bmax,nbmin,nbmax,bouB,bouC;
global superhden;
# For 1 <= i <= 29329 superhden[i][1] contains a prime q_i.
# For all primes pk satisfying pk < q_{i+1}
# this procedure calculates b_{sigma_k} with sigma_k=pi_1(pk),
# bptauk=b_{sigma_k} + tau_k and muk defined by
# bptauk=2/3+c +muk*log(log(sigma_{k+1}))/log(sigma_{k+1}).
# Thisprocedure calculates the minimum "bmin" and the maximum "bmax"
# of b_{sigma_k} for sigma_k >= 100.
# This procedure also displays the largest values of p_k such that
# bptauk > 1.044 and muk > 0.77.
cm23:=2/3-0.046117644421509023827;
cp23:=2/3+0.046117644421509023827;
bmin:=infinity; bmax:=-infinity;
lis:=superhden[i+1];
pkp1:=lis[1]; sigkp1:=lis[2]; thkp1:=lis[3];
logsigkp1:=log(evalf(sigkp1));
lim1:=Lim1Newtonx0(sigkp1,sigkp1*logsigkp1,0); raclikp1:=sqrt(lim1);
pk:=prevprime(pkp1);
bouB:=true; bouC:=true;
while pk > 2 do
   sigk:=sigkp1-pkp1; thk:=thkp1-log(evalf(pkp1));
   logsigk:=log(evalf(sigk)); sigklog14:=(sigk*logsigk)^(1/4);
   lim1:=Lim1Newtonx0(sigk,lim1,0); raclik:=sqrt(lim1);
   bsigk:=(raclik-thk)/sigklog14;
   tauk:=(raclikp1-raclik)/sigklog14; bptauk:=bsigk+tauk;
   muk:=solve(bptauk=cp23+t*log(logsigkp1)/logsigkp1,t);
   if sigk >= 100 then
      if bsigk < bmin then bmin:=bsigk; nbmin:= sigk;
     fi;
      if bsigk > bmax then bmax:=bsigk; nbmax:= sigk; fi;
   fi;
   if bptauk > 1.044 and bouB then
      print("sigk=",sigk,"sigkp1=",sigkp1,"bptauk=",evalf(bptauk,10)); bouB:=false;
   fi;
   if muk > 0.77 and bouC then
      print("sigk=",sigk,"sigkp1=",sigkp1,"muk=",evalf(muk,10)); bouC:=false;
   fi;
   sigkp1:=sigk; logsigkp1:=logsigk; thkp1:=thk; pkp1:=pk; raclikp1:=raclik;
   pk:=prevprime(pk);
od;
print("min of b_{sigma_k} (for sigma_k >= 100) =",evalf(bmin,10),"for sigma_k=",nbmin);
print("max of b_{sigma_k} (for sigma_k >= 100) =", evalf(bmax,10),"for sigma_k=",nbmax);
end;


##########################################################

testr:=proc(i,imp) local A,B,C,pA,pB,pC,nA,nB,nC,lis,pmin,cp23,cm23,pk,sigk,
thk,pkp1,sigkp1,thkp1,raclik,raclikp1,bsigk,tauk,muk,
sigklog14,logsigk,xinit,lim1,bptauk,bouA,bouB,bouC;
# For 1 <= i <= 29329 superhden[i][1] contains a prime q_i.
# For all primes p satisfying q_i <= p < q_{i+1}
# this procedure calculates b_sigma with sigma=pi_1(p)
# and the minimum A of these b_sigma.
global superhden;
cm23:=2/3-0.046117644421509023827;
cp23:=2/3+0.046117644421509023827;
lis:=superhden[i+1]; pmin:=superhden[i][1];
pkp1:=lis[1]; sigkp1:=lis[2]; thkp1:=lis[3];
lim1:=Lim1Newtonx0(sigkp1,sigkp1*log(evalf(sigkp1)),0); raclikp1:=sqrt(lim1);
pk:=prevprime(pkp1);
A:=infinity; B:=-infinity; C:=-infinity;
bouA:=true; bouB:=true; bouC:=true;
while pkp1 > pmin do
   sigk:=sigkp1-pkp1; thk:=thkp1-log(evalf(pkp1));
   logsigk:=log(evalf(sigk)); sigklog14:=(sigk*logsigk)^(1/4);
   lim1:=Lim1Newtonx0(sigk,lim1,0); raclik:=sqrt(lim1);
   bsigk:=(raclik-thk)/sigklog14;
   if bsigk < A then A:=bsigk; pA:=pk; nA:=sigk; fi;
   tauk:=(raclikp1-raclik)/sigklog14; bptauk:=bsigk+tauk;
   if bptauk > B then B:=bptauk; pB:=pk; nB:=sigk;fi;
   muk:=solve(bptauk=cp23+t*log(logsigk)/logsigk,t);
   if muk > C then C:=muk; pC:=pk; nC:=sigk; fi;
   if imp=0 then
      if bsigk < cm23 and bouA then
         print("sigk=",sigk,"sigkp1=",sigkp1,"bsigk=",bsigk); bouA:=false;
      fi;
      if bptauk > 1.044 and bouB then
         print("sigk=",sigk,"sigkp1=",sigkp1,"bptauk=",bptauk); bouB:=false;
      fi;
      if muk > 0.77 and bouC then
         print("sigk=",sigk,"sigkp1=",sigkp1,"muk=",muk); bouC:=false;
      fi;
   fi;
   sigkp1:=sigk; thkp1:=thk; pkp1:=pk; raclikp1:=raclik;
   pk:=prevprime(pk);
od;
[A,B,C,pA,pB,pC,nA,nB,nC];
end;

testh:=proc(imin,imax,pas) local i,lis,ti; global lisA,lisB,lisC;
ti:=time();
lisA:=[]; lisB:=[]; lisC:=[];
for i from imin to imax do
   lis:=testr(i,1);
   lisA:=[[lis[1],i,lis[4],lis[7]],op(lisA)];
   lisB:=[[lis[2],i,lis[5],lis[8]],op(lisB)];
   lisC:=[[lis[3],i,lis[6],lis[9]],op(lisC)];
   if irem(i,pas)=0 then
      print("i=",i,superhden[i][1..2],"temps=",time()-ti);
   fi;
od;
end;