######################################################################
## WZpairs.txt Save this file as WZpairs.txt  to use it,             #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `WZpairs.txt` <Enter>                        #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Robert Dohgherty-Bliss and  Doron Zeilberger,          #
#Rutgers University ,                                                #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: April 2021: tested for Maple 2018 `):
print(`Version  : April 2021:  `):
print():
print(`This is WZpairs.txt, A Maple package`):
print(`accompanying Robert Dougherty-Bliss  and Doron Zeilberger's  article: `):
print(`" Experimenting with Apery Limits and WZ pairs" `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/AperyLimits.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 


ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:`):
print(` Comps, CompsAd, GenPol, IsWZ, WZpairsN `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(`The MAIN procedures are:`):
print(` WZpairs `):

elif nargs=1 and args[1]=Comps then
print(`Comps(n,d): The set of vectors of length n with non-negative integers that add-up to d. Try:`):
print(`Comps(3,4);`):

elif nargs=1 and args[1]=CompsAd then
print(`CompsAd(n,d): The set of vectors of length n with non-negative integers that add-up to <= d. Try:`):
print(`CompsAd(3,4);`):

elif nargs=1 and args[1]=GenPol then
print(`GenPol(n,d,x,a): inputs a positive integer n, a positive integer d, a symbol x and a symbol a`):
print(`outputs a generic homog. polynomial in the n variable x[1],...,x[n] of degree d with coefficients of the form a[i]`):
print(`followed by the set of coefficients. Try:`):
print(`GenPol(4,2,x,a);`):


elif nargs=1 and args[1]=IsWZ then
print(`IsWZ(P,n,k): Is the pair P a WZ pair?`):

elif nargs=1 and args[1]=WZpairs then
print(`WZpairs(K,n,k,d): inputs a closed-form kernel K in (n,k), variables n and k, and a non-negative integer d`):
print(`outputs a base for pairs [P1,P2] such that [K*Pol1,K*Pol2] is a WZ pair, where Pol1, Pol2 are polynomials in k,n, of degree <=d.`):
print(`Try:`):
print(`WZpairs((-1)^k*n!*k!/(n+k+1)!/x^k/(1+x)^(n+1),n,k,1);`):
print(`WZpairs((-1)^(n+k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1),n,k,1);`):
print(`WZpairs((-1)^(k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1)^2/(k+1),n,k,3);`):

elif nargs=1 and args[1]=WZpairsN then
print(`WZpairsN(K,n,k,d): inputs a closed-form kernel K in (n,k), variables n and k, and a non-negative integer d`):
print(`outputs a base for pairs [P1,P2] such that [K*Pol1,K*Pol2] is a WZ pair, where P1, P2 are polynomials in k,n, of degree <=d.`):
print(`Try:`):
print(`WZpairsN((-1)^k*n!*k!/(n+k+1)!/x^k/(1+x)^(n+1),n,k,1)`):
print(`WZpairsN((-1)^(n+k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1),n,k,1)`):
print(`WZpairsN((-1)^(k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1)^2/(k+1),n,k,3)`):

else
 print(`There is no such thing as`, args):

fi:


end:


ez:=proc(): print(`  WZpairs(K,n,k,d), , WZpairsN(K,n,k,d) `): 
print(`WZpairs((-1)^k*n!*k!/(n+k+1)!/x^k/(1+x)^(n+1),n,k,1);`):
print(`WZpairs((-1)^(n+k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1),n,k,1);`):
print(`WZpairs((-1)^(k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1)^2/(k+1),n,k,3);`):
end:

#Comps(n,d): The set of vectors of length n with non-negative integers that add-up to d. Try:
#Comps(3,4);
Comps:=proc(n,d) local gu,d1,mu,mu1:
option remember:

if n=0 then
 if d=0 then
  RETURN({[]}):
else
 RETURN({}):
 fi:
fi:

gu:={}:

for d1 from 0 to d do
 mu:=Comps(n-1,d-d1):
 gu:=gu union {seq([op(mu1),d1], mu1 in mu)}:
od:
gu:
end:

#CompsAd(n,d): The set of vectors of length n with non-negative integers that add-up to <= d. Try:
#CompsAd(3,4);

CompsAd:=proc(n,d) local d1:
{seq(op(Comps(n,d1)),d1=0..d)}:
end:

#GenPol(n,d,x,a): inputs a positive integer n, a positive integer d, a symbol x and a symbol a
#outputs a generic homog. polynomial in the n variable x[1],...,x[n] of degree d with coefficients of the form a[i]
#followed by the set of coefficients. Try:
#GenPol(4,2,x,a);
GenPol:=proc(n,d,x,a) local gu,gu1,P,mu,co,i1:
gu:=CompsAd(n,d):
P:=0:
mu:={}:
co:=0:
 for gu1 in gu do
 co:=co+1:
  P:=P+a[co]*mul(x[i1]^gu1[i1],i1=1..n):
  mu:=mu union {a[co]}:
 od:
P,mu:
end:


#WZpairsN(K,n,k,d): inputs a closed-form kernel K in (n,k), variables n and k, and a non-negative integer d
#outputs a base for pairs [P1,P2] such that [K*Pol1,K*Pol2] is a WZ pair, where P1, P2 are polynomials in k,n, of degree <=d.
#Try:
#WZpairsN((-1)^k*n!*k!/(n+k+1)!/x^k/(1+x)^(n+1),n,k,1)
#WZpairsN((-1)^(n+k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1),n,k,1)
#WZpairsN((-1)^(k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1)^2/(k+1),n,k,3)
WZpairsN:=proc(K,n,k,d) local P1,P2,a,b,gu,var,eq,i,j,var1,var2,ka,kv,tovim,P:
P1:=GenPol(2,d,[n,k],a):
P2:=GenPol(2,d,[n,k],b):

var:=P1[2] union P2[2]:
P1:=P1[1]:
P2:=P2[1]:

gu:=subs(n=n+1,P1)*normal(simplify(subs(n=n+1,K)/K))-P1-subs(k=k+1,P2)*normal(simplify(subs(k=k+1,K)/K))+P2:
gu:=numer(normal(gu)):

eq:={}:

for i from 0 to degree(gu,n) do
 for j from 0 to degree(coeff(gu,n,i),k) do
  eq:=eq union {coeff(coeff(gu,n,i),k,j)}:
 od:
od:

var1:=solve(eq,var):

var2:={}:

for ka in var1 do
 if op(1,ka)=op(2,ka) then
  var2:=var2 union {op(1,ka)}:
 fi:
od:


P1:=subs(var1,P1):
P2:=subs(var1,P2):

kv:={seq([expand(coeff(P1,ka,1)),expand(coeff(P2,ka,1))],ka in var2)}:

tovim:={}:

for P in kv do
 if type(sum(P[1]*K,k),function) then
  tovim:=tovim union {P}:
 fi:
od:
 
tovim: 

end:



#WZpairs(K,n,k,d): inputs a closed-form kernel K in (n,k), variables n and k, and a non-negative integer d
#outputs a base for pairs [P1,P2] such that [K*Pol1,K*Pol2] is a WZ pair, where Pol1, Pol2 are polynomials in k,n, of degree <=d.
#Try:
#WZpairs((-1)^k*n!*k!/(n+k+1)!/x^k/(1+x)^(n+1),n,k,1);
#WZpairs((-1)^(n+k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1),n,k,1);
#WZpairs((-1)^(k)*k!^2*(n-k-1)!/(n+k+1)!/(n+1)^2/(k+1),n,k,3);
WZpairs:=proc(K,n,k,d) local gu,P,ku,gu1,ku1:
gu:=WZpairsN(K,n,k,d):


ku:={}:

for gu1 in gu do:
ku1:=[factor(gu1[1])*K,factor(gu1[2])*K]:


if IsWZ(ku1,n,k)=0 then

ku:=ku union {ku1}:

fi:
od:
ku:
end:



#IsWZ(P,n,k): Is the pair P a WZ pair?
IsWZ:=proc(P,n,k) local F,G:
F:=P[1]:
G:=P[2]:
normal(simplify(subs(n=n+1,F)/F)-1+simplify(G/F)-simplify(subs(k=k+1,G)/F)):
end: