#John Kim
#Use whatever you like

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

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



#F=binomial(n,k)^4: yields -(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: yields (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: wouldn't run in reasonable time.

#SumTools[Hypergeometric][Zeilberger](binomial(n,k)^4,n,k,N)[1];
#(n^3+6*n^2+12*n+8)*N^2+(-12*n^3-54*n^2-82*n-42)*N-60-188*n-192*n^2-64*n^3 (same as EmpiricalZeilberger)

#SumTools[Hypergeometric][Zeilberger](binomial(n,k)^2*binomial(n+k,k)^2,n,k,N)[1];
#(n^3+6*n^2+12*n+8)*N^2+(-34*n^3-153*n^2-231*n-117)*N+1+3*n+3*n^2+n^3 (same as EmpiricalZeilberger)

#SumTools[Hypergeometric][Zeilberger](binomial(n,k)^3*binomial(n+k,k)^2,n,k,N)[1];
#takes too long to run



#FirstMterms(F,n,k,M1,M2):
#First finds the first M1 terms of the sequence
#a(n1):=Sum(F(n1,k),k=0..n1);
#guesses a recurrence, and uses it to output the first M2 terms of the same sequence. 

FirstMterms:=proc(F,n,k,M1,M2) local ope,ORD,ini:

ope:=EmpiricalZeilberger(F,n,k,N,M1):
ORD:=degree(ope,N):
ini:=[seq(sum(subs({n=n1,k=k1},F),k1=0..n1),n1=1..ORD)]:

HtoSeq(ope,n,N,ini,M2):

end:


#FirstMterms(binomial(n,k),n,k,20,1000)[1000]; outputs:
#10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

#FirstMterms(binomial(n,k)^2,n,k,20,1000)[1000]; 
#2048151626989489714335162502980825044396424887981397033820382637671748186202083755828932994182610206201464766319998023692415481798004524792018047549769261578563012896634320647148511523952516512277685886115395462561479073786684641544445336176137700738556738145896300713065104559595144798887462063687185145518285511731662762536637730846829322553890497438594814317550307837964443708100851637248274627914170166198837648408435414308177859470377465651884755146807496946749238030331018187232980096685674585602525499101181135253534658887941966653674904511306110096311906270342502293155911108976733963991149120

#FirstMterms(binomial(n,k)^4,n,k,20,1000)[1000];
#105810444912123795746231188557055866350850558228633681091593665418370891494186275699363465887734946955573051849447639465543296792349845753260795072348672281435244339445003035470626981308227036571425935313621568747047661757164581052168975616577819084529082449013216134436958891223182786100322813203762145710894630450188413914920375083595574265703359026647551016680914812050327692150253932517553517682754328737817391626069875857890121286766293006090656371389265055516027461061014824692725429540752107600345915254195328546681146533848930777668240326735505658783193375415535837739494931868932224703232279412455484853533268921741004604852764742330532840806455047020202486696009448291836958936259156477010359581246025046926254548470278218976807808650103588050982636548700042842293937137848051390625872102175960282257212742133546253368076535931728158394219114048806842179920660399068531367371391904541836071070807649588481957206221580439672161712514061734033231350845783017471436709709958415079454130204505015942750969106879666330319919439079121140789444100878283595892393949012927409847571897262243425940707180177270366024748975695858386958047762220607527688791414832542037054612552185462380343555329821248



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

#print(ope,eq,var):

  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(L,n,N,q): Inputs a sequence L, of polynomials or rational functions in q and tries to conjecture a recurrence operator
#ope(q^n,q,N) annihilating L. 

qGuessH:=proc(L,n,N,q) 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:=qGuessH1(L,ORD,DEG,n,N,q):
         if ope<>FAIL then
            RETURN(ope):
          fi:
   
 od:
 od:

  FAIL:

end:


#OTN(S,N): Inputs a set of non-negative integers, S, and outputs the list of non-neg. integers, 
#let's call it L, such that L[n] is the number of ordered trees with n vertices where every vertex 
#must have a number of children that is a member of S. For example,
#OTN({0,2},7);
#should give
#[1,0,1,0,2,0,5]
#and OTN({0,1,2},7);
#should give the first seven terms of the Motzkin numbers sequence. 

OTN:=proc(S,N) local P,i,u:

P:=add(u^S[i],i=1..nops(S)):
LIF(P,u,N):

end:


#TreeH(S,N1,n,N): inputs a set of non-negative integers, S, and conjectures a linear recurrence operator, 
#ope(n,N) annihilating the enumerating sequence for the number of ordered trees with n vertices where 
#every vertex must have a number of children that is a member of S. N1 is a positive integer indicating 
#how many terms of the sequence to use for guessing. If it fails, the procedure should return FAIL. 
#For example
#TreeH({1,2},40,n,N);
#should return the recurrence operator annihilating the Motzkin number sequence. 

TreeH:=proc(S,N1,n,N) local L:

L:=OTN(S union {0},N1):
GuessH(L,n,N):

end:



#DiagRec(F,x,y,z,M,n,N): inputs a rational function F of the variables x,y,z 
#(whose denominator does not vanish at (x,y,z)=(0,0,0) so it has a Taylor series around the origin, 
#and a positive integer M (indicating how many terms to use for guessing), and symbols n and N and 
#outputs a linear recurrence operator with polynomial coefficients, ope(n,N), annihilating the sequence
#a(n):=coeff of xnynzn in the Taylor expansion of F around the origin.
#For example,
#DiagRec(1/(1-x-y-z),x,y,z,40,n,N);
#should return (I hope, I did it in my head)
#(n+1)^2*N-3*(3*n+1)*(3*n+2).

DiagRec:=proc(F,x,y,z,M,n,N) local L,i:

L:=[seq(coeftayl(F,[x,y,z]=[0,0,0],[i,i,i]),i=1..M)]:
GuessH(L,n,N):

end:



#F=1/(1-x-y-z+4*x*y*z) yields:
#-8*(1+n)^2+(-16-21*n-7*n^2)*N+(n+2)^2*N^2
#Franel number a(n) = Sum C(n,k)^3, k=0..n. 

#F=1/(1-x-y-z+x*y+x*z+y*z) yields:
#(3*(3*n+4))*(3*n+2)+(n+2)^2*N^2
#Number of paths of length 3n in an n X n X n grid from (0,0,0) to (n,n,n).