#VERSION OF April 13, 1995
#First version: March 3, 1995; A package that implements Sylvester's
#enhancement of Tchebychev's method for finidng upper and lower
#bounds for psi(x)/x where psi(x) is the sum of log(p) over all
#prime powers p^m <=x
print(`Version of April 13, 1995`):
print(`SYLVESTER: A Maple package written by Doron Zeilberger`):
print(`For general help, and a list of functions, type "ezra()"`):
print(`For specific help, type ezra(procedure_name)`):
 
 
ezra:=proc()
if args=NULL then
print(`This is SYLVESTER, a Maple package implementing Sylvester's`):
print( `obsolete(?) method given, in final form in his 1892 paper`):
print(` In CW IV 687-731`):  
print(`Contains the following procedures: stigser,parsepos,parseneg,zugotpos`):
print(`'zugotneg' ,'matri' ,'Iter' `):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`Iter` then
print(` Iter(scheme,A_l,A_u) inputs a scheme  and the current`):
print(`lower and upper bounds for psi(x)/x, x large `):
print(`and outputs improved ones`):
print(`see Sylvester's CW IV, p. 717`):
fi:
 
if nops([args])=1 and op(1,[args])=`matri` then
print(` matri(scheme,cutof) inputs a scheme  and a numerical cutoff`):
print(` finds the 2 by 2 matrix [[a,b],[c,d]] such that `):
print(`u[i+1]=M+a*u[i]-b*v[i] ; v[i+1]=M+c*v[i]-d*u[i]`):
print(`see Sylvester's CW IV, p. 713`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`stigser` then
print(`stigser(scheme) gives the stignatic series corresponding to the`):
print(`scheme given as a list of lists [[a_1,...,a_r];[b_1,...,b_s]]`):
print(`Corresponding to forming Sum_{i=0}^r T(x/a_i)-Sum_{i=0}^s T(x/b_i)`):
print(`where T(x)=log(x!); The stigmatic series is the coding of the`):
print(`periodic expression in psi(x). It is given as a list,each `):
print(`entry of which is a list of two elements [a,b] where b is +1 or -1`):
print(`and it represents the term b*psi(x/a), in the first period`):
fi:
 
if nops([args])=1 and op(1,[args])=`parsepos` then
print(`parsepos(stig) inputs a stimatics series, in the form that is`):
print(` outputted from stigser(), and outputs the decomposition into`):
print(`an initial segment, folowed by "Irred. ballot paths", i.e. hills`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`parseneg` then
print(`parseneg(stig) inputs a stimatics series, in the form that is`):
print(` outputted from stigser(), and outputs the decomposition into`):
print(`an initial segment, folowed by "inverted Irr.ballot paths",valleies`):
 
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`zugotpos` then
print(`zugotpos(IP): given a list of 2-element lists [n,e] such that`):
print(`the e's form an "irreducible upper ballot sequence", pairs up`):
print(` like parentheses the corresponding pairs. `):
print(`The output is given as a list of pairs of the n s`):
fi:
 
if nops([args])=1 and op(1,[args])=`zugotneg` then
print(`zugotneg(IP): given a list of 2-element lists [n,e] such that`):
print(`the e's form an "irreducible lower ballot sequence", pairs up`):
print(` like parentheses the corresponding pairs. `):
print(`The output is given as a list of pairs of the n s`):
fi:
 
end:
 
 
#pol is a poly. ra is a rational function
 
ple:=proc(pol,ra,x)
local gu,i,pol1,gun,gud:
 
pol1:=expand(pol):
 
gu:=0:
 
for i from 1 to degree(pol1,x) do
gu:=gu+coeff(pol1,x,i)*subs(x=x^i,ra):
gu:=normal(gu):
 
od:
 
gun:=expand(numer(gu)*(1-x)):
gud:=expand(denom(gu)*(1-x)):
gun:=taylor(gun,x=0,degree(gun)+1):
gun,gud:
end:
 
#stigesr: Given a scheme [[a_1, .., a_r],[b_1,..., b_s]] finds
#the `Stigmatic Series'
 
stigser:=proc(scheme)
local lu,schpos,schneg,i,j,pol, ra,gu,gun,gud,resh,x :
 
schpos:=op(1,scheme):
schneg:=op(2,scheme):
ra:=x/(1-x):
pol:=0:
 
 
for i from 1 to nops(schpos) do
 pol:=pol+x^op(i,schpos):
od:
 
 
for i from 1 to nops(schneg) do
 pol:=pol-x^op(i,schneg):
od:
 
gu:=ple(pol,ra,x):
 
gun:=gu[1]:
gud:=gu[2]:
 
resh:=[]:
 
for i from 1 to degree(gud,x) do
 lu:=coeff(gun,x,i):
 
 if lu>0 then
    for j from 1 to lu do
       resh:=[op(resh),[i,1]]:
    od:
 elif lu<0 then
       for j from 1 to -lu do
        resh:=[op(resh),[i,-1]]:
       od:
 fi: 
od:
 
resh:
 
end:
 
 
#parsepos(stig) intputs a stigmatics series given as a list of
#2-element list [a,e] where e=+-1 and a in an integer that corresponds
# to e*psi(x/a), and divides them into groups, just like in
#Sylvester's CW  IV, p. 707
 
parsepos:=proc(stig)
local makh,i,gu,mu,parsum,hakhi,orekh,makom:
#makh is the `makhzor, the period; parsum is the partial sum, hakhi is
#minimum of the partial sums ; makom is the place where the recod takes
#place
 
orekh:=nops(stig):
makh:=op(1,op(orekh,stig)):
parsum:=0:
hakhi:=0:
makom:=0:
for i from 1 to orekh do
parsum:=parsum+op(2,op(i,stig)):
 
if parsum<hakhi then
   hakhi:=parsum:
   makom:=i:
 fi:
od:
 
mu:=[[op(1..makom,stig)]]:
 
parsum:=0:
 
lu:=[]:
 
 
for i from makom+1 to orekh do
 parsum:=parsum+op(2,op(i,stig)):
  lu:=[op(lu),op(i,stig)]:
 
  if parsum=0 then
     mu:=[op(mu),lu]:
     lu:=[]:
  fi:
od:
 
 
 
for i from 1 to makom do
 parsum:=parsum+op(2,op(i,stig)):
  lu:=[op(lu),[op(1,op(i,stig))+makh,op(2,op(i,stig))]]:
 
  if parsum=0 then
     mu:=[op(mu),lu]:
     lu:=[]:
  fi:
od:
 
if op(1,mu)=[] then
 mu:=[op(2..nops(mu),mu)]:
fi:
mu:
   
end:
 
 
 
#parseneg(stig) intputs a stigmatics series given as a list of
#2-element list [a,e] where e=+-1 and a in an integer that corresponds
# to e*psi(x/a), and divides them into groups, just like in
#Sylvester's CW  IV, p. 707; the second way, where it is negatively oriented
 
parseneg:=proc(stig)
local makh,i,i1,lu,lu1,mu,parsum,hakhi,orekh,makom:
#makh is the `makhzor, the period; parsum is the partial sum, hakhi is
#minimum of the partial sums ; makom is the place where the recod takes
#place
 
orekh:=nops(stig):
makh:=op(1,op(orekh,stig)):
parsum:=0:
hakhi:=0:
makom:=0:
for i from 1 to orekh do
parsum:=parsum+op(2,op(i,stig)):
 
if parsum>hakhi then
   hakhi:=parsum:
   makom:=i:
 fi:
od:
 
mu:=[[op(1..makom,stig)]]:
 
parsum:=0:
 
lu:=[]:
 
 
for i1 from makom+1 to orekh do
 parsum:=parsum+op(2,op(i1,stig)):
  lu:=[op(lu),op(i1,stig)]:
 
  if parsum=0 then
     mu:=[op(mu),lu]:
     lu:=[]:
  fi:
od:
 
 
 
for i from 1 to makom do
 parsum:=parsum+op(2,op(i,stig)):
  lu:=[op(lu),[op(1,op(i,stig))+makh,op(2,op(i,stig))]]:
  if parsum=0 then
    
    lu1:=[]:
     for i1 from 1 to nops(lu) do
  
       if op(1,op(i1,lu))<makh then
               lu1:=[op(lu1),[op(1,op(i1,lu))+makh,op(2,op(i1,lu))]]:
         else
               lu1:=[op(lu1),[op(1,op(i1,lu)),op(2,op(i1,lu))]]:
       fi:
 
    od:
     mu:=[op(mu),lu1]:
     lu:=[]:
  fi:
od:
 
mu:
   
end:
 
 
#zugotpos(IP): given a list of 2-element lists [n,e] such that
#the e's form an `irreducible upper ballot sequence', pairs up like parentheses
#the corresponding pairs. The output is given as a list of pairs of the n's
 
zugotpos:=proc(IP)
local rishon, akharon,gu,i,IPk,lu,lu1:
 
if nops(IP)=0 then
 RETURN([]):
fi:
 
rishon:=op(1,IP):
akharon:=op(nops(IP),IP):
 
gu:=[[op(1,rishon),op(1,akharon)]]:
 
IPk:=[op(2..nops(IP)-1,IP)]:
 
if nops(IPk)=0 then
 RETURN(gu):
fi:
 
lu:=parsepos(IPk):
 
for i from 1 to nops(lu) do
 lu1:=op(i,lu):
 gu:=[op(gu),op(zugotpos(lu1))]
od:
 
gu:
end:
 
 
#zugotneg(IP): given a list of 2-element lists [n,e] such that
#the e's form an `irreducible upper ballot sequence', pairs up like parentheses
#the corresponding pairs. The output is given as a list of pairs of the n's
 
zugotneg:=proc(IP)
local rishon, akharon,gu,i,IPk,lu,lu1:
 
if nops(IP)=0 then
 RETURN([]):
fi:
 
rishon:=op(1,IP):
akharon:=op(nops(IP),IP):
 
gu:=[[op(1,rishon),op(1,akharon)]]:
 
IPk:=[op(2..nops(IP)-1,IP)]:
 
if nops(IPk)=0 then
 RETURN(gu):
fi:
 
lu:=parseneg(IPk):
 
for i from 1 to nops(lu) do
 lu1:=op(i,lu):
 gu:=[op(gu),op(zugotneg(lu1))]
od:
 
gu:
end:
 
  
#matri: Given a scheme and cutof, finds the 2 by 2 matrix
#[[a,b],[d,c]] such that (see Sylvester's CW IV, p. 713)
#u[i+1]=M+a*v[i]-b*u[i] ; v[i+1]=M+c*v[i]-d*u[i]
 
matri:=proc(scheme,cutof)
local stigser1,parsepos1,parseneg1,a,b,c,d,mu,mu1,mu11,zugotpos1,zugotneg1,
i1,zug,makhzor,i,j,lista,listb,listc,listd,schel,scher,M,E,F,ku1,rho,lu:
 
M:=0:
 
schel:=op(1,scheme):
scher:=op(2,scheme):
 
 
for i from 1 to nops(schel) do
 M:=M-log(op(i,schel))/op(i,schel):
od:
 
 
for i from 1 to nops(scher) do
 M:=M+log(op(i,scher))/op(i,scher):
od:
 
print(`The stigmatic   multiplier is`):
print(M):
print(`Which is roughly`,evalf(M)):
 
 
a:=0:
b:=0:
c:=0:
d:=0:
stigser1:=stigser(scheme):
makhzor:=op(1,op(nops(stigser1),stigser1)):
parseneg1:=parseneg(stigser1):
parsepos1:=parsepos(stigser1):
 
mu:=op(1,parsepos1):
 
mu1:=op(1,mu):
 
if not (op(1,mu1)=1 and op(2,mu1)=1 ) then
  ERROR(`Something went wrong`):
fi:
 
lista:=[]:listb:=[]:listc:=[]:listd:=[]:
for i from 2 to nops(mu) do
 mu11:=op(i,mu):
 
 if op(2,mu11)=-1 then
   a:=a+1/op(1,mu11):
    lista:=[op(lista),1/op(1,mu11)]:
 elif op(2,mu11)=1 then
    b:=b+1/op(1,mu11):
      listb:=[op(listb),1/op(1,mu11)]:
 else
   ERROR(`Something is wrong`):
 fi:
od:
 
for i from 2 to nops(parsepos1) do
 mu:=op(i,parsepos1):
 zugotpos1:=zugotpos(mu);
 
 for j from 1 to nops(zugotpos1) do
   zug:=op(j,zugotpos1):
 
   if op(2,zug)/op(1,zug)>cutof then
    a:=a+1/op(2,zug):
    lista:=[op(lista),1/op(2,zug)]:
    b:=b+1/op(1,zug):
    listb:=[op(listb),1/op(1,zug)]:
 
      if (op(2,zug)+makhzor)/(op(1,zug)+makhzor) >cutof then
       for i1 from 1 while (op(2,zug)+i1*makhzor)/(op(1,zug)+i1*makhzor)>
                                                              cutof do
               a:=a+1/(op(2,zug)+i1*makhzor):
       lista:=[op(lista),1/(op(2,zug)+i1*makhzor)]:
               b:=b+1/(op(1,zug)+i1*makhzor):
       listb:=[op(listb),1/(op(1,zug)+i1*makhzor)]:
       od:
      fi:
  fi:
 od:
od:
 
mu:=op(1,parseneg1):
 
mu1:=op(1,mu):
 
if not (op(1,mu1)=1 and op(2,mu1)=1 ) then
  ERROR(`Something went wrong`):
fi:
 
for i from 2 to nops(mu) do
 mu11:=op(i,mu):
 
 if op(2,mu11)=-1 then
   c:=c+1/op(1,mu11):
   listc:=[op(listc),1/op(1,mu11)]:
 elif op(2,mu11)=1 then
    d:=d+1/op(1,mu11):
   listd:=[op(listd),1/op(1,mu11)]:
 else
   ERROR(`Something is wrong`):
 fi:
od:
 
for i from 2 to nops(parseneg1) do
 mu:=op(i,parseneg1):
 zugotneg1:=zugotneg(mu);
 
 for j from 1 to nops(zugotneg1) do
   zug:=op(j,zugotneg1):
 
   if op(2,zug)/op(1,zug)>cutof then
    c:=c+1/op(1,zug):
   listc:=[op(listc),1/op(1,zug)]:
    d:=d+1/op(2,zug):
   listd:=[op(listd),1/op(2,zug)]:
 
 
      if (op(2,zug)+makhzor)/(op(1,zug)+makhzor) >cutof then
       for i1 from 1 while (op(2,zug)+i1*makhzor)/(op(1,zug)+i1*makhzor)>
                                                              cutof do
               c:=c+1/(op(1,zug)+i1*makhzor):
         listc:=[op(listc),1/(op(1,zug)+i1*makhzor)]:
               d:=d+1/(op(2,zug)+i1*makhzor):
         listd:=[op(listd),1/(op(2,zug)+i1*makhzor)]:
       od:
      fi:
  fi:
 od:
od:
 
ku1:=fsolve(rho^2-(a+c)*rho+(a*c-b*d)=0,rho):
print(`The latent roots are`, ku1):
 
E:=(1-b-c)/((1-a)*(1-c)-b*d)*M:
F:=(1-a-d)/((1-a)*(1-c)-b*d)*M:
print(`It follows from the Sylvester-Tchebycheff method that`):
print(`If we already know that there exist A_u and A_l such that`):
print(`for large x, psi(x)<=A_u x and psi(x)>=A_l x, with A_u/A_l larger `):
print(`than`,cutof):
print(`Then the new A_u is`):
print(E,`=`,evalf(E),`...`):
 
print(`and the new A_l is`):
print(F,`=`,evalf(F),`...`):
 
[[a,b],[d,c]],lista,listb,listd,listc,evalf(E),evalf(F),evalf(E/F):
 
 
end:
    
 
#Iter: Given a scheme and current best known upper and lower 
#bounds A_u and A_l finds new, possibly better ones
 
Iter:=proc(scheme,A_l,A_u)
local stigser1,parsepos1,parseneg1,a,b,c,d,mu,mu1,mu11,zugotpos1,zugotneg1,
i1,zug,makhzor,i,j,lista,listb,listc,listd,schel,scher,M,E,F,cutof,lu:
 
cutof:=A_u/A_l:
M:=0:
 
schel:=op(1,scheme):
scher:=op(2,scheme):
 
 
for i from 1 to nops(schel) do
 M:=M-log(op(i,schel))/op(i,schel):
od:
 
 
for i from 1 to nops(scher) do
 M:=M+log(op(i,scher))/op(i,scher):
od:
 
print(`The schematic multiplier is`):
print(M):
print(`Which is roughly`,evalf(M)):
 
 
a:=0:
b:=0:
c:=0:
d:=0:
stigser1:=stigser(scheme):
makhzor:=op(1,op(nops(stigser1),stigser1)):
parseneg1:=parseneg(stigser1):
parsepos1:=parsepos(stigser1):
mu:=op(1,parsepos1):
mu1:=op(1,mu):
 
if not (op(1,mu1)=1 and op(2,mu1)=1 ) then
  ERROR(`Something went wrong`):
fi:
 
lista:=[]:listb:=[]:listc:=[]:listd:=[]:
for i from 2 to nops(mu) do
 mu11:=op(i,mu):
 
 if op(2,mu11)=-1 then
   a:=a+1/op(1,mu11):
    lista:=[op(lista),1/op(1,mu11)]:
 elif op(2,mu11)=1 then
    b:=b+1/op(1,mu11):
      listb:=[op(listb),1/op(1,mu11)]:
 else
   ERROR(`Something is wrong`):
 fi:
od:
 
 
for i from 2 to nops(parsepos1) do
 
 mu:=op(i,parsepos1):
 zugotpos1:=zugotpos(mu);
 
 
 for j from 1 to nops(zugotpos1) do
   zug:=op(j,zugotpos1):
   if op(2,zug)/op(1,zug)>evalf(cutof) then
    a:=a+1/op(2,zug):
 
 
    lista:=[op(lista),1/op(2,zug)]:
    b:=b+1/op(1,zug):
    listb:=[op(listb),1/op(1,zug)]:
 
      if (op(2,zug)+makhzor)/(op(1,zug)+makhzor) >evalf(cutof) then
       for i1 from 1 while (op(2,zug)+i1*makhzor)/(op(1,zug)+i1*makhzor)>
                                                      evalf(cutof) do
               a:=a+1/(op(2,zug)+i1*makhzor):
       lista:=[op(lista),1/(op(2,zug)+i1*makhzor)]:
               b:=b+1/(op(1,zug)+i1*makhzor):
       listb:=[op(listb),1/(op(1,zug)+i1*makhzor)]:
       od:
      fi:
  fi:
 od:
od:
 
mu:=op(1,parseneg1):
 
mu1:=op(1,mu):
 
if not (op(1,mu1)=1 and op(2,mu1)=1 ) then
  ERROR(`Something went wrong`):
fi:
 
for i from 2 to nops(mu) do
 mu11:=op(i,mu):
 
 if op(2,mu11)=-1 then
   c:=c+1/op(1,mu11):
   listc:=[op(listc),1/op(1,mu11)]:
 elif op(2,mu11)=1 then
    d:=d+1/op(1,mu11):
   listd:=[op(listd),1/op(1,mu11)]:
 else
   ERROR(`Something is wrong`):
 fi:
od:
 
for i from 2 to nops(parseneg1) do
 mu:=op(i,parseneg1):
 zugotneg1:=zugotneg(mu);
 
 for j from 1 to nops(zugotneg1) do
   zug:=op(j,zugotneg1):
 
   if op(2,zug)/op(1,zug)>evalf(cutof) then
    c:=c+1/op(1,zug):
   listc:=[op(listc),1/op(1,zug)]:
    d:=d+1/op(2,zug):
   listd:=[op(listd),1/op(2,zug)]:
 
 
      if (op(2,zug)+makhzor)/(op(1,zug)+makhzor) >evalf(cutof) then
       for i1 from 1 while (op(2,zug)+i1*makhzor)/(op(1,zug)+i1*makhzor)>
                                                              evalf(cutof) do
               c:=c+1/(op(1,zug)+i1*makhzor):
         listc:=[op(listc),1/(op(1,zug)+i1*makhzor)]:
               d:=d+1/(op(2,zug)+i1*makhzor):
         listd:=[op(listd),1/(op(2,zug)+i1*makhzor)]:
       od:
      fi:
  fi:
 od:
od:
 
E:=M+a*A_u-b*A_l:
F:=M+c*A_l-d*A_u:
 
print(`It follows from the Sylvester-Tchebycheff method that`):
print(`by taking the scheme`,scheme):
print(`if we already know that there exist A_u and A_l such that`):
print(`for large x, psi(x)<=A_u x and psi(x)>=A_l x`):
print(`with A_u=`,evalf(A_u),`and A_l=`,evalf(A_l)):
print(`Then the new A_u can be taken  to be `):
print(E,`=`,evalf(E),`...`):
 
 
if evalf(E)<evalf(A_u) then
 print(`Since this is a better upper bound, we adopt it`):
else
 print(`Since this is a worse upper bound, we do not dopt it`):
 E:=A_u:
fi:
 
 
print(`and the new A_v can be taken  to be `):
print(F,`=`,evalf(F),`...`):
 
 
if evalf(F)>evalf(A_l) then
 print(`Since this is a better lower bound, we adopt it`):
else
 print(`Since this is a worse lower bound, we do not adopt it`):
 F:=A_l:
fi:
 
 
F,E,[[a,b],[d,c]],lista,listb,listd,listc,evalf(F),evalf(E):
 
 
end: