# Matthew Russell
# Experimental Math
# Homework 14
# I give permission to post this

read `public_html/em12/C13.txt`:
read `public_html/em12/C14.txt`:

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

# By using GuessH, and by first generating the first M terms of a binomial coefficient sum,
# write a procedure,
# EmpiricalZeilberger(F,n,k,N,M)
# that 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);

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:

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

# Try EmpiricalZeilbegrer out (with M large enough) and
#   * F1=binomial(n,k)^4
#   * F2=binomial(n,k)^2*binomial(n+k,k)^2
#   * F3=binomial(n,k)^3*binomial(n+k,k)^2 

# Output for F1:
# -4*(4*n+5)*(4*n+3)*(n+1)-2*(2*n+3)*(3*n^2+9*n+7)*N+(n+2)^3*N^2
# Output for F2:
# (n+1)^3-(2*n+3)*(17*n^2+51*n+39)*N+(n+2)^3*N^2
# Output for F3:
# (Omitted due to size; it returned a degree-9, order-7 recurrence)

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

# Compare the above outputs with the Maple procedure
# SumTools[Hypergeometric][Zeilberger](F,n,k,N)[1] ;

# Output for F1:
# (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, as shown by simplify)
#
# Output for F2:
# (n^3+6*n^2+12*n+8)*N^2+(-34*n^3-153*n^2-231*n-117)*N+3*n+1+n^3+3*n^2
# (same as EmpiricalZeilberger, as shown by simplify)
#
# Output for F3:
# (3156425287*n^16+166415220936*n^15+4074844116325*n^14+61498559727260*n^13+640172471751394*n^12
# +4872692486746748*n^11+28047431954557238*n^10+124511211463021580*n^9+430737539643208151*n^8
# +1164813743532739780*n^7+2453578828663815173*n^6+3982505680032003304*n^5+4881992669157267744*n^4
# +4368384221599323880*n^3+2690154674576001600*n^2+1018445595942672000*n+178527210071040000)*N^5
# +(-1211619464770238154*n^10-216815927488303005*n^11-200267172803497*n^14-8485425629186*n^15-167290540211*n^16
# -2918670844608321*n^13-97172055884733806784*n^2-160533108500200656104*n^3-182772812605736794696*n^4
# -29392033464083602*n^12-152109839143693573272*n^5-95747398763994516869*n^6-46512766707912198813*n^7
# -5230904376648202499*n^9-36206611716191253120*n-6254313816674893824-17627951730325162043*n^8)*N^4
# +(-1751689608347831*n^14-568997303095543701024*n^2-966728401405852173360*n^3-1133131123746236452196*n^4
# -971853830099819479672*n^5-1584525494074*n^16-9132026529819485757*n^10-238220965947661415*n^12
# -24563723416900436*n^13-206374286524503995136*n-1693830244246437584*n^11-631076437791894778768*n^6
# -316567196755095106096*n^7-77202338933576*n^15-34739924439108507840-124009230668323369495*n^8-38071512294985967932*n^9)*N^3
# +(-4679715452127809001*n^12-7583607653514254290848*n^2-2661171465870805396032*n-660535209206907173697*n^9
# -1720462823050524*n^15-2066356262067961347465*n^8-434027162953117248192-37429006976899209*n^14
# -31912775772437938311*n^11-165067834351062788639*n^10-16198267945257350167676*n^4-14414016117751691640372*n^5
# -9721310346166804361868*n^6-5069692873453886555985*n^7-13335403697005011927408*n^3-36822857398142*n^16-503234007303291111*n^13)*N^2
# +(-3510763960094640000-133850617599135363264*n^3-9280145415246130*n^13-2595828050407894058*n^10-70898837136517466112*n^2
# -779637045889*n^16-724425035650330*n^14-66574810344510743354*n^7-174328360572932286696*n^4-166070581377844837892*n^5
# -119692925523077495756*n^6-23149740187448631360*n-28880410286523533735*n^8-9803942551520133938*n^9-82031184499920736*n^12
# -34867463204080*n^15-530477077097799022*n^11)*N
# -42127403291374464*n-269701936132*n^15-5265994423350528-518344649519172976*n^4-565668751525906904*n^5
# -461224250488658984*n^6-286876109392130180*n^7-137694180628217774*n^8-51220159882437968*n^9-14731237157676456*n^10
# -3244281056259256*n^11-536817719656228*n^12-64568218913944*n^13-5328162498400*n^14-154187771896707840*n^2-342788693239830912*n^3-6312850574*n^16
# This is actually different than the one found by EmpiricalZeilberger

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

# By combining EmpiricalZeilberger with
# HtoSeq(ope,n,N,Ini,M):
# write a program,
# FirstMterms(F,n,k,M1,M2);
# That 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.
# What are
#   * FirstMterms(binomial(n,k),n,k,20,1000)[1000];
#   * FirstMterms(binomial(n,k)^2,n,k,20,1000)[1000];
#   * FirstMterms(binomial(n,k)^4,n,k,20,1000)[1000]; 

FirstMterms:=proc(F,n,k,M1,M2) local Ini,EZ,d:

EZ:=EmpiricalZeilberger(F,n,k,N,M1):
d:=degree(EZ,N):
Ini:=[seq(expand(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..d)]:
return HtoSeq(EZ,n,N,Ini,M2):

end:

# We find:
# FirstMterms(binomial(n,k),n,k,20,1000)[1000];
# gives
# 107150860718626732094842504906000181056140481170553360744375038837035105112493612249319837881=
# 569585812759467291755314682518714528569231404359845775746985748039345677748242309854210746050=
# 623711418779541821530464749835819412673987675591655439460770629145711964776865421676604298316=
# 52624386837205668069376,
# as it should (this is 2^1000).
#
# FirstMterms(binomial(n,k)^2,n,k,20,1000)[1000];
# gives
# 204815162698948971433516250298082504439642488798139703382038263767174818620208375582893299418=
# 261020620146476631999802369241548179800452479201804754976926157856301289663432064714851152395=
# 251651227768588611539546256147907378668464154444533617613770073855673814589630071306510455959=
# 514479888746206368718514551828551173166276253663773084682932255389049743859481431755030783796=
# 444370810085163724827462791417016619883764840843541430817785947037746565188475514680749694674=
# 923803033101818723298009668567458560252549910118113525353465888794196665367490451130611009631=
# 1906270342502293155911108976733963991149120,
# as it should (this is the 1000th central binomial coefficient).
#
# FirstMterms(binomial(n,k)^4,n,k,20,1000)[1000]; 
# gives
# 105810444912123795746231188557055866350850558228633681091593665418370891494186275699363465887=
# 734946955573051849447639465543296792349845753260795072348672281435244339445003035470626981308=
# 227036571425935313621568747047661757164581052168975616577819084529082449013216134436958891223=
# 182786100322813203762145710894630450188413914920375083595574265703359026647551016680914812050=
# 327692150253932517553517682754328737817391626069875857890121286766293006090656371389265055516=
# 027461061014824692725429540752107600345915254195328546681146533848930777668240326735505658783=
# 193375415535837739494931868932224703232279412455484853533268921741004604852764742330532840806=
# 455047020202486696009448291836958936259156477010359581246025046926254548470278218976807808650=
# 103588050982636548700042842293937137848051390625872102175960282257212742133546253368076535931=
# 728158394219114048806842179920660399068531367371391904541836071070807649588481957206221580439=
# 672161712514061734033231350845783017471436709709958415079454130204505015942750969106879666330=
# 319919439079121140789444100878283595892393949012927409847571897262243425940707180177270366024=
# 748975695858386958047762220607527688791414832542037054612552185462380343555329821248

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

# A sequence is q-holonomic if it satisfies a linear recurrence equation with coefficients that
# are polynomials in q^n and q, for example the sequence
# qfac(n)=(1-q)*(1-q^2)*....*(1-q^n)
# where q is a symbol. Write a procedure
# 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.
#
# Hint, think of ope(q^n,q,n) as ope(q^n,N) over the field of rational functions in q.
# When you solve a linear system, now your system has coefficients that depend on q,
# so it takes longer for Maple to solve the system. 

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

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:

return `FAIL`:

end:

qGuessH1:=proc(L,ord,deg,n,N,q) local ope,a,i,j,var,var1,eq:

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+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
	return `FAIL`:
else
	return ope:
fi:

end:

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

# By using LIF(P,u,N1), write a procedure
# TreeH(S,N1,n,N)
# that 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 s,P,L:

P:=1+add(u^s,s in S):
L:=LIF(P,u,N1):
return GuessH(L,n,N):

end:

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

# Write a procedure
# DiagRec(F,x,y,z,M,n,N)
# that 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 x^ny^nz^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 (I hope, I did it in my head)
# (n+1)^2*N-3*(3*n+1)*(3*n+2) .
# Hint: Either use mtaylor or iterated taylor.

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

m:=mtaylor(F,[x,y,z],3*M+1):
L:=[seq(coeff(coeff(coeff(m,x,i),y,i),z,i),i=1..M)]:
print(L):
return GuessH(L,n,N):

end:

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

# Apply DiagRec(F,x,y,z,M,n,N) to
#   * F1=1/(1-x-y-z+4*x*y*z)
#   * F2=1/(1-x-y-z+x*y+x*z+y*z)
#     Are these sequences in OEIS? Are their recurrences mentioned? 

# For F1, we find
# -8*(1+n)^2+(-16-21*n-7*n^2)*N+(n+2)^2*N^2
# The sequence is 
# 1, 2, 10, 56, 346, 2252, 15184, 104960, 739162, 5280932, 38165260, 278415920,...
# which is A000172, the Franel numbers. a(n) = Sum C(n,k)^3, k=0..n
# This recurrence is mentioned.
#
# For F2, we find
# 3*(3*n+4)*(3*n+2)+(n+2)^2*N^2
# The sequence is
# 1, 0, -6, 0, 90, 0, -1680, 0, 34650, 0, -756756, 0, 17153136, 0, -399072960,...
# which is not in OEIS.
# However, the sequence
# 1, 6, 90, 1680, 34650, 756756, 17153136, 399072960, ...
# is A006480, De Bruijn's S(3,n): (3n)!/(n!)^3. 
# Its corresponding recurrence is given by 
# -3*(3*n-1)*(3*n-2)+n^2*N
# This actually does not appear to be in OEIS.
#####################################################################################