######################################################################### 
# ContMarkovWZ: save this file as ContMarkovgWZ. To use it, stay in the##
# same directory, get into Maple (by typing: maple <Enter> )           ##
# and then type:  read ContMarkovWZ : <Enter>                          ##
# Then follow the instructions given there                             ##
#                                                                      ## 
# Written by Mohamud Mohammed and Doron Zeilberger,                    ##
# Rutgers University ,                                                 ## 
#########################################################################


print():
print(`This is ContMarkovWZ, A Maple package`):
print(` by M. Mohammed and D. Zeilberger`):
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg .`):
print(`Please report all bugs to: zeilberg at math `):
print(`dot rutgers dot edu.`): 
print():
print(`For general help, and a list of the available functions,`): 
print(`type "ezra();". For specific help type "ezra(procedure_name)" `):
print():

ezra:=proc()

if args=NULL then
print(`ContMarkovWZ `):
print(`A Maple package for discovering Cont. Markov-WZ Pairs. `):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`For a list of all procedures, type ezra(); `):
print(`The Main  procedure is: `):
print(`  Markov `):

print():

elif nargs=1 and args[1]=Gosper1 then
print(`Gosper1(POL,F,k): Given a polynomial POL(k), and a hypergeometric`):
print(`F(k), and a variable k, decides whther there exist a hypergeometric`): 
print(`z(k) such that z(k+1)-z(k)=POL(k)*F(k),or FAIL.`):
print(`For example, try: Gosper(1,k*k!,k)`):

elif nargs=1 and args[1]=CheckZeil then
print(` Cheks if F and G satifiey (d/dx)F-(G(k+1)-G(k))=0`):

elif nargs=1 and args[1]=GosperP then
print(` GosperP(POL,F,k,VAR1): Given a polynomial POL(k), `):
print(` featuring coefficients from a list VAR1, and a hypergeometric`):
print(` F(k), and a variable k, finds the general assignment for VAR1 `):
print(`that will be it Gosperable together with the the anti-difference`): 
print(` z(k) divided by F `):
print(`z(k) such that z(k+1)-z(k)=POL(k)*F(k)`):
print(`For example, try: GosperP(a0+a1*k,k!,k,[a0,a1]); `):

elif nargs=1 and args[1]=GosperP1 then
print(` GosperP1(POL,F,k,VAR1): Like GosperP `):
print(`but only returns the solution set and the the z(k) `):
print(`For example, try: GosperP1(a0+a1*k,k!,k,[a0,a1]); `):


elif nargs=1 and args[1]=Markov then
print(` Markov(H,k,x,para): Given a hypergemetric term H in discrete `):
print(` variables k and continous variable x, returns the WZ-Markov pair.`):
print(`The first component is the F part in the WZ-pair`):
print(`The second component is a transition matrix A such that `):
print(` [(d/dx)a[0](x), (d/dx)a[1](x), ... ]= A times` ):
print(` [a[0](x),a[1](x), ... ] `):
print(`The third component is the coefficients of POL `):
print(`The last component is the G part `):
print(`i.e. we have d/dx(F(x,k)POL(x,k)-F(n,k)*POl(x,k)=G(n,k+1)-G(n,k) `):
print(`Try: Markov((-1)^k*x^(2*k+1)/(2*k+1)!,k,x,k)`):

else
 print(`There is no Ezra for: `, convert( args, string ) ):
fi:

end:

#Gosper1(POL,F,k): Given a polynomial POL(k), and a hypergeometric 
#F(k),and a variable k, decides whther there exist a hypergeometric 
#z(k) suchthat z(k+1)-z(k)=POL(k)*F(k)
Gosper1:=proc(POL,F,k)
local r,i,a,b,c,deg,x,eq,var,e,lu:
r:=normal(expand(simplify(subs(k=k+1,F)/F))):
a:=factor(numer(r)): b:=factor(denom(r)):
c:=POL:


lu:=MakeBest(a,b,c,k):
a:=lu[1]: b:=lu[2]:c:=lu[3]:


deg:=FindDeg(a,b,c,k):
if deg=FAIL then
 RETURN(FAIL):
fi:

x:=GenPol(k,deg,e):
var:=x[2]:
x:=x[1]:

lu:=expand(a*subs(k=k+1,x)-subs(k=k-1,b)*x-c):
eq:={seq(coeff(lu,k,i),i=0..degree(lu,k))}:
var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

x:=subs(var,x):
normal(subs(k=k-1,b)*POL*x/c):
end:


#FindDeg(a,b,c,k): Finds the degree for x(k)
FindDeg:=proc(a,b,c,k)
local D1,D2,i1,A,B:
if degree(a,k)<>degree(b,k) or
   coeff(a,k,degree(a,k))<>coeff(b,k,degree(b,k)) then
    RETURN(degree(c,k)-max(degree(a,k),degree(b,k))):
fi:

A:=coeff(a,k,degree(a,k)-1):
B:=coeff(b,k,degree(b,k)-1):

D1:={degree(c,k)-degree(a,k)+1,(B-A)/coeff(a,k,degree(a,k))}:
D2:={}:
for i1 from 1 to nops(D1) do
 if type(D1[i1],integer) and D1[i1] >=0 then
    D2:=D2 union {D1[i1]}:
 fi:
od:
if D2={} then
  RETURN(FAIL):
else
 RETURN(max(op(D2))):
fi:
end:


#GenPol(n,deg,a): A generic polynomial of degree deg in n
#with coeffs. indexed by a, followed by the set of coeffs.
GenPol:=proc(n,deg,a) local i:
convert([seq(a[i]*n^i,i=0..deg)],`+`), {seq(a[i],i=0..deg)}:
end:

#GosperP(POL,F,k,VAR1): Given a polynomial POL(k),
#featuring coefficients from a list VAR1, and a hypergeometric
#F(k), and a variable k, finds the general assignment for VAR1
#that will be it Gosperable together with the the anti-difference
#z(k) divided by F
#z(k) such that z(k+1)-z(k)=POL(k)*F(k)
GosperP:=proc(POL,F,k,VAR1)
local r,i,a,b,c,deg,x,eq,var,e,lu,var1,tov,pu:
r:=normal(expand(simplify(subs(k=k+1,F)/F))):
a:=factor(numer(r)): b:=factor(denom(r)):

c:=POL:


lu:=MakeBest(a,b,c,k):
a:=lu[1]: b:=lu[2]:c:=lu[3]:

deg:=FindDeg(a,b,c,k):
if deg=FAIL then
 RETURN(FAIL):
fi:

x:=GenPol(k,deg,e):
var:=x[2]:
x:=x[1]:

lu:=expand(a*subs(k=k+1,x)-subs(k=k-1,b)*x-c):
eq:={seq(coeff(lu,k,i),i=0..degree(lu,k))}:
var1:=solve(eq,var union convert(VAR1,set)):

if var1=NULL then
 RETURN(FAIL):
fi:

tov:={}:
for i from 1 to nops(var1) do
 if op(1,var1[i])=op(2,var1[i])    then
   tov:=tov union {op(1,var1[i])}:
  fi:
od:

x:=subs(var1,x):
pu:=[seq([VAR1[i],subs(var1,VAR1[i])],i=1..nops(VAR1))]:
normal(subs(k=k-1,b)*POL*x/c),pu, tov,var1:
end:


#Gormim(P,k): all the roots a of P(k)=a such that in
#the fact
#it has non-linear factors
Gormim:=proc(P,k) local gu,mu,i,gu1:

if degree(P,k)=0 then
 RETURN({}):
fi:

mu:={}:
if degree(P,k)=1 then
RETURN({solve(P=0,k)}):
fi:

gu:=factor(P):

if type(gu,`*`) and degree(op(1,gu),k)=0 then
  gu:=factor(gu/op(1,gu)):
fi:

if type(gu,`^`) then
 gu1:=op(1,gu):
 if degree(gu1,k)<>1 then
   RETURN(FAIL):
 else
   RETURN({solve(gu1,k)}):
 fi:
fi:

for i from 1 to nops(gu) do
gu1:=op(i,gu):
 if type(gu1,`^`) then
  gu1:=op(1,gu1):
 fi:

 if degree(gu1,k)<>1 then
    RETURN(FAIL):
  else
    mu:=mu union {solve(gu1,k)}:
 fi:
od:
mu:
end:


Findh:=proc(a,b,k) local h,g,gu1,gu2,lu,i,j,lu1:

gu1:=Gormim(a,k):
gu2:=Gormim(b,k):
if gu1<>FAIL and gu2<>FAIL then
lu:={seq(seq(gu2[i]-gu1[j],j=1..nops(gu1)),i=1..nops(gu2))}:
 for i from 1 to nops(lu) do
  lu1:=op(i,lu):
   if type(lu1,integer) and lu1>0 then
     RETURN(lu1):
   fi:
 od:
  RETURN(FAIL):
else
print(`Using empirical approach`):
fi:
for h from 1 to 100 do
g:=gcd(a,subs(k=k+h,b)):
if degree(g,k)>0 then
RETURN(h):
fi:
od:
FAIL:
end:

MakeBetter:=proc(a,b,c,k) local h,g,a1,b1,c1,i1:
a1:=a: b1:=b:c1:=c:

h:=Findh(a,b,k):

if h<>FAIL then
g:=gcd(a1,subs(k=k+h,b1)):
 c1:=c1*convert([seq(subs(k=k-i1,g),i1=1..h)],`*`):
 a1:=normal(a1/g):
 b1:=normal(b1/subs(k=k-h,g)):
fi:

a1,b1,c1:
end:

MakeBest:=proc(a,b,c,k) local lu,lu1:
lu:=a,b,c:
lu1:=MakeBetter(lu,k):

while lu1<>lu do
lu:=lu1:
lu1:=MakeBetter(lu1,k):
od:
lu1:
end:


#GosperP1(POL,F,k,VAR1): Given a polynomial POL(k),
#featuring coefficients from a list VAR1, and a hypergeometric
#F(k), and a variable k, finds the general assignment for VAR1
#that will be it Gosperable together with the the anti-difference
#z(k) divided by F
#z(k) such that z(k+1)-z(k)=POL(k)*F(k)
GosperP1:=proc(POL,F,k,VAR1)
local r,i,a,b,c,deg,x,eq,var,e,lu,var1:
r:=normal(expand(simplify(subs(k=k+1,F)/F))):
a:=factor(numer(r)): b:=factor(denom(r)):

c:=POL:


lu:=MakeBest(a,b,c,k):
a:=lu[1]: b:=lu[2]:c:=lu[3]:

deg:=FindDeg(a,b,c,k):
if deg=FAIL then
 RETURN(FAIL):
fi:

x:=GenPol(k,deg,e):
var:=x[2]:
x:=x[1]:

lu:=expand(a*subs(k=k+1,x)-subs(k=k-1,b)*x-c):
eq:={seq(coeff(lu,k,i),i=0..degree(lu,k))}:
var1:=solve(eq,var union convert(VAR1,set)):

if var1=NULL then
  RETURN(FAIL):
fi:

x:=subs(var1,x):


var1,factor(normal(subs(k=k-1,b)*POL*x/c)) :

end:

#CheckZeil(F,x,z,cert): checks the output to Markov

CheckZeil:=proc(F,x,z,cert) local gu:

gu:=diff(F,x)-subs(k=k+1,F*cert)+F*cert:

simplify(normal(expand(gu))):

end:


Markov1:=proc(F,k,n,deg,para)
local gu,POL,POL1,F1,a,b,Var2, i,lu,i1,mu,mu1,cert:
option remember:

POL:=convert([seq(a[i]*para^i,i=0..deg)], `+`):

gu:=(diff(F,n)/F*POL)+convert([seq(b[i]*para^i,i=0..deg)], `+`):

gu:=simplify(normal(expand(gu))):

POL1:=numer(gu):
F1:=F/denom(gu):
Var2:=[seq(b[i],i=0..deg)]:

lu:=GosperP1(POL1,F1,k,Var2):
if lu=FAIL then
 RETURN(FAIL):
fi:
cert:=lu[2]/denom(gu):

lu:=lu[1]:
mu:=[]:

for i from 0 to deg do
mu1:=expand(subs(lu,b[i])):
mu:=[op(mu),[seq(factor(normal(coeff(mu1,a[i1],1))),i1=0..deg)]]: 
od:
POL1:=subs(lu,POL):

POL1:=[seq(coeff(POL1,a[i1],1),i1=0..deg)]:
cert:=[seq(coeff(cert,a[i1],1),i1=0..deg)]:

[F,mu,POL1,cert]:

end:


Markov:=proc(F,k,n,para) 
local deg,gu:

for deg from 0 do

gu:=Markov1(F,k,n,deg,para):

if gu<>FAIL then
 RETURN(gu):
 fi:
 
 od:
end: