######################################################################
##qFindRec.txt: Save this file as  qFindRec.txt                      #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  qFindRec.txt <Enter>                         #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: March , 2017
 
print(`Created:  March 2017`):
print(` This is qFindRec.txt `):
print(`That guesses linear recurrences with coefficients that are polynomials in q^n and q for sequencds of polynomials in q`):
print(`In other words, the 1D q-holonomic ansatz, or q-P-recursive ansatz. `):

print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: qbin, qfac, qFR1 `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: qFR, SumSeq, SumSeqBi `):
 print(` `):



elif nops([args])=1 and op(1,[args])=qbin then
print(`qbin(a1,b1,q): the q-binomial  coefficient for numeric a1 and b1. Try:`):
print(` qbin(6,3,q); `):

elif nops([args])=1 and op(1,[args])=qfac then
print(`qfac(n1,q): (1-q)*...*(1-q^n1)/(1-q)^n1. Try:`):
print(`qfac(5,q);`):


elif nops([args])=1 and op(1,[args])=qFR then
print(`qFR(L,q,n,N,C): inputs a list L of polynomials (or rational functions) in q, and non-negative integers`):
print(` tries to find an operator ope(q,q^n,N) of degree deg1 and order ord1 such that deg1+ord1<=C`):
print(`annihilating L. or returns FAIL. L must be at least of length (C+3)^2/4+4.  Try`):
print(`qFR([seq(q^(i1^2),i1=1..15)],q,n,N,3);`):

elif nops([args])=1 and op(1,[args])=qFR1 then
print(`qFR1(L,deg1,ord1,q,n,N): inputs a list L of polynomials (or rational functions) in q, and non-negative integers`):
print(`deg1 and ord1 tries to find an operator ope(q,q^n,N) annihilating L. Try:`):
print(` qFR1([seq(q^(i1^2),i1=1..20)],2,1,q,n,N); `):


elif nops([args])=1 and op(1,[args])=SumSeq then
print(`SumSeq(F,N0): inputs  a summand F phrased in terms of the qbinomial coefficients (called qbin), or q-factorial (called qfac),  that`):
print(`depends on k and n, outputs the first N0 terms (starting at n=1) of the Sum(F,k=0..n). Try: `):
print(`SumSeq((n,k)->q^(k*(k-1)/2)*qbin(n,k,q),10); `):
print(`SumSeq((n,k)->q^(k^2)/qfac(k,q)/qfac(n-k,q),10); `):

elif nops([args])=1 and op(1,[args])=SumSeqBi then
print(`SumSeqBi(F,N0): inputs  a summand F phrased in terms of the qbinomial coefficients (called qbin) or q-factorial (called qfac) that`):
print(`depends on k and n, outputs the first N0 terms (starting at n=1) of the Sum(F,k=-n..n). Try: `):
print(`SumSeqBi((n,k)->(-1)^k*q^((5*k^2+k)/2)/qfac(n-k,q)/qfac(n+k,q),10); `):


else
print(`There is no ezra for`,args):
fi:
 
end:

#qFR1(L,deg1,ord1,q,n,N): inputs a list L of polynomials (or rational functions) in q, and non-negative integers
#deg1 and ord1 tries to find an operator ope(q,q^n,N) annihilating L. Try
#qFR1([seq(q^(i1^2),i1=1..10)],1,2,q,n,N);
qFR1:=proc(L,deg1,ord1,q,n,N) local a,ope,var,i,j,eq,eq1,hal,i1,n1,lu,hal1,lu1:

if nops(L)<= (ord1+1)*(deg1+2)+3  then
 print(`The list is too short, it has `, nops(L), `terms but it should have at least`, (ord1+1)*(deg1+2)+3, `terms .` ):
   RETURN(FAIL):
fi:

ope:=add(add(a[i,j]*(q^(n*i))*N^j,i=0..deg1),j=0..ord1):

var:={seq(seq(a[i,j],i=0..deg1),j=0..ord1)}:

eq:=expand({seq(add(subs(n=n1,coeff(ope,N,i1))*L[n1+i1],i1=0..ord1),n1=1..(ord1+1)*(deg1+1)+3  ) }):

eq1:=subs(q=2,eq):

hal:=solve(eq1,var):

if expand(subs(hal,ope))=0 then
  RETURN(FAIL):
else
print(`There is hope for a recurrence of order`, ord1, `and degree`, deg1):
fi:


hal:=solve(eq,var):
if expand(subs(hal,ope))=0 then
  RETURN(FAIL):
fi:

lu:={}:

for hal1 in hal do
if op(1,hal1)=op(2,hal1) then
 lu:=lu union {op(1,hal1)}:
fi:
od:


if nops(lu)>1 then
 print(`Non unique solution`):
 RETURN({seq(coeff(ope,lu1,1),lu1 in lu)}):
fi:

ope:=subs(hal,ope):

ope:=subs(lu[1]=1,ope):

if normal(expand({seq(add(subs(n=n1,coeff(ope,N,i1))*L[n1+i1],i1=0..ord1),n1=(ord1+1)*(deg1+1)+3..nops(L)-ord1   ) }))<>{0}  then
 print(ope, `did not work out`):
 RETURN(FAIL):
fi:

add(factor(coeff(ope,N,i1))*N^i1,i1=0..ord1):


end:




#qFR(L,q,n,N,C): inputs a list L of polynomials (or rational functions) in q, and non-negative integers
# tries to find an operator ope(q,q^n,N) of degree deg1 and order ord1 such that deg1+ord1<=C
#annihilating L. or returns FAIL. L must be at least of length (C+3)^2/4+4.  Try
#qFR([seq(q^(i1^2),i1=1..10)],q,n,N,3);
qFR:=proc(L,q,n,N,C) local ope, ord1,deg1,C1:

if nops(L)< (C+3)^2/4+4 then
 print(`L must be at least of length`,  trunc((C+3)^2/4+4)):
 RETURN(FAIL):
fi:

for C1 from 1 to C do
 for ord1 from 1 to C1 do
    deg1:=C1-ord1:
    ope:=qFR1(L,deg1,ord1,q,n,N):
    if ope<>FAIL then
      RETURN(ope):
   fi:
 od:
od:

FAIL:

end:

#qfac(n1,q): (1-q)*...*(1-q^n1)/(1-q)^n1. Try:
#qfac(5,q);
qfac:=proc(n1,q) local i:
expand(mul(1-q^i,i=1..n1)):
end:


#qbin(a1,b1,q): the q-binomial  coefficient for numeric a1 and b1. try:
#qbin(6,3,q);
qbin:=proc(a1,b1,q)  local i:
expand(normal(mul(1-q^(a1-i+1),i=1..b1)/mul(1-q^i,i=1..b1))):
end:

#SumSeq(F,N0): inputs  a summand F phrased in terms of the qbinomial coefficients (called qbin) or q-factorial (called qfac) that
#depends on k and n, outputs the first N0 terms (starting at n=1) of the Sum(F,k=0..n). Try
#SumSeq((n,k)->q^(k*(k-1)/2)*qbin(n,k,q),10);
SumSeq:=proc(F,N0) local n1,k1:
normal([seq(expand(add(F(n1,k1),k1=0..n1)),n1=1..N0)]):
end:


#SumSeqBi(F,N0): inputs  a summand F phrased in terms of the qbinomial coefficients (called qbin) or q-factorial (called qfac) that
#depends on k and n, outputs the first N0 terms (starting at n=1) of the Sum(F,k=-n..n). Try
#SumSeqBi((n,k)->(-1)^k*q^((5*k^2+k)/2)/qfac(n-k,q)/qfac(n+k,q),10);
SumSeqBi:=proc(F,N0) local n1,k1:
normal([seq(expand(add(F(n1,k1),k1=-n1..n1)),n1=1..N0)]):
end: