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

#Created: March 2001, previos version: March 2001
#This version:  Dec. 2018
#LinDiophantus.txt: 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.rutgers.edu

with(combinat):

print(`LinDiophantus.txt: 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.rutgers.edu`):

lprint(``):
print(`Created: March 2001`):
print(`This version: Dec. 2018`):
lprint(``):
print(`Written by Doron Zeilberger, DoronZeil at gmail dot com`):
lprint(``):
print(`Please report bugs to  DoronZeil at gmail dot com `):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/`):

print(`----------------------------------------------------`):
 print(`For a list of the Empirical procedures type ezraE(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):

print(`----------------------------------------------------`):
 print(`For a list of the Story procedures type ezraS(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):

print(`----------------------------------------------------`):
 print(`For a list of Supporting procedures type ezra1(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):

print(`----------------------------------------------------`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):


ezra1:=proc()
if nops([args])=0  then
print(` The supporting procedures are: `):
print(`CondorcetIE `):
else
ezra(args):
fi:
end:

ezraS:=proc()
if nops([args])=0  then
print(` The STORY procedures are: `):
print(`BPs, CondorcetS, GFtS , LHPs, SeferBPs`):
else
ezra(args):
fi:
end:

ezraE:=proc()
if nops([args])=0  then
print(` The procedures for the EMPIRICAL approach (written June 2012) are: `):
print(` GFtE, NuMagicSq, NuMagicSqwd, NuP, Pitronot, WtP, RecE `):
else
ezra(args):
fi:
end:

ezra:=proc()
print(``):
if nops([args])=0  then
print(`contains procedures: `):
print(`  BoS, BP,Condorcet, CondorcetG, `):
print(` CondorcetW, CT1, GFIx, GFt, GFtW, GFtoQP`):
print(` GFx, GFxt, LHP, LHPwt, MagicSq, MagicSqwd, MagicSqG, MagicSqGwd,`):
print(` MagicRect, MagicRectG, PPars, Simpson, ToUm , OMEGA`): 



elif nops([args])=1 and op(1,[args])=BoS then
print(`BoS(m): the Boolean lattice of dimension m,`):
print(`try: BoS(3);`):

elif nops([args])=1 and op(1,[args])=BP then
print(`BP(k,r,t): the generating function for partiotions of n`):
print(`p_1>=p_2>=...>=p_k>=0 such that`):
print(`the sum of the r lowest entries is >=the sum of the `):
print(`top k-r entries`):
print(`Try:`):
print(`BP(3,2,t);`):

elif nops([args])=1 and op(1,[args])=BPempir then
print(`BPempir(k,r,t,K): Like BP(k,r,t) but empirically`):
print(`using K terms`):
print(`Try:`):
print(`BPempir(3,2,t,40);`):

elif nops([args])=1 and op(1,[args])=BPs then
print(`BPs(k,r,t): verbose form of BP(k,r,t) (q.v.)`):
print(`Try:`):
print(`BPs(3,2,t);`):

elif nops([args])=1 and op(1,[args])=Condorcet then
print(`Condorcet(t): the generating function in t`):
print(`whose coeff. of t^n (in its Maclaurin expansion)`):
print(`gives you the number of possible "score vectors"`):
print(`for all six total rankings of the three candidates`):
print(`among n votes, that yields the (strict) Condorcet`):
print(`scenatrio. Try:`):
print(`Condorcet(t);`):

elif nops([args])=1 and op(1,[args])=CondorcetG then
print(`CondorcetG(t,c12,c23,c31): the generating function in t`):
print(`whose coeff. of t^n (in its Maclaurin expansion)`):
print(`gives you the number of possible "score vectors"`):
print(`for all six total rankings of the three candidates`):
print(`among n votes, that yields the (strict) Condorcet`):
print(`scenatrio with margins `):
print(` 1 beats 2 with a margin of c12`):
print(` 2 beats 3 with a margin of c23`): :
print(` 3 beats 1 with a margin of c31 `):
print(` Try:`):
print(`CondorcetG(t,2,2,2);`):

elif nops([args])=1 and op(1,[args])=CondorcetIE then
print(`CondorcetIE(p123,p132,p213,p231,p312,p321): `):
print(`The Condorcet inequalities`):
print(`Try:`):
print(`CondorcetIE(p123,p132,p213,p231,p312,p321);`):

elif nops([args])=1 and op(1,[args])=CondorcetS then
print(`CondorcetS(t): verbose form of Condorcet(t) (q.v.)`):
print(`Try:`):
print(`CondorcetS(t);`):

elif nops([args])=1 and op(1,[args])=CondorcetW then
print(`CondorcetW(t): the generating function in t`):
print(`whose coeff. of t^n (in its Maclaurin expansion)`):
print(`gives you the number of possible "score vectors"`):
print(`for all six total rankings of the three candidates`):
print(`among n votes, that yields the (weak) Condorcet`):
print(`scenatrio. Try:`):
print(`CondorcetW(t);`):

elif 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 CT1(1/(1-x*t)/(y-t),[x,y],t);`):


elif nops([args])=1 and op(1,[args])=GFt then
print(`GFt(a,Ineqs,t): Inputs a list of discrete variables`):
print(`a and a list of affine-linear expressions, Ineqs,`):
print(`finds the rational generating function, in t,`):
print(`whose coefficient of t^n (in the Maclaurin expansion)`):
print(`is the number of vectors of length nops(a) that`):
print(`add up to n and obey the inequalities, For example try:`):
print(`GFt([a,b,c],[a-b,b-c],t);`):

elif nops([args])=1 and op(1,[args])=GFtS then
print(`GFtS(a,Ineqs,t): verbose form of GFt(a,Ineqs,t) (q.v.)`):
print(`For example try:`):
print(`GFtS([a,b,c],[a-b,b-c],t);`):

elif nops([args])=1 and op(1,[args])=GFtW then
print(`GFtW(a,Ineqs,x,t): `):
print(`Inputs a list of discrete variables`):
print(`a and a list of affine-linear expressions, Ineqs,`):
print(`and variable x (for indexed variable) and variable t`):
print(`finds the rational generating function, in t,`):
print(`and x[a[1]], ..., x[a[nops(a)]],`):
print(`whose coefficient of t^n (in the Maclaurin expansion)`):
print(`is the weight-enumerator of vectors of length nops(a) that`):
print(`add up to n and obey the inequalities, For example try:`):
print(`GFtW([a,b,c],[a-b,b-c],x,t);`):

elif nops([args])=1 and op(1,[args])=GFtoQP then
print(`GFtoQP(f,q,K,n): converts a rational function`):
print(`whose denominators is a product of (1-t^i)'s`):
print(`into a sum of quasi-polynomials`):
print(`For example, try:`):
print(`GFtoQP(1/(1-q)/(1-q^2),q,20,n);`):

elif nops([args])=1 and op(1,[args])=GFtE then
print(`GFtE(a,Ineqs,t,K): an empirical approach to`):
print(`GFt(a,Ineqs,t) (q.v.), where K is the max. number of data`):
print(`points to gather. If they do not suffice, it returns FAIL.`):
print(`For example try:`):
print(`GFtE([a,b,c],[a-b,b-c],t,40);`):


elif 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); `):


elif 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);`):


elif 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);`):

elif nops([args])=1 and op(1,[args])=LHP then
print(`LHP(j,y,t): The weight-enumerator of all Lecture Hall Partitions`):
print(`of length j, b[j]+...+b[1] (Lecture Hall partitions, defined`):
print(`by  Mireille Bousquet-Melou and Kimmo Eriksson means that`):
print(` b[j]/(j-1)>=b[j-1]/(j-1)>= ... >=b[1] >=0)`):
print(`according to the weight `):
print(` t^(sum of entries)*y^(b[j]-b[j-1]+(-1)^(j-1)*b[1]) `):
print(`Try: `):
print(`LHP(j,y,t); `):

elif nops([args])=1 and op(1,[args])=LHPs then
print(`LHPs(J): tells the story of Lecture Hall Partitions up to J parts`):

elif nops([args])=1 and op(1,[args])=LHPwt then
print(`LHPwt(j,x,t): The weight-enumerator of all Lecture Hall Partitions`):
print(`of length j, b[j]+...+b[1] (Lecture Hall partitions, defined`):
print(`by  Mireille Bousquet-Melou and Kimmo Eriksson means that`):
print(`b[j]/(j-1)>=b[j-1]/(j-1)>= ... >=b[1] >=0)`):
print(`according to the weight `):
print(`t^(sum of entries)*x[1]^b[1]*...*x[j]^b[j] `):
print(`Try:`):
print(`LHPwt(j,x,t);`):

elif nops([args])=1 and op(1,[args])=MagicSq then
print(`MagicSq(r,z): The weight-enumerator of ALL r by r`):
print(`arrays of non-negative integers all whose line-sums are`):
print(`the same (by line I mean row or column)`):
print(`according to the weight z^LineSum`):
print(`Try:MagicSq(2,z);`):

elif nops([args])=1 and op(1,[args])=MagicSqwd then
print(`MagicSqwd(r,z): The weight-enumerator of ALL r by r`):
print(`arrays of non-negative integers all whose line-sums are`):
print(`the same (by line I mean row or column)`):
print(`AND also  the two main diagonals`):
print(`according to the weight z^LineSum`):
print(`Try:MagicSqwd(2,z);`):

elif nops([args])=1 and op(1,[args])=MagicSqG then
print(`MagicSqG(r,x,z): The weight-enumerator of ALL r by r`):
print(`arrays of non-negative integers all whose line-sums are`):
print(`the same (by line I mean row or column)`):
print(`according to the weight x[i,j]^M[i,j]*z^LineSum`):
print(`Try:MagicSqG(2,x,z);`):

elif nops([args])=1 and op(1,[args])=MagicSqGwd then
print(`MagicSqGwd(r,x,z): The weight-enumerator of ALL r by r`):
print(`arrays of non-negative integers all whose line-sums are`):
print(`the same (by line I mean row or column)`):
print(`AND also the two MAIN DIAGONALS add up to the same`):
print(`according to the weight x[i,j]^M[i,j]*z^LineSum`):
print(`Try:MagicSqGwd(2,x,z);`):

elif nops([args])=1 and op(1,[args])=MagicRect then
print(`MagicRect(r,s,z,w): The weight-enumerator of ALL r by s`):
print(`arrays of non-negative integers all whose row sums are`):
print(`the same and all whose column-sums are the same`):
print(`according to the weight z^RowSum*w^ColumnSum`):
print(`Try:MagicRect(2,2,z,w);`):

elif nops([args])=1 and op(1,[args])=MagicRectG then
print(`MagicRectG(r,s,x,z,w): The weight-enumerator of ALL r by s`):
print(`arrays of non-negative integers all whose row sums are`):
print(`the same and all whose column-sums are the same`):
print(`according to the weight x[i,j]^M[i,j]*z^RowSum*w^ColumnSum`):
print(`Try:MagicRectG(2,2,x,z,w);`):


elif nops([args])=1 and op(1,[args])=NuMagicSq then
print(`NuMagicSq(r,N): The number of r by r semi-magic squares`):
print(`with line-sum N`):
print(`Try: NuMagicSqG(2,4);`):

elif nops([args])=1 and op(1,[args])=NuMagicSqwd then
print(`NuMagicSqwd(r,N): The number of r by r semi-magic squares`):
print(`with line-sum N and the two MAIN DIAGONALS`):
print(`Try: NuMagicSqGwd(2,4);`):

elif nops([args])=1 and op(1,[args])=NuP then
print(`NuP(a,Ineqs,n): `):
print(`INPUTS: a list of discrete variables`):
print(`a and a list of affine-linear expressions, Ineqs,`):
print(`and a positive integer n.`): 
print(`OUTPUTS: the NUMBER`):
print(`of vectors of length nops(a), with non-negative integer`):
print(`entries that adds up to n and makes all the`):
print(`members of Ineqs non-negative.`):
print(`For example try:`):
print(`NuP([a,b,c],[a-b,b-c],10);`):

elif nops([args])=1 and op(1,[args])=Pitronot then
print(`Pitronot(a,Ineqs,n): `):
print(`INPUTS: a list of discrete variables`):
print(`a and a list of affine-linear expressions, Ineqs,`):
print(`and a positive integer n.`):
print(` OUTPUS: the SET`):
print(`of vectors of length nops(a), with non-negative integer`):
print(`entries that adds up to n and makes all the`):
print(`members of Ineqs non-negative.`):
print(`For example try:`):
print(`Pitronot([a,b,c],[a-b,b-c],10);`):


elif nops([args])=1 and op(1,[args])=PPars then
print(`PPars(m,P,t): Inputs a pos. integer n, and a poset P`):
print(`with vertices {1, ..., m}`):
print(`given a set of edges of its Hasse diagram,`):
print(`outputs the generating function for the number of`):
print(`P-partitions of n. Try:`):
print(`PPars(4,{[1,2],[1,3],[2,4],[3,4]},t);`):

elif 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);`):

elif nops([args])=1 and op(1,[args])=RecE then
print(`RecE(a,Ineqs,n,N,K): inputs a list of`):
print(`discrete variables a, a list if affine-linear`):
print(`expressions in them and tries to find`):
print(`a recurrence operator ope(n,N) with ORD+DEG<=K, annihilating`):
print(`the sequence n!*sum(Wt(a), |a|=n, a satisfies the ineqalities`):
print(`of Ineqs), using K data points. If it fails,`):
print(`it returns the terms of the sequence that it found`):
print(`For example try:`):
print(`RecE([a,b],[a-b],n,N,50);`):

elif nops([args])=1 and op(1,[args])=SeferBPs then
print(`SeferBPs(K,t): inputs a positive integer and`):
print(`outputs all the theorems about partitions where`):
print(`the rich are not much richer than the poor.`):
print(`Try SeferBPs(5,t);`):

elif nops([args])=1 and op(1,[args])=Simpson then
print(`Simpson(a,b,t): Inputs positve inetegers, a and b, and a variable t.`):
print(`Outputs the generating function, in t,  for Simpson's paradox where`):
print(` t^n is the number of 2 by 2 matrices of non-negative integers, [[x11,x12],[x21,x22]] such that `):
print(`x11<=a*n, x12<=b*n, x21<=b*n, x22<=a*n, x11/a<x12/b, x21/b<x22/a, yet x11+x21-x12-x22>0`):
print(`Try: `):
print(`Simpson(2,1,t); `):

elif 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`):

elif nops([args])=1 and op(1,[args])=WtP then
print(`WtP(a,Ineqs,n): `):
print(`INPUTS: a list of discrete variables`):
print(`a and a list of affine-linear expressions, Ineqs,`):
print(`and a positive integer n.`): 
print(`OUTPUTS: the SUM of 1/(a1!*a2!...) `):
print(`of all vectors [a1,a2, ...]`):
print(`of length nops(a), with non-negative integer`):
print(`entries that adds up to n and makes all the`):
print(`members of Ineqs non-negative.`):
print(`For example try:`):
print(`WtP([a,b,c],[a-b,b-c],10);`):

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



############Start guessing procedures
###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L,ope1:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
   RETURN(ope):

      fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
###end Findrec

#from Cfinite
#SeqFromRec(C,N): the first N terms of the C-finite
#sequence defined by C, try:
#SeqFromRec([[-1,-1],[1,1]],10);
SeqFromRec:=proc(C,N) local gu,C1,C2,i,i1:
C1:=C[1]: C2:=C[2]:
gu:=C2:

for i from 1 to N-nops(C2) do
gu:=[op(gu),-add(C1[i1]*gu[nops(gu)+1-i1],i1=1..nops(C1))]:
od:

gu:

end:

#GFfromRec(C,t): the generating function of the C-finite
#sequence C (starting with n=0)
#Try:
#GFfromRec([[-1,-1],[1,1]],t);
GFfromRec:=proc(C,t)  local mu,lu,N,lu1,i,i1:
N:=5*nops(C[1])+5*nops(C[2]):
mu:=SeqFromRec(C,N):
lu:=add(mu[i]*t^(i-1),i=1..nops(mu))*(1+add(C[1][i1]*t^i1,i1=1..nops(C[1]))):
lu:=expand(lu):
lu:=add(coeff(lu,t,i)*t^i,i=0..N-nops(C[2])):

if degree(lu,t)>=2*nops(C[1]) then
print(lu):
RETURN(FAIL):
fi:
lu:=lu/(1+add(C[1][i1]*t^i1,i1=1..nops(C[1]))):
N:=2*(nops(C[1])+5*nops(C[2])):
lu1:=taylor(lu,t=0,N+1):
mu:=SeqFromRec(C,N):

if [seq(coeff(lu1,t,i),i=0..N-1)]<>mu then
 RETURN(FAIL):
fi:

factor(lu):

end:


#GuessRec1(L,r): inputs a sequence of numbers, L, and a pos.
#integer r, and tries to guess a linear recurrence equation
#of order r satisfied by it
#L(n)+c1*L(n-1)+...+cr*L(n-r)=0 for all n within the range
#The  output would be given  as a list of r numbers
#[[c1,c2,...,cr],InitialConditions]
#For example, the Fibonacci sequence is
#[[-1,-1], [1,1]] 
GuessRec1:=proc(L,r) local c,i,var,eq,n:

if nops(L)<=2*r+4 then
 RETURN(FAIL):
fi:


var:={seq(c[i],i=1..r)}:


eq:={seq( L[n]+add(c[i]*L[n-i],i=1..r)=0 , n=r+1..nops(L))}:

var:=solve(eq,var):

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

[[seq(subs(var, c[i]),i=1..r)], L[1..r]   ] ;

end:


#GuessRec(L): inputs a sequence of numbers, L, 
# and tries to guess a linear recurrence equation
#of some order, r, satisfied by it
#L(n)+c1*L(n-1)+...+cr*L(n-r)=0 for all n within the range
#The  output would be given  as a list of r numbers
#[[c1,c2,...,cr],InitialConditions]
#For example, the Fibonacci sequence is
#[[-1,-1], [1,1]] 
GuessRec:=proc(L) local r,ans:

for r from 1 to (nops(L)-4)/2 do
 ans:=GuessRec1(L,r):
 if ans<>FAIL then
  RETURN(ans):
 fi:
od:

FAIL:
end:

#End Cfinite
############End guessing procedures
#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:
 




#GFt(a,Ineqs,t): 
#Inputs a list of discrete variables
#a and a list of affine-linear expressions, Ineqs,
#finds the rational generating function, in t,
#whose coefficient of t^n (in the Maclaurin expansion)
#is the number of vectors of length nops(a) that
#add up to n and obey the inequalities, For example try:
#GFt([a,b,c],[a-b,b-c],t);
GFt:=proc(a,Ineqs,t) local var,Eqs,i,x,y,gu:
var:=a:
Eqs:=[]:

for i from 1 to nops(Ineqs) do
 var:=[op(var),y[i]]:
 Eqs:=[op(Eqs),Ineqs[i]-y[i]]:
od:

gu:=GFx(var,Eqs,x);

gu:=subs({seq(x[a[i]]=t,i=1..nops(a)),seq(x[y[i]]=1,i=1..nops(Ineqs))},gu):
factor(gu):
end:



#Pitornot(a,Ineqs,n): 
#Inputs a list of discrete variables
#a and a list of affine-linear expressions, Ineqs,
#and a positive integer n, outputs the set
#of vectors of length nops(a), with non-negative integer
#entries that adds up to n and makes all the
#members of Ineqs non-negative.
#For example try:
#Pitronot([a,b,c],[a-b,b-c],10);
Pitronot:=proc(a,Ineqs,n) local gu,lu,i,k,ak,gu1,var,Ineqs1,gu11:
option remember:

if nops(a)=1 then

lu:=subs(a[1]=n,Ineqs):

 if {seq(evalb(lu[i]>=0),i=1..nops(lu))}={true} then
   RETURN({[n]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:


k:=nops(a):


for ak from 0 to n do
var:=[op(1..k-1,a)]:
Ineqs1:=subs(a[k]=ak,Ineqs):
gu1:=Pitronot(var,Ineqs1,n-ak):
gu:=gu union {seq([op(gu11),ak],gu11 in gu1)}:
od:

gu:

end:




#NuP(a,Ineqs,n): 
#Inputs a list of discrete variables
#a and a list of affine-linear expressions, Ineqs,
#and a positive integer n, outputs the NUMBER 
#of vectors of length nops(a), with non-negative integer
#entries that adds up to n and makes all the
#members of Ineqs non-negative.
#For example try:
#NuP([a,b,c],[a-b,b-c],10);
NuP:=proc(a,Ineqs,n) local gu,lu,i,k,ak,var,Ineqs1:
option remember:

if nops(a)=1 then

lu:=subs(a[1]=n,Ineqs):

 if {seq(evalb(lu[i]>=0),i=1..nops(lu))}={true} then
   RETURN(1):
 else
  RETURN(0):
 fi:
fi:

gu:=0:


k:=nops(a):


for ak from 0 to n do
 var:=[op(1..k-1,a)]:
 Ineqs1:=subs(a[k]=ak,Ineqs):
 gu:=gu+NuP(var,Ineqs1,n-ak):
od:

gu:

end:

  




#GFtE(a,Ineqs,t,K): an empirical approach to
#GFt(a,Ineqs,t) (q.v.), K is the largest
#term of the sequence to go to. If none found, it
#returns fail
#For example try:
#GFtE([a,b,c],[a-b,b-c],t,50);
GFtE:=proc(a,Ineqs,t,K) local gu,n1,lu,i:


for n1 from 10 to K by 10 do
gu:=[seq(NuP(a,Ineqs,i,K),i=0..n1)]:
lu:=GuessRec(gu):
if lu<>FAIL then
 RETURN(GFfromRec(lu,t)):
fi:
od:
print(K, `is not big enough, make last argument larger`):
FAIL:
end:



#WtP(a,Ineqs,n): 
#Inputs a list of discrete variables
#a and a list of affine-linear expressions, Ineqs,
#and a positive integer n, outputs the sum of the weights 
#of vectors of length nops(a), with non-negative integer
#entries that adds up to n and makes all the
#members of Ineqs non-negative,
#where the weight is 1/(a[1]!*...*a[nops(a)]!)
#For example try:
#WtP([a,b,c],[a-b,b-c],10);
WtP:=proc(a,Ineqs,n) local gu,lu,i,k,ak,var,Ineqs1:
option remember:

if nops(a)=1 then

lu:=subs(a[1]=n,Ineqs):

 if {seq(evalb(lu[i]>=0),i=1..nops(lu))}={true} then
   RETURN(1/n!):
 else
  RETURN(0):
 fi:
fi:

gu:=0:


k:=nops(a):


for ak from 0 to n do
 var:=[op(1..k-1,a)]:
 Ineqs1:=subs(a[k]=ak,Ineqs):
 gu:=gu+WtP(var,Ineqs1,n-ak)/ak!:
od:

gu:

end:



#RecE(a,Ineqs,n,N,K): inputs a list of
#discrete variables a, a list if affine-linear
#expressions in them and tries to find
#a recurrence operator ope(n,N) with ORD+DEG<=K, annihilating
#the sequence n!*sum(Wt(a), |a|=n, a satisfies the ineqalities
#of Ineqs), using K data points. If it fails,
#it returns the terms of the sequence that it found
#it just returns the first K value
#GFt(a,Ineqs,t) (q.v.), K is the largest
#term of the sequence to go to. If none found, it
#returns fail
#For example try:
#RecE([a,b],[a-b],n,N,50);
RecE:=proc(a,Ineqs,n,N,MaxC) local gu,ope,L,i,DEG,ORD,ope1:


for L from 4 to MaxC do

gu:=[seq(WtP(a,Ineqs,i),i=0..trunc((L/2+3/2)^2+5))]:
gu:=[seq(i!*gu[i],i=1..nops(gu))]:
ope:=Findrec(gu,n,N,L):

if ope<>FAIL then
DEG:=degree(numer(ope),n)+2:
ORD:=degree(numer(ope),N)-1:

gu:=[seq(WtP(a,Ineqs,i),i=0..(1+DEG)*(1+ORD)+5+ORD)]:
gu:=[seq(i!*gu[i],i=1..nops(gu))]:

ope1:=findrec(gu,DEG,ORD,n,N):

 if ope1<>FAIL then
 RETURN(ope1,[op(1..degree(ope1,N),gu)]):
  else
  RETURN(ope,[op(1..degree(ope,N),gu)]):
 fi:
fi:
od:

gu:

end:



######Begin GFtoQP#########

#GuessPol1(L,d,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i+1] for i=1..nops(L)
#For example, try: 
#GuessPol1([seq(i,i=1..10)],1,n);
GuessPol1:=proc(L,d,n) local P,i,a,eq,var:
if d>nops(L)-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=1..nops(L))}:

var:=solve(eq,var):

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

subs(var,P):

end:

#GuessQP11(L,d,r,n): Inputs a list L and outputs a list of polynomials
#of degree d
#of length r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r, For example try:
#GuessQP11:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],0,2,n);
GuessQP11:=proc(L,d,r,n) local M,i,L1,gu,m:
M:=[]:
for i from 1 to r do
L1:=[seq(L[i+(m-1)*r],m=1..trunc((nops(L)-i)/r)+1)]:
gu:=GuessPol1(L1,d,m):
if gu=FAIL then
  RETURN(FAIL):
fi:
gu:=expand(subs(m=(n-i)/r+1,gu)):
M:=[op(M),gu]:
end:
M:
end:


#GuessQP2(L,d,n): Inputs a list L and outputs a list of polynomials
#of degree d, the list being
#length r, for the smallest possible
#r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r, 
#For example try:
#GuessQP2:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],0,n);
GuessQP2:=proc(L,d,n) local r,gu:

for r from 1 to trunc(nops(L)/(d+2)) do
gu:=GuessQP11(L,d,r,n):
if gu<>FAIL then
 RETURN(gu):
fi:
od:
FAIL:
end:

#FindExp(P): Given an expression of the form
#P=C*pol(q)^i finds i, for example, try:
#FindExp(6*(1+q+q^2)^3);
FindExp:=proc(P) local P1:

if type(P,`^`) then
 RETURN(op(2,P)):
fi:

if type(P,`*`) then
  P1:=op(2,P):
  if type(P1,`^`) then
    RETURN(op(2,P1)):
  else
    RETURN(1):
  fi:
fi:
1:
end:

#KhaberLists(S): given a set of lists of the same size, adds them
#up, returns FAIL otherwise. For example, type:
#KhaberLists({[1,2],[2,4]});
KhaberLists:=proc(S) local s,kv,k,i1:

kv:={seq(nops(s), s in S)}:

if nops(kv)<>1 then
 RETURN(FAIL):
fi:
k:=kv[1]:

[seq(sort(add(s[i1],s in S)),i1=1..k)]:

end:

#EvalQP(Q,n,n0): evaluates the quasi-polynomial Q of n at n0
#For example try:
#EvalQP([2*n^3,n^3],n,5);
EvalQP:=proc(Q,n,n0) local i:

i:=n0 mod nops(Q):

if i=0 then i:=nops(Q): fi:

subs(n=n0,Q[i]):
end:



#EvalQPS(L,n,n0): evaluates the sum of the 
#quasi-polynomials in the list L of n at n0
#For example try:
#EvalQPS([[n^2],[2*n^3,n^3],[1,2,6]],n,5);
EvalQPS:=proc(L,n,n0) local i:
add(EvalQP(L[i],n,n0),i=1..nops(L)):
end:

#GFtoQP(f,q,K,n): converts a rational function
#whose denominators is a product of (1-t^i)'s
#into a sum of quasi-polynomials
#For example, try:
#GFtoQP(1/(1-q)/(1-q^2),q ,20,n
GFtoQP:=proc(f,q,K,n) local gu,i,lu1,d,H,i1,gu1,T,makh,mu:
option remember:
gu:=convert(f,parfrac):


H:={}:

for i from 1 to nops(gu) do
 gu1:=op(i,gu):



 d:=FindExp(denom(gu1))-1:

 lu1:=taylor(gu1,q=0,K+1):
 lu1:=[seq(coeff(lu1,q,i1),i1=1..K)]:

 lu1:= GuessQP2(lu1,d,n):

   if lu1=FAIL then
     RETURN(FAIL):
   else
    H:=H union {lu1}:
   fi:
od:

makh:={seq(nops(H[i1]),i1=1..nops(H))}:
makh:=sort(convert(makh,list)):

for i1 from 1 to nops(makh) do
 T[makh[i1]]:= {}:
od:

for i1 from 1 to nops(H) do
T[nops(H[i1])]:=T[nops(H[i1])] union {H[i1]}:
od:

H:=[seq(T[makh[i1]],i1=1..nops(makh))]:

H:=[seq(KhaberLists(H[i1]),i1=1..nops(H))]:

gu:=taylor(f,q=0,2*K+9):
gu:=[seq(coeff(gu,q,i),i=1..2*K+8)]:
mu:=[seq(EvalQPS(H,n,i),i=1..2*K+8)]:
if gu<>mu then
print(H, `didn't work out, too bad`):
RETURN(FAIL):
fi:
H:
end:

####End GFtoQP#########



#CondorcetIE(p123,p132,p213,p231,p312,p321): 
#The Condorcet inequalities
#Try:
#CondorcetIE(p123,p132,p213,p231,p312,p321);
CondorcetIE:=proc(p123,p132,p213,p231,p312,p321)
[p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321-1,
 p123+p213+p231-p132-p312-p321-1,
p312+p321+p231-p132-p123-p213-1]:
end:

#Condorcet(t): the generating function in t
#whose coeff. of t^n (in its Maclaurin expansion)
#gives you the number of possible "score vectors"
#for all six total rankings of the three candidates
#among n votes, that yields the (strict) Condorcet
#scenatrio. Try:
#Condorcet(t);
Condorcet:=proc(t) local p123,p132,p213,p231,p312,p321:
GFt([p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321-1,
 p123+p213+p231-p132-p312-p321-1,
p312+p321+p231-p132-p123-p213-1],t):
end:

#CondorcetE(t): Like Condorcet(t) (q.v.) but with the
#empirical approach, try:
#CondorcetE(t);
CondorcetE:=proc(t) local p123,p132,p213,p231,p312,p321:
GFtE([p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321-1,
 p123+p213+p231-p132-p312-p321-1,
p312+p321+p231-p132-p123-p213-1],t,30):
end:

#CondorcetW(t): the generating function in t
#whose coeff. of t^n (in its Maclaurin expansion)
#gives you the number of possible "score vectors"
#for all six total rankings of the three candidates
#among n votes, that yields the (weak) Condorcet
#scenatrio. Try:
#CondorcetW(t);
CondorcetW:=proc(t) local p123,p132,p213,p231,p312,p321:
GFt([p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321,
 p123+p213+p231-p132-p312-p321,
p312+p321+p231-p132-p123-p213],t):
end:

#CondorcetWE(t): Like CondorcetW(t) (q.v.) but with the
#empirical approach, try:
#CondorcetWE(t);
CondorcetWE:=proc(t) local p123,p132,p213,p231,p312,p321:
GFtE([p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321,
 p123+p213+p231-p132-p312-p321,
p312+p321+p231-p132-p123-p213],t,30):
end:



#BP(k,r,t): the generating function for partiotions of n
#p_1>=p_2>=...>=p_k>=0 such that
#the sum of the r lowest entries is >=the sum of the 
#top k-r entries
BP:=proc(k,r,t) local i,p:
GFt([seq(p[i],i=1..k)],
[seq(p[i]-p[i+1],i=1..k-1),p[k]-1,
add(p[i],i=k-r+1..k)-add(p[i],i=1..k-r)-1],t);
end:

#BPempir(k,r,t,K): the generating function for partiotions of n
#p_1>=p_2>=...>=p_k>=0 such that using guessing
#the sum of the r lowest entries is >=the sum of the 
#top k-r entries.
#Try:
#BPempir(3,2,t,40);
BPempir:=proc(k,r,t,K) local i,p:
GFtE([seq(p[i],i=1..k)],
[seq(p[i]-p[i+1],i=1..k-1),p[k]-1,
add(p[i],i=k-r+1..k)-add(p[i],i=1..k-r)-1],t,K);
end:



####Start story procedures

#GFtS(a,Ineqs,t): verbose form of GFtS(a,Ineqs,t)
GFtS:=proc(a,Ineqs,t) local lu,f,n,i,t0,mu,K,ku1,ku2:

t0:=time():
print(`On the Number of Non-Negative Integer Sequencs that Add-up to n`):
print(`and that obey a certain set of inequalities`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

lu:=GFt(a,Ineqs,t):

print(`Theorem: Let `, f(n), ` be the number of vectors of non-negative integers`):
print(a):
print(`that add-up to`,n, `and in addition satisfy the following`):
print(` inequalities: `):

for i from 1 to nops(Ineqs) do
print(Ineqs[i]>=0):
od:

print(`Then we have the generating function `):
print(Sum(f(n)*t^n,n=0..infinity)=lu):
print(``):
print(`In Maple input notation  the right side is:`):
lprint(lu):

K:=2*degree(denom(lu),t)+6:
mu:=GFtoQP(lu,t,K,n):

if mu<>FAIL then
print(f(n), `itself can be given explicitly as a sum of the following`):
print(`quasi-polynomials, where a list of polynomials of length k`):
print(`indicates a quasi-polynomial whose value when n is i (mod k)`):
print(`is given by the i-th entry, for i=1..k .`):

print(f(n)=mu):
print(`In Maple input notation  the right side is:`):
lprint(mu):
ku1:=[seq(coeff(taylor(lu,t=0,41),t,i),i=1..40)]:
ku2:=[seq(EvalQPS(mu,n,i),i=1..40)]:
if ku1<>ku2 then
lprint(lu,mu):
print(ku1,ku2):
 print(`Something bad happended`):
else
 print(`For the sake of Sloane, the first 40 terms (starting with n=1) are:`):
print(ku1):
fi:

fi:
print(`-------------------------------------------------------`):
print(`This took`, time()-t0, `seconds. `):

end:





#CondorcetS(t): 
#verbose form of Condorcet(t) (q.v.)
#CondorcetS(t);
CondorcetS:=proc(t) local p123,p132,p213,p231,p312,p321,lu,
ku1,ku2,i,K,t0,mu,n,f,luW,muE,muO:

t0:=time():
print(`The Exact Number of Ways that the Condorcet Paradox`):
print(`Can Happen `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

lu:=GFt([p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321-1,
 p123+p213+p231-p132-p312-p321-1,
p312+p321+p231-p132-p123-p213-1],t):

luW:=GFtW([p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321-1,
 p123+p213+p231-p132-p312-p321-1,
p312+p321+p231-p132-p123-p213-1],x,t):


print(`Theorem:, Suppose that you have n voters and three candidates`):
print(`and each voter is asked to totally-order the three candidates`):
print(`and then the votes given to each of the six rankings is counted`):
print(`resulting in a certain composition of n into six parts.`):
print(`Of course the total number of such compositions is`):
print(factor(binomial(n+5,6))):
print(`Let `, f(n), `be the number of such compositions that give`):
print(`rise to the famous Condorcet paradox`):
print(`i.e. a majority of the voters prefer`):
print(` candidate 1 to candidate 2,`):
print(`candidate 2 to candidate 3,`):
print(`candidate 3 to candidate 1,`):
print(`or vice-versa `):
print(`Then we have the generating function `):
print(Sum(f(n)*t^n,n=0..infinity)=lu):
print(``):
print(`In Maple input notation  the right side is:`):
lprint(lu):

print(`More informatively, the weighted generating function`):
print(`whose coeff. of t^n gives all the score vectors is: `):
print(luW):
print(`and in Maple input notation`):
lprint(luW):

K:=2*degree(denom(lu),t)+6:
mu:=GFtoQP(lu,t,K,n):
muE:=factor(subs(n=2*n,mu[1][1]+mu[2][2])):
muO:=factor(subs(n=2*n+1,mu[1][1]+mu[2][1])):

mu:=factor(mu[1][1]+1/2*mu[2][1]+1/2*mu[2][2])+
factor((-mu[2][1]+mu[2][2])/2)*(-1)^n:
print(`In fact, we have the following explicit expression`):
print(f(n)=mu):
print(``):
print(`Alternatively:`):
print(f(2*n)=muE):
print(``):
print(f(2*n+1)=muO):
print(`In Maple input notation  the right side is:`):
lprint(mu):

ku1:=[seq(coeff(taylor(lu,t=0,41),t,i),i=0..40)]:
ku2:=[seq(subs(n=i,mu),i=0..40)]:
if ku1<>ku2 then
 print(`Something terrible happened`):
 RETURN(FAIL):
else
 print(`For the sake of Sloane, the first 40 terms are:`):
print(ku1):
fi:
print(`-------------------------------------------------------`):
print(`This took`, time()-t0, `seconds. `):


end:




#BPs(k,r,t): verbose form of BP(k,r,t)
BPs:=proc(k,r,t) local i,lu,n,K,ku1,ku2,mu,t0:

t0:=time():
print(`The Number of Partitions of n into`, k, `non-negative parts`):
print(`where the sum of the`, r,  `lowest entries is At Least the sum of the `):
print(`top ` ,  k-r, ` entries `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

lu:=BP(k,r,t):

print(`Theorem: Let `, f(n), ` be the number of non-increasing`):
print(`arrays of non-negative integers of length`, k ):
print(`such that the sum of the`, r ,`lowest entries is at Least the sum of the `):
print(`top ` ,  k-r, ` entries. `):
print(`We have the generating function `):
print(Sum(f(n)*t^n,n=0..infinity)=lu):
print(``):
print(`In Maple input notation  the right side is:`):
lprint(lu):

K:=2*degree(denom(lu),t)+6:
mu:=GFtoQP(lu,t,K,n):

if mu<>FAIL then
print(f(n), `itself can be given explicitly as a sum of the following`):
print(`quasi-polynomials, where a list of polynomials of length k`):
print(`indicates a quasi-polynomial whose value when n is i (mod k)`):
print(`is given by the i-th entry, for i=1..k .`):

print(f(n)=mu):
print(`In Maple input notation  the right side is:`):
lprint(mu):
ku1:=[seq(coeff(taylor(lu,t=0,41),t,i),i=1..40)]:
ku2:=[seq(EvalQPS(mu,n,i),i=1..40)]:
if ku1<>ku2 then
lprint(lu,mu):
print(ku1,ku2):
 print(`Something bad happended`):
else
 print(`For the sake of Sloane, the first 40 terms (starting with n=1) are:`):
print(ku1):
fi:

fi:
print(`-------------------------------------------------------`):
print(`This took`, time()-t0, `seconds. `):

end:



#PPars(m,P,t): Inputs a pos. integer n, and a poset P
#with vertices {1, ..., m}
#given a set of edges of its Hasse diagram,
#outputs the generating function for the number of
#P-partitions of n. Try:
#PPars(4,{[1,2],[1,3],[2,4],[3,4]},t);
PPars:=proc(m,P,t) local i,j, Ine,x,e,a:

if P minus {seq(seq([i,j],j=i+1..m),i=1..m)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:


a:=[seq(x[i],i=1..m)]:
Ine:=[seq(x[e[2]]-x[e[1]],e in P)]:
GFt(a,Ine,t):
end:



#BoS(m): the Boolean lattice of dimension m,
#try: BoS(3);
BoS:=proc(m) local gu1,gu2,co,T,E,i,j,k1,k2:
co:=0:

E:={seq([1,i],i=2..m+1)}:

T[[]]:=1:
for i from 1 to m do
 T[[i]]:=i+1:
od:

gu1:=choose(m,1):
co:=m+1:

for i from 2 to m do
 gu2:=choose(m,i):

for j from 1 to nops(gu2) do
 co:=co+1:
 T[gu2[j]]:=co:
od:


for k1 from 1 to nops(gu1) do
 for k2 from 1 to nops(gu2) do
  if convert(gu1[k1],set) minus convert(gu2[k2],set)={} then
    E:=E union {[T[gu1[k1]],T[gu2[k2]]]}:
  fi:
 od:
od:
gu1:=gu2:
od:
co,E:
end:



#BPs1(k,r,t,co): verbose form of BP(k,r,t) with theorem number
BPs1:=proc(k,r,t,co) local i,lu,n,K,ku1,ku2,mu:
lu:=BP(k,r,t):
print(`Theorem Number`, co, `: Let `, f(n), ` be the number of non-increasing`):
print(`arrays of non-negative integers of length`, k ):
print(`such that the sum of the`, r , `lowest entries is at Least the sum of the `):
print(`top ` ,  k-r, ` entries. `):
print(`We have the generating function `):
print(Sum(f(n)*t^n,n=0..infinity)=lu):
print(``):
print(`In Maple input notation  the right side is:`):
lprint(lu):

K:=2*degree(denom(lu),t)+6:
mu:=GFtoQP(lu,t,K,n):

if mu<>FAIL then
print(f(n), `itself can be given explicitly as a sum of the following`):
print(`quasi-polynomials, where a list of polynomials of length k`):
print(`indicates a quasi-polynomial whose value when n is i (mod k)`):
print(`is given by the i-th entry, for i=1..k .`):

print(f(n)=mu):
print(`In Maple input notation  the right side is:`):
lprint(mu):
ku1:=[seq(coeff(taylor(lu,t=0,41),t,i),i=1..40)]:
ku2:=[seq(EvalQPS(mu,n,i),i=1..40)]:
if ku1<>ku2 then
lprint(lu,mu):
print(ku1,ku2):
 print(`Something bad happended`):
else
 print(`For the sake of Sloane, the first 40 terms (starting with n=1) are:`):
print(ku1):
fi:

fi:
print(`-------------------------------------------------------`):


end:

#SeferBPs(K,t): inputs a positive integer and
#outputs all the theorems about partitions where
#the rich are not much richer than the poor.
#Try SeferBPs(5,t);
SeferBPs:=proc(K,t) local k,r,co,t0:

t0:=time():
print(`Enuerating Equitable Partitions of length up to`,K):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`In this article we will give generating functions, as well`):
print(`explicit expressions as sums of quasi-polynomials`):
print(`that enumerate partitions where the sum of the bottom entries`):
print(`exceeds the sum of the top entries, whenever feasible`):
co:=1:
for k from 2 to K do
for r from trunc(k/2)+1 to k-1 do
 BPs1(k,r,t,co):
  co:=co+1:
od:
od:

print(`It took`, time()- t0, `seconds to generate this article.`):
end:











#GFtW(a,Ineqs,x,t): 
#Inputs a list of discrete variables
#a and a list of affine-linear expressions, Ineqs,
#and variable x (for indexed variable) and variable t
#finds the rational generating function, in t,
#and x[a[1]], ..., x[a[nops(a)]],
#whose coefficient of t^n (in the Maclaurin expansion)
#is the weight-enumerator of vectors of length nops(a) that
#add up to n and obey the inequalities, For example try:
#GFtW([a,b,c],[a-b,b-c],t);
GFtW:=proc(a,Ineqs,x,t) local var,Eqs,i,y,gu:
var:=a:
Eqs:=[]:

for i from 1 to nops(Ineqs) do
 var:=[op(var),y[i]]:
 Eqs:=[op(Eqs),Ineqs[i]-y[i]]:
od:

gu:=GFx(var,Eqs,x);

gu:=subs({seq(x[a[i]]=t*x[a[i]],i=1..nops(a)),seq(x[y[i]]=1,i=1..nops(Ineqs))},gu):
factor(gu):
end:


#MagicRectG(r,s,x,z,w): The weight-enumerator of ALL r by s
#arrays of non-negative integers all whose row sums are
#the same and all whose column-sums are the same
#according to the weight x[i,j]^M[i,j]*z^RowSum*w^ColumnSum
#Try:MagicRectG(2,2,x,z,w);
MagicRectG:=proc(r,s,x,z,w) local m,n,i,j,b,Eqs,lu,a:
option remember:
a:=[seq(seq(b[i,j],j=1..s),i=1..r),n,m]:
Eqs:=[seq(add(b[i,j],j=1..s)-n,i=1..r),
seq(add(b[i,j],i=1..r)-m,j=1..s)]:
lu:=GFx(a,Eqs,x);
lu:=subs({seq(seq(x[b[i,j]]=x[i,j],j=1..s),i=1..r)},lu):
lu:=subs({x[n]=z,x[m]=w},lu):
end:


#MagicRect(r,s,z,w): The weight-enumerator of ALL r by s
#arrays of non-negative integers all whose row sums are
#the same and all whose column-sums are the same
#according to the weight z^RowSum*w^ColumnSum
#Try:MagicRect(2,2,z,w);
MagicRect:=proc(r,s,z,w) local gu,x,i,j:
gu:=MagicRectG(r,s,x,z,w):
gu:=subs({seq(seq(x[i,j]=1,j=1..s),i=1..r)},gu):
factor(gu):
end:


#MagicSqG(r,x,z): The weight-enumerator of ALL r by r
#arrays of non-negative integers all whose row sums are
#the same and all whose column-sums are the same
#and they equal to each other
#according to the weight x[i,j]^M[i,j]*z^LineSum
#Try:MagicSqG(2,2,x,z);
MagicSqG:=proc(r,x,z) local n,i,j,b,Eqs,lu,a:
option remember:
a:=[seq(seq(b[i,j],j=1..r),i=1..r),n]:
Eqs:=[seq(add(b[i,j],j=1..r)-n,i=1..r),
seq(add(b[i,j],i=1..r)-n,j=1..r)]:
lu:=GFx(a,Eqs,x);
lu:=subs({seq(seq(x[b[i,j]]=x[i,j],j=1..r),i=1..r)},lu):
lu:=subs({x[n]=z},lu):
end:


#MagicSq(r,z): The weight-enumerator of ALL r by r
#arrays of non-negative integers all whose line-sums are
#the same 
#according to the weight z^LineSum (line=row, column)
#Try:MagicSq(2,z);
MagicSq:=proc(r,z) local gu,x,i,j:
option remember:
gu:=MagicSqG(r,x,z):
gu:=subs({seq(seq(x[i,j]=1,j=1..r),i=1..r)},gu):
factor(gu):
end:

#MagicSqwd(r,z): The weight-enumerator of ALL r by r
#arrays of non-negative integers all whose line-sums are
#the same AND the two main diagonals too
#according to the weight z^LineSum (line=row, column)
#Try:MagicSqwd(2,z);
MagicSqwd:=proc(r,z) local gu,x,i,j:
option remember:
gu:=MagicSqGwd(r,x,z):
gu:=subs({seq(seq(x[i,j]=1,j=1..r),i=1..r)},gu):
factor(gu):
end:


#LHPwt(j,x,t): The weight-enumerator of all Lecture Hall Partitions
#of length j, b[j]+...+b[1] (Lecture Hall partitions, defined
#by  Mireille Bousquet-Melou and Kimmo Eriksson means that
#b[j]/(j-1)>=b[j-1]/(j-1)>= ... >=b[1] >=0)
#according to the weight 
#t^(sum of entries)*x[1]^b[1]*...*x[j]^b[j] 
#Try:
#LHPwt(j,x,t);
LHPwt:=proc(j,x,t) local i,a,Ine,b,lu:
option remember:
a:=[seq(b[i],i=1..j)]:
Ine:=[seq((i-1)*b[i]-i*b[i-1],i=2..j)]:
lu:=GFtW(a,Ine,x,t);
lu:=subs({seq(x[b[i]]=x[i],i=1..j)},lu):
end:



#LHP(j,y,t): The weight-enumerator of all Lecture Hall Partitions
#of length j, b[j]+...+b[1] (Lecture Hall partitions, defined
#by  Mireille Bousquet-Melou and Kimmo Eriksson means that
#b[j]/(j-1)>=b[j-1]/(j-1)>= ... >=b[1] >=0)
#according to the weight 
#t^(sum of entries)*y^(b[1]-b[2]+(-1)^(j-1)) 
#Try:
#LHP(j,y,t);
LHP:=proc(j,y,t) local gu,x,i:
gu:=LHPwt(j,x,t):
gu:=subs({seq(x[j-2*i]=y,i=0..trunc(j/2)),
seq(x[j-2*i-1]=1/y,i=0..trunc((j-1)/2))},gu):
factor(gu):
end:





#NuPe(a,Eqs,N): 
#Inputs a list of discrete variables
#a and a list of affine-linear expressions, Eqs,
#and a positive integer n, outputs the NUMBER 
#of vectors of length nops(a), with non-negative integer
#entries less than N
#the members of Eqs 0
#For example try:
#NuPe([a,b,n],[a+b-n],10);
NuPe:=proc(a,Eqs,N) local gu,i,k,ak,var,Eqs1,co:
option remember:

if nops(a)=1 then
co:=0:
for i from 0 to N do
 if convert(subs(a[1]=i,Eqs),set)={0} then
  co:=co+1:
 fi:
od:
RETURN(co):
fi:

gu:=0:


k:=nops(a):


for ak from 0 to N do

 var:=[op(1..k-1,a)]:

 Eqs1:=subs(a[k]=ak,Eqs):


 gu:=gu+NuPe(var,Eqs1,N):
od:

gu:

end:



#NuMagicSq(r,N): The number of r by r semi-magic squares
#with line-sum N
#Try: NuMagicSqG(2,4);
NuMagicSq:=proc(r,N) local n,i,j,b,Eqs,lu,a:
option remember:
a:=[seq(seq(b[i,j],j=1..r),i=1..r),n]:
Eqs:=[seq(add(b[i,j],j=1..r)-n,i=1..r),
seq(add(b[i,j],i=1..r)-n,j=1..r)]:

NuPe(a,Eqs,N):
end:



###
#CondorcetG(t,c12,c23,c31): the generating function in t
#whose coeff. of t^n (in its Maclaurin expansion)
#gives you the number of possible "score vectors"
#for all six total rankings of the three candidates
#among n votes, that yields the (strict) Condorcet
#scenatrio with margins 
#1 beats 2 with a margin c12
#2 beats 3 with a margin c23
#3 beats 1 with a margin c31
CondorcetG:=proc(t,c12,c23,c31) local p123,p132,p213,p231,p312,p321:
GFt([p123,p132,p213,p231,p312,p321],
[p123+p132+p312-p213-p231-p321-c12,
 p123+p213+p231-p132-p312-p321-c23,
p312+p321+p231-p132-p123-p213-c31],t):
end:



#LHPs(J): tells the story of Lecture Hall Partitions up to J parts
LHPs:=proc(J) local x,y,t,lu,luW,t0,j:

t0:=time():
print(`The Generating Function for Lecture-Hall partitions with`, j, ` parts `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`We will derive generating functions for the so-called`):
print(`Lecture Hall Partitions first defined (and proved!) `):
print(`by Mireille Bousquet-Melou and Kimmo Eriksson `):
print(``):

for j from 2 to J do
 print(`Theorem No. `, j-1, `: The GLOBAL generating function for`):
 print(`for Lecture Hall Partitions with`, j, `parts is:`):
 luW:=LHPwt(j,x,t):
 print(luW):
 print(`and in Maple input notation: `):
 lprint(luW):

 print(`where the weight of a partition [b[j],b[j-1], ..., b[1]] is`):
 print(` x[j]^b[j]*...*x[1]^b[1] `):
 print(`the generating function according to the weight`):
 print(`t^n*y^(b[j]-b[j-1]+b[j-2]+ ... ) is:`):
 lu:=LHP(j,y,t):
 print(lu):
 print(`and in Maple input notation: `):
 lprint(lu):
od:

print(`It took`, time()-t0, `seconds to generate this fascinating web-book`):

end:


#MagicSqGwd(r,x,z): The weight-enumerator of ALL r by r
#arrays of non-negative integers all whose row sums are
#the same and all whose column-sums are the same
#and they equal to each other
#and the two main diagonals are the same
#according to the weight x[i,j]^M[i,j]*z^LineSum
#Try:MagicSqGwd(2,2,x,z);
MagicSqGwd:=proc(r,x,z) local n,i,j,b,Eqs,lu,a:
option remember:
a:=[seq(seq(b[i,j],j=1..r),i=1..r),n]:
Eqs:=[seq(add(b[i,j],j=1..r)-n,i=1..r),
seq(add(b[i,j],i=1..r)-n,j=1..r),
add(b[i,i],i=1..r)-n,add(b[i,r+1-i],i=1..r)-n
]:

lu:=GFx(a,Eqs,x);
lu:=subs({seq(seq(x[b[i,j]]=x[i,j],j=1..r),i=1..r)},lu):
lu:=subs({x[n]=z},lu):
end:

#NuMagicSqwd(r,N): The number of r by r magic squares
#with line-sum N
#Try: NuMagicSwd(2,4);
NuMagicSqwd:=proc(r,N) local n,i,j,b,Eqs,lu,a:
option remember:
a:=[seq(seq(b[i,j],j=1..r),i=1..r),n]:
Eqs:=[seq(add(b[i,j],j=1..r)-n,i=1..r),
seq(add(b[i,j],i=1..r)-n,j=1..r),
add(b[i,i],i=1..r)-n,add(b[i,n+1-i],i=1..r)-n
]:
NuPe(a,Eqs,N):
end:

#Simpson(a,b,t): Inputs positve inetegers, a and b, and a variable t.
#Outputs the generating function, in t,  for Simpson's paradox where
#t^n is the number of 2 by 2 matrices of non-negative integers, [[x11,x12],[x21,x22]] such that
#x11<=a*n, x12<=b*n, x21<=b*n, x22<=a*n, x11/a<x12/b, x21/b<x22/a, yet x11+x21-x12-x22>0
Simpson:=proc(a,b,t) local x11,x12,x21,x22,z,Eqs,gu:
Eqs:=[a*n-x11,b*n-x12,b*n-x21,a*n-x22,a*x12-b*x11-1,b*x22-a*x21-1, x11+x21-x12-x22-1]:
gu:=GFIx([x11,x12,x21,x22,n],Eqs,X):
gu:=subs({X[x11]=1,X[x12]=1,X[x21]=1,X[x22]=1,X[n]=t},gu);
end: