Pictographs for intertwiners of SU(n): O-blades

O-BLADES ASSOCIATED WITH AN ADMISSIBLE TRIPLE OF IRREDUCIBLE REPRESENTATIONS OF SU(n)

Version 2.0
Robert Coquereaux
CPT, CNRS, Luminy, Marseille

WHAT THE PROGRAM DOES

A triple of irreducible representations (ω1, ω2, ω3) of SU(n) specified by their highest weights is admissible if ω1 ⊗ ω2 ⊗ ω3 contains ℂ (the scalars), or equivalently if the decomposition of ω1 ⊗ ω2 into a direct sum contains (OverBar[ω3]) (the conjugate of ω3). By definition an admissible triple has non-zero multiplicity,  i.e. the representation  (OverBar[ω3]) may occur once or several times in the decomposition of ω1 ⊗ ω2 into a direct sum.
If the multiplicity is equal to m one can associate m diagrams, generically called pictographs to the chosen triple.

There are several versions of pictographs: 
Chronologically the first version was introduced by A.D. Berenstein and A.Z. Zelevinsky in 1992; they are called BZ-triangles. 
Then we had the honeycombs introduced by A. Knutson and T. Tao in 1992; they are called KT-honeycomb  (actually they were introduced for the group GL(n)).
A kind of dual version of the latter was later introduced by A. Ocneanu (unpublished); they are called O-blades, or oblades, and were presented and used, together with an adapted SU(n) version of the KT-honeycomb (called SU(n) honeycombs, isometric honeycombs, or just honeycombs, for short)  in an article by R. Coquereaux and J.-B. Zuber,  in 2014.

The edges of a honeycomb carry integers : one can display this information with regular honeycombs, with edges labelled by the corresponding non-negative integers (but they may vanish).
One can also display the same information by using isometric honeycombs (also presented in the same paper by R.C. and J.-B. Zuber) where edges have precisely a length equal to the integer that they carry.
In the latter case some edges may have length 0 so that the hexagons of an isometric honeycomb can be degenerated.

The precise references can be found in the subsection DETAILS AND REFERENCES (see below).

Given an admissible triple, the present program determines its multiplicity and draws the corresponding O-blades.

A companion program (called honeycombs)  draws the associated regular honeycombs, and also draws the isometric honeycombs.

The present program can also determine the decomposition of a tensor product of irreducible representations (irreps) into a direct sum of irreps, although this is not is prime objective, as this is something that can be done separately by the function tensorproduct contained in contained in the notebook tensorProdIrrep.nb  (the function tensorproduct is defined in the package shorttensorProdIrrepPackage).  

A product of two irreps can contain a huge quantity of terms in its decomposition, whereas the functions multiplicity and plotOblades, defined below, assume that their three arguments are provided by the user, so that they do not need to use the function tensorproduct (but the user may call this latter function if needed).

Constraints obeyed by the relevant pictographs (here O-blades) are summarized in the  subsection “The SU(N) game”.

Here, like in the program tensorProdIrrep, the decomposition of a tensor product of irreducible representations is done by using an improved hive technique developed by the author: in order to speed up the calculation we use a basis of elementary pictographs where syzygies (constraints) are eliminated before the actual reduction of the system of inequalities takes place.

Note: the package obladePackage given below contains many more functions than those needed to draw the O-blades.
The user should have a look at the list of exported functions at the beginning of the package obladePackage and read their documentation (see function::usage).

What the program does can be explained in a purely combinatorial way, read the next subsection :  “The SU(N) game”.

CAPTION

Given an admissible triple, the program determines its multiplicity and draws the corresponding O-blades.
This can be explained in terms of representation theory (see the subsection Details And References), or in a purely combinatorial way, see the next subsection :  “The SU(N) game”.

THE SU(N) GAME

External vertices (sides of the triangle, excluding corners) are given: they are labelled with N-1 arbitrary non-negative integers.
Internal edges carry non-negative integers.

Problem: find all the pictographs obeying the following constraints: 

	Conservation of integers at each external vertex (its label is equal to the sum of the two internal edges ending there).
	Equality of opposite angles at each internal vertex (the value of an acute angle being defined as the sum of its adjacent edges).
	
A triple (three sides) is called admissible if the above problem admits a solution.   Arbitrary triples are usually not admissible. 
Given two sides, the set of admissible triples is finite : the third side has to be chosen in a specific list (the one returned by the  tensorproduct command).
A chosen admissible triple may have multiplicity (several pictographs with the same sides) : the various solutions differ by the values of internal edges.
The solutions are given by the pictographs.

Example
	
{3,1,2} ⊗ {2,5,1} ⟶ {3,2,5}
Multiplicity is 4
Here are the four O-blades :

DETAILS AND REFERENCES

Irreducible representations (irreps) of the Lie group SU(n) are denoted by their highest weights. 
Their  n-1 components are given in a  particular basis called the Dynkin basis.

One starts from three irreps building an admissible triple.
Warning about the orientation of the last member of the triple... the triple obtained by conjugating the last member of an admissible triple is usually not admissible!
If the user only knows the first two members of the triple, he has  first to decomposes their tensor product ω1⊗ω2 as a sum of irreps, with multiplicities.

The user chooses one of them ((OverBar[ω3])), the program then draws the pictographs (O-blade version) associated with the coupling ω1⊗ω2 ⟶ (OverBar[ω3]).
Or, equivalently, with the coupling ω1⊗ω2 ⊗ω3 ⟶  1

The multiplicity of the coupling (aka Littlewood-Richardson coefficient) is equal to the number of pictographs.

Components of highest weights (non-negative integers) are read as follows:

ω1 : external vertices belonging to the left side of the triangle, bottom to top.
ω2 : external vertices belonging to the right side of the triangle, top to bottom.
ω3 : external vertices belonging to the bottom side of the triangle, right to left (equivalently (OverBar[ω3]) : left to right)

Internal edges carry non-negative integers.
They obey the following constraints : 

	Conservation of integers at each external vertex, 
	Equality of opposite angles at each internal vertex (the value of an angle being defined as the sum of its edges).

The displayed collection of pictographs, for a given triple  (ω1, ω2 ,  ω3), solves the above constraints.

Several types of pictographs can be used to display the same space of intertwiners (Subscript[Hom, SU4])(ω1⊗ω2, (OverBar[ω3])) :  BZ-triangles, KT-honeycombs, O-blades.
BZ stands for Berenstein-Zelevinski , KT for Knutsen-Tao, and O for Ocneanu. 

Here we display O-blades.

O-blades (also called “oblades”) are introduced and documented in [1][2][3]

[1] MouseOver
[2] MouseOver
[3] MouseOver

The decomposition of the tensor product of irreducible representations is done by using an improved hive technique developed by the author and inspired from [4], [5], [6]. 
In order to speed up the calculation we use a basis of elementary pictographs where syzygies (constraints) are eliminated before the actual reduction of the system of inequalities takes place.

[4] MouseOver
[5] MouseOver
[6] MouseOver

CAVEATS AND COMMENTS

Web version:  The computation is terminated if it takes more than 1.5 seconds

Constraint on the values of weight components when using the  tensorproduct command :
In order to run the command  tensorproduct[la, mu], the components of the weights la and mu should obey the constraint Or[max(la) < 10, max(mu)] < 10 == True
This value of 10 comes from the default setting of the global system options for Reduce in Mathematica version 11.1
Strangely enough, the constraint was much less drastic -- it was 10^4 rather than 10^1 -- in Mathematica 6.0
This behavior will be “repaired” (hope so) in a later version of Mathematica.
Meanwhile,  If you need to manipulate simultaneously two weights that have components bigger than 10 (but smaller than 20, for instance),  quit the kernel, start it again, run the command 
SetSystemOptions[“ReduceOptions” -> {“DiscreteSolutionBound” -> 20}]
and evaluate the two packages adjmatPackage and obladePackage. 
You can also run the command 
SetSystemOptions[“ReduceOptions”->{“DiscreteSolutionBound”->Infinity}]
and evaluate the two packages adjmatPackage and obladePackage.

On the iPad the cdf programs usually work  when the tablet is in portrait (vertical) mode, but not always when it is in landscape ( horizontal) mode !
The graphs obtained with the tablet often appear to be badly scaled, although they look nice on the computer screen.
The code uses a single Manipulate wrapped into a dynamic module, in order to avoid a code using nested manipulates, something forbidden (as of Nov 2017) by the IOS cdf Reader.

The dynamic interface

The main function is plotOblade -- see also the function multiplicity,  and, if needed, the function tensorproduct.
They are defined in the main program : the obladePackage, that, in turns, calls the adjmatPackage package (exported functions adjmatfct and adjmatgraph).
The user can directly call these functions from a Mathematica session.
See the section Examples.

The section “Code for interface” contains commands that generate a dynamical interface (and a cdf file).

Why two versions for the CDF Player in version 1.0  (years < 2019) ?
The first version of the interface code can be used in Mathematica, and it was used to generate a cdf file that runs both in Mathematica and with the Wolfram CDF Player version for computers.
Unfortunately, some commands do not work well with the IOS version (iPad 2017 version) of the CDF Player.
Moreover the graphics uses GraphPlot that, in turns, uses an undocumented value of RandomSeed in order to choose a global orientation of the graphs; unfortunately the defaults for RandomSeed are not the same on different machines. This feature has changed in later releases of Mathematica and in various versions of the cdf Player; in that case the a-list nntoseed that works up to n=12 and appears both in adjmatPackageIOS and in adjmatPackage had to be adapted.  This was the reason for a second version specially written for tablets.
This is no longer the case.

Update 2020 : 
With Mathematica Version 12, the command GraphPlot uses an undocumented value of RandomSeed that makes the previous code (years < 2018) use a bad global orientation to display graphs, and the meaning of several options has changed. For this reason the code is now re-written without using GraphPlot : it uses the command AdjacencyGraph instead.

To do

Warning: because of the recent change of behavior of Reduce when the nb of solutions  > 10 we use Solve in the tensorProdIrrep program, but one should find a replacement for componentsEdgesAndnuc[{lam__},{mu__},{nuc__}] with three arguments.
 This is not used anyway in the shorttensorProdIrrepPackage, the latter being an adapted (because it uses Solve) and shortened version (because it only contains the tensorproduct code) of the /Users/robert/Library/Mathematica/Applications/oblade.m package.
 
When we use not only tensorproduct with two arguments, but also tensorproductCC, or the fct multiplicity, with three arguments, ie when we use the full tensorProdIrrepPackage we should worry when two lists  of components of weights involve integers larger than 10. To solve this problem one can use Reduce (instead of Solve) but modify its DiscreteSolutionBound option.

On the iPad the cdf programs usually work  when the tablet is in portrait (vertical) mode, but not always when it is in landscape ( horizontal) mode !

The graphs obtained with the tablet often appear to be badly scaled, although they look nice on the computer screen.

The code uses a single Manipulate wrapped into a dynamic module, in order to avoid a code using nested manipulates, something forbidden (as of Nov 2017) by the IOS cdf Reader.

tensorproduct[mat_]:= Apply[tensorproduct,mat];(* not used in the program, but it may be useful *)

HISTORY

The package was written in 2014.
The code for the interface (dynamic module, cdf file) was written in August 2017.
It was made public on the R.C. web site in September 2017.
Mathematica version 12.0 has modified the global orientation of graphs produced by GraphPlot (undocumented variable RandomSeed), our code for O-blades was re-written in 2019 to use AdjacencyGraph instead of GraphPlot.
The same version now works for Mathematica and the CDF player both on computer and on IOS tablets (iPad).

The decomposition of the tensor product of irreducible representations is done by using an improved hive technique. 
In order to speed up the calculation we use a basis of elementary pictographs where syzygies (constraints) are eliminated before the actual reduction of the system of inequalities takes place. 
See the subsection DETAILS.

Other related programs : see the program tensorProdIrrep (that only determines the reduction of a product of irreps into a direct sum), and the program SU(n) honeycombs.

LINKS

Related programs :

Notebooks (nb files, they also include the source code for the cdf version) 
(HyperlinkURL["Tensor product of SU(n) irreducible representations (Littlewood-Richardson coefficients)", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/nb/tensorProdIrrep.nb])
(HyperlinkURL["Pictographs: oblades for SU(n)", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/nb/obladesSUn.nb])
(HyperlinkURL["Pictographs: honeycombs for SU(n)", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/nb/honeycombsSUn.nb])
(HyperlinkURL["Dimensions of irreps of SU(n)", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/nb/dimensionsOfIrrepsSUn.nb])

CDF files (standalone)
for the Wolfram CDF Player
(HyperlinkURL["Tensor product of SU(n) irreducible representations (Littlewood-Richardson coefficients). CDF file", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/tensorProdIrrep.cdf])
(HyperlinkURL["Pictographs: oblades for SU(n). CDF file", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/obladesSUn.cdf])
(HyperlinkURL["Pictographs: honeycombs for SU(n). CDF file", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/honeycombsSUn.cdf])
(HyperlinkURL["Dimensions of irreps of SU(n). CDF file", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/dimensionsOfIrrepsSUn.cdf])

CDF files adapted to the first version of the CDF Player on IOS 
(when the latter is updated,  one should use the previous files instead):
(HyperlinkURL["Pictographs: oblades for SU(n). CDF file, IOS.", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/obladesSUnForIOS.cdf])
(HyperlinkURL["Pictographs: honeycombs for SU(n). CDF file, IOS.", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/honeycombsSUnForIOS.cdf])

Files for the cdf web plug-in 
(dying technology: the plug-in is incompatible with most recent browsers, these files will not be updated):
(HyperlinkURL["Tensor product of SU(n) irreducible representations (Littlewood-Richardson coefficients). CDF file for web plug-in", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/tensorProdIrrepWeb.cdf])
(HyperlinkURL["Pictographs: oblades for SU(n). CDF file for web plug-in", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/obladesSUnWeb.cdf])
(HyperlinkURL["Pictographs: honeycombs for SU(n). CDF file for web plug-in", http://www.cpt.univ-mrs.fr/~coque/Computer_programs/cdf/honeycombsSUnWeb.cdf])

CITATION

If you use this program in a scientific publication or talk, please give proper academic credit to the author :

MouseOver
MouseOver
September 2017

LICENCE

MIT License
Copyright (c) 2017 RobertCoquereaux
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sub license, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHERLIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Packages for Mathematica and code for a user interface (for Mathematica or for the CDF Player on a computer)

Use : Evaluate Initialization Cells.
The Manipulate window appears when outputObladeManipulate is evaluated (no argument).

Two packages (to be evaluated)

These two packages are local copies taken from robert/Library/Mathematica.Applications/*

adjmatPackage

BeginPackage["adjmatPackage`"];
(*usages go here*)

adjmatfct::usage="adjmatfct[n] returns the adjacency matrix of SU(3) at level n"

adjmatgraph::usage="adjmatgraph[n] returns the graph of the adjacency matrix of SU(3) at level n"

vertexcoordtriangle::usage=
"vertexcoordtriangle[n] specifies the coordinates of the inner vertices of the graph drawn by the function adjmatgraph"

(* nntoseed::usage = "nntoseed is association list and nntoseed[NN] gives the RandomSeed to be used in GraphPlot (pre 2018 version) in order to display correctly (orientation) the graph of adjacency matrices. The list nntoseed should be completed manually." *)

Begin["Private`"];
(*Function definitions go here*)
adjmatfct[NN_]:=adjmatfct[NN]=Module[{levelplus,rk,One,iii,SparseAdjTens,normalAdjTens,lab,partlab,adjacencyMatrices},
levelplus = NN +1;
rk=(NN + 1)(NN +2)/2;
One = IdentityMatrix[rk];
SparseAdjTens=SparseArray[{{a_,b_,c_,d_} /;  c == a+1 &&  d == b ||c == a &&  d == b - 1 || c == a-1 && d == b + 1 -> 1},{levelplus,levelplus,levelplus,levelplus}];
normalAdjTens= Normal[SparseAdjTens];
lab={{0,0},{1,0},{0,1}};
Do[Do[iii=i+1;partlab[iii]={j-i,i},{i,0,j}];Do[lab=Append[lab,partlab[i]],{i,1,j+1}],{j,2,levelplus-1}];
adjacencyMatrices=Module[{v,crit0, crit1, crit2, indices0,indices1,indices2,indices,orderedindices,smallAdjTens, matbid,trisort,tritolab,permutritolab,permutritolabinverse,AdjMatTriality,AdjMatLevel},
indices = Flatten[Table[v[i,j],{i,1,levelplus},{j,1,levelplus-i+1}]];
crit0[v[p_,q_]] := Mod[p + 2q - 3,3] ==0; 
crit1[v[p_,q_]] := Mod[p + 2q - 3,3] ==1; 
crit2[v[p_,q_]] := Mod[p + 2q - 3,3] ==2; 
indices0=Select[indices,crit0];
indices1=Select[indices,crit1];
indices2=Select[indices,crit2];
orderedindices=Join[indices0,indices1,indices2];
smallAdjTens= Table[normalAdjTens[[i,j,k,l]],{i,1,levelplus},{j,1,levelplus-i+1},{k,1,levelplus},{l,1,levelplus-k+1}];
matbid[v[a_,b_],v[c_,d_]] :=smallAdjTens[[a,b,c,d]];
AdjMatTriality=Table[matbid[orderedindices[[i]],orderedindices[[j]]],{i,1,rk},{j,1,rk}];
(* Generator [1,0] with vertices sorted by triality blocks *)
trisort=(orderedindices /. {v->List})-ConstantArray[{1,1},rk];
tritolab=Table[Flatten[List[s,Position[lab,trisort[[s]]]]],{s,1,rk}];
permutritolab=SparseArray[MapThread[Rule,{tritolab, ConstantArray[1,rk]}],{rk,rk}];
permutritolabinverse=SparseArray[Inverse[permutritolab]];
 AdjMatLevel = permutritolabinverse.AdjMatTriality.permutritolab ;
(* Generator [1,0] with vertices sorted by increasing levels, like lab *)
{AdjMatTriality,AdjMatLevel}];
adjacencyMatrices[[2]]]

(* 
nntoseed=<|1->22,2->8,3->1,4->2,5->3,6->0,7->4,8->2,9->0,10->10,11->1,12->9|>;
(* nntoseed : to be completed *)
 
adjmatgraph[NN_]:= If[MissingQ[nntoseed[NN]]==True, Print["The a-list nntoseed should be completed for this argument"],
GraphPlot[adjmatfct[NN], VertexLabeling->True,Method->{Automatic,"RandomSeed"->nntoseed[NN]}]]
old version *)

vertexcoordtriangle[s_]:=Module[{nw={-1/2,Sqrt[3]/2}},Flatten[Table[{p,0}+q nw,{p,0,s},{q,0,p}],1]//N]

(* adjmatgraph[NN_]:= 
GraphPlot[adjmatfct[NN], VertexLabeling->True,VertexCoordinateRules->vertexcoordtriangle[NN]] *)
(*  March 2019 - Before Mathematica version 12 *)

(* New Avril 2020 : Mathematica Version 12 : lack of backward compatibility for GraphPlot. Rewriting  code with AdjacencyGraph *)

(* The test If[$VersionNumber<12, GraphPlot[... etc fails on iPad -  I keep it ... *)

(* adjmatgraph[NN_]:= If[$VersionNumber<12,
GraphPlot[adjmatfct[NN], VertexLabeling->True,VertexCoordinateRules->vertexcoordtriangle[NN]],
AdjacencyGraph[adjmatfct[NN],VertexLabels->All,VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False]] *)

adjmatgraph[NN_]:= AdjacencyGraph[adjmatfct[NN],VertexLabels->All,VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False]

End[];(*private*)(*protect what you want*)EndPackage[];

Names["adjmatPackage`*"]

obladePackage (calling adjmatPackage)

BeginPackage["obladePackage`", {"adjmatPackage`"}];
(*usages go here*)

setGroupSU::usage=
"Evaluation of setGroupSU[n] initializes all relevant data for SU(n). It returns n and set the variable currentgroup to n.";

(* old name for setGroupSU[n] was oblade[n] *)

currentgroup::usage = "currentgroup returns the argument n of the group SU(n) currently in use.";

tensorproduct::usage="tensorproduct[lambda,mu] returns the decomposition into irreducible components nu of a tensor product of two irreps lambda, mu, of the Lie group SU(n), given by the components of their highest weights along the basis of fundamental weights.  The last entry of each term is its multiplicity.
One can use this command with three arguments (tensorproduct[lambda,mu,nu]), and only some of the components of the last argument, nu, can be given, in which case all the possible solutions are listed.
It uses a hive-type algorithm (oblades).
If the group is SU(3) one can use the command SU3MultiplicityDynkin[lambda, mu] which is faster and uses the semi-magic algorithm.
In the case of SU(3) again, one can use the command tensorproductGenFun which uses a generating function.
";

SU3MultiplicityDynkin::usage="SU3MultiplicityDynkin[lambda, mu] returns the decomposition into irreducible components nu of a tensor product of two irreps lambda, mu, of the Lie group SU(3), given by the components of their highest weights along the basis of fundamental weights.
It uses the semi-magic algorithm (this only works for SU(3)). The general command is tensorproduct.";

multiplicity::usage="multiplicity[lambda,mu,nu] gives the multiplicity of the term nu in the decomposition into irreducible components of a tensor product lambda ⊗ mu of two irreps of the Lie group SU(n), given by the components of their highest weights along the basis of fundamental weights.
This integer is also called Littlewood-Richardson coefficient.  Equivalently it is the multiplicity of the trivial representation in the triple tensor product: lambda ⊗ mu ⊗ Conjugate[nu].
This function admits optional arguments.
By default it calls the function multiplicityHive, that uses the hive (oblade) method, and works for all SU(n).
With the option Kostant->True multiplicity calls the function multiplicityKostant that only works for those Lie groups for which the Kostant partition function is known.
With the option Magic->True multiplicity calls the function multiplicityMagic that uses the semi-magic algorithm and only works for SU(3). Very fast.
With the option Volume->True multiplicity calls the function volumeJ that uses the volume algorithm and only works for SU(3).Very fast.
With the option GenFun->True multiplicity calls the function multiplicityGenFun that uses a generating function that only works for SU(3). Slow.
With the option LiE->True multiplicity calls the function multiplicityLiE that uses the external program LiE.
";

multiplicityHive::usage="multiplicityHive[lambda,mu,nu] gives multiplicity[lambda,mu,nu] by using the hive (oblade) method, which is the default. See the documentation of the command multiplicity.";

multiplicityKostant::usage="multiplicityKostant[lambda,mu,nu] gives multiplicity[lambda,mu,nu] by using the Kostant method.
It is equivalent to multiplicity[lambda,mu,nu, Kostant->True].
It only works for those Lie groups for which the Kostant partition function is known. See the documentation of the command multiplicity.
Private: it is called lr in other notebooks.";

Kostant::usage ="An option for multiplicity";
Magic::usage ="An option for multiplicity";
GenFun::usage ="An option for multiplicity";
LiE::usage ="An option for multiplicity";

multiplicityMagic::usage="multiplicityMagic[lambda,mu,nu] gives multiplicity[lambda,mu,nu] by using the semi-magic method (this only works for SU(3).
It is equivalent to multiplicity[lambda,mu,nu, Magic->True].
It only works for SU(3). See the documentation of the command multiplicity.";

multiplicityLiE::usage="multiplicityLiE[lambda,mu,nu] gives multiplicity[lambda,mu,nu] by calling the external program LiE.
It is equivalent to multiplicity[lambda,mu,nu, LiE->True]. The group is assumed to be of type SU.
See the documentation of the command multiplicity.";

fundamentalcomponentsLabels::usage = 
"Prints the chosen basis of fundamental intertwiners (triplets of related external vertices, using a clockwise ordering of the external sides of an SU(n) oblade)" ;

pictographsCoordinates::usage =
"Assuming that the chosen triple has multiplicity m, the function pictographsCoordinates[lambda,mu,nu] gives the matrix
 of components of the m different pictographs along a basis of fundamental intertwiners. 
For SU(n) there are 3(n-1)+ (n-1)(n-2) fundamental intertwiners. A basis has only 3(n-1)+ (n-1)(n-2)/2 elements because of relations between fundamental intertwiners (syzygies). 
The rank of the matrix, minus 1, gives the dimension of the hive polytope, a rational polytope (not always integral) whose integral points are precisely the different pictographs. 
The dimension of the hive polytope is at most (n-1)(n-2)/2, with equality for the so-called generic triples." ;  

plotOblades::usage ="plotOblades[lambda,mu,nu] returns the list of drawings of all oblades with sides defined by the three highest weights lambda, mu, nu and describes the branching la ⊗ mu -> nu. 
To display this list as a tableform (with no braces) one can do as follows:  li = plotOblades[lambda,mu,nu];  GraphicsRow[li,ImageSize-> Length[li]*300]. 
The command admits and optional argument pos which is an integer chosen between 1 and multiplicityHive[lambda, mu, nu]; 
if specified, plotOblades[lambda,mu,nu, pos] displays only the particular oblade numbered pos.";

edgesEWvalues::usage = "edgesEWvalues[lambda,mu,nucc] returns the list of East-West edges values for the branching rule lambda ⊗ mu ⊗ nucc -> 1 or la ⊗ mu -> nu with nuc = Reverse[nu].";
edgesNESWvalues::usage = "edgesEWvalues[lambda,mu,nucc] returns the list of NorthEast-SouthWest edges values for the branching rule lambda ⊗ mu ⊗ nucc -> 1 or la ⊗ mu -> nu with nuc = Reverse[nu].";
edgesNWSEvalues::usage = "edgesEWvalues[lambda,mu,nucc] returns the list of NorthWest-SouthEast edges values for the branching rule lambda ⊗ mu ⊗ nucc -> 1 or la ⊗ mu -> nu with nuc = Reverse[nu].";

pictographsCoordinatesToWeights::usage= "pictographsCoordinatesToWeights[list_of_components] returns the Dynkin labels of the three highest weights determined by the components of the chosen pictograph (intertwiner)";

pictographsCoordinatesToDynkinAndLeaves::usage = "pictographsCoordinatesToDynkinAndLeaves[list_of_components] returns a pair. 
The first member of the pair give the Dynkin labels of the three highest weights determined by the components of the chosen pictograph (intertwiner).
The second member is a list of the values of the leaves of the pictograph, in the three directions.
";

pictographsCoordinatesToHoneycomb::usage = "pictographsCoordinatesToHoneycomb[list_of_components] draws the honeycomb determined by the components of the chosen pictograph (intertwiner). 
This is done without re-calculating the whole list of pictographs determined by the associated external three weights.";

pictographsCoordinatesToOblade::usage = "pictographsCoordinatesToOblade[list_of_components] draws the oblade determined by the components of the chosen pictograph (intertwiner).
This is done without re-calculating the whole list of pictographs determined by the associated external three weights.";

nondegeneratedOblades::usage = "List of non degenerated pictographs for a given tensor product (the command accepts two or three arguments)";

hivePolytope::usage = "hivePolytope[lambda, mu, nu] returns a pair (vars, ineqs) where vars are independent leaves, and ineqs are inequations in the variables vars defining the hive polytope associated with the three highest weights lambda, mu, nu.";

multiplicityGeneratingFunctionSU3::usage = "multiplicityGeneratingFunctionSU3[L1,L2,M1,M2,N1,N2] is a generating function for the multiplicity of L⊗M⊗N --> 1. 
It is only implemented for SU(3). It can be used to determine the tensor product of two irreps (see tensorproductGenFun) or the multiplicity of a triple (see multiplicityGenFun)";

tensorproductGenFun::usage = "tensorproductGenFun[{l1,l2},{m1,m2}] returns the tensor product of the two irreps (with multiplicities).
It uses the generating function for multiplicities given by multiplicityGeneratingFunctionSU3. This is a SU3 specific command.";

multiplicityGenFun::usage = "multiplicityGenFun[lambda,mu,nu] returns the multiplicity of  lambda⊗mu -> nu.
It uses the generating function for multiplicities given by multiplicityGeneratingFunctionSU3. This is a SU3 specific command.";

Begin["Private`"];
(*Function definitions go here*)
 
(* Sides,edges,components and equations *)
(* Sides *)
(* p will run from 1 to NN-1 *)
lambdalabel[n_,p_]:=1/2 (2+3 p+p^2);
mulabel[n_,p_]:=(n+1)(n+2)/2-p;
nulabel[n_,p_]:=1/2 (n+n^2-2 n p+(-1+p) p+2);

SetSystemOptions["ReduceOptions"->{"DiscreteSolutionBound"->Infinity}];

setGroupSU[var_]:= With[{NN=var},

currentgroup=NN;
 
(* Equations on sides *)
(* conservation equations at the external vertices (weights) *)
lambdaequations = And[
Apply[And,Table[lambda[p,0,0]== Sum[ j[p,NN-(s+p),s] ,{s,0,NN-p}],{p,1,NN-1}]],
Apply[And,Table[lambda[0,p,0]== Sum[ j[s,p,NN-(s+p)] ,{s,0,NN-p}],{p,1,NN-1}]],
Apply[And,Table[lambda[0,0,p]== Sum[ j[s,NN-(s+p),p] ,{s,0,NN-p}],{p,1,NN-1}]]];
(* Edges values (in terms of fundamental pictographs)*)
edgesEW[NN]= Table[FoldList[Plus][Table[j[t,NN-s,s-t],{s,t,NN-1}]],{t,1,NN-1}];
edgesNESW[NN]= Table[FoldList[Plus][Table[j[s-p,p,NN-s],{s,p,NN-1}]],{p,1,NN-1}];
edgesNWSE[NN]= Table[FoldList[Plus][Table[j[NN-s,s-p,p],{s,p,NN-1}]],{p,1,NN-1}];
(* Edge labelling doublets *)
edgesEWlabels[lev_]:=Table[Take[Table[1/2 (4+p+p^2)+(lev-1),{p,lev,NN}],{z,z+1}],{z,1,NN-lev}];
edgesNESWlabels[lev_]:=Table[Take[Table[1/2 (3p+p^2)+1/2 (lev-1) (lev)+ (lev-1) p,{p,1,NN}],{z,z+1}],{z,1,NN-lev}];
edgesNWSElabels[lev_]:=Table[Take[Range[1+1/2 (-lev+NN) (1-lev+NN),1-lev+NN+1/2 (-lev+NN) (1-lev+NN)],{z,z+1}],{z,1,NN-lev}];

(* Fundamental pictographs and positivity constraints *)
components = Join[
Table[Table[ j[p,NN-(s+p),s] ,{s,0,NN-p}],{p,1,NN-1}],
Table[Table[ j[s,p,NN-(s+p)] ,{s,0,NN-p}],{p,1,NN-1}],
Table[Table[ j[s,NN-(s+p),p] ,{s,0,NN-p}],{p,1,NN-1}]];
fundamentalcomponents = DeleteDuplicates[Flatten[components]];
positivityconstraints=And[
Apply[And,Map[# >= 0 &,Flatten@edgesEW[NN]]],
Apply[And,Map[# >= 0 &,Flatten@edgesNESW[NN]]],
Apply[And,Map[# >= 0 &,Flatten@edgesNWSE[NN]]]
];
buildlambdainput[{l1__},{l2__}]:=Flatten[{
MapThread[Rule,{Table[lambda[p,0,0],{p,1,NN-1}],{l1}}],
MapThread[Rule,{Table[lambda[0,p,0],{p,1,NN-1}],{l2}}]}];
buildlambdainput[{l1__},{l2__},{l3__}]:=Flatten[{
MapThread[Rule,{Table[lambda[p,0,0],{p,1,NN-1}],{l1}}],
MapThread[Rule,{Table[lambda[0,p,0],{p,1,NN-1}],{l2}}],
MapThread[Rule,{Table[lambda[0,0,p],{p,1,NN-1}],{l3}}]}];
fundamentalcomponentsLabels=ToString[fundamentalcomponents];

(* Values of edges if only lambda and mu are given *) 
componentsEdgesAndnuc[{lam__},{mu__}]:= componentsEdgesAndnuc[{lam},{mu}]= Module[{redsol,locrules,nuclist,fundcompval,edgesEWValues,edgesNESWValues,edgesNWSEValues},
redsol=Reduce[And[lambdaequations/. buildlambdainput[{lam},{mu}],positivityconstraints],fundamentalcomponents,Integers];
locrules=ToRules[redsol];
fundcompval=fundamentalcomponents//. {locrules} ;
nuclist=Table[lambda[0,0,p],{p,1,NN-1}]//. {locrules};
edgesEWValues = Table[edgesEW[NN]/.MapThread[Rule,{fundamentalcomponents, fundcompval[[s]]}],{s,1,Length[fundcompval]}];
edgesNESWValues = Table[edgesNESW[NN]/.MapThread[Rule,{fundamentalcomponents, fundcompval[[s]]}],{s,1,Length[fundcompval]}];
edgesNWSEValues = Table[edgesNWSE[NN]/.MapThread[Rule,{fundamentalcomponents, fundcompval[[s]]}],{s,1,Length[fundcompval]}];
{fundcompval,edgesEWValues ,edgesNESWValues,edgesNWSEValues,nuclist}];

(* Values of edges if lambda, mu and nu are given - Warning: nuc is reverse nu *) 
componentsEdgesAndnuc[{lam__},{mu__},{nuc__}]:= componentsEdgesAndnuc[{lam},{mu},{nuc}]= Module[{redsol,locrules,nuclist,fundcompval,edgesEWValues,edgesNESWValues,edgesNWSEValues},
redsol=Reduce[And[lambdaequations/. buildlambdainput[{lam},{mu},{nuc}],positivityconstraints],fundamentalcomponents,Integers];
redsol = And[redsol, Apply[And,Table[lambda[0,0,p]=={nuc}[[p]],{p,1,NN-1}]]];
(* If[redsol=!=False,locrules= ToRules[redsol], Print["Not an admissible intertwiner"];Return[{0,0,0,0,0}]];*)
If[redsol=!=False,locrules= ToRules[redsol],Return[{0,0,0,0,0}]];
fundcompval=fundamentalcomponents//. {locrules} ;
nuclist=Table[lambda[0,0,p],{p,1,NN-1}]//. {locrules};
edgesEWValues = Table[edgesEW[NN]/.MapThread[Rule,{fundamentalcomponents, fundcompval[[s]]}],{s,1,Length[fundcompval]}];
edgesNESWValues = Table[edgesNESW[NN]/.MapThread[Rule,{fundamentalcomponents, fundcompval[[s]]}],{s,1,Length[fundcompval]}];
edgesNWSEValues = Table[edgesNWSE[NN]/.MapThread[Rule,{fundamentalcomponents, fundcompval[[s]]}],{s,1,Length[fundcompval]}];
{fundcompval,edgesEWValues ,edgesNESWValues,edgesNWSEValues,nuclist}];

(* Plot oblades *)

AdjMatbis=adjmatfct[NN]; (* this function is defined in the dependent package adjmatPackage  *)

plotSUNoblade[mytripleofweights_,loccoupling_]:=Module[{adj,vert=mytripleofweights,grid,gr2,grla,grmu,grnu,grEW,grNESW,grNWSE},
(* Warning: mytripleofweights is a list of three weights *)
(* loccoupling is a positive integer that refers to the position in the list of coupling with fixed nu. *)
adj=AdjMatbis;
fullinfo=componentsEdgesAndnuc[mytripleofweights[[1]],mytripleofweights[[2]], Reverse[mytripleofweights[[3]]]];
grid=AdjacencyGraph[adj,VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False];
grla=AdjacencyGraph[adj,VertexShapeFunction->( Which@@Flatten[Table[{#2==lambdalabel[NN,s],{Text[Style[vert[[1,s]],25,Blue],# ,Offset[{20,-30}](*, {Right,Bottom}*)], Disk[#,.03],{Black}}},{s,1,NN-1}],1]&)(* lambda side *),
VertexCoordinates->vertexcoordtriangle[NN]
,DirectedEdges->False];
grmu=AdjacencyGraph[adj,VertexShapeFunction->( Which@@Flatten[Table[{#2==mulabel[NN,s],{Text[Style[vert[[2,s]],25,Blue],#,Offset[{-30,-30}](*{Left,Bottom}*)],Disk[#,.03],{Black}}},{s,1,NN-1}],1]&)(* mu side *),
VertexCoordinates->vertexcoordtriangle[NN]
,DirectedEdges->False];
grnu=AdjacencyGraph[adj,VertexShapeFunction->( Which@@Flatten[Table[{#2==nulabel[NN,s],{Text[Style[vert[[3,NN-s]],25,Blue],#,Offset[{-5,+25}](*{Right,Top}*)],Disk[#,.03],{Black}}},{s,1,NN-1}],1]&)(* reversed nu side *),
VertexCoordinates->vertexcoordtriangle[NN]
,DirectedEdges->False];
grEW=AdjacencyGraph[adj,EdgeShapeFunction->
( Which@@Flatten[Table[Flatten[Table[{First[#2]==First[u] && Last[#2]==Last[u],{Text[Style[Reverse[fullinfo[[2,loccoupling]][[f]]][[First@Flatten[Position[edgesEWlabels[f],u]]]],20],(#1[[1]]+#1[[2]])/2]}},
{u,edgesEWlabels[f]}],1],{f,1,NN-1}],1]&),(* edges  EW *)
VertexCoordinates->vertexcoordtriangle[NN],
DirectedEdges->False
];

grNESW=AdjacencyGraph[adj,EdgeShapeFunction->
( Which@@Flatten[Table[Flatten[Table[{First[#2]==First[u] && Last[#2]==Last[u],{Text[Style[fullinfo[[3,loccoupling]][[f]][[First@Flatten[Position[edgesNESWlabels[f],u]]]],20],(#1[[1]]+#1[[2]])/2]}},
{u,edgesNESWlabels[f]}],1],{f,1,NN-1}],1]&),(* edges  NE-SW *)
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False];

grNWSE=AdjacencyGraph[adj,EdgeShapeFunction->
( Which@@Flatten[Table[Flatten[Table[{First[#2]==First[u] && Last[#2]==Last[u],{Text[Style[Reverse[fullinfo[[4,loccoupling]][[f]]][[First@Flatten[Position[edgesNWSElabels[f],u]]]],20],(#1[[1]]+#1[[2]])/2]}},
{u,edgesNWSElabels[f]}],1],{f,1,NN-1}],1]&),(* edges  NW-SE *)
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False];

Show[grid,grla,grmu,grnu,grEW,grNESW,grNWSE]];

plotAllOblades[myweights_]:= 
Table[plotSUNoblade[myweights, pos],{pos,1,Apply[multiplicity,myweights]}];

(* we add a CC at the end of tensorproduct to say that the rhs of the tensor product is reversed (complex conjugated) *)
componentsAlongFundPict[{lam__},{mu__}]:=componentsEdgesAndnuc[{lam},{mu}][[1]];
edgesEWvalues[{lam__},{mu__}]:=componentsEdgesAndnuc[{lam},{mu}][[2]];
edgesNESWvalues[{lam__},{mu__}]:=componentsEdgesAndnuc[{lam},{mu}][[3]];
edgesNWSEvalues[{lam__},{mu__}]:=componentsEdgesAndnuc[{lam},{mu}][[4]];
tensorproductCC[{lam__},{mu__}]:=componentsEdgesAndnuc[{lam},{mu}][[5]];
(* Getting only the weights nuc (ie conjugated nu),tallied,and their multiplicities *)

tensorproduct[{lam__},{mu__}]:=Sort@Switch[
ValueQ[currentgroup],False,setGroupSU[Length[{lam}]+1];tensorproduct[{lam},{mu}],
True,If[
Length[{lam}]== currentgroup-1, 
SortBy[Tally[Map[Reverse,tensorproductCC[{lam},{mu}]]],Last], 
setGroupSU[Length[{lam}]+1];
tensorproduct[{lam},{mu}]]];

tensorproduct[{lam__},{mu__},{nu__}]:=Switch[
ValueQ[currentgroup],False,setGroupSU[Length[{lam}]+1];tensorproduct[{lam},{mu},{nu}],
True,If[
Length[{lam}]== currentgroup-1, 
SortBy[Tally[Map[Reverse,tensorproductCC[{lam},{mu},Reverse[{nu}]]]],Last], 
setGroupSU[Length[{lam}]+1];
tensorproduct[{lam},{mu}, {nu}]]];

(* Information for a given admissible triple:components along primitive pictographs,edges,multiplicity,etc. *)
(* Warning:the last component of the next command is a nu_conjugate *)
(* Getting part of the above information: *)
(* we add a CC at the end of tensorproduct to say that the rhs of the tensor product is reversed (complex conjugated) *)
(* Getting only the weights nuc (ie conjugated nu),tallied,and their multiplicities *)

componentsAlongFundPict[{lam__},{mu__},{nuc__}]:=componentsEdgesAndnuc[{lam},{mu},{nuc}][[1]];
edgesEWvalues[{lam__},{mu__},{nuc__}]:=componentsEdgesAndnuc[{lam},{mu},{nuc}][[2]];
edgesNESWvalues[{lam__},{mu__},{nuc__}]:=componentsEdgesAndnuc[{lam},{mu},{nuc}][[3]];
edgesNWSEvalues[{lam__},{mu__},{nuc__}]:=componentsEdgesAndnuc[{lam},{mu},{nuc}][[4]];
tensorproductCC[{lam__},{mu__},{nuc__}]:=componentsEdgesAndnuc[{lam},{mu},{nuc}][[5]];

multiplicityHive[{lam__},{mu__},{nu__}]:=Switch[
ValueQ[currentgroup],False,setGroupSU[Length[{lam}]+1];multiplicityHive[{lam},{mu},{nu}],
True,If[Length[{lam}]== currentgroup-1, 
Length[tensorproductCC[{lam},{mu},Reverse[{nu}]]], 
setGroupSU[Length[{lam}]+1];
multiplicityHive[{lam},{mu},{nu}]]];
(* multiplicityCC[{lam__},{mu__},{nuc__}]:=Length[tensorproductCC[{lam},{mu},{nuc}]]; *)

pictographsCoordinates[{lam__},{mu__},{nu__}]:=Switch[
ValueQ[currentgroup],False,setGroupSU[Length[{lam}]+1];pictographsCoordinates[{lam},{mu},{nu}],
True,If[Length[{lam}]== currentgroup-1, 
componentsAlongFundPict[{lam},{mu},Reverse[{nu}]], 
setGroupSU[Length[{lam}]+1];
pictographsCoordinates[{lam},{mu},{nu}]]];(* Warning: the third argument nu is not reversed *)

plotOblades[{lam__},{mu__},{nu__}]:=
Switch[ValueQ[currentgroup],
False,setGroupSU[Length[{lam}]+1];plotOblades[{lam},{mu},{nu}],
True,If[Length[{lam}]==currentgroup-1,plotAllOblades[{{lam},{mu},{nu}}],
setGroupSU[Length[{lam}]+1];plotOblades[{lam},{mu},{nu}]]];

plotOblades[{lam__},{mu__},{nu__},pos_]:=
Switch[ValueQ[currentgroup],
False,setGroupSU[Length[{lam}]+1];plotOblades[{lam},{mu},{nu},pos],
True,If[Length[{lam}]==currentgroup-1,plotSUNoblade[{{lam},{mu},{nu}},pos],
setGroupSU[Length[{lam}]+1];plotOblades[{lam},{mu},{nu},pos]]];

(* Improve the way of introducing the global variable currentgroup *)

pictographsCoordinatesToWeights[{li__}]:=Module[
{locjjpol = {li}.ToExpression[fundamentalcomponentsLabels],jj},
locjjpol=locjjpol /. { Private`j-> jj};
jj[u_,v_,w_]:= jj[la,u]+jj[mu,v]+jj[nu,w];
locjjpol=locjjpol /.{ jj[la,0]-> 0, jj[mu,0]->0,jj[nu,0]->0};
With[{pol =locjjpol},
{Table[Coefficient[pol,jj[la,p]], {p,1,currentgroup-1}],
Table[Coefficient[pol,jj[mu,p]], {p,1,currentgroup-1}],Table[Coefficient[pol,jj[nu,p]], {p,1,currentgroup-1}]}]];

(* The following command pictographsCoordinatesToWeightsBis re-calculates the external weights from the given list of coordinates by decomposing the intertwiner on the basis of fundamental ones.
It is not used in the sequel but it should be kept. 
*)

pictographsCoordinatesToWeightsBis[{li__}]:=Module[{check,p,n,nbdegfree},
check=Reduce[Length[{li}]==1/2 (-1+p) (4+p)&&p>0,p, Integers];
If[check==False,Print["The number of components should be equal to(-1+n) (4+n)/2 for some positive integer n"];Return[{}],
n= p/. ToRules[check]];nbdegfree=Length[{li}];
Apply[Plus,{li}*Map[pictographsCoordinatesToWeights, Table[ UnitVector[nbdegfree,s],{s,1,nbdegfree}]]]];

pictographsCoordinatesToDynkinAndLeaves[{li__}]:=Module[{unitvecs,edges,leavesEW,leavesNESW,leavesNWSE},
unitvecs=Module[{n=currentgroup,nbdegfree},nbdegfree = 3(n-1)+(n-1)(n-2)/2; Table[ UnitVector[nbdegfree,s],{s,1,nbdegfree}]];
edges=Apply[Plus,{li}*Map[pictographsCoordinatesToWeights,unitvecs]];
leavesEW=Apply[Plus,{li}*Apply[edgesEWvalues,Map[pictographsCoordinatesToWeights,unitvecs],1]];
leavesNESW=Apply[Plus,{li}*Apply[edgesNESWvalues,Map[pictographsCoordinatesToWeights,unitvecs],1]];
leavesNWSE=Apply[Plus,{li}*Apply[edgesNWSEvalues,Map[pictographsCoordinatesToWeights,unitvecs],1]];
{edges,{leavesEW[[1]],leavesNESW[[1]],leavesNWSE[[1]]}}];

pictographsCoordinatesToHoneycomb[{li__}]:=
With[{u = {Sqrt[3]/2,1/2},v = {-(Sqrt[3]/2),1/2},w={0,-1}},
 Module[{weightsandleaves,weights, leaves, la,mu,nu,nn,ptg,gra1,gra2, tailsize},
weightsandleaves=pictographsCoordinatesToDynkinAndLeaves[{li}];
weights=First@weightsandleaves;
la=weights[[1]];mu=weights[[2]];nu=Reverse@weights[[3]];
leaves=Last@weightsandleaves;
nn=currentgroup;(* for SU(NN) *)
(* Points coordinates for the floor f of the honeycomb *)
ptg[1,1]={0,0}; (* starting position: left bottom corner *)
ptg[1,s_]:=ptg[1,s-1]+ u leaves[[3,-s/2 ,-1]]  /; EvenQ[s];(* -1 special *)
ptg[1,s_]:=ptg[1,s-1]-v leaves[[2,(s-1)/2,1]] /; OddQ[s];

ptg[f_,1]:=ptg[f-1,2] - w leaves[[1,f-1,nn-(f-1)]];
ptg[f_,s_]:=ptg[f,s-1]+u  leaves[[3,-(s/2+f-1),s/2]]/; EvenQ[s];
ptg[f_,s_]:=ptg[f,s-1]- v leaves[[2,(s-1)/2,f]]/; OddQ[s];
(* Internal: points and edges*)
tailsize=Round[Max[{Total[la],Total[mu], Total[nu]}]/20];
gra1=Graphics[Join[
Table[{Black,Point[Table[ptg[f,s],{s,1,2nn- 2f +1}]]},{f,1,nn}],
Table[{Blue, Line[Table[{ptg[f,j],ptg[f,j+1]},{j,1,2(nn-f)}]]},{f,1,nn-1}],
Table[{Blue,Line[Table[{ptg[f,j],ptg[f+1,j-1]},{j,2,2(nn-f),2}]]},{f,1,nn-1}],
Table[{Red,Table[Private`affineSegment[ptg[f,j],tailsize*w],{j,1,2 nn-f,2}]},{f,1,1}],
Table[{Red,Table[Private`affineSegment[ptg[f,j],tailsize*v],{j,1,1}]},{f,1,nn}],
Table[{Red,Table[Private`affineSegment[ptg[f,j],tailsize*u],{j,2(nn-f)+1,2(nn-f)+1}]},{f,1,nn}]
],AspectRatio->Automatic];
(* External: outgoing edges *)
gra2=Graphics[{
Table[Table[Private`segmentLabel[leaves[[1,f,-j]],Black,ptg[f,2j],ptg[f+1,2j-1],{5,0}],{j,1,nn-f}],{f,1,nn-1}],
Table[Table[Private`segmentLabel[leaves[[2,f, j]],Black,ptg[j,2f],ptg[j,2f+1],{-6,-6}],{j,1,nn-f}],{f,1,nn-1}],
Table[Table[Private`segmentLabel[leaves[[3,f, -j]],Black,ptg[j,2nn-2j-2f+1],ptg[j,2nn - 2j - 2f+1+1],{0,-6}],{j,1,nn-f}],{f,1,nn-1}],

Table[Private`segmentLabel[la[[f]],Blue,Bold,ptg[f,1]-w,ptg[f+1,1]-w,{25,-10}],{f,1,nn-1}],
Table[Private`segmentLabel[mu[[nn-f]],Blue,Bold,ptg[f,2nn-2f+1]-w,ptg[f+1,2nn-2f-1]-w,{-40,-10}],{f,1,nn-1}],
Table[Private`segmentLabel[nu[[nn-f]],Blue,Bold,ptg[1,2nn-2f-1]-w,ptg[1,2nn-2f+1]-w,{-10,30}],{f,1,nn-1}]
}, AspectRatio->Automatic];
Show[gra1,gra2]]];
(* The test If[$VersionNumber<12, GraphPlot[... etc fails on iPad - I replaced it by False so the first clause is useless but I keep it ... *)

(* for $VersionNumber ≥ 12 : *)
pictographsCoordinatesToOblade[{li__}]:=
Module[{weightsandleaves,weights,la,mu,nu,leaves,adj,grid,grla,grmu,grnu, grEW,grNESW,grNWSE},

weightsandleaves=pictographsCoordinatesToDynkinAndLeaves[{li}];
weights=First@weightsandleaves;
la=weights[[1]];mu=weights[[2]];nu=Reverse@weights[[3]];
leaves=Last@weightsandleaves;
adj=AdjMatbis;

grid=AdjacencyGraph[adj,VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False];

grla=AdjacencyGraph[adj,VertexShapeFunction->( Which@@Flatten[Table[{#2==lambdalabel[NN,s],{Text[Style[la[[s]],Large,Blue],#,{Right,Bottom}],Disk[#,.05],{Black}}},{s,1,NN-1}],1]&)(* lambda side *),
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False
];

grmu=AdjacencyGraph[adj,VertexShapeFunction->( Which@@Flatten[Table[{#2==mulabel[NN,s],{Text[Style[mu[[s]],Large,Blue],#,{Left,Bottom}],Disk[#,.05],{Black}}},{s,1,NN-1}],1]&)(* mu side *),
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False
];

grnu=AdjacencyGraph[adj,VertexShapeFunction->( Which@@Flatten[Table[{#2==nulabel[NN,s],{Text[Style[nu[[NN-s]],Large,Blue],#,{Right,Top}],Disk[#,.05],{Black}}},{s,1,NN-1}],1]&)(* reversed nu side *),
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False
];

grEW=AdjacencyGraph[adj,EdgeShapeFunction->
( Which@@Flatten[Table[Flatten[Table[{First[#2]==First[u] && Last[#2]==Last[u],{Text[Style[Reverse[leaves[[1,f]]][[First@Flatten[Position[edgesEWlabels[f],u]]]],20],(#1[[1]]+#1[[2]])/2]}},
{u,edgesEWlabels[f]}],1],{f,1,NN-1}],1]&),(* edges  EW *)
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False
];

grNESW=AdjacencyGraph[adj,EdgeShapeFunction->
( Which@@Flatten[Table[Flatten[Table[{First[#2]==First[u] && Last[#2]==Last[u],{Text[Style[leaves[[2,f]][[First@Flatten[Position[edgesNESWlabels[f],u]]]],20],(#1[[1]]+#1[[2]])/2]}},
{u,edgesNESWlabels[f]}],1],{f,1,NN-1}],1]&),(* edges  NE-SW *)
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False];

grNWSE=AdjacencyGraph[adj,EdgeShapeFunction->
( Which@@Flatten[Table[Flatten[Table[{First[#2]==First[u] && Last[#2]==Last[u],{Text[Style[Reverse[leaves[[3,f]]][[First@Flatten[Position[edgesNWSElabels[f],u]]]],20],(#1[[1]]+#1[[2]])/2]}},
{u,edgesNWSElabels[f]}],1],{f,1,NN-1}],1]&),(* edges  NW-SE *)
VertexCoordinates->vertexcoordtriangle[NN],DirectedEdges->False]; 

Show[grid,grla,grmu,grnu,grEW,grNESW,grNWSE]];

(* Inequations defining the hive polytope *)
hivePolytope[la_,mu_,nu_]:=Module[
{externalConstraints,eliminatedLeaves,remainingLeaves},
externalConstraints = lambdaequations/. buildlambdainput[la,mu,Reverse@nu];
eliminatedLeaves = Solve[externalConstraints];
remainingLeaves = Union@Flatten[Map[Variables,Flatten[fundamentalcomponents/. eliminatedLeaves]]];
Print["The hive polytope is defined by the following pair {vars, ineqs} where vars are ", Length[remainingLeaves]," ≤ ",(NN-1)(NN-2)/2,
" independent parameters obeying the inequations ineqs:"];
{remainingLeaves,Simplify[positivityconstraints/. eliminatedLeaves]}]; 

(* List of non degenerated pictographs in a given tensor product *)

nondegeneratedOblades[λ_,μ_]:=Switch[ValueQ[currentgroup],
False,setGroupSU[Length[λ]+1];nondegeneratedOblades[λ,μ],
True,If[Length[λ]==currentgroup-1,nondegeneratedObladesMain[λ,μ],
setGroupSU[Length[λ]+1];nondegeneratedOblades[λ,μ]]];

nondegeneratedObladesMain[λ_,μ_]:= Module[
{tp=tensorproduct[λ,μ],len},
len=Length@tp;
Position[Table[
Map[FreeQ[#,0]&,Map[Flatten,Transpose@{Map[Flatten,edgesEWvalues[λ,μ,Reverse[tp[[j,1]]]]],Map[Flatten,edgesNESWvalues[λ,μ,Reverse[tp[[j,1]]]]],Map[Flatten,edgesNWSEvalues[λ,μ,Reverse[tp[[j,1]]]]]}]],{j,1,len}],True]
];

nondegeneratedOblades[λ_,μ_,ν_]:=Switch[ValueQ[currentgroup],
False,setGroupSU[Length[λ]+1];nondegeneratedOblades[λ,μ,ν],
True,If[Length[λ]==currentgroup-1,nondegeneratedObladesMain[λ,μ,ν],
setGroupSU[Length[λ]+1];nondegeneratedOblades[λ,μ,ν]]];

nondegeneratedObladesMain[λ_,μ_,ν_]:= 
Position[
Map[FreeQ[#,0]&,Map[Flatten,Transpose@{Map[Flatten,edgesEWvalues[λ,μ,Reverse[ν]]],
Map[Flatten,edgesNESWvalues[λ,μ,Reverse[ν]]],Map[Flatten,edgesNWSEvalues[λ,μ,Reverse[ν]]]}]],True];



(* multiplicity using Kostant for SUn *)
(* multiplicityKostant is called lr in other notebooks *)

Options[multiplicity]={Kostant->False, Magic->False, Volume -> False, GenFun -> False, LiE -> False};
multiplicity[lambda_,mu_,nu_, OptionsPattern[]]:= Which[
OptionValue[Kostant]==True, multiplicityKostant[lambda,mu,nu], 
OptionValue[Magic]==True, multiplicityMagic[lambda,mu,nu],
OptionValue[Volume]==True, multiplicityVolume[lambda,mu,nu],
OptionValue[GenFun]==True, multiplicityGenFun[lambda,mu,nu],
OptionValue[LiE]==True, multiplicityLiE[lambda,mu,nu],
True, multiplicityHive[lambda,mu,nu]];

multiplicityLiE[λ1_,λ2_,λ3_,diagram_]:=
 Block[{group,V1, V2,V3,  listofLieCommands ,res},
group=ToString[diagram];
V1 = StringJoin["[", StringRiffle[ToString/@λ1,", "], "]"];
V2 = StringJoin["[", StringRiffle[ToString/@λ2,", "], "]"];
V3 = StringJoin["[", StringRiffle[ToString/@λ3,", "], "]"];
listofLieCommands = {StringJoin["maxobjects ","1500000"],StringJoin["setdefault ",group ], StringJoin["tensor(",V1,",",V2,",",V3,")"]};
SetDirectory[$HomeDirectory];
Export["filequestion.dat",listofLieCommands];
res=ReadList["!source LiEfromM < filequestion.dat",Expression];
Return[First@res ]; 
ResetDirectory[]];(* warning: code to be changed for non-A because of difference of convention for Dynkin labels ordering *)

multiplicityLiE[lambda_,mu_,nu_]:= multiplicityLiE[lambda,mu,nu,StringJoin["A", ToString[currentgroup-1]]];

multiplicityMagic[lambda_,mu_,nu_]:= 
If[currentgroup =!= 3, Print["This algorithm works only for SU(3). I continue using the hive method"]; multiplicityHive[lambda,mu,nu],
magicmult[lambda,mu,Reverse[nu]]];

multiplicityKostant[lambda_,mu_,nu_]:=
Which[
currentgroup==3,
Kostant[{x_,y_}]:=If[x<0 || y< 0, 0,  1+ Min[x,y]];
dynkintokac[{a_,b_}]:= a {2/3,1/3}+ b {1/3,2/3};
weyldot[s_][{la1_,la2_}]:=weyl[s][{la1,la2}+{1,1}]-{1,1};
weyl[1][{la1_,la2_}]:={la1,la2};
weyl[2][{la1_,la2_}]:={-la1,la1+la2};
weyl[3][{la1_,la2_}]:={la1+la2,-la2};
weyl[4][{la1_,la2_}]:={la2,-(la1+la2)};
weyl[5][{la1_,la2_}]:={-(la1+la2),la1};
weyl[6][{la1_,la2_}]:={-la2,-la1};
sign[weyl[s_]]:={1,-1,-1,1,1,-1}[[s]];

If[Nand@@Map[IntegerQ,dynkintokac[lambda+mu-nu]],Print["Not an admissible intertwiner"];0,
Sum[sign[weyl[v]]sign[weyl[w]]Kostant[dynkintokac[weyldot[v][lambda]+weyldot[w][mu]-nu]],{v,1,3!},{w,1,3!}]],

currentgroup==4,
Kostant[{k1_,k2_, k3_}]:=
Piecewise[{
{0,k1<0 || k2< 0 || k3 < 0},

 { (k2+1)(k2+2)(k2+3)/6, k2 <= k1 && k2 <= k3},

 {(k1+1)(k1+2)(3k2-2k1+3)/6, k1<=k2<= k3},
{(k1+1)(k1+2)(3k3-k1+3)/6, k1<= k3<= k1+k3<= k2},
{ (k1+1)(k1+2)(3k3-k1+3)/6 - Binomial[(k3+k1-k2+2),3], k1<=k3<= k2<= k1+k3},

{(k3+1)(k3+2)(3k2-2k3+3)/6, k3<=k2<= k1},
{(k3+1)(k3+2)(3k1-k3+3)/6, k3<= k1<= k3+k1<= k2},
{ (k3+1)(k3+2)(3k1-k3+3)/6 - Binomial[(k1+k3-k2+2),3], k3<=k1<= k2<= k3+k1}

}];(* There are seven cases. Note that the last three cases are obtained from the previous three by exchanging k1 and k3.
Giving explicitly the case 0 is useless if Kostant is defined using Piecewise (this is the default: True).
Warning: k1,k2,k3 refer to Kac labels (components on the basis of simple roots) *)
dynkintokac[{a_,b_,c_}]:= a {3/4,1/2,1/4}+ b {1/2,1,1/2}+c{1/4,1/2,3/4};
weyldot[s_][{la1_,la2_, la3_}]:=weyl[s][{la1,la2,la3}+{1,1,1}]-{1,1,1}; (* warning: acts on weight labels *)
sign[weyl[s_]]:={1,-1,-1,-1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,-1,-1,-1,1}[[s]];
weyl[1][{la1_,la2_, la3_}]:={la1,la2, la3};
weyl[2][{la1_,la2_, la3_}]:={-la1,la1+la2, la3};
weyl[3][{la1_,la2_,la3_}]:={la1+la2,-la2, la2+la3};
weyl[4][{la1_,la2_,la3_}]:={la1,(la2+la3),-la3};
weyl[5][{la1_,la2_,la3_}]:={la2,-(la1+la2),la1+la2+la3};
weyl[6][{la1_,la2_,la3_}]:={-la1,la1+la2+la3,-la3};
weyl[7][{la1_,la2_,la3_}]:={-la1-la2,la1,la2+la3};
weyl[8][{la1_,la2_,la3_}]:={la1+la2,la3,-la2-la3};
weyl[9][{la1_,la2_,la3_}]:={la1+la2+la3, -la2-la3, la2};
weyl[10][{la1_,la2_,la3_}]:={-la2,-la1, la1+la2+la3};
weyl[11][{la1_,la2_,la3_}]:={la2,la3,-( la1+la2+la3)};
weyl[12][{la1_,la2_,la3_}]:={la2+la3, -( la1+la2+la3), la1+la2};
weyl[13][{la1_,la2_,la3_}]:={-la1-la2, ( la1+la2+la3), -(la2+la3)};
weyl[14][{la1_,la2_,la3_}]:={( la1+la2+la3), -la3, -la2};
weyl[15][{la1_,la2_,la3_}]:={-( la1+la2+la3), la1, la2};
weyl[16][{la1_,la2_,la3_}]:={-la2,la2+la3,-( la1+la2+la3)};
weyl[17][{la1_,la2_,la3_}]:={la2+la3,-la3,-la1-la2};
weyl[18][{la1_,la2_,la3_}]:={-la2-la3,-la1,la1+la2};
weyl[19][{la1_,la2_,la3_}]:={la3,-( la1+la2+la3),la1};
weyl[20][{la1_,la2_,la3_}]:={-( la1+la2+la3), la1+la2, -la2};
weyl[21][{la1_,la2_,la3_}]:={la3,-la2-la3, -la1};
weyl[22][{la1_,la2_,la3_}]:={-la2-la3, la2, -la1-la2};
weyl[23][{la1_,la2_,la3_}]:={-la3,- la1-la2, la1};
weyl[24][{la1_,la2_,la3_}]:={-la3, -la2, -la1};

If[Nand@@Map[IntegerQ,dynkintokac[lambda+mu-nu]],Print["Not an admissible intertwiner"];0,
Sum[sign[weyl[v]]sign[weyl[w]]Kostant[dynkintokac[weyldot[v][lambda]+weyldot[w][mu]-nu]],{v,1,4!},{w,1,4!}]],

currentgroup > 4, Print["Kostant partition function not implemented, I continue using the hive method"]; multiplicityHive[lambda,mu,nu]];

(*  Generating Function multiplicityGeneratingFunctionSU3  for L@M@N --> 1.   It is only implemented for SU(3)
It can be used to determine the tensor product of two irreps L@M, or the multiplicity of a triple lambda@mu-> nu*)
(* It is called GSU3 in other notebooks *)

multiplicityGeneratingFunctionSU3[L1_,L2_,M1_,M2_,N1_,N2_]:=(1−L1 L2 M1 M2 N1 N2)/(((1−L1 N2)(1−L1 M2)(1−L2 M1)(1−L2 N1))((1−M2 N1)(1−M1 N2)(1−L1 M1 N1)(1−L2 M2 N2)));

tensorproductGenFunPol[{l1_,l2_},{m1_,m2_}]:= Coefficient[
Expand[Normal[Series[multiplicityGeneratingFunctionSU3[L1,L2,M1,M2,N2,N1],
{L1,0,l1},{L2,0,l2},{M1,0,m1},{M2,0,m2},{N1,0,l1+m1+Min[l2,m2]},{N2,0,l2+m2+Min[l1,m1]}]]],L1^l1 L2^l2 M1^m1 M2^m2 ];
(* warning: we have reversed N1 N2 in the arguments of GSU3 (conjugate) because we are interested in L M --> N  *)

tensorproductGenFun[lambda_,mu_]:= 
Sort@Which[
currentgroup==3,  CoefficientRules[tensorproductGenFunPol[lambda,mu],{N1,N2}]/. {Rule -> List}//Sort,
currentgroup > 3, Print["Generating function not implemented, I continue using the hive method"]; tensorproduct[lambda,mu]];

multiplicityGenFun[lambda_,mu_,nu_]:= Coefficient[tensorproductGenFunPol[lambda,mu], N1^nu[[1]] N2^nu[[2]]];

(* The following family of commands work for SU(3) only. They are copied from the package Sympol$Package.wl *)

HornFull[α_,β_]:=
{{ Max[α[[1]]+β[[3]],α[[2]]+β[[2]],α[[3]]+β[[1]]],α[[1]] +β[[1]] },
{ Max[α[[2]]+β[[3]],α[[3]]+β[[2]]],Min[α[[1]]+β[[2]],α[[2]]+β[[1]]]},
{α[[3]]+β[[3]],Min[α[[1]]+β[[3]],α[[2]]+β[[2]],α[[3]]+β[[1]]]}
};

magicmult[la_,mu_,nuc_]:=
Module[{irr={la,mu,nuc},irrtra,redirrtra,sig,k,magic1, magic2},
irrtra=Transpose[irr]; sig=Map[Total,irrtra];
 k=Abs[sig[[1]]-sig[[2]]]/3; If[Not[IntegerQ[k]],Return[0]];redirrtra=irrtra- If[sig[[1]]>=sig[[2]],{ ConstantArray[k,3],0},{0,  ConstantArray[k,3]}];
{magic1,magic2}=Map[Total,redirrtra]; 
Min[Flatten[{redirrtra,ConstantArray[magic1,3]-redirrtra[[1]]-redirrtra[[2]]}]]+1];(* semi-magic method. The last argument is complex conjugated *)

SchurProdPairMagicMethod[α_,β_]:=Module[{la,mu,horngam,gamma1min,gamma1max,gamma2min,gamma2max,γ3},
la={α[[1]]-α[[2]],α[[2]]-α[[3]]};mu={β[[1]]-β[[2]],β[[2]]-β[[3]]};
horngam=HornFull[α,β];
gamma1min= horngam[[1,1]];gamma1max= horngam[[1,2]];gamma2min= horngam[[2,1]];gamma2max= horngam[[2,2]];
Map[{Take[#,3],Last[#]}& ,DeleteCases[Flatten[Table[
γ3= Total[α]+Total[β] - (γ1+γ2) ;
{γ1,γ2,γ3,Ramp@magicmult[la,mu,Reverse@{γ1-γ2,γ2- γ3}]},
{γ1,gamma1min, gamma1max},{γ2,gamma2min,gamma2max}],1],{_,_,_,0},1]]];

SU3MultiplicityYoung[la_,mu_]:=
Transpose@With[{li=SchurProdPairMagicMethod[{la[[1]]+la[[2]],la[[2]],0},{mu[[1]]+mu[[2]],mu[[2]],0}]},
With[{n=3},{Map[PadRight[#,n-1]&,Map[SUreductionOnlists[#,n]&,Map[First,li]]],Map[Last,li]}]];

SU3MultiplicityDynkin[la_,mu_]:=With[{finpar=SU3MultiplicityYoung[la,mu]},
Transpose[{Map[youngPartitionToHighestWeight[#,3]&,Transpose[finpar][[1]]],Transpose[finpar][[2]]}]];

youngPartitionToHighestWeight[li_,cox_]:=Module[{lifull=PadRight[li,cox]},
Table[lifull[[i]]-lifull[[i+1]],{i,1,cox-1}]];

SUreductionOnlists[lis_, n_] := 
 Which[Length[lis] > n, Print["illegal input"], 
  Length[lis] == n && CountDistinct[lis] == 1, {0}, True, 
  transposeYoungPartition[
   SUreductionTr[transposeYoungPartition[lis], n]]];
   
transposeYoungPartition[li_] := 
 Apply[Plus, 
  Transpose[
   PadRight[Table[ConstantArray[1, li[[s]]], {s, 1, Length[li]}]]], 1];
   (* We have already defined ConjugatePartition! This should be the same! *)

SUreductionTr[lisT_, n_] := 
 Which[First[lisT] > n, 0, First[lisT] == n, 
  SUreductionTr[Drop[lisT, 1], n], True, lisT];
  (* assume that the lists have a length which is not bigger than n. todo *) 
  
SUreductionOnlists[lis_, n_] := 
 Which[Length[lis] > n, Print["illegal input"], 
  Length[lis] == n && CountDistinct[lis] == 1, {0}, True, 
  transposeYoungPartition[
   SUreductionTr[transposeYoungPartition[lis], n]]];
   
   
   (* Calculate multiplicity from the volume function J3 renamed volumeJ *)
   
  volumeJ[{la1_,la2_},{mu1_,mu2_},{nu1_,nu2_}]:=4 (la1+la2+mu1+mu2+nu1+nu2-Abs[la1-la2+mu1-mu2-nu1+nu2]
-2 (({{\[Piecewise], ({{la1+la2-mu1-nu1, la1+2 la2+nu2>=2 mu1+mu2+nu1&&2 la1+la2+mu2>mu1+2 nu1+nu2}, {-la2+mu1+mu2-nu2, la1+2 la2+nu2<2 mu1+mu2+nu1&&la1+mu1+2 mu2>=la2+nu1+2 nu2}, {-la1-mu2+nu1+nu2, 2 la1+la2+mu2<=mu1+2 nu1+nu2&&la1+mu1+2 mu2<la2+nu1+2 nu2}, {0, True}})}}))
-2 (({{\[Piecewise], ({{-la1+mu1+mu2-nu1, 2 la1+la2+nu1<=mu1+2 mu2+nu2&&la1+2 nu1+nu2<la2+2 mu1+mu2}, {la1+la2-mu2-nu2, 2 la1+la2+nu1>mu1+2 mu2+nu2&&la1+2 la2+mu1>=mu2+nu1+2 nu2}, {-la2-mu1+nu1+nu2, la1+2 nu1+nu2>=la2+2 mu1+mu2&&la1+2 la2+mu1<mu2+nu1+2 nu2}, {0, True}})}})));
multiplicityVolume[{la1_,la2_},{mu1_,mu2_},{nu1_,nu2_}]:= 
If[currentgroup =!= 3, Print["This algorithm works only for SU(3). I continue using the hive method"]; 
multiplicityHive[lambda,mu,nu],
volumeJ[{la1,la2},{mu1,mu2},{nu1,nu2}]/24+1];
 
   ]  
(* end of the big With NN *)
End[];(*private*)
(* protect what you want*)
EndPackage[];

Names["obladePackage`*"]

The command setGroupSU[n] should be evaluated at least once (initialization cell), for some value of n. See examples. 
This is not always needed because several functions  (not all : code to be improved !)  adjust automatically the value of n to the number of their arguments.

setGroupSU[3]

Examples for Mathematica

Example with n=3 :  λ={3,1};μ={2,5};

Although this is not its principal interest, the present program can list all the possible ν  (with multiplicity) that appear in the decomposition of the product λ and μ.

setGroupSU[3]

λ={3,4};μ={2,5};

tensorproduct[λ,μ](* The irreducible representations that appear in the decomposition of λ⊗μ into a direct sum, with their multiplicity (the last entry of each sublist) *)

The program recognizes that the current group is SU(3)

currentgroup

The following example assumes that ν is given (it is chosen from the above list).

ν={3,7};

multiplicity[λ,μ,ν](* For :  la  mu -> nu. We knew the result already. *)

The multiplicity being 3, there are 3 oblades associated with this space of intertwiners (ie with this branching rule) :

plotOblades[λ,μ,ν]

plotOblades[λ,μ,ν, 2](* We only display the second one *)

Clear[λ,μ,ν]

Example with n=6:  λ={2,1,1,1,1};μ={1,1,2,1,1};ν={2,2,3,2,1};

setGroupSU[6]

λ={2,1,1,1,1};μ={1,1,2,1,1};

tensorproduct[λ,μ]

currentgroup

The following example assumes that ν is given (it is chosen in the above list).

ν={2,2,3,2,1};

multiplicity[λ,μ,ν](* For :  la  mu -> nu. We knew the result already. *)

The multiplicity being 16, there are 16 oblades associated with this space of intertwiners (ie with this branching rule):

?plotOblades

plotOblades[λ,μ,ν]

plotOblades[λ,μ,ν, 5](* We only display the fifth one *)

νcc=Reverse[ν](* The conjugate of ν *)

?pictographsCoordinates

pictographsCoordinates[λ,μ,ν]
(* Components of the 16 oblades along a basis of 25 fundamental pictographs *)

Dimensions[%]

MatrixRank[%%]-1 (* Dimension of the hive polytope: the triple is therefore generic. *)

Other functions : value of edges. Example with λ={2,1,1,1,1};μ={1,1,2,1,1};ν={2,2,3,2,1}

The following three functions give the list of values of (internal) edges in the directions EW, NE-SW and NW-SE.
Warning: one has to conjugate the last argument in the function call

edgesEWvalues[λ, μ,Reverse[ν]]

edgesNESWvalues[λ, μ,Reverse[ν]]

edgesNWSEvalues[λ, μ,Reverse[ν]]

An application (using  n = 4)

setGroupSU[4]

λ={3,4,2};μ={2,5,3};

tensorproduct[λ,μ]

ν={5,4,7};

{ λ, μ, ν}
multiplicity[λ,μ,ν]

tb=Table[multiplicity[s λ,s μ,s ν],{s,0,7}](* Warning: be patient if you increase the scaling factor too much ! *)

pictographsCoordinates[λ,μ,ν]

deg=MatrixRank[%]-1

pol=InterpolatingPolynomial[tb,s+1]//Expand  (* the Littlewood Richardson polynomial for the chosen triple *)

(pol/. s-> 1 )==multiplicity[λ,μ,ν](* check for :  la  mu -> nu *)

With[{t=8}, {(pol /. s-> t), multiplicity[t λ,t μ,t ν]}](* direct check *)

100 { λ, μ, ν}  (* a triple of large weights *)
pol/. s-> 100   (* Its multiplicity *)

An example with n = 5

setGroupSU[5]

multiplicity[{3,4,3,5},{4,3,5,4},{3,4,4,7}]

We shall not display the 1427 oblades associated with this triple (!),  but we can display for instance the one which is in position 799

plotOblades[{3,4,3,5},{4,3,5,4},{3,4,4,7}, 799]

An example with n = 6

setGroupSU[6]

multiplicity[{3,2,3,2,2},{2,3,2,2,3},{4,5,4,2,5}]

We shall not display the 592 oblades associated with this triple (!),  but we can display for instance the one which is in position 5

plotOblades[{3,2,3,2,2},{2,3,2,2,3},{4,5,4,2,5}, 5]

A non admissible case

setGroupSU[3]

The result of multiplicity is 0

multiplicity[{1,0},{2,1},{3,0}]

Then plotOblades returns an empty list (maybe one should modify this behaviour and make plotOblades return a point)

plotOblades[{1,0},{2,1},{3,0}]

currentgroup

Clear variables

Clear[λ,μ,ν,tb]

User interface (Manipulate) for Mathematica or for the CDF player on a computer : Code and output

Messages (initialization cells) : size limits for the web version of Manipulate

limitlenval=800;
limittimeval=1.5;

Print[Style["Web version: CHANGE The number of inequivalent irreps generated by the program should be smaller than ",12,Bold], Style[limitlenval,12,Bold]]

Print[Style["More time needed: the computation is terminated if it takes more than ",12,Bold],
Style[limittimeval,12,Bold],Style[" seconds",12,Bold]]

limittext="Web version: CHANGE the computation is terminated if it takes more than "1.5`" seconds";

limittime="More time needed: the computation is terminated if it takes more than "1.5`" seconds";

abouttabletext={"A triple of irreducible representations (ω1, ω2, ω3) specified by their higest weights is admissible if ω1 ⊗ ω2 ⊗ ω3 contains ℂ (the scalars)," , "or equivalently if the decomposition of ω1 ⊗ ω2 into a direct sum contains (OverBar[ω3]) (the conjugate of ω3).","Lines of the above table are the components of three highest weights (ω1, ω2; (OverBar[ω3])) in the basis of fundamental weights (Dynkin basis).",
"The triple (ω1, ω2, ω3) should be admissible.","","FILL THE ABOVE TABLE (USE TABS) AND HIT \"Update Components\".",
"If needed, the decomposition of ω1 ⊗ ω2 as a direct sum is obtained by clicking the button \"Tensor product decomposition of ω1 ⊗ ω2 \"; the last entry of each term returned by this command is its multiplicity; warning: the output can be quite long. "
  };

aboutobladetext={"A triplet of multiplicity m is associated with m distinct couplings represented by pictographs.",
 "There are several types of pictographs: BZ-triangles, KT-honeycombs, O-blades. Here we display O-blades (or oblades).",
"Display all pictographs, one at a time, by increasing or descreasing the count index (buttons Next pictograph, Previous pictograph, or the popup menu)."};

notadmissibletext=StringJoin[
"(ω1, ω2, ω3) is not an admissible triple: Modify the components and press \"Update\" \n
To see all possible choices for (OverBar[ω3]), using the displayed values of ω1 and ω2 (the list may be long):\n
 Press the button \"Tensor product decomposition of ω1 ⊗ ω2\" "];

Code for Manipulate (the user should enable dynamic updating)

outputObladeManipulate=With[{fieldSize=3,n=3, limitlen=limitlenval, limittime=limittimeval},
DynamicModule[{rows=3,cols=n-1,
mat=With[{mat12=Table[RandomInteger[7],{2},{n-1}]},Append[mat12,First@RandomChoice[Apply[tensorproduct,mat12]]]],
coupling=1, mul, myRange, myPart,tp=0,tpval},
myRange[s_]:= Range[s];
myPart[expr_,s__]:= Part[expr,s];
mul=Dynamic[Apply[multiplicity,mat]];
Dynamic@Manipulate[
Which[
tp==0,(* Draw pictograph *)
Row[{Invisible[doit],
CheckAbort[If[Length[#]>limitlen, Abort[],#],Print[limittext, limitlen]]&[TimeConstrained[If[Setting@mul==0,notadmissibletext,
TableForm[{PopupMenu[Dynamic[coupling],Range[1,Setting@mul]],Show[plotOblades[myPart[mat,1],myPart[mat,2],myPart[mat,3], coupling],ImageSize-> Full (* Large *)]}]],limittime,Print[limittime, limittime, " seconds"]]]}],

tp==1,(* Calculate the decomposition of ω1 ⊗ ω2 *)
tpval=tensorproduct[myPart[mat,1],myPart[mat,2]],

tp==2,(* a nice SU(6) example *)
tp=0;
mat={{3,2,3,2,2},{2,3,2,2,3},{4,5,4,2,5}};
TableForm[{PopupMenu[Dynamic[coupling],Range[1,Setting@mul]],Show[plotOblades[myPart[mat,1],myPart[mat,2],myPart[mat,3], coupling],ImageSize-> Full]}]
],

Evaluate@With[{makeRow=Function[{rowIndex}, Map[Function[{colIndex},InputField[Dynamic[mat[[rowIndex,colIndex]]],Number,FieldSize->fieldSize]],myRange[cols]],HoldAll]},
Row[{TextGrid[{{"ω1",""},{},{"ω2",""},{},{OverBar["ω3"],""}}],Grid[Map[makeRow,myRange[rows]],Dividers->{{False,False,False},{False,False,Red}}],Invisible["xxxxxx"],
Button[Style["Update components",FontSize->14,Bold, Blue],mat;coupling=1;cols++;cols--;tp=0],
Column[{
Button[Style["Previous pictograph",Bold],If[coupling>1, coupling--; ]],
Button[Style["Next pictograph",Bold], If[coupling<Setting[mul], coupling++; ]]}],
Column[abouttabletext]}]],
Column[{Panel[Style[cols+1,FontSize->14,Bold],Style["SU(n), with n =",Bold, FontSize->14],Left,BaselinePosition->Bottom],Panel[Style[Column[{"Pictographs for admissible triples","of irreducible representations", "of SU(n)","O-blade version"},Center], "Title"]]},Invisible["xxxxxx"]],
Panel[mul,Style["Multiplicity is",FontSize->14,Bold],Left],
Panel[coupling,Style["The count index of this pictograph is",Bold],Left],
Column[aboutobladetext],
Button[Style["SU(n): decrease n",FontSize->14,Bold],If[rows>1,
mat=Map[(Drop[#,-1])&,mat];
coupling=1;
cols--;]],
Button[Style["SU(n): increase n",FontSize->14,Bold],
mat=Transpose[Append[Transpose[mat],Table[0,{3}]]];
coupling=1;
cols++;],
Button[Style["Tensor product decomposition of ω1 ⊗ ω2 ",FontSize->14,Bold, Blue],coupling=1; 
tp=1;cols++;cols--;],
Button[Style["An SU(6) example with multiplicity 592",FontSize->10,Bold, Brown,TextAlignment->Left],coupling=1;cols=6-1;tp=2; cols++;cols--; ],
Button[Style["Reset random SU(3) example",FontSize->10,Bold, Brown,TextAlignment->Left],
cols=2;mat=With[{mat12=Table[RandomInteger[7],{2},{3-1}]},Append[mat12,First@RandomChoice[Apply[tensorproduct,mat12]]]];
coupling=1;cols++;cols--;tp=0],
{{doit,0,""},{0},ControlType->None},
TrackedSymbols:> {doit}(* ,ContentSize->{650,650}*)],
SaveDefinitions->True, SynchronousUpdating-> False]];

Manipulate window

CellPrint[ExpressionCell[outputObladeManipulate,CellLabel->""]]