# Alex Gentul hw 27 OK to post

info1:=proc() print(`SSe(Se,N)`); end:
info1();

	# Finish Procedure 
	# SSe(Se,N) 
	# that we didn't have a chance to do in class.

	# SSe(Se,N): inputs a set of non-negative integers Se, that must include 0, and a pos. integer k, and outputs
	# the sequence consisting of the number of unlabelled rooted trees with outdegree belonging to Se, with n vertices for n=1..N. For
	# example 
	# SSe({0,1,2,3},10); 
	# should output the same as 
	# S3(10);

SSe:=proc(Se,N) local s,i,P,p,A,sus,f,j,a;

P:=x;
A:=CIPseq(max(Se),s);
print(A);

for i in Se minus {0} do
	P:=P+A[i];
end;

for i from 1 to N do
	sub:={s[1]=f};
		for j from 2 to max(Se) do
			sub:=sub union {s[j]=subs(x=x^j,f};
		end:
	f:=expand(x+x*(subs(sub,P)));
	f:=add(coeff(f,x,a,)*x^a,a=0..N);
end:
[seq(ceoff(f,x,a),a=1..N)];
end:


# Find as many sets Se (that contain 0) for which 
# SSe(Se,20) 
# is in Sloane.
















































####### c27.txt ########
Help:=proc(): print(`Wt(pi,S), CIP(G,s), S2(N), s3(N) `): 
print( ` CIPseq(N,s) `):
end:

with(combinat):
with(group):

###From C26.txt
#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:
###End C26.txt

#T(x)=x+x*(T(x)+T(x)^2/2+T(x^2)/2)

#s2(N): inputs a pos. integer N and outputs
#the sequence consisting of the number of
#unlabelled rooted trees with outdegree in {0,1,2}
#with n vertices for n=1..N
S2:=proc(N) local f,i,brian:
f:=x:

for i from  1 to N do
f:=expand(x+x*(f+f^2/2+ subs(x=x^2,f)/2)):
f:=add(coeff(f,x,brian)*x^brian,brian=0..N):
od:

[seq(coeff(f,x,brian),brian=1..N)]:


end:




# T(x)=x+x*( T(x)+ 1/2*T(x)^2+T(x^2)/2+ 
# (1/6)*T(x)^3+1/2*T(x^2)*T(x)+1/3*T(x^3)



#s3(N): inputs a pos. integer N and outputs
#the sequence consisting of the number of
#unlabelled rooted trees with outdegree in {0,1,2,3}
#with n vertices for n=1..N
S3:=proc(N) local f,i,brian:
f:=x:

for i from  1 to N do
f:=expand(x+x*(f+f^2/2+ subs(x=x^2,f)/2
+f^3/6+1/2*subs(x=x^2,f)*f+1/3*subs(x=x^3,f))
):
f:=add(coeff(f,x,brian)*x^brian,brian=0..N):
od:
print(f);

[seq(coeff(f,x,brian),brian=1..N)]:


end:


#exp(s[1]*t+s[2]*t^2/2+s[3]*t^3/3+...)
#CIPseq(N,s): The Polya cycle-index polynomials
#for S_n for n=1..N using s[1], ... , s[n]
CIPseq:=proc(N,s) local n,t,f:
f:=exp(add(s[n]*t^n/n,n=1..N)):
f:=taylor(f,t=0,N+1):
[seq(expand(coeff(f,t,n)),n=1..N)]:
end: