#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 16
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.
Help := proc(): print(`LargestAB(L), FindMu(ope, n, N), HtoSeq(ope,n,N,Ini,M), FindTheta(ope, n, N, ini),
FindC(ope, n, N, ini, M:=1000), Asy(ope, n, N, ini), EmpiricalZeilberger(F,n,k,N,M), AsyBS(F, n, k, L)`): end:
#############
# Problem 1 #
#############

LargestAB := proc(L) local i, mc:
    mc := 0:
    for i in L do:
        if evalf(abs(i)) >= evalf(abs(mc)) then:
            mc := i:
        fi:
    od:
    return mc:
end:

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

FindMu := proc(ope, n, N) local roots, largest:
    roots := convert({solve(lcoeff(ope, n) = 0)}, list):
    largest := LargestAB(roots):
    if Im(largest) <> 0 then:
        return FAIL:
    fi:
    if LargestAB(convert(convert(roots, set) minus {largest}, list)) = largest then:
        return FAIL:
    fi:
    if evalf(largest) < 0 then:
        return FAIL:
    fi:
    return largest:
end:

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

# We need HtoSeq from class 14:

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:

##

FindTheta := proc(ope, n, N, ini) local L, mu, absr2, m, cvgts:
    mu := FindMu(ope, n, N):
    if mu = FAIL then:
        return FAIL:
    fi:

    # To decide how many terms to compute, let r2 be the
    # second-largest root of lcoeff(ope, n, N). We'll continue until
    # (abs(r2)/abs(mu))^m < 0.001, where 0.001 is just a constant
    # chosen arbitrarily but is hopefully good enough.

    absr2 := max(seq(abs(i), i in {solve(lcoeff(ope, n))} minus {mu})):
    if mu = 0 then:
        # We're in trouble!
        return FAIL:
    fi:
    if evalf(absr2) <= 0 then:
        m := 100 # aribtrary
    else:
        m := ceil(log(0.001)*(log(absr2)-log(abs(mu)))):
    fi:
    L := HtoSeq(ope, n, N, ini, m):
    if L = FAIL then:
        return FAIL:
    fi:
    convert(log(L[m]/L[m-1]/mu)/log((m+1)/m), confrac, cvgts):
    if nops(cvgts) = 1 then:
        return cvgts[1]:
    else:
        return cvgts[2]:
    fi:
end:

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

# M = number of terms to compute in approximating C.
# The call to identify is very slow and the user can provide it
# manually if desired.
FindC := proc(ope, n, N, ini, M:=1000) local L, mu, theta, L2:
    L := HtoSeq(ope, n, N, ini, M):
    mu := FindMu(ope, n, N):
    theta := FindTheta(ope, n, N, ini):
    return evalf(L[M]/mu^M/M^theta):
end:

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

Asy := proc(ope, n, N, ini):
    return [FindMu(ope, n, N), FindTheta(ope, n, N, ini), FindC(ope, n, N, ini)];
end:

###############################################################################
# From HW 14:

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)];
    return GuessH(L,n,N):
end:

# From Class 14:

  #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:

###############################################################################

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

AsyBS := proc(F, n, k, L) local N, ope, ord:
    ope := EmpiricalZeilberger(F,n,k,N,L):
    printf("Found operator %s...\n", convert(ope, string)):
    ord := degree(ope, N):
    return Asy(ope, n, N, [seq(sum(subs({n=i,k=j}, F), j=0..i), i=1..ord)]):
end:

# AsyBS(binomial(n,k)^2,n,k,50) = [4, -1/2, 0.2820595321]
# (that is, 0.28206*4^n/sqrt(n)
# identify is not able to determine a formula for 0.2820595321, but it
# is easy to get (from Stirling's formula) that the constant should be
# 1/(2sqrt(pi)), which is 0.2820947918.

# AsyBS(binomial(n,k)^3,n,k,50) = [8, -1, 0.3674300930]
# AsyBS(binomial(n,k)^4,n,k,50) = [16, -3/2, 0.2538317153]
# AsyBS(binomial(n+k,k)^2*binomial(n,k)^2,n,k,50) : I had to time it
# out after 10 minutes. (I can't figure out what is so slow!)