# Ok to post homework
# Lucy Martinez, 04-27-2025, Assignment 25


######################
#Problem 1: Briefly describe the progress on your team's project
# My team has now written a few procedures for the algorithm
# that we plan to implement. We have met a few times this week





#Problem 3: According to Polya theory (see above essay):
# The number of sequences of length n with entries in {1,...,c}
# up to the action of the symmetric group
#(i.e. the number of equivalence classes)
# is obtained by plugging in x[i]=c for all i=1..n.
# Write a procedure NuEqC(n,c) that uses CIPsn(n,x) to find that number

# What is
# [seq(NuEqC(i,2),i=1..12)]=[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]

# Triangular numbers:
# [seq(NuEqC(i,3),i=1..12)]=[3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91]

# Tetrahedral numbers:
# [seq(NuEqC(i,4),i=1..12)]=[4, 10, 20, 35, 56, 84, 120, 165, 220, 286, 364, 455]

NuEqC:=proc(n,c) local poly,x,i:
poly:=CIPsn(n,x):
subs(seq({x[i]=c},i=1..n),poly):
end:

####################################from previous classes:
#C25.txt: April 24, 2025
Help25:=proc():
print(`EC(pi,i), PtoC(pi),WtP(pi,x),CIP(G,x),Mult(a)`):
print(`Pars1(n,m), CIPsn(n,x)`):
end:

with(combinat):

#You have a set S of objects on n "point", with a natural action
# of a group of permutations (subset of S_n)
# and you want to count the number of isomorphism classes
# for each s in S, Images(s)={g(s), g in G},
# you want to count how many such isomorphism classes
# Def: For a permutation pi decompose into cycles 
# pi=C1*C2*...*Ck, let X be a formal variable
# Wt(pi)=X[nops(C1)]*X[nops(C2)]*...*X[nops(Ck)]
# Cycle-index polynomial of a group of permutations G (a subset of S_n)
# is Sum(Wt(pi), pi in G)/nops(G)


#EC(pi,i): inputs a permutation pi and finds the cycle that belongs
# to i
EC:=proc(pi,i) local a,C,m,j:
C:=[i]:
a:=pi[i]:
while a<>i do
 C:=[op(C),a]:
 a:=pi[a]:
od:
m:=min(C):
for j from 1 to nops(C) while C[j]<>m do od:
[op(j..nops(C),C),op(1..j-1,C)]:
end:


#PtoC(pi): inputs a permutation in one-line notation, outputs
# its cycle decomposition
PtoC:=proc(pi) local S,C,c:
C:={}:
S:={op(pi)}:
while S<>{} do
 c:=EC(pi,S[1]):
 S:=S minus {op(c)}:
 C:=C union {c}:
od:
C:
end:


#WtP(pi,x): the weight of the permutation pi according to
#x[1]^NumberOfOneCycles*x[2]^NumberOfTwoCycles*...*x[n]^NumberOfnCycles
WtP:=proc(pi,x) local C,c:
C:=PtoC(pi):
mul(x[nops(c)],c in C):
end:


#CIP(G,x): the Polya cycle-index polynomial of the group G of
# permutations
CIP:=proc(G,x) local g:
add(WtP(g,x),g in G)/nops(G):
end:

#For the symmetric group of permutations:
#Every monomial x1^a1*...xn^a correspond to
# 1*a1+2*a2+...+n*an=n al add-up to n*an
# integer patition of n in FREQUENCY notation

# Why is it that if there are a1 1-cycle, a2 2-cycles, ... ,
# an n-cycles then 1*a1+2*a2+...+n*an=n*an
# Because you are partitioning the permutation

# What is the coefficient of x1^a1*x2^a2*...*xn^an in
# CIP(permute(n),x)
# How many ways can you fill-in:
# -,-,-,... a1 times
# [-,-],[-,-],... a2 times
# and so on
# Example: suppose we have p:=[3,0,2,1] where
# there are p[1] 1-cycle
# there are p[2] 2-cycles
# there are p[3] 3-cycles
# there are p[4] 4-cycles
# [[],[],[],[,,],[,,],[,,,]]
# [1^3, 1^0, 3^2, 4^1]
# n!/(1!^a1*2!^a2*...*n!^an)/(a1!*a2!*...*an!)*0!^a1*1!^a2*...*(n-1)!^an


#Mult(a): inputs an integer partition a in frequency-notation
# 1^a1...n^an [a1,...,an] of non-negative integers
# n=nops(a)
Mult:=proc(a) local n,i:
n:=nops(a):
if n<>add(a[i]*i,i=1..n) then
 return(FAIL):
fi:
n!/(mul(i!^a[i],i=1..n))/(mul(a[i]!,i=1..n))*(mul((i-1)!^a[i],i=1..n)):
#n!/(1!^a1*2!^a2*...*n!^an)/(a1!*a2!*...*an!)*0!^a1*1!^a2*...*(n-1)!^an
end:


#Added after class

#Pars1(n,m): The set of [a1,a2,...,an] such that 1*a1+...+n*an=m
Pars1:=proc(n,m) local S,s1,i,S1:
option remember:
if m=0 then
 RETURN({[0$n]}):
fi:

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

S:={}:

for i from 0 to trunc(m/n) do
S1:=Pars1(n-1,m-i*n):
S:=S union {seq([op(s1),i], s1 in S1)}:
od:
S:
end:

#CIPsn(n,x):CIP(permute(n),x), done cleverly, using Mult(a):
CIPsn:=proc(n,x) local S,s,i:
S:=Pars1(n,n):
add(Mult(s)*mul(x[i]^s[i],i=1..n), s in S)/n!:
end: