read "superch2.m": 
# contains the list "listsuperchE2" and the variable "indicemax" (= the number
# of terms of the list "listsuperchE2")
# which are described and used in the procedure "superchamp" below.

# The procedure "functLandau(n,imp)" computes the standard factorization 
# into primes of g(n) according to the scheme given in § 1.3 of the paper.
# To get the numerical value, use "numvalueofg(n)".
# n should satisfy 7 <= n <= 1073484873411753. If n <= 165, the algorithm
# sometimes does not work.
# "imp" is a printing variable between 0 and 9. If imp = 0, nothing is printed;
# if imp = 1, are printed N = the superchampion preceding g(n), n, ell(g(n)),
# the value  of the fraction g(n)/N and its factorisation into primes.
# For 2 <= imp <= 4, more parameters are printed. The values of imp 
# betwen 5 and 9 are used to print intermediate results.

functLandau:=proc(n,imp)  
local
temps,tempar,tN,tp1,tp2,tp3,tb1,tb2,tn,tgr,  # time variable
np1,np2,np3,nn1,nn2, # count the lengths of the lists
i,gofnoverN,lis,Boverho,PI,PIf,pif,m,mf,nmlgofn,pom,gg,suff,omf,
Bprime,   # guessed value for B; called B' in the paper, see § 7.3
Bexact,   # value of ben(g(n))+n-ell(g(n))
pmin;     # the smallest prime dividing the superchampion N with v_p(N)=1
global
rho,      # parameter s.t. ell(N_rho) <= n < ell(N_rho^+)  (cf. § 5.1)
x1,       # implicitely defined by x1/log(x1)=rho (cf. formula (28)
x2,       # implicitely defined by (x2^2-x2)/log(x2)=rho (cf. formula (28)
pk,       # the greatest prime factor of N_rho (noted p_k in the paper)
listexpo, # the list where each term [p,alpha] gives the exponent alpha of the
          # prime p in the superchampion N_rho (for all p with alpha >= 2)
lofN,     # =ell(N_rho) which is  <= n (by (32))
lofNprime,# =ell(N'_rho) > n; often, lofNprime=lofN+nextprime(pk)
lprime,   # lofNprime-lofN
pprime,   # quotient of Nprime_rho by N_rho; often, pprime=nextprime(pk)
B,        # an upper bound for ben(g(n)) + n -ell(g(n)) (cf. § 7.3 and (42))
Baug,     # a slightly increased value of B to avoid rounding errors
B1,       # We assume Baug <= B1 in order to apply Proposition 2, § 6.1
listpref1,listpref2, # provisional values of "listpref"
listpref, # list of the possible plain prefixes pi (ben(pi*N) <= B) 
          # listpref is noted {curly D}(B) in the article.
          # Each plain prefix pi is represented by a list of 3 terms
          # [pi,ell(N*pi)-ell(N),ben(N*pi)]
listprefnorm1, # provisional value of "listprefnorm" before "fight" § 7.7
listprefnorm; # list of possible normalized prefixes after "fight"
              # (noted {curly N}) § 7.9
          # Each prefix PI is represented by a list of 6 numbers: 
          # [PI,ell(N*PI)-ell(N),ben(N*PI),omega,p_{k+omega},pi]
 temps:=time();tempar:=temps;
################
# First step, determination of the superchampion N (cf. § 5.2)
################
superchamp(n); # computes N satisfying (32) : ell(N) <= n < ell(N')
 tN:=time()-tempar;tempar:=time(); # tN is the time of computation of N
if imp >= 1 then
   pmin:=nextprime(listexpo[nops(listexpo)][1]);
   print(`N=`,op(impression(listexpo,pmin,pk)));
fi;
if imp >= 2 then
   print(`n=`,n,`ell(N)=`,lofN,`n-ell(N)=`,n-lofN);
   print(`rho=`,rho,`prime number corresponding to rho=`,pprime);
fi;
################
# Second step : first determination of plain prefixes.
# Bprime represents B' as described in § 7.3
################
x2:=fsolve((t^2-t)/log(t)=rho,t=2...x1); # formula (28)
B1:=min(x2^2-2.*x2,x1/2.-sqrt(x1)); # cf. (34) 
if n >= 10^10 then Bprime:=0.5*rho;  # cf. § 7.3 for the choice of Bprime
elif n >= 2485 then Bprime:=rho; 
else Bprime:=0.9999*B1-0.0001; # a little less than B1
fi;
listpref1:=prefix(Bprime,imp);
 tp1:=time()-tempar;tempar:=time(); np1:=nops(listpref1);
################
# Third step: first upper bound for the benefit of g(n).
################
B:=majben(listpref1,n,imp); # cf. § 7.3
Baug:=1.0001*B+0.0001; # to avoid rounding errors
 tb1:=time()-tempar;tempar:=time(); 
 tp2:=-1; tb2:=-1; np2:=-1; # initialization of these variables in case
                            # of Baug <= Bprime
if imp >= 3 then
   print(`time for N=`,tN,`time of 1st computation of prefixes=`,tp1,
          `number of prefixes=`,np1);
   Boverho:=B/rho;
   print(`first computation of B=`,B,`B/rho=`,evalf(Boverho,3),`time=`,tb1);
fi;
################
# Fourth step : if the upper bound B is not good enough, we start again
################
if Baug > B1 then 
   if imp >= 1 then
      print(`it does not work for n=`,n,`since B=`,B,`Baug=`,Baug,
            `are greater than B1=`,B1,`try an n > 165`);
   fi;
   return [false,`n=`,n,`B=`,B,`Baug=`,Baug,`B1=`,B1]; 
fi;
# Now, we have Baug <= B1.
if Baug > Bprime then  # we have to compute again prefixes, cf. § 7.3
   listpref2:=prefix(Baug,imp);
    tp2:=time()-tempar;tempar:=time(); np2:=nops(listpref2);
   B:=majben(listpref2,n,imp);
   Baug:=1.0001*B+0.0001; # to avoid rounding errors
    tb2:=time()-tempar;tempar:=time(); 
   listpref:=shorten(listpref2,Baug); # removing the prefixes pi of "listpref2"
                                      #  s.t. ben(N*pi) > Baug
else
   listpref:=shorten(listpref1,Baug); # removing the prefixes pi of "listpref1"
                                      #  s.t. ben(N*pi) > Baug
fi;
 tp3:=time()-tempar;tempar:=time(); np3:=nops(listpref);
if imp >= 3 then
   if tp2 < 0 then 
      print(`no second computation of prefixes`); 
   else
 print(`time of 2nd computation of prefixes=`,tp2,`number of prefixes=`,np2);
 print(`second computation of B=`,B,`time=`,tb2);
   fi:
   Boverho:=B/rho; 
   print(`B=`,B,`B/rho=`,evalf(Boverho,3),`time of shortening`,tp3,
          `number of prefixes=`,np3);
fi;
################
# Fifth step : construction of normalized prefixes (§ 7.6 and 7.7).
################
listprefnorm1:=prefnorm(listpref,n,imp); nn1:=nops(listprefnorm1);
listprefnorm:=fight(listprefnorm1,pk,n,imp); # cf. § 7.9
 tn:=time()-tempar; tempar:=time(); nn2:=nops(listprefnorm); 
if imp >= 3 then
   print(`time of computation of normalized prefixes=`,tn);
   print(`number of normalized prefixes=`,nn1,`after fight=`,nn2);
fi;
################
# Sixth and last step : determination of the suffix and of g(n)
################
gofnoverN:=0; # cf. § 7.8
for i from 1 to nops(listprefnorm) do
   lis:=listprefnorm[i]; PI:=lis[1]; m:=n-lofN-lis[2]; pom:=lis[5];
   gg:=Gpkm(pom,m,imp);
   if imp >= 4 then
      print(`Function G(p,m); omega=`,lis[4],`p=`,pom,`m=`,m,`G=`,
             ifactor(gg[1]),gg[2]);
   fi; 
   if PI*gg[1] > gofnoverN then 
      gofnoverN:=PI*gg[1]; PIf:=PI; pif:=lis[6]; 
      suff:=PIf/pif*gg[1]; nmlgofn:=m-ellfrac(gg[1]); omf:=lis[4]; mf:=m;
   fi;
od;
 tgr:=time()-tempar;
if imp >= 4 then
   print(`time of computation of the suffix=`,tgr);
   Bexact:= evalf(n-lofN)-rho*log(evalf(gofnoverN)); Boverho:=Bexact/rho;
   print(`Bexact=ben(g(n))+n-ell(g(n))=`,Bexact,
          `Bexact/rho=`,evalf(Boverho,3));
fi;
if imp >= 1 then
   print(`n=`,n,`g(n)/N=`,gofnoverN,`g(n)/N=`,ifactor(gofnoverN),
         `n-ell(g(n))=`,nmlgofn,`total time=`,time()-temps);
fi;
if imp >= 2 then
   print(`m=`,mf,`plain prefix=`,ifactor(pif),`omega=`,omf,`PI/pi=`,
         ifactor(PIf/pif),`suffix=`,ifactor(suff));
fi;
gofnoverN;
end; # of functLandau

ell:=proc(p,k) 
if k=0 then 0 elif k > 0 then p^k 
       else error "the variable k= %1 is non-positive for p= %2", k,p; 
fi;
end;

ellfrac:=proc(M) local S,liste,i;
# computes the function "ell" for the fraction M
if type(M,fraction) then ellfrac(numer(M))-ellfrac(denom(M));
elif type(M,integer) then
   liste:=op(2,ifactors(M));
   S:=0;
   for i from 1 to nops(liste) do
      S:=S+ell(op(liste[i]));
   od;
   S;
else error "bad type of M=%1",M; 
fi;
end;

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;


impression:=proc(liste,p1,p2) local i,nb,com,T,lisres,mot;
# This procedure allows the printing of the factorization of superchampions.
nb:=nops(liste); mot:=` ... `;
if nb < 5 or liste[nb-2][2] > 2 then
   if p1=p2 then
      lisres:=[op(liste),p1];
   elif p2=nextprime(p1) then 
      lisres:=[op(liste),p1,p2];
   else
      lisres:=[op(liste),p1,mot,p2];
   fi;
else
   com:=0; i:=1;
   while liste[i][2] > 2 do
      if liste[i][2] > liste[i+1][2] then
         com:=com+1; T[com]:=liste[i]; i:=i+1;
      elif liste[i][2] > liste[i+2][2] then
         com:=com+1; T[com]:=liste[i]; 
         com:=com+1; T[com]:=liste[i+1]; i:=i+2;
      else
         com:=com+1; T[com]:=liste[i]; 
         com:=com+1; T[com]:=mot;        
         while liste[i][2]=liste[i+1][2] do i:=i+1; od;
         com:=com+1; T[com]:=liste[i]; i:=i+1;
      fi;
   od;
        com:=com+1; T[com]:=liste[i]; 
        com:=com+1; T[com]:=mot; 
        com:=com+1; T[com]:=liste[nb]; 
        com:=com+1; T[com]:=p1; 
        com:=com+1; T[com]:=mot; 
        com:=com+1; T[com]:=p2;
        lisres:=[seq(T[icom],icom=1..com)];
fi;
lisres;
end; 




superchamp:=proc(n) local i,aile,ii,q,alfaq,maxq,n0,ltab,pkp1,t,p,oldp,pmin;
global pk,pprime,alfaprime,lofN,lofNprime,lprime,rho,alfa,listexpo,x1,
nmaximum,  # =1073484873411753 max value of n for this algo to compute g(n)
indicemax,# =1360, number of terms of the list "listsuperchE2" 
listsuperchE2; # each element of this list is of the form [q,j,p,ell(N_\rho^+];
          # rho=(q^{j+1}-q^j)/log q and p is the greatest prime s.t.
          # p/log p < rho. cf. The superchampion algorithm, § 5.2
# This procedure  determines the two consecutive superchampion 
# numbers N and N' such that ell(N) <= n < ell(N') (formula (32). 
nmaximum:=listsuperchE2[indicemax][4]-1;
if n < 7 or n > nmaximum then 
   error "IMPOSSIBLE,n=%1 is not in the interval [ 7 .. %2]",n,nmaximum;
fi;
# We determine the two consecutive elements of listsuperchE2 sandwiching n
i:=cherchez(listsuperchE2,n,1,indicemax,4);
q:=listsuperchE2[i+1][1]; alfaq:=listsuperchE2[i+1][2];
aile:=listsuperchE2[i+1][4]-q^(alfaq-1)*(q-1);
if aile <= n then # here, n is just a little smaller than listsuperchE2[i+1][4]
   rho:=evalf(q^(alfaq-1)*(q-1)/log(q)); # and rho belongs to the set E' 
   pk:=listsuperchE2[i+1][3]; pkp1:=nextprime(pk);
   pprime:=q;
   alfaprime:=alfaq-1;
   lofN:=aile;
   lofNprime:=listsuperchE2[i+1][4]; lprime:=lofNprime-lofN;
   x1:=fsolve(t/log(t)=rho,t=0.9999*evalf(pk)..1.0001*evalf(pkp1));
else
   n0:=listsuperchE2[i][4]; oldp:=listsuperchE2[i][3]; p:=nextprime(oldp);
   n0:=n0+p;
   while n0 <= n do
      oldp:=p; p:=nextprime(p);
      n0:=n0+p;
   od;
   lofNprime:=n0; lofN:=n0-p; lprime:=lofNprime-lofN;
   pk:=oldp; pprime:=p;
   x1:=evalf(pprime);
   rho:=x1/log(x1);
   alfaprime:=0;
fi;
maxq:=0; # computation of "listexpo"
for ii from 1 to i do
   q:=listsuperchE2[ii][1]; alfa[q]:=listsuperchE2[ii][2];
   if q > maxq then maxq:=q; fi;
od;
q:=2; ii:=0;
while q <= maxq do
   ii:=ii+1;
   ltab[ii]:=[q,alfa[q]];
   q:=nextprime(q);
od;
pmin:=q;
listexpo:=[ltab[iii] $iii=1..ii];
[n,lofN,lofNprime,n-lofN,lprime,pmin,pk,rho,x1,pprime,listexpo];
end;

evalN:=proc() local p,alpha,S,N,palpha,k;
global pk,listexpo;
# computes (for small n's) the numerical value of N, the greatest superchampion
# such that ell(N) <= n. This proc is not used in the computation of g(n)
N:=1; S:=0;
for k from 1 to nops(listexpo) do
   p:=listexpo[k][1]; alpha:=listexpo[k][2]; palpha:=p^alpha;
   N:=N*palpha; S:=S+palpha;
od;
while p < pk do
   p:=nextprime(p); N:=N*p; S:=S+p;
od;
[S,N];
end;

numvalueofg:=proc(n) local gofnoverN;
# computes the number  g(n), for small n's
gofnoverN:=functLandau(n,0);
evalN()[2]*gofnoverN;
end;

superchampp:=proc(x) local i,aile,aileprime,S,p,p1,q,alfaq;
global listsuperchE2;
# Let p be the greatest prime <= x (supposed to be an integer).
# This proc computes S=ell(N), where N is the greatest superchampion
# such that its greatest prime factor is p. 
# This proc is not used in the computation of g(n)
i:=cherchez(listsuperchE2,x,1,nops(listsuperchE2),3);
aile:=listsuperchE2[i][4];
q:=listsuperchE2[i+1][1];alfaq:=listsuperchE2[i+1][2];
aileprime:=listsuperchE2[i+1][4]-q^(alfaq-1)*(q-1);
p1:=listsuperchE2[i][3];
p:=nextprime(p1);
S:=aile;
while p <= x do
   S:=S+p;
   p:=nextprime(p);
od:
S;
end;

elaguage:=proc(liste) local comp,icom,a,maxl,i,T;
# elaguage means pruning. Each term of "liste" is a list of at least 2
# elements : [a[1],a[2],...]. "liste" is supposed to be ordered so that
# the first element of its terms are in increasing order.
# This procedure removes of "list" the terms a st there exists an other term
# b satisfying a[1] < b[1]  a[2] >= b[2].
comp:=nops(liste);
if comp <= 1 then return liste fi;
T[comp] :=liste[comp]; # the greatest term cannot be taken off
icom:=comp; maxl:=T[comp][2];
for i from comp -1 by -1 to 1 do
   a:=liste[i];
   if a[2] < maxl then # here the term a is kept
      icom:=icom-1;
      T[icom]:=a;
      maxl:=a[2];
   fi;
od;
[seq(T[scom],scom=icom..comp)];
end;

fusionelaguage:=proc(na,A,nb,B) local ic,C,ia,ib,maxl,nc,res;
# The terms of the two lists A and B of length na and nb have the same form
# [a_i,a'_i,...] and [b_j,b'_j,...]. In each of these two lists
# the two first elements are ordered in increasing order, i.e. they satisfy
# a_i < a_{i+1}, a'_i < a'_{i+1}, b_j < b_{j+1} and b'_j < b'_{j+1}.
# One supposes that the first elements a_i and b_j are all distinct.
# The lists A and B are merged, yielding a new list C=[...,[c_k,c'_k,...],...]
# with c_k < c_{k+1} . Further, the list C is pruned, i.e. the
# element C[i]=[c_i,c'_i,...] for which there exists C[j] with i < j,
# c_i < c_j and c'_i >= c'_j is taken off.
if na=0 then return B; fi;
if nb=0 then return A; fi;
ia:=na; ib:=nb; 
nc:=na+nb; ic:=nc+1;
maxl:=infinity;
while ia > 0 and ib > 0 do
   if A[ia][1] < B[ib][1] then
      if B[ib][2] < maxl then
         ic:=ic-1;
         C[ic]:=B[ib];
         maxl:=B[ib][2];
      fi;
      ib:=ib-1;
   else
      if A[ia][2] < maxl then
         ic:=ic-1;
         C[ic]:=A[ia];
         maxl:=A[ia][2];
      fi;
      ia:=ia-1;
   fi;
od;
# at the end of the above loop, ia and ib cannot vanish simultaneously 
if ia=0 then
   while B[ib][2] >= maxl and ib >= 1 do ib:=ib-1;od;
   res:= [seq(B[ibb],ibb=1..ib),seq([icc],icc=ic..nc)];
   # if ib = 0, seq(B[ibb],ibb=1..ib) is empty.
else   # here, ib=0
   while A[ia][2] >= maxl and ia >= 1 do ia:=ia-1;od;
   res:= [seq(A[iaa],iaa=1..ia),seq(C[icc],icc=ic..nc)];
fi;
res; 
end;


prefix:=proc(Bprime,imp) local p,j,U,compD,compU,icom,beta,pbeta,
alpha,palpha,b,delta,listel,lisD,lisU,rhologp,p2mp,lisUp,ti;
global pk,listexpo,rho,x1;
# This procedure computes the set D(B') as explained in the paper § 7.2.
# All the computed prefixes have the form pi=U/V satisfying (44)
# Note that more prefixes can exist, if Bprime is > B1. 
# Each prefix is saved in the form [g,lg,ben]; g is an element of U_j
# or D_j, lg=ell(Ng)-ell(N), and ben=ben(gN) <= Bprime.
# B' is supposed to be >= 0. If imp >= 9, partial results are printed.
ti:=time();
lisD:=[[1,0,0.]];   # 1 is the only element of D_0
compD:=1;   # compD counts the elements of D_0, D_1,...
for j from 1 to nops(listexpo) do # here we consider the primes < x_2
   listel:=listexpo[j]; 
   p:=listel[1]; # p is equal to p_j
   alpha:=listel[2]; # alpha is v_p(N)
   palpha:=p^alpha; rhologp:=rho*log(evalf(p));
   compU:=compD; # compD is the number of elements of D_{j-1}
      # First, we put all the elements of D_{j-1} in U_j so that 
      # compU is equal to compD.
      # Considering the exponents greater than alpha
   beta:=alpha+1; pbeta:=palpha*p;
   b:=evalf(pbeta-palpha)-rhologp; #b is ben(Np^gamma) with gamma=beta-alpha
   while b <= Bprime do # if b > Bprime, then (54) cannot hold
      for icom from 1 to compD do
         delta:=lisD[icom]; # delta runs over D_j
         if b+delta[3] <= Bprime then # this is formula (54)
            compU:=compU+1; # we add a new element to set U_j
            U[compU]:=[pbeta/palpha*delta[1],delta[2]+pbeta-palpha,b+delta[3]];
         fi;
      od;
      beta:=beta+1; pbeta:=pbeta*p;
      b:=evalf(pbeta-palpha)-evalf(beta-alpha)*rhologp;
         # by Lemma 3, b is increasing on beta
   od; 
         # Now, the exponents smaller than alpha
   beta:=alpha-1; pbeta:=palpha/p;
   b:=evalf(pbeta-palpha)+rhologp; #b is ben(Np^gamma) with gamma=beta-alpha <0
   while b <= Bprime and beta >= 0 do
      for icom from 1 to compD do
         delta:=lisD[icom]; # delta runs over D_j
           if b+delta[3] <= Bprime then
            compU:=compU+1; # we add a new element to set U_j
      U[compU]:=[pbeta/palpha*delta[1],delta[2]+ell(p,beta)-palpha,b+delta[3]];
           fi;
      od;
      beta:=beta-1; pbeta:=pbeta/p;
      b:=evalf(ell(p,beta)-palpha)-evalf(beta-alpha)*rhologp;
   od; # the computation of U_j is over
   if compU > compD then # computation of D_j
      lisU:=[op(lisD), seq(U[icom],icom=compD+1..compU)]; # lisU is the set U_j
      lisU:=sort(lisU,proc(a,b) evalb(a[1] < b[1]) end);
      lisD:=elaguage(lisU); # lisD is the set D_j after pruning U_j
      compD:=nops(lisD); # compD is the number of elements of D_j
      if imp >= 6 then 
         print(`p=`,p,`alpha=`,alpha,`number of terms of lists=`,
                compU,compD, `time=`, time()-ti); ti:=time();
      fi; 
   fi;
od;
# Now we consider the primes between x_2 and sqrt(x_1) (with exponent 1 in N); 
# thanks to the condition Bprime < B1, one does not divide and one multiplies 
# only one time by these primes. (cf. Proposition 4, (2))
p:=nextprime(p); # p is the smallest prime s.t. alpha=v_p(N)=1
rhologp:=rho*log(evalf(p)); p2mp:=p^2-p; 
b:=evalf(p2mp)-rhologp; # b is increasing on p cf. § 7.5
while b <= Bprime and p <= pk do
   compU:=0;
      # here, we construct the set U_j in two parts. The first part contains
      # the multiples of p. The second part is D_{j-1}. These two parts are 
      # sorted by merging and in the same time pruned by the procedure
      # "fusionelaguage".
   for icom from 1 to compD do
      if b+lisD[icom][3] <= Bprime then
         compU:=compU+1;
         U[compU]:=lisD[icom];
      fi;
   od;
   lisUp:=[seq(U[icom],icom=1..compU)];
   lisUp:=map(x->[x[1]*p,x[2]+p2mp,x[3]+b],lisUp); # the first part of U_j
   lisD:=fusionelaguage(compD,lisD,compU,lisUp);      
   compD:=nops(lisD);
   if imp >= 6 then 
      print(`p=`,p,`alpha=`,1,`nb of terms of lists=`,compU,compD, 
            `time=`, time()-ti); ti:=time();
   fi;
   p:=nextprime(p); rhologp:=rho*log(evalf(p)); p2mp:=p^2-p; 
   b:=evalf(p2mp)-rhologp;
od;
lisD;
end;




majben:=proc(liste,n,imp) local i,imax,omega,nmlofN,lm,Pom,Som,benom,
delta,ldelta,bendelta,benM,lM,pf,pfi,omegamin,omegamax,bi,imaxB;
global pk,lofN,x1,B,rho;
# this procedure computes B such that ben(g(n))+n-ell(g(n)) <= B
# by the method described in the paper, § 7.3. "liste" is a list of prefixes
# containing [1,0,0.]. First we compute the products
# of the primes close to p_k.
nmlofN:=n-lofN; 
# First the non-negative omega's
omega:=0; Pom[0]:=1; Som[0]:=0; benom[0]:=0.; 
pfi:=pk; pf:=evalf(pfi); lm:=liste[1][2];
# liste[1] is the smallest prefix of "liste". "omegamax" is the greatest 
# integer s.t. ell(N*delta_omegamax) <= n.
while lm +Som[omega] <= nmlofN do
   omega:=omega+1; 
   pfi:=nextprime(pfi); pf:=evalf(pfi); 
   Pom[omega]:=Pom[omega-1]*pfi; # Pom[omega]=p_{k+1}...p_{k+omega}
   Som[omega]:=Som[omega-1]+pfi; # Som[omega]=ell(Pom[omega])
   benom[omega]:=benom[omega-1]+pf-rho*log(pf); # =ben(N*Pom[omega])
od;
omegamax:=omega-1; # if delta=liste[1] is the smallest prefix, 
# "omegamax" is the greatest index s.t. ell(N*delta_omegamax) <= n.  
# Now, the negative omega's.
omega:=0; pfi:=pk; pf:=evalf(pfi); lm:=liste[nops(liste)][2];
while lm +Som[omega] > nmlofN and pf >= 1.0001*sqrt(x1) do 
   # the 2nd condition ensures that the prime factors of any prefix are smaller
   # than the prime factors of delta/delta_omega. if sqrt(x1) <= pf < 1.0001
   # sqrt(x1), we may forget a prefix, but no importance. 
   omega:=omega-1; 
   Pom[omega]:=Pom[omega+1]/pfi; # =1/(p_k p_{k-1}...p_{k+omega+1}
   Som[omega]:=Som[omega+1]-pfi; # Som[omega]=ell(Pom[omega])=S_omega
   benom[omega]:=benom[omega+1]+rho*log(pf)-pf; # =ben(N*Pom[omega])
   pfi:=prevprime(pfi); pf:=evalf(pfi); 
od;
omegamin:=omega; # the smallest omega to be used to construct delta_omega
imax:=cherchez(liste,nmlofN-Som[omegamin],1,nops(liste),2);
# imax is the index of the greatest prefix delta s.t. 
# ell(delta*N) <= ell(N/(p_k...p_{k+omegamin+1};
# imax is equal to nops(liste) unless we exit from the "while" loop
# with "pfi" small and lm+Som[omega] > nmlofN. This is very unlikely to 
# happen. If imax < nops(liste), we forget the elements of "liste" 
# of index > imax in the computation of the benefit.
# End of precomputation; now, computation of B
omega:=omegamax;
# the initial value of B is ben M +n-ell(M) with M=delta_omega where
# delta is the smallest prefix of "liste". It always exists provided that
# the prefix [1,0,0.] is in "liste".
B:=liste[1][3]+benom[omega]+evalf(nmlofN-liste[1][2]-Som[omega]); imaxB:=1;
for i from 2 to imax do
   delta:=liste[i][1]; # delta is the current prefix
   ldelta:=liste[i][2]; # ldelta=ell(delta*N)-ell(N)
   bendelta:=liste[i][3]; # bendelta=ben(delta*N);
   # computation of omega
   while ldelta+Som[omega] > nmlofN do omega:=omega-1; od;
   benM:=bendelta+benom[omega]; # by formula (6.5), benM=ben M.
   lM:=ldelta+Som[omega]; # lM=ell(M)-ell(N) 
   bi:=benM+evalf(nmlofN-lM);
     # bi is equal to ben M +n-ell(M) with M=delta_omega
   if bi < B then B:=bi; imaxB:=i; fi;
od;
if imp >= 8 then 
   print(`B=`,B, `imaxB=`,imaxB,liste[imaxB],`nmlofN=`,nmlofN);
fi;
B;
end;

shorten:=proc(liste,maj) local i,com,T;
# "liste" is a list of plain prefixes. Each term of "liste" is a list of 3 
# elements [g,lg,ben]. The terms of "liste" with ben > maj are crossed out.
com:=0;
for i from 1 to nops(liste) do
   if liste[i][3] <= maj then
      com:=com+1; T[com]:=liste[i];
   fi;
od;
[seq(T[icom],icom=1..com)];
end;


prefnorm:=proc(liste,n,imp) local t1dim,sqrtx1,di,i,ii,omega,nmlofN,lm,Pom,
Som,benom,pfi,pf,unmoinsrhosurt1,omegamin,omegamax,comPI,T,PI,lPI,
benPI,diPI,pi,lipn,pfin,sqrtx1aug;
global pk,lofN,x1,B,Baug,rho,t1;
# "liste" is the list of plain prefixes D(Baug) whose benefit is <= Baug. 
# We apply (7.36) which asserts that ell(N*pi) is at most
# n-Som[omega] and at least n-Som[omega]-di with di=B/(1-rho/t1).
# (replacing B by Baug and t1 by t1dim will eventually increase {curly N})
# For each possible omega, we look for the elements "pi" of "liste"
# satisfying the above inequalities.
# First we compute the sums S_omega almost in the same
# way than in procedure "majben". Each possible normalized prefix 
# is represented by a list of 6 elements 
# [PI,ell(PI*N)-ell(N),ben(PI*N),omega,p_{k+omega},pi]
t1:=fsolve(B+t11-rho*log(t11)=0,t11=0.9999*rho..1.0001*x1); # cf. § 7.6, (69)
t1dim:=0.9999*t1-0.0001; # to avoid round off errors
sqrtx1aug:=1.0001*sqrt(x1);
if t1dim <= sqrtx1aug then # very unlikely to happen
   error "t1dim=%1 is smaller than sqrt(x1)=%2; t1=%3",t1dim,sqrtx1,t1;
fi;
unmoinsrhosurt1:=1.-rho/t1dim;
di:=ceil(Baug/unmoinsrhosurt1);
nmlofN:=n-lofN;
# computation of omegamax
omega:=0; Pom[0]:=1; Som[0]:=0; benom[0]:=0.;  
pfi:=pk; pf:=evalf(pfi); pfin[0]:=pk; lm:=liste[1][2];
while lm +Som[omega] <= nmlofN do
   omega:=omega+1; pfi:=nextprime(pfi); 
   pfin[omega]:=pfi; pf:=evalf(pfi);# pfin[omega]=p_{k+omega}
   Pom[omega]:=Pom[omega-1]*pfi; # Pom[omega]=p_{k+1}...p_{k+omega}
   Som[omega]:=Som[omega-1]+pfi; # Som[omega]=ell(Pom[omega])=S_omega
   benom[omega]:=benom[omega-1]+pf-rho*log(pf); # =ben(N*Pom[omega])
od;
omegamax:=omega-1;  
# computation of omegamin
omega:=0; pfi:=pk; pf:=evalf(pfi); lm:=liste[nops(liste)][2];
while lm +Som[omega] > nmlofN -di and pf >= t1dim do # formula (77)    
   omega:=omega-1; 
   Pom[omega]:=Pom[omega+1]/pfi; # =1/(p_k p_{k-1}...p_{k+omega+1}
   Som[omega]:=Som[omega+1]-pfi; # Som[omega]=ell(Pom[omega])
   benom[omega]:=benom[omega+1]+rho*log(pf)-pf; # =ben(N*Pom[omega])
   pfi:=prevprime(pfi); pfin[omega]:=pfi; pf:=evalf(pfi);
od; 
omegamin:=omega;   
comPI:=0;
for omega from omegamin to omegamax do
   i:=cherchez(liste,nmlofN-Som[omega],1,nops(liste),2);
   ii:=i;
   while ii >= 1 and liste[ii][2] >= nmlofN-Som[omega]-di do
      lPI:=liste[ii][2]+Som[omega]; # lPI=ell(PI*N)-ell(N)
      benPI:=liste[ii][3]+benom[omega]; # benPI=ben(PI*N)
      diPI:=ceil((Baug-benPI)/unmoinsrhosurt1);
      if lPI >= nmlofN-diPI then # this condition comes from formula (77). 
         pi:=liste[ii][1];
         PI:=pi*Pom[omega]; # PI is a new normalized prefix
         comPI:=comPI+1; 
         T[comPI]:=[PI,lPI,benPI,omega,pfin[omega],pi];
         if pfin[omega+1]-(nmlofN-lPI) < sqrtx1aug then 
            error "there is a problem with the normalized prefix %1",T[comPI]; 
         fi; # very unlikely to happen, cf. (80) and § 7.8.
      fi;
      ii:=ii-1;
   od;
od;
lipn:=[seq(T[icom],icom=1..comPI)];
if imp >= 8 then 
   print(`list of normalized prefixes=`,lipn,`card=`,nops(lipn));
fi;
lipn;
end;

fight:=proc(liste,pk,n,imp) local nmlofN,sui,omax,omin,om,pfin,i,m,p,PI,T,TT,
lisT,lisTord,mini,maxi,bou,compt,lisTT;
global lofN;
# this procedure eliminates some normalized prefixes as explained in § 7.9
# "liste" is a list of normalized prefixes. Each prefix is represented
# by a list of 6 elements [PI,ell(PI*N)-ell(N),ben(PI*N),omega,pfin[omega],pi]
if nops(liste) <= 1 then return liste; fi;
nmlofN:=n-lofN;
# construction of pfin
sui:=op(map(x->x[4],liste)); omin:=min(sui); omax:=max(sui);
pfin[0]:=pk;
for om from 0 to omax do pfin[om+1]:=nextprime(pfin[om]); od;
for om from 0 by -1 to omin+2  do pfin[om-1]:=prevprime(pfin[om]); od;
# ordering of "liste" on the increasing values of "maxi"
for i from 1 to nops(liste) do
   PI:=liste[i][1]; om:=liste[i][4]; p:=pfin[om+1]; m:=nmlofN-liste[i][2];
   maxi:=PI*p/(p-m); # this is the upper bound g''(PI,n) of § 7.9 and (90)
   T[i]:=[op(liste[i]),maxi];
od;
lisT:=[seq(T[ii],ii=1..nops(liste))];
lisTord:=sort(lisT,proc(a,b) evalb(a[7] > b[7]) end);
PI:=lisTord[1][1]; om:=lisTord[1][4]; p:=pfin[om+1]; m:=nmlofN-lisTord[1][2];
mini:=PI*p/nextprime(p-m-1);
bou:=true;
compt:=1; TT[compt]:=lisTord[1]; # the 1st element of lisTord cannot be killed
for i from 2 to nops(liste) while bou do
   if lisTord[i][7] < mini then
      bou:=false;
   else
      compt:=compt+1; TT[compt]:=lisTord[i];
      PI:=lisTord[i][1]; om:=lisTord[i][4]; p:=pfin[om+1]; 
      m:=nmlofN-lisTord[i][2];
      mini:=max(mini,PI*p/nextprime(p-m-1));
   fi;
od;
lisTT:=[seq(TT[icompt],icompt=1..compt)];
lisTT:=map(x->[x[ii] $ ii=1..6],lisTT);
if imp >= 8 then 
   print(`ordered list of normalized prefixes=`,lisTord,evalf(lisTord),
          `after fight`,lisTT);
fi;
lisTT;
end;

mj:=proc(j,K,r) local a; global S;
# This short procedure calculates m_j(P_r) defined in the article, (103).
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, § 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. (8.7))
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 (107)
         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;

listediffprem:=[1,3,2,5,4,11,6,29,8,97,14,127,18,541,20,907,22,1151,34,1361,
36,9587,44,15727,52,19661,72,31469,86,156007,96,360749,112,370373,114,492227,
118,1349651,132,1357333,148,2010881,154,4652507,180,17051887,210,20831533,
220,47326913,222,122164969,234,189695893,248,191913031,250,387096383,
282,436273291,288,1294268779,292,1453168433,320,2300942869,
336,3842611109,354,4302407713,382,10726905041,384,20678048681,
394,22367085353,456,25056082543,464,42652618807,468,127976335139,
474,182226896713,486,241160624629,490,297501076289,500,303371455741,
514,304599509051,516,416608696337,532,461690510543,534,614487454057,
540,738832928467,582,1346294311331,588,1408695494197,602,1968188557063,
652,2614941711251];

# The number following 468 is wrong in the table of Riesel's book 

diffprem:=[];
for i from 1 by 2 to nops(listediffprem) do
   diffprem:=[op(diffprem),[listediffprem[i],listediffprem[i+1]]];
od:

deltaprem:=proc(x) local a; global diffprem;
# returns the largest gap between two consecutive primes not exceeding x
if not type(x,integer) then error "bad choice of x=%1",x fi;
a:=cherchez(diffprem,x,1,54,2);
diffprem[a][1];
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 > pkp1-3 or m < 0 or not type(m,integer) then error "m=%1 is wrong",m; fi;
if m < pkp1-pk then 
   return [1,[`m is smaller than p_{k+1}-p_k`]];
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;