######################################################################
##PiDay2019.txt: Save this file as  PiDay2019.txt                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read PiDay2019.txt<Enter>                          #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: March 14, 2019
 
print(`Created:  March 14, 2019`):
print(` This is PiDay2019.txt `):
print(`It is one of the packages that accompany the talk `):
print(` What is Pi and What it is Not`):
print(`by  Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AsySe `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Asy, AsyBin, PtoSeq, RecBin, SeferBin `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Asy then
print(`Asy(rec,K,M): inputs a P-recursive sequence rec=[INI, [n, [P0,...,Pk]]], pos. integer K and a a large pos. integer M`):
print(`and outputs the asympotic expansion to order K such L[n]=C*mu^n*n^theta (1+a1/n+...+aK/n^K).`):
print(` M is used to estimate the constant C. Try: `):
print(`Asy([[1/2],[n,[2*n,1-2*n]]],4,1000);`):

elif nops([args])=1 and op(1,[args])=AsyBin then
print(`AsyBin(F,k,n,K,M): the asympototic expression to order K of the binomial coefficients sum Sum(F(n,k),k=0..n) with the natural limits.`):
print(`Try: `):
print(`AsyBin(binomial(n,k)^3,k,n,3,1000);`):

elif nops([args])=1 and op(1,[args])=AsySe then
print(`AsySe(rec,K,M): Like Asy(rec,K,M) but the output is separated into [F,mu,theta,f]. Try: `):
print(`AsySe([[1/2],[n,[2*n,1-2*n]]],4,1000);`):



elif nops([args])=1 and op(1,[args])=PtoSeq then
print(`PtoSeq(P,N): inputs a P-recursive sequence in the format [INI,[n,[P0,.., Pk]] where INI is the list of initial values for n=1, ...,k and`):
print(`P0,..., Pk are polynomials in n, representing the coefficients of the linear recurrence with positive coefficients`):
print(`P0(n)*x(n)+ ...+ Pk(n)*x(n-k)=0 with those initial conditions, and outputs the first N terms of the sequence`):
print(`(written by Victoria Chayes)`):
print(`Try:`):
print(`PtoSeq([[1],[n,[1,-n]]],10)`);


elif nops([args])=1 and op(1,[args])=RecBin then
print(`RecBin(F,k,n): the P-recursive represenation of Sum(F(n,k),k=0..n) where F(n,k) is a bi-hypergoemetric sum`):
print(`RecBin(binomial(n,k)^3,k,n);`):
print(``):

elif nops([args])=1 and op(1,[args])=SeferBin then
print(`SeferBin(G,n,k,K,Ga): inputs a positive integer G, symbols n and k, pos. integer K, and a large pos. integer Ga, outputs`):
print(`an article about the asympotics ofsum of the r-th power of the binomial coefficients . Try:`):
print(`SeferBin(4,n,k,6,5000);`):

else
print(`There is no ezra for`,args):
fi:
 
end:

Z:=SumTools[Hypergeometric][Zeilberger]:



#PtoSeq(P,N): inputs a P-recursive sequence in the format [INI,[n,[P0,.., Pk]] where INI is the list of initial values for n=1, ...,k and
#P0,..., Pk are polynomials in n, representing the coefficients of the linear recurrence with positive coefficients
#P0(n)*x(n)+ ...+ Pk(n)*x(n-k)=0 with those initial conditions, and outputs the first N terms of the sequence
#(stolen from Victoria Chayes' hw12.txt)
PtoSeq:=proc(P,N) local L,n, nt, i, l,k,P0:
L:=P[1]: n:=P[2][1]: l:=nops(P[2][2])-1: k:=nops(L)+1: P0:=-1/P[2][2][1]:
while k<=N do   nt:=subs(n=k,P0)*add(subs(n=k, P[2][2][i+1])*L[-i],i=1..l):   L:=[op(L),nt]:  k:=k+1:od:
L:
end:

#The max. in absolute value
MaxAbs:=proc(L) local i,cha,rec:
cha:=L[1]:
rec:=evalf(abs(cha)):

for i from 2 to nops(L) do
  if evalf(abs(L[i]))>rec then
    cha:=L[i]:
    rec:=evalf(abs(cha)):
 fi:
od:
cha:
end:

 #Asy(rec,K,M): inputs a P-recursive sequence rec=[INI, [n, [P0,...,Pk]]], pos. integer K and a a large pos. integer M
 #and outputs the asympotic expansion to order K such L[n]=C*mu^n*n^theta (1+a1/n+...+aK/n^K).
 #M is used to estimate the constant C. Try:
 #Asy([[1/2],[n,[2*n,1-2*n]],4,1000);
Asy:=proc(rec,K,M) local n,rec0,d,i,rec00,mu,p,mu0,LS,x,theta,a,f,gu,eq,var,C:
 n:=rec[2][1]:
 rec0:=rec[2][2]:
  d:=max(seq(degree(rec0[i],n),i=1..nops(rec0))):
 rec00:=[seq(coeff(rec0[i],n,d),i=1..nops(rec0))]:
 p:=numer(add(rec00[i+1]/mu^i,i=0..nops(rec00)-1)):
 mu0:=MaxAbs([solve(p,mu)]):
  if not evalf(mu0)> 0 then
       RETURN(FAIL):
   fi:
 rec0:=[seq(rec0[i+1]/mu0^i,i=0..nops(rec0)-1)]:
 LS:=add(rec0[i+1]*(1-i/n)^theta/n^d,i=0..nops(rec0)-1):
 LS:=normal(subs(n=1/x,LS)):
 LS:=taylor(LS,x=0,3);
  if coeff(LS,x,0)<>0 then
     RETURN(FAIL):
   fi:
 LS:=coeff(LS,x,1):

theta:=solve(LS,theta):
 f:=1+add(a[i]/n^i,i=1..K):
   var:={seq(a[i],i=1..K)}:
 LS:=add(rec0[i+1]*(1-i/n)^theta/n^d*subs(n=n-i,f),i=0..nops(rec0)-1):
 LS:=normal(subs(n=1/x,LS)):
 gu:=taylor(LS,x=0,K+2):
  eq:={seq(coeff(gu,x,i),i=2..K+1)}:
 var:=solve(eq,var):
  f:=subs(var,f):
f:=mu0^n*n^theta*f:

C:=evalf(PtoSeq(rec,M)[M]/ subs(n=M,f)):
C:=identify(C):

  C*f:
 end:

#RecBin(F,k,n): the P-recursive represenation of Sum(F(n,k),k=0..n) where F(n,k) is a bi-hypergoemetric sum
#RecBin(binomial(n,k)^3,k,n);
RecBin:=proc(F,k,n) local ope,N,ord,rec,INI,i,k1,n1:
if eval(subs({n=11,k=-1},F))<>0 or eval(subs({n=11,k=12},F))<>0  then
 print(`bad input`):
 RETURN(FAIL):
fi:

ope:=Z(F,n,k,N)[1]:
ord:=degree(ope,N):
ope:=expand(subs(n=n-ord,ope)/N^ord):
ope:=add(factor(coeff(ope,N,-i))*N^(-i),i=0..ord):

INI:=[seq(add(eval(subs({n=n1,k=k1},F)), k1=0..n1),n1=1..ord)]:
rec:=[INI,[n,[seq(coeff(ope,N,-i),i=0..ord)]]]:
rec:
end:




#AsyBin(F,k,n,K,M): the asympototic expression to order K of the binomial coefficients sum Sum(F(n,k),k=0..n) with the natural limits.
#Try:
#AsyBin(binomial(n,k)^3,k,n,K,M);
AsyBin:=proc(F,k,n,K,M) local ope,N,ord,rec,INI,i,k1,n1:

rec:=RecBin(F,k,n):
Asy(rec,K,M):
end:


 #AsySe(rec,K,M): Like Asy(rec,K,M) but the output is split into [C,mu,theta,f]. Try:
 #AsySe([[1/2],[n,[2*n,1-2*n]],4,1000);
AsySe:=proc(rec,K,M) local n,rec0,d,i,rec00,mu,p,mu0,LS,x,theta,a,f,gu,eq,var,C,f1:
 n:=rec[2][1]:
 rec0:=rec[2][2]:
  d:=max(seq(degree(rec0[i],n),i=1..nops(rec0))):
 rec00:=[seq(coeff(rec0[i],n,d),i=1..nops(rec0))]:
 p:=numer(add(rec00[i+1]/mu^i,i=0..nops(rec00)-1)):
 mu0:=MaxAbs([solve(p,mu)]):
  if not evalf(mu0)> 0 then
       RETURN(FAIL):
   fi:
 rec0:=[seq(rec0[i+1]/mu0^i,i=0..nops(rec0)-1)]:
 LS:=add(rec0[i+1]*(1-i/n)^theta/n^d,i=0..nops(rec0)-1):
 LS:=normal(subs(n=1/x,LS)):
 LS:=taylor(LS,x=0,3);
  if coeff(LS,x,0)<>0 then
     RETURN(FAIL):
   fi:
 LS:=coeff(LS,x,1):

theta:=solve(LS,theta):
 f:=1+add(a[i]/n^i,i=1..K):
   var:={seq(a[i],i=1..K)}:
 LS:=add(rec0[i+1]*(1-i/n)^theta/n^d*subs(n=n-i,f),i=0..nops(rec0)-1):
 LS:=normal(subs(n=1/x,LS)):
 gu:=taylor(LS,x=0,K+2):
  eq:={seq(coeff(gu,x,i),i=2..K+1)}:
 var:=solve(eq,var):
  f:=subs(var,f):
f1:=mu0^n*n^theta*f:

C:=evalf(PtoSeq(rec,M)[M]/ subs(n=M,f1)):
C:=identify(C):

[C,mu0,theta,f]:
 end:

#SeferBin(G,n,k,K,Ga): inputs a positive integer G, symbols n and k, pos. integer K, and a large pos. integer Ga, outputs
#an article about the asympotics ofsum of the r-th power of the binomial coefficients . Try:
#SeferBin(4,n,k,6,5000);

SeferBin:=proc(G,n,k,K,Ga) local g,A,rec,INI,rec1,as1,i:

Digits:=20:

print(`Recurrences and Asymptotic for the Sum of powers of the Rows of Pascal's Triangle`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
for g from 1 to G do
print(``):
print(`Theorem No. `, g, `Let `):
print(``):
print(A(n)=Sum(binomial(n,k)^g,k=0..n) ):
print(``):
rec:=RecBin(binomial(n,k)^g,k,n):


print(``):
INI:=rec[1]:
rec1:=rec[2][2]:
print(`Then A(n) satisfies the recurrence `):
print(``):
print(add(rec1[i+1]*A(n-i),i=0..nops(rec1)-1) =0 ):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(A(i)=INI[i],i=1..nops(INI))):
print(``):
print(`In Maple notation:`):
print(``):
lprint(add(rec1[i+1]*A(n-i),i=0..nops(rec1)-1) =0 ):
print(``):
print(`The first 30 terms are`):
print(``):
print(PtoSeq(rec,30)):
print(``):
as1:=Asy(rec,K,Ga):
print(``):
print(`A(n) is asymptotic to `):
print(``):
print(as1):
print(``):
print(`In Maple notation:`):
print(``):
lprint(as1):
od:

print(``):

print(`----------------------------------`):
print(``):
print(`This took`, time(), `seconds. `):
end: