Help:=proc()
if (args=NULL) then
print(`Functions included are:`):
print(`RecSum(myI,J,F,n,k,N,T,d,eqT,f,a)`):
print(`ReduceT(T,eqT,N,shift,base,n)`):
elif(args = RecSum) then
print(`myI and J are upperbounds on the offsets of n and k respectively`):
print(`n,k,N,T,f,a are all left just as symbols`):
print(`Tries to compute a recurrence that is satisfied f which is`):
print(`the sum over all k of F*g(k)^d`):
print(`where g is the solution to the recurrence eqT`):
print(`for formats of eqT, see the help for procedure ReduceT`):
print(`For Example: RecSum(2,2,binomial(n,k),n,k,N,T,1,1+N,f,a);`):
print(`will return f[n+2] = -f[n]+3*f[n+1]`):
print(`which is the recurrence that is satisfied by the sequence`):
print(`A001906 a_n = \sum_k binomial(n,k)*fibonacci(k)`):
print(`If this procedure fails it will usually output a useless tautology, try increasing myI or J`):
elif(args = ReduceT) then
print(`Format for eqT is that the coefficient i of N is T(k+i), and the whole thing is equal to T(k+degree(eqT,N)+1)`):
print(`n is the value of k in T(k), it should appear in eqT`):
print(`The reduced value will contain expressions T[k-i] for i at most the degree (in N) of eqT`):
print(`For Example: seq(subs({T[1]=1,T[0]=1},(ReduceT(T,((2*n-1)*N+3*(n-1))/n,N,k,2,n)^1)),k=-1..8);`):
print(`gives the first 10 central trinomial coefficients`):
fi:
end:

#F is an expression in n and k
RecSum:=proc(myI,J,F,n,k,N,T,d,eqT,f,a) local Texpr,i,j,eqs,mysum,cNum,sols,comps,comp,myCoeff,CelineRec,CelineRec2:

mysum:=add(add(simplify(expand(subs({n=n+i,k=k+j},F)/F))*a[i][j]*(ReduceT(T,eqT,N,j-1,k,n))^d,i=0..myI),j=0..J):
#print(mysum);
eqs:={}:
comps:= combinat[composition](d+degree(eqT,N)+1,degree(eqT,N)+1):
for comp in comps do
  myCoeff:=mysum:
  for i from 1 to nops(comp) do
    myCoeff:=coeff(myCoeff,T[k-i],comp[i]-1):
  od:
  cNum:= numer(simplify(myCoeff)):
 # print(`coff`,cNum):
  eqs:=eqs union {seq(coeff(cNum,k,i),i=0..degree(cNum,k))}:
od:
sols:=solve(eqs,{seq(seq(a[i][j],i=0..myI),j=0..J)}):
#print(sols):
CelineRec:=add(add(subs({n=n+i,k=k+j},F)*subs(sols,a[i][j])*(ReduceT(T,eqT,N,0,k+j+1,n))^d,i=0..myI),j=0..J):
#print(simplify(CelineRec)):
CelineRec2:=simplify(add(add(f[n+i]*subs(sols,a[i][j]),i=0..myI),j=0..J)):
#print(CelineRec2):
f[n+myI] - simplify(solve(CelineRec2,f[n+myI]));
end:

#N is shift operator
ReduceT:=proc(T,eqT,N,shift,base,n) local myexpr,bigpart,i:
if(shift<0) then
RETURN(T[shift+base]):
fi:
normal(add(subs(n=base+shift,coeff(eqT,N,i))*ReduceT(T,eqT,N,shift - degree(eqT,N)+i-1,base,n),i=0..degree(eqT,N))):
end:

ReduceBinom:=proc(n,k,i,j) local expr,N1,N2,l:
expr:=N1^i*N2^j*binomial(n,k):
if(i>j) then
  for l from 1 to i-j do
    expr:=expr*N2*(k+degree(expr,N2))/(n+degree(expr,N1)+1-k-degree(expr,N2)):
  od:
fi:
if(j>i) then
  for l from 1 to j-i do
    expr:=expr/N2/(k+degree(expr,N2))*(n+degree(expr,N1)+1-k-degree(expr,N2)):
  od:
fi:
for l from 1 to i do
  expr:=expr/N2/N1*(n+degree(expr,N1))/(k+degree(expr,N2)):
od:
expr:
end:


RecToSeq:=proc(Rec,f,Ini,n,numterms,d) local i,ret,testrec,degree,myexpr:
testrec:=Rec:
#for i from 0 to nops(Ini)-1 do
#print(testrec, coeff(testrec,f[n+i],1)):
#testrec:=testrec - coeff(testrec,f[n+i],1)*f[n+i]:
#od:
#testrec:=testrec:
if(nops(Ini)<d) then
  print(`give longer Ini`):
#  print(`unknown terms are`,testrec)
  RETURN(FAIL):
fi:
ret:=Ini:
i:=1:
while(nops(ret)<numterms) do
myexpr:=subs(n=i,Rec):
myexpr:=subs({seq(f[i+j]=ret[i+j],j=0..d-1)},myexpr):
#print(myexpr):
ret:=[op(ret),solve(myexpr,f[i+d])]:
#print(ret):
i:=i+1:
od:
RETURN(ret):
end: