####################################################################################3
# Making sure I understand how basic plethysm works in Sage
####################################################################################3

# p3[p2] = p6
print("p[3].plethysm(p[2]) should equal p6: %s" % (p[3].plethysm(p[2])))

# p1 should be the identity
print("s[3,2].plethysm(p[1]) should be s[3,2]: %s" % (s[3,2].plethysm(p[1])))
print("p[1].plethysm(s[3,2]) should be s[3,2]: %s" % (s(p[1].plethysm(s[3,2]))))

# pn should be in the center
print("p[3].plethysm(s[3,2])-s[3,2].plethysm(p[3]) should be 0: %s" % (p[3].plethysm(s[3,2])-s[3,2].plethysm(p[3])))

# See Example 6 in Expose
# q, t, z should be all be treated as generators of the K-algebra; pm[q] = q**m, pm[t] = t**m, pm[z] = z**m
print("p3[q*p3] = %s, p4[t*p2] = %s, p5[z*p4] = %s"%\
                  (p[3].plethysm(q*p[3]), p[4].plethysm(t*p[2]), p[5].plethysm(z*p[4])))
                  
# We can instruct Sage to leave certain generators fixed if we wish
# The default is to treat them all as degree-one generators that are not fixed by plethysm with p_n
print(p[3].plethysm(t*q*z*p[0]))
print(p[3].plethysm(t*q*z*p[0],exclude=[z]))
print(p[3].plethysm(t*q*z*p[0],exclude=[z,q,t]))

####################################################################################3
# Code for checking that we understand the various negatives
####################################################################################3

# All of these should be negative
print("p[3].plethysm(-p[3]): %s" % (p[3].plethysm(-p[4])))
print("p[3].plethysm(-p[4]): %s" % (p[3].plethysm(-p[4])))
print("p[4].plethysm(-p[3]): %s" % (p[3].plethysm(-p[4])))
print("p[4].plethysm(-p[4]): %s" % (p[3].plethysm(-p[4])))

# Theorem 8 of Expose: g[-A] = (-1)**n*omega(g)[A] if g homogeneous of degree n
ans1 = s[3,1].plethysm(-p[2])
ans2 = (-1)**4*s[3,1].omega().plethysm(p[2])
print("Difference in sides of formula for s[3,1].plethysm(-p[2]): %s" % (ans1-ans2))

ans1 = s[3,2].plethysm(-p[3])
ans2 = (-1)**5*s[3,2].omega().plethysm(p[3])
print("Difference in sides of formula for s[3,2].plethysm(-p[3]): %s" % (ans1-ans2))

ans1 = h[2,1,1].plethysm(-p[2])
ans2 = (-1)**4*h[2,1,1].omega().plethysm(p[2])
print("Difference in sides of formula for h[2,1,1].plethysm(-p[2]): %s" % (ans1-ans2))

# Note that f[epsilon X] = - omega(f)[X]
# In particular, pulling epsilon out front of pn[X] plethysm leads to (-1)**n in front
# This is not true for pulling out from other symmetric functions, as third example shows
ans1 = -p[3].omega().plethysm(p[2])
ans2 = (-1)**3*p[3].plethysm(p[2])
print("Diff: %s" % (ans1-ans2))

ans1 = -p[4].omega().plethysm(p[2])
ans2 = (-1)**4*p[4].plethysm(p[2])
print("Diff: %s" % (ans1-ans2))

ans1 = -s[3,1].omega().plethysm(p[2])
ans2 = (-1)**4*s[3,1].plethysm(p[2])
print("Diff: %s" % (ans1-ans2))