#John Kim
#Use whatever you like

#read "C:/Users/John Y. Kim/Documents/Maple/hw26.txt":



Help:=proc(): print(` A81(N), A55(N), Wt(pi,s) ,CIP(G,s), NuCol(G,c) `):
print(`Cn(n) , WtCol(G,c,x) `): 
end:
with(combinat):
with(group):

###From C23.txt
#Counting unlabelled rooted trees
#A[n]=number of rooted unalabelled trees with n vertices
#A[1]=1, A[2]=1, A[3]=2,
#Removing the root, leads to a MULTISET of smaller rooted trees
#HOW MANY with one vertex? 1/(1-x)^A[1]
#HOW MANY with two vertices 1/(1-x^2)^A[2]
#...
#How MANY with i vertices? 1/(1-x^i)^A[i]
#A[1]*x+A[2]*x^2+A[3]*x^3+A[4]*x^4+...=
#x/(1-x)^A[1]/(1-x^2)^A[2]/(1-x^3)^A[3]/.....

#A81(N): the first N terms of the enumerating 
#sequence for the  numnber of UNLABELLED rooted trees
#with n vertices
A81:=proc(N) local A,f,n,x,i:

A:=[1]:

for n from 2 to N do
f:=x*mul(1/(1-x^i)^A[i],i=1..nops(A)):
f:=taylor(f,x=0,n+1):
A:=[op(A),coeff(f,x,n)]:
od:

A:

end:

###end from C23.txt



#Let R(x) be the (ORDINARY) generating function for
#Unlabelled rooted trees (counted by URT(N))
#R(x)=add(URT(infinity)[i]*x^i,i=0..infinity);
#Let A(x) be the generating function for unlabelled trees
#A(x)=R(x)-R(x)^2/2+R(x^2)/2

#A55(N): the first N terms in OEIS A000055, Neil Sloane's
#favorite sequence (default, viewed April 23, 2012 and many
#years before, possibly from the start)
#the number of unlabelled trees
A55:=proc(N) local R,A,L,x,i:

L:=A81(N):

R:=add(L[i]*x^i,i=1..N):

A:=expand(R-R^2/2+subs(x=x^2,R)/2):

[seq(coeff(A,x,i), i=1..N)]:

end:



#Wt(pi,s): the weight (according to Polya) of
#the permutation pi
#For example : Wt([1,2,3],s)=s[1]^3
Wt:=proc(pi,s) local L,n,i,y,W:
n:=nops(pi):
L:=convert(pi, `disjcyc`);

W:=mul( s[ nops(L[i])], i=1..nops(L)):

W*s[1]^n/subs({seq(s[i]=s[1]^i,i=2..n)},W):

end:

#CIP(G,s): inputs a group (or even a set) of permutations
#and outputs the Polya Cycle-index polynomial
#1/|G| *( Sum of the weights of all pi in G)
#weight(pi)=prod(s[i]^(#of cycles of length i)
#s is a letter
#For example,
#wt([2,1,4,5,3])=s[2]^1*s[3]^1=s[2]*s[3]
#wt([1,2,3,4,5,6])=s[1]^6
#wt([2,3,4,5,6,1])=s[6]
CIP:=proc(G,s) local pi:
add(Wt(pi,s), pi in G)/nops(G):
end:

#NuCol(G,c): the number of ways of coloring a combinatorial
#object (up to equivalence) with the group being G
#with c colors
#For example
#NuCol(permute(4),c);
NuCol:=proc(G,c) local s,P,n:

n:=nops(G[1]):

P:=CIP(G,s):

factor(subs( { seq( s[i]=c, i=1..n)}, P)):

end:


#Cn(n): the group of rotations of a circle with n beads
Cn:=proc(n) local i,j:

#[i,i+1,...,n,1,..., i-1]

{ seq( [seq(j,j=i..n), seq(j,j=1..i-1)], i=1..n)}:

end:

##WtCol(G,c,x): The generating function of all colorings
#with the weigt
#x[1]^NumberOFBeadsColoredColor1*
#x[2]^NumberOFBeadsColoredColor2* ....

WtCol:=proc(G,c,x) local i,n,P:


n:=nops(G[1]):

P:=CIP(G,s):


#s[1]=x[1]+...+x[c]
#s[2]=x[1]^2+ ...+x[c]^2
#....

expand(subs({seq( s[i]=add(x[j]^i,j=1..c),i=1..n)},P)):

end:




#Mul(P,Q) that inputs two permutations (written as lists) of the same lengths and outputs their product. For example,
#Mul([2,3,1,4],[2,1,4,3]);
#should output
#[1,4,2,3] 

Mul:=proc(P,Q) local i:

[seq(Q[P[i]],i=1..nops(P))]:

end:



#Gp(S) that inputs a set of permutations of the same length 
#(output FAIL if S is not a set or it is not a set of lists 
#of the same length, etc.) and outputs the set of all 
#permutations that belong to the group generated by the 
#members of S. (Don't forget the identity permutation!). 
#For example
#Gp({[2,1]});
#should output
#{[1,2],[2,1]}.
#Gp({[2,1,3,4],[1,2,4,3]});
#should output
#{[1,2,3,4],[2,1,3,4],[1,2,4,3],[2,1,4,3]}   . 

Gp:=proc(S) local G1,G,pi,g:

G1:={}:
G:=S:

while G1<>G do
 G1:=G:
 for pi in S do
  for g in G do
   G:=G union {Mul(g,pi)}:
  od:
 od:
od:

G:

end:




#Find the subgroup of S6 generated by the three rotations 
#(about the x,y, and z axes) of a cube. Using Nathaniel Shar's 
#suggestion, let's stick to the usual convention:
#Front: 1 ; Back: 6; Left:3; Right: 4; Bottom:2; Top=5. 

#Gp({[1,4,2,5,3,6],[5,1,3,4,6,2],[4,2,1,6,5,3]});








#How many necklaces are there with 15 beads colored Red, 
#23 beads colored Blue, and 12 beads colored Green? 

#coeff(coeff(coeff(WtCol(Cn(50),3,x),x[1]^15),x[2]^23),x[3]^12);
#                      37564175809042384320


#In how many ways can you color the faces of a cube 
#(up to rotational symmetry) with 3 faces colored Red, 
#2 faces colored Blue, and one face colored Green? 

#WtCol(Gp({[1,4,2,5,3,6],[5,1,3,4,6,2],[4,2,1,6,5,3]}),3,x);
#coeff(coeff(coeff(%,x[1]^3),x[2]^2),x[3]);
#                      3