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

# Nathaniel Shar
# HW 14
# Experimental Mathematics

###############################################################################
# Stuff from Class 14
###############################################################################

Help:=proc(): print( ` GuessH(L,n,N) , SpellOut(ope,n,N,f)`):

print(`HtoSeq(ope,n,N,Ini,M)`):
end:

#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 ORD, DEG, ope:

for ORD from 0 to (nops(L)-7)/2 do
for DEG from 0 to (nops(L)-ORD-6)/(1+ORD)-1 do

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:

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

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

# We trust the user to provide an F of the appropriate form.

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:

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

# F=binomial(n,k)^4, M=50:
#     -4*(4*n+5)*(4*n+3)*(1+n)-2*(2*n+3)*(3*n^2+9*n+7)*N+(n+2)^3*N^2
# F=binomial(n,k)^2*binomial(n+k,k)^2, M=50:
#     (1+n)^3-(2*n+3)*(17*n^2+51*n+39)*N+(n+2)^3*N^2
# F=binomial(n,k)^3*binomial(n+k,k)^2, M=100:
#     something too big to include here

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

# Everything seems to check out.

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

FirstMterms := proc(F,n,k,M1,M2) local deg, oper, k1, n1, L:
    oper := EmpiricalZeilberger(F,n,k,N,M1):
    deg := degree(oper, N):
    L := [seq(sum(subs({n=n1,k=k1},F),k1=0..n1), n1=1..deg)]:
    return HtoSeq(oper, n, N, L, M2):
end:

#####################
# Problem 4, part 2 #
#####################

# 1. It is 2^1000
# 2. It is binomial(2000,1000)
# 3. It is the following large number:
# 105810444912123795746231188557055866350850558228633681091593665418370891494186275699363465887734946955573051849447639465543296792349845753260795072348672\
#    28143524433944500303547062698130822703657142593531362156874704766175716458105216897561657781908452908244901321613443695889122318278610032281320376214\
#    57108946304501884139149203750835955742657033590266475510166809148120503276921502539325175535176827543287378173916260698758578901212867662930060906563\
#    71389265055516027461061014824692725429540752107600345915254195328546681146533848930777668240326735505658783193375415535837739494931868932224703232279\
#    41245548485353326892174100460485276474233053284080645504702020248669600944829183695893625915647701035958124602504692625454847027821897680780865010358\
#    80509826365487000428422939371378480513906258721021759602822572127421335462533680765359317281583942191140488068421799206603990685313673713919045418360\
#    60355721577725808747721971089839654056098465944678697156913341899313960925460348733483095665973245923739996021793575411414458867062515938685156211869\
#    40230347966132461912478194198877296683489110154798652502713079538244280770829709082779250381614354461311759493714599986939418272053283799350634966175\
#    1872

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

qGuessH1:=proc(L,ORD,DEG,n,N,q) 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]*(q^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(expand(ope),N,j))*L[n1+j],j=0..ORD)=0,
             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:

qGuessH := proc(L, n, N,q) local oper, imax, i, ord, deg:
    imax := floor((nops(L)-7)/2):
    for i from 1 to imax do:
        for ord from 1 to i do:
            for deg from 0 to i-ord do:
                if nops(L) >= (1+ord)*(1+deg)+ord+6 then:
                    oper := qGuessH1(L,ord,deg,n,N,q):
                fi:
                if oper <> FAIL then:
                    return oper:
                fi:
            od:
        od:
    od:
    return FAIL:
end:

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

TreeH := proc(S,N1,n,N) local i, u:
    return GuessH(LIF(1+add(u^i, i in S), u, N1), n, N):
end:

#############
# Problem 7 #
#############

# Note: I think I used a slightly different convention for the
# starting index.

# Also note: This is slow. It could be made faster by using an
# iterated approach. If I ever have to use this function for something
# else, I will do that, but for now, this seems good enough.

DiagSeq := proc(F,x,y,z,M) local T:
    T := mtaylor(F, [x,y,z], 3*M):
    return [seq(coeff(coeff(coeff(T,z,i),y,i),x,i), i=0..M-1)]:
end:

DiagReq := proc(F,x,y,z,M,n,N) local L:
    L := DiagSeq(F,x,y,z,M):
    return GuessH(L,n,N):
end:

#############
# Problem 8 #
#############

# 1. Yes, it is in OEIS as A000172, and the recurrence is mentioned
# 2. The sequence contains 0s and negative terms. After removing the
# 0s and taking absolute value, it is in the OEIS (A0006480). The
# recurrence for that sequence does not seem to be mentioned, probably
# because it is trivial from the explicit formula that is provided.