Help:=proc()
if (args=NULL) then
print(`Functions included are:`):
print(`RecSum(maxI,maxJ,timeout,D,eqT,F,n,k,N,T,f,a)`):
print(`ReduceT(T,eqT,N,shift,base,n)`):
print(`RecToSeq(Rec,f,Ini,n,numterms,d)`):
elif(args = RecSum) then
print(`maxI and maxJ are upperbounds on the offsets of n and k respectively`):
print(`increasing maxI will allow looking for higher order recurrences`):
print(`increasing maxJ will take more time looking for the lower order recurrences`):
print(`timeout is the time in seconds that if a case longer than, the procedure will fail`):
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,10,1,1+N,binomial(n,k),n,k,N,T,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 least 0 and 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`):
elif (args= RecToSeq) then
print(`Rec is the recurrence output from RecSum`):
print(`f and n are left as symbols`):
print(`Ini is the initial segment to start the recurrence off`):
print(`d is the degree of the recurrence`):
print(`for example RecToSeq(f[n+3] - f[n+1] - f[n]-f[n+2],f,[0,0,1],n,10,3)`):
print(`will give you the first 10 tribonacci numbers`):
fi:
end:

RecSum:=proc(maxI,maxJ,timeout,d,eqT,F,n,k,N,T,f,a) local myI,myJ,ret,t:
for myI from 7 to maxI do
for myJ from 7 to maxJ do
print(`i,j,time is`,myI,myJ,time()):
t:=time();
ret:=RecSum1(myI,myJ,F,n,k,N,T,d,eqT,f,a):
if(ret<>FAIL) then
RETURN(ret):
elif (time()-t > timeout) then
RETURN(FAIL):
fi:
od:
od:
FAIL:
end:

RecTribonacci:=((2*(n)-1)*N+3*((n)-1))/(n):
RecMotzkin:=(3*n-3)/(n+2) + (2*n+1)/(n+2)*N:

#F is an expression in n and k
RecSum1:=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,k,n))^d,i=0..myI),j=0..J):
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+1],comp[i]-1):
  od:
  cNum:= numer(simplify(myCoeff)):
  eqs:=eqs union {seq(coeff(cNum,k,i),i=0..degree(cNum,k))}:
od:
#print(eqs):
if(QuickCheckSystem(eqs,{seq(seq(a[i][j],i=0..myI),j=0..J)},3,{n}) = false) then
#print(`caughtquick`):
RETURN(FAIL):
fi:
sols:=solve(eqs,{seq(seq(a[i][j],i=0..myI),j=0..J)}):
CelineRec2:=simplify(add(add(N^(n+i)*subs(sols,a[i][j]),i=0..myI),j=0..J)):
myCoeff:=N^(n+myI) - simplify(solve(CelineRec2,N^(n+myI)));
#myCoeff:=CelineRec2:
myCoeff:=simplify(myCoeff*denom(simplify(solve(CelineRec2,N^(n+myI))))):
if(myCoeff=0) then
  FAIL:
else
  myCoeff:
fi:
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:

#Obsolete, apparently maple just does this automatically
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:

#maple won't evaluate degree without this for some reason
foodeg:=proc(expr,x)
degree(expr,x):
end:

RecToSeq:=proc(Rec,N,Ini,n,numterms) local d,i,ret,myRec,testrec,degree,myexpr:
myRec:=simplify(Rec/N^n):
myRec:=simplify(myRec/N^(ldegree(myRec,N))):
d:=foodeg(myRec,N):
if(nops(Ini)<d) then
  print(`give longer Ini`):
  RETURN(FAIL):
fi:
ret:=Ini:
i:=1:
while(nops(ret)<numterms) do
  myexpr:=subs(n=i-1,myRec):
  myexpr:=solve(myexpr,N^d):
  myexpr:=subs({seq(N^j=ret[i+j],j=0..d-1)},myexpr):
  ret:=[op(ret),myexpr]:
  i:=i+1:
od:
RETURN(ret):
end:

QuickCheckSystem:=proc(eqs,vars,NumCheck,plugvars) local i,myEqs,eq,sol:
for i from 1 to NumCheck do
myEqs:=subs({seq(plugvars[j] = 1/(i*nops(plugvars)+j),j=1..nops(plugvars))},eqs):
sol:=solve(myEqs,vars):
if(subs(sol,vars)={0}) then
RETURN(false):
fi:
od:
RETURN(true):
end:

SumToSeq:=proc(numTerms,d,Tini,eqT,F,n,k,N,T) local i,j,ii:
seq(
add(
subs({seq(T[1-i] = Tini[nops(Tini)-i],i=0..nops(Tini)-1)},
ReduceT(T,eqT,N,j,1,n))^d*subs({k=j,n=ii},F)
,j=0..ii),ii=1..numTerms):
end: