#######################################################################
## GuessRat: Save this file as GuessRat. To use it, stay in the       #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read GuessRat: <Enter>                             #
## Then follow the instructions given there                           #
##                                                                    #
## Written by Drew Sills and                                          #
## Doron Zeilberger, Rutgers University ,                             #
##  sills at math dot rutgers dot edu ,                               #
##  zeilberg at math dot rutgers dot edu.                             # 
#######################################################################


print(`First Written: Feb. 19 , 2004: tested for Maple 8 and 9 `):
print(`Version of Feb. 19, 2004:  `):
print():
print(`This is GuessRat, A Maple package`):
print(`that guesses explicit expressions, as rational functions of `):
print(`several variables, for empirically derived functions `):
print(``):
print(`It accompanyies the article: Variations on Dyson's Ex-Conjecture `):
print(` by Drew Sills and Doron  Zeilberger`):
print(`But is of very general Applicability. `):
print():
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 nargs=0 then
print(` GuessRat `):
print(` A Maple package for `):
print(` help with a specific procedure, type "ezra(procedure_name);"`):
print(`The MAIN procedure is: GuessRat `):
print(` the supporting procedures are : GenPol, GenRat, `):
print():

elif nargs=1 and args[1]=GenPol then
print(`GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a,`):
print(`and a non-negative integer deg, outputs a generic polynomial`):
print(`of total degree deg,in the variables Resh, indexed by the letter`):
print(` followed by the set of coeffs. For example, try `):
print(` GenPol([x,y],a,1); `):


elif nargs=1 and args[1]=GenRat then
print(`GenRat(Resh,a,b,degT,degB): The generic form of a rational function`):
print(`in the list of variables Resh indexed by a on top and b at bottom`):
print(`of total degree degT on top and degB at bottom`):
print(`for example, try:  GenRat([x,y],a,b,1,2)`):

elif nargs=1 and args[1]=GuessRat then
print(`GuessRat(F,Resh,deg,M): Given a function or procedure R with`):
print(`nops(Resh)  arguments (that are integers), guesses a rational`):
print(` function, using data that starts at M`):
print(`with the sum of total degrees of top and bottom <=deg`):
print(`that seems fits it, or returns FAIL`):
print(`For example, try: GuessRat(a->a/(a+1),[a],2,1);`):


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


#GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a,
#and a non-negative integer deg, outputs a generic polynomial
#of total degree deg,in the variables Resh, indexed by the letter a of 
#followed by the set of coeffs. For example, try:
#GenPol([x,y],a,1);
GenPol:=proc(Resh,a,deg)
local var,POL1,i,x,Resh1,POL:

if not type(Resh,list) or nops(Resh)=0 then
 ERROR(`Bad Input`):
fi:

x:=Resh[1]:
if nops(Resh)=1 then
RETURN(convert([seq(a[i]*x^i,i=0..deg)],`+`),{seq(a[i],i=0..deg)}):
fi:

Resh1:=[op(2..nops(Resh),Resh)]:
var:={}:
POL:=0:

for i from 0 to deg do
POL1:=GenPol(Resh1,a[i],deg-i):
var:=var union POL1[2]:
POL:=expand(POL+POL1[1]*x^i):
od:
POL,var:
end:

#GenRat(Resh,a,b,degT,degB): The generic form of a rational function
#in the list of variables Resh indexed by a on top and b at bottom
#of total degree degT on top and degB at bottom
#for example, try:  GenRat([x,y],a,b,1,2)
GenRat:=proc(Resh,a,b,degT,degB)
local TOP,BOT:
TOP:=GenPol(Resh,a,degT):
BOT:=GenPol(Resh,b,degB):
TOP[1]/BOT[1],TOP[2] union BOT[2]:
end:


#Simp(k,L): all the weakly-decreasing k-vectors of non-negative integers with 
#entries<=L, for example Simp(2,4);
Simp:=proc(k,L)
local gu,mu,i,j,i1:
if k=1 then
RETURN({seq([i],i=0..L)}):
fi:
gu:={}:

for i from 0 to L do
mu:=Simp(k-1,i):
 gu:=gu union {seq(seq([op(mu[i1]),j],j=i..L),i1=1..nops(mu))}:
od:
gu:
end:

#Simp1(k,M,L): all the weakly-decreasing k-vectors of non-negative 
#integers with 
#M<=entries<=L, for example Simp(2,3,4);
Simp1:=proc(k,M,L)
local gu,mu,i,j,i1:
if k=1 then
RETURN({seq([i],i=M..L)}):
fi:
gu:={}:

for i from M to L do
mu:=Simp1(k-1,M,i):
 gu:=gu union {seq(seq([op(mu[i1]),j],j=i..L),i1=1..nops(mu))}:
od:
gu:
end:


#GuessRat1(F,Resh,degT,degB): Given a function or procedure R with
#nops(Resh)  arguments (that are integers), guesses a rational function
#with total degree degT and degB for numerator and denom. resp.
#that seems fits it, using data of integer k-tuples that start at M
GuessRat1:=proc(F,Resh,degT,degB,M)
local eq,var,Rat,a,b,i1,j1,mu,k,L:
k:=nops(Resh):
Rat:=GenRat(Resh,a,b,degT,degB):
var:=Rat[2]:
Rat:=Rat[1]:
eq:={}:

for L from M while nops(Simp1(k,M,L))<nops(var)+10 do
od:
mu:=Simp1(k,M,L):

eq:=
{seq(
numer(normal(F(op(mu[i1]))-subs({seq(Resh[j1]=mu[i1][j1],j1=1..k)},Rat))),
i1=1..nops(mu))}:

eq:=eq minus {0}:

if nops(eq)<nops(var)+5 then
print(`Too many zeroes, try to raise L `):
RETURN(FAIL):
fi:

var:=solve(eq,var):
if subs(var,denom(Rat))=0 or var=NULL then
  RETURN(FAIL):
else
 RETURN(normal(subs(var,Rat))):
fi:


end:


#GuessRat(F,Resh,deg): Given a function or procedure R with
#nops(Resh)  arguments (that are integers), guesses a rational function
#with total degrees of sums of degrees of top and bottom up to
#deg that seems fits it, or  RETURNS FAIL
#For example, try: GuessRat1((a,b)->a/(a+b+1),[a,b],2);

GuessRat:=proc(F,Resh,deg,M)
local lu,L:

for L from 0 to deg do
lu:=GuessRat1(F,Resh,L,deg-L,M):
if lu<>FAIL then
  RETURN(lu):
fi:
od:
FAIL:
end:



Dyson3:=proc(b1,b2,a1,a2,a3)
local k,L:
L:=max(abs(b1),abs(b2)):
(a1,a2,a3)->
convert([seq((-1)^k*binomial(a1+a2,a2+k)*binomial(a1+a3,a1+k-b1)*
binomial(a2+a3,a3+k+b2),k=-max(a1,a2,a3)-L..max(a1,a2,a3)+L)],`+`)/
((a1+a2+a3)!/(a1!*a2!*a3!)):
end:

Dyson3a:=proc(b1,b2,a1,a2,a3)
local k,L:
L:=max(abs(b1),abs(b2)):
convert([seq((-1)^k*binomial(a1+a2,a2+k)*binomial(a1+a3,a1+k-b1)*
binomial(a2+a3,a3+k+b2),k=-max(a1,a2,a3)-L..max(a1,a2,a3)+L)],`+`)/
((a1+a2+a3)!/(a1!*a2!*a3!)):
end:


#VerifyMN(Rat,Resh):Given a rational function Rat in the variables
#in the list Resh, verifies that Rat*MultiNomial(Resh) satisfies the
#same the famous Pascal-Like recurrence of the Mult-nomial coeff.
VerifyMN:=proc(Rat,Resh) local i,su:
su:=convert(Resh,`+`):
normal(su*Rat-
convert([seq(subs(Resh[i]=Resh[i]-1,Rat)*Resh[i],i=1..nops(Resh))],`+`)):
end: