######################################################################
##BABUSHKAS: Save this file as  BABUSHKAS                            #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read BABUSHKAS<Enter>                              #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Aug. 2009
 
print(`Created:  Aug. 2009`):
print(` This is BABUSHKAS `):
print(`It accompanies the article `):
print(`The Number of Ways of Breaking Up identical Russian Dolls`):
print(`by  Doron Zeilberger`):
print(`dedicated to Thotsaporn and Thotsaphon Thanatipanonda `):

print(`and also available from Zeilberger's website, and arxiv.org .`):
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(`For a list of the supporting procedures type: ezra1(); `):
print(`For procedures about random generation,type: ezraRand(); `):

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

with(combinat):

ezraRand:=proc()

if args=NULL then
 print(` The procedures about random generation are: LD, RandSPn, RandSPnk `):

else
ezra(args):
fi:
end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Bdirect `):
 print(` Brn, Cdirect, Comps, CompsP(n,k), CompToStates1, DefVec, DefVec1, `):
 print(` F, gammaku,`):
 print(` HafelOper, HafelOperM, Khaber, MonoState, MonoStateC, `):
 print(`Oper, OperC, Pars1, Pars1Ad, ParToM,  ParToM1, ParToU, `):
 print(` SeqBrnDJ, SeqHrnDJ ,States1`):
 print(` Sugim, Wt1, Wt, Yeladim , Yeladim0, Yeladim1, YeladimS `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: B, C, MSSP  `):
 print(` OperD, RatGF, SeqBrn, SeqCrn, SeqC1rn, SeqGrn, SeqHrn,SipurBrn, `):
 print(` SipurRatGF, SSP, Srnk `):



elif nops([args])=1 and op(1,[args])=B then
print(`B(M): Given a multiset 1^m[1]2^m[2]... `):
print(`represented in terms of the multiplicty vector`):
print(`M=[m[1], m[2], ...], finds the number of`):
print(`multisets of sets whose union equals to it.`):
print(`For example, try:`):
print(`B([3,3]);`):

elif nops([args])=1 and op(1,[args])=Bdirect then
print(`Bdirect(M): B(M) by direct counting, only use`):
print(`for checking purposes, where B(M) is not too large`):
print(`For example, try:`):
print(`Bdirect([2,2]);`):

elif nops([args])=1 and op(1,[args])=Brn then
print(`Brn(r,n): the generalized Bell numbers`):
print(`with  r identical Russian Dolls each with n layers.`):
print(`Warning: very slow, only to be used for checking`):
print(`For example, try:`):
print(`Brn(2,5);`):

elif nops([args])=1 and op(1,[args])=C then
print(`C(M): Given a multiset 1^m[1]2^m[2]... `):
print(`represented in terms of the multiplicty vector`):
print(`M=[m[1], m[2], ...], finds the number of`):
print(`sets of sets whose union equals to it.`):
print(`For example, try:`):
print(`C([3,3]);`):

elif nops([args])=1 and op(1,[args])=Cdirect then
print(`Cdirect(M): B(M) by direct counting, only use`):
print(`for checking purposes, where B(M) is not too large`):
print(`For example, try:`):
print(`Cdirect([2,2]);`):

elif nops([args])=1 and op(1,[args])=Comps then
print(` Comps(n,k): the set of lists of length k of non-negative `):
print(` integers that add up to n. `):
print(` For example, try: `):
print(` Comps(4,2);`):


elif nops([args])=1 and op(1,[args])=CompToStates1 then
print(`CompToStates1(C): given a composition C`):
print(`finds all the states corresponding to it`):
print(`(vectors of the same length of partions of C[i], i=1..k) (k=nops(C)):`):
print(`For example, try:`):
print(`CompToStates1([1,1,3,2]);`):


elif nops([args])=1 and op(1,[args])=DefVec then
print(`DefVec(S): the deficieny vector for state S.`):
print(`For example, try:`):
print(`DefVec([0, [[1], 0, [2]]]);`):

elif nops([args])=1 and op(1,[args])=DefVec1 then
print(`DefVec1(k,a,j): the effect on a k-vector`):
print(` (describing the current state)`):
print(`of placing j copies of (n+1) into a bunch of a identical sets`):
print(`For example, try:`):
print(`DefVec1(3,2,1); `):

elif nops([args])=1 and op(1,[args])=F1 then
print(`F1(P,a,n,a0): The value of F1(n;a0) given by the`):
print(`evolution operator P (phrased in terms of a[i]) .`):
print(`For example, try:`):
print(`F1({[1,[-1]]},a,5,[5]);`):


elif nops([args])=1 and op(1,[args])=gammaku then
print(`gammaku(k,u): gamma_k(u) in Theorem 2.1 of `):
print(` J.S. Devitt, and D.M.Jackson's `):
print(` "The enumeration of covers of a finite set"  `):
print(`(J. London Math. Soc. (2), 25(1982), 1-6), for numeric k and`):
print(`symbolic u[1],u[2], ..., u[k]. For example`):
print(`try gammaku(2,u);`):

elif nops([args])=1 and op(1,[args])=GuessRat then
print(`GuessRat(L,t): given a sequence guesses a rational generating`):
print(` function for it. For example, try: `):
print(`GuessRat([1$18],t);`):

elif nops([args])=1 and op(1,[args])=HafelOper then
print(`HafelOper(P,z,D1,k,f): Given a  differential operator`):
print(`P in terms of z[1], ..., z[k] and D1[1], ..., D1[k]`):
print(`applies it to the polynomial f in z[1], ..., z[k] .`):
print(`For example, try:`):
print(`HafelOper(z[1]*D1[1]+z[2]*D1[2],z,D1,2,z[1]*z[2]);`):

elif nops([args])=1 and op(1,[args])=HafelOperM then
print(`HafelOperM(M,z,D1,k,f): Given a monomial differential operator`):
print(`in terms of z[1], ..., z[k] and D1[1], ..., D1[k]`):
print(`applies it to the polynomial f in z[1], ..., z[k]`):
print(`For example, try`):
print(`HafelOperM(z[1]*D1[1],z,D1,1,z[1]^5); `);

elif nops([args])=1 and op(1,[args])=Khaber then
print(`Khaber(v1,v2): addint two vecors, for example, try`):
print(`Khaber([1,2],[3,4]);`):

elif nops([args])=1 and op(1,[args])=LD then
print(`LD(L): loaded die with L`):

elif nops([args])=1 and op(1,[args])=MonoState then
print(`MonoState(S,a,A): the contribution to the recurrence operator`):
print(`coming from state S in terms of a[1], ..., a[k] and`):
print(`the corresonding shift operators A[1], ..., A[k]. For`):
print(`example, try:`):
print(`MonoState([0, [[1], 0, [2]]],a,A); `):

elif nops([args])=1 and op(1,[args])=MSSP then
print(`MSSP(L): given a list of elements finds the set of`):
print(`all multisets of sets whose union (as a multiset, i.e.`):
print(`counting multiplicities) is the multiset given by L.`):
print(`For example, try:`):
print(`MSSP([1,2,3]);`):

elif nops([args])=1 and op(1,[args])=Oper then
print(`Oper(k,a,A): The partial recurrence operator in a[1], ... , a[k]`):
print(`(and their respective shift operators) A[1], ..., A[k] relating.`):
print(`The number of ways of arranging k identical size n Russian Dolls`):
print(`with specification [a[1], ..., a[k]] (i.e. where there are`):
print(`exactly a[i] sets that only show up i times (i=1..k)`):
print(`For example, try:`):
print(`Oper(2,a,A);`):

elif nops([args])=1 and op(1,[args])=OperC then
print(`OperC(k,a): like Oper(k,a,A), but in more convenient form `):
print(`for Maple. `):
print(`For example, try:`):
print(`OperC(2,a);`):

elif nops([args])=1 and op(1,[args])=OperD then
print(`OperD(k,z,D1): The differential operator repsonsible `):
print(`for the evlution of k identical Russian dolls`):
print(`in terms of indexed variables z[1],z[2], ...`):
print(`and the corresponding differential operators D1[1],D1[2], ...`):
print(`For example, try:`):
print(`OperD(2,z,D1);`):

elif nops([args])=1 and op(1,[args])=Pars1 then
print(`Pars1(n,k): all the partitions of n with largest part k`):
print(`For example, try: Pars1(5,2);`):

elif nops([args])=1 and op(1,[args])=Pars1Ad then
print(`Pars1Ad(n,k): all the partitions of n with largest part<=k`):
print(`For example, try: Pars1Ad(5,2);`):

elif nops([args])=1 and op(1,[args])=ParToM then
print(`ParToM(L): the integer partition L in multiplicity notation.`):
print(`For example, try: `):
print(`ParToM([2,2,2,1,1,1,1]);`):

elif nops([args])=1 and op(1,[args])=ParToM1 then
print(`ParToM1(L,r): given a partition L with largest part<=r`):
print(`translates it to a vector [a[1],..., a[r]]`):
print(`where a[i] is the number of times i shows up`):
print(`For example, try: `):
print(`ParToM1([2,2,2,1,1,1,1],2);`):

elif nops([args])=1 and op(1,[args])=ParToU then
print(`ParToU(L,k): the integer partition L in multiplicity notation`):
print(`1^u[1]...k^u[k] only returns [u[1],..., u[k]].`):
print(`For example, try: `):
print(`ParToU([2,2,2,1,1,1,1],3);`):

elif nops([args])=1 and op(1,[args])=RandSPn then
print(`RandSPn(n): a (uniformly) random set-partition of {1,...,n}.`):
print(` For example, try: `):
print(`RandSPn(3);`):
print(`add(x[RandSPn(3)],i=1..500);`):

elif nops([args])=1 and op(1,[args])=RandSPnk then
print(` RandSPnk(n,k): picks a randon set partition `):
print(` of {1, ..., n} with k sets. `):
print(`For example, try: RandSPnk(4,2); `):

elif nops([args])=1 and op(1,[args])=RatGF then
print(`RatGF(r,N,t): the ordinary generating functions in z`):
print(`for the number of ways of rearranging r identical Russian Dolls`):
print(`into exactly k other (possibly repeated) Russian Dolls`):
print(`for k=1,2,..., guessing from the data up to N terms.`):
print(`Note: this is empirical, but it can be easily proved rigorously`):
print(`(but who cares?!, we have better things to do!)`):
print(`For example, try:`):
print(`RatGF(2,30,z);`):

elif nops([args])=1 and op(1,[args])=SeqBrn then
print(`SeqBrn(r,n): the sequence of length n+1 whose`):
print(`(i+1)^th item is the number of ways of breaking-up`):
print(`r identical Russian Dolls, each consisting of i layers`):
print(`into smaller Russian Dolls. For example, try:`):
print(`SeqBrn(2,10);`):

elif nops([args])=1 and op(1,[args])=SeqBrnDJ then
print(`SeqBrnDJ(r,n): same as SeqBrn(r,n,t)`):
print(`but using the Devitt-Jackson formula .`):
print(` For example, try:`):
print(`SeqBrnDJ(2,10);`):

elif nops([args])=1 and op(1,[args])=SeqCrn then
print(`SeqCrn(r,n): the sequence of length n+1 whose`):
print(`(i+1)^th item is the number of ways of breaking-up`):
print(`r identical Russian Dolls, each consisting of i layers`):
print(`into smaller ALL DIFFERENT Russian Dolls. For example, try`):
print(`SeqCrn(2,10);`):

elif nops([args])=1 and op(1,[args])=SeqC1rn then
print(`SeqC1rn(r,n,t): the sequence of length n+1 whose`):
print(`(i+1)^th item is the polynomial in t, whose coeff. of t^j is`):
print(`the number  of ways of breaking-up`):
print(`r identical Russian Dolls, each consisting of i layers`):
print(`into j DIFFERENT Russian Dolls. For example, try`):
print(`SeqC1rn(2,10,t);`):

elif nops([args])=1 and op(1,[args])=SeqGrn then
print(`SeqGrn(r,n,t): the sequence of length n+1 whose`):
print(`(i+1)^th item is the polynomial in t `):
print(`whose coefficient of t^j is the number of ways`):
print(` of breaking-up`):
print(`r identical Russian Dolls, each consisting of i layers`):
print(`into j Different kinds of Russian Dolls (i.e. you don't count`):
print(`duplications). For example, try:`):
print(`SeqGrn(2,10,t);`):

elif nops([args])=1 and op(1,[args])=SeqHrn then
print(`SeqHrn(r,n,t): the sequence of length n+1 whose`):
print(`(i+1)^th item is the polynomial in t `):
print(`whose coefficient of t^j is the number of ways`):
print(` of breaking-up`):
print(`r identical Russian Dolls, each consisting of i layers`):
print(`into j smaller Russian Dolls. `):
print(` For example, try:`):
print(`SeqHrn(2,10,t);`):

elif nops([args])=1 and op(1,[args])=SeqHrnDJ then
print(`SeqHrnDJ(r,n,t): same as SeqHrn(r,n,t)`):
print(`but using the Devitt-Jackson formula .`):
print(` For example, try:`):
print(`SeqHrnDJ(2,10,t);`):


elif nops([args])=1 and op(1,[args])=SipurBrn then
print(`SipurBrn(R,N): all the sequenecs, up to N+1 terms`):
print(`of the total number of ways of breaking up r identical`):
print(`Russian Dolls into smaller dolls, for r=1..N .`):
print(`For example, try:`):
print(`SipurBrn(3,30);`):

elif nops([args])=1 and op(1,[args])=SipurRatGF then
print(`SipurRatGF(R,N,t): the story for the ordinary generating `):
print(` functions in z `):
print(`for the number of ways of rearranging r identical Russian Dolls`):
print(`into exactly k other (possibly repeated) Russian Dolls`):
print(`(for r=1..R)`):
print(`for k=r,r+1,..., guessing from the data up to N terms.`):
print(`Note: this is empirical, but it can be easily proved rigorously`):
print(`(but who cares?!, we have better things to do!).`):
print(`For example, try:`):
print(`SipurRatGF(2,30,z);`):

elif nops([args])=1 and op(1,[args])=Srnk then
print(`Srnk(r,n,k): the generalized Stirling numbers of the second kind`):
print(`with  r identical Russian Dolls each with n layers and forming`):
print(`k resulting dolls. For example, try:`):
print(`Srnk(2,5,3);`):

elif nops([args])=1 and op(1,[args])=SSP then
print(`SSP(L): given a list of elements finds the set of`):
print(`all sets of sets whose union (as a multiset, i.e.`):
print(`counting multiplicities) is the multiset given by L.`):
print(`For example, try:`):
print(`SSP([1,2,3]);`):

elif nops([args])=1 and op(1,[args])=States1 then
print(`States1(k): all the possible states with k copies of {n}`):
print(`For example, try:`):
print(`States1(2);`):


elif nops([args])=1 and op(1,[args])=Sugim then
print(`Sugim(k,r): all vectors of non-negative integers`):
print(`such that a[1]+2*a[2]+...+r*a[r]=k. For example try:`):
print(`Sugim(4,2);`):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(a,S): The weight of state S as a monomial in a[1], ..., a[k].`):
print(`For example, try:`):
print(`Wt(a,[0, [[1], 0, [2]]]);`):

elif nops([args])=1 and op(1,[args])=Wt1 then
print(`Wt1(a,p): The atomic contribution of integer-partition p`):
print(`to a compartment of (symbolic) size a.`):
print(`For example, try: `):
print(`Wt1(a,[2,2,1]);`):

elif nops([args])=1 and op(1,[args])=Yeladim then
print(`Yeladim(S,a): Given a multiset of sets S `):
print(`finds all multisets of sets obtained by adding`):
print(`another a to it, either by itself, or in an already existing set.`):
print(`For example, try:`):
print(`Yeladim({[{1,2},2]},3);`):

elif nops([args])=1 and op(1,[args])=Yeladim0 then
print(`Yeladim0(S,a): Given a multiset of sets S `):
print(`and an element a, returns the multiset of`):
print(`sets obtained by adding the singleton {a}.`):
print(`For example, try:`):
print(`Yeladim0({[{1,2},2]},2);`):

elif nops([args])=1 and op(1,[args])=Yeladim1 then
print(`Yeladim1(S,s,a): Given a multiset of sets S `):
print(`and one of the sets in S, and an element a, returns the multiset of`):
print(`sets obtained by putting a in one of the s's in S .`):
print(`For example, try:`):
print(`Yeladim1({[{1,2},2]},{1,2},3);`):

elif nops([args])=1 and op(1,[args])=YeladimS then
print(`YeladimS(G,a): Given a set of multisets G`):
print(`finds the union of all multisets of sets obtained by adding`):
print(`another {a} to it for each and every single member of G.`):
print(`For example, try:`):
print(`YeladimS({{[{1,2},2]}},3);`):

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


#####Stuff from Findrec
###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 (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
RETURN(FAIL):
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:

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:

#Findrec0(f,N): 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);
Findrec0:=proc(f,N) local ORDER,n,ope:
for ORDER from 0 to nops(f)-8 do
 ope:=findrec(f,0,ORDER,n,N):
 if ope<>FAIL then
  RETURN(ope):
 fi:
od:
FAIL:

end:

#GuessRat(L,t): given a sequence guesses a rational generating function
#for it. For example, try: GuessRat([1$18],t);
GuessRat:=proc(L,t) local L1,i,N,mekh,ope,lu,mone,lu1:

if L=[0$(nops(L))] then
 RETURN(0):
fi:

for i from 1 while L[i]=0 do od:

L1:=[op(i..nops(L),L)]:


ope:=Findrec0(L1,N):

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

mekh:=numer(normal(subs(N=1/t,ope))):

lu:=add(L[i+1]*t^i,i=0..nops(L)-1):
mone:=expand(lu*mekh):
mone:=add(coeff(mone,t,i)*t^i,i=0..nops(L)-3):
lu:=mone/mekh:
lu1:=taylor(lu,t=0,nops(L)+1):
if not [seq(coeff(lu1,t,i),i=0..nops(L)-1)]=L then
  RETURN(FAIL):
else
 RETURN(factor(lu)):
fi:

end:


#####end stuff from Findrec

#Comps(n,k): the set of lists of length k of non-negative
#integers that add up to n
#For example, try:
#Comps(4,2);
Comps:=proc(n,k) local gu,ak,lu,L:
option remember:

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

gu:={}:

for ak from 0 to n do
lu:=Comps(n-ak,k-1):
gu:=gu union {seq([op(L),ak],L in lu)}:
od:
gu:
end:




#CompsP(n,k): the set of pairs [num,C], with num+sum(C)=n and |C|=k
#CompsP(4,2);
CompsP:=proc(n,k) local n1,gu,lu,L:
option remember:
gu:={}:
for n1 from 0 to n do
 lu:=Comps(n-n1,k):
 gu:=gu union {seq([n1,L],L in lu)}:
od:
gu:
end:

#Pars1(n,k): all the partitions of n with largest part k
Pars1:=proc(n,k) local lu,L:
option remember:
if n<1 then
  RETURN({}):
fi:

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

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

if k>n then
 RETURN({}):
fi:

lu:=Pars1Ad(n-k,k):
{seq([k,op(L)],L in lu)}:
end:

#Pars1Ad(n,k): all the partitions of n with largest part <=k
Pars1Ad:=proc(n,k) local k1:
option remember:
{seq(op(Pars1(n,k1)),k1=1..k)}:
end:


#CompToStates1(C): given a composition C
#finds all the states corresponding to it
#(vectors of the same length of partions of C[i], i=1..k) (k=nops(C)):
#For example, try:
#CompToStates1([1,1,3,2]);
CompToStates1:=proc(C) local gu,k,C1,g,mu,m:
k:=nops(C):
if k=1 then
 if C[1]=0 then
  RETURN({[0]}):
 else
  RETURN({[[1$C[1]]]}):
 fi:
fi:

C1:=[op(1..k-1,C)]:

gu:=CompToStates1(C1):

if C[k]=0 then
 RETURN({seq([op(g),0],g in gu)}):
fi:

mu:=Pars1Ad(C[k],k):

{seq(seq([op(g),m],g in gu),m in mu)}:
end:

#States1(k): all the possible states with k copies of {n}
#For example, try:
#States1(2);
States1:=proc(k) local mu,gu,m,lu,L:
mu:=CompsP(k,k):
gu:={}:
for m in mu do
 lu:=CompToStates1(m[2]):
 gu:=gu union {seq([m[1],L],L in lu)}:
od:
gu:
end:


#DefVec1(k,a,j): the effect on a k-vector (describing the current state)
#of placing j copies of (n+1) into a bunch of a identical sets
DefVec1:=proc(k,a,j) local i1,T:
for i1 from 1 to k do
T[i1]:=0:
od:

T[a]:=T[a]-1:
T[j]:=T[j]+1:
if a>j then
T[a-j]:=T[a-j]+1:
fi:

[seq(T[i1],i1=1..k)]:
end:


#Khaber(v1,v2): addint two vecors, for example, try
#Khaber([1,2],[3,4]);
Khaber:=proc(v1,v2) local i: 
if nops(v1)<>nops(v2) then RETURN(FAIL): fi:
[seq(v1[i]+v2[i],i=1..nops(v1))]:end:



#DefVec(S): the deficieny vector for state S
#For example, try:
#DefVec([0, [[1], 0, [2]]]);
DefVec:=proc(S) local a0,C,i,gu,p,j,k:
a0:=S[1]:
C:=S[2]:
k:=nops(C):
gu:=[0$k]:

for i from 1 to k do
p:=C[i]:

 if p<>0 then
  for j from 1 to nops(p) do
     gu:=Khaber(gu,DefVec1(k,i,p[j])):
  od:
 fi:

od:

if a0>=1 then
gu:=Khaber(gu,[0$(a0-1),1,0$(k-a0)]):
fi:
gu:

end:


#ParToM(L): the integer partition L in multiplicity notation
#For example, try: 
#ParToM([2,2,2,1,1,1,1]);
ParToM:=proc(L) local i,L1:
option remember:
if L=[] then 
 RETURN([]):
fi:

for i from 2 to nops(L) while L[i]=L[1] do od:
i:=i-1:
L1:=[op(i+1..nops(L),L)]:
[ [L[1],i], op(ParToM(L1))]:
end:


#Wt1(a,p): The atomic contribution of integer-partition p
#to a compartment of (symbolic) size a
#For example, try
#Wt1(a,[2,2,1]);
Wt1:=proc(a,p) local p1,i:
p1:=ParToM(p):
p1:=[seq(p1[i][2],i=1..nops(p1))]:
simplify(a!/(a-convert(p1,`+`))!/mul(p1[i]!,i=1..nops(p1))):
end:



#Wt(a,S): The weight of state S as a monomial in a[1], ..., a[k]
#For example, try:
#Wt(a,[0, [[1], 0, [2]]]);
Wt:=proc(a,S) local C,k,mu,i:
C:=S[2]:
k:=nops(C):
mu:=1:
for i from 1 to k do
 if C[i]<>0 then
  mu:=mu*Wt1(a[i],C[i]):
 fi:
od:
mu:
end:



#MonoState(S,a,A): the contribution to the recurrence operator
#coming from state S in terms of a[1], ..., a[k] and
#the corresonding shift operators A[1], ..., A[k]. For
#example, try:
#MonoState([0, [[1], 0, [2]]],a,A);
MonoState:=proc(S,a,A) local v,w,i,k:
k:=nops(S[2]):
v:=DefVec(S):
w:=Wt(a,S):
w:=subs({seq(a[i]=a[i]-v[i],i=1..k)},w):
w*mul(A[i]^(-v[i]),i=1..k):
end:



#MonoStateC(S,a): the contribution to the recurrence operator
#coming from state S in terms of a[1], ..., a[k] and
#the corresonding shift operators A[1], ..., A[k]. For
#example, try:
#MonoStateC([0, [[1], 0, [2]]],a);
MonoStateC:=proc(S,a) local v,w,i,k:
k:=nops(S[2]):
v:=DefVec(S):
w:=Wt(a,S):
w:=subs({seq(a[i]=a[i]-v[i],i=1..k)},w):
[w,[seq(-v[i],i=1..k)]]:
end:

#Oper(k,a,A): The partial recurrence operator in a[1], ... , a[k]
#(and their respective shift operators) A[1], ..., A[k] relating
#The number of ways of arranging k identical size n Russian Dolls
#with specification [a[1], ..., a[k]] (i.e. where there are
#exactly a[i] sets that only show up i times (i=1..k)
#For example, try:
#Oper(2,a,A);
Oper:=proc(k,a,A) local gu,g:
option remember:
gu:=States1(k):
add(MonoState(g,a,A), g in gu):
end:

#OperC(k,a): Like Oper(k,a,A) but in a more convenient form
#For example, try:
#OperC(2,a);
OperC:=proc(k,a) local gu,g:
option remember:
gu:=States1(k):
{seq(MonoStateC(g,a), g in gu)}:
end:


#F1(P,a,n,a0): The value of F1(n;a0) given by the
#evolution operator P (phrased in terms of a[i], i=1..r)
#For example, try:
#F1({[1,[-1]]},a,5,[5]);
F1:=proc(P,a,n,a0) local r,p,i:
option remember:
r:=nops(a0):
if n=0 then
 if a0=[0$r] then
   RETURN(1):
  else
    RETURN(0):
 fi:
fi:

add(subs({seq(a[i]=a0[i],i=1..r)},p[1])*F1(P,a,n-1,Khaber(a0,p[2])),
p in P):

end:




#ParToM1(L,r): given a partition L with largest part<=r
#translates it to a vector [a[1],..., a[r]]
#where a[i] is the number of times i shows up
#For example, try: 
#ParToM1([2,2,2,1,1,1,1],2);
ParToM1:=proc(L,r) local L1,T,i:
option remember:
L1:=ParToM(L):

for i from 1 to r do
T[i]:=0:
od:

for i from 1 to nops(L1) do
T[L1[i][1]]:=L1[i][2]:
od:

[seq(T[i],i=1..r)]:

end:


#Sugim(k,r): all vectors of non-negative integers
#such that a[1]+2*a[2]+...+r*a[r]=k. For example try:
#Sugim(4,2);
Sugim:=proc(k,r) local gu,g:
gu:=Pars1Ad(k,r):

{seq(ParToM1(g,r), g in gu)}:
end:


#Srnk(r,n,k): the generalized Stirling numbers of the second kind
#with  r identical Russian Dolls each with n layers and forming
#k resulting dolls. For example, try:
#Srnk(2,5,3);
Srnk:=proc(r,n,k) local gu,P1,a,g:
option remember:
P1:=OperC(r,a):
gu:=Sugim(k,r):
add(F1(P1,a,n,g), g in gu):
end:


#Brn(r,n): the generalized Bell numbers
#with  r identical Russian Dolls each with n layers and forming
#For example, try:
#Brn(2,5);
Brn:=proc(r,n) local k:
option remember:
add(Srnk(r,n,k),k=1..n*r):
end:




#MonoStateD(S,z,D1): the contribution to the differential operator
#coming from state S in terms of z[1], ..., z[k] and
#the corresonding diff. operators D1[1], ..., D1[k]. For
#example, try:
#MonoStateD([0, [[1], 0, [2]]],z,D1);
MonoStateD:=proc(S,z,D1) local w,i,k,lu,j,P1,j1:
k:=nops(S[2]):
w:=z[S[1]]:
for i from 1 to k do
lu:=S[2][i]:
 if lu<>0 then
 for j from 1 to nops(lu) do
  w:=w*z[i-lu[j]]*z[lu[j]]*D1[i]:
 od:
  P1:=ParToM(lu):
  w:=w/mul(P1[j1][2]!,j1=1..nops(P1)):
 fi:
od:
w:=subs(z[0]=1,w):

end:

#OperD(k,z,D1): The differential operator for the evolution
#of k identical Russian Dolls
#OperD(2,z,D1);
OperD:=proc(k,z,D1) local gu,g:
option remember:
gu:=States1(k):
add(MonoStateD(g,z,D1), g in gu):
end:



#HafelOperM(M,z,D1,k,f): Given a monomial differential operator
#in terms of z[1], ..., z[k] and D1[1], ..., D1[k]
#applies it to the polynomial f in z[1], ..., z[k]
#For example, try:
#HafelOperM(z[1]*D1[1],z,D1,1,z[1]^5);
HafelOperM:=proc(M,z,D1,k,f) local mu,i,d1,f1:
f1:=f:
mu:=M:

for i from 1 to k do
 d1:=degree(M,D1[i]):
 if d1>0 then
  f1:=diff(f1,z[i]$d1):
  mu:=normal(mu/D1[i]^d1):
 fi:
od:
mu:=normal(mu):
expand(mu*f1):
end:


#HafelOper(P,z,D1,k,f): Given a  differential operator
#P in terms of z[1], ..., z[k] and D1[1], ..., D1[k]
#applies it to the polynomial f in z[1], ..., z[k]
#For example, try:
#HafelOper(z[1]*D1[1]+z[2]*D1[2],z,D1,2,z[1]*z[2]);
HafelOper:=proc(P,z,D1,k,f) local i:
if type(P,`+`) then
 expand(add(HafelOperM(op(i,P),z,D1,k,f),i=1..nops(P))):
else
HafelOperM(P,z,D1,k,f):
fi:
end:

 


#SeqBrn(r,n): the sequence of length n+1 whose
#(i+1)^th item is the number of ways of breaking-up
#r identical Russian Dolls, each consisting of i layers
#into smaller Russian Dolls
#SeqBrn(2,10);
SeqBrn:=proc(r,n) local gu,mu,z,D1,P,i,j:
P:=OperD(r,z,D1):
mu:=1:
gu:=[1]:

for i from 1 to n do
 mu:=HafelOper(P,z,D1,r,mu):
 gu:=[op(gu),subs({seq(z[j]=1,j=1..r)},mu)]:
od:
gu:
end:



SeqGrn:=proc(r,n,t) local gu,mu,z,D1,P,i,j:
P:=OperD(r,z,D1):
mu:=1:
gu:=[1]:

for i from 1 to n do
 mu:=HafelOper(P,z,D1,r,mu):
 gu:=[op(gu),subs({seq(z[j]=t,j=1..r)},mu)]:
od:
gu:
end:

#SeqHrn(r,n,t): the sequence of length n+1 whose`):
#(i+1)^th item is the polynomial in t `):
#whose coefficient of t^j is the number of ways`):
#of breaking-up`):
#r identical Russian Dolls, each consisting of i layers`):
#into j smaller Russian Dolls. `):
#For example, try:`):
#SeqHrn(2,10,t);`):
SeqHrn:=proc(r,n,t) local gu,mu,z,D1,P,i,j:
option remember:
P:=OperD(r,z,D1):
mu:=1:
gu:=[1]:

for i from 1 to n do
 mu:=HafelOper(P,z,D1,r,mu):
 gu:=[op(gu),subs({seq(z[j]=t^j,j=1..r)},mu)]:
od:
gu:
end:




#SipurBrn(R,N): all the sequenecs, up to N+1 terms
#of the total number of ways of breaking up r identical
#Russian Dolls into smaller dolls, for r=1..N .
#For example, try:
#SipurBrn(3,30);
SipurBrn:=proc(R,N) local r:
print(`This is the story about the number of ways of rearanging`):
print(`r identical Russian Dolls , each with n atomic dolls, into`):
print(`smaller Russian Dolls, for r between 1 and `, R):
print(`and for n  between 0 and `, N):

for r from 1 to 1 do
 print(`The enumerating sequence for the number of ways of breaking-up`):
 print(`one Russian Doll with i atomic dolls for i between 0 `):
 print(`and `, N, `equals `):
 print(SeqBrn(r,N)):
 print(`Note that this rings a Bell!`):
od:

for r from 2 to R do
 print(`The enumerating sequence for the number of ways of breaking-up`):
 print(r, `identical Russian Dolls each with i atomic dolls, for i between 0`):
 print(`and `, N, `equals `):
 print(SeqBrn(r,N)):
od:

end:


#RatGF(r,N,t): the ordinary generating functions in z
#for the number of ways of rearranging r identical Russian Dolls
#into exactly k other (possibly repeated) Russian Dolls
#for k=r,r+1,..., guessing from the data up to N terms.
#Note: this is empirical, but it can be easily proved rigorously
#(but who cares?!, we have better things to do!)
#For example, try:
#RatGF(2,30,z);
RatGF:=proc(r,N,z) local gu,mu,i,lu,i1,t,bu:
lu:=SeqHrn(r,N,t):
gu:=[]:
for i from r do
mu:=[seq(coeff(lu[i1],t,i),i1=1..nops(lu))]:
bu:=GuessRat(mu,z):
 if bu=FAIL then
   RETURN(gu):
else
 gu:=[op(gu),bu]:
 fi:
od:
end:






#SipurRatGF(R,N,t): the story for the ordinary generating functions in z
#for the number of ways of rearranging r identical Russian Dolls
#into exactly k other (possibly repeated) Russian Dolls
#(for r=1..R)
#for k=r,r+1,..., guessing from the data up to N terms.
#Note: this is empirical, but it can be easily proved rigorously
#(but who cares?!, we have better things to do!)
#For example, try:
#SipurRatGF(2,30,z);
SipurRatGF:=proc(R,N,z) local r,mu:
print(`This is the story for the ordinary generating functions`):
print(`for the number of ways of breaking up r identical Russian Dolls`):
print(`into k smaller Russian Dolls, according to the number`):
print(`of atomic dolls, for k small.`):
print(`Note, these were empirically guessed, but it can all be proved`):
print(`rigorously, but frankly, we don't care!`):
for r from 1 to R do
mu:=RatGF(r,N,z):
print(`The first few generating functions for breaking`, r, `identical`):
print(`Russian Dolls into k smaller Russian Dolls, for k from 1 to`, nops(mu)):
print(`using the variable`, z, `happen to be `):
print(mu):
od:

end:


#SeqCrn(r,n): the sequence of length n+1 whose
#(i+1)^th item is the number of ways of breaking-up
#r identical Russian Dolls, each consisting of i layers
#into smaller ALL DIFFERENT Russian Dolls. For example, try
#SeqCrn(2,10);
SeqCrn:=proc(r,n) local gu,mu,z,D1,P,i,j:
P:=OperD(r,z,D1):
mu:=1:
gu:=[1]:

for i from 1 to n do
 mu:=HafelOper(P,z,D1,r,mu):
 gu:=[op(gu),subs({z[1]=1,seq(z[j]=0,j=2..r)},mu)]:
od:
gu:
end:

#SeqC1rn(r,n,t): the sequence of length n+1 whose
#(i+1)^th item is the polynomial in t whose coeff. of t^j is
#the number of ways of breaking-up
#r identical Russian Dolls, each consisting of i layers
#into j DIFFERENT Russian Dolls. For example, try
#SeqC1rn(2,10,t);
SeqC1rn:=proc(r,n,t) local gu,mu,z,D1,P,i,j:

mu:=1:
gu:=[1]:
P:=OperD(r,z,D1):
for i from 1 to n do
 mu:=HafelOper(P,z,D1,r,mu):
 gu:=[op(gu),subs({z[1]=t,seq(z[j]=0,j=2..r)},mu)]:
od:
gu:
end:



#B(M): Given a multiset 1^m[1]2^m[2]... 
#represented in terms of the multiplicty vector
#M=[m[1], m[2], ...], finds the number of
#multisets of sets whose union equals to it
#For example, try:
#B([3,3]);
B:=proc(M) local mu,r,z,D1,j,P,i:

mu:=1:

for r from 1 to nops(M) do
P:=OperD(r,z,D1):
for i from 1 to M[r] do
 mu:=HafelOper(P,z,D1,r,mu):
od:
od:
subs({seq(z[j]=1,j=1..r)},mu):
end:


#Yeladim0(S,a): Given a multiset of sets S 
#and an element a, returns the multiset of
#sets obtained by adding the singleton {a}
#For example, try:
#Yeladim0({[{1,2},2]},2);
Yeladim0:=proc(S,a) local T,gu,i,g:
gu:={}:
for i from 1 to nops(S) do
 T[S[i][1]]:=S[i][2]:
 gu:=gu union {S[i][1]}:
od:

if not member({a},gu) then
 T[{a}]:=0:
fi:

gu:=gu union {{a}}:

T[{a}]:=T[{a}]+1:

{seq([g,T[g]],g in gu)}:

end:





#Yeladim1(S,s,a): Given a multiset of sets S 
#and one of the sets in S, and an element a, returns the multiset of
#sets obtained by putting a in one of the s's in S .
#For example, try:
#Yeladim1({[{1,2},2]},{1,2},3);
Yeladim1:=proc(S,s,a) local T,gu,i,g,s1,mu:

if member(a,s) then
 RETURN(FAIL):
fi:

gu:={}:
for i from 1 to nops(S) do
 T[S[i][1]]:=S[i][2]:
 gu:=gu union {S[i][1]}:
od:

if not member(s,gu) then
 RETURN(FAIL):
fi:

s1:= s union {a}:

if not member(s1,gu) then
 T[s1]:=0:
 gu:=gu union {s1}:
fi:

T[s1]:=T[s1]+1:
T[s]:=T[s]-1:
mu:={}:

for g in gu do
 if T[g]>0 then
   mu:=mu union {[g,T[g]]}:
 fi:
od:

mu:
end:

#Yeladim(S,a): Given a multiset of sets S 
#finds all multisets of sets obtained by adding
#another {a} to it
#Yeladim({[{1,2},2]},3);
Yeladim:=proc(S,a) local mu,gu,s,g,lu:

gu:={seq(s[1],s in S)}:

mu:={Yeladim0(S,a)}:

for g in gu do
lu:=Yeladim1(S,g,a):
 if lu<>FAIL then
  mu:=mu union {lu}:
 fi:
od:
mu:
end:



#YeladimS(G,a): Given a set of multisets G
#finds the union of all multisets of sets obtained by adding
#another {a} to it for each and every single member of G.
#For example, try:
#YeladimS({{[{1,2},2]}},3);
YeladimS:=proc(G,a) local S:
{seq(op(Yeladim(S,a)), S in G)}:
end:


#MSSP(L): given a list of elements finds the set of
#all multisets of sets whose union (as a multiset, i.e.
#counting multiplicities) is the multiset given by L.
#For example, try:
#MSSP([1,2,3]);
MSSP:=proc(L) local i,gu:
gu:={{}}:
for i from 1 to nops(L) do
 gu:=YeladimS(gu,L[i]):
od:
gu:
end:



#Bdirect(M): B(M) computed directly, by actually constructing
#the object being counted. Only used as a check
#For example, try:
#Bdirect([2,2]);
Bdirect:=proc(M) local mu,i,c,j:
mu:=[]:
c:=0:
for i from 1 to nops(M) do
 for j from 1 to M[i] do
  c:=c+1:
  mu:=[op(mu),c$i]:
 od:
od:
nops(MSSP(mu)):
end:


#C(M): Given a multiset 1^m[1]2^m[2]... 
#represented in terms of the multiplicty vector
#M=[m[1], m[2], ...], finds the number of
#sets of sets whose union equals to it
#For example, try:
#C([3,3]);
C:=proc(M) local mu,r,z,D1,j,P,i:

mu:=1:

for r from 1 to nops(M) do
P:=OperD(r,z,D1):
for i from 1 to M[r] do
 mu:=HafelOper(P,z,D1,r,mu):
od:
od:
subs({z[1]=1,seq(z[j]=0,j=2..r)},mu):
end:


#IsProper(M): is the multi-set set partition M proper?
IsProper:=proc(M) local m:
evalb({seq(evalb(m[2]=1),m in M)}={true}):
end:



#SSP(L): given a list of elements finds the set of
#all sets of sets whose union (as a multiset, i.e.
#counting multiplicities) is the multiset given by L.
#For example, try:
#SSP([1,2,3]);
SSP:=proc(L) local mu,gu,m:
mu:=MSSP(L):
gu:={}:
for m in mu do
if IsProper(m) then
  gu:=gu union {m}:
fi:
od:
gu:
end:


#Cdirect(M): C(M) computed directly, by actually constucting
#the object being counted. Only used as a check
#For example, try:
#Cdirect([2,2]);
Cdirect:=proc(M) local mu,i,c,j:
mu:=[]:
c:=0:
for i from 1 to nops(M) do
 for j from 1 to M[i] do
  c:=c+1:
  mu:=[op(mu),c$i]:
 od:
od:
nops(SSP(mu)):
end:


#gammaku(k,u): gamma_k(u) in Theorem 2.1 of J.S. Devitt, and D.M.Jackson's
#"The enumeration of covers of a finite set" 
#(J. London Math. Soc. (2), 25(1982), 1-6), for numeric k and
#symbolic u[1],u[2], ..., u[k]. For example
#try gammaku(2,u);
gammaku:=proc(k,u) local i,z,gu:
gu:=mul((1+z^i)^u[i],i=1..k):
gu:=taylor(gu,z=0,k+1):
expand(coeff(gu,z,k)):
end:








#ParToU(L,k): the integer partition L in multiplicity notation
#1^u[1]...k^u[k] only returns [u[1],..., u[k]]
#For example, try: 
#ParToU([2,2,2,1,1,1,1],3);
ParToU:=proc(L,k) local mu,T,i:
mu:=ParToM(L):

for i from 1 to k do
T[i]:=0:
od:

for i from 1 to nops(mu) do
 T[mu[i][1]]:=mu[i][2]:
od:

[seq(T[i],i=1..k)]:

end:






#SeqHrnDJ(r,n,t): 
#Like SeqHrn(r,n,t) but using the generating function
#of J.S.Devitt and D.M. Jackson
#For example, try:
#SeqHrnDJ(2,10,t);`):
SeqHrnDJ:=proc(r,n,t) local y,gam,u,mu,k,lu,kv,pa1,ua1,vu,i:
gam:=gammaku(r,u):
lu:=exp(-add(t^i/i,i=1..r)):
mu:=0:
for k from 0 to r*n do
 kv:=Pars1Ad(k,r):
  vu:=0:
 for pa1 in kv do
  ua1:=ParToU(pa1,r):
   vu:=vu+1/mul(i^ua1[i],i=2..r)/mul(ua1[i]!,i=1..r)*
    exp(y*subs({seq(u[i]=ua1[i],i=1..r)},gam)):
  od:
 mu:=mu+t^k*vu:
od:
lu:=lu*mu:




lu:=taylor(lu,t=0,n*r+1):
lu:=add(coeff(lu,t,i)*t^i,i=0..n*r):
lu:=taylor(lu,y=0,n+1):
[1,seq(expand(i!*coeff(lu,y,i)),i=1..n)]:
end:



#SeqBrnDJ(r,n): 
#Like SeqBrn(r,n,t) but using the generating function
#of J.S.Devitt and D.M. Jackson
#For example, try:
#SeqBrnDJ(2,10);`):
SeqBrnDJ:=proc(r,n) local t:
subs(t=1,SeqHrnDJ(r,n,t)):
end:



#LD(L): loaded die with L
LD:=proc(L) local s,ra,i:
s:=convert(L,`+`):
ra:=rand(1..s)():
for i from 1 to nops(L) do
 if ra<=convert([op(1..i,L)],`+`) then
   RETURN(i):
 fi:
od:
FAIL:
end:

#RandSPnk(n,k): picks a randon set partition
RandSPnk:=proc(n,k) local L,ra,lu:
if n=1 then
 if k=1 then
  RETURN({{1}}):
 else
   RETURN({}):
 fi:
fi:

L:=[stirling2(n-1,k-1),stirling2(n-1,k)$k]:

ra:=LD(L):

if ra=1 then
 lu:=RandSPnk(n-1,k-1):
  RETURN({{n}} union lu):
else
  lu:=RandSPnk(n-1,k):
  ra:=ra-1:
  RETURN({op(1..ra-1,lu),lu[ra] union {n},op(ra+1..k,lu)}):
fi:

end:


#RandSPn(n): a (uniformly) random set-partition of {1,...,n}
#For example, try:
#RandSPn(3);
#add(x[RandSPn(3)],i=500);
RandSPn:=proc(n) local L,k,ra:
L:=[seq(stirling2(n,k),k=1..n)]:
ra:=LD(L):
RandSPnk(n,ra):
end:



#SeqB2n(n): the sequence of length n+1 whose
#(i+1)^th item is the number of ways of breaking-up
#2 identical Russian Dolls, each consisting of i layers
#into smaller Russian Dolls
#SeqB2n(10);
SeqB2n:=proc(n) local gu,mu,z1,z2,P,i:
P:=
f->expand(z2*diff(f,z2)+z1^4*diff(f,z2$2)/2+
z1^3*diff(diff(f,z1),z2)+z1^2*diff(f,z1$2)/2+
z1^3*diff(f,z2)+z1^2*diff(f,z1)+z2*f):

mu:=1:
gu:=[1]:

for i from 1 to n do
 mu:=P(mu):
 gu:=[op(gu),subs({z1=1,z2=1},mu)]:
od:

gu:

end: