1InstallMethod( ProjectionMatrix,
2"for two semisimple Lie algebras", true, [ IsLieAlgebra, IsLieAlgebra ], 0,
3
4function( L, K )
5
6    local H, HK, N, K0, hK, hL, b;
7
8    H:= CartanSubalgebra(L);
9    HK:= Intersection( H, K );
10
11    N:= LieNormalizer( K, HK );
12    if N <> HK then
13       Error("The Cartan subalgebras do not match.");
14    fi;
15
16    if HasCartanSubalgebra(K) and CartanSubalgebra(K) <> HK then
17       K0:= Subalgebra( L, BasisVectors( Basis(K) ), "basis" );
18       SetCartanSubalgebra(K0,HK);
19    else
20       K0:= K;
21    fi;
22
23    hK:= CanonicalGenerators( RootSystem(K0) )[3];
24    hL:= CanonicalGenerators( RootSystem(L) )[3];
25
26    b:= Basis( Subspace( L, hL ), hL );
27
28    return List( hK, x -> Coefficients( b, x ) );
29
30
31end );
32
33
34InstallOtherMethod( ProjectionMatrix,
35"for two semisimple Lie algebras", true, [ IsLieAlgebra, IsLieAlgebra, IsList ], 0,
36
37function( L, K, cc )
38
39    local H, HK, N, K0, hK, hL, b;
40
41    H:= CartanSubalgebra(L);
42    HK:= Intersection( H, K );
43
44    N:= LieNormalizer( K, HK );
45    if N <> HK then
46       Error("The Cartan subalgebras do not match.");
47    fi;
48
49    if not ForAll( cc, x -> x in H ) then
50       Error("The toral elements do noy lie in the Cartan.");
51    fi;
52
53    if not HasCartanSubalgebra(K) then
54       K0:= K;
55       SetCartanSubalgebra(K0,HK);
56    elif  CartanSubalgebra(K) <> HK then
57       K0:= Subalgebra( L, BasisVectors( Basis(K) ), "basis" );
58       SetCartanSubalgebra(K0,HK);
59    else
60       K0:= K;
61    fi;
62
63    hK:= Concatenation( CanonicalGenerators( RootSystem(K0) )[3], cc );
64    hL:= CanonicalGenerators( RootSystem(L) )[3];
65
66    b:= Basis( Subspace( L, hL ), hL );
67
68    return List( hK, x -> Coefficients( b, x ) );
69
70
71end );
72
73
74
75
76InstallMethod( Branching,
77"for two semisimple Lie algebras and a weight", true,
78[ IsLieAlgebra, IsLieAlgebra, IsList ], 0,
79
80function( L, K, wt )
81
82
83    local P, ch, R, W, w0, m0, i, o, mu, pos, S, sim, b, ord, ww, mm, mult,
84          w, isg;
85
86
87    P:= ProjectionMatrix(L,K);
88    ch:= DominantCharacter( L, wt );
89    R:= RootSystem(L);
90    W:= WeylGroup(R);
91    w0:= [ ]; m0:= [ ];
92    for i in [1..Length(ch[1])] do
93        o:= WeylOrbitIterator( W, ch[1][i] );
94        while not IsDoneIterator(o) do
95             mu:= P*NextIterator( o );
96             if ForAll( mu, u -> u >= 0 ) then
97                pos:= Position( w0, mu );
98                if pos = fail then
99                   Add( w0, mu );
100                   Add( m0, ch[2][i] );
101                else
102                   m0[pos]:= m0[pos]+ch[2][i];
103                fi;
104             fi;
105        od;
106    od;
107
108    S:= RootSystem(K);
109    sim:= SimpleRootsAsWeights(S);
110    b:= Basis( VectorSpace( Rationals, sim ), sim );
111
112    ord:= function( mu, nu )
113
114        # true if mu < nu
115
116        local cf;
117
118        if mu = nu then return fail; fi;
119        cf:= Coefficients( b, nu-mu );
120        if ForAll( cf, IsInt ) then
121           if ForAll( cf, x -> x >= 0 ) then
122              return true;
123           elif ForAll( cf, x -> x <= 0 ) then
124              return false;
125           fi;
126        fi;
127
128        return fail;
129
130    end;
131
132    ww:=[ ]; mm:= [ ];
133
134    while Length(w0) > 0 do
135
136        # look for a maximum in w0...
137
138        w:= w0[1];
139        mult:= m0[1];
140        for i in [1..Length(w0)] do
141            isg:= ord( w, w0[i] );
142            if isg <> fail and isg then
143                 w:= w0[i]; mult:= m0[i];
144            fi;
145        od;
146
147        Add( ww, w ); Add( mm, mult );
148        ch:= DominantCharacter( K, w );
149        for i in [1..Length(ch[1])] do
150            pos:= Position( w0, ch[1][i] );
151            m0[pos]:= m0[pos] - mult*ch[2][i];
152        od;
153        for i in [1..Length(w0)] do
154            if m0[i] = 0 then
155               Unbind( m0[i] );
156               Unbind( w0[i] );
157            fi;
158        od;
159        m0:= Filtered( m0, x -> IsBound(x) );
160        w0:= Filtered( w0, x -> IsBound(x) );
161    od;
162
163    return [ww,mm];
164
165end );
166
167
168InstallOtherMethod( Branching,
169"for two semisimple Lie algebras and a weight", true,
170[ IsLieAlgebra, IsLieAlgebra, IsList, IsList ], 0,
171
172function( L, K, cc, wt )
173
174
175    local P, ch, R, W, w0, m0, i, o, mu, pos, S, sim, b, ord, ww, mm, mult,
176          w, isg, semsimrk;
177
178
179
180    P:= ProjectionMatrix( L, K, cc );
181    semsimrk:= Dimension( CartanSubalgebra(K) );
182
183    ch:= DominantCharacter( L, wt );
184    R:= RootSystem(L);
185    W:= WeylGroup(R);
186    w0:= [ ]; m0:= [ ];
187    for i in [1..Length(ch[1])] do
188        o:= WeylOrbitIterator( W, ch[1][i] );
189        while not IsDoneIterator(o) do
190             mu:= P*NextIterator( o );
191             if ForAll( [1..semsimrk], k -> mu[k] >= 0 ) then
192                pos:= Position( w0, mu );
193                if pos = fail then
194                   Add( w0, mu );
195                   Add( m0, ch[2][i] );
196                else
197                   m0[pos]:= m0[pos]+ch[2][i];
198                fi;
199             fi;
200        od;
201    od;
202
203    S:= RootSystem(K);
204    sim:= SimpleRootsAsWeights(S);
205    b:= Basis( VectorSpace( Rationals, sim ), sim );
206
207    ord:= function( mu, nu )
208
209        # true if mu < nu
210
211        local cf, mu0, nu0;
212
213        mu0:= mu{[1..semsimrk]};
214        nu0:= nu{[1..semsimrk]};
215        if mu0 = nu0 then return fail; fi;
216        cf:= Coefficients( b, nu0-mu0 );
217        if ForAll( cf, IsInt ) then
218           if ForAll( cf, x -> x >= 0 ) then
219              return true;
220           elif ForAll( cf, x -> x <= 0 ) then
221              return false;
222           fi;
223        fi;
224
225        return fail;
226
227    end;
228
229    ww:=[ ]; mm:= [ ];
230
231    while Length(w0) > 0 do
232
233        # look for a maximum in w0...
234
235        w:= w0[1];
236        mult:= m0[1];
237        for i in [1..Length(w0)] do
238            isg:= ord( w, w0[i] );
239            if isg <> fail and isg then
240                 w:= w0[i]; mult:= m0[i];
241            fi;
242        od;
243
244        Add( ww, w ); Add( mm, mult );
245        ch:= DominantCharacter( K, w{[1..semsimrk]} );
246        for i in [1..Length(ch[1])] do
247	    mu:= Concatenation( ch[1][i], w{[semsimrk+1..Length(w)]} );
248            pos:= Position( w0, mu );
249            m0[pos]:= m0[pos] - mult*ch[2][i];
250        od;
251        for i in [1..Length(w0)] do
252            if m0[i] = 0 then
253               Unbind( m0[i] );
254               Unbind( w0[i] );
255            fi;
256        od;
257        m0:= Filtered( m0, x -> IsBound(x) );
258        w0:= Filtered( w0, x -> IsBound(x) );
259    od;
260
261    return [ww,mm];
262
263end );
264
265
266
267
268InstallMethod( DirectSumDecomposition,
269"for a module over a semisimple Lie algebra",
270true, [ IsAlgebraModule ], 0,
271
272# IN FACT: we assume that it is a left algebra module!
273
274function( V )
275
276   local L, R, e, sp, x, m, c, sp0, K, h, U, es;
277
278   L:= LeftActingAlgebra( V );
279
280   # We assume that L is reductive, and that the action is "algebraic".
281
282   if LieDerivedSubalgebra( L ) <> L then
283      K:= LieDerivedSubalgebra(L);
284   else
285      K:= L;
286   fi;
287
288   R:= RootSystem(K);
289   if R = fail then return fail; fi;
290   e:= CanonicalGenerators( R )[1];
291
292   sp:= V;
293   for x in e do
294       m:= MatrixOfAction( Basis(V), x );
295       m:= TransposedMatDestructive(m);
296       c:= NullspaceMat(m);
297       sp0:= Subspace( V, List( c, u -> u*Basis(V) ), "basis" );
298       sp:= Intersection( sp, sp0 );
299   od;
300
301   # in sp find a basis of weight vectors
302   h:= CanonicalGenerators( R )[3];
303   sp:= [ sp ];
304   for x in h do
305       sp0:= [ ];
306       for U in sp do
307           m:= List( Basis(U), a -> Coefficients( Basis(U), x^a ) );
308           es:= Eigenspaces(LeftActingDomain(L),m);
309           Append( sp0, List( es, M -> Subspace( V, List( Basis(M), v ->
310                                          v*Basis(U) ) ) ) );
311       od;
312       sp:= sp0;
313   od;
314
315   sp0:= [ ];
316   for U in sp do
317       Append( sp0, List( Basis(U), u -> SubAlgebraModule( V, [u] ) ) );
318   od;
319   return sp0;
320
321
322end );
323
324
325
326InstallMethod( CharacteristicsOfStrata,
327"for a semisimple Lie algebra and a highest weight",
328true, [ IsLieAlgebra, IsList ], 0,
329
330function( L, hw ) # L: semisimple Lie algebra, hw: highest weight
331
332    local ch, W, wts, w, o, hh, KM, hwts, perms, rk, i, R, B, rho, exprs,
333          p, a, sa, mults, H, BH, ip, posR, csa, inconvexhull, hdimstrata,
334          weylorbs, j, A, sim, id;
335
336# first we have three functions...
337
338weylorbs:= function( W, wts, A )
339
340     # wts is a set of weights, in usual notation, permuted
341     # by the Weyl group, we find subsets of length up to the rank
342     # containing reps of each orbit of the Weyl group, but
343     # probably rather redundant.
344
345     local sets, zeros, rk, s1, z1, i, j, k, s, z, s0, z0, inds, len, cont,
346           r, l, OO, dt, dets, mat, rows, row, cols, col;
347
348     sets:= [ ];
349     zeros:= [ ];
350     rk:= Length( wts[1] );
351
352     s1:= [ ];
353     z1:= [ ];
354     for i in [1..Length(wts)] do
355         if ForAll( wts[i], m -> m >= 0 ) then
356            Add( s1, [i] );
357            Add( z1, Filtered( [1..rk], m -> wts[i][m] = 0 ) );
358         fi;
359     od;
360     Add( sets, s1 );
361     Add( zeros, z1 );
362
363     for len in [2..rk] do
364         s1:= [ ]; z1:= [ ]; dets:= [ ]; rows:= [ ]; #cols:= [ ];
365         for i in [1..Length(sets[len-1])] do
366             s:= sets[len-1][i];
367             z:= zeros[len-1][i];
368             for j in [1..Length(wts)] do
369                 if not j in s then
370                    if ForAll( z, m -> wts[j][m] >= 0 ) then
371                       s0:= ShallowCopy(s); Add( s0, j );
372                       z0:= [ ];
373                       for k in [1..rk] do
374                           if k in z and wts[j][k]=0 then
375                              Add( z0, k );
376                           fi;
377                       od;
378                       cont:= false;
379                       mat:= List( [1..Length(s0)], s -> [ ] );
380                       for k in [1..Length(s0)] do
381                           for l in [k..Length(s0)] do
382                               mat[k][l]:= wts[s0[k]]*A*wts[s0[l]];
383                               mat[l][k]:= mat[k][l];
384                           od;
385                       od;
386                       dt:= Permanent(mat);
387                       row:= List( mat, Sum );
388                       Sort( row );
389
390                       #col:= List( TransposedMat(mat), Sum );
391                       #Sort(col);
392                       for l in [1..Length(s1)] do
393                           if dt = dets[l] and row = rows[l] #and col = cols[l]
394                                                             then
395                              r:= RepresentativeAction( W, s0, s1[l], OnSets );
396                              if r <> fail then
397                                 cont:= true;
398                                 break;
399                              fi;
400                           fi;
401                       od;
402
403                       if not cont then
404                          Add( s1, s0 ); Add( z1, z0 );
405                          Add( dets, dt ); Add( rows, row ); #Add( cols, col );
406                       fi;
407                    fi;
408                 fi;
409             od;
410         od;
411         Add( sets, s1 ); Add( zeros, z1 );
412     od;
413
414     return sets;
415
416end;
417
418
419inconvexhull:= function( B, S0, p0, dist, ip )
420                                            # S set of vecs in R^m (rat coords),
421                                            # p a point in R^m, is p\in S?
422                                            # dist: distance fct
423
424    local m, i, one, eps, dists, pos, v, pp, k, j, u, t, S, p;
425
426    S:= List( S0, x -> Coefficients( B, x ) );
427    p:= Coefficients( B, p0 );
428    one:= 1.000000000000000000000;
429    S:= List( S, u -> u*one );
430    p:= p*one;
431
432    eps:= one*1/10; # small eps, but not very small, the algorithm for
433                    # the strata is rather sensitive to the choice of eps,
434                    # smaller eps, longer runtime...
435
436    dists:= List( S, s -> dist(s,p) );
437    pos:= Position( dists, Minimum( dists ) );
438    v:= S[pos];
439    pp:= S[pos];
440
441    while true do
442       if dist(p,pp) < eps*dist(p,v) then
443          return [ pp, true ];
444       else
445          k:= 0;
446          for j in [1..Length(S)] do
447              if dist(pp,S[j]) >= dist(p,S[j]) then
448                 k:= j; break;
449              fi;
450          od;
451          if k > 0 then
452             v:= S[k];
453          else
454             return [ pp, false ];
455          fi;
456       fi;
457
458       u:= pp-v;
459       t:= -ip(u,p-pp)/ip(u,u);
460       pp:= pp+t*(-u);
461
462    od;
463end;
464
465hdimstrata:= function( R, H, BH, ip, rk, posR, csa, KM, h_wts, mults, wts, W, exps, A, bool )
466
467    # R: the root system of the big Lie alg
468    # H: Cartan subalgebra of the "big" Lie alg
469    # BH: basis of H, wrt Chevalley basis vecs
470    # ip: inner product function on H (fct of two arguments)
471    # rk: rank of the input Lie algebra (in general a reductive subalgebra)
472    # posR: pos roots of input Lie alg, as weights,
473    # csa : basis of the Cartan of the input Lie alg
474    # h_wts: List of wts of the rep, as elts in the CSA,
475    # mults: their multiplicities
476    # wts: the weights in usual form
477    # W : list of perms on [1..Length(wts)], giving the action on the
478    # Weyl group
479    # exps: reduces expressions of the reflections corr to the rts in posR
480    # bool: a boolean, if true then we do not consider only the stata of
481    #       highest dim.
482
483    local dist, out, i, G, k, XX, OO, orbs, st, hst, diffs, V, j, hset, h,
484          mat, sol, bas, B, dims, dim, cfs, p, inds, pR, wt0, ex0,
485          h0, cs, hw0, mu0, perms, res, hs, dist0, w, ip0, bcsa, BC, bigger,
486          sums, delt, xx, yy, hh, conjdomh, Onew, d, r, u, v, rep, len, N,
487          totlen, dets, dt, orbs0, q;
488
489    dist:= function(u,v) return ip(u-v,u-v); end;
490
491    dist0:= function(u,v) return (u-v)*KM*(u-v); end;
492    ip0:= function(u,v) return u*KM*v; end;
493
494    # detect the trivial rep, ie all weights are 0
495    w:= List( h_wts, v -> List( csa, u -> ip(u,v) ) );
496    if ForAll( w, IsZero ) then
497       # no nilpotent elements
498       return [ [], [] ];
499    fi;
500
501    if Length(posR)=0 then  #torus
502       out:= [ [], [] ];
503       for i in [1..Length(h_wts)] do
504           w:= List( csa, u -> ip(u,h_wts[i]) ); # so the weight corr.
505           if not IsZero(w) then
506              Add( out[1], h_wts[i] );
507              Add( out[2], mults[i] );
508           fi;
509       od;
510       return out;
511    fi;
512
513    # compute simple roots, Weyl group etc.
514    # to later erase W-conjugate h-s...
515
516    sums:= [ ];
517    for i in [1..Length(posR)] do
518        for j in [i+1..Length(posR)] do
519            Add( sums, posR[i]+posR[j] );
520        od;
521    od;
522    delt:= Filtered( posR, x -> not x in sums );
523
524    inds:= List( delt, x -> Position( PositiveRootsAsWeights(R), x ) );
525    xx:= PositiveRootVectors(R){inds};
526    yy:= NegativeRootVectors(R){inds};
527    hh:= List( [1..Length(inds)], i -> xx[i]*yy[i] );
528
529    conjdomh:= function( h )
530           # the conjugate dominant elt in H to h.
531
532       local h0, a, c, pos;
533
534       h0:= h;
535       while true do
536          c:= Coefficients( BH, h0 );
537          a:= List( delt, u -> u*c );
538          pos:= PositionProperty( a, x -> x<0 );
539          if pos = fail then
540             return h0;
541          else
542             h0:= h0-a[pos]*hh[pos];
543          fi;
544       od;
545
546    end;
547
548
549       G:= Group(W);
550       hs:= [ ];
551
552       if bool then orbs0:= weylorbs( G, wts, A ); fi;
553
554       for k in [1..rk] do
555           if bool then
556              OO:= orbs0[k];
557           else
558              N:= NrCombinations( [1..Length(h_wts)], k );
559
560              if N <= 500000 then # some bound for memory reasons...
561                 XX:= Combinations( [1..Length(h_wts)], k );
562                 OO:= OrbitsDomain( G, XX, OnSets );
563                 OO:= List( OO, u -> u[1] );
564
565              else
566                 Onew:= [ ];
567                 totlen:= 0;
568                 dets:= [ ];
569                 for u in OO do
570                     d:= u[Length(u)];
571                     v:= ShallowCopy(u);
572                     len:= Length(v)+1;
573                     for i in [d+1..Length(h_wts)] do
574
575                         v[len]:= i;
576                         dt:= Permanent( List( v, i -> List( v,
577                                 j -> ip( h_wts[i], h_wts[j] ) ) ) );
578
579                         rep:= false;
580                         for j in [1..Length(Onew)] do
581                             if dets[j] = dt then
582                                r:= RepresentativeAction( G, v, Onew[j], OnSets );
583                                if r <> fail then
584                                   rep:= true;
585                                   break;
586                                fi;
587                             fi;
588                         od;
589                         if not rep then
590                            Add( Onew, ShallowCopy(v) );
591                            Add( dets, dt );
592                            totlen:= totlen + OrbitLength(G,v,OnSets);
593                         fi;
594                         if totlen=N then break; fi;
595                     od;
596                 od;
597                 OO:= Onew;
598
599              fi;
600           fi;
601
602           orbs:= Filtered( OO, v -> not 0*h_wts[1] in
603                       h_wts{v} );
604           for st in orbs do
605               # find the affine space generated by st...
606
607               hst:= h_wts{st};
608               diffs:= [ ];
609               for i in [1..Length(hst)] do
610                   for j in [i+1..Length(hst)] do
611                       Add( diffs, hst[j]-hst[i] );
612                   od;
613               od;
614               V:= Subspace( H, diffs );
615
616               # so the affine space is A=hst[1]+V...
617               # see if there are more elms in this space
618
619               hset:= ShallowCopy( hst );
620               #bigger:= false;
621               for h in h_wts do
622                   if h-hst[1] in V and not h in hset then
623                      Add( hset, h );
624                      #bigger:= true;
625                   fi;
626               od;
627
628               #if not bigger or k=rk then
629
630                 # find point closest to 0 on A:
631                 if Length(st)>1 then
632                    mat:= List( [1..Dimension(H)], i ->
633                            List( [1..Dimension(V)], j ->
634                                     ip( Basis(H)[i], Basis(V)[j] ) ) );
635                    sol:= List( NullspaceMat( mat ), u -> u*Basis(H) );
636
637                   # so this is a basis of the space orthogonal to V...
638                   # so sol and V span the whole space, therefore, finding the
639                 # intersection point amounts to expressing hst[1] on this basis
640
641                    bas:= -ShallowCopy( Basis(V) );
642                    Append( bas, sol );
643                    B:= Basis( H, bas );
644                    cfs:= Coefficients( B, hst[1] );
645                    h0:= hst[1]+cfs{[1..Dimension(V)]}*Basis(V);
646                 else
647                    h0:= h_wts[ st[1] ];
648                 fi;
649
650                 # see whether it is contained
651                 # in the convex hull of hst.
652
653                 if not IsZero(h0) then
654                    if Length(hst)=1 or
655                           inconvexhull( BH, hst, h0, dist0, ip0 )[2] then
656                       h0:= (2/ip(h0,h0))*h0;
657                       h0:= conjdomh(h0);
658                       if not h0 in hs then
659                          Add( hs, h0 );
660                       fi;
661                    fi;
662                 #fi;
663              fi;
664           od;
665       od;
666
667       # now for each elt of hs we compute the dimension of the corr stratum
668       # (in case it is a characteristic).
669
670       dims:= [ ];
671       for h0 in hs do
672          dim:= 2*Length(posR)+Length(csa);
673          for i in [1..Length(h_wts)] do
674              if ip(h0,h_wts[i]) >= 2 then dim:= dim+mults[i]; fi;
675          od;
676          cfs:= Coefficients( BH, h0 );
677          for p in posR do
678              if cfs*p <> 0 then # so either the pos or the neg root contrbs
679                 dim:= dim-1;
680              else
681                 # now both the pos and the neg root contribs
682                 dim:= dim-2;
683              fi;
684          od;
685          dim:= dim-Length(csa);
686          Add( dims, dim );
687       od;
688
689       if not bool then
690          dim:= 0; # compute maximum dim, less than total dim of V...
691          for i in [1..Length(dims)] do
692              if dims[i] > dim and dims[i] <= Sum( mults ) then
693                 dim:= dims[i];
694              fi;
695          od;
696          inds:= Filtered( [1..Length(hs)], i -> dims[i] = dim );
697          hs:= hs{inds};
698       else
699          dim:= Sum( mults );
700          inds:= Filtered( [1..Length(hs)], i -> dims[i] <= dim );
701          hs:= hs{inds};
702       fi;
703
704       # now for each elt of hs, do a recursive call...
705
706       out:= [ [], [] ]; # first the chars, then the dims of the corr stratum.
707
708       for k in [1..Length(hs)] do
709
710           h0:= hs[k];
711           # have to compute tilde z(h)...
712           # first its posR...
713
714           pR:= [ ];
715           ex0:= [ ];
716           wt0:= [ ];
717           cfs:= Coefficients( BH, h0 );
718           for i in [1..Length(posR)] do
719               if cfs*posR[i] = 0 then
720                  Add( pR, posR[i] );
721                  Add( ex0, exps[i] );
722               fi;
723           od;
724
725
726           # next the csa:
727           mat:= List( csa, h1 -> [ip(h1,h0)] );
728           cs:= List( NullspaceMat(mat), u -> u*csa );
729
730           bcsa:= ShallowCopy(cs); Add( bcsa, h0 );
731           BC:= Basis( Subspace(H,bcsa), bcsa );
732           # now the h_wts of V_2(h0):
733           # those are the elts h of h_wts such that ip(h0,h)=2...
734           hw0:= [ ];
735           mu0:= [ ];
736           for i in [1..Length(h_wts)] do
737               if ip(h0,h_wts[i]) = 2 then
738                  # project h_wts[i] on cs:
739                  if Length(cs) > 0 then
740                     cfs:= Coefficients( BC, h_wts[i] );
741                     Add( hw0, cfs{[1..Length(cs)]}*cs );
742                  else
743                     Add( hw0, Zero(H) );
744                  fi;
745                  Add( mu0, mults[i] );
746                  Add( wt0, wts[i] );
747               fi;
748           od;
749
750           # and finally the Weyl group:
751           perms:= [ ];
752           for i in [1..Length(pR)] do
753               Add( perms, PermList(
754                 List( wt0, w -> Position( wt0,
755                               ApplyWeylElement(WeylGroup(R),w,ex0[i]) ))));
756           od;
757           res:= hdimstrata( R, H, BH, ip, rk-1, pR, cs, KM, hw0, mu0, wt0,
758                             perms, ex0, A, false );
759
760           # is there a stratum in the recursive output of dimension equal
761           # to the dimension of V_2(h)? If yes, then h is not a char.
762
763           if not Sum( mu0 ) in res[2] then
764              # so h0 is a char
765              if bool then
766                 Add( out[1], h0 ); Add( out[2], dims[k] );
767              else
768                 Add( out[1], h0 ); Add( out[2], dim );
769              fi;
770           fi;
771
772       od;
773       return out;
774
775end;
776
777   # Now we start with the main function.
778
779    # first we compute the character of the module,
780    # then map the weights to the CSA, where we use the Killing form.
781
782    R:= RootSystem(L);
783    ch:= DominantCharacter( L, hw );
784    W:= WeylGroup( R );
785    wts:= [ ];
786    mults:= [ ];
787    for w in [1..Length(ch[1])] do
788        o:= WeylOrbitIterator( W, ch[1][w] );
789        while not IsDoneIterator( o ) do
790            Add( wts, NextIterator(o) );
791            Add( mults, ch[2][w] );
792        od;
793    od;
794
795    hh:= ChevalleyBasis(L)[3];
796    rk:= Length(hh);
797    KM:= List( hh, h1 -> List( hh, h2-> TraceMat( AdjointMatrix(Basis(L),h1)*
798           AdjointMatrix(Basis(L),h2) ) ) );
799
800    hwts:= [ ];
801    for w in wts do
802        Add( hwts, SolutionMat( KM, w )*hh );
803    od;
804
805    # next we get permutations generating the permutation action of W on the
806    # weights:
807
808    perms:= [ ];
809    for i in [1..rk] do
810        Add( perms, PermList(
811              List( wts, w -> Position( wts, ApplyWeylElement(W,w,[i]) ))));
812    od;
813
814    # find reduced expressions for the reflections corr the pos rts:
815    B:= BilinearFormMatNF(R);
816    rho:= List( [1..rk], i -> 1);
817    exprs:= [ ];
818    for i in [1..Length(PositiveRoots(R))] do
819        p:= PositiveRootsNF(R)[i];
820        a:= Sum( [1..rk], j -> p[j]*B[j][j] );
821        sa:= rho -a/(p*B*p)*PositiveRootsAsWeights(R)[i];
822        w:= ConjugateDominantWeightWithWord( W, sa );
823        Add( exprs, w[2] );
824    od;
825
826    H:= CartanSubalgebra(L);
827    BH:= Basis(H, hh );
828
829    ip:= function(h1,h2) return Coefficients(BH,h1)*KM*Coefficients(BH,h2);
830         end;
831    posR:= PositiveRootsAsWeights(R);
832    csa:= ShallowCopy( hh );
833
834    sim:= SimpleRootsAsWeights( R );
835    sim:= Basis( VectorSpace( Rationals, sim ), sim );
836    id:= IdentityMat(rk);
837    A:= List( [1..rk], i -> [ ] );
838    for i in [1..rk] do
839        for j in [i..rk] do
840            A[i][j]:= Coefficients(sim,id[i])*B*Coefficients(sim,id[j]);
841            A[j][i]:= A[i][j];
842        od;
843
844    od;
845
846    return  hdimstrata( R, H, BH, ip, rk, posR, csa, KM, hwts, mults,
847                 wts, perms, exprs, A, true );
848
849end );
850
851
852
853