######################################################################
##LinDiophantus: Save this file as LinDiophantus. To use it,         # 
##stay in the    same directory, get into Maple                      #  
#(by typing: maple <Enter> )                                         #
##and then type:  read LinDiophantus : <Enter>                       #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Temple University ,                   #
#zeilberg@math.temple.edu.                                           # 
#######################################################################

#Created: March 2001
#This version: March 2001
#LinDiophantus: A Maple package that finds generating functions
#representating solutions of systems of Linear Diophantine Equations
#It accompanies Doron Zeilberger's article:
#A Constant Term Method for Solving Systems of Linear Diophantine
#Equations
#Please report bugs to zeilberg@math.temple.edu


print(`LinDiophantus: A Maple package that finds generating functions`):
print(`representating solutions of systems of Linear Diophantine Equations`):
print(`It accompanies Doron Zeilberger's article:`):
print(`A Constant Term Method for Solving Systems of Linear Diophantine`):
print(` Equations `):
print(`Please report bugs to zeilberg@math.temple.edu`):

lprint(``):
print(`Created: March 2001`):
print(`This version: March 2001`):
lprint(``):
print(`Written by Doron Zeilberger, zeilberg@math.temple.edu`):
lprint(``):
print(`Please report bugs to zeilberg@math.temple.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.temple.edu/~zeilberg/`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):

ezra:=proc()
print(``):
if nops([args])=0  then
print(`contains procedures: `):
print(`  CT1, GFIx, GFx, GFxt , ToUm , OMEGA`): 
fi:

if nops([args])=1 and op(1,[args])=CT1 then
print(`CT1(f,xSet,t): Given a rational`):
print(`function f of t and a set, xSet, of  x-variables`):
print(`finds the constant term in t (i.e. the coeff. of t^0)`):
print(`of the Laurent expansion of t, viewed as a formal Laurent`):
print(`series in the x's of xSet`):
print(`For example try CT(1/(1-x*t)/(y-t),[x,y],t);`):
fi:

if nops([args])=1 and op(1,[args])=GFx then
print(`GFx(a,Elist,x): Given a list of variables, a, `):
print(`of type non-negative integer, and a list`):
print(`of linear-diophantine equations in these variables, Elist `):
print(`and a continuous indexed-variable x, outputs the`):
print(`generating function of all solutions `):
print(`where x[1]^a[1]*x[2]^a[2]*... shows up iff (a[1],a[2], ...)`):
print(`is a solution`):
print(`For example, try GFx([a,b],[a-b],x); `):
fi:

if nops([args])=1 and op(1,[args])=GFIx then
print(`GFIx(a,Elist,x): Given a list of variables, a, `):
print(`of type non-negative integer, and a list`):
print(`of linear-diophantine inequalities in these variables`):
print(`and a continuous indexed-variable x, outputs the`):
print(`generaitng function of all solutions `):
print(`where x[1]^a[1]*x[2]^a[2]*... shows up iff (a[1],a[2], ...)`):
print(`is a solution. Elist is is a list of affine`):
print(`linear expressions in the members of a and the inequalities`):
print(`are Elist[1]>=0, Elist[2]>=0, ...`):
print(`For example, type: GFIx([b1,b2,b3],[b2-2*b1,2*b3-3*b2],x);`):
fi:

if nops([args])=1 and op(1,[args])=GFxt then
print(`GFxt(a,E1,x,t): Given a list of symbols, a`):
print(`denoting generic (symbolic) non-negative integers`):
print(`and a list, E1, of affine linear expressions with integer`):
print(`coeffs. in them , and symbols x and t,`):
print(`finds the generating function:`):
print(`sum of  x[1]^a[1]*...*x[n]^a[n]*t[1]^E1[1]*...*t[r]^E1[r]`):
print(`(where n=nops(a) and r=nops(E1)), and the sum is `):
print(`over 0<=a[1]<infinity, ..., 0<=a[n]<infinity, ...,`):
print(`For example, try: GFxt([a,b,c],[a-b,a-c],x,t);`):
fi:

if nops([args])=1 and op(1,[args])=OMEGA then
print(`OMEGA(f,xSet,t): Given a rational`):
print(`function f of t and a set, xSet, of  x-variables`):
print(`Applies MacMahon's operation Omega w.r.t to the`):
print(`variable t (i.e. only takes that part of f that involves)`):
print(`non-negative powers of t`):
print(`Here f is viewed as a formal Laurent in t and a`):
print(`formal power series in the x's of xSet`):
print(`For example try OMEGA(1/(1-x*t)/(y-t),[x,y],t);`):
fi:

if nops([args])=1 and op(1,[args])=ToUm then
print(`ToUm(R,xList): Given a rational function R in the`):
print(`variables of xList`):
print(`tries to find an umbra that sends f(xList) to Constant`):
print(`term of R(xList)f(1/xList), it only treats the case where`):
print(`the denominator of R(x) factors completely over the `):
print(`rationals and there are no multiplicity and`):
print(`the degree of the numerator is less than the degree of`):
print(`of the denominator for `):
print(`if this is not the case, it returns 0`):
print(`it also avoids the case when the partial umbra depends on`):
print(`the newly introduced variable`):
fi:


end:

#GFxt(a,E1,x,t): Given a list of symbols, a
#denoting generic (symbolic) non-negative integers
#and a list, E1, of affine linear expressions with integer
#coeffs. in them , and symbols x and t,
#finds the generating function:
#sum of  x[1]^a[1]*...*x[n]^a[n]*t[1]^E1[1]*...*t[r]^E1[r]
#(where n=nops(a) and r=nops(E1)), and the sum is over 0<=a[1]<infinity, ...,
#0<=a[n]<infinity, ...,
GFxt:=proc(a,E1,x,t)
local i,n,gu,r:
gu:=1:
n:=nops(a):
r:=nops(E1):

for i from 1 to n do
gu:=gu*x[i]^a[i]:
od:

for i from 1 to r do
gu:=gu*t[i]^E1[i]:
od:

for i from 1 to n do
gu:=sum(gu,a[i]=0..infinity):
gu:=normal(gu):
od:
gu:
end:

#WhoAreYout(f,t): Given a Laurent polynomial f in t
#decides it is of the form P(t), where P is a polynomial in t
#or P(1/t)/t and returns 1 and -1 respectively
WhoAreYout:=proc(f,t)
local f1:
f1:=expand(f):
if degree(f1,t) >0 and ldegree(f1,t)>=0 then 1:
elif degree(f1,t) <0 then -1:
else 0:
fi:
end:

#WhoAreYouxt(f,x,t): Given a polynomial f, in the variables
#x and t, looks at the `normalized' polynomial (in x)
#obtained by dividing by the coeff. of x^0
#and seeing whether all the  coeffs. of the powers of x
#(that are Laurent polynomials in t) have consistently
#coeffs that are polynomials in t, or consistently coeffs.
#that are polynomials in 1/t, or none of the above
#and returns 1, -1, 0 respectively

WhoAreYouxt:=proc(f,x,t)
local f1,lu,i1,i:
f1:=expand(f):

if degree(f1,t)=0 then
RETURN(1):
fi:

if degree(f1,t)<>0 and normal(f1-subs(t=1,f1)*t^degree(f1,t))=0 then
RETURN(-1):
fi:


if coeff(f1,x,0)<>0 then
f1:=expand(f1/coeff(f1,x,0)):
else
f1:=f:
fi:

for i1 from 1 while expand(coeff(f1,x,i1))=0 do od:
lu:=WhoAreYout(coeff(f1,x,i1),t):

for i from i1 to degree(f1,x) do
if WhoAreYout(coeff(f1,x,i),t)<>lu then
 RETURN(0):
fi:
od:
lu:
end:


#isbad2(evar,xSet): given an expression in the variables
#in the set xSet, decides whether evar has positive degree
#in any of the variables of xSet

isbad2:=proc(evar,xSet)
local i:

for i from 1 to nops(xSet) do
if degree(evar,xSet[i])>0 then
RETURN(true):
fi:
od:
false:
end:

#CT1old(f,xSet,t): Given a rational
#function f of t and a set, xSet, of  x-variables
#finds the constant term in t (i.e. the coeff. of t^0)
#of the Laurent expansion of t, viewed as a formal Laurent
#series in the x's of xSet
#For example try CT1old(1/(1-x*t)/(y-t),[x,y],t);
CT1old:=proc(f,xSet,t)
local f1,i,lu,g,mu,BAD1:
mu:=denom(normal(f)):
mu:={solve(mu,t)}:

BAD1:={}:

for i from 1 to nops(mu) do
if isbad2(mu[i],xSet) or mu[i]=0 then
 BAD1:=BAD1 union {mu[i]}:
fi:
od:

f1:=convert(f,parfrac,t):

g:=f1:
for i from 1 to nops(f1) do
lu:=op(i,f1):
mu:=denom(lu):
if isbad3(mu,BAD1,t) then
g:=g-lu:
fi:

od:


factor(normal(subs(t=0,g))):
end:
 


#isbad3(evar,kv,t): if any of the substitutions of t=kv[i]
#into evar result in 0
isbad3:=proc(evar,kv,t)
local i:

for i from 1 to nops(kv) do

 if expand(subs(t=kv[i],evar))=0 then
  RETURN(true):
 fi:

od:

false:
end:


####################NEW STUFF

isbad:=proc(evar,t,xSet)
local i,j,coe:

for i from 0 to degree(evar,t)-1 do
 coe:=coeff(evar,t,i):
  for j from 1 to nops(xSet) do
    if degree(coe,xSet[j])>0 then
     RETURN(true):
    fi:
  od:
od:
false:
end:

#isbada(evar,BAD1,t): Given a denominator term (coming
#from a single term of a prtial-fraction decomposition w.r.t t)
#decides whether it is of the form BAD1[i]^(some power) for
#some bad term of BAD1
isbada:=proc(evar,BAD1,t)
local j,gu,d1:
if normal(diff(evar,t))=0 then
 RETURN(false):
fi:
for j from 1 to nops(BAD1) do
gu:=BAD1[j]:
d1:=degree(evar,t)/degree(gu,t):
if type(d1,integer) then
  gu:=normal(evar/gu^d1):
  if normal(diff(gu,t))=0 then
   RETURN(true):
 fi:
fi:
od:
false:
end:

#CT1(f,xSet,t): Given a rational
#function f of t and a set, xSet, of  x-variables
#finds the constant term in t (i.e. the coeff. of t^0)
#of the Laurent expansion of t, viewed as a formal Laurent
#series in the x's of xSet
#For example try CT1(1/(1-x*t)/(y-t),[x,y],t);
CT1:=proc(f,xSet,t)
local f1,i,lu,g,mu,BAD1,mu1,mu11,mu11a:
mu:=denom(normal(f)):

BAD1:={}:
if type(mu,`*`) then
for i from 1 to nops(mu) do
mu1:=op(i,mu):
if type(mu1, `^`) then
 mu11:=op(1,mu1):
else
mu11:=mu1:
fi:

mu11a:=expand(mu11/coeff(mu11,t,degree(mu11,t))):

if mu11a=t or isbad(mu11a,t,xSet) then
 BAD1:=BAD1 union {mu11}:
fi:
od:

else

mu1:=mu:
if type(mu1, `^`) then
 mu11:=op(1,mu1):
else
mu11:=mu1:
fi:

mu11a:=expand(mu11/coeff(mu11,t,degree(mu11,t))):

if mu11a=t or isbad(mu11a,t,xSet) then
 BAD1:=BAD1 union {mu11}:
fi:
fi:


f1:=convert(f,parfrac,t):

g:=f1:
for i from 1 to nops(f1) do
lu:=op(i,f1):
mu:=denom(lu):
if isbada(mu,BAD1,t) then
g:=g-lu:
fi:

od:


factor(normal(subs(t=0,g))):
end:
 

#GFx(a,Elist,x): Given a list of variables, a, 
#of type non-negative integer, and a list
#of linear-diophantine equations in these variables
#and a continuous indexed-variable x, outputs the
#generating function of all solutions 
#where x[1]^a[1]*x[2]^a[2]*... shows up iff (a[1],a[2], ...)
#is a solution
GFx:=proc(a,Eqs,x)
local f,i,t,xSet,g:
f:=GFxt(a,Eqs,x,t):
xSet:=[seq(x[i],i=1..nops(a))]:
for i from 1 to nops(Eqs) do
f:=CT1(f,xSet,t[i]):
od:
g:=f:
for i from 1 to nops(a) do
g:=subs(x[i]=x[a[i]],g):
od:
g:
end:


#GFIx(a,Elist,x): Given a list of variables, a, 
#of type non-negative integer, and a list
#of linear-diophantine inequalities in these variables
#and a continuous indexed-variable x, outputs the
#generaitng function of all solutions 
#where x[1]^a[1]*x[2]^a[2]*... shows up iff (a[1],a[2], ...)
#is a solution. Elist is is a list of affine
#linear expressions in the members of a and the inequalities
#are Elist[1]>=0, Elist[2]>=0, ...
GFIx:=proc(a,Eqs,x)
local f,i,t,xSet,g:
f:=GFxt(a,Eqs,x,t):
xSet:=[seq(x[i],i=1..nops(a))]:
for i from 1 to nops(Eqs) do
f:=OMEGA(f,xSet,t[i]):
od:
g:=f:
for i from 1 to nops(a) do
g:=subs(x[i]=x[a[i]],g):
od:
g:
end:


#numfacs1(mekh,x): given a product, mekh, of terms
#counts how many of them depend on x
numfacs1:=proc(mekh,x)
local co,i:
co:=0:
for i from 1 to nops(mekh) do
if degree(op(i,mekh),x)>0 then
 co:=co+1:
fi:
od:
co:
end:


#ToUm1(R,x): Given a rational function R and a variable x
#tries to finds an umbra that sends f(x) to Constant
#term of R(x)f(1/x), it only treats the case where
#the denominator of R(x) factors completely over the 
#rationals and there are no multiplicity and
#the degree of the numerator is less than the degree of
#of the denominator
#if this is not the case, it returns 0
ToUm1:=proc(R,x)
local mekh,gu,i,mu,p:
mekh:=factor(denom(R)):

if degree(numer(R),x)>=degree(denom(R),x) then
 RETURN(0):
fi:

if type(mekh, `*`) and numfacs1(mekh,x)<>degree(mekh,x) then
 RETURN(0):
fi:

gu:={solve(mekh,x)}:

mu:={}:
readlib(residue):
for i from 1 to nops(gu) do
p:=normal(-residue(R,x=gu[i])/gu[i]):
mu:=mu union {[p,[0],[1/gu[i]]]}:
od:
mu:
end:




ToUm:=proc(R,xList)
local gu,mu,i,xList1,mui,mui1,mui2,mui3,x,lu,j:

if nops(xList)=1 then
RETURN(ToUm1(R,xList[1])):
fi:

x:=xList[nops(xList)]:
xList1:=[op(1..nops(xList)-1,xList)]:
mu:=ToUm(R,xList1):

gu:={}:

for i from 1 to nops(mu) do
mui:=mu[i]:
mui1:=mui[1]:
mui2:=mui[2]:
mui3:=mui[3]:

for j from 1 to nops(mui3) do
if diff(mui3[j],x)<>0 then
 RETURN(0):
fi:
od:

lu:=ToUm1(mui1,x):

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

for j from 1 to nops(lu) do
gu:=gu union {[lu[j][1],[op(mui2),0],[op(mui3),op(lu[j][3])]]} :
od:
od:
gu:
end:

#OMEGA(f,xSet,t): Given a rational
#function f of t and a set, xSet, of  x-variables
#Applies MacMahon's operation Omega w.r.t to the
#variable t (i.e. only takes that part of f that involves)
#non-negative powers of t
#Here f is viewed as a formal Laurent in t and a
#formal power series in the x's of xSet
#For example try OMEGA(1/(1-x*t)/(y-t),[x,y],t);
OMEGA:=proc(f,xSet,t)
local f1,i,lu,g,mu,BAD1,mu1,mu11,mu11a:
mu:=denom(normal(f)):

BAD1:={}:
if type(mu,`*`) then
for i from 1 to nops(mu) do
mu1:=op(i,mu):
if type(mu1, `^`) then
 mu11:=op(1,mu1):
else
mu11:=mu1:
fi:

mu11a:=expand(mu11/coeff(mu11,t,degree(mu11,t))):

if mu11a=t or isbad(mu11a,t,xSet) then
 BAD1:=BAD1 union {mu11}:
fi:
od:

else

mu1:=mu:
if type(mu1, `^`) then
 mu11:=op(1,mu1):
else
mu11:=mu1:
fi:

mu11a:=expand(mu11/coeff(mu11,t,degree(mu11,t))):

if mu11a=t or isbad(mu11a,t,xSet) then
 BAD1:=BAD1 union {mu11}:
fi:
fi:


f1:=convert(f,parfrac,t):

g:=f1:
for i from 1 to nops(f1) do
lu:=op(i,f1):
mu:=denom(lu):
if isbada(mu,BAD1,t) then
g:=g-lu:
fi:

od:
g:=subs(t=1,g):
normal(g):
end: