Lie-algebra representation powers - plethysms

lpetrich

I've done some more work on my SemisimpleLieAlgebras.zip package, adding decompositions of representation powers. It was only recently that I discovered that they are called "plethysms".

Rep powers can be decomposed by symmetry types, a feature that can be important in their applications. For power p, the types are symmetric (p >=1), antisymmetric (p >= 2), and mixed (p >= 3), with more than one mixed type for p >= 4.

These are related to the irreps in (SU(n) fundamental rep)^p, and they share multiplicities. The symmetry types and irreps can be denoted with Young diagrams having p boxes.

With a lot of experimentation and trial and error, I have worked out how to do rep-power decompositions, but the algorithm is hideously complicated, Though I've successfully tested it for low powers, I don't have a rigorous proof for it, and I've been unsuccessful in finding an algorithm in the professional literature. I've also had to create versions for reducible reps and algebra-product reps, because rep-power decompositions do not have any simple decompositions of those in terms of the rpd's of the component irreps of the component algebras.

-

As an example, consider flavor and spin symmetries of the three light quarks: SU(3) and SU(2) respectively. They can be combined into SU(6), and one can undo this combination into SU(3)*SU(2) again with

Mathematica: SubalgMultSU[ldqksp,{3,2}]; res = DoBranching[ldqksp,{1,0,0,0,0}]
Python: ldqksp = SubalgMultSU((3,2)); res = ldqksp.DoBranching((1,0,0,0,0))

Breakdown of cube of fundamental rep of SU(6):

Mathematica: MakeLieAlgebra[laqksp,{1,5}]; res = DecomposeRepPower[laqksp,{1,0,0,0,0},3]
Python: res = DecomposeRepPower((1,5), (1,0,0,0,0), 3)

Count, Irrep:
S: (1, (3, 0, 0, 0, 0))
M: (1, (1, 1, 0, 0, 0))
A: (1, (0, 0, 1, 0, 0))

Breakdown of cube of corresponding rep of SU(3)*SU(2):

Mathematica: MakeLieAlgebra[laqk, {1, 2}]; MakeLieAlgebra[lasp, {1, 1}]; res = DecomposeAlgProdRepPower[{laqk, lasp}, {{1, 0}, {1}}, 3]
Python: res = DecomposeAlgProdRepPower(((1,2),(1,1)), ((1,0),(1,)), 3)

Count, SU(3) Irrep, SU(2) Irrep
S: (1, ((3, 0), (3,))) . (1, ((1, 1), (1,)))
M: (1, ((1, 1), (3,))) . (1, ((3, 0), (1,))) . (1, ((1, 1), (1,))) . (1, ((0, 0), (1,)))
A: (1, ((1, 1), (1,))) . (1, ((0, 0), (3,)))

-

I'll now translate these results into typical physicist language, using rep multiplicities and spins:

Each quark: (3,1/2)
S: (10,3/2) + (8,1/2)
M: (8,3/2) + (10,1/2) + (8,1/2) + (1,1/2)
A: (8,1/2) + (1,3/2)

The symmetric one is the one that matches the light-baryon spectrum, and not the mixed or the antisymmetric ones. This is contrary to Fermi-Dirac statistics, and it took QCD to resolve that discrepancy. The quarks are antisymmetric in color, which gives them overall antisymmetry, thus preserving FD statistics.

Simon_Tyler

Hi Ipetrich,

Glad to see you're still working on this - it's looking good (although it's getting beyond the stuff that I've played with). I like how you're maintaining the Mathematica and python code and always giving examples from QFT.

Branching rules from plethysms are implemented in sage.
It's described here (with references):
http://www.sagemath.org/doc/thematic...html#plethysms
and the code can be found here:
http://hg.sagemath.org/sage-main/fil...cters.py#l2048

This should provide some good comparisons/test cases for your code.

lpetrich

Thanx. I checked on that, and it does not seem to list plethysm decompositions by symmetry type.

That was part of the algorithms for branching rules, and I'd implemented all of them:

• Levi -- root demotion. Sage's code does not seem to calculate the resulting U(1) factors, however.
• Dynkin-diagram extension splitting.
• Division of orthogonal and symplectic groups: SO(m+n)->SO(m)*SO(n), Sp(2(m+n)) -> Sp(2m)*Sp(2n). All but SO(even) -> SO(odd)*SO(odd) can be derived from root demotion and extension splitting, however, using SO(2) ~ U(1).
• Symmetric subgroups: SU(n) -> SO(n), SU(2n) -> Sp(2n), SO(2n) -> SO(2n-1), E(6) -> F(4)
• Tensor products: SU(m*n) -> SU(m)*SU(n), SO/Sp(m*n) -> SO/Sp(m)*SO/Sp(n) (SO -> SO*SO or Sp*Sp, Sp -> SO*Sp or Sp*SO)
• Symmetric powers: I implement what seems like a superset of their special case: anything -> SU(2). Root height -> SU(2) root value.

The plethysm rule seemed like a subcase of one of Dynkin's theorems. I dug up a translation of a collection of his papers, and I found various useful theorems.

I've also implemented all the maximal subalgebras for the exceptional ones; Slansky lists all the maximal subalgebras for algebras with rank <= 8.


These special cases do not encompass all the subalgebras of classical algebras, but Dynkin has some theorems that supply constraints on the remaining ones:

* At least one "simple" irrep must map onto a single term in the subalgebra, but that term need not have a "simple" irrep. Here are the "simple" ones:

A(n) / SU(n-1): vector and conjugate vector (dim 2n+1)
B(n) / SO(2n+1): vector (dim 2n+1) and spinor (dim 2^n)
C(n) / Sp(2n): vector (dim 2n)
D(n) / SO(2n): vector (dim 2n) and the two spinors (dim 2^(n-1))

* The algebra's adjoint rep must map onto subalgebra's adjoint rep, though usually with extra irreps. The adjoint can be constructed with the help of powers and symmetry types (plethysms) of the defining irreps for the original and the subalgebra.

* The original algebra's irrep and the subalgebra's one must have the same "reality":

• Real -- self-conjugate with even height -- adjoint in antisymmetric part of square
• Pseudoreal or quaternionic -- self-conjugate with odd height -- adjoint in symmetric part of square
• Complex -- not self-conjugate -- adjoint in higher power

In general, (irrep)*(its conjugate) always contains the adjoint rep.

I could not find any general algorithm for calculating projection matrices in this case, or any subcases, so I left unimplemented most of these ones that Slansky had listed. Some examples of this case:

SU(6) -> SU(3): (1,0,0,0,0) 6 -> (2,0) 6
SO(7) -> G(2): (1,0,0) 7 -> (0,1) 7

lpetrich

I'd found branching-rule projection matrices in ad hoc fashion; part of the problem is finding which subalgebra rep root an original-algebra rep root maps onto.

Turning to Sage's D7 -> G2 plethysm example, I note that D7's vector rep (1,0,0,0,0,0,0) maps onto G2's adjoint rep (1,0), both with dimension 14. For G2, Sage uses reverse order from my code.

In the example's final output, the first five entries were all antisymmetric products, and they agree with what my code produces. I tried my code on D7 also, and the antisymmetric products were all single irreps with one weight value 1 and the others zero. Going to the 6th power yielded (0,0,0,0,0,1,1) -- the product of the two spinors. So mapping an individual spinor onto G2 will be tricky.


In my code, I identify the symmetry types with Young diagrams, and these map directly onto the irreps in powers of the fundamental rep of SU(n) / A(n+1).
Power p symmetric: (p)
Power p antisymmetric: (p 1's)
Power p mixed: (other values that add up to p)

Symmetric: all index permutations equally weighted
Antisymmetric: like symmetric, but all odd permutations have reverse sign
Mixed: hard to find any simple formulas


Getting back to Young diagrams, one can easily find SU(n) rep dimensions with them, and there's a moderately complicated algorithm called the Littlewood-Richardson rule for rep products. Are there similar techniques for SO(n) and Sp(2n)? For SO(n), product irreps of the vector rep have to be traceless, and one has to take that into account. Like:

SU(n) -- vector * vector = (symmetric tensor) + (antisymmetric tensor)
SO(n) -- the symmetric tensor becomes (symmetric traceless tensor) + (scalar)
Sp(2n) -- the antisymmetric tensor becomes (antisymmetric sort-of-traceless tensor) + (scalar)

I say sort-of-traceless, because it's an operation that involves subtracting out the symplectic form for the algebra:
J = {{0, I}, {-I, 0}}
in the fashion of tracelessness, subtracting out the identity matrix.

lpetrich

I have updated my Lie-algebra code yet again, this time to support U(1) factors in the algebra-list versions of my rep powers. U(1) factors can be any integers or rational numbers; they follow the max-weight vectors. It's possible to have only U(1) factors if one so desires. Here, I'll give some results in Python form - it should be easy to translate ino Mma form.

DecomposeAlgProdRepPower( ((1,1),), ((1,),42), 2)
yields
((1, ((2,), 84)),)
((1, ((0,), 84)),)
(adding two spin-1/2's to make a spin 1 (symmetric) and a spin 0 (antisymmetric))

Now to raise vectors in various dimensions in various powers. The outputs are the symmetry-type decomposition lists, whose members have format (count, rep).

2D vector: SO(2) = U(1) -- helicity. Demo of using no semisimple algebras, just U(1)'s
DecomposeAlgProdRepListPower,( (), ((1,),(-1,)), p)
p = 1
((1, (1,)), (1, (-1,)))
p = 2
((1, (2,)), (1, (0,)), (1, (-2,)))
((1, (0,)),)
p = 3
((1, (3,)), (1, (1,)), (1, (-1,)), (1, (-3,)))
((1, (1,)), (1, (-1,)))
()
p = 4
((1, (4,)), (1, (2,)), (1, (0,)), (1, (-2,)), (1, (-4,)))
((1, (2,)), (1, (0,)), (1, (-2,)))
((1, (0,)),)
()
()

3D vector: SO(3) = B(1)
DecomposeRepPower( (2,1), (2,), p)
p = 1
((1, (2,)),)
p = 2
((1, (4,)), (1, (0,)))
((1, (2,)),)
p = 3
((1, (6,)), (1, (2,)))
((1, (4,)), (1, (2,)))
((1, (0,)),)
p = 4
((1, (8,)), (1, (4,)), (1, (0,)))
((1, (6,)), (1, (4,)), (1, (2,)))
((1, (4,)), (1, (0,)))
((1, (2,)),)
()

4D vector: SO(4) = D(2)
DecomposeRepPower( (4,2), (1,1), p)
p = 1
((1, (1, 1)),)
p = 2
((1, (2, 2)), (1, (0, 0)))
((1, (0, 2)), (1, (2, 0)))
p = 3
((1, (3, 3)), (1, (1, 1)))
((1, (1, 3)), (1, (3, 1)), (1, (1, 1)))
((1, (1, 1)),)
p = 4
((1, (4, 4)), (1, (2, 2)), (1, (0, 0)))
((1, (2, 4)), (1, (4, 2)), (1, (2, 2)), (1, (2, 0)), (1, (0, 2)))
((1, (2, 2)), (1, (4, 0)), (1, (0, 4)), (1, (0, 0)))
((1, (2, 2)), (1, (2, 0)), (1, (0, 2)))
((1, (0, 0)),)

The second one for p = 2 is the antisymmetric 2-tensor.
For 2D; scalar: 0
For 3D: vector: (2,)
For 4D: self-dual and anti-self-dual parts: (0, 2) and (2, 0)
Duality: (1/2)*e.A = +/- A

For the Faraday tensor, the self-dual and anti-self-dual parts are, in 3+1 notation, E+i*B and E-i*B

The third one for p = 4 has the symmetries of the Riemann tensor. Its Young diagram is (2,2)
For 2D, it's given by the Ricci scalar: 0
For 3D, it's given by the Ricci tensor (symmetric 2-tensor): (4,) and (0,)
For 4D it's given by the Ricci tensor + the Weyl tensor:
(2, 2) and (0, 0)
(4, 0) and (0, 4)
The latter has self-dual and anti-self-dual parts.

lpetrich

I've fixed a bug in the Python version, SemisimpleLieAlgebras.py -- DecomposeRepProduct had been returning empty products, but it's now returning correct ones. I've also added counted-list inputs to rep powers, those being lists of (count,maxwts) instead of just maxwts.

I also have some more conjectures in addition to my conjectured rep-power-decomposition (plethysm) formula. I've looked in the literature, but I can't find whether they've been anticipated.

A rep-product one:

Take representation R and some power, and take the part with the symmetry type denoted by Young diagram Y. I will call this part R^Y (R to power Y).

The rep product R^Y1 * R^Y2 = sum_Y c_Y1,Y2,Y R^Y

where c_Y1,Y2,Y are the Littlewood-Richardson coefficients for how many copies of Y are in the product of Y1 and Y2.

This is correct if R is the fundamental representation of SU(n), and I've verified it in some other cases. I've looked in the literature for a general proof, without success.

A rep-sum one:

For representations R1 and R2 of a Lie algebra, the sum rep R1 + R2 satisfies
(R1 + R2)^Y = sum_Y1,Y2 c_Y1,Y2,Y R1^Y1 * R2^Y2

I've verified that one also in some cases. I'm in the design phase of implementing it more generally, and I should be able to implement it before long.

Testing cannot be complete, of course, but it should be feasible for the lower-rank algebras and lower-multiplicity reps.

Others?

I've found counterexamples to (R^Y1)^Y2 = (R^Y2)^Y1, and I've been unable to guess what cx is in

(R1,R2)^Y = sum_Y1,Y2 cx_Y1,Y2,Y (R1^Y1,R2^Y2)

other than Y1 and Y2 having the same order as Y for cx to be nonzero.

lpetrich

I've updated SemisimpleLieAlgebras.zip to include code for testing my previous post's conjectures in the Mathematica version. I've run it on several examples, and it has been successful so far. I still haven't found proofs, however.

But both Mathematica and Python versions now have code for doing purely-symmetric or purely-antisymmetric representation-power decompositions (plethysms).

Example in Python:
The Sage results that Simon_Tyler had linked to I could duplicate with (Mma)

MakeLieAlgebra[lag2, {7,2}]
DecomposeRepPwrSym[lag2, {1,0}, pwr, -1]

with pwr = 1,2,3,4,5
Sage's code uses the reverse order of roots/weights from mine, which is Robert Cahn's order.

lpetrich

I've added some stuff on Schur functions and related ones to my Young-diagram Mma notebook in that archive. I've done so because Schur functions are related to the irreps of SU(n) / A(n-1) / SL(n,R). They are even specified in the same way, with partitions / Young diagrams. In general:

Schur function: S(part / YD, vector)

The smallest ones:
S({},{x1,x2,x3}) = 1
S({1},{x1,x2,x3}) = x1+x2+x3
S({2},{x1,x2,x3}) = x1^2+x1*x2+x2^2+x1*x3+x2*x3+x3^2
S({1,1},{x1,x2,x3}) = x1*x2+x1*x3+x2*x3

They form a basis for the symmetric polynomials in their vector arg.

Products of them give linear combinations of them specified with the Littlewood-Richardson rule, familiar from products of SU(n) irreps.

Plethysms on Schur functions are defined with this algorithm:

To calculate S(a)[S(b)] the plethysm of S(a) over S(b), do
For a vector x, decompose S(b,x) as a sum of monic monomials in products and powers of the x components: x(1)^2*x(2), etc.
Make them into a vector, and evaluate S(a,(this monic-monomial vector))

You can now decompose this value into a sum of S(c,x) values. In short,
S(a)[S(b)] = sum over c of pl(a,b,c) S(c)

One can replace S(a) with some other sort of symmetric function on a vector, and still get (Schur function) -> (combination of Schur functions).

However, it's been difficult for me to find an algorithm for calculating pl(a,b,c) comparable to the Littlewood-Richardson rule. One can do it with the algorithm here, or with A(n) Lie-algebra irreps, but I wish to avoid these algorithms, because they can require a very large amount of calculation.

lpetrich

I've succeeded in finding an algorithm for pl(a,b,c) when a and b are partitions / Young diagrams for symmetric or antisymmetric SU(n) reps. That is, {n} or {n 1's} / {1^n}

This algorithm sidesteps calculating reps or working with Schur-function vector args, thus making it much faster. It's essentially the algorithm in EUROCAL '85. European Conference on ... - Bob F. Caviness - Google Books

Some special cases of Schur functions:
S({n},x) = homogeneous symmetric polynomial in x of order n
Terms have the form (x1^p1)*(x2^p2)*... powers summing to n
S({1^n},x) = elementary symmetric polynomial in x of order n
Terms have the form x1*x2*... powers no greater than 1 and summing to n
All coefficients are 1

Calculate S({n},x^m) and S({1^n},x^m} from power sums of x using (power-sum to homogeneous) and (power-sum to elementary) formulas.

Calculate the plethysm by treating the above S's as power sums in the (power-sum to homogeneous) and (power-sum to elementary) formulas.

Along the way, I found a formula for power sums in terms of Schur functions:
sum(x^n) = sum over k: 0 to n-1 of (-1)^k * S({n-k, k 1's},x)

I was unable to prove it, though I was able to verify it for small n's.

lpetrich

Various identities involving plethysms: rep-power decompositions.

http://www.sbfisica.org.br/bjp/files/v32_641.pdf
The Plethysm Technique Applied to the Classification of Nuclear States, by Alcaras, Tambergs, Krasta, Ruza, Katkevicius

plt(a,x) is the plethysm of a on x, with a the plethysmer and x the plethysmed one, an irrep. a is a partition / Young diagram. A power p decomposes into diagrams a with length p.

Inner product: a.b = sum of ipc(a,b,c) * c

ipc(a,b,c) = 1/N sum over diagrams r of N(r)*X(a,r)*X(b,r)*X(c,r)
X(a,r) = characters of symmetric group for irrep a and class b (always integer), N(r) is order of class with part/diag r, N is total order.

Outer product: a*b = sum of opc(a,b,c) * c

opc(a,b,c) one calculates uses the Littlewood-Richardson rule.

Reps x,y have sum (sum of coefficients of irreps in it), difference (likewise with differences), product (rep product)

The identities:

plt(a +- b, x) = plt(a,x) +- plt(b,x)
plt(a, x + y) = sum over b,c of opc(b,c,a) * plt(b,x)*plt(c,y)
plt(a, x - y) = sum over b,c of opc(b,c,a) * (-1)^length(c) * plt(b,x)*plt(c,y)
plt(a*b,x) = plt(a,x)*plt(b,x)
plt(a,x*y) = sum over b,c of ipc(b,c,a) * plt(b,x)*plt(c,y)
conjg(plt(a,x)) = plt(a,conjg(x))

I earlier went through a lot of trouble to find formulas for plt(a,x+y) and plt(a,x*y), and while I successfully guessed the former one by extrapolating from small examples, I had much more difficulty with the latter one.

lpetrich

I now have plethysms / rep-power decompositions working on Young diagrams, without calculating Lie-algebra reps or Schur functions of vectors.

SelYDPlethysm with args plethysmer selector, order, plethysmed one selector, order
Selectors:
sym = {order}, sym all = {p} up to {order}
ats = {1^order}, ats all = {1^p} up to {1^order}
gen = (all YD's with size = order), gen all = (all YD's with size <= order)

As far as I can tell, it works.

lpetrich

conjg(plt(a,x)) = plt(a,conjg(x)) is correct for conjugate representations, but for conjugate Schur functions / Young diagrams, one finds this:

conjg(plt(a,x)) =
plt(a,conjg(x)) if len(x) is even
plt(conjg(a),conjg(x)) if len(x) is odd

len(x) is the length or order or total power of x.

Another identity: plt(a,plt(b,x)) = plt(plt(a,b),x)

Some special cases: symmetric and antisymmetric, {n} and {1^n}. I call them that because of the corresponding SU(N) / SL(N) irreps and how they are composed from the fundamental one: the symmetrized or the antisymmetrized product.

plt({n},x+y) = sum over k of plt({n-k},x)*plt({k},y)
plt({1^n},x+y) = sum over k of plt({1^(n-k)},x)*plt({1^k},y)

plt({n},x*y) = sum over a of plt(a,x)*plt(a,y) where len(a) = n
plt({1^n},x*y) = sum over a of plt(a,x)*plt(conjg(a),y) where len(a) = n

lpetrich

I think that I ought to give the algorithm that I conjecture for plethysms of Lie-algebra representations. I start with a rep given by a list of sets (n,v) where n is the multiplicity/degeneracy and v is a root vector. The root vectors are all distinct in the list.

For power p, the results are lists of multiplicity-vector sets (N(n_i,1,n_i,2,...n_i,p), v_i,1 + v_i,2 + ... + v_i,p)
I then find irreps from it in standard fashion, by repeatedly looking for the highest weight and subtracting out the irrep with that highest weight.


The problem is, of course, finding the N functions. For the pure symmetric and antisymmetric cases, it is easy. First, find the multiplicities m of each n value, a list of
(m_k,n_k)

and then one finds
N = product over the (m,n)'s of n*(n+s)*(n+2*s)...(n+(m-1)*s)/m!
where s = 1 for symmetric and -1 for antisymmetric. It's rather easy to find this result from symmetry.

But the mixed-symmetry cases are much more difficult. They are plethysms of Young diagrams y on the rep, and I calculate them for all the y's that add up to the power. The y's correspond to irreps of SU(n) that one finds from powers of SU(n)'s fundamental rep, and a plethysm of y on the fundamental rep yields the rep with Young-diagram y.

That's rather easy to show from the correspondence between Schur functions and the characters of SU(n) Lie algebras. In general, a Lie-algebra character for a rep is the sum of n*exp(c.v) for vector c and rep pairs (n,v).

So I calculate all the N(y,{n}) for Young diagrams y and vectors of n: {n}.

-

I first calculate NW(yw,{n}) the number for the SU(n) Weyl orbit with max weight yw, then find
N(y,{n}) = sum over yw of Kostka(y,yw)*NW(yw,{n})

where

(SU(n) irrep for y) = sum over yw of Kostka(y,yw) * (SU(n) Weyl orbit for yw)

or in Schur functions,

S(y,x) = sum over yw of Kostka(y,yw) * SW(yw,x)

where SW(yw,x) is a Schur-Weyl function, the sum over all terms x_i1^r1*x_i2^r2*...*x_i(n)^r(n) for yw = {r1,r2,...,r(n)} and selections x_i's from the vector x.

The Kostka matrix I find using the algorithm in "Determinantal Expression and Recursion for Jack Polynomials", by L. Lapointe, A. Lascoux, J. Morse. Are there any recursive algorithms for this matrix? Algorithms that calculate each power's Kostka matrix from the previous ones.

-

NW(yw,{n}) turns out to be relatively simple. From {n}, find {(m,n)}, and sort the list to place the m's in non-increasing order. The m's form the Young diagram ym.

Find all nestings of yw in ym. Row lengths: yw = {rw}, ym = {rm}. These are all possibilities of
rm(i) = sum over k of rw(j(i,k))
where each rw value appears once and only once in the sums and nestings with different arrangements of repeated values are counted as a single nesting.

For each nesting, find the product over the rm(i)'s in ym of:
using power = count of rw(j(i,k))'s
using repeat counts = counts of repeats of rw(j(i,k)) values.
power!/(product of count! for count in repeat counts) * n(i)*(n(i)-1)*...*(n(i)-power+1)/power!

Add up over all nestings. No nestings possible yields zero. That gives NW(yw,{n}). Find for all yw's. Multiply by the Kostka matrix for the overall power to get the N(y,{n})'s.

This algorithm essentially calculates the plethysm of SU(n) Weyl orbits on the rep character, then uses the appropriate Kostka matrix to find the plethysm of SU(n) irreps on the rep character.

This suggests a proof: turn the rep character into a sum of monomials exp(c.v(i)) where each such term is repeated n(i) times. Then make it a vector of exp(c.v(i)) values and use it as the vector arg in the Schur or Schur-Weyl function. I'd have to work out a lot of combinatorics, but it seems like it could be done in the Schur-Weyl case.

Simon_Tyler

Hey lpetrich,

I don't have much to say about your recent posts and work, since it's out of my range of experience and I don't have time to study it at the moment.
So I'll just say that it sounds good and keep going!

lpetrich

Thanx.

I worked out the combinatorics for the Schur-Weyl case. Let's consider the case where x1, x2, ..., x(n) have the same root value v, the same value of exp(c.v). Let's also consider them divided among the Schur-Weyl row lengths r1, r2, ..., r(p): x1^r1*x2^r2*...*x(p)^r(p) and all permutations. There are thus n*(n-1)*...*(n-p+1) such terms. However, sets of equal r's have the permutations of x's in them counted only once, which means dividing the number of terms by the product (number of r equal to some value)! for all values that the r's equal.

If the full set of terms contains something like x1*x2*x3 + x1*x3*x2 + x2*x1*x3 + x2*x3*x1 + x3*x1*x2 + x3*x2*x1, one must divide by (# of x's here)! = 3! to count only x1*x2*x3.

Which is the final step of the proof.

I note as an aside the Weyl character formula for SU(2) irreps. Using quantum-mechanical notation with angular momentum j, the highest weight becomes {2j}, and the character

X(j,c) = exp((-j/2)*c) + exp((-j/2+1)*c) + ... + exp((j/2)*c)
= sinh((j+1)*c/2)/sinh(c/2)

lpetrich

Not surprisingly, large reps take a long time to process. I tried the square plethysm of the fundamental irrep of E8, the "248", and I found about 3 minutes for Mathematica and 1 minute for Python on my 2-GHz Intel Core 2 Duo iMac

starting
(0, 0, 0, 0, 0, 0, 1, 0): 248
symmetric part of square
1 of (0, 0, 0, 0, 0, 0, 2, 0): 27000
1 of (1, 0, 0, 0, 0, 0, 0, 0): 3875
1 of (0, 0, 0, 0, 0, 0, 0, 0): 1
antisymmetric part of square
1 of (0, 0, 0, 0, 0, 1, 0, 0): 30380
1 of (0, 0, 0, 0, 0, 0, 1, 0): 248
Using Robert Cahn's root-numbering convention

I've been considering alternatives to both Mma and Python in the hope of getting a speedup. They have varying advantages and disadvantages.

Here are some benchmarking sites:
Which programming languages are fastest? | Computer*Language*Benchmarks*Game, the most useful that I've found
Programming Languages Benchmarks
The winners are Fortran, C, C++, Ada, ATS, and Java.

Neither Ada nor ATS is very well-known, so I'll ignore them. The next question is whether the languages have convenience features for some of the more common tasks. Associative arrays are good for finding whether a vector is already in a list of vectors, a very common task in my code, and operator overloading is good for implementing fractions and abstract-algebra entities.

Programming Language Comparison has a section on operator overloading, and of the fast languages, only C++ has it. Looking in Wikipedia, I find that of those languages, only C++ and Java have them.

So I choose C++. With its templates and operator overloading, one can implement fractions, and with its Standard Template Library template "map", one can implement associative arrays.

Simon_Tyler

Have you tried just using compilation to C in both Mathematica and Python?
http://reference.wolfram.com/mathema...ionTarget.html
http://docs.cython.org/src/quickstart/overview.html
The generated C code can sometimes be almost as fast or faster than a hand coded C/C++
http://fredrik-j.blogspot.com/2009/0...e-days-15.html <-- this blog is worth reading
http://carlo-hamalainen.net/blog/200...vs-c-improved/
It also makes parallelization easier:
http://www.walkingrandomly.com/?p=3635
http://reference.wolfram.com/mathema...fficiency.html

The advantage of using these is that you only need to compile the key steps in your code and only need to make minimal modifications to your existing code.

The Mathematica Compile[..., CompilationTarget->"C"] function will fall back to normal virtual machine compilation if no C compiler exists.
And for Cython, you can use python decorators to achieve a similar thing
http://docs.cython.org/src/tutorial/pure.html
This is even easier if you're using sage
http://docs.cython.org/src/quickstar...highlight=sage

The other thing to try is to use (integer) packed arrays and make sure that you don't accidentally unpack them.
For Mathematica, see
http://www.wolfram.com/technology/guide/PackedArrays/
http://stackoverflow.com/questions/4...in-mathematica
http://stackoverflow.com/questions/8...for-using-them
For pure Python see
http://docs.python.org/library/array.html
but maybe it would be easier and faster to use Numpy
http://docs.scipy.org/doc/numpy/user/basics.types.html
http://stackoverflow.com/questions/8...arrays-so-fast

Of course, if you do completely rewrite your code in C or C++ then you can create interfaces for Mathematica and Python (and most other languages). However, your code will also become less accessible to most people.

Anyway, sorry for another link dump. Hopefully you find something useful in it!