#Yusra Naqvi: HW #14
#OK to post

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

#1

###

#GuessH1(L,ord,deg,n,N): inputs a sequence L, non-negative integers ord and deg,
#and symbols n and N, and outputs a linear recurrence operator with polynomial 
#coefficients conjecturally annihilating the sequence L

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(cat(`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+ord)*(1+deg)+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
 return(FAIL):
else
 ope:
fi:

end:

###

#GuessH(L,n,N): inputs a sequence L, and symbols n and N, and 
#outputs a minimal (in the sense of order+degree) linear recurrence operator
#with polynomial coefficients conjecturally annihilating the sequence L
#with nops(L)>=(1+ord)*(1+deg)+ord+6

GuessH:=proc(L,n,N) local K,ope,ord,deg:

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:

###

#EmpiricalZeilberger(F,n,k,N,M): inputs a positive integer M, a symbol N, and
#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 are expressible as products/quotients of binomial coefficients of 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 

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

L:=[]:
for i from 1 to M do
 L:=[op(L),expand(add(subs({n=i,k=j},F),j=0..i))];
od:

GuessH(L,n,N):

end:

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

#2
#F=binomial(n,k)^4;
#-(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 
#(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 
#Too long to list here: has degree 9 and order 7.
#Coefficients involve large numbers.

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

#3
#This gives the same output in all the suggested example cases except for:
#F=binomial(n,k)^3*binomial(n+k,k)^2 
#where the output in this case has degree 16 and order 5.
#This is also much faster!

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

#4

###

#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, and
#outputs the first M terms of the sequence. 
#For example:
#HtoSeq(n+1-N,n,N,[1],4);
#should yield 
#[1,2,6,24]

HtoSeq:=proc(ope,n,N,ini,M) local ord,ope1,L,n1,k,i:

ord:=degree(ope,N):

if nops(ini)<>ord then
 return(FAIL):
fi:

ope1:=expand(subs(n=n-ord,ope)/N^ord):

L:=ini:

for n1 from nops(ini)+1 to M do
 k:=subs(n=n1,coeff(ope1,N,0)):
 if k=0 then
  print(cat(`Blows up at `,n1,`.`)):
  print(cat(`So far, we have `,L)):
  return(FAIL):
 fi:
 
 L:=[op(L),-add(subs(n=n1,coeff(ope1,N,-i))*L[n1-i],i=1..ord)/k]:
od:

L:

end:

###

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

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

ope:=EmpiricalZeilberger(F,n,k,N,M1):

ini:=[]:
for i from 1 to degree(ope,N) do
 ini:=[op(ini),expand(add(subs({n=i,k=j},F),j=0..i))];
od:

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

end:

###

#This gives us:

#FirstMterms(binomial(n,k),n,k,20,1000)[1000]; 
#107150860718626732094842504906000181056140481170553360744375038837035105112493
#612249319837881569585812759467291755314682518714528569231404359845775746985748
#039345677748242309854210746050623711418779541821530464749835819412673987675591
#65543946077062914571196477686542167660429831652624386837205668069376

#FirstMterms(binomial(n,k)^2,n,k,20,1000)[1000]; 
#204815162698948971433516250298082504439642488798139703382038263767174818620208
#375582893299418261020620146476631999802369241548179800452479201804754976926157
#856301289663432064714851152395251651227768588611539546256147907378668464154444
#533617613770073855673814589630071306510455959514479888746206368718514551828551
#173166276253663773084682932255389049743859481431755030783796444370810085163724
#827462791417016619883764840843541430817785947037746565188475514680749694674923
#803033101818723298009668567458560252549910118113525353465888794196665367490451
#1306110096311906270342502293155911108976733963991149120

#FirstMterms(binomial(n,k)^4,n,k,20,1000)[1000]; 
#105810444912123795746231188557055866350850558228633681091593665418370891494186
#275699363465887734946955573051849447639465543296792349845753260795072348672281
#435244339445003035470626981308227036571425935313621568747047661757164581052168
#975616577819084529082449013216134436958891223182786100322813203762145710894630
#450188413914920375083595574265703359026647551016680914812050327692150253932517
#553517682754328737817391626069875857890121286766293006090656371389265055516027
#461061014824692725429540752107600345915254195328546681146533848930777668240326
#735505658783193375415535837739494931868932224703232279412455484853533268921741
#004604852764742330532840806455047020202486696009448291836958936259156477010359
#581246025046926254548470278218976807808650103588050982636548700042842293937137
#848051390625872102175960282257212742133546253368076535931728158394219114048806
#842179920660399068531367371391904541836071070807649588481957206221580439672161
#712514061734033231350845783017471436709709958415079454130204505015942750969106
#879666330319919439079121140789444100878283595892393949012927409847571897262243
#425940707180177270366024748975695858386958047762220607527688791414832542037054
#612552185462380343555329821248

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

#5
#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,ord,deg:

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:

###

#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) with
#order ord and degree deg, annihilating L.

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
 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+ord)*(1+deg)+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
 return(FAIL):
else
 ope:
fi:

end:

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

#6

### 

#LIF(P,u,N): inputs an expression P(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 equation u(t)=t*P(u(t)).
#For example:
#LIF(1+u^2,u,12);
#should return
#[1,0,1,0,2,0,5,0,14,0,42,0]

LIF:=proc(P,u,N) local n:

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

end:

###

#OTN(S,N): inputs a set of positive integers, S, and outputs the list L of
#length N of non-neg. integers 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 or 0. 
#Note that this differs from the statement of problem in that we always take 0
#to be one of the possible numbers of children in order to get finite trees.

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

S1:=S union {0}:

P:=0:
for s in S1 do
 P:=P+u^s:
od:

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 union {0}.

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

L:=OTN(S,N1):
GuessH(L,n,N):

end:

###

#We get:
#TreeH({1,2},40,n,N);
#-3*n+(-3-2*n)*N+(3+n)*N^2

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

#7
#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)), 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 x^n*y^n*z^n 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
#(n+1)^2*N-3*(3*n+1)*(3*n+2) 

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

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

end:

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

#8
#F=1/(1-x-y-z+4*x*y*z); 
#-8*(n+1)^2+(-16-21*n-7*n^2)*N+(n+2)^2*N^2
#Sequence A172 in OEIS
#The recurrence is also given.

#F=1/(1-x-y-z+x*y+x*z+y*z);
#(3*(3*n+4))*(3*n+2)+(n+2)^2*N^2
#If you remove the alternating 0's in the sequence, and take absolute values,
#this is Sequence A6480 in OEIS.

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