Ticket #5794: trac_5794-continued.patch

File trac_5794-continued.patch, 23.9 KB (added by bump, 13 years ago)
  • sage/combinat/root_system/all.py

    # HG changeset patch
    # User Daniel Bump <bump@match.stanford.edu>
    # Date 1242139111 25200
    # Node ID dabb54d76d8fb9ace494b1a1ca8b48c298d7df94
    # Parent  73f7ba09b8f420d93f86b8caa8ea16ad5c68aa14
    [mq]: trac_5794-continued.patch
    
    diff --git a/sage/combinat/root_system/all.py b/sage/combinat/root_system/all.py
    a b  
    44from coxeter_matrix import coxeter_matrix
    55from root_system import RootSystem, WeylDim
    66from weyl_group import WeylGroup, WeylGroupElement
    7 from weyl_characters import WeylCharacter, WeylCharacterRing, branch_weyl_character, WeightRing, WeightRingElement
     7from weyl_characters import WeylCharacter, WeylCharacterRing, branch_weyl_character, WeightRing, WeightRingElement, branching_rule_from_plethysm
  • sage/combinat/root_system/weyl_characters.py

    diff --git a/sage/combinat/root_system/weyl_characters.py b/sage/combinat/root_system/weyl_characters.py
    a b  
    1919from sage.combinat.root_system.root_system import RootSystem
    2020from sage.combinat.root_system.cartan_type import CartanType
    2121from sage.modules.free_module import VectorSpace
     22from sage.modules.free_module_element import vector
    2223from sage.structure.element import is_Element
     24from sage.matrix.constructor import matrix
    2325from sage.rings.all import ZZ, QQ
    2426from sage.misc.misc import repr_lincomb
    2527from sage.misc.functional import is_odd, is_even
     
    347349        """
    348350        return [[k,m] for k,m in self._hdict.iteritems()]
    349351
     352    def is_irreducible(self):
     353        """
     354        Returns True if self is an irreducible character.
     355
     356        EXAMPLES::
     357       
     358            sage: B3 = WeylCharacterRing(['B',3])
     359            sage: [B3(x).is_irreducible() for x in B3.fundamental_weights()]
     360             [True, True, True]
     361            sage: sum(B3(x) for x in B3.fundamental_weights()).is_irreducible()
     362             False
     363        """
     364        h = self.hlist()
     365        return len(h) is 1 and h[0][1] is 1
     366
     367    def symmetric_square(self):
     368        """
     369        Returns the symmetric square of the character.
     370
     371        EXAMPLES::
     372       
     373            sage: A2 = WeylCharacterRing("A2",style="coroots")
     374            sage: A2(1,0).symmetric_square()
     375             A2(2,0)
     376        """
     377
     378        cmlist = self.mlist()
     379        mdict = {}
     380        for j in range(len(cmlist)):
     381            for i in range(j+1):
     382                t = cmlist[i][0]+cmlist[j][0]
     383                if i < j:
     384                    coef = cmlist[i][1]*cmlist[j][1]
     385                else:
     386                    coef = cmlist[i][1]*(cmlist[i][1]+1)/2
     387                if t in mdict:
     388                    mdict[t] += coef
     389                else:
     390                    mdict[t] = coef
     391        hdict = self._parent.char_from_weights(mdict)
     392        return WeylCharacter(self._parent, hdict, mdict)
     393
     394    def exterior_square(self):
     395        """
     396        Returns the exterior square of the character.
     397
     398        EXAMPLES::
     399       
     400            sage: A2 = WeylCharacterRing("A2",style="coroots")
     401            sage: A2(1,0).exterior_square()
     402             A2(0,1)
     403        """
     404        cmlist = self.mlist()
     405        mdict = {}
     406        for j in range(len(cmlist)):
     407            for i in range(j+1):
     408                t = cmlist[i][0]+cmlist[j][0]
     409                if i < j:
     410                    coef = cmlist[i][1]*cmlist[j][1]
     411                else:
     412                    coef = cmlist[i][1]*(cmlist[i][1]-1)/2
     413                if t in mdict:
     414                    mdict[t] += coef
     415                else:
     416                    mdict[t] = coef
     417        hdict = self._parent.char_from_weights(mdict)
     418        return WeylCharacter(self._parent, hdict, mdict)
     419
     420    def frobenius_schur_indicator(self):
     421        """
     422        Returns:
     423
     424        1 if the representation is real (orthogonal)
     425        -1 if the representation is quaternionic (symplectic)   
     426        0 if the representation is complex (not self dual)
     427
     428        The Frobenius-Schur indicator of a character 'chi'
     429        of a compact group G is the Haar integral over the
     430        group of 'chi(g^2)'. Its value is 1,-1 or 0. This
     431        method computes it for irreducible characters of
     432        compact Lie groups by checking whether the symmetric
     433        and exterior square characters contain the trivial
     434        character.
     435
     436        EXAMPLES::
     437       
     438            sage: B2 = WeylCharacterRing("B2",style="coroots")
     439            sage: B2(1,0).frobenius_schur_indicator()
     440             1
     441            sage: B2(0,1).frobenius_schur_indicator()
     442             -1
     443        """
     444        if not self.is_irreducible():
     445            raise ValueError, "Frobenius-Schur indicator is only valid for irreducible characters"
     446        z = self._parent.space().zero()
     447        if z in self.symmetric_square()._hdict:
     448            return 1
     449        elif z in self.exterior_square()._hdict:
     450            return -1
     451        else:
     452            return 0
     453
    350454    def mlist(self):
    351455        """
    352456        Returns a list of weights in self with their multiplicities.
     
    895999    King, Branching rules for classical Lie groups using tensor and
    8961000    spinor methods. J. Phys. A 8 (1975), 429-449, Howe, Tan and
    8971001    Willenbring, Stable branching rules for classical symmetric pairs,
    898     Trans. Amer. Math. Soc. 357 (2005), no. 4, 1601-1626 and McKay and
     1002    Trans. Amer. Math. Soc. 357 (2005), no. 4, 1601-1626, McKay and
    8991003    Patera, Tables of Dimensions, Indices and Branching Rules for
    900     Representations of Simple Lie Algebras (Marcel Dekker, 1981).
     1004    Representations of Simple Lie Algebras (Marcel Dekker, 1981),
     1005    and Fauser, Jarvis, King and Wybourne, New branching rules induced
     1006    by plethysm. J. Phys. A 39 (2006), no. 11, 2611--2655.
    9011007   
    9021008    INPUT:
    9031009
     
    9261032    We will list give predefined rules that cover most cases where the
    9271033    branching rule is to a maximal subgroup. For convenience, we
    9281034    also give some branching rules to subgroups that are not maximal.
    929     For example, a Levi subgroup may or may not be maximal, but you
    930     can usually branch to it.
     1035    For example, a Levi subgroup may or may not be maximal.
    9311036   
    9321037    LEVI TYPE. These can be read off from the Dynkin diagram. If
    9331038    removing a node from the Dynkin diagram produces another Dynkin
     
    10021107    (e.g. B2xB3 -> D6) cannot be explained this way but for
    10031108    uniformity is implemented under rule="extended".
    10041109
    1005     You can therefore get any branching rule
     1110    Using rule="extended" you can get any branching rule
     1111    SO(n) => SO(a) x SO(b) x SO(c) x ... where n = a+b+c+ ...
     1112    Sp(2n) => Sp(2a) x Sp(2b) x Sp(2c) x ... where n = a+b+c+ ...
     1113    where O(a) = ['D',r] (a=2r) or ['B',r] (a=2r+1)
     1114    and Sp(2r)=['C',r].
    10061115
    1007     O(n) => O(a)xO(b)xO(c)x ... where n = a+b+c+ ...
    1008     Sp(2n) => Sp(2a)xSp(2b)xSp(2c)x ... where n = a+b+c+ ...
    1009 
    1010     where O(a) = ['D',r] (a=2r) or ['B',r] (a=2r+1)
    1011     and Sp(2r)=['C',r] using rule="extended".   
    1012 
    1013     TENSOR: The branching rule
     1116    TENSOR: There are branching rules
    10141117
    10151118    ['A', rs-1] => ['A',r-1] x ['A',s-1]
    1016 
    1017     corresponding to the tensor product homomorphism GL(r) x GL(s) -> GL(rs)
    1018     is implemented using rule="tensor". There are also branching rules
    1019     for types B,C and D. These are:
    1020 
    10211119    ['B',2rs+r+s] => ['B',r] x ['B',s]
    10221120    ['D',2rs+s] => ['B',r] x ['D',s]
    10231121    ['D',2rs] => ['D',r] x ['D',s]
     
    10251123    ['C',2rs+s] => ['B',r] x ['C',s]
    10261124    ['C',2rs] => ['C',r] x ['D',s].
    10271125
    1028     These are not implemented yet though you may obtain them using
    1029     handwritten rules.
     1126    corresponding to the tensor product homomorphism. For type
     1127    A, the homomorphism is GL(r) x GL(s) -> GL(rs). For the
     1128    classical types, the relevant fact is that if V,W are
     1129    orthogonal or symplectic spaces, that is, spaces endowed
     1130    with symmetric or skew-symmetric bilinear forms, then V
     1131    tensor W is also an orthogonal space (if V and W are both
     1132    orthogonal or both symplectic) or symplectic (if one of
     1133    V and W is orthogonal and the other symplectic).
     1134   
     1135    The corresponding branching rules are obtained using rule="tensor".
     1136
    10301137
    10311138    SYMMETRIC POWER: The k-th symmetric and exterior power homomorphisms
    10321139    map GL(n) --> GL(binomial(n+k-1,k)) and GL(binomial(n,k)). The
     
    10381145    ['B',r] => A1
    10391146    ['C',r] => A1
    10401147
    1041     and these may be obtained using the rule "symmetric power".
     1148    and these may be obtained using the rule "symmetric_power".
    10421149
    10431150    MISCELLANEOUS: Use rule="miscellaneous" for the following rule,
    10441151    which does not fit into the above framework.
    10451152   
    10461153    B3 => G2
    10471154   
     1155    BRANCHING RULES FROM PLETHYSMS
     1156
     1157    Nearly all branching rules G => H where G is of type A,B,C or D
     1158    are covered by the preceding rules. The function
     1159    branching_rules_from_plethysm covers the remaining cases.
     1160
    10481161    ISOMORPHIC TYPE: Although not usually referred to as a branching
    10491162    rule, the effects of the accidental isomorphisms may be handled
    10501163    using rule="isomorphic"
     
    12191332         A2xA1(1,0,2) + A2xA1(0,2,0),
    12201333         A2xA1(0,1,1)]
    12211334
    1222     Although tensor products for the classical types are not implemented
    1223     as hard-coded rules, you may still write your own rule for any given
    1224     case. For example the tensor product homomorphism SO(3)xSO(3)-->SO(9)
    1225     gives a branching rule B4=>B1xB1 which may be computed as follows:
    1226 
    12271335        sage: B4=WeylCharacterRing("B4",style="coroots")
    12281336        sage: B1xB1=WeylCharacterRing("B1xB1",style="coroots")
    1229         sage: [B4(hwv).branch(B1xB1, rule=lambda x : [x[0]-x[2]+x[3],x[0]+x[1]+x[2]]) for hwv in B4.fundamental_weights()]
     1337        sage: [B4(f).branch(B1xB1,rule="tensor") for f in B4.fundamental_weights()]
    12301338        [B1xB1(2,2),
    12311339         B1xB1(0,2) + B1xB1(2,0) + B1xB1(2,4) + B1xB1(4,2),
    12321340         B1xB1(0,2) + B1xB1(0,6) + B1xB1(2,0) + B1xB1(2,2) + B1xB1(2,4) + B1xB1(4,2) + B1xB1(4,4) + B1xB1(6,0),
    12331341         B1xB1(1,3) + B1xB1(3,1)]
    12341342
     1343        sage: D4=WeylCharacterRing("D4",style="coroots")
     1344        sage: C2xC1=WeylCharacterRing("C2xC1",style="coroots")
     1345        sage: [D4(f).branch(C2xC1,rule="tensor") for f in D4.fundamental_weights()]
     1346        [C2xC1(1,0,1),
     1347         C2xC1(0,0,2) + C2xC1(0,1,2) + C2xC1(2,0,0),
     1348         C2xC1(1,0,1),
     1349         C2xC1(0,0,2) + C2xC1(0,1,0)]
     1350
     1351        sage: C3=WeylCharacterRing("C3",style="coroots")
     1352        sage: B1xC1=WeylCharacterRing("B1xC1",style="coroots")
     1353        sage: [C3(f).branch(B1xC1,rule="tensor") for f in C3.fundamental_weights()]
     1354        [B1xC1(2,1), B1xC1(2,2) + B1xC1(4,0), B1xC1(0,3) + B1xC1(4,1)]
     1355
    12351356    EXAMPLES: (Symmetric Power)
    12361357
    12371358        sage: A1=WeylCharacterRing("A1",style="coroots")
     
    12881409        sage: [D2(fw).branch(A1xA1,rule="isomorphic") for fw in D2.fundamental_weights()]
    12891410        [A1xA1(1,0), A1xA1(0,1)]
    12901411
     1412    EXAMPLES: (Branching rules from plethysms)
     1413
     1414    This is a general rule that includes any branching rule
     1415    from types A,B,C or D as a special case. Thus it could be
     1416    used in place of the above rules and would give the same
     1417    results.
     1418
     1419    We consider a homomorphism H --> G where G is one of
     1420    SL(r+1), SO(2r+1), Sp(2r) or SO(2r). The function
     1421    branching_rule_from_plethysm produces the corresponding
     1422    branching rule. The main ingredient is the character
     1423    chi of the representation of H that is the homomorphism
     1424    to GL(r+1), GL(2r+1) or GL(2r).
     1425
     1426    This rule is so powerful that it contains the other
     1427    rules implemented above as special cases. First let
     1428    us consider the symmetric fifth power representation
     1429    of SL(2).
     1430
     1431        sage: A1=WeylCharacterRing("A1",style="coroots")
     1432        sage: chi=A1([5])
     1433        sage: chi.degree()
     1434         6
     1435        sage: chi.frobenius_schur_indicator()
     1436        -1
     1437
     1438    This confirms that the character has degree 6 and
     1439    is symplectic, so it corresponds to a homomorphism
     1440    SL(2) --> Sp(6), and there is a corresponding
     1441    branching rule C3 => A1.
     1442
     1443        sage: C3 = WeylCharacterRing("C3",style="coroots")
     1444        sage: sym5rule = branching_rule_from_plethysm(chi,"C3")
     1445        sage: [C3(hwv).branch(A1,rule=sym5rule) for hwv in C3.fundamental_weights()]
     1446        [A1(5), A1(4) + A1(8), A1(3) + A1(9)]
     1447
     1448    This is identical to the results we would obtain using
     1449    rule="symmetric_power". The next example gives a branching
     1450    not available by other standard rules.
     1451
     1452        sage: G2 = WeylCharacterRing("G2",style="coroots")
     1453        sage: D7 = WeylCharacterRing("D7",style="coroots")
     1454        sage: ad=G2(0,1); ad.degree(); ad.frobenius_schur_indicator()
     1455         14
     1456         1
     1457        sage: spin = D7(0,0,0,0,0,1,0); spin.degree()
     1458         64
     1459        sage: spin.branch(G2, rule=branching_rule_from_plethysm(ad, "D7"))
     1460         G2(1,1)
     1461
     1462    We have confirmed that the adjoint representation of G2
     1463    gives a homomorphism into SO(14), and that the pullback
     1464    of the one of the two 64 dimensional spin representations
     1465    to SO(14) is an irreducible representation of G2.
     1466
    12911467    BRANCHING FROM A REDUCIBLE ROOT SYSTEM
    12921468
    12931469    If you are branching from a reducible root system, the rule is
     
    13101486
    13111487    WRITING YOUR OWN RULES
    13121488
    1313     The rules that are already coded in SAGE are extensive but
    1314     not exhaustive. If you need a rule that is not implemented,
    1315     you may write it yourself.
    1316 
    13171489    Suppose you want to branch from a group G to a subgroup H.
    13181490    Arrange the embedding so that a Cartan subalgebra U of H is
    13191491    contained in a Cartan subalgebra T of G. There is thus
     
    13951567    s = Stype.rank()
    13961568    rdim = Rtype.root_system().ambient_space().dimension()
    13971569    sdim = Stype.root_system().ambient_space().dimension()
     1570    if Stype.is_reducible():
     1571        stypes = Stype.component_types()
    13981572    if rule == "default":
    13991573        return lambda x : x
    14001574    elif rule == "levi":
     
    14021576            raise ValueError, "Incompatible ranks"
    14031577        if Rtype[0] == 'A':
    14041578            if Stype.is_reducible():
    1405                 if all(ct[0]=='A' for ct in Stype.component_types()) \
     1579                if all(ct[0]=='A' for ct in stypes) \
    14061580                       and rdim == sdim:
    14071581                    return lambda x : x
    14081582                else:
     
    14171591                    return  lambda x : x
    14181592                elif Stype[0] == Rtype[0]:
    14191593                    return lambda x : list(x)[1:]
    1420             elif Stype.component_types()[-1][0] == Rtype[0] and all(t[0] == 'A' for t in Stype.component_types()[:-1]):
     1594            elif stypes[-1][0] == Rtype[0] and all(t[0] == 'A' for t in stypes[:-1]):
    14211595                return lambda x : x
    14221596            else:
    14231597                raise ValueError, "Rule not found"
     
    14721646            raise ValueError, "Rule not found"
    14731647    elif rule == "extended":
    14741648        if Stype.is_reducible():
    1475             if Rtype[0] in ['B','D'] and all(t[0] in ['B','D'] for t in Stype.component_types()):
     1649            if Rtype[0] in ['B','D'] and all(t[0] in ['B','D'] for t in stypes):
    14761650                if Rtype[0] == 'D':
    14771651                    rdeg = 2*r
    14781652                else:
    14791653                    rdeg = 2*r+1
    14801654                sdeg = 0
    1481                 for t in Stype.component_types():
     1655                for t in stypes:
    14821656                    if t[0] == 'D':
    14831657                        sdeg += 2*t[1]
    14841658                    else:
     
    14881662                else:
    14891663                    raise ValueError, "Rule not found"
    14901664            elif Rtype[0] == 'C'  and s == r:
    1491                 if all(t[0] == Rtype[0] for t in Stype.component_types()):
     1665                if all(t[0] == Rtype[0] for t in stypes):
    14921666                    return lambda x : x
    14931667            elif Rtype[0] == 'G' and s == r:
    1494                 if all(t[0] == 'A' and t[1] == 1 for t in Stype.component_types()):
     1668                if all(t[0] == 'A' and t[1] == 1 for t in stypes):
    14951669                    return lambda x : [(x[1]-x[2])/2,-(x[1]-x[2])/2, x[0]/2, -x[0]/2]
    14961670            else:
    14971671                raise ValueError, "Rule not found"
     
    15081682    elif rule == "isomorphic":
    15091683        if r != s:
    15101684            raise ValueError, "Incompatible ranks"
    1511         if Rtype[0] == 'B' and r == 2 and Stype[0] == 'C':
     1685        if Rtype == Stype:
     1686            return lambda x : x
     1687        elif Rtype[0] == 'B' and r == 2 and Stype[0] == 'C':
    15121688            def rule(x) : [x1, x2] = x; return [x1+x2, x1-x2]
    15131689            return rule
    1514         if Rtype[0] == 'B' and r == 1 and Stype[0] == 'A':
     1690        elif Rtype[0] == 'B' and r == 1 and Stype[0] == 'A':
    15151691            return lambda x : [x[0],-x[0]]
    15161692        elif Rtype[0] == 'C' and r == 2 and Stype[0] == 'B':
    15171693            def rule(x) : [x1, x2] = x; return [(x1+x2)/2, (x1-x2)/2]
    15181694            return rule
     1695        elif Rtype[0] == 'C' and r == 1 and Stype[0] == 'A':
     1696            return lambda x : [x[0]/2,-x[0]/2]
    15191697        elif Rtype[0] == 'A' and r == 3 and Stype[0] == 'D':
    15201698            def rule(x): [x1, x2, x3, x4] = x; return [(x1+x2-x3-x4)/2, (x1-x2+x3-x4)/2, (x1-x2-x3+x4)/2]
    15211699            return rule
    15221700        elif Rtype[0] == 'D' and r == 2 and Stype.is_reducible() and \
    1523                  all(t[0] == 'A' for t in Stype.component_types()):
     1701                 all(t[0] == 'A' for t in stypes):
    15241702            def rule(x): [t1, t2] = x; return [(t1-t2)/2, -(t1-t2)/2, (t1+t2)/2, -(t1+t2)/2]
    15251703            return rule
    15261704        elif Rtype[0] == 'D' and r == 3 and Stype[0] == 'A':
     
    15281706            return rule
    15291707        else:
    15301708            raise ValueError, "Rule not found"
    1531     elif rule == "tensor":
     1709    elif rule == "tensor" or rule == "tensor-debug":
    15321710        if not Stype.is_reducible():
    15331711            raise ValueError, "Tensor product requires more than one factor"           
    1534         if len(Stype.component_types()) is not 2:
     1712        if len(stypes) is not 2:
    15351713            raise ValueError, "Not implemented"
     1714        if Rtype[0] is 'A':
     1715            nr = Rtype[1]+1
     1716        elif Rtype[0] is 'B':
     1717            nr = 2*Rtype[1]+1
     1718        elif Rtype[0] in ['C', 'D']:
     1719            nr = 2*Rtype[1]
     1720        else:
     1721            raise ValueError, "Rule not found"
     1722        [s1, s2] = [stypes[i][1] for i in range(2)]
     1723        ns = [s1, s2]
     1724        for i in range(2):
     1725            if stypes[i][0] is 'A':
     1726                ns[i] = ns[i]+1
     1727            if stypes[i][0] is 'B':
     1728                ns[i] = 2*ns[i]+1
     1729            if stypes[i][0] in ['C','D']:
     1730                ns[i] = 2*ns[i]
     1731        if nr != ns[0]*ns[1]:
     1732            raise ValueError, "Ranks don't agree with tensor product"
    15361733        if Rtype[0] == 'A':
    1537             if all(t[0] == 'A' for t in Stype.component_types()):
    1538                 [s1,s2] = [Stype.component_types()[i][1]+1 for i in range(2)]       
    1539                 if not s1*s2 == r+1:
    1540                     raise ValueError, "Ranks don't agree with tensor product"
     1734            if all(t[0] == 'A' for t in stypes):
    15411735                def rule(x):
    1542                     ret = [sum(x[i*s2:(i+1)*s2]) for i in range(s1)]
    1543                     ret.extend([sum(x[s2*j+i] for j in range(s1)) for i in range(s2)])
     1736                    ret = [sum(x[i*ns[1]:(i+1)*ns[1]]) for i in range(ns[0])]
     1737                    ret.extend([sum(x[ns[1]*j+i] for j in range(ns[0])) for i in range(ns[1])])
    15441738                    return ret
    15451739                return rule
    15461740            else:
    15471741                raise ValueError, "Rule not found"
    1548         else:
    1549             raise ValueError, "Not implemented"
     1742        elif Rtype[0] == 'B':
     1743            if not all(t[0] == 'B' for t in stypes):
     1744                raise ValueError, "Rule not found"
     1745        elif Rtype[0] == 'C':
     1746            if stypes[0][0] in ['B','D'] and stypes[1][0] is 'C':
     1747                pass
     1748            elif stypes[1][0] in ['B','D'] and stypes[0][0] is 'C':
     1749                pass
     1750            else:
     1751                raise ValueError, "Rule not found"
     1752        elif Rtype[0] == 'D':
     1753            if stypes[0][0] in ['B','D'] and stypes[1][0] is 'D':
     1754                pass
     1755            elif stypes[1][0] is 'B' and stypes[0][0] is 'D':
     1756                pass
     1757            elif stypes[1][0] is 'C' and stypes[0][0] is 'C':
     1758                pass
     1759            else:
     1760                raise ValueError, "Rule not found"
     1761        rows = []
     1762        for i in range(s1):
     1763            for j in range(s2):
     1764                nextrow = (s1+s2)*[0]
     1765                nextrow[i] = 1
     1766                nextrow[s1+j] = 1
     1767                rows.append(nextrow)
     1768        if stypes[1][0] == 'B':
     1769            for i in range(s1):
     1770                nextrow = (s1+s2)*[0]
     1771                nextrow[i] = 1
     1772                rows.append(nextrow)
     1773        for i in range(s1):
     1774            for j in range(s2):
     1775                nextrow = (s1+s2)*[0]
     1776                nextrow[i] = 1
     1777                nextrow[s1+j] = -1
     1778                rows.append(nextrow)
     1779        if stypes[0][0] == 'B':
     1780            for j in range(s2):
     1781                nextrow = (s1+s2)*[0]
     1782                nextrow[s1+j] = 1
     1783                rows.append(nextrow)
     1784        mat = matrix(rows).transpose()
     1785        if rule == "tensor-debug":
     1786            print mat
     1787        return lambda x : tuple(mat*vector(x))
    15501788    elif rule == "symmetric_power":
    15511789        if Stype[0] == 'A' and s == 1:
    15521790            if Rtype[0] == 'B':
     
    15701808            raise ValueError, "Rule not found"
    15711809
    15721810
     1811def branching_rule_from_plethysm(chi, cartan_type, return_matrix = False):
     1812    """
     1813    INPUT:
     1814
     1815    - ``chi`` - the character of an irreducible representation pi of a group G
     1816    - ``cartan_type`` - a classical Cartan type (A,B,C or D).
     1817
     1818    It is assumed that the image of the irreducible representation pi
     1819    naturally has its image in the group G.
     1820
     1821    Returns a branching rule for this plethysm.
     1822
     1823    EXAMPLE:
     1824
     1825    The adjoint representation SL(3) --> GL(8) factors
     1826    through SO(8). The branching rule in question will
     1827    describe how representations of SO(8) composed with
     1828    this homomorphism decompose into irreducible characters
     1829    of SL(3).
     1830
     1831        sage: A2 = WeylCharacterRing("A2")
     1832        sage: A2 = WeylCharacterRing("A2", style="coroots")
     1833        sage: ad = A2(1,1)
     1834        sage: ad.degree()
     1835         8
     1836        sage: ad.frobenius_schur_indicator()
     1837         1
     1838
     1839    This confirms that ad has degree 8 and is orthogonal,
     1840    hence factors through SO(8)=D4.
     1841
     1842        sage: br = branching_rule_from_plethysm(ad,"D4")
     1843        sage: D4 = WeylCharacterRing("D4")
     1844        sage: [D4(f).branch(A2,rule = br) for f in D4.fundamental_weights()]
     1845         [A2(1,1), A2(1,1) + A2(0,3) + A2(3,0), A2(1,1), A2(1,1)]
     1846    """
     1847    ct = CartanType(cartan_type)
     1848    if ct[0] not in ["A","B","C","D"]:
     1849        raise ValueError, "not implemented for type %s"%ct[0]
     1850    if ct[0] is "A":
     1851        ret = []
     1852        for [v,n] in chi.mlist():
     1853            ret.extend(n*[v.to_vector()])
     1854        M = matrix(ret).transpose()
     1855        if len(M.columns()) != ct[1] + 1:
     1856            raise ValueError, "representation has wrong degree for type %s"%ct.__repr__()
     1857        return lambda x : tuple(M*vector(x))
     1858    if ct[0] in ["B","D"]:
     1859        if chi.frobenius_schur_indicator() != 1:
     1860            raise ValueError, "character is not orthogonal"
     1861    if ct[0] is "C":
     1862        if chi.frobenius_schur_indicator() != -1:
     1863            raise ValueError, "character is not symplectic"
     1864    if ct[0] is "B":
     1865        if is_even(chi.degree()):
     1866            raise ValueError, "degree is not odd"
     1867    if ct[0] is ["C","D"]:
     1868        if is_odd(chi.degree()):
     1869            raise ValueError, "degree is not even"
     1870    ret = []
     1871    for [v,n] in chi.mlist():
     1872        v = v.to_vector()
     1873        if all(x==0 for x in v):
     1874            if ct[0] is "B":
     1875                n = (n-1)/2
     1876            else:
     1877                n = n/2
     1878        elif [x for x in v if x !=0][0] < 0:
     1879            continue
     1880        ret.extend(n*[v])
     1881    M = matrix(ret).transpose()           
     1882    if len(M.columns()) != ct.root_system().ambient_space().dimension():
     1883        raise ValueError, "representation has wrong degree for type %s"%ct.__repr__()
     1884    if return_matrix:
     1885        return M
     1886    else:
     1887        return lambda x : tuple(M*vector(x))
     1888
     1889
    15731890class WeightRingElement(AlgebraElement):
    15741891    """
    15751892    A class for weights, and linear combinations of weights. See