######################################################################
##NaiveCounting.txt: Save this file as  NaiveCounting.txt            #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read NaiveCounting.txt<Enter>                      #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created:  Feb. 16, 2016
 
print(`Created:  Feb. 16, 2016`):
print(` This is  NaiveCounting.txt `):
print(`It is one of the packages that accompany the Class room note `):
print(` Enumerating Patitions Satisfying Linear Inequalties by the Most Naive Method: Counting and Guessing`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 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 Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

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

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: GuessRec, CtoR, Pars  `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: GenF, GenFv, SeqC `):
 print(` `):


elif nops([args])=1 and op(1,[args])=CtoR then
print(`CtoR(S,t), the rational function, in t, whose coefficients`):
print(`are the members of the C-finite sequence S. For example, try:`):
print(`CtoR([[1,1],[1,1]],t);`):

elif nops([args])=1 and op(1,[args])=GenF then
print(`GenF(k,a,C,t,N): inputs a positive integer k, a symbol C, a set of linear inequalities in a[1], ..., a[k], and a variable t.`):
print(`It also inputs a positive integer N (for maximum number of terms to guess)`):
print(`Outputs the ordinary `):
print(`generating function, in the variable t, whose coefficient of t^n is the number of vectors a[1]>=a[2]>=...>=a[k]>=0 that add-up to n`):
print(`and satisfy the set of conditions C. If it can't find anything (because N is not large enough) it returns FAIL. Try:`):
print(`GenF(3,a,{a[1]<=a[2]+a[3]},t,30);`):
print(``):

elif nops([args])=1 and op(1,[args])=GenFv then
print(`GenFv(k,a,C,t,N): Verbose form of GenF(k,a,C,t,N) (q.v.). Try:`):
print(`GenFv(3,a,{a[1]<=a[2]+a[3]},t,30);`):
print(``):

elif nops([args])=1 and op(1,[args])=GuessRec then
print(`GuessRec(L): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant coefficients  `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):

elif nops([args])=1 and op(1,[args])=Pars then
print(`Pars(n,k): the set of vectors [a1,a2,..., ak] such that a1>=a2>=...>=ak>=0 and a1+...+ak=n`):


elif nops([args])=1 and op(1,[args])=SeqC then
print(`SeqC(N,k,a,C): the first N terms of the sequence enumerating the sets of partitions of n (with 0 allowed) of size k,`):
print(`obeying the set of conditions in the set C, phrased in terms of the indexed variable a. For example,try:`):
print(`SeqC(20,3,a,{a[1]<=a[2]+a[3]});`):

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


#CountPartitions


#ez:=proc(): print(` Pars1(n,a1,k), Pars(n,k), IsGood1(p,a,Cond), IsGood(p,a,SetConditions), ParsC(n,k,a,C), SeqC(N,k,a,C) `): 
#print(` Pars1Cslow(n,g,k,a,C), ParsCslow(n,k,a,C) ,  Pars1Cn(n,g,k,a,C), ParsCn(n,k,a,C), GenF(k,a,C,t,N) `):


#######start from Cfinite
#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

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

end:

#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(factor(f)):
fi:
end:


#RtoC(R,t): Given a rational function R(t) finds the
#C-finite description of its sequence of coefficients
#For example, try:
#RtoC(1/(1-t-t^2),t);
RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a:
R1:=normal(R):
D1:=denom(R1):
d:=degree(D1,t):
if degree(numer(R1),t)>=d then
print(`The degree of the top must be less than the degree of the bottom`):
 RETURN(FAIL):
fi:

a:=coeff(D1,t,0):
S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]:
f1:=taylor(R,t=0,d+1):
S1:=[seq(coeff(f1,t,i),i=0..d-1)]:
[S1,S2]:
end:

#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:
#######end from Cfinite

#Pars1(n,g,k): all vectors [a1, ..., ak] with a1>=a2>=..ak>=0 that sum up to n
#and a1=g
Pars1:=proc(n,g,k) local gu,lu,lu1,a2:
option remember:

if k=1 then
 if n=g then
  RETURN({[g]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:

for a2 from 0 to g do
 lu:=Pars1(n-g,a2,k-1):
 gu:=gu union {seq([g,op(lu1)],lu1 in lu)}:
od:

gu:
end:

#Pars(n,k): the set of vectors [a1,a2,..., ak] such that a1>=a2>=...>=ak>=0 and a1+...+ak=n
Pars:=proc(n,k) local g:
option remember:
{seq(op(Pars1(n,g,k)),g=0..n)}:
end:

#IsGood1(p,a,Cond):inputs a partition p of size k, say (so k=nops(p)) and a symbolic condition
#phrased in terms of a[1], ..., a[k] decides whether the partition p satisfies it.
#Try:
#IsGood([8,6,3],a,a[1]<=a[2]+a[3]);
IsGood1:=proc(p,a,Cond)  local k,i:
k:=nops(p):

evalb(subs({seq(a[i]=p[i],i=1..k)},Cond)):

end:


#IsGood(p,a,Conds):inputs a partition p of size k, say (so k=nops(p)) and a set of symbolic conditions
#phrased in terms of a[1], ..., a[k] decides whether the partition p satisfies it.
#Try:
#IsGood([8,6,3],a,{a[1]<=a[2]+a[3]});
IsGood:=proc(p,a,S)  local c:
for c in S do
if not IsGood1(p,a,c) then
 RETURN(false):
fi:
od:
true:
end:



#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

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

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:



#ParsC(n,k,a,C) : the set of weakly decreasing non-negative integers of size k that add-up to n that
#satisy the conditions given in the set C. Try:
#ParsC(5,3,a,{a[1]<=a[2]+a[3]});
ParsC:=proc(n,k,a,C) local gu,mu,mu1:
option remember:
mu:=Pars(n,k):

gu:={}:

for mu1 in mu do
 if IsGood(mu1,a,C) then
  gu:=gu union {mu1}:
 fi:
od:

gu:
end:


#SeqC(N,k,a,C): the first N terms of the sequence enumerating the sets of partitions of n (with 0 allowed) of size k,
#obeying the set of conditions in the set C, phrased in terms of the indexed variable a. For example,try:
#SeqC(20,3,a,{a[1]<=a[2]+a[3]});
SeqC:=proc(N,k,a,C) local n:
[seq(nops(ParsC(n,k,a,C)),n=0..N)]:
end:


#Pars1Cslow(n,g,k,a,C): all vectors [a1, ..., ak] with a1>=a2>=..ak>=0 that sum up to n
#and satisfying the conditions in C
#and a1=g
Pars1Cslow:=proc(n,g,k,a,C) local gu,lu,lu1,a2,C1,i:
option remember:

if k=1 then
 if n=g and IsGood([g],a,C) then
     RETURN({[g]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:

C1:=subs(a[1]=g,C):
C1:=subs({seq(a[i]=a[i-1],i=2..k)},C1):
for a2 from 0 to g do
 lu:=Pars1Cslow(n-g,a2,k-1,a,C1):
 gu:=gu union {seq([g,op(lu1)],lu1 in lu)}:
od:

gu:
end:

#ParsCslow(n,k,a,C): all vectors [a1, ..., ak] with a1>=a2>=..ak>=0 that sum up to n
#and satisfying the conditions in C
ParsCslow:=proc(n,k,a,C) local g:
{seq(op(Pars1Cslow(n,g,k,a,C)),g=0..n)}:
end:

#Pars1Cn(n,g,k,a,C): The number of vectors [a1, ..., ak] with a1>=a2>=..ak>=0 that sum up to n
#and satisfying the conditions in C
#and a1=g
Pars1Cn:=proc(n,g,k,a,C) local a2,C1,i:
option remember:

if k=1 then
 if n=g and IsGood([g],a,C) then
     RETURN(1):
 else
  RETURN(0):
 fi:
fi:


C1:=subs(a[1]=g,C):
C1:=subs({seq(a[i]=a[i-1],i=2..k)},C1):

add(Pars1Cn(n-g,a2,k-1,a,C1),a2=0..g):



end:

#ParsCn(n,k,a,C): all vectors [a1, ..., ak] with a1>=a2>=..ak>=0 that sum up to n
#and satisfying the conditions in C
ParsCn:=proc(n,k,a,C) local g:
add(Pars1Cn(n,g,k,a,C),g=0..n):
end:


#GenF(k,a,C,t,N): inputs a positive integer k, a symbol C, a set of linear inequalities in a[1], ..., a[k], and a variable t.
#It also inputs a positive integer N (for maximum number of terms to guess)
#Outputs the ordinary 
#generating function, in the variable t, whose coefficient of t^n is the number of vectors a[1]>=a[2]>=...>=a[k]>=0 that add-up to n
#and satisfy the set of conditions C. If it can't find anything (because N is not large enough) it returns FAIL. Try:
#GenF(3,a,{a[1]<=a[2]+a[3]},t,30);
GenF:=proc(k,a,C,t,N) local gu,N0,lu:

for N0 from 10 to N by 10 do

gu:=SeqC(N0,k,a,C):

lu:=GuessRec(gu):

if lu<>FAIL then
  RETURN(CtoR(lu,t)):
fi:

od:

FAIL:

end:

 



#GenFv(k,a,C,t,N):  Verbose form of GenFv(k,a,C,t,N) (q.v.) . Try:
#GenFv(3,a,{a[1]<=a[2]+a[3]},t,30);
GenFv:=proc(k,a,C,t,N) local gu,i,c,n,t0:
t0:=time():
gu:=GenF(k,a,C,t,N):

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

print(`Theorem :`):
print(`Let c(n) be the number of weakly-decreasing arrays of non-negative integers of length`, k, [seq(a[i],i=1..k)]):
print(`satisfing the set of conditions `):
print(C):
print(`Then the ordinary generating function of c(n) is given by the following rational  function`):
print(``):
print(Sum(c(n)*t^n,n=0..infinity)=gu):
print(``):
print(`and in Maple foramt it is`):
print(``):
lprint(gu):
print(``):
print(`--------------------------------`):
print(`This ends this theorem that took`, time()-t0, `seconds. `):
end: