1InstallMethod( FiniteOrderInnerAutomorphisms,
2"for string, integer and integer", true, [ IsString, IsInt, IsInt ], 0,
3function( type, rank, m )
4
5    # finite order auts of the simple Lie algebra,
6    # of order m, that correspond to untwisted diagrams,
7    # in other words, they correpond to the identity
8    # diagram automorphism.
9
10    local L, w, cc, ch, g, a, ss, good, s, i, auts, g0, t, G, f, stack,
11          stack0, j, u, p1, p2, list, n;
12
13    n:= rank;
14    if type = "A" then
15       list:= [2..n+1]; Add( list, 1 );
16       p1:= PermList( list );
17       p2:= ();
18       i:= 1;
19       while 1+i < n+2-i do
20           p2:= p2*(1+i,n+2-i);
21           i:= i+1;
22       od;
23       G:= Group( [ p1, p2 ] );
24    elif type = "B" then
25       if rank = 2 then
26          G:= Group([(1,3)]);
27       else
28          G:= Group([(1,2)]);
29       fi;
30    elif type = "C" then
31       p2:= ();
32       i:= 1;
33       while i < n+2-i do
34           p2:= p2*(i,n+2-i);
35           i:= i+1;
36       od;
37       G:= Group( [p2] );
38    elif type = "D" and rank > 4 then
39       p2:= ();
40       i:= 1;
41       while i < n+2-i do
42           p2:= p2*(i,n+2-i);
43           i:= i+1;
44       od;
45       G:= Group( [(1,2),(n,n+1),p2] );
46    elif type = "D" and rank = 4 then
47       G:= Group( [ (1,2,4,5), (1,2) ] );
48    elif type = "E" and rank = 6 then
49       G:= Group( [ (1,2,7)*(3,4,6), (1,2)*(3,4)] );
50    elif type = "E" and rank = 7 then
51       G:= Group( [ (1,8)*(2,7)*(4,6) ] );
52    else
53       G:= Group( [ () ] );
54    fi;
55
56    G:= Elements( G );
57    G:= Filtered( G, x -> x <> x^0 );
58
59    L:= SimpleLieAlgebra( type, rank, CF(m) );
60    w:= E(m);
61
62    cc:= ExtendedCartanMatrix( RootSystem(L) );
63
64    ch:= ChevalleyBasis(L);
65    g:= [ ch[2][ Length(ch[2]) ] ];
66    Append( g, ch[1]{[1..rank]} );
67
68    a:= cc.labels;
69    ss:= [ ];
70
71    stack:= [ List( g, x -> 0 ) ];
72
73    for i in [1..Length(g)] do
74
75        stack0:= [ ];
76        for s in stack do
77            u:= a*s;
78            if u = m and Gcd(s) = 1 then
79               good:= true;
80               for p1 in G do
81                   t:= Permuted( s, p1 );
82                   if t in ss then
83                      good:= false; break;
84                   fi;
85               od;
86               if good then
87                  Add( ss, s );
88               fi;
89            elif u < m then
90               for j in [0..m-u] do
91                   t:= ShallowCopy(s);
92                   t[i]:= j;
93                   Add( stack0, t );
94               od;
95            fi;
96        od;
97        stack:= stack0;
98    od;
99
100    for s in stack do
101        u:= a*s;
102        if u = m and Gcd(s) = 1 then
103           good:= true;
104           for p1 in G do
105               t:= Permuted( s, p1 );
106               if t in ss then
107                  good:= false; break;
108               fi;
109           od;
110           if good then
111              Add( ss, s );
112           fi;
113        fi;
114    od;
115
116
117    auts:= [ ];
118
119    for s in ss do
120        g0:= List( [1..Length(g)], i -> w^s[i]*g[i] );
121        f:= AlgebraHomomorphismByImagesNC( L, L, g, g0 );
122        SetOrder(f,m);
123        SetKacDiagram( f, rec( CM:= cc.ECM, labels:= cc.labels, weights:= s ) );
124        Add( auts, f );
125    od;
126
127    return auts;
128
129
130end );
131
132
133
134InstallMethod( FiniteOrderOuterAutomorphisms,
135"for string, and three integers", true, [ IsString, IsInt, IsInt, IsInt ], 0,
136function( type, rank, m, d )
137
138
139     # corresponding to the diagram automorphism of
140     # order d.
141
142     local phi, L, w, R, cg, cg0, sim, i, pos, f, mat, mat0, sol, cK, H, K,
143           V, y, rt, sp, h, g, rts, B, C, j, v, a, ss, done, s, auts, g0, G,
144           t, n, p2, good, u, p1, stack, stack0;
145
146     phi:= function( rt )
147
148        local r0;
149        if type = "A"  then
150           return Reversed(rt);
151        elif type = "D" and d=2 then
152           r0:= ShallowCopy( rt );
153           r0[ rank ]:= rt[rank-1];
154           r0[rank-1]:= rt[rank];
155        elif type = "D" and d = 3 then
156           if rank <> 4 then
157              Error( "only D_4 has a diagram automorphism of order 3");
158           fi;
159           r0:= ShallowCopy( rt );
160           r0[1]:= rt[4]; r0[3]:= rt[1]; r0[4]:= rt[3];
161        else
162           r0:= ShallowCopy(rt);
163           r0[1]:= rt[6]; r0[6]:= rt[1];
164           r0[3]:=rt[5]; r0[5]:= rt[3];
165        fi;
166        return r0;
167
168     end;
169
170    if type = "A" and rank = 1 then return [ ]; fi;
171
172    if type ="D" and d = 2 then
173       n:= rank-1;
174       p2:= ();
175       i:= 1;
176       while i < n+2-i do
177           p2:= p2*(i,n+2-i);
178           i:= i+1;
179       od;
180       G:= Group( [ p2 ] );
181       G:= Elements( G );
182       G:= Filtered( G, x -> x <> x^0 );
183    elif type = "A" and IsOddInt(rank) then
184       G:= [ (1,2) ];
185    else
186       G:= [ ];
187    fi;
188
189
190    if d=2 then
191       L:= SimpleLieAlgebra( type, rank, CF(m) );
192    else
193       L:= SimpleLieAlgebra( type, rank, CF(3*m) );
194    fi;
195
196    w:= E(m);
197
198    R:= RootSystem(L);
199    cg:= CanonicalGenerators( R );
200    cg0:= [ [], [], [] ];
201
202    sim:= SimpleSystemNF( R );
203    for i in [1..Length(sim)] do
204        pos:= Position( sim, phi(sim[i]) );
205        Add( cg0[1], cg[1][pos] );
206        Add( cg0[2], cg[2][pos] );
207        Add( cg0[3], cg[3][pos] );
208    od;
209
210    f:= AlgebraHomomorphismByImagesNC( L, L, Flat(cg), Flat(cg0) );
211
212    mat:= [ ];
213    for i in [1..Dimension(L)] do
214        Add( mat, Coefficients( Basis(L), Image( f, Basis(L)[i] ) ) );
215    od;
216
217    mat0:= mat- IdentityMat( Dimension(L) );
218
219    sol:= NullspaceMat( mat0 );
220
221    K:= Subalgebra( L, List( sol, x -> LinearCombination(Basis(L),x) ) );
222
223    cK:= CanonicalGenerators( RootSystem(K) );
224
225    if d=2 then
226       mat0:= mat+IdentityMat( Dimension(L) );
227    else
228       mat0:= mat-E(3)*IdentityMat( Dimension(L) );
229    fi;
230
231    sol:= NullspaceMat( mat0 );
232    V:= LeftAlgebraModuleByGenerators( K, function(x,v) return x*v; end,
233           List( sol, x -> LinearCombination( Basis(L), x ) ) );
234
235    y:= List( cK[2], x -> MatrixOfAction( Basis(V), x ) );
236
237    # get simultaneous kernel...
238
239    mat:= y[1];
240    for i in [2..Length(y)] do Append( mat, y[i] ); od;
241
242    sol:= NullspaceMat( TransposedMatDestructive(mat) );
243
244    g:= [ LinearCombination( Basis(V), sol[1] )![1] ];
245
246    sp:= Basis( VectorSpace( LeftActingDomain(L), g ), g );
247
248    rt:= [ ];
249    for h in cK[3] do
250        Add( rt, Coefficients( sp, h*g[1] )[1] );
251    od;
252
253    sim:= SimpleRootsAsWeights( RootSystem(K) );
254    sp:= Basis( VectorSpace( Rationals, sim ), sim );
255    rts:= [ Coefficients( sp, rt ) ];
256    Append( rts, SimpleSystemNF(RootSystem(K) ) );
257
258    B:= BilinearFormMatNF( RootSystem(K) );
259    C:= NullMat( Length(rts), Length(rts) );
260    for i in [1..Length(rts)] do
261        for j in [1..Length(rts)] do
262            C[i][j]:= 2*( rts[i]*(B*rts[j]) )/( rts[j]*(B*rts[j]) );
263        od;
264    od;
265
266    v:= NullspaceMat(C)[1];
267    a:= Lcm( List( v, DenominatorRat ) );
268    v:= a*v;
269
270    Append( g, cK[1] );
271
272    ss:= [ ];
273
274    stack:= [ List( g, x -> 0 ) ];
275
276    for i in [1..Length(g)] do
277
278        stack0:= [ ];
279        for s in stack do
280            u:= d*(v*s);
281            if u = m and Gcd(s) = 1 then
282               good:= true;
283               for p1 in G do
284                   t:= Permuted( s, p1 );
285                   if t in ss then
286                      good:= false; break;
287                   fi;
288               od;
289               if good then
290                  Add( ss, s );
291               fi;
292            elif u < m then
293               for j in [0..m-u] do
294                   t:= ShallowCopy(s);
295                   t[i]:= j;
296                   Add( stack0, t );
297               od;
298            fi;
299        od;
300        stack:= stack0;
301    od;
302    for s in stack do
303        u:= d*(v*s);
304        if u = m and Gcd(s) = 1 then
305           good:= true;
306           for p1 in G do
307               t:= Permuted( s, p1 );
308               if t in ss then
309                  good:= false; break;
310               fi;
311           od;
312           if good then
313              Add( ss, s );
314           fi;
315        fi;
316    od;
317
318    auts:= [ ];
319    for s in ss do
320        g0:= List( [1..Length(g)], i -> w^s[i]*g[i] );
321        f:= AlgebraHomomorphismByImagesNC( L, L, g, g0 );
322        SetOrder(f,m);
323        SetKacDiagram(f,rec( CM:= C, labels:= v, weights:= s ));
324        Add( auts, f );
325    od;
326
327    return auts;
328
329end );
330
331
332InstallOtherMethod( Grading,
333"for a finite order automorphism", true, [ IsGeneralMapping ], 0,
334function( f )
335
336    local L, m, w, mat, id, spaces, i, sp;
337
338    L:= Source(f);
339    m:= Order(f);
340    w:= E(m);
341
342    mat:= List( Basis(L), x -> Coefficients( Basis(L), Image(f,x) ) );
343    id:= mat^0;
344
345    spaces:= [ ];
346    for i in [0..m-1] do
347        sp:= NullspaceMat( mat - w^i*id );
348        Add( spaces, List( sp, x -> LinearCombination(Basis(L),x) ) );
349    od;
350
351    return spaces;
352
353end );
354
355
356
357SLAfcts.nil_orbs_inner:= function( L, gr0, gr1, gr2 )
358
359     # Here L is a simple graded Lie algebra; gr0 a basis of the
360     # elts of degree 0, gr1 of degree 1, and gr2 of degree -1.
361     # We find the nilpotent G_0-orbits in g_1.
362     # We assume that the given CSA of L is also a CSA of g_0.
363
364     local F, g0, s, r, HL, Hs, R, Ci, hL, hl, C, rank, posRv_L, posR_L,
365           posR, i, j, sums, fundR, inds, tr, h_candidates, BH, W, h,
366           c_h, ph, stb, v, w, is_rep, h0, wr, Omega, good_h, g1, g2, h_mats1,
367           h_mats2, mat, sl2s, id1, id2, V, e, f, bb, ef, found, good, co, x,
368           C_h0, sp, sp0, y, b, bas, c, Cs, B, k, sol, info;
369
370     F:= LeftActingDomain(L);
371
372     g0:= Subalgebra( L, gr0, "basis" );
373
374     s:= LieDerivedSubalgebra( g0 );
375     r:= LieCentre(g0);
376
377     HL:= CartanSubalgebra(L);
378     Hs:= Intersection( s, HL );
379     SetCartanSubalgebra( s, Hs );
380
381     R:= RootSystem(L);
382     Ci:= CartanMatrix( R )^-1;
383     hL:= ChevalleyBasis(L)[3];
384     hl:= List( NilpotentOrbits(L), x -> (Ci*WeightedDynkinDiagram(x))*hL );
385
386     for i in [1..Length(hl)] do
387         if hl[i] = 0*hl[i] then
388            Unbind( hl[i] );
389         fi;
390     od;
391     hl:= Filtered( hl, x -> IsBound(x) );
392
393     C:= CartanMatrix( R );
394     rank:= Length(C);
395
396     # we have to compute a root system of s such that its
397     # positive roots are also positive in L...
398     # Note that since the CSA of s is a subset of the CSA of L,
399     # the roots of s are a subset of the roots of L; also:
400     # the part of the CSA of L that is not in s, commutes with s,
401     # the coordinates of the roots of s with respect to those h-s
402     # is zero (if you understand what I mean...).
403
404     posRv_L:= PositiveRootVectors(R);
405     posR_L:= PositiveRootsNF(R);
406     posR:= [ ];
407     for i in [1..Length(posRv_L)] do
408         if posRv_L[i] in s then
409            Add( posR, posR_L[i] );
410         fi;
411     od;
412
413     sums:= [ ];
414     for i in [1..Length(posR)] do
415         for j in [i+1..Length(posR)] do
416             Add( sums, posR[i]+posR[j] );
417         od;
418     od;
419     fundR:= Filtered( posR, x -> not x in sums );
420     inds:= List( fundR, x -> Position( posR_L, x ) );
421     tr:= WeylTransversal( R, inds );
422
423info:= "Constructed a Weyl transversal of ";
424Append( info, String(Length(tr)) );
425Append( info, " elements.");
426
427Info(InfoSLA,2,info);
428
429     h_candidates:= [ ];
430     BH:= Basis( VectorSpace( F, hL ), hL );
431     W:= WeylGroup(R);
432     for h in hl do
433
434         # first we get the indices of the simple reflections that
435         # stabilise h...
436         c_h:= Coefficients( BH, h );
437         ph:= C*c_h;
438         stb:= Filtered( [1..rank], k -> ph[k] = 0 );
439
440         for w in tr do
441             is_rep:= true;
442             for i in stb do
443                 # see whether there is an expression for w ending with i
444                 v:= ShallowCopy(w); Add( v, i );
445                 if LengthOfWeylWord( W, v ) < Length(w) then
446                    is_rep:= false;
447                    break;
448                 fi;
449             od;
450             if is_rep then
451                h0:= ShallowCopy(c_h);
452                wr:= Reversed(w);
453                for i in wr do
454                    h0[i]:= h0[i] - (C[i]*h0);
455                od;
456                AddSet( h_candidates, h0 );
457             fi;
458         od;
459     od;
460
461
462info:=  "Constructed ";
463Append( info, String( Length(h_candidates) ) );
464Append( info, " Cartan elements to be checked." );
465Info( InfoSLA, 2, info );
466
467     # now we need to compute sl_2 triples wrt the h-s found...
468
469     Omega:= [0..Dimension(L)];
470     good_h:= [ ];
471
472     g1:= Basis( Subspace( L, gr1 ), gr1 );
473     g2:= Basis( Subspace( L, gr2 ), gr2 );
474
475     # the matrices of hL[i] acting on g1
476     h_mats1:= [ ];
477     for h0 in hL do
478         mat:= [ ];
479         for i in [1..Length(g1)] do
480             Add( mat, Coefficients( g1, h0*g1[i] ) );
481         od;
482         Add( h_mats1, mat );
483     od;
484
485     # those of wrt g2...
486     h_mats2:= [ ];
487     for h0 in hL do
488         mat:= [ ];
489         for i in [1..Length(g1)] do
490             Add( mat, Coefficients( g2, h0*g2[i] ) );
491         od;
492         Add( h_mats2, mat );
493     od;
494
495     sl2s:= [ ];
496     id1:= IdentityMat( Length(g1) );
497     id2:= IdentityMat( Length(g2) );
498     for h in h_candidates do
499
500         mat:= h*h_mats1;
501         mat:= mat - 2*id1;
502         V:= NullspaceMat( mat );
503         e:= List( V, v -> v*gr1 );
504
505         mat:= h*h_mats2;
506         mat:= mat + 2*id2;
507         V:= NullspaceMat( mat );
508         f:= List( V, v -> v*gr2 );
509
510         # check whether h0 in [e,f]....
511         bb:= [ ];
512         for x in e do
513             for y in f do
514                 Add( bb, x*y );
515             od;
516         od;
517         ef:= Subspace( L, bb );
518
519         h0:= h*hL;
520
521         if h0 in ef then  #otherwise we can just discard h...
522            found:= false;
523            good:= false;
524            while not found do
525
526                co:= List( e, x -> Random(Omega) );
527                x:= co*e;
528                sp:= Subspace( L, List( f, y -> x*y) );
529
530                if Dimension(sp) = Length(e) and h0 in sp then
531
532                   # look for a nice one...
533                   for i in [1..Length(co)] do
534                       k:= 0;
535                       found:= false;
536                       while not found do
537                           co[i]:= k;
538                           x:= co*e;
539                           sp:= Subspace( L, List( f, y -> x*y) );
540
541                           if Dimension(sp) = Length(e) and h0 in sp then
542                              found:= true;
543                           else
544                              k:= k+1;
545                           fi;
546                       od;
547                   od;
548
549                   mat:= List( f, u -> Coefficients( Basis(sp), x*u ) );
550                   sol:= SolutionMat( mat, Coefficients( Basis(sp), h0 ) );
551
552                   Add( good_h, h );
553                   Add( sl2s, [sol*f,h0,x] );
554
555                   found:= true;
556
557
558                else
559                   C_h0:= LieCentralizer( g0, Subalgebra( g0, [h0] ) );
560                   sp0:= Subspace( L, List( Basis(C_h0), y -> y*x ) );
561                   if Dimension(sp0) = Length(e) then
562                      found:= true;
563                      good:= false;
564                   fi;
565                fi;
566
567            od;
568
569         fi;
570     od;
571
572     # Now we compute a set of canonical generators of s...
573     inds:= List( fundR, x -> Position( posR_L, x ) );
574
575     x:= PositiveRootVectors( R ){inds};
576     y:= NegativeRootVectors( R ){inds};
577     for i in [1..Length(x)] do
578         V:= VectorSpace( F, [ x[i] ] );
579         b:= Basis( V, [x[i]] );
580         c:= Coefficients( b, (x[i]*y[i])*x[i] )[1];
581         y[i]:= y[i]*2/c;
582     od;
583     bas:= List( [1..Length(x)], i -> x[i]*y[i] );
584
585     Append( bas, BasisVectors( Basis(r) ) );
586     b:= Basis( Subspace( L, bas ), bas );
587
588     # Cartan matrix of s...
589     Cs:= NullMat( Length(fundR), Length(fundR) );
590     B:= BilinearFormMatNF(R);
591     for i in [1..Length(fundR)] do
592         for j in [1..Length(fundR)] do
593             Cs[i][j]:= 2*( fundR[i]*(B*fundR[j]) )/( fundR[j]*(B*fundR[j]) );
594         od;
595     od;
596
597     return rec( hs:= good_h, sl2:= sl2s, chars:= List( good_h, x ->
598                   Cs*( Coefficients( b, x*hL ){[1..Length(x)]} ) ) );
599
600end;
601
602SLAfcts.loop_W:= function( C, h_lst, func )
603
604
605     # C: Cartan matrix
606     # h_lst: list of initial elements of H (given as coefficient vectors,
607     # rel to basis of Chevalley type).
608     # func: function H --> true, false,
609     # if func(orb elt) = true, then orb elt is included...
610
611     local rank, sim, path, h_orb, h, r, i, j, idone, nu, ispos, wrd, hs0;
612
613     rank:= Length( C );
614     sim:= ShallowCopy(C);
615
616     path:= [ rec( wt:= List( [1..rank], x -> 1 ),
617                            word:= [ ],
618                            hs:= h_lst,
619                            ind:= 0 ) ];
620     h_orb:= [ ];
621     for h in h_lst do
622         if func(h) then Add( h_orb, h ); fi;
623     od;
624
625     while Length(path) > 0 do
626
627          r:= path[ Length(path) ];
628          i:= r.ind+1;
629          idone:= false;
630          while i <= rank and not idone do
631                if r.wt[i] <= 0 then
632                   i:= i+1;
633                else
634
635                   nu:= r.wt - r.wt[i]*sim[i];  # i.e. s_i(r.wt)
636                   ispos:= true;
637                   for j in [i+1..rank] do
638                       if nu[j] < 0 then
639                          ispos:= false;
640                          break;
641                       fi;
642                   od;
643                   if ispos then
644                      path[Length(path)]:= rec( wt:= r.wt,
645                                                word:= r.word,
646                                                hs:= r.hs,
647                                                ind:= i );
648                      wrd:= [ i ]; Append( wrd, r.word );
649                      hs0:= ShallowCopy(r.hs);
650                      for j in [1..Length(hs0)] do
651                          h:= ShallowCopy(hs0[j]);
652                          h[i]:= h[i] - C[i]*h;  # i.e., s_i(h)
653                          hs0[j]:= h;
654                      od;
655
656                      Add( path, rec( wt:= nu,
657                                      word:= wrd,
658                                      hs:= hs0,
659                                      ind:= 0 ) );
660                      for h in hs0 do
661                          if func( h ) then
662                             if not h in h_orb then
663                                Add( h_orb, h );
664                             fi;
665                          fi;
666                      od;
667                      idone:= true;
668                   else
669                      i:= i+1;
670                   fi;
671                fi;
672          od;
673          if not idone then  # get rid of last elt as it is searched through
674             Unbind( path[Length(path)] );
675          fi;
676
677     od;
678
679     return h_orb;
680
681end;
682
683
684SLAfcts.nil_orbs_outer:= function( L, gr0, gr1, gr2 )
685
686     # Here L is a simple graded Lie algebra; gr0 a basis of the
687     # elts of degree 0, gr1 of degree 1, and gr2 of degree -1.
688     # We find the nilpotent G_0-orbits in g_1.
689     # We *do not* assume that the given CSA of L is also a CSA of g_0.
690
691     local F, g0, s, r, HL, Hs, R, Ci, hL, hl, C, rank, posRv_L, posR_L,
692           posR, i, j, sums, fundR, inds, tr, h_candidates, BH, W, h,
693           c_h, ph, stb, v, w, is_rep, h0, wr, Omega, good_h, g1, g2, h_mats1,
694           h_mats2, mat, sl2s, id1, id2, V, e, f, bb, ef, found, good, co, x,
695           C_h0, sp, sp0, y, b, bas, c, Cs, B, Rs, nas, b0, ranks, in_weylch,
696           charact, k, sol, info;
697
698     F:= LeftActingDomain(L);
699
700     g0:= Subalgebra( L, gr0, "basis" );
701
702     s:= LieDerivedSubalgebra( g0 );
703     r:= LieCentre(g0);
704
705     HL:= CartanSubalgebra(L);
706     Hs:= Intersection( s, HL );
707     SetCartanSubalgebra( s, Hs );
708
709     R:= RootSystem(L);
710     Ci:= CartanMatrix( R )^-1;
711     hL:= ChevalleyBasis(L)[3];
712
713     hl:= List( NilpotentOrbits(L), x -> Ci*WeightedDynkinDiagram(x) );
714     for i in [1..Length(hl)] do
715         if hl[i] = 0*hl[i] then
716            Unbind( hl[i] );
717         fi;
718     od;
719     hl:= Filtered( hl, x -> IsBound(x) );
720
721     C:= CartanMatrix( R );
722     rank:= Length(C);
723
724     Rs:= RootSystem(s);
725     Cs:= CartanMatrix( Rs );
726     ranks:= Length( Cs );
727
728     bas:= ShallowCopy( CanonicalGenerators(Rs)[3] );
729     Append( bas, BasisVectors( Basis(r) ) );
730     b0:= Basis( VectorSpace( F, bas ), bas );
731
732
733     in_weylch:= function( h )
734
735          local cf, u;
736
737          u:= h*hL;
738          if not u in g0 then return false; fi;
739          cf:= Coefficients( b0, u ){[1..ranks]};
740          if ForAll( Cs*cf, x -> x >= 0 ) then
741             return true;
742          else
743             return false;
744          fi;
745
746     end;
747
748     charact:= function( h )
749
750          local cf;
751
752          cf:= Coefficients( b0, h ){[1..ranks]};
753          return Cs*cf;
754
755     end;
756
757     h_candidates:= SLAfcts.loop_W( C, hl, in_weylch );
758
759info:= "Constructed ";
760Append( info, String(Length(h_candidates)) );
761Append( info, " Cartan elements to be checked.");
762
763Info(InfoSLA,2,info);
764
765     # now we need to compute sl_2 triples wrt the h-s found...
766
767     Omega:= [0..Dimension(L)];
768     good_h:= [ ];
769
770     g1:= Basis( Subspace( L, gr1 ), gr1 );
771     g2:= Basis( Subspace( L, gr2 ), gr2 );
772
773     # the matrices of hL[i] acting on g1
774     h_mats1:= [ ];
775     for h0 in bas do
776         mat:= [ ];
777         for i in [1..Length(g1)] do
778             Add( mat, Coefficients( g1, h0*g1[i] ) );
779         od;
780         Add( h_mats1, mat );
781     od;
782
783     # those of wrt g2...
784     h_mats2:= [ ];
785     for h0 in bas do
786         mat:= [ ];
787         for i in [1..Length(g1)] do
788             Add( mat, Coefficients( g2, h0*g2[i] ) );
789         od;
790         Add( h_mats2, mat );
791     od;
792
793     sl2s:= [ ];
794     id1:= IdentityMat( Length(g1) );
795     id2:= IdentityMat( Length(g2) );
796     for h in h_candidates do
797
798         c_h:= Coefficients( b0, h*hL );
799
800         mat:= c_h*h_mats1;
801         mat:= mat - 2*id1;
802         V:= NullspaceMat( mat );
803         e:= List( V, v -> v*gr1 );
804
805         mat:= c_h*h_mats2;
806         mat:= mat + 2*id2;
807         V:= NullspaceMat( mat );
808         f:= List( V, v -> v*gr2 );
809
810         # check whether h0 in [e,f]....
811         bb:= [ ];
812         for x in e do
813             for y in f do
814                 Add( bb, x*y );
815             od;
816         od;
817         ef:= Subspace( L, bb );
818
819         h0:= h*hL;
820
821         if h0 in ef then  #otherwise we can just discard h...
822            found:= false;
823            good:= false;
824            while not found do
825
826                co:= List( e, x -> Random(Omega) );
827                x:= co*e;
828                sp:= Subspace( L, List( f, y -> x*y) );
829
830                if Dimension(sp) = Length(e) and h0 in sp then
831
832                   # look for a nice one...
833                   for i in [1..Length(co)] do
834                       k:= 0;
835                       found:= false;
836                       while not found do
837                           co[i]:= k;
838                           x:= co*e;
839                           sp:= Subspace( L, List( f, y -> x*y) );
840
841                           if Dimension(sp) = Length(e) and h0 in sp then
842                              found:= true;
843                           else
844                              k:= k+1;
845                           fi;
846                       od;
847                   od;
848
849                   mat:= List( f, u -> Coefficients( Basis(sp), x*u ) );
850                   sol:= SolutionMat( mat, Coefficients( Basis(sp), h0 ) );
851
852                   Add( good_h, h0 );
853                   Add( sl2s, [sol*f,h0,x] );
854
855                   found:= true;
856
857                else
858                   C_h0:= LieCentralizer( g0, Subalgebra( g0, [h0] ) );
859                   sp0:= Subspace( L, List( Basis(C_h0), y -> y*x ) );
860                   if Dimension(sp0) = Length(e) then
861                      found:= true;
862                      good:= false;
863                   fi;
864                fi;
865
866            od;
867
868         fi;
869     od;
870
871     return rec( hs:= good_h, sl2:= sl2s, chars:= List( good_h, charact ) );
872
873end;
874
875SLAfcts.roots_and_vecs:= function( f )
876
877  # we return the roots and corresponding vectors of g_0, and g_1;
878  # the output is a list with two records the first describing
879  # g0, the second describing g1. In the case of g0 the roots are
880  # split in positive/negative.
881
882  local L, R, posR, posRv, negRv, m, vv, g0, g1, pr0, pv0, nr0, nv0,
883        r1, rv1, i, w, m0, gm, rm, rvm, ord2;
884
885  if Order(f) = 2 then ord2:= true; else ord2:= false; fi;
886
887  L:= Source(f);
888  w:= E( Order(f) );
889  R:= RootSystem(L);
890  posR:= PositiveRootsNF(R);
891  posRv:= PositiveRootVectors( R );
892  negRv:= NegativeRootVectors( R );
893
894  m:= List( Basis(L), x -> ShallowCopy( Coefficients( Basis(L), Image(f,x))));
895  m0:= m - IdentityMat( Dimension(L) );
896
897  vv:= NullspaceMat( m0 );
898  vv:= List( vv, x -> LinearCombination( Basis(L), x ) );
899  g0:= Subspace( L, vv, "basis" );
900
901  m0:= m - w*IdentityMat( Dimension(L) );
902  vv:= NullspaceMat( m0 );
903  vv:= List( vv, x -> LinearCombination( Basis(L), x ) );
904  g1:= Subspace( L, vv, "basis" );
905
906  m0:= m - w^(Order(f)-1)*IdentityMat( Dimension(L) );
907  vv:= NullspaceMat( m0 );
908  vv:= List( vv, x -> LinearCombination( Basis(L), x ) );
909  gm:= Subspace( L, vv, "basis" );
910
911
912  pr0:= [ ]; pv0:= [ ];
913  nr0:= [ ]; nv0:= [ ];
914
915  r1:= [ ]; rv1:= [ ];
916
917  rm:= [ ]; rvm:= [ ];
918
919  for i in [1..Length(posR)] do
920      if posRv[i] in g0 then
921         Add( pr0, posR[i] );
922         Add( pv0, posRv[i] );
923         Add( nr0, -posR[i] );
924         Add( nv0, negRv[i] );
925         if not negRv[i] in g0 then Print("OOOOOOOPS!!!!\n"); fi;
926      elif posRv[i] in g1 then
927         Add( r1, posR[i] );
928         Add( rv1, posRv[i] );
929      elif posRv[i] in gm then
930         Add( rm, posR[i] );
931         Add( rvm, posRv[i] );
932      fi;
933      if negRv[i] in g1 then
934         Add( r1, -posR[i] );
935         Add( rv1, negRv[i] );
936      elif negRv[i] in gm then
937         Add( rm, -posR[i] );
938         Add( rvm, negRv[i] );
939      fi;
940  od;
941
942  if ord2 then # g_{-1} = g_{1}....
943      return [ rec( pr0:= pr0, pv0:= pv0, nr0:= nr0, nv0:= nv0 ),
944           rec( r1:= r1, rv1:= rv1 ), rec( rm:= r1, rvm:= rv1 ) ];
945  else
946      return [ rec( pr0:= pr0, pv0:= pv0, nr0:= nr0, nv0:= nv0 ),
947           rec( r1:= r1, rv1:= rv1 ), rec( rm:= rm, rvm:= rvm ) ];
948  fi;
949
950
951end;
952
953SLAfcts.root_basis:= function( B, posr )
954
955  local inds, i, j, pos, bas, C, tp, subs, sub, s, rrr, R, pi, posRw,
956        rts, concs, news, r;
957
958  inds:=[ ];
959  for i in [1..Length(posr)] do
960      for j in [i+1..Length(posr)] do
961          pos:= Position( posr, posr[i]+posr[j] );
962          if pos <> fail then AddSet( inds, pos ); fi;
963      od;
964  od;
965
966  bas:=[ ];
967  for i in [1..Length(posr)] do
968      if not i in inds then
969         Add( bas, posr[i] );
970      fi;
971  od;
972
973  C:=List( bas, x -> [ ] );
974  for i in [1..Length(bas)] do
975      for j in [1..Length(bas)] do
976          C[i][j]:= 2*bas[i]*( B*bas[j] )/( bas[j]*(B*bas[j]) );
977      od;
978  od;
979
980  tp:= CartanType( C );
981
982  subs:=[ ];
983  for i in [1..Length(tp.types)] do
984
985      rrr:= bas{tp.enumeration[i]};
986      R:= RootSystem( tp.types[i] );
987      pi:= SLAfcts.pi_systems( R );
988      sub:= [ ];
989      posRw:= PositiveRootsAsWeights( R );
990      for j in [1..Length( pi.types )] do
991          rts:= pi.roots[j];
992          s:= [ ];
993          for r in rts do
994              pos:= Position( posRw, r );
995              if pos <> fail then
996                 Add( s, PositiveRootsNF(R)[pos]*rrr );
997              else
998                 pos:= Position( posRw, -r );
999                 Add( s, -PositiveRootsNF(R)[pos]*rrr );
1000              fi;
1001          od;
1002          Add( sub, s );
1003      od;
1004      Add( subs, sub );
1005  od;
1006
1007  concs:= [ [ ] ];
1008  for i in [1..Length(subs)] do
1009      news:= [ ];
1010      for s in concs do
1011          for j in [1..Length(subs[i])] do
1012              sub:= ShallowCopy( s );
1013              Append( sub, subs[i][j] );
1014              Add( news, sub );
1015          od;
1016      od;
1017      concs:= news;
1018  od;
1019
1020  return concs;
1021
1022
1023end;
1024
1025
1026
1027SLAfcts.zero_systems:= function( B, posr )
1028
1029  local inds, i, j, pos, bas, C, tp, subs, sub, s, rrr, R, pi, posRw,
1030        rts, concs, news, r;
1031
1032  inds:=[ ];
1033  for i in [1..Length(posr)] do
1034      for j in [i+1..Length(posr)] do
1035          pos:= Position( posr, posr[i]+posr[j] );
1036          if pos <> fail then AddSet( inds, pos ); fi;
1037      od;
1038  od;
1039
1040  bas:=[ ];
1041  for i in [1..Length(posr)] do
1042      if not i in inds then
1043         Add( bas, posr[i] );
1044      fi;
1045  od;
1046
1047  C:=List( bas, x -> [ ] );
1048  for i in [1..Length(bas)] do
1049      for j in [1..Length(bas)] do
1050          C[i][j]:= 2*bas[i]*( B*bas[j] )/( bas[j]*(B*bas[j]) );
1051      od;
1052  od;
1053
1054  tp:= CartanType( C );
1055
1056  subs:=[ ];
1057  for i in [1..Length(tp.types)] do
1058
1059      rrr:= bas{tp.enumeration[i]};
1060      R:= RootSystem( tp.types[i] );
1061      pi:= SLAfcts.sub_systems( R );
1062      sub:= [ [ ] ];
1063      posRw:= PositiveRootsAsWeights( R );
1064      for j in [1..Length( pi.types )] do
1065          rts:= pi.roots[j];
1066          s:= [ ];
1067          for r in rts do
1068              pos:= Position( posRw, r );
1069              if pos <> fail then
1070                 Add( s, PositiveRootsNF(R)[pos]*rrr );
1071              else
1072                 pos:= Position( posRw, -r );
1073                 Add( s, -PositiveRootsNF(R)[pos]*rrr );
1074              fi;
1075          od;
1076          Add( sub, s );
1077      od;
1078      Add( subs, sub );
1079  od;
1080
1081  concs:= [ [ ] ];
1082  for i in [1..Length(subs)] do
1083      news:= [ ];
1084      for s in concs do
1085          for j in [1..Length(subs[i])] do
1086              sub:= ShallowCopy( s );
1087              Append( sub, subs[i][j] );
1088              Add( news, sub );
1089          od;
1090      od;
1091      concs:= news;
1092  od;
1093
1094  return rec( bas:= bas, subs:= concs );
1095
1096
1097end;
1098
1099
1100
1101SLAfcts.my_are_conjugate_0:= function( W, R, B, mus, lams )
1102
1103
1104   # R is the big root system, B the bilin form mat wrt weights,
1105   # mus and lams are lists of weights, we determine whether
1106   # there exists w in W wich w(mus[i]) = lams[i], all i.
1107
1108   local sim, i, j, k, a, b, w, mu, rmu;
1109
1110   sim:= SimpleRootsAsWeights( R );
1111
1112   for i in [1..Length(mus)] do
1113
1114       rmu:= List( W.roots, x -> 2*mus[i]*( B*x )/( x*(B*x) ) );
1115       a:= SLAfcts.conj_dom_wt( mus[i], rmu, W );
1116
1117       rmu:= List( W.roots, x -> 2*lams[i]*( B*x )/( x*(B*x) ) );
1118       b:= SLAfcts.conj_dom_wt( lams[i], rmu, W );
1119
1120       if a[1] <> b[1] then return false; fi;
1121
1122       w:= Reversed( b[3] );
1123       Append( w, a[3] );
1124       w:= Reversed(w);
1125
1126       for k in [i..Length(mus)] do
1127
1128           mu:= ShallowCopy(mus[k]);
1129           rmu:= List( W.roots, x -> 2*mu*( B*x )/( x*(B*x) ) );
1130
1131           for j in w do
1132
1133               mu:= mu -rmu[j]*W.roots[j];
1134               rmu:= rmu - rmu[j]*W.wgts[j];
1135
1136           od;
1137
1138           mus[k]:= mu;
1139
1140       od;
1141
1142       W:= SLAfcts.stabilizer( lams[i], B, W );
1143
1144   od;
1145
1146   return true;
1147
1148
1149end;
1150
1151
1152SLAfcts.my_are_conjugate:= function( W, R, B, mus, lams )
1153
1154
1155   # same as previous function, but now we also permute
1156   # the mus, lams. We do assume that they arrive in an
1157   # order that defines the same Cartan matrix...
1158
1159   local C, perms, i, newperms, p, q, good, j, k, l, nus;
1160
1161   # however,... first we try the identity permutation...
1162
1163   if SLAfcts.my_are_conjugate_0( W, R, B, mus, lams ) then
1164      return true;
1165   fi;
1166
1167   # The Cartan matrix:
1168   C:= List( mus, x -> List( mus, y -> 2*x*(B*y)/( y*(B*y) ) ) );
1169
1170   # Now we determine all permutations of the mus that leave C invariant:
1171
1172   perms:= List( [1..Length(mus)], x -> [x] );
1173   # i.e., the first element can be mapped to one of the other elts.
1174   # now we see whether this can be extended...
1175
1176   for i in [2..Length(mus)] do
1177
1178       newperms:= [ ];
1179       for p in perms do
1180           for j in [1..Length(mus)] do
1181               # see whether p can be extended by adding j...
1182               if not j in p then
1183                  q:= ShallowCopy(p);
1184                  Add( q, j );
1185                  good:= true;
1186                  for k in [1..i] do
1187                      if not good then break; fi;
1188                      for l in [1..i] do
1189                          if not good then break; fi;
1190                          if C[k][l] <> C[ q[k] ][ q[l] ] then
1191                             good:= false;
1192                          fi;
1193                      od;
1194                  od;
1195                  if good then Add( newperms, q ); fi;
1196               fi;
1197           od;
1198       od;
1199       perms:= newperms;
1200   od;
1201
1202   perms:= Filtered( perms, x -> x <> [1..Length(mus)] ); # already tried it
1203
1204   # now we see whether there is a permutation mapping
1205   # a permuted mus to lams...
1206
1207   for p in perms do
1208
1209       nus:= [ ];
1210       for i in [1..Length(p)] do
1211           nus[p[i]]:= mus[i];
1212       od;
1213
1214       if SLAfcts.my_are_conjugate_0( W, R, B, nus, lams ) then
1215          return true;
1216       fi;
1217   od;
1218
1219   return false;
1220
1221end;
1222
1223
1224
1225SLAfcts.inner_orbits_carrier:= function( f )
1226
1227   # we give a list of all flat Z-graded subalgebras of the
1228   # graded Lie algebra corresponding to f.
1229
1230
1231   local L, R, B, ch, posR, N, rts, rr, pi, r1, zero, stack, res, r,
1232         start, rrr, ips, i, vv, u, h, C, CT, pi_0, pi_1, t, s, pos,
1233         ct, eqns, rhs, eqn, j, sol, h0, psi0, psi1, good, x, y, es, fs,
1234         valmat, val, chars, u0, v, done, gr1, gr2, g1, g2, h_mats1, h_mats2,
1235         mat, sl2s, id1, id2, Omega, V, e, ff, found, co, k, sp, extended,
1236         zz, bas, sim, Bw, W0, types, weights, wrts, tp, a, c, comb, hZ, hs,
1237         info;
1238
1239
1240   L:= Source(f);
1241
1242   ch:= ChevalleyBasis(L);
1243
1244   R:= RootSystem(L);
1245
1246   posR:= PositiveRootsNF(R);
1247   N:= Length( posR );
1248   rts:= ShallowCopy(posR);
1249   Append( rts, -posR );
1250
1251   B:= BilinearFormMatNF(R);
1252
1253   rr:= SLAfcts.roots_and_vecs( f );
1254
1255   zz:= SLAfcts.zero_systems( B, rr[1].pr0 );
1256   pi:= zz.subs;
1257
1258   # now see how we can extend each element in pi with roots of
1259   # weight 1... and compute the maximal ones first!
1260
1261   bas:= zz.bas;
1262   sim:= [ ];
1263   for a in bas do
1264       pos:= Position( posR, a );
1265       Add( sim, PositiveRootsAsWeights( R )[pos] );
1266   od;
1267
1268   Bw:= SLAfcts.bilin_weights( R );
1269   W0:= rec( roots:= sim, wgts:= List( sim, x -> List( sim, y ->
1270                   2*x*(Bw*y)/( y*(Bw*y) ) ) ) );
1271
1272
1273   r1:= rr[2].r1;
1274
1275   zero:= 0*r1[1];
1276
1277   res:= [ ];
1278   for k in [1..Length(pi)] do
1279
1280       types:= [ ];
1281       weights:= [ ];
1282
1283       stack:= [ rec( rts0:= pi[k], rts1:= [ ], start:= 0,
1284                      sp:= VectorSpace( Rationals, pi[k], zero ) ) ];
1285       while Length(stack) > 0 do
1286           r:= stack[Length(stack)];
1287           RemoveElmList( stack, Length(stack) );
1288           start:= r.start+1;
1289           rrr:= Concatenation( r.rts0, r.rts1 );
1290           extended:= false;
1291           for i in [start..Length(r1)] do
1292               ips:= List( rrr, x -> x - r1[i] );
1293               if ForAll( ips, x -> not ( x in rts ) ) and
1294                           not r1[i] in r.sp then
1295                  vv:= ShallowCopy( BasisVectors( Basis(r.sp) ) );
1296                  Add( vv, r1[i] );
1297                  u:= ShallowCopy( r.rts1 );
1298                  Add( u, r1[i] );
1299                  Add( stack, rec( rts0:= r.rts0, rts1:= u, start:= i,
1300                          sp:= VectorSpace( Rationals, vv ) ) );
1301                  extended:= true;
1302               fi;
1303           od;
1304           if not extended then # see whether we can extend by
1305                                # adding something "smaller"
1306              for i in [1..start-1] do
1307                  if not r1[i] in rrr then
1308                     ips:= List( rrr, x -> x - r1[i] );
1309                     if ForAll( ips, x -> not ( x in rts ) ) and
1310                                    not r1[i] in r.sp then
1311                        extended:= true; break;
1312                     fi;
1313                  fi;
1314              od;
1315           fi;
1316
1317           if not extended then
1318              C:= List( rrr, x -> List( rrr, y -> 2*x*(B*y)/(y*(B*y)) ) );
1319              tp:= CartanType( C );
1320              SortParallel( tp.types, tp.enumeration );
1321              wrts:= [ ];
1322              for i in [1..Length(tp.enumeration)] do
1323                  for j in tp.enumeration[i] do
1324                      pos:= Position( rts, rrr[j] );
1325                      if pos <= N then
1326                         Add( wrts, PositiveRootsAsWeights(R)[pos] );
1327                      else
1328                         Add( wrts, -PositiveRootsAsWeights(R)[pos-N] );
1329                      fi;
1330                  od;
1331              od;
1332              found:= false;
1333              if tp.types in types then
1334                 for i in [1..Length(types)] do
1335                     if tp.types = types[i] then
1336                        if SLAfcts.my_are_conjugate( W0, R, Bw, wrts, weights[i] ) then
1337                           found:= true;
1338                           break;
1339                        fi;
1340                     fi;
1341                 od;
1342              fi;
1343              if not found then
1344                 Add( types, tp.types );
1345                 Add( weights, wrts );
1346                 Add( res, r );
1347              fi;
1348           fi;
1349       od;
1350
1351   od;
1352
1353   stack:= [ ];
1354   for r in res do
1355
1356       comb:= Combinations( [1..Length(r.rts1)] );
1357       comb:= Filtered( comb, x -> x <> [ ] );
1358       for c in comb do
1359           Add( stack, rec( rts0:= r.rts0, rts1:= r.rts1{c} ) );
1360       od;
1361
1362   od;
1363
1364   res:= stack;
1365
1366info:= "Constructed ";
1367Append( info, String(Length(res)) );
1368Append( info, " root bases of possible flat subalgebras, now checking them...");
1369Info( InfoSLA, 2, info );
1370
1371   h:= BasisVectors( Basis( CartanSubalgebra(L) ) );
1372
1373   C:= CartanMatrix(R);
1374   CT:= TransposedMat( C );
1375
1376   # HERE we assume inner!
1377
1378   good:= [ ];
1379   for r in res do
1380
1381       pi_0:= r.rts0;
1382       pi_1:= r.rts1;
1383
1384       t:= [ ];
1385       pi:= Concatenation( pi_0, pi_1 );
1386       for s in pi do
1387           pos:= Position( rts, s );
1388           if pos <= N then
1389              Add( t, ch[1][pos]*ch[2][pos] );
1390           else
1391              Add( t, ch[2][pos-N]*ch[1][pos-N] );
1392           fi;
1393       od;
1394
1395       t:= BasisVectors( Basis( Subspace( L, t ) ) );
1396
1397       ct:= List( t, x -> Coefficients( Basis(CartanSubalgebra(L)), x ) );
1398
1399       # i.e. t is a Cartan subalgebra of s
1400
1401       # find h0 in t such that a(h0)=1 for all a in pi_1, a(h0)=0
1402       # for all a in pi_0
1403
1404       eqns:=[ ];
1405       rhs:= [ ];
1406       for j in [1..Length(pi_0)] do
1407           eqn:= [ ];
1408           for i in [1..Length(t)] do
1409               eqn[i]:= pi_0[j]*( C*ct[i] );
1410           od;
1411           Add( eqns, eqn ); Add( rhs, 0 );
1412       od;
1413       for j in [1..Length(pi_1)] do
1414           eqn:= [ ];
1415           for i in [1..Length(t)] do
1416               eqn[i]:= pi_1[j]*( C*ct[i] );
1417           od;
1418           Add( eqns, eqn ); Add( rhs, 1 );
1419       od;
1420
1421       sol:= SolutionMat( TransposedMat(eqns), rhs );
1422       h0:= sol*t;
1423
1424       # Find a basis of the subspace of h consisting of u with
1425       # a(u) = 0, for a in pi = pi_0 \cup pi_1.
1426
1427       eqns:= [ ];
1428       for i in [1..Length(h)] do
1429           eqns[i]:= [ ];
1430           for j in [1..Length(pi_0)] do
1431               Add( eqns[i], pi_0[j]*CT[i] );
1432           od;
1433           for j in [1..Length(pi_1)] do
1434               Add( eqns[i], pi_1[j]*CT[i] );
1435           od;
1436       od;
1437       sol:= NullspaceMat( eqns );
1438       hZ:= List( sol, u -> u*h );
1439
1440       # Now we compute |Psi_0| and |Psi_1|...
1441
1442       psi0:= [ ];
1443       for a in rr[1].pv0 do
1444           if h0*a = 0*a and ForAll( hZ, u -> u*a = 0*a ) then
1445              Add( psi0, a );
1446           fi;
1447       od;
1448
1449       psi1:= [ ];
1450       for a in rr[2].rv1 do
1451           if h0*a = a and ForAll( hZ, u -> u*a = 0*a ) then
1452              Add( psi1, a );
1453           fi;
1454       od;
1455
1456       if Length(pi_0)+Length(pi_1) + 2*Length(psi0) = Length(psi1) then
1457
1458          if not 2*h0 in good then
1459             Add( good, 2*h0 );
1460          fi;
1461
1462       fi;
1463   od;
1464
1465info:= "Obtained ";
1466Append( info, String( Length(good) ) );
1467Append( info, " Cartan elements, weeding out equivalent copies...");
1468Info(InfoSLA,2,info);
1469
1470# NEXT can be obtained from Kac diagram!!
1471
1472   x:= ChevalleyBasis(L)[1];
1473   y:= ChevalleyBasis(L)[2];
1474   es:= [ ];
1475   fs:= [ ];
1476   if Image( f, y[Length(y)] ) = y[Length(y)] then
1477      Add( fs, x[Length(x)] );
1478      Add( es, y[Length(y)] );
1479   fi;
1480   for i in [1..Length(CartanMatrix(R))] do
1481       if Image( f, x[i] ) = x[i] then
1482          Add( es, x[i] );
1483          Add( fs, y[i] );
1484       fi;
1485   od;
1486   hs:= List( [1..Length(es)], i -> es[i]*fs[i] );
1487
1488   valmat:= [ ];
1489   for i in [1..Length(hs)] do
1490       val:= [ ];
1491       for j in [1..Length(hs)] do
1492           Add( val, Coefficients( Basis( Subspace(L,[es[j]]), [es[j]] ),
1493                       hs[i]*es[j] )[1] );
1494       od;
1495       Add( valmat, val );
1496   od;
1497
1498
1499   chars:= [ ];
1500   for i in [1..Length(good)] do
1501
1502       u0:= good[i];
1503       v:= List( es, z -> Coefficients( Basis(Subspace(L,[z]),[z]), u0*z )[1] );
1504       done:= ForAll( v, z -> z >= 0 );
1505       while not done do
1506           pos:= PositionProperty( v, z -> z < 0 );
1507           u0:= u0 - v[pos]*hs[pos];
1508           v:= v - v[pos]*valmat[pos];
1509           done:= ForAll( v, z -> z >= 0 );
1510       od;
1511
1512       if not u0 in chars then
1513          Add( chars, u0 );
1514       fi;
1515   od;
1516
1517
1518   gr1:= rr[2].rv1;
1519   gr2:= rr[3].rvm;
1520
1521     g1:= Basis( Subspace( L, gr1 ), gr1 );
1522     g2:= Basis( Subspace( L, gr2 ), gr2 );
1523
1524     # the matrices of hL[i] acting on g1
1525     h_mats1:= [ ];
1526     for h0 in h do
1527         mat:= [ ];
1528         for i in [1..Length(g1)] do
1529             Add( mat, Coefficients( g1, h0*g1[i] ) );
1530         od;
1531         Add( h_mats1, mat );
1532     od;
1533
1534     # those of wrt g2...
1535     h_mats2:= [ ];
1536     for h0 in h do
1537         mat:= [ ];
1538         for i in [1..Length(g1)] do
1539             Add( mat, Coefficients( g2, h0*g2[i] ) );
1540         od;
1541         Add( h_mats2, mat );
1542     od;
1543
1544     sl2s:= [ ];
1545     id1:= IdentityMat( Length(g1) );
1546     id2:= IdentityMat( Length(g2) );
1547     Omega:= [1..Dimension(L)];
1548     for h0 in chars do
1549
1550         ch:= Coefficients( Basis( CartanSubalgebra(L) ), h0 );
1551         mat:= ch*h_mats1;
1552         mat:= mat - 2*id1;
1553         V:= NullspaceMat( mat );
1554         e:= List( V, v -> v*gr1 );
1555
1556         mat:= ch*h_mats2;
1557         mat:= mat + 2*id2;
1558         V:= NullspaceMat( mat );
1559         ff:= List( V, v -> v*gr2 );
1560
1561         found:= false;
1562         while not found do
1563
1564             co:= List( e, x -> Random(Omega) );
1565             x:= co*e;
1566             sp:= Subspace( L, List( ff, y -> x*y) );
1567
1568             if Dimension(sp) = Length(e) and h0 in sp then
1569
1570                # look for a nice one...
1571                for i in [1..Length(co)] do
1572                    k:= 0;
1573                    found:= false;
1574                    while not found do
1575                        co[i]:= k;
1576                        x:= co*e;
1577                        sp:= Subspace( L, List( ff, y -> x*y) );
1578
1579                        if Dimension(sp) = Length(e) and h0 in sp then
1580                           found:= true;
1581                        else
1582                           k:= k+1;
1583                        fi;
1584                    od;
1585                od;
1586
1587                mat:= List( ff, u -> Coefficients( Basis(sp), x*u ) );
1588                sol:= SolutionMat( mat, Coefficients( Basis(sp), h0 ) );
1589
1590                Add( sl2s, [sol*ff,h0,x] );
1591
1592                found:= true;
1593
1594
1595             fi;
1596
1597         od;
1598
1599     od;
1600
1601   return sl2s;
1602
1603end;
1604
1605
1606
1607
1608InstallMethod( NilpotentOrbitsOfThetaRepresentation,
1609"for a finite order automorphism", true, [ IsGeneralMapping ], 0,
1610function( f )
1611
1612   local g, L, rank, r, meth, kd, C, inds, i, w, tr;
1613
1614   g:= Grading(f);
1615
1616   if g[2] = [ ] then return [ ]; fi;
1617
1618   meth:= ValueOption( "method" );
1619
1620   L:= Source(f);
1621   rank:= Length( CartanMatrix( RootSystem(L) ) );
1622   if Length( KacDiagram( f ).weights ) = rank +1 then
1623
1624      if meth = fail then
1625         kd:= KacDiagram( f );
1626         C:= kd.CM;
1627         inds:= [ ];
1628         for i in [1..Length(kd.weights)] do
1629             if kd.weights[i] = 0 then Add( inds, i ); fi;
1630         od;
1631         w:= SizeOfWeylGroup( CartanType( C{inds}{inds} ).types );
1632         tr:= SizeOfWeylGroup( RootSystem(L) )/w;
1633         if tr > 8000 then
1634            meth:= "Carrier";
1635         else
1636            meth:= "WeylOrbit";
1637         fi;
1638      fi;
1639
1640      if meth = "WeylOrbit" then
1641         Info(InfoSLA,2,"Selected Weyl orbit method.");
1642         r:= SLAfcts.nil_orbs_inner( L, g[1], g[2], g[Length(g)] ).sl2;
1643      else
1644         Info(InfoSLA,2,"Selected carrier algebra method.");
1645         r:= SLAfcts.inner_orbits_carrier( f );
1646      fi;
1647   else
1648      r:= SLAfcts.nil_orbs_outer( L, g[1], g[2], g[Length(g)] ).sl2;
1649   fi;
1650
1651   return r;
1652
1653
1654
1655end );
1656
1657SLAfcts.CartanMatrixToPositiveRoots:= function( C )
1658
1659        local   rank,  posr,  ready,  ind,  le,  i,  a,  j,  ej,  r,  b,
1660                q, CT;
1661
1662        rank:= Length( C );
1663        CT:= TransposedMat(C);
1664
1665        # posr will be a list of the positive roots. We start with the
1666        # simple roots, which are simply unit vectors.
1667
1668        posr:= IdentityMat( rank );
1669
1670        ready:= false;
1671        ind:= 1;
1672        le:= rank;
1673        while ind <= le  do
1674
1675            # We loop over those elements of posR that have been found in
1676            # the previous round, i.e., those at positions ranging from
1677            # ind to le.
1678
1679            le:= Length( posr );
1680            for i in [ind..le] do
1681                a:= posr[i];
1682
1683                # We determine whether a+ej is a root (where ej is the j-th
1684                # simple root.
1685                for j in [1..rank] do
1686                    ej:= posr[j];
1687
1688                    # We determine the maximum number r such that a-r*ej is
1689                    # a root.
1690                    r:= -1;
1691                    b:= ShallowCopy( a );
1692                    while b in posr do
1693                        b:= b-ej;
1694                        r:=r+1;
1695                    od;
1696                    q:= r-LinearCombination( CT[j], a );
1697                    if q>0 and (not a+ej in posr ) then
1698                        Add( posr, a+ej );
1699                    fi;
1700                od;
1701            od;
1702            ind:= le+1;
1703            le:= Length( posr );
1704        od;
1705
1706        return posr;
1707    end;
1708
1709
1710
1711SLAfcts.sub_systems_Delta:= function( R )
1712
1713   # simple root system..., we give reps of all orbits of
1714   # sub root systems that have a basis which is a subset of the basis of R,
1715   # under the Weyl group
1716
1717   local pis, B, roots, types, tps, rts, mus, pos, found, i, j, k, comb,
1718         r0, c, C, r1, tp, e, u, t1, rank;
1719
1720   tp:= CartanType( CartanMatrix( R ) );
1721   pis:= rec( types:= [tp.types], roots:= [SimpleRootsAsWeights(R){tp.enumeration[1]}] );
1722   B:= SLAfcts.bilin_weights( R );
1723
1724   roots:= [ ];
1725   types:= [ ];
1726
1727   rank:= Length(B);
1728   comb:= Combinations( [1..rank] );
1729   comb:= Filtered( comb, x -> (x <> [] and Length(x) <> rank ) );
1730
1731
1732   for i in [1..Length(pis.types)] do
1733       tps:= pis.types[i];
1734       rts:= pis.roots[i];
1735
1736       Add( roots, rts );
1737       Add( types, tps );
1738
1739
1740       for c in comb do
1741
1742           r0:= rts{c};
1743           # find its type in normal enumeration...
1744
1745           C:= List( r0, x -> List( r0, y -> 2*x*(B*y)/(y*(B*y)) ) );
1746           tp:= CartanType( C );
1747
1748           e:= tp.enumeration;
1749           r1:= [ ];
1750           for j in [1..Length(e)] do
1751               u:= [ ];
1752               for k in e[j] do
1753                   Add( u, r0[k] );
1754               od;
1755               Add( r1, u );
1756           od;
1757
1758           t1:= tp.types;
1759           SortParallel( t1, r1 );
1760
1761           mus:= Concatenation( r1 );
1762
1763           pos:= Position( types, t1 );
1764           if pos = fail then
1765              Add( types, t1 );
1766              Add( roots, mus );
1767           else
1768
1769              found:= false;
1770              for j in [pos..Length(types)] do
1771                  if types[j] = t1 then
1772
1773                     if SLAfcts.are_conjugate( R, B, mus, roots[j] ) then
1774                        found:= true; break;
1775                     fi;
1776
1777                  fi;
1778              od;
1779              if not found then
1780                 Add( types, t1 );
1781                 Add( roots, mus );
1782              fi;
1783
1784           fi;
1785       od;
1786   od;
1787
1788   return rec( types:= types, roots:= roots );
1789
1790
1791end;
1792
1793
1794
1795
1796SLAfcts.roots_and_vecs_Z:= function( L, g0,g1,gm )
1797
1798  # we return the roots and corresponding vectors of g_0, and g_1;
1799  # the output is a list with two records the first describing
1800  # g0, the second describing g1. In the case of g0 the roots are
1801  # split in positive/negative.
1802
1803  local R, posR, posRv, negRv, m, vv, pr0, pv0, nr0, nv0,
1804        r1, rv1, i, rm, rvm;
1805
1806  R:= RootSystem(L);
1807  posR:= PositiveRootsNF(R);
1808  posRv:= PositiveRootVectors( R );
1809  negRv:= NegativeRootVectors( R );
1810
1811  pr0:= [ ]; pv0:= [ ];
1812  nr0:= [ ]; nv0:= [ ];
1813
1814  r1:= [ ]; rv1:= [ ];
1815
1816  rm:= [ ]; rvm:= [ ];
1817
1818  for i in [1..Length(posR)] do
1819      if posRv[i] in g0 then
1820         Add( pr0, posR[i] );
1821         Add( pv0, posRv[i] );
1822         Add( nr0, -posR[i] );
1823         Add( nv0, negRv[i] );
1824         if not negRv[i] in g0 then Print("OOOOOOOPS!!!!\n"); fi;
1825      elif posRv[i] in g1 then
1826         Add( r1, posR[i] );
1827         Add( rv1, posRv[i] );
1828      elif posRv[i] in gm then
1829         Add( rm, posR[i] );
1830         Add( rvm, posRv[i] );
1831      fi;
1832      if negRv[i] in g1 then
1833         Add( r1, -posR[i] );
1834         Add( rv1, negRv[i] );
1835      elif negRv[i] in gm then
1836         Add( rm, -posR[i] );
1837         Add( rvm, negRv[i] );
1838      fi;
1839  od;
1840
1841  return [ rec( pr0:= pr0, pv0:= pv0, nr0:= nr0, nv0:= nv0 ),
1842           rec( r1:= r1, rv1:= rv1 ), rec( rm:= rm, rvm:= rvm ) ];
1843
1844
1845end;
1846
1847
1848SLAfcts.zero_systems_Z:= function( B, posr )
1849
1850  local inds, i, j, pos, bas, C, tp, subs, sub, s, rrr, R, pi, posRw,
1851        rts, concs, news, r;
1852
1853  inds:=[ ];
1854  for i in [1..Length(posr)] do
1855      for j in [i+1..Length(posr)] do
1856          pos:= Position( posr, posr[i]+posr[j] );
1857          if pos <> fail then AddSet( inds, pos ); fi;
1858      od;
1859  od;
1860
1861  bas:=[ ];
1862  for i in [1..Length(posr)] do
1863      if not i in inds then
1864         Add( bas, posr[i] );
1865      fi;
1866  od;
1867
1868  C:=List( bas, x -> [ ] );
1869  for i in [1..Length(bas)] do
1870      for j in [1..Length(bas)] do
1871          C[i][j]:= 2*bas[i]*( B*bas[j] )/( bas[j]*(B*bas[j]) );
1872      od;
1873  od;
1874
1875  tp:= CartanType( C );
1876
1877  subs:=[ ];
1878  for i in [1..Length(tp.types)] do
1879
1880      rrr:= bas{tp.enumeration[i]};
1881      R:= RootSystem( tp.types[i] );
1882      pi:= SLAfcts.sub_systems_Delta( R );
1883      sub:= [ [ ] ];
1884      posRw:= PositiveRootsAsWeights( R );
1885      for j in [1..Length( pi.types )] do
1886          rts:= pi.roots[j];
1887          s:= [ ];
1888          for r in rts do
1889              pos:= Position( posRw, r );
1890              if pos <> fail then
1891                 Add( s, PositiveRootsNF(R)[pos]*rrr );
1892              else
1893                 pos:= Position( posRw, -r );
1894                 Add( s, -PositiveRootsNF(R)[pos]*rrr );
1895              fi;
1896          od;
1897          Add( sub, s );
1898      od;
1899      Add( subs, sub );
1900  od;
1901
1902  concs:= [ [ ] ];
1903  for i in [1..Length(subs)] do
1904      news:= [ ];
1905      for s in concs do
1906          for j in [1..Length(subs[i])] do
1907              sub:= ShallowCopy( s );
1908              Append( sub, subs[i][j] );
1909              Add( news, sub );
1910          od;
1911      od;
1912      concs:= news;
1913  od;
1914
1915  return rec( bas:= bas, subs:= concs );
1916
1917
1918end;
1919
1920
1921# NOTE: basis of simple roots in g0 directly from grading-diagram!
1922
1923
1924SLAfcts.zgrad_orbits_carrier:= function( L, grading )
1925
1926   # L: Lie algebra, gr: grading (0,1,-1 components).
1927   #
1928
1929
1930   local R, B, ch, posR, N, rts, rr, pi, r1, zero, stack, res, r,
1931         start, rrr, ips, i, vv, u, h, C, CT, pi_0, pi_1, t, s, pos,
1932         ct, eqns, rhs, eqn, j, sol, h0, psi0, psi1, good, x, y, es, fs,
1933         valmat, val, chars, u0, v, done, gr1, gr2, g2, h_mats1, h_mats2,
1934         mat, sl2s, id1, id2, Omega, V, e, ff, found, co, k, sp, extended,
1935         zz, bas, sim, Bw, W0, types, weights, wrts, tp, a, c, comb, hZ, hs,
1936         info, posRv, negRv, g0, g1, gm, CM, rr0, l0, l1, gr, deg;
1937
1938
1939   ch:= ChevalleyBasis(L);
1940
1941   R:= RootSystem(L);
1942
1943   posR:= PositiveRootsNF(R);
1944   posRv:= PositiveRootVectors(R);
1945   negRv:= NegativeRootVectors(R);
1946   N:= Length( posR );
1947   rts:= ShallowCopy(posR);
1948   Append( rts, -posR );
1949
1950   B:= BilinearFormMatNF(R);
1951
1952   rr:= [ rec( pr0:= [ ], pv0:= [ ], nv0:= [] ), rec( r1:= [ ], rv1:= [ ] ), rec( rvm:= [ ] ) ];
1953   for i in [1..Length(posR)] do
1954         v:= posR[i]*grading;
1955         if v = 0 then
1956            Add( rr[1].pr0, posR[i] );
1957            Add( rr[1].pv0, posRv[i] );
1958            Add( rr[1].nv0, negRv[i] );
1959         elif v = 1 then
1960            Add( rr[2].r1, posR[i] );
1961            Add( rr[2].rv1, posRv[i] );
1962            Add( rr[3].rvm, negRv[i] );
1963         fi;
1964   od;
1965
1966   zz:= SLAfcts.zero_systems_Z( B, rr[1].pr0 );
1967   pi:= zz.subs;
1968
1969   # now see how we can extend each element in pi with roots of
1970   # weight 1... and compute the maximal ones first!
1971
1972   bas:= zz.bas;
1973   sim:= [ ];
1974   for a in bas do
1975       pos:= Position( posR, a );
1976       Add( sim, PositiveRootsAsWeights( R )[pos] );
1977   od;
1978
1979   Bw:= SLAfcts.bilin_weights( R );
1980   W0:= rec( roots:= sim, wgts:= List( sim, x -> List( sim, y ->
1981                   2*x*(Bw*y)/( y*(Bw*y) ) ) ) );
1982
1983
1984   r1:= rr[2].r1;
1985
1986   zero:= 0*r1[1];
1987
1988   res:= [ ];
1989   for k in [1..Length(pi)] do
1990
1991       types:= [ ];
1992       weights:= [ ];
1993
1994       stack:= [ rec( rts0:= pi[k], rts1:= [ ], start:= 0,
1995                      sp:= VectorSpace( Rationals, pi[k], zero ) ) ];
1996       while Length(stack) > 0 do
1997           r:= stack[Length(stack)];
1998           RemoveElmList( stack, Length(stack) );
1999           start:= r.start+1;
2000           rrr:= Concatenation( r.rts0, r.rts1 );
2001           extended:= false;
2002           for i in [start..Length(r1)] do
2003               ips:= List( rrr, x -> x - r1[i] );
2004               if ForAll( ips, x -> not ( x in rts ) ) and
2005                           not r1[i] in r.sp then
2006                  vv:= ShallowCopy( BasisVectors( Basis(r.sp) ) );
2007                  Add( vv, r1[i] );
2008                  u:= ShallowCopy( r.rts1 );
2009                  Add( u, r1[i] );
2010                  Add( stack, rec( rts0:= r.rts0, rts1:= u, start:= i,
2011                          sp:= VectorSpace( Rationals, vv ) ) );
2012                  extended:= true;
2013               fi;
2014           od;
2015           if not extended then # see whether we can extend by
2016                                # adding something "smaller"
2017              for i in [1..start-1] do
2018                  if not r1[i] in rrr then
2019                     ips:= List( rrr, x -> x - r1[i] );
2020                     if ForAll( ips, x -> not ( x in rts ) ) and
2021                                    not r1[i] in r.sp then
2022                        extended:= true; break;
2023                     fi;
2024                  fi;
2025              od;
2026           fi;
2027
2028           if not extended then
2029              C:= List( rrr, x -> List( rrr, y -> 2*x*(B*y)/(y*(B*y)) ) );
2030              tp:= CartanType( C );
2031              SortParallel( tp.types, tp.enumeration );
2032              wrts:= [ ];
2033              for i in [1..Length(tp.enumeration)] do
2034                  for j in tp.enumeration[i] do
2035                      pos:= Position( rts, rrr[j] );
2036                      if pos <= N then
2037                         Add( wrts, PositiveRootsAsWeights(R)[pos] );
2038                      else
2039                         Add( wrts, -PositiveRootsAsWeights(R)[pos-N] );
2040                      fi;
2041                  od;
2042              od;
2043              found:= false;
2044              if tp.types in types then
2045                 for i in [1..Length(types)] do
2046                     if tp.types = types[i] then
2047                        if SLAfcts.my_are_conjugate( W0, R, Bw, wrts, weights[i] ) then
2048                           found:= true;
2049                           break;
2050                        fi;
2051                     fi;
2052                 od;
2053              fi;
2054              if not found then
2055                 Add( types, tp.types );
2056                 Add( weights, wrts );
2057                 Add( res, r );
2058              fi;
2059           fi;
2060       od;
2061
2062   od;
2063
2064   stack:= [ ];
2065   for r in res do
2066
2067       comb:= Combinations( [1..Length(r.rts1)] );
2068       comb:= Filtered( comb, x -> x <> [ ] );
2069       for c in comb do
2070           Add( stack, rec( rts0:= r.rts0, rts1:= r.rts1{c} ) );
2071       od;
2072
2073   od;
2074
2075   res:= stack;
2076
2077info:= "Constructed ";
2078Append( info, String(Length(res)) );
2079Append( info, " root bases of possible flat subalgebras, now checking them...");
2080Info( InfoSLA, 2, info );
2081
2082   h:= BasisVectors( Basis( CartanSubalgebra(L) ) );
2083
2084   C:= CartanMatrix(R);
2085   CT:= TransposedMat( C );
2086
2087   good:= [ ];
2088   for r in res do
2089
2090       pi_0:= r.rts0;
2091       pi_1:= r.rts1;
2092       pi:= Concatenation( pi_0, pi_1 );
2093
2094       CM:= List( pi, x -> List( pi, y -> 2*x*(B*y)/( y*(B*y) ) ) );
2095       rr0:= SLAfcts.CartanMatrixToPositiveRoots( CM );
2096       l0:= 0; l1:= 0;
2097       gr:= Concatenation( List( pi_0, x -> 0 ), List( pi_1, x -> 1 ) );
2098       for s in rr0 do
2099           deg:= s*gr;
2100           if deg=0 then
2101              l0:= l0+1;
2102           elif deg=1 then
2103              l1:= l1+1;
2104           fi;
2105       od;
2106
2107       if 2*l0+Length(pi) = l1 then
2108
2109          t:= [ ];
2110          for s in pi do
2111              pos:= Position( rts, s );
2112              if pos <= N then
2113                 Add( t, ch[1][pos]*ch[2][pos] );
2114              else
2115                 Add( t, ch[2][pos-N]*ch[1][pos-N] );
2116              fi;
2117          od;
2118
2119          t:= BasisVectors( Basis( Subspace( L, t ) ) );
2120
2121          ct:= List( t, x -> Coefficients( Basis(CartanSubalgebra(L)), x ) );
2122
2123          # i.e. t is a Cartan subalgebra of s
2124
2125          # find h0 in t such that a(h0)=1 for all a in pi_1, a(h0)=0
2126          # for all a in pi_0
2127
2128          eqns:=[ ];
2129          rhs:= [ ];
2130          for j in [1..Length(pi_0)] do
2131              eqn:= [ ];
2132              for i in [1..Length(t)] do
2133                  eqn[i]:= pi_0[j]*( C*ct[i] );
2134              od;
2135              Add( eqns, eqn ); Add( rhs, 0 );
2136          od;
2137          for j in [1..Length(pi_1)] do
2138              eqn:= [ ];
2139              for i in [1..Length(t)] do
2140                  eqn[i]:= pi_1[j]*( C*ct[i] );
2141              od;
2142              Add( eqns, eqn ); Add( rhs, 1 );
2143          od;
2144
2145          sol:= SolutionMat( TransposedMat(eqns), rhs );
2146          h0:= sol*t;
2147
2148          # Find a basis of the subspace of h consisting of u with
2149          # a(u) = 0, for a in pi = pi_0 \cup pi_1.
2150
2151          eqns:= [ ];
2152          for i in [1..Length(h)] do
2153              eqns[i]:= [ ];
2154              for j in [1..Length(pi_0)] do
2155                  Add( eqns[i], pi_0[j]*CT[i] );
2156              od;
2157              for j in [1..Length(pi_1)] do
2158                  Add( eqns[i], pi_1[j]*CT[i] );
2159              od;
2160          od;
2161          sol:= NullspaceMat( eqns );
2162          hZ:= List( sol, u -> u*h );
2163
2164          # Now we compute |Psi_0| and |Psi_1|...
2165
2166          psi0:= [ ];
2167          for a in rr[1].pv0 do
2168              if h0*a = 0*a and ForAll( hZ, u -> u*a = 0*a ) then
2169                 Add( psi0, a );
2170              fi;
2171          od;
2172
2173          psi1:= [ ];
2174          for a in rr[2].rv1 do
2175              if h0*a = a and ForAll( hZ, u -> u*a = 0*a ) then
2176                 Add( psi1, a );
2177              fi;
2178          od;
2179
2180          if Length(pi_0)+Length(pi_1) + 2*Length(psi0) = Length(psi1) then
2181
2182             if not 2*h0 in good then
2183                Add( good, 2*h0 );
2184             fi;
2185
2186          fi;
2187       fi;
2188   od;
2189
2190info:= "Obtained ";
2191Append( info, String( Length(good) ) );
2192Append( info, " Cartan elements, weeding out equivalent copies...");
2193Info(InfoSLA,2,info);
2194
2195# NEXT can be obtained from Kac diagram!!
2196
2197   x:= ChevalleyBasis(L)[1];
2198   y:= ChevalleyBasis(L)[2];
2199   es:= [ ];
2200   fs:= [ ];
2201   g0:= Subspace( L, Concatenation( Basis(CartanSubalgebra(L)), rr[1].pv0, rr[1].nv0 ) );
2202
2203   for i in [1..Length(CartanMatrix(R))] do
2204       if x[i] in g0 then
2205          Add( es, x[i] );
2206          Add( fs, y[i] );
2207       fi;
2208   od;
2209   hs:= List( [1..Length(es)], i -> es[i]*fs[i] );
2210
2211   valmat:= [ ];
2212   for i in [1..Length(hs)] do
2213       val:= [ ];
2214       for j in [1..Length(hs)] do
2215           Add( val, Coefficients( Basis( Subspace(L,[es[j]]), [es[j]] ),
2216                       hs[i]*es[j] )[1] );
2217       od;
2218       Add( valmat, val );
2219   od;
2220
2221
2222   chars:= [ ];
2223   for i in [1..Length(good)] do
2224
2225       u0:= good[i];
2226       v:= List( es, z -> Coefficients( Basis(Subspace(L,[z]),[z]), u0*z )[1] );
2227       done:= ForAll( v, z -> z >= 0 );
2228
2229       while not done do
2230           pos:= PositionProperty( v, z -> z < 0 );
2231           u0:= u0 - v[pos]*hs[pos];
2232           v:= v - v[pos]*valmat[pos];
2233           done:= ForAll( v, z -> z >= 0 );
2234       od;
2235
2236       if not u0 in chars then
2237          Add( chars, u0 );
2238       fi;
2239   od;
2240
2241   gr1:= rr[2].rv1;
2242   gr2:= rr[3].rvm;
2243
2244     g1:= Basis( Subspace( L, gr1 ), gr1 );
2245     g2:= Basis( Subspace( L, gr2 ), gr2 );
2246
2247     # the matrices of hL[i] acting on g1
2248     h_mats1:= [ ];
2249     for h0 in h do
2250         mat:= [ ];
2251         for i in [1..Length(g1)] do
2252             Add( mat, Coefficients( g1, h0*g1[i] ) );
2253         od;
2254         Add( h_mats1, mat );
2255     od;
2256
2257     # those of wrt g2...
2258     h_mats2:= [ ];
2259     for h0 in h do
2260         mat:= [ ];
2261         for i in [1..Length(g1)] do
2262             Add( mat, Coefficients( g2, h0*g2[i] ) );
2263         od;
2264         Add( h_mats2, mat );
2265     od;
2266
2267     sl2s:= [ ];
2268     id1:= IdentityMat( Length(g1) );
2269     id2:= IdentityMat( Length(g2) );
2270     #Omega:= [1..Dimension(L)];
2271     Omega:= [-1,0,1,1];
2272     for h0 in chars do
2273
2274         ch:= Coefficients( Basis( CartanSubalgebra(L) ), h0 );
2275         mat:= ch*h_mats1;
2276         mat:= mat - 2*id1;
2277         V:= NullspaceMat( mat );
2278         e:= List( V, v -> v*gr1 );
2279
2280         mat:= ch*h_mats2;
2281         mat:= mat + 2*id2;
2282         V:= NullspaceMat( mat );
2283         ff:= List( V, v -> v*gr2 );
2284
2285         found:= false;
2286         while not found do
2287
2288             co:= List( e, x -> Random(Omega) );
2289             x:= co*e;
2290             sp:= Subspace( L, List( ff, y -> x*y) );
2291
2292             if Dimension(sp) = Length(e) and h0 in sp then
2293
2294                # look for a nice one...
2295                for i in [1..Length(co)] do
2296                    k:= 0;
2297                    found:= false;
2298                    while not found do
2299                        co[i]:= k;
2300                        x:= co*e;
2301                        sp:= Subspace( L, List( ff, y -> x*y) );
2302
2303                        if Dimension(sp) = Length(e) and h0 in sp then
2304                           found:= true;
2305                        else
2306                           k:= k+1;
2307                        fi;
2308                    od;
2309                od;
2310
2311                mat:= List( ff, u -> Coefficients( Basis(sp), x*u ) );
2312                sol:= SolutionMat( mat, Coefficients( Basis(sp), h0 ) );
2313
2314                Add( sl2s, [sol*ff,h0,x] );
2315
2316                found:= true;
2317
2318
2319             fi;
2320
2321         od;
2322
2323     od;
2324
2325   return sl2s;
2326
2327end;
2328
2329###############################################################################################
2330#
2331#  method based on Weyl group action...
2332#
2333
2334SLAfcts.nil_orbits_weyl:= function( L, grading )
2335
2336     # grading is a list with the degree of each simple root..., required to be
2337     # non-negative.
2338
2339     local R, posR, posRv, negRv, g0, g1, gm, R1, D0, rank, inds0, v, i, perm,
2340           wrep, rts, w, N, p, D, P0, P1, j, es, fs, hs, valmat, val, chars,
2341           done, pos, u0, sg1, sgm, h_mats1, h_mats2, mat, sl2s, id1, id2, Omega,
2342           ch, V, e, ff, found, co, x, sp, k, c0, c1, s0, s1, pi_0, pi_1, t, pi,
2343           s, ct, eqns, rhs, C, CT, h, good, sol, h0, hZ, psi0, psi1, a, g00, eqn, info,
2344           orth, B, U, pU, CM, rr0, l0, l1, gr, deg;
2345
2346     R:= RootSystem(L);
2347     posR:= PositiveRootsNF(R);
2348     posRv:= PositiveRootVectors(R);
2349     negRv:= NegativeRootVectors(R);
2350
2351     g0:= ShallowCopy( BasisVectors( Basis( CartanSubalgebra(L) ) ) );
2352     g1:= [ ]; gm:= [ ];
2353     g00:= [ ];
2354
2355     R1:= [ ];
2356     D0:= [ ];
2357
2358     rank:= Length( CartanMatrix(R) );
2359     inds0:=[ ];
2360
2361     for i in [1..Length(posR)] do
2362         v:= posR[i]*grading;
2363         if v = 0 then
2364            Add( g0, posRv[i] );
2365            Add( g0, negRv[i] );
2366            Add( g00, posRv[i] );
2367            if i <= rank then Add( D0, posR[i] ); Add( inds0, i ); fi;
2368         elif v = 1 then
2369            Add( g1, posRv[i] );
2370            Add( gm, negRv[i] );
2371            Add( R1, posR[i] );
2372         fi;
2373     od;
2374
2375     perm:= SLAfcts.perms(R);
2376     wrep:= WeylTransversal( R, inds0 );
2377
2378info:= "Constructed a Weyl transversal of ";
2379Append( info, String(Length(wrep)) );
2380Append( info, " elements.");
2381
2382Info(InfoSLA,2,info);
2383
2384
2385     N:= Length( posR );
2386     rts:= Concatenation( posR, -posR );
2387
2388     ch:= ChevalleyBasis(L);
2389     h:= BasisVectors( Basis( CartanSubalgebra(L) ) );
2390     C:= CartanMatrix(R);
2391     CT:= TransposedMat( C );
2392     B:= BilinearFormMatNF( R );
2393     U:=[];
2394
2395     good:= [ ];
2396     for w in wrep do
2397
2398         p:= Product( perm{w} ); p:= p^-1;
2399         D:= rts{ List( [1..rank], i -> i^p ) };
2400         P0:= Intersection( D, D0 );
2401         P1:= Intersection( D, R1 );
2402
2403         if Length(P1) > 0 then
2404
2405            c0:= Combinations( [1..Length(P0)] );
2406
2407            for s0 in c0 do
2408
2409                    pi_0:= P0{s0};
2410                    pi_1:= P1;
2411
2412                    orth:= true;
2413                    for s in P0 do
2414                        if not s in pi_0 then
2415                           if ForAny( pi_0, x -> not IsZero( x*(B*s) ) ) or
2416                              ForAny( pi_1, x -> not IsZero( x*(B*s) ) ) then
2417                              orth:= false;
2418                              break;
2419                           fi;
2420                        fi;
2421                    od;
2422
2423                    if orth then
2424                       Sort( pi_0 ); Sort( pi_1 );
2425
2426                       pi:= Concatenation( pi_0, pi_1 );
2427                       CM:= List( pi, x -> List( pi, y -> 2*x*(B*y)/( y*(B*y) ) ) );
2428                       rr0:= SLAfcts.CartanMatrixToPositiveRoots( CM );
2429                       l0:= 0; l1:= 0;
2430                       gr:= Concatenation( List( pi_0, x -> 0 ), List( pi_1, x -> 1 ) );
2431                       for s in rr0 do
2432                           deg:= s*gr;
2433                           if deg=0 then
2434                              l0:= l0+1;
2435                           elif deg=1 then
2436                              l1:= l1+1;
2437                           fi;
2438                       od;
2439
2440                       if 2*l0+Length(pi) = l1 then
2441
2442                          t:= [ ];
2443                          for s in pi do
2444                              pos:= Position( posR, s );
2445                              Add( t, ch[1][pos]*ch[2][pos] );
2446                          od;
2447                          t:= BasisVectors( Basis( Subspace( L, t ) ) );
2448                          ct:= List( t, x -> Coefficients( Basis(CartanSubalgebra(L)), x ) );
2449
2450                          # i.e. t is a Cartan subalgebra of s
2451                          # find h0 in t such that a(h0)=1 for all a in pi_1, a(h0)=0
2452                          # for all a in pi_0
2453
2454                          eqns:=[ ];
2455                          rhs:= [ ];
2456                          for j in [1..Length(pi_0)] do
2457                              eqn:= [ ];
2458                              for i in [1..Length(t)] do
2459                                  eqn[i]:= pi_0[j]*( C*ct[i] );
2460                              od;
2461                              Add( eqns, eqn ); Add( rhs, 0 );
2462                          od;
2463                          for j in [1..Length(pi_1)] do
2464                              eqn:= [ ];
2465                              for i in [1..Length(t)] do
2466                                  eqn[i]:= pi_1[j]*( C*ct[i] );
2467                              od;
2468                              Add( eqns, eqn ); Add( rhs, 1 );
2469                          od;
2470
2471                          sol:= SolutionMat( TransposedMat(eqns), rhs );
2472                          h0:= sol*t;
2473
2474                          if not 2*h0 in good then
2475
2476                             # Find a basis of the subspace of h consisting of u with
2477                             # a(u) = 0, for a in pi = pi_0 \cup pi_1.
2478
2479                             eqns:= [ ];
2480                             for i in [1..Length(h)] do
2481                                 eqns[i]:= [ ];
2482                                 for j in [1..Length(pi_0)] do
2483                                     Add( eqns[i], pi_0[j]*CT[i] );
2484                                 od;
2485                                 for j in [1..Length(pi_1)] do
2486                                     Add( eqns[i], pi_1[j]*CT[i] );
2487                                 od;
2488                             od;
2489                             sol:= NullspaceMat( eqns );
2490                             hZ:= List( sol, u -> u*h );
2491
2492                             # Now we compute |Psi_0| and |Psi_1|...
2493
2494                             psi0:= [ ];
2495                             for a in g00 do
2496                                 if h0*a = 0*a and ForAll( hZ, u -> u*a = 0*a ) then
2497                                    Add( psi0, a );
2498                                 fi;
2499                             od;
2500
2501                             psi1:= [ ];
2502                             for a in g1 do
2503                                 if h0*a = a and ForAll( hZ, u -> u*a = 0*a ) then
2504                                    Add( psi1, a );
2505                                 fi;
2506                             od;
2507
2508                             if Length(pi_0)+Length(pi_1) + 2*Length(psi0) = Length(psi1) then
2509                                Add( good, 2*h0 );
2510                             fi;
2511                          fi;
2512                       fi;
2513                    fi;
2514                #od;
2515            od;
2516
2517         fi;
2518     od;
2519
2520info:= "Obtained ";
2521Append( info, String( Length(good) ) );
2522Append( info, " Cartan elements, weeding out equivalent copies...");
2523Info(InfoSLA,2,info);
2524
2525   es:= ChevalleyBasis(L)[1]{inds0};
2526   fs:= ChevalleyBasis(L)[2]{inds0};
2527   hs:= List( [1..Length(es)], i -> es[i]*fs[i] );
2528
2529   valmat:= [ ];
2530   for i in [1..Length(hs)] do
2531       val:= [ ];
2532       for j in [1..Length(hs)] do
2533           Add( val, Coefficients( Basis( Subspace(L,[es[j]]), [es[j]] ),
2534                       hs[i]*es[j] )[1] );
2535       od;
2536       Add( valmat, val );
2537   od;
2538
2539
2540   chars:= [ ];
2541   for i in [1..Length(good)] do
2542
2543       u0:= good[i];
2544       v:= List( es, z -> Coefficients( Basis(Subspace(L,[z]),[z]), u0*z )[1] );
2545       done:= ForAll( v, z -> z >= 0 );
2546       while not done do
2547           pos:= PositionProperty( v, z -> z < 0 );
2548           u0:= u0 - v[pos]*hs[pos];
2549           v:= v - v[pos]*valmat[pos];
2550           done:= ForAll( v, z -> z >= 0 );
2551       od;
2552
2553       if not u0 in chars then
2554          Add( chars, u0 );
2555       fi;
2556   od;
2557
2558     sg1:= Basis( Subspace( L, g1 ), g1 );
2559     sgm:= Basis( Subspace( L, gm ), gm );
2560
2561     # the matrices of hL[i] acting on g1
2562     h_mats1:= [ ];
2563     for h0 in h do
2564         mat:= [ ];
2565         for i in [1..Length(g1)] do
2566             Add( mat, Coefficients( sg1, h0*g1[i] ) );
2567         od;
2568         Add( h_mats1, mat );
2569     od;
2570
2571     # those of wrt gm...
2572     h_mats2:= [ ];
2573     for h0 in h do
2574         mat:= [ ];
2575         for i in [1..Length(gm)] do
2576             Add( mat, Coefficients( sgm, h0*gm[i] ) );
2577         od;
2578         Add( h_mats2, mat );
2579     od;
2580
2581     sl2s:= [ ];
2582     id1:= IdentityMat( Length(g1) );
2583     id2:= IdentityMat( Length(gm) );
2584     #Omega:= [1..Dimension(L)];
2585     Omega:= [-1,0,1,1];
2586     for h0 in chars do
2587
2588         ch:= Coefficients( Basis( CartanSubalgebra(L) ), h0 );
2589         mat:= ch*h_mats1;
2590         mat:= mat - 2*id1;
2591         V:= NullspaceMat( mat );
2592         e:= List( V, v -> v*g1 );
2593
2594         mat:= ch*h_mats2;
2595         mat:= mat + 2*id2;
2596         V:= NullspaceMat( mat );
2597         ff:= List( V, v -> v*gm );
2598
2599         found:= false;
2600
2601         while not found do
2602
2603             co:= List( e, x -> Random(Omega) );
2604             x:= co*e;
2605             sp:= Subspace( L, List( ff, y -> x*y) );
2606
2607             if Dimension(sp) = Length(e) and h0 in sp then
2608
2609                # look for a nice one...
2610                for i in [1..Length(co)] do
2611                    k:= 0;
2612                    found:= false;
2613                    while not found do
2614                        co[i]:= k;
2615                        x:= co*e;
2616                        sp:= Subspace( L, List( ff, y -> x*y) );
2617
2618                        if Dimension(sp) = Length(e) and h0 in sp then
2619                           found:= true;
2620                        else
2621                           k:= k+1;
2622                        fi;
2623                    od;
2624                od;
2625
2626                mat:= List( ff, u -> Coefficients( Basis(sp), x*u ) );
2627                sol:= SolutionMat( mat, Coefficients( Basis(sp), h0 ) );
2628
2629                Add( sl2s, [sol*ff,h0,x] );
2630
2631                found:= true;
2632
2633
2634             fi;
2635
2636         od;
2637
2638     od;
2639
2640   return sl2s;
2641
2642
2643end;
2644
2645
2646InstallOtherMethod( NilpotentOrbitsOfThetaRepresentation,
2647"for a Lie algebra and a grading diagram", true, [ IsLieAlgebra, IsList ], 0,
2648function( L, d )
2649
2650   local meth, rank, C, inds, i, w, tr, r;
2651
2652   meth:= ValueOption( "method" );
2653
2654   rank:= Length( d );
2655   C:= CartanMatrix( RootSystem(L) );
2656   if meth = fail then
2657      inds:= [ ];
2658      for i in [1..Length(d)] do
2659          if d[i] = 0 then Add( inds, i ); fi;
2660      od;
2661      w:= SizeOfWeylGroup( CartanType( C{inds}{inds} ).types );
2662      tr:= SizeOfWeylGroup( RootSystem(L) )/w;
2663      if tr > 8000 then
2664         meth:= "Carrier";
2665      else
2666         meth:= "WeylOrbit";
2667      fi;
2668   fi;
2669
2670   if meth = "WeylOrbit" then
2671      Info(InfoSLA,2,"Selected Weyl orbit method.");
2672      r:= SLAfcts.nil_orbits_weyl( L, d );
2673   else
2674      Info(InfoSLA,2,"Selected carrier algebra method.");
2675      r:= SLAfcts.zgrad_orbits_carrier( L, d );
2676   fi;
2677
2678   return r;
2679
2680end );
2681
2682#### functions for computing the Hasse diagram of the nil orbs of a theta group ##########
2683
2684# First: functions for dealing with the type of Weyl orbits that we need.
2685#
2686SLAfcts.NextIterator_WeylOrbitH:= function( it )
2687    local   output,  mu,  rank,  bound,  foundsucc,
2688            pos,  i,  nu,  a, B, v, lam, sims, stack, simples, stabs;
2689
2690    if it!.isDone then Error("the iterator is exhausted"); fi;
2691
2692    output:= it!.currentWeight;
2693
2694    B:= it!.B;
2695    v:= it!.v;
2696    lam:= it!.otherWt;
2697    stabs:= it!.stabinds;
2698    sims:= it!.simples;
2699
2700    #calculate the successor of curweight
2701
2702    mu:= ShallowCopy(it!.currentWeight);
2703    rank:= Length( mu );
2704    stack:= it!.stack;
2705    bound:= 1;
2706    foundsucc:= false;
2707    while not foundsucc do
2708
2709        pos:= fail;
2710        for i in [bound..rank] do
2711            if mu[i]>0 then
2712               nu:= ShallowCopy(mu);
2713               nu:= nu -nu[i]*sims[i];
2714               if ForAll( nu{[i+1..rank]}, x -> x >= 0 ) then
2715                  if Concatenation( nu, it!.c0 )*(B*Concatenation( lam, it!.c1 )) >= v then
2716                     pos:= i; break;
2717                  fi;
2718               fi;
2719            fi;
2720        od;
2721
2722        if pos <> fail then
2723            Add( stack, [ mu, pos ] );
2724            if not pos in stabs and ForAll( stabs, x -> nu[x] >= 0 ) then
2725               foundsucc:= true;
2726            else
2727               mu:= nu;
2728               bound:= 1;
2729            fi;
2730        else
2731
2732            if mu = it!.root then
2733
2734                # we cannot find a sucessor of the root: we are done
2735
2736                it!.isDone:= true;
2737                nu:= [];
2738                foundsucc:= true;
2739            else
2740                a:= stack[Length(stack)];
2741                mu:= a[1]; bound:= a[2]+1;
2742                RemoveElmList( stack, Length(stack) );
2743            fi;
2744
2745        fi;
2746
2747    od;
2748
2749    it!.stack:= stack;
2750    it!.currentWeight:= nu;
2751
2752    return output*it!.dualBas;
2753
2754end;
2755
2756SLAfcts.ClosureWeylOrbit:= function( L, K, H, basH, C, dualBas, B, h0, h1, c1, c0 )
2757
2758    # K: semsim Lie algebra,
2759    # H: CSA of K,
2760    # basH: basis of H, consisting of alpha_i^\vee
2761    # C: the Cartan matrix wrt the alpha_i
2762    # h0: element of H, of which we compute the orbit; assumed to be dominant.
2763    # h1: the other element, used for the Killing value...
2764
2765    local   mu,  nu, rank, Ci, i, j, c, v, bas, pos, sims, nu1;
2766
2767    rank:= Length(C);
2768
2769    bas:= Basis( H, dualBas );
2770
2771    mu:= Coefficients( bas, h0 );
2772    nu:= Coefficients(bas,h1);
2773
2774    nu1:= Concatenation( nu, c1 );
2775    v:= nu1*B*nu1;
2776
2777    sims:= List( basH, x -> Coefficients( bas, x ) );
2778
2779    return IteratorByFunctions( rec(
2780               IsDoneIterator := IsDoneIterator_WeylOrbit,
2781               NextIterator   := SLAfcts.NextIterator_WeylOrbitH,
2782
2783               ShallowCopy:= function( iter )
2784                      return rec( root:= ShallowCopy( iter!.root ),
2785                        currentWeight:= ShallowCopy( iter!.currentWeight ),
2786                        stack:= ShallowCopy( iter!.stack ),
2787                        otherWt:= ShallowCopy( iter!.otherWt ),
2788                        stabinds:= ShallowCopy( iter!.stabinds ),
2789                        B:= ShallowCopy( iter!.B ),
2790                        v:= iter!.v,
2791                        c0:= iter!.c0,
2792                        c1:= iter!.c1,
2793                        dualBas:= iter!.dualBas,
2794                        simples:= ShallowCopy( iter!.simples ),
2795                        isDone:= iter!.isDone );
2796                     end,
2797                        root:= mu,
2798                        currentWeight:= mu,
2799                        stack:= [ ],
2800                        otherWt:= nu,
2801                        stabinds:= Filtered( [1..Length(nu)], i -> nu[i]=0 ),
2802                        B:= B,
2803                        v:= v,
2804                        c0:= c0,
2805                        c1:= c1,
2806                        dualBas:= dualBas,
2807                        simples:= sims,
2808                        isDone:= false ) );
2809end;
2810
2811
2812##########################################################################################
2813# Next: the main procedure, starting with some subfunctions...
2814#
2815SLAfcts.cartanorb_prop:= function( L, K0, H0, basH, C, dualBas, B, h, oh, c1, c0, func )
2816
2817
2818    # h in H; its orbit under the Weyl group...
2819
2820    local o, hh, count;
2821
2822    o:= SLAfcts.ClosureWeylOrbit( L, K0, H0, basH, C, dualBas, B, h, oh, c1, c0 );
2823
2824    while not IsDoneIterator(o) do
2825        hh:= NextIterator(o);
2826        if func(hh) then
2827           return true;
2828        fi;
2829    od;
2830
2831    return false;
2832
2833end;
2834
2835
2836SLAfcts.normaliser:= function( L, K, U )
2837
2838   # L: Lie algebra
2839   # K: subalgebra
2840   # U: basis of a subspace of L
2841   # return: subalgebra of K, of all x such that [x,U] subset U.
2842
2843   local n, sp, m, s, eqns, k, d, r, j, i, sol, eqn;
2844
2845   n:= Dimension(K);
2846   sp:= Subspace(L,U);
2847   U:= BasisVectors( Basis(sp) );
2848   m:= Dimension(sp);
2849   s:= Dimension(L);
2850
2851   eqns:= NullMat( n+m^2, s*m );
2852   for k in [1..m] do
2853       d:= Coefficients( Basis(L), U[k] );
2854       for r in [1..s] do
2855           if not IsZero(d[r]) then
2856              for j in [1..m] do
2857                 eqns[n+k+(j-1)*m][j+(r-1)*m]:= eqns[n+k+(j-1)*m][j+(r-1)*m]+d[r];
2858              od;
2859           fi;
2860       od;
2861   od;
2862
2863   for i in [1..n] do
2864       for j in [1..m] do
2865           d:= Coefficients( Basis(L), Basis(K)[i]*U[j] );
2866           for r in [1..s] do
2867               if not IsZero(d[r]) then
2868                  eqns[i][j+(r-1)*m]:= eqns[i][j+(r-1)*m]+d[r];
2869               fi;
2870           od;
2871       od;
2872   od;
2873
2874   sol:= NullspaceMatDestructive( eqns );
2875   return List( sol, x -> x{[1..n]}*Basis(K) );
2876
2877end;
2878
2879
2880SLAfcts.inc:= function( sl2, domh, L, K, GM, G1, H0, basH, C, dualBas, B, i0, j0, file, matlist )
2881
2882   # see whether orbit i0 is inclued in j0.
2883
2884   local hi, hj, gh, W, ww, V2, U, wd, k, v0, f_chk, R, hh, BH, KH, Ci, Um, q,
2885         K0, h0, b0, sp0, h_start, c0, kk, KL, v, t, mats, oh, sp, kval, adh, Uinds, ord,
2886         maxU, i, j, mats0, c1, inconvexhull, kappamat, kapinv, invnu, dist0, ip0;
2887
2888
2889inconvexhull:= function( B, S0, p0, dist, ip, eps0 )
2890                                            # S set of vecs in R^m (rat coords),
2891                                            # p a point in R^m, is p\in S?
2892                                            # dist: distance fct
2893
2894    local m, i, one, eps, dists, pos, v, pp, k, j, u, t, S, p;
2895
2896    S:= List( S0, x -> Coefficients( B, x ) );
2897    p:= Coefficients( B, p0 );
2898    one:= 1.00000000000000000000000000;
2899    S:= List( S, u -> u*one );
2900    p:= p*one;
2901
2902    eps:= one*eps0;
2903
2904    dists:= List( S, s -> dist(s,p) );
2905    pos:= Position( dists, Minimum( dists ) );
2906    v:= S[pos];
2907    pp:= S[pos];
2908
2909    while true do
2910       if dist(p,pp) < eps*dist(p,v) then
2911          return [ pp, true ];
2912       else
2913          k:= 0;
2914          for j in [1..Length(S)] do
2915              if dist(pp,S[j]) >= dist(p,S[j]) then
2916                 k:= j; break;
2917              fi;
2918          od;
2919          if k > 0 then
2920             v:= S[k];
2921          else
2922             return [ pp, false ];
2923          fi;
2924       fi;
2925
2926       u:= pp-v;
2927       t:= -ip(u,p-pp)/ip(u,u);
2928       pp:= pp+t*(-u);
2929
2930    od;
2931end;
2932
2933   kappamat:= List( Basis(H0), h1 -> List( Basis(H0), h2 ->
2934        TraceMat( AdjointMatrix( Basis(L), h1 )*AdjointMatrix(Basis(L),h2)) ));
2935   kapinv:= kappamat^-1;
2936
2937   invnu:= function(x)  # x a root vec in g1, compute invnu of corr root.
2938
2939        local sp, b, u;
2940
2941        sp:= Basis( Subspace( L, [x] ), [x] );
2942        b:= List( Basis(H0), h1 -> Coefficients( sp, h1*x )[1] );
2943        return (kapinv*b)*Basis(H0);
2944   end;
2945
2946    dist0:= function(u,v) return (u-v)*kappamat*(u-v); end;
2947    ip0:= function(u,v) return u*kappamat*v; end;
2948
2949
2950
2951   ord:= function( i1, i2 )
2952
2953      if Length(i1) <> Length(i2) then
2954         return Length(i1) < Length(i2);
2955      else
2956         return i1 < i2;
2957      fi;
2958   end;
2959
2960   hi:= domh[i0];
2961   hj:= domh[j0];
2962
2963   kval:=  TraceMat( AdjointMatrix( Basis(L), hi )*AdjointMatrix(Basis(L),hj));
2964   kval:= kval/TraceMat(AdjointMatrix(Basis(L),hi)*AdjointMatrix(Basis(L),hi));
2965
2966   if kval < 1 then
2967      return false;
2968   fi;
2969
2970   ww:= [ ];
2971   for v in Basis(G1) do
2972       sp:= Basis( Subspace( L, [v] ), [v] );
2973       k:= Coefficients( sp, hi*v )[1];
2974       if k = 2 then Add( ww, v ); fi;
2975   od;
2976   V2:= ww;
2977
2978   K0:= LieDerivedSubalgebra(K);
2979   h0:= BasisVectors(Basis( LieCentre(K) ));
2980
2981   b0:= Concatenation( basH, h0 );
2982   sp0:= Basis( Subspace(L,b0), b0 );
2983
2984   R:= RootSystem(K0);
2985
2986   hh:= CanonicalGenerators(R)[3];
2987
2988   BH:= Basis( Subspace( L, hh ), hh );
2989   KH:= List( hh, x -> [] );
2990   adh:= List( hh, x -> AdjointMatrix( Basis(L), x ) );
2991   for i in [1..Length(hh)] do
2992       for j in [i..Length(hh)] do
2993           kval:= TraceMat( adh[i]*adh[j] );
2994           KH[i][j]:= kval;
2995           if i <> j then KH[j][i]:= kval; fi;
2996       od;
2997   od;
2998
2999   Um:= Filtered( Basis(GM), x -> hi*x = -2*x );
3000   q:= LieCentralizer( K, Subalgebra(K,[hi]) );
3001
3002   mats:= [ ];
3003   mats0:= [ ];
3004
3005   Uinds:= [ ];
3006   maxU:= [ ];
3007
3008   f_chk:= function( h, c0, t )
3009
3010      local gh, U, k, v0, u1, A, sol, P, M, sp, B, q0, inds, p,
3011            i, j, m, n, x, u, c, cf, R, a, s, v, d, r, V, t0, Uind, pos, cf0, found, Om,
3012            matrc, colinds, Bt, kval0, kval1, hhh, kv, h00, eps, bl;
3013
3014
3015      if Length(c0) > 0 then
3016         u1:= h+c0*h0;
3017      else
3018         u1:= h;
3019      fi;
3020
3021      U:= [ ];
3022      Uind:= [ ];
3023      for i in [1..Length(V2)] do
3024          v:= V2[i];
3025          sp:= Basis( Subspace( L, [v] ), [v] );
3026          k:= Coefficients( sp, u1*v )[1];
3027          if k >= 2 then Add( U, v ); Add( Uind, i ); fi;
3028      od;
3029
3030      if Length(U) = 0 then return false; fi;
3031
3032      pos:= PositionSorted( Uinds, Uind, ord );
3033      if IsBound(Uinds[pos]) and Uinds[pos] = Uind then
3034         return false;
3035      else
3036         Add( Uinds, Uind, pos );
3037      fi;
3038
3039      for i in [1..Length(maxU)] do
3040          if IsSubset( maxU[i], Uind ) then
3041             return false;
3042          fi;
3043      od;
3044      #
3045
3046      M:= SLAfcts.normaliser( L, K, U );
3047      Om:= [0..Dimension(L)];
3048      for k in [1..5] do
3049
3050          cf0:= List( U, x -> Random( Om ) );
3051          v0:= cf0*U;
3052          A:= List( Um, x -> Coefficients( Basis(K), v0*x ) );
3053          if Length(A) = 0 then
3054             sol:= fail;
3055          else
3056             sol:= SolutionMat( A, Coefficients( Basis(K), hi ) );
3057          fi;
3058          if sol <> fail then
3059             return true;
3060          else
3061             if Dimension( Subspace( L, v0*M ) ) = Length(U) then
3062                Add( maxU, Uind );
3063                return false;
3064             fi;
3065          fi;
3066
3067      od;
3068
3069      hhh:= List( U, x -> invnu(x) );
3070      kv:= TraceMat(AdjointMatrix(Basis(L),hi)*AdjointMatrix(Basis(L),hi));
3071      h00:= 2*hi/kv;
3072      eps:= 1/10;
3073      while eps > 1/10000000 do
3074         bl:= inconvexhull( Basis(H0), hhh, h00, dist0, ip0, eps )[2];
3075         if bl then
3076            eps:= eps/10;
3077         else
3078            return false;
3079         fi;
3080      od;
3081
3082      if Length(file) > 0 then
3083
3084         V:= Subspace( L, V2 );
3085
3086         m:= Length(U);
3087         if m = 0 then return false; fi;
3088         n:= Dimension(q);
3089         s:= Dimension(V);
3090         x:= BasisVectors( Basis(q) );
3091         v:= BasisVectors( Basis(V) );
3092
3093         c:= List( [1..n], r -> [ [ ] ] );
3094         for i in [1..n] do
3095             for j in [1..s] do
3096                 c[i][j]:= Coefficients( Basis( V ), x[i]*v[j] );
3097             od;
3098         od;
3099
3100         d:= List( U, r -> Coefficients( Basis(V), r ) );
3101
3102         P:= PolynomialRing( Rationals, m );
3103         a:= IndeterminatesOfPolynomialRing( P );
3104         A:= List( [1..s], r -> [ ] );
3105         for r in [1..s] do
3106             for i in [1..n] do
3107                 cf:= Zero(P);
3108                 for k in [1..m] do
3109                     for j in [1..s] do
3110                         cf:= cf + a[k]*d[k][j]*c[i][j][r];
3111                     od;
3112                 od;
3113                 A[r][i]:= cf;
3114             od;
3115         od;
3116
3117         if Dimension(q)-Length(A) > t then
3118            return false;
3119         fi;
3120
3121         if ForAny(A, IsZero ) then
3122            return false;
3123         fi;
3124
3125         matrc:= rec( inds:= [i0,j0], numindets:= Length(a), fullmat:= A );
3126
3127         q0:= BasisVectors( CanonicalBasis( Intersection(
3128                                     q, KappaPerp( L, Subalgebra(L,[hi]) ) ) ) );
3129         B:= [ ];
3130         # first we find a max set of lin indep rows, for a random v0\in U
3131         cf0:= List( U, x -> Random( [1..Dimension(L)^2] ) );
3132         v0:= cf0*U;
3133         for i in [1..Length(q0)] do
3134             Add( B, Coefficients( Basis(V), q0[i]*v0 ) );
3135         od;
3136         inds:= [ 1 ];
3137         sp:= MutableBasis( LeftActingDomain(L), [ ], List( Basis(V), x -> 0 ) );
3138         for i in [1..Length(B)] do
3139             if not IsContainedInSpan( sp, B[i] ) then
3140                CloseMutableBasis( sp, B[i] );
3141                Add( inds, i+1 );
3142             fi;
3143         od;
3144
3145         colinds:= [ ];
3146         Bt:= TransposedMat(B{[1..Length(B)]});
3147         sp:= MutableBasis( LeftActingDomain(L), [ ], List( [1..Length(q0)], x -> 0 ) );
3148         for i in [1..Length(Bt)] do
3149             if not IsContainedInSpan( sp, Bt[i] ) then
3150                CloseMutableBasis( sp, Bt[i] );
3151                Add( colinds, i );
3152             fi;
3153         od;
3154
3155         # now make the "real" matrix...
3156         x:= Concatenation( [hi], q0 );
3157         A:= List( inds, zz -> [ ] );
3158         c:= List( [1..n], r -> [ [ ] ] );
3159         for i in [1..n] do
3160             for j in [1..m] do
3161                 c[i][j]:= Coefficients( Basis( V ), x[i]*U[j] );
3162             od;
3163         od;
3164
3165         for i in [1..Length(inds)] do
3166             for j in [1..s] do
3167                 p:= Zero(P);
3168                 for k in [1..m] do
3169                     p:= p + a[k]*c[inds[i]][k][j];
3170                     A[i][j]:= p;
3171                 od;
3172             od;
3173         od;
3174
3175         matrc.redmat:= A;
3176         matrc.colinds:= colinds;
3177
3178         Add( matlist, matrc );
3179
3180      else
3181        Add( matlist, rec( inds:= [i0,j0] ) );
3182      fi;
3183
3184      return false;
3185
3186   end;
3187
3188   t:= Dimension( Intersection( LieCentralizer( K, Subalgebra(K,[sl2[i0][2]]) ),
3189                LieCentralizer(L,Subalgebra(L,[sl2[i0][3]])) ) );
3190
3191   Ci:= CartanMatrix(R)^-1;
3192
3193   h_start:= Coefficients( sp0, hj ){[1..Length(basH)]}*basH;
3194   c0:= Coefficients( sp0, hj ){[ Length(basH)+1..Length(b0) ]};
3195   c1:= Coefficients( sp0, hi ){[ Length(basH)+1..Length(b0) ]};
3196
3197   oh:= Coefficients( sp0, hi ){[1..Length(basH)]}*basH;
3198
3199   v:= SLAfcts.cartanorb_prop( L, K0, H0, basH, C, dualBas, B, h_start, oh, c1, c0,
3200                        h -> f_chk(h,c0,t) );
3201
3202   return v;
3203
3204end;
3205
3206
3207SLAfcts.is_included_in:= function( diag, j )
3208
3209    local list, done, list0, d, pos;
3210
3211    list:= [ j ];
3212    done:= false;
3213    while not done do
3214
3215        list0:= [ ];
3216        for d in diag do
3217            if d[1] in list then
3218               AddSet( list0, d[2] );
3219            fi;
3220        od;
3221        list0:= Set( Concatenation( list, list0 ) );
3222        if Length(list0) = Length(list) then
3223           done:= true;
3224        fi;
3225        list:= list0;
3226    od;
3227
3228    pos:= Position( list, j );
3229    RemoveElmList( list, pos );
3230
3231    return list;
3232
3233end;
3234
3235SLAfcts.minimize_diag:= function( N, diag )
3236
3237    # if [i,j] and [j,k] and [i,k] then get rid of [i,k]...
3238
3239    local edges, s, i, j, k, path2, p1, p2;
3240
3241    Sort(diag);
3242
3243    edges:= [ ];
3244    for s in diag do
3245
3246        i:= s[1]; j:= s[2];
3247        path2:= false;
3248        for k in [1..N] do
3249            if k <> i and k <> j then
3250               p1:= PositionSorted( diag, [i,k] );
3251               p2:= PositionSorted( diag, [k,j] );
3252               if IsBound(diag[p1]) and IsBound(diag[p2]) and
3253                  diag[p1] = [i,k] and diag[p2] = [k,j] then
3254               #if [i,k] in diag and [k,j] in diag then
3255                  path2:= true;
3256               fi;
3257            fi;
3258            if path2 then break; fi;
3259        od;
3260
3261        if not path2 then
3262           Add( edges, s );
3263        fi;
3264    od;
3265
3266    return edges;
3267
3268end;
3269
3270
3271
3272SLAfcts.hasse_diag:= function( L, grad, sl2 )
3273
3274   local K, GM, G1, dim, dim1, d1, d2, diag, i, j, k, incs, b, file, K0, R, posRv,
3275         posR, negRv, fundR, sums, inds, basH, H0, rank, C, B, posR_L, dualBas,
3276         Ci, c, g0, g1, gm, gsp, Cu, domh, sims, pos, bas, mu, h0, b0, c0, sp0,
3277         m, m0, r, numvar, matlist, info, set, magmaprog, dualBas0;
3278
3279   g0:= grad[1]; gm:= grad[ Length(grad) ];
3280   if Length( grad ) > 1 then
3281      g1:= grad[2];
3282   else
3283      g1:= grad[1];
3284   fi;
3285
3286   file:= ValueOption( "filenm" );
3287   if file = fail then file:= ""; fi;
3288
3289   K:= Subalgebra( L, g0 );
3290   GM:= Subspace( L, gm );
3291   G1:= Subspace( L, g1 );
3292
3293   K0:= LieDerivedSubalgebra(K);
3294
3295   R:= RootSystem(K0);
3296   posR:= PositiveRootsNF(R);
3297   fundR:= SimpleSystemNF(R);
3298   inds:= List( fundR, x -> Position( posR, x ) );
3299   posRv:= PositiveRootVectors(R);
3300   negRv:= NegativeRootVectors(R);
3301   basH:= List( inds, i -> posRv[i]*negRv[i] );
3302
3303   H0:= Subalgebra( K0, basH );
3304
3305   rank:= Length( basH );
3306   C:= CartanMatrix(R);
3307
3308    dualBas:= [ ];
3309    Ci:= C^-1;
3310    for i in [1..rank] do
3311        c:= 0*[1..rank];
3312        c[i]:= 1;
3313        c:= Ci*c;
3314        Add( dualBas, c*basH );
3315    od;
3316
3317    dualBas0:= Concatenation( dualBas, BasisVectors(Basis( LieCentre(K) ) ) );
3318    B:= List( [1..Length(dualBas0)], x -> [] );
3319    for i in [1..Length(dualBas0)] do
3320        for j in [i..Length(dualBas0)] do
3321            B[i][j]:= TraceMat( AdjointMatrix( Basis(L), dualBas0[i] )*
3322                                AdjointMatrix( Basis(L), dualBas0[j] ) );
3323            B[j][i]:= B[i][j];
3324        od;
3325    od;
3326
3327   dim:= function( u ) return Dimension( Subspace( L, Basis(K)*u[3] ) ); end;
3328
3329   d1:= List( sl2, dim );
3330
3331   SortParallel( d1, sl2 );
3332   d1:= Reversed(d1);
3333   sl2:= Reversed(sl2);
3334
3335   # now map all h-s to dominant Weyl chamber.
3336   h0:= BasisVectors(Basis( LieCentre(K) ));
3337   b0:= Concatenation( basH, h0 );
3338   sp0:= Basis( Subspace(L,b0), b0 );
3339
3340   domh:= [ ];
3341   bas:= Basis( H0, dualBas );
3342   sims:= List( basH, x -> Coefficients( bas, x ) );
3343   for i in [1..Length(sl2)] do
3344
3345       h0:= Coefficients( sp0, sl2[i][2] ){[1..Length(basH)]}*basH;
3346       c0:= Coefficients( sp0, sl2[i][2] ){[ Length(basH)+1..Length(b0) ]};
3347
3348       mu:= Coefficients( bas, h0 );
3349       pos:= PositionProperty( mu, x -> x < 0 );
3350       while pos <> fail do
3351          mu:= mu -mu[pos]*sims[pos];
3352          pos:= PositionProperty( mu, x -> x < 0 );
3353       od;
3354       if Length(c0) = 0 then
3355          Add( domh, mu*bas );
3356       else
3357          Add( domh, mu*bas + c0*b0{[Length(basH)+1..Length(b0) ]} );
3358       fi;
3359   od;
3360
3361   gsp:= List( grad, x -> Subspace( L, x ) );
3362   d2:= [ ];
3363   for i in [1..Length(sl2)] do
3364       Cu:= LieCentralizer( L, Subalgebra( L, [sl2[i][3]] ) );
3365       Add( d2, List( gsp, x -> Dimension( Intersection( x, Cu ) ) ) );
3366   od;
3367
3368   diag:= [ ];
3369
3370   # [ i, j ] in diag means orbit i is included in the closure of
3371   # orbit j.
3372
3373   matlist:= [ ];
3374
3375   for i in [2..Length(sl2)] do
3376       for j in [i-1,i-2..1] do
3377
3378           if (not [i,j] in diag) and d1[i] < d1[j] and ForAll( [1..Length(grad)], k ->
3379              d2[i][k] >= d2[j][k] ) then
3380
3381              b:= SLAfcts.inc( sl2, domh, L, K, GM, G1, H0, basH, C, dualBas, B, i, j, file, matlist );
3382              if b then
3383                 Add( diag, [i,j] );
3384                 incs:= SLAfcts.is_included_in( diag, j );
3385                 for k in incs do Add( diag, [i,k] ); od;
3386              fi;
3387           fi;
3388       od;
3389   od;
3390
3391   matlist:= Filtered( matlist, x -> not x.inds in diag );
3392
3393   if Length(matlist) = 0 then
3394      Info( InfoSLA,2,"All (non-) inclusions proved!");
3395   else
3396      info:= "For the following pairs of orbits the first could be included in the\nclosure of the second\
3397 (but this is unlikely):\n";
3398      set:= Set( List( matlist, r -> r.inds ) );
3399      for i in  [1..Length(set)] do
3400          Append( info, String(set[i]) );
3401          if i < Length( set ) then
3402             Append( info, ", " );
3403          fi;
3404      od;
3405      Info( InfoSLA,2, info );
3406   fi;
3407
3408   if Length(file) > 0 and Length(matlist) > 0 then
3409
3410magmaprog:=
3411"minors:= function( m )\n\n\
3412   \/\/ m and rxs matrix with s >= r, compute all rxr minors,\n\
3413   \/\/ return false if a nonzero found; otherwise return true.\n\n\
3414   r:= NumberOfRows(m);\n\
3415   s:= NumberOfColumns(m);\n\
3416   if r eq s then return Determinant(m) eq 0; end if;\n\
3417   rows:= [1..r];\n\
3418   done:= false;\n\
3419   len:= s-r;\n\
3420   exc:= [1..len];\n\
3421   while not done do\n\
3422        cols:= [ ];\n\
3423        for i in [1..s] do\n\
3424            if not i in exc then Append( ~cols, i ); end if;\n\
3425        end for;\n\
3426        d:= Minor( m, rows, cols ); \n\
3427        if not IsZero(d) then return false; end if;\n\
3428        pos:= len;\n\
3429        found:= (exc[pos] lt s); \n\
3430        while not found do \n\
3431            pos:= pos-1; \n\
3432            if pos eq 0 then // done... \n\
3433               done:= true;\n\
3434               found:= true;\n\
3435            else\n\
3436               found:= (exc[pos] lt s-(len-pos));\n\
3437            end if;\n\
3438        end while;\n\
3439        if not done then \n\
3440           l:= exc[pos]+1;\n\
3441           for i in [pos..len] do \n\
3442               exc[i]:= l+i-pos; \n\
3443           end for; \n\
3444        end if; \n\
3445   end while; \n\
3446   return true;\n\
3447end function;\n\n\n\n\
3448minors0:= function( m, cols )\n\n\
3449   \/\/ m and rxs matrix with s > r \n\
3450   \/\/ and cols are indices of lin indep columns.\n\
3451   \/\/ check whether first row is in span of last r-1 rows...\n\n\
3452   r:= NumberOfRows(m);\n\
3453   s:= NumberOfColumns(m);\n\
3454   rows:= [1..r]; \n\
3455   for i in [1..s] do \n\
3456       if not i in cols then \n\
3457          columns:= cols; \n\
3458          Append( ~columns, i ); \n\
3459          Sort( ~columns ); \n\
3460          d:= Minor( m, rows, columns ); \n\
3461          if not IsZero(d) then return false; end if; \n\
3462       end if; \n\
3463   end for; \n\
3464   return true; \n\
3465end function; \n\n";
3466
3467      AppendTo( file, magmaprog );
3468
3469      numvar:= 0;
3470      for r in matlist do
3471          if r.numindets > numvar then numvar:= r.numindets; fi;
3472      od;
3473      AppendTo( file, "F<" );
3474      for i in [1..numvar] do
3475          AppendTo( file, "x_" );
3476          AppendTo( file, i );
3477          if i < numvar then
3478             AppendTo(file,",");
3479          fi;
3480      od;
3481      AppendTo( file, ">:= RationalFunctionField( Rationals(), ");
3482      AppendTo( file, numvar );
3483      AppendTo( file, ");\n\n");
3484
3485      for r in matlist do
3486          AppendTo( file, "print \"inclusion: orbit \", ",r.inds[1], ", \" in orbit \",",r.inds[2],";\n");
3487          m:= r.fullmat; m0:= r.redmat;
3488          AppendTo( file, "m:= \n ",m,";\n\n", "m:= Matrix(m);\n\n",
3489                          "m0:= \n ",m0,";\n\n", "m0:= Matrix(m0);\n\n" );
3490          AppendTo( file, "cols:= ",r.colinds,";\n\n" );
3491          if Length(m[1])-Length(m) < Length(m0[1])-Length(m0) then
3492             AppendTo( file, "minors(m);\n\n" );
3493          else
3494             AppendTo( file,
3495                    "if not minors0(m0,cols) then minors(m); else true; end if;\n\n");
3496          fi;
3497      od;
3498   fi;
3499
3500   return rec( diag:= diag, sl2:= sl2 );
3501
3502end;
3503
3504InstallMethod( ClosureDiagram,
3505"for Lie algebra, list or automorphism, list of sl2 triples", true, [ IsLieAlgebra, IsObject, IsList ], 0,
3506function( L, obj, sl2 )
3507
3508   local d, g0, g1, gm, R, posR, posRv, negRv, i, v, r, diag, f, g;
3509
3510   if IsList( obj ) then
3511      d:= obj;
3512      g0:= ShallowCopy( ChevalleyBasis(L)[3] );
3513      g1:= [ ]; gm:= [ ];
3514      R:= RootSystem(L);
3515      posR:= PositiveRootsNF(R);
3516      posRv:= PositiveRootVectors(R);
3517      negRv:= NegativeRootVectors(R);
3518
3519      for i in [1..Length(posR)] do
3520          v:= posR[i]*d;
3521          if v = 0 then
3522              Add( g0, posRv[i] );
3523              Add( g0, negRv[i] );
3524          elif v = 1 then
3525            Add( g1, posRv[i] );
3526            Add( gm, negRv[i] );
3527         fi;
3528      od;
3529      if Length(g1) > 0 then
3530         r:= SLAfcts.hasse_diag( L, [g0,g1,gm], sl2 );
3531      else
3532         r:= SLAfcts.hasse_diag( L, [g0], sl2 );
3533      fi;
3534   elif IsMapping( obj ) then
3535      f:= obj;
3536      g:= Grading(f);
3537      r:= SLAfcts.hasse_diag( L, g, sl2 );
3538   else
3539      Error( "the second argument has to be an automorphism or a list giving a Z-grading");
3540   fi;
3541
3542   diag:= SLAfcts.minimize_diag( Length(sl2), r.diag );
3543   return rec( diag:= diag, sl2:= r.sl2 );
3544
3545end );
3546
3547SLAfcts.carrZ:= function( L, diag, e )
3548
3549   # L Z-graded by diag, e in L(1); get its carrier algebra.
3550
3551   local R, posR, pRv, nRv, g0, i, v, K, h, sp, gp, gn, lams, b, C, c, eigensp, gpos, gneg,
3552         gzero, h1, good;
3553
3554   R:= RootSystem(L);
3555   posR:= PositiveRootsNF(R);
3556   pRv:= PositiveRootVectors(R);
3557   nRv:= NegativeRootVectors(R);
3558
3559   gpos:= [ ];
3560   gneg:= [ ];
3561
3562   gzero:= ShallowCopy( CanonicalGenerators(R)[3] );
3563   for i in [1..Length(posR)] do
3564       v:= posR[i]*diag;
3565       if v = 0 then
3566          Add( gzero, pRv[i] );
3567          Add( gzero, nRv[i] );
3568       elif v > 0 then
3569          if not IsBound( gpos[v] ) then
3570             gpos[v]:= [ pRv[i] ];
3571             gneg[v]:= [ nRv[i] ];
3572          else
3573             Add( gpos[v], pRv[i] );
3574             Add( gneg[v], nRv[i] );
3575          fi;
3576       fi;
3577   od;
3578   K:= Subalgebra( L, gzero );
3579
3580   h1:= BasisVectors( CanonicalBasis( Intersection( LieNormalizer( L, Subalgebra(L,[e]) ),
3581             CartanSubalgebra(L) ) ) );
3582   h:= BasisVectors(CanonicalBasis( ( CartanSubalgebra(LieNormalizer(L,Subalgebra(L,[e]))))));
3583
3584
3585   if Length(h) = Length(h1) then
3586      good:= true;
3587   else
3588      good:= false;
3589   fi;
3590
3591   lams:= [ ];
3592   sp:= Basis( Subspace( L, [e] ), [e] );
3593   for i in [1..Length(h)] do
3594       Add( lams, Coefficients( sp, h[i]*e )[1] );
3595   od;
3596
3597   gp:= [ ]; gn:= [ ];
3598   if good then
3599
3600      g0:= ShallowCopy( CanonicalGenerators(R)[3] );
3601      b:= List( h, x -> Coefficients( CanonicalBasis(CartanSubalgebra(L)), x ) );
3602      C:= CartanMatrix( R );
3603
3604      for i in [1..Length(posR)] do
3605          v:= posR[i]*diag;
3606          c:= List( b, x ->  posR[i]*C*x );
3607
3608          if v = 0 then
3609             if c = 0*c then
3610                Add( g0, pRv[i] );
3611                Add( g0, nRv[i] );
3612             fi;
3613          else
3614             if c = v*lams then
3615                if not IsBound( gp[v] ) then
3616                   gp[v]:= [];
3617                   gn[v]:= [];
3618                fi;
3619                Add( gp[v], pRv[i] );
3620                Add( gn[v], nRv[i] );
3621             fi;
3622          fi;
3623      od;
3624
3625   else
3626
3627      eigensp:= function( uu, t )
3628
3629         local m, s, sp, eqns, i, j, k, c, sol;
3630
3631         m:= Length(h);
3632         s:= Length(uu);
3633         sp:= Basis( Subspace( L, uu ), uu );
3634         eqns:= NullMat( s, s*m );
3635         for j in [1..m] do
3636             for i in [1..s] do
3637                 c:= Coefficients( sp, h[j]*uu[i] );
3638                 for k in [1..s] do
3639                     eqns[i][(k-1)*m+j]:= c[k];
3640                 od;
3641             od;
3642         od;
3643         for k in [1..s] do
3644             for j in [1..m] do
3645                 eqns[k][(k-1)*m+j]:= eqns[k][(k-1)*m+j]-t*lams[j];
3646             od;
3647         od;
3648
3649         sol:= NullspaceMat( eqns );
3650         return List( sol, x -> x*uu );
3651      end;
3652
3653      g0:= eigensp( gzero, 0 );
3654      for i in [1..Length(gpos)] do
3655          if IsBound( gpos[i] ) then
3656             gp[i]:= eigensp( gpos[i], i );
3657             gn[i]:= eigensp( gneg[i], -i );
3658          fi;
3659      od;
3660
3661   fi;
3662
3663   K:= Subalgebra(L,Concatenation(g0,Flat(gp),Flat(gn)));
3664   K:= LieDerivedSubalgebra(K);
3665
3666   sp:= Intersection( Subspace(L,g0), K );
3667   g0:= BasisVectors( Basis( sp ) );
3668   for i in [1..Length(gp)] do
3669       if IsBound( gp[i] ) then
3670          sp:= Intersection( Subspace(L,gp[i]), K );
3671          gp[i]:= BasisVectors( Basis( sp ) );
3672       else
3673          gp[i]:= [ ];
3674       fi;
3675   od;
3676   for i in [1..Length(gn)] do
3677       if IsBound( gn[i] ) then
3678          sp:= Intersection( Subspace(L,gn[i]), K );
3679          gn[i]:= BasisVectors( Basis( sp ) );
3680       else
3681          gn[i]:= [ ];
3682       fi;
3683   od;
3684
3685   return rec( g0:= g0, posdeg:= gp, negdeg:= gn );
3686
3687end;
3688
3689
3690SLAfcts.carrZm:= function( L, f, e )
3691
3692   local h, lams, sp, i, gp, gn, eigensp, g0, g1, gm, m, gr, K, k, dim,t0;
3693
3694   gr:= Grading(f);
3695   sp:= Subspace( L, gr[1] );
3696   h:= BasisVectors(CanonicalBasis( Intersection( sp, CartanSubalgebra(LieNormalizer(L,
3697                            Subalgebra(L,[e]))))));
3698   lams:= [ ];
3699   sp:= Basis( Subspace( L, [e] ), [e] );
3700   for i in [1..Length(h)] do
3701       Add( lams, Coefficients( sp, h[i]*e )[1] );
3702   od;
3703
3704   gp:= [ ]; gn:= [ ];
3705
3706    eigensp:= function( uu, t )
3707
3708         local m, s, sp, eqns, i, j, k, c, sol;
3709
3710         m:= Length(h);
3711         s:= Length(uu);
3712         sp:= Basis( Subspace( L, uu ), uu );
3713         eqns:= NullMat( s, s*m );
3714         for j in [1..m] do
3715             for i in [1..s] do
3716                 c:= Coefficients( sp, h[j]*uu[i] );
3717                 for k in [1..s] do
3718                     eqns[i][(k-1)*m+j]:= c[k];
3719                 od;
3720             od;
3721         od;
3722         for k in [1..s] do
3723             for j in [1..m] do
3724                 eqns[k][(k-1)*m+j]:= eqns[k][(k-1)*m+j]-t*lams[j];
3725             od;
3726         od;
3727
3728         sol:= NullspaceMat( eqns );
3729         return List( sol, x -> x*uu );
3730      end;
3731
3732   m:= Length(gr);
3733   g0:= eigensp( gr[1], 0 );
3734   g1:= eigensp( gr[2], 1 );
3735   gm:= eigensp( gr[ m ], -1 );
3736
3737   K:= LieDerivedSubalgebra( Subalgebra( L, Concatenation( gm, g0, g1 ) ) );
3738
3739   g0:= BasisVectors( Basis( Intersection( Subspace( L, g0 ), K ) ) );
3740
3741   dim:= Length(g0);
3742   k:= 1;
3743   while dim < Dimension(K) do
3744      g1:= BasisVectors( Basis( Intersection( Subspace( L,
3745              eigensp( gr[ (k mod m) +1 ], k ) ), K ) ) );
3746      Add( gp, g1 );
3747      dim:= dim+Length(g1);
3748      gm:= BasisVectors( Basis( Intersection( Subspace( L,
3749              eigensp( gr[ (-k mod m) +1 ], -k ) ), K ) ) );
3750      Add( gn, gm );
3751      dim:= dim+Length(gm);
3752      k:= k+1;
3753   od;
3754
3755   return rec( g0:= g0, gp:= gp, gn:= gn );
3756
3757end;
3758
3759InstallMethod( CarrierAlgebra,
3760"for Lie algebra, list or automorphism, a nilpotent element", true,
3761[ IsLieAlgebra, IsObject, IsObject ], 0,
3762function( L, obj, e )
3763
3764   local d, g0, g1, gm, R, posR, posRv, negRv, i, v, r, diag, f, g;
3765
3766   if IsList( obj ) then
3767      d:= obj;
3768      r:= SLAfcts.carrZ( L, d, e );
3769   elif IsMapping( obj ) then
3770      f:= obj;
3771      r:= SLAfcts.carrZm( L, f, e );
3772   else
3773      Error( "the second argument has to be an automorphism or a list giving a Z-grading");
3774   fi;
3775
3776   return r;
3777
3778end );
3779
3780
3781SLAfcts.jord_dec:= function ( mat )
3782
3783    local  F, p, B, f, g, fac, ff, h, w;
3784
3785    F := DefaultFieldOfMatrix( mat );
3786    if F = fail  then
3787        return TRY_NEXT_METHOD;
3788    fi;
3789    p := Characteristic( F );
3790    f := CharacteristicPolynomial( F, F, mat );
3791    fac := Set( Factors( f ) );
3792    g:= Product( fac );
3793    if f = g  then
3794        return [ mat, 0 * mat ];
3795    fi;
3796    w := GcdRepresentation( g, Derivative( g ) )[2];
3797    w := w * g;
3798    B := ShallowCopy( mat );
3799    while Value( g, B ) <> 0 * B  do
3800        B := B - Value( w, B );
3801    od;
3802    return [ B, mat - B ];
3803end;
3804
3805SLAfcts.semsim_part:= function( L, x, b, B )
3806
3807  # x elm of semsim Lie alg, L, b: basis of Lie alg
3808  # B: basis of space spanned by ad L, such that i-th
3809  # basis elements is ad b_i.
3810
3811  local adx, cf;
3812
3813  adx:= AdjointMatrix( Basis(L), x );
3814  cf:= Coefficients( B, SLAfcts.jord_dec( adx )[1] );
3815  return cf*b;
3816
3817end;
3818
3819InstallMethod( CartanSubspace,
3820"for a finite order automorphism", true, [ IsGeneralMapping ], 0,
3821function( f )
3822
3823  local F, sp0, sp1, g, s, r, Omega, s1, b, ad, B, found, x, xs, s0, V,
3824        i, cf, cf1, x1, xs1, adx, found0, L, g0, g1;
3825
3826  L:= Source(f);
3827  g0:= Grading(f)[1];
3828  g1:= Grading(f)[2];
3829
3830  F:= LeftActingDomain(L);
3831
3832  sp0:= Subspace( L, g0 );
3833  sp1:= Subspace( L, g1 );
3834
3835  g:= L;
3836  s:= L;
3837  r:= Subalgebra( L, [] );
3838
3839  Omega:= [ 0, 1, 1, 1 ];
3840
3841  while true do
3842
3843     s1:= Intersection( s, sp1 );
3844     if Dimension( s1 ) = 0 then
3845        return Intersection( g, sp1 );
3846     fi;
3847
3848     b:= BasisVectors( Basis(s) );
3849     ad:= List( b, x -> AdjointMatrix( Basis(s), x ) );
3850     B:= Basis( VectorSpace( F, ad ), ad );
3851
3852     found:= false;
3853
3854     while not found do
3855
3856          cf:= List( Basis(s1), u -> Random( Omega) );
3857          x:= cf*BasisVectors( Basis(s1) );
3858          adx:= AdjointMatrix( Basis(s), x );
3859
3860          if not IsZero(adx^Dimension(s)) then
3861
3862             found0:= false;  # we find one with maybe better coefficients
3863             while not found0 do
3864               cf:= List( Basis(s1), u -> Random( [0,1] ) );
3865               x:= cf*BasisVectors( Basis(s1) );
3866               adx:= AdjointMatrix( Basis(s), x );
3867               found0:= not IsZero( adx^Dimension(s) );
3868             od;
3869
3870             for i in [1..Length(cf)] do
3871                 if not IsZero(cf[i]) then
3872                    cf1:= ShallowCopy(cf);
3873                    cf1[i]:= 0;
3874                    x1:= cf1*BasisVectors( Basis(s1) );
3875                    adx:= AdjointMatrix( Basis(s), x1 );
3876                    if not IsZero( adx^Dimension(s) ) then
3877                       x:= x1; cf:= cf1;
3878                    fi;
3879                 fi;
3880             od;
3881
3882             for i in [1..Length(cf)] do
3883                 if not IsZero(cf[i]) then
3884                    cf1:= ShallowCopy(cf);
3885                    cf1[i]:= 0;
3886                    x1:= cf1*BasisVectors( Basis(s1) );
3887                    adx:= AdjointMatrix( Basis(s), x1 );
3888                    if not IsZero( adx^Dimension(s) ) then
3889                       x:= x1; cf:= cf1;
3890                    fi;
3891                 fi;
3892             od;
3893
3894             xs:= SLAfcts.semsim_part( s, x, b, B );
3895
3896             g:= LieCentralizer( g, Subalgebra( g, [xs] ) );
3897             s:= LieDerivedSubalgebra(g);
3898             r:= LieCentre(g);
3899
3900             found:= true;
3901          else
3902
3903             s0:= Intersection( s, sp0 );
3904             V:= VectorSpace( F, List( Basis(s0), u -> u*x ) );
3905             if Dimension(V) = Dimension(s1) then
3906                return Intersection( r, sp1 );
3907             fi;
3908          fi;
3909
3910     od;
3911
3912  od;
3913
3914end );
3915