### Homework 16 - Patrick Devlin - March 23 2012 ###
## I have no preference about keeping my work private ##

Help:=proc():
 print(`LargestAB(L), FindMu(ope, n, N), FindTheta(ope, n, N, Ini), FindC(ope, n, N, Ini)`):
 print(`Asy(ope, n, N, Ini), AsyBS(F, n, k, L), BetterAsy(ope, n, N, Ini, r)`):
 print(`-- Type  HelpAll()  to see all functions --`):
end:

HelpAll:=proc():
 print(`HtoSeq(ope, n, N, Ini, M), EmpiricalZeilberger(F, n, k, N, L), GuessH(L, n, N)`):
end:


#############
# Problem 1 #
#############

LargestAB:=proc(L) local biggestModulus, biggestIndex, i:
  biggestIndex:=1:  biggestModulus:=evalf(abs(L[1])):
  for i from 1 to nops(L) do
    if (evalf(abs(L[i])) > biggestModulus) then
      biggestModulus:=evalf(abs(L[i])):
      biggestIndex:=i:
    fi:
  od:
  return L[biggestIndex]:
end proc:


#############
# Problem 2 #
#############

FindMu:=proc(ope, n, N) local L, largest:
  L:=[fsolve(lcoeff(ope,n)=0, N)]:
  try:
    largest:=LargestAB(L):
    if(largest>0) then return largest: fi:
  catch:
    return FAIL:
  end try:
  return FAIL:
end proc:


#############
# Problem 3 #
#############

FindTheta:=proc(ope, n, N, Ini, big:=10000):
  return FindThetaGivenMu(ope, n, N, Ini, FindMu(ope, n, N), big):
end proc:

FindThetaGivenMu:=proc(ope, n, N, Ini, mu, big:=10000) local L, th:
  if(mu = FAIL) then return FAIL: fi:
  L:=HtoSeq(ope, n, N, Ini, big):
  try:
   th:=convert( log(L[big]/(L[big-1] *mu))/log(big/(big-1)), confrac, cvgsts):
   return th[1] + 1/th[2]:
  catch:
   return FAIL:
  end try:
  
end proc:


#############
# Problem 4 #
#############

FindC:=proc(ope, n, N, Ini, big:=10000) local mu:
  mu:=FindMu(ope, n, N):
  FindCGivenMuAndTheta(ope, n, N, Ini, mu, FindThetaGivenMu(ope, n, N, Ini,
mu, big)):
end proc:


FindCGivenMuAndTheta:=proc(ope, n, N, Ini, mu, theta, big:=10000) local L, C:
  if( (mu = FAIL) or (theta = FAIL)) then return FAIL: fi:
  L:=HtoSeq(ope, n, N, Ini, big):
  try:
   return identify(L[big]/( (mu^big) * (big ^ theta))):
  catch:
   return FAIL:
  end try:
end proc:

#############
# Problem 5 #
#############

Asy:=proc(ope, n, N, Ini, big:=100000) local mu, theta, C, m:
  mu:=FindMu(ope, n, N):
  theta:=FindThetaGivenMu(ope, n, N, mu, big):
  C:=FindCGivenMuAndTheta(ope, n, N, Ini, mu, theta, big):
  if((mu = FAIL) or (theta = FAIL) or (C = FAIL)) then return FAIL: fi:
  return (C* (mu^m) * (m ^theta)):
end proc:


#############
# Problem 6 #
#############

AsyBS:=proc(F, n, k, L) local Ini, n1, k1, ope:
  Ini:=[seq(expand(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..L)]:
  ope:=EmpiricalZeilberger(F, n, k, N, L):
  if(ope = FAIL) then return FAIL: fi:
  return Ays(ope, n, N, Ini):
end proc:


####################
# From other stuff #
####################

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

L:=[seq(expand(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..M)]:
return GuessH(L,n,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:





#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:
 return L:
end: