#John Kim
#Use whatever you like

#read "C:/Users/John Y. Kim/Documents/Maple/hw16.txt":
 


#old stuff from C13.txt
Help13:=proc(): print(` LIF(P,u,N), GuessH1(L,ORD,DEG,n,N)`): end:
 #Lagrange Inversion applied to tree-enumeration
   #of ordered
 #Starting the holonomic ansatz (a.k.a P-recursive)

 #LIF: Inputs an expression P of the variable 
 #u that possesses
 #a Maclaurin expansion in u, and outputs the
  #first N terms of the series expansion of the
  #unique f.p.s u(t) satisfying the functional eq.
  #u(t)=t*P(u(t)). For example,
  #LIF(1+u^2,u,10); should return something with
 #Catalan numbers
 LIF:=proc(P,u,N) local n:

 [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]:
 
  end:

  #GuessH1(L,ORD,DEG,n,N): Inputs a sequence of numbers L
  #and pos. integer ORD and non-neg. integer DEG and letters
  #n and N and outputs an linear recurrence operator with
  #poly. coeffs. (conj.) annihilating the sequence L
  #nops(L)>=(1+ORD)*(1+DEG)+ORD+6  
  GuessH1:=proc(L,ORD,DEG,n,N) local ope,a,i,j,var,eq,n1,var1:
  
  if nops(L)<(1+ORD)*(1+DEG)+ORD+6 then
   print(`We need at least`, (1+ORD)*(1+DEG)+ORD+6, `terms in L`):
   RETURN(FAIL): 
  fi:

  ope:=add(add(a[i,j]*n^i*N^j,i=0..DEG),j=0..ORD):

   
  var:={seq(seq(a[i,j],i=0..DEG),j=0..ORD)}:

   eq:={seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD),
    n1=1..(1+DEG)*(1+ORD)+6)}:

  var1:=solve(eq,var):

  ope:=subs(var1,ope):

   if ope=0 then
     RETURN(FAIL):
    fi:

  ope:=subs({seq(v=1, v in var)}, ope):

  ope:=add(factor(coeff(ope,N,j))*N^j, j=0..degree(ope,N)):


eq:=expand({seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD),
    n1=1..nops(L)-ORD)}):

  if eq<>{0} then
   print(ope, `didn't work out`):
    RETURN(FAIL):
   else
     RETURN(ope):
   fi:


   end:


#end old stuff


 #start new stuff
 #GuessH(L,n,N): Inputs a sequence of numbers L
 # and symbols n and N
  #and finds a minimal (in the sense of ord+deg)
  # linear recurrence operator with
  #poly. coeffs. (conj.) annihilating the sequence L with
  #nops(L)>=(1+ORD)*(1+DEG)+ORD+6  
  GuessH:=proc(L,n,N) local K,ope,DEG,ORD:

  ORD:=1:
   DEG:=0:

  for K from 1 while nops(L)>=(1+ORD)*(1+DEG)+ORD+6 do
    for  ORD from 1 to K do
         DEG:=K-ORD:
         ope:=GuessH1(L,ORD,DEG,n,N):
         if ope<>FAIL then
            RETURN(ope):
          fi:
   
 od:
 od:

  FAIL:

   end:
   
#SpellOut(ope,n,N,f)
#writes the fact ope(n,N)f(n)=0 in common languae 
SpellOut:=proc(ope,n,N,f) local i:

add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0:

end:



#HtoSeq(ope,n,N,Ini,M): Inputs a linear recurrence operator
#ope phrased in terms of the discrete variable n and
#its corresponding shift operator N, and ORD=degree(ope,N)
#values for the sequence outputs the first M terms of
#the sequence. For example,
#HtoSeq(n+1-N,n,N,[1],4); should yield
#[1,2,6,24]
#N^(-ORD)*ope(N,n)
HtoSeq:=proc(ope,n,N,Ini,M)  local L,n1,ORD,BenMiles,kirsten,i:

ORD:=degree(ope,N):

if nops(Ini)<>ORD then
  RETURN(FAIL):
fi:

BenMiles:=expand(subs(n=n-ORD,ope)/N^ORD):

L:=Ini:


for n1 from nops(Ini)+1 to M do

kirsten:=subs(n=n1,coeff(BenMiles,N,0)):

if kirsten=0 then
print(`blows up at`, n1):
print(`so far we have`, L):

  RETURN(FAIL):
fi:

L:=[op(L),
-add(subs(n=n1,coeff(BenMiles,N,-i))*L[n1-i],i=1..ORD)/kirsten]:

od:

L:

end:



#EmpiricalZeilberger(F,n,k,N,M):
#inputs, in addition to the positive integer M, an expression F in the letters n and k, that is a summand of a binomial coefficient sum of the form
#binomial(n,k)^r*(ProductsOrQuotientOfTermsOfTheForm (a*n+b*k+c)!)
#(a,b integers, c number)
#(that be expressible as products/quotients of binomial coefficients if the form binomial(a*n+b*k,a1*n+b1(k))
#and outputs a conjenctured recurrence satisfied by the sequence
#a(n1):=Sum(F(n1,k),k=0..n1);
#For example,
#EmpiricalZeilberger(binomial(n,k)*2^k,n,k,N,30);
#should return
#N-3
#and EmpiricalZeilberger(binomial(n,k)^2,n,k,N,30);
#(n+1)*N-(4*n+2) 

EmpiricalZeilberger:=proc(F,n,k,N,M) local L,n1,k1:

L:=[seq(sum(subs({n=n1,k=k1},F),k1=0..n1),n1=1..M)]:
GuessH(L,n,N):

end:



Help16:=proc(): print(`LIFtoS(F,P,n), NRope(F,P,n,N,K)`):
print(`Rope(F,P,n,N) `):
end:


#LIFtoS(F,P,n): 
#Inputs an expression F in the variable P.
#  and a pos. integern
#outputs the list of the first n coefficients in the formal power
#series that is the UNIQUE solution of the FUNCTIONAL equatation
#P(x)=x*F(P(x))
#For example, to get the the enumerating  sequence for the
#number of complete binary trees with i vertices, up to 10
#do LIFtoS(1+P^2,P,x,10);
#P(x)=x*(1+P(x)^2) 
LIFtoS:=proc(F,P,n) local x,Y,i,p:


Y:=coeff(F,P,0)*x:


for i from 1 to n do

 Y:=x*subs(P=Y,F):

 Y:=add(coeff(Y,x,p)*x^p,p=0..n):

od:

[seq(coeff(Y,x,i),i=1..n)]:

end:



############stuff from AZ.txt##############


#AZ.txt: The Almkvist-Zeilberger algorithm
#taken from
#http://www.math.rutgers.edu/~zeilberg/EKHAD
HelpAZ:=proc(): print(` AZd(A,y,n,N) `): end:

pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:


goremd:=proc(N,ORDER)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:

duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1,apc,apd:
KAMA:=10:
gorem:=goremd(N,ORDER):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
 
 
shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
 
RETURN(shad[1],S):
end:

#AZd(A,y,n,N): inputs a hyperexponential expression A of the
#contintuous variable y and the discrete variable n and
#outputs a linear recurrence operator annihilating the
#integral of A(y,n) along an interval where
#A vanishes at the endpoints (including infinite intervals)
#or a closed circle around the orign. For example try
#AZd((1+y+y^2)^n/y^(n+1),y,n,N);
AZd:= proc(A,y,n,N)
local ORDER,gu,KAMA:
KAMA:=200:
 
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
FAIL:
end:


#####End of stuff from AZ.txt




#NRope(F,P,n,N,K): Inputs an expression F in variable P and
#symbols n and N and  a pos. integer K for
#guessing, outputs a non-rogorously guessed
#linear recurrence operator in the shift operator N
#
#annihilating the sequence of coeffs. of the unique f.p.s.
#solution of
#P(x)=x*F(P(x)).For example,
#NRope(1+P+P^2,P,n,N,30);
NRope:=proc(F,P,n,N,K) local L:

L:=LIFtoS(F,P,K): GuessH(L,n,N):



end:

#Rope(F,P,n,N): Inputs an expression F in variable P and
#symbols n and N and  a pos. integer K for
#guessing, outputs a rigorously derived (using the 
#not-as-famous-as-it-should-be Almkvist-Zeilberger algorithm
#AND Lagrange-Inversion, a
#linear recurrence operator in the shift operator N
#
#annihilating the sequence of coeffs. of the unique f.p.s.
#solution of
#P(x)=x*F(P(x)).For example,
#Rope(1+P+P^2,P,n,N,30);
Rope:=proc(F,P,n,N):

#LIF: the n-th term is the coeff. of P^(n-1) in
#F(P)^n divided by n

AZd((F/P)^n/n,P,n,N)[1]:


end:











#######################
######START HW 16######
#######################

#LargestAB(L): inputs a list of complex numbers, and finds the one with the largest absolute value.  
#(Use the evalf and evalc Maple commands) For example,
#LargestAB([5,-3,-4,6,1]);
#should give 6 and
#LargestAB([5,-3,-4,6,-7,1]);
#should give -7. 

LargestAB:=proc(L) local x,norm,max,maxNorm,numMax:

if nops(L)=0 then
 RETURN(FAIL):
fi:

norm:=0:
max:=0:
maxNorm:=0:
for x in L do:
 norm:=evalf(abs(evalc(x))):
 if norm > maxNorm then
  maxNorm:=norm:
  max:=x:
 fi:
od:

numMax:=0:
for x in L do:
 norm:=evalf(abs(evalc(x))):
 if norm = maxNorm then
  numMax:=numMax+1:
 fi:
od:

if numMax>1 then
 RETURN(FAIL):
fi:

max:

end:


#FindMu(ope,n,N): inputs a linear recurrence operator with polynomial coefficients, ope(n,N), 
#and outputs the root of the largest absolute value of lcoeff(ope,n) (a certain polynomial in N), 
#if it is positive, or returns FAIL if it is negative, or if there are ties (i.e. both a and -a are roots), 
#or if the root(s) with maximal absolute value are complex. For example,
#FindMu((n+2)*N^2-(3*n+3)*N+(2*n+5),n,N);
#should return: 2. 

FindMu:=proc(ope,n,N) local mu, ope1,sol:

ope1:=lcoeff(ope,n):
sol:=[solve(ope1,N)]:
mu:=LargestAB(sol):

if mu=FAIL or evalf(mu)<0 then
 RETURN(FAIL):
fi:

mu:

end:



#FindTheta(ope,n,N,Ini): inputs a linear recurrence operator ope, phrased in terms of n and the shift N, 
#and Ini, a list of length degree(ope,N) indicating the initial values of the sequence, and outputs a 
#guessed RATIONAL number such that the sequence defined by
#a(m):=HtoSeq(ope,n,N,Ini,m)[m] is given by:
#a(m) IS ASYMPOTIC TO   C*mu^m * m^theta
#log(a(m+1)/(mu*a(m))) = theta*log((m+1)/m)
#for some constant C, and the above mu that you found, using FindMu(ope,n,N). 

FindTheta:=proc(ope,n,N,Ini) local mu,a,thetaGuess,theta:

mu:=FindMu(ope,n,N):
if mu=FAIL then
 RETURN(mu):
fi:

a:=HtoSeq(ope,n,N,Ini,10000):
thetaGuess:=evalf(log(a[10000]/(mu*a[9999]))/log(10000/9999)):
convert(thetaGuess,confrac,'cvgts'):
print(cvgts):
theta:=cvgts[3]:

theta:

end:



#FindC(ope,n,N,Ini): estimates the C in 
#a(m) IS ASYMPTOTIC TO   C*mu^m* m^theta

FindC:=proc(ope,n,N,Ini) local mu,theta,a:

mu:=FindMu(ope,n,N):
theta:=FindTheta(ope,n,N,Ini):
a:=HtoSeq(ope,n,N,Ini,10000):

identify(evalf(a[10000]/mu^10000/10000^theta)):

end: 


#Asy(ope,n,N,Ini): inputs a linear recurrence operator ope
#in n and the shift operator N, and the initial values
#(of length the order of ope, i.e. degree(ope,N) and
#outputs (parially non-rig.) the asymptic expression
#for the sequence in the form 
#C*mu^n*n^theta

Asy:=proc(ope,n,N,Ini) local mu,theta,C:

mu:=FindMu(ope,n,N):
theta:=FindTheta(ope,n,N,Ini):
C:=FindC(ope,n,N,Ini):

C*mu^n*n^theta:

end:



#AsyBS(F,n,k,L): inputs a binomial coefficient expression phrased in terms of n and k, 
#and a pos. integer L (for guessing, made large enough) and conjectures asymptotics for
#a(n):=add(F(n,k),k=0..n).

AsyBS:=proc(F,n,k,L) local ope,ORD,Ini,m,j:

ope:=EmpiricalZeilberger(F,n,k,N,L):
ORD:=degree(ope,N):
Ini:=[seq(evalf(add(subs({n=m,k=j},F),j=0..m)) ,m=1..ORD)]:

Asy(ope,n,N,Ini):

end:


#AsyBS(binomial(n,k)^2,n,k,50);
#.56418253122204200526*4^n/sqrt(n)

#AsyBS(binomial(n,k)^3,n,k,50); 
#.36754034533076489914*8^n/n

#AsyBS(binomial(n,k)^4,n,k,50);
#.25396025799131365077*16^n/n^(3/2)

#AsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50);
#.22003455751284834012*(17+12*sqrt(2))^n/n^(3/2)