# Alex Gentul Homework 14 OK to share

info1:=proc() print(`EmpiricalZeilberger(expression F,symbol n, symbol k,symbol N,integer M); EZTest(); FirstMterms(expression F,symbol n,symbol k,integer M1,integer M2); FMTest(); qGuessH(sequence L,symbol n,symbol N,variable q); TreeH(set S,integer N1,symbol n,symbol N); DiagRec(function F,variable x,variable y,variable z,integer M,symbol n,symbol N)  `); end:
info1();

	# (1) #

	# EmpiricalZeilberger(expression F,symbol n, symbol k,symbol N,integer M)

EmpiricalZeilberger:=proc(F,n,k,N,M) local L;
with(combinat);
L:=eval([seq(eval(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..M)]);

GuessH(L,n,N);

end:

	# (2) and (3) #

	# EZTest()

EZTest:=proc()

print(`Emp`,expand(EmpiricalZeilberger(binomial(n,k)^4,n,k,N,50))); 
print(`Maple`,expand(SumTools[Hypergeometric][Zeilberger](binomial(n,k)^4,n,k,N)[1])); 
print(`Emp`,expand(EmpiricalZeilberger(binomial(n,k)^2*binomial(n+k,k)^2,n,k,N,50))); 
print(`Maple`,expand(SumTools[Hypergeometric][Zeilberger](binomial(n,k)^2*binomial(n+k,k)^2,n,k,N)[1])); 
print(`Emp`,expand(EmpiricalZeilberger(binomial(n,k)^3*binomial(n+k,k)^2,n,k,N,160))); 
print(`Maple`,expand(SumTools[Hypergeometric][Zeilberger](binomial(n,k)^3*binomial(n+k,k)^2,n,k,N)[1])); 

end:


	# (4) #

	# FirstMterms(expression F,symbol n,symbol k,integer M1,integer M2)  
FirstMterms:=proc(F,n,k,M1,M2) local ope,L;

ope:=EmpiricalZeilberger(F,n,k,N,M1);
L:=eval([seq(eval(add(subs({n=n1,k=k1},F),k1=0..n1)),n1=1..degree(ope,N))]);
HtoSeq(ope,n,N,L,M2);

end:

	# FMTest()

FMTest:=proc()

print(FirstMterms(binomial(n,k),n,k,20,1000)[1000]);
print(FirstMterms(binomial(n,k)^2,n,k,20,1000)[1000]);
print(FirstMterms(binomial(n,k)^4,n,k,20,1000)[1000]);

end:


	# (5) #

	# qGuessH(sequence L,symbol n,symbol N,variable q)

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

  FAIL:

   end:

qGuessH1:=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]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
   print(ope, `didn't work out`):
    RETURN(FAIL):
   else
     RETURN(ope):
   fi:


   end:


	# (6) #

	# TreeH(set S,integer N1,symbol n,symbol N) 

TreeH(S,N1,n,N) local A,a,L;

A:=add(a^S[i],i=1..nops(S));
L:=LIF(A,a,N);
GuessH(L,n,N);

end:


	# (7) #

	# DiagRec(function F,variable x,variable y,variable z,integer M,symbol n,symbol N) 

DiagRec(F,x,y,z,M,n,N) 








































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