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

print(`First Written: May 26, 2004: tested for Maple 8 and 9 `):
print(`Version of May 26, 2004:  `):
print():
print(`This is MarkovAZ, One of the Maple packages`):
print(`accompanying the article: The  Markov-WZ Method. `):
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 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(`MarkovAZ: `):
print(`A Maple package for discovering Markov-Almkvist-Zeilberger Pairs. `):
print(`written by Mohamud Mohammed and Doron Zeilberger. `):
print(`It is one of the Maple packages accompaanying their article`):
print(`The Markov-WZ Method. `):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`For a list of all procedures, type ezra1(); `):
print(`The procedures are: `):
print(` AZd,AZG,AZGwP,MarkovAZ`):
print(` `):
print():

elif nargs=1 and args[1]=AZd then
print(`AZd(F,x,n,N): The Almkvist-Zeilberger algorithm, `):
print(`also found in the Maple package EKHAD. `):
print(`Given an expression F in the cont. variable`):
print(` x and the discrete variable n, and a letter N (denoting`):
print(` the shift operator in n: N g(n):=g(n+1)) `):
print(` that is hypergeometric, `):
print(`i.e. F(x,n+1)/F(x,n) and (D_xF)/F are rational functions of`):
print(` of n and x , finds a linear recurrence opeator ope(N,n) `):
print(`and a rational function R(n,x) such that `):
print(`ope(N,n)F(n,x)=D_x (R(n,x) F(n,x) ) `):
print(`For example, try: AZd((1+x)^(2*n)/x^(n+1),x,n,N); `):

elif nargs=1 and args[1]=AZG then
print(`AZG(A,POL,y): Given a hyperexponential A(y) and a polynomial P(y)`):
print(`decides whether there exists a hyperexponetial S(y) such`):
print(`that S'(y)=A(y)P(y). If it fails, it returns FAIL, otherwise`):
print(`it returns the rational function Rat(y) such that`):
print(`the derivative of (A(y)*Rat(y)) w.r.t. y equals POL(y)*A(y)`):
print(`For example, try: AZG(exp(x),1,x);`):

elif nargs=1 and args[1]=AZGwP then
print(`AZGwP(A,P,y,Pars): Given a hyperexponential A(y)`):
print(`and a polynomial P(y) depending linearly on the list of`):
print(` parameters in Pars decides whether there exists a hyperexponetial `):
print(` S(y) such that S'(y)=A(y)P(y). If it fails, it returns FAIL, `):
print(`otherwise it returns the rational function Rat(y) such that`):
print(`(A(y)*Rat(y))' equals POL(y)*A(y), followed by the succesful POL`):
print(`For example, try: AZGwP(exp(x^2),a+b*x,x,[a,b]); `):

elif nargs=1 and args[1]=MarkovAZ then
print(`MarkovAZ(H,x,n): An AZ-Markov Pair with kernel H(x,n) where`):
print(`x is a continuous variable and n is a discrete variable`):
print(`and H(x,n) is hyper-exp.-geom. For example, try:`):
print(`The output consists of four parts`):
print(`The first part is the kernel H itself`):
print(`The second part is the transition matrix`):
print(`between a_0(n), a_1(n), ..., a_d(n) `):
print(` and a_0(n+1), a_1(n+1), ..., a_d(n+1) `):
print(` which is a matrix of rational functions `):
print(`The third part is always [1,x,x^2, ..., x^d]`):
print(`reminding us that to get the F-part of the MAZ pair`):
print(`one multiplies the kernel F(x,n) by the dot product of`):
print(` [a_0(n), ..., a_d(n)] and [1,x,x^2, ..., x^d ] `):
print(`The fourth (and last) part is a vector of length d+1 `):
print(`of rational functions in x and d`):
print(`from which the G-part of the MAZ pair is formed by `):
print(`taking the dot product of it with`):
print(` [a_0(n), ..., a_d(n)]`):
print(`Recall that (F(x,n),G(x,n)) is a MAZ-pair if F is hyperexponential`):
print(`in x and hypergemetric in n, and `):
print(` We have F(x,n+1)-F(x,n)=D_x G(n,x) .`):
print(` For example, try: `):
print(`MarkovAZ(x^n*(1-x)^b,x,n); `):



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

end:



#AZG: Given a hyperexponential A(y) and a polynomial P(y)
#decides whether there exists a hyperexponetial S(y) such
#that S'(y)=A(y)P(y). If it fails, it returns FAIL, otherwise
#it returns the rational function Rat(y) such that
#(A(y)*Rat(y))' equals POL(y)*A(y)
AZG:= proc(A,POL,y) local Q,R,P,KAMA,yakhas,j1,g,Q1,l1,l2,k1,mumu,fu,
Q2,meka1,mekb1,gugu,degg,eq,var,P1,i: 
KAMA:=20:
yakhas:=normal(diff(A,y)/A):
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
P:=1:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:

P1:=POL*P:

P1:=expand(P1):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P1,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P1,y)-l2):
  else
   k1:=degree(P1,y)-l2:
  fi:
fi:

if k1 < 0 then
RETURN(FAIL):
fi:

fu:=convert([seq(a[i]*y^i,i=0..k1)],`+`):

gugu:=expand(Q1*fu+R*diff(fu,y)-P1):
degg:=degree(gugu,y):
eq:={seq(coeff(gugu,y,i),i=0..degg)} minus {0}:
var:={seq(a[i],i=0..degg)}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:
fu:=subs(var,fu):
if fu=0  then
 RETURN(FAIL):
fi:
fu:=normal(fu):
 
R*fu/P:
end:

CheckAZG:=proc(f,x)
local mu:
mu:=diff(f,x):

normal(mu*AZG(mu,1,x)):
end:


CheckAZG1:=proc(f,POL,x) local mu:
mu:=normal(diff(f*POL,x)/f):
normal(AZG(f/denom(mu),numer(mu),x)*denom(mu)),POL;
end:

#AZGwP(A,POL,y,Pars): Given a hyperexponential A(y) and a polynomial P(y)
#depending linearly on the list of parameters in Pars
#decides whether there exists a hyperexponetial S(y) such
#that S'(y)=A(y)P(y). If it fails, it returns FAIL, otherwise
#it returns the rational function Rat(y) such that
#(A(y)*Rat(y))' equals POL(y)*A(y), followed by the succesful POL
AZGwP:= proc(A,POL,y,Pars) local Q,R,P,KAMA,yakhas,j1,g,Q1,l1,l2,k1,mumu,fu,
Q2,meka1,mekb1,gugu,degg,eq,var,P1,a,i: 
KAMA:=20:
yakhas:=normal(diff(A,y)/A):
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
P:=1:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:

P1:=POL*P:

P1:=expand(P1):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P1,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P1,y)-l2):
  else
   k1:=degree(P1,y)-l2:
  fi:
fi:

if k1 < 0 then
RETURN(FAIL):
fi:

fu:=convert([seq(a[i]*y^i,i=0..k1)],`+`):

gugu:=expand(Q1*fu+R*diff(fu,y)-P1):
degg:=degree(gugu,y):
eq:={seq(coeff(gugu,y,i),i=0..degg)} minus {0}:
var:={seq(a[i],i=0..degg)}:
var:=solve(eq,var union {op(Pars)}):
if var=NULL then
 RETURN(FAIL):
fi:

fu:=subs(var,fu):
if fu=0  then
 RETURN(FAIL):
fi:


fu:=normal(fu):
 
R*fu/P, subs(var,POL), subs(var,Pars),Pars,var:
end:





AZd1:=proc(F,x,n,N,ORDER)
local gu,var,i,gu2,pu,POL,ope,cert,a:
var:=[seq(a[i],i=0..ORDER)]:
gu:=0:
for i from 0 to ORDER do
gu:=gu+a[i]*normal(expand(simplify(subs(n=n+i,F)/F))):
od:

POL:=numer(gu):
gu2:=denom(gu):

pu:=AZGwP(F/gu2,POL,x,var):


if pu=FAIL then
 RETURN(0):
fi:
if pu[1]=0 then
 RETURN(0):
fi:


ope:=convert([seq(factor(pu[3][i])*N^(i-1),i=1..ORDER+1)],`+`):
cert:=normal(pu[1]/gu2):
ope:=normal(ope):
gu:=denom(ope):
cert:=normal(cert*gu):
ope:=numer(ope):

ope:=factor(ope):

ope:=yafe(ope,N):

ope,cert:
end:






GetRid:=proc(p)
local p1,i,q:
p1:=factor(p):
q:=p:

for i from 1 to nops(p1) do
if type(op(i,p1),indexed) then
  q:=normal(q/op(i,p1)):
fi:
od:
q:
end:


yafe:=proc(ope,N) local i:
sort(convert([seq(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N))],`+`)):end:


AZd:=proc(F,x,n,N) local i,gu:
for i from 0 do
gu:=AZd1(F,x,n,N,i):
if gu[1]<>0 then
 RETURN(yafe(GetRid(gu[1]),N),factor(GetRid(gu[2]))):
fi:
od:
end:






#AZGwP1(A,POL,y,Pars): Given a hyperexponential A(y) and a polynomial P(y)
#depending linearly on the list of parameters in Pars
#decides whether there exists a hyperexponetial S(y) such
#that S'(y)=A(y)P(y). If it fails, it returns FAIL, otherwise
#it returns the rational function Rat(y) such that
#(A(y)*Rat(y))' equals POL(y)*A(y), followed by the succesful POL
AZGwP1:= proc(A,POL,y,Pars) local Q,R,P,KAMA,yakhas,j1,g,Q1,l1,l2,k1,mumu,fu,
Q2,meka1,mekb1,gugu,degg,eq,var,P1,a,i: 
KAMA:=20:
yakhas:=normal(diff(A,y)/A):
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
P:=1:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:

P1:=POL*P:

P1:=expand(P1):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P1,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P1,y)-l2):
  else
   k1:=degree(P1,y)-l2:
  fi:
fi:

if k1 < 0 then
RETURN(FAIL):
fi:

fu:=convert([seq(a[i]*y^i,i=0..k1)],`+`):

gugu:=expand(Q1*fu+R*diff(fu,y)-P1):
degg:=degree(gugu,y):
eq:={seq(coeff(gugu,y,i),i=0..degg)} minus {0}:
var:={seq(a[i],i=0..degg)}:
var:=solve(eq,var union {op(Pars)}):
if var=NULL then
 RETURN(FAIL):
fi:

fu:=subs(var,fu):
if fu=0  then
 RETURN(FAIL):
fi:


fu:=normal(fu):

var,R*fu/P:
end:

MarkovAZ1:=proc(F,x,n,deg)
local gu,POL,POL1,F1,a,Var2, i,lu,i1,mu,mu1,cert:
option remember:
POL:=convert([seq(a[i](n)*x^i,i=0..deg)], `+`):
gu:=normal(normal(simplify(expand(subs(n=n+1,F)/F)))*subs(n=n+1,POL)-POL):

POL1:=numer(gu):

F1:=F/denom(gu):

Var2:=[seq(a[i](n+1),i=0..deg)]:


lu:=AZGwP1(F1,POL1,x,Var2):

if lu=FAIL then
 RETURN(FAIL):
fi:


cert:=normal(lu[2]/denom(gu)):

lu:=lu[1]:

mu:=[]:

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

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

[F,mu,POL1, cert]:

end:


#MarkovAZ(F,x,n): An AZ-Markov Pair with kernel F(x,n) where
#x is a continuous variable and x is a discrete variable
#and F(x,n) is hyper-exp.-geom. For example, try:
#MarkovAZ(x^n,x,n):
MarkovAZ:=proc(F,x,n) local deg:
option remember:
for deg from 0 while MarkovAZ1(F,x,n,deg)=FAIL do od:
MarkovAZ1(F,x,n,deg):end: