1#############################################################################
2##
3##  This file is part of GAP, a system for computational discrete algebra.
4##  This file's authors include Thomas Breuer, Ansgar Kaup.
5##
6##  Copyright of GAP belongs to its developers, whose names are too numerous
7##  to list here. Please refer to the COPYRIGHT file for details.
8##
9##  SPDX-License-Identifier: GPL-2.0-or-later
10##
11##  This file contains functions that mainly deal with lattices in the
12##  context of character tables.
13##
14
15
16#############################################################################
17##
18#F  LLL( <tbl>, <characters>[, <y>][, \"sort\"][, \"linearcomb\"] )
19##
20InstallGlobalFunction( LLL, function( arg )
21
22    local tbl,         # character table, first argument
23          characters,  # list of virtual characters, second argument
24          sorted,      # characters sorted by degree
25          L,           # lattice gen. by the virtual characters
26          i,           # loop variable
27          lllrb,       # result of `LLLReducedBasis'
28          lll,         # result
29          v,           # loop over the LLL reduced basis
30          perm,        # permutation arising from sorting characters
31          y,           # optional argument <y>
32          scpr;        # scalar product
33
34    # 1. Check the arguments.
35    if   Length( arg ) < 2 or Length( arg ) > 5
36       or not IsNearlyCharacterTable( arg[1] ) then
37      Error( "usage: ",
38             "LLL( <tbl>, <chars> [,<y>][,\"sort\"][,\"linearcomb\"] )" );
39    fi;
40
41    # 2. Get the arguments.
42    tbl:= arg[1];
43    characters:= arg[2];
44    if "sort" in arg then
45      sorted:= SortedCharacters( tbl, characters, "degree" );
46      perm:= Sortex( ShallowCopy( sorted ) )
47             / Sortex( ShallowCopy( characters ) );
48      characters:= sorted;
49    fi;
50    if IsBound( arg[3] ) and IsRat( arg[3] ) then
51      y:= arg[3];
52    else
53      y:= 3/4;
54    fi;
55
56    # 3. Call the LLL algorithm.
57    L:= AlgebraByGenerators( Rationals, [ TrivialCharacter( tbl ) ] );
58    if "linearcomb" in arg then
59      lllrb:= LLLReducedBasis( L, characters, y, "linearcomb" );
60    else
61      lllrb:= LLLReducedBasis( L, characters, y );
62    fi;
63
64    # 4. Make a new result record.
65    lll:= rec( irreducibles := [],
66               remainders   := [],
67               norms        := [] );
68
69    # 5. Sort the relations and transformation if necessary.
70    if IsBound( lllrb.relations ) then
71      lll.relations      := lllrb.relations;
72      lll.transformation := lllrb.transformation;
73      if IsBound( perm ) then
74        lll.relations      := List( lll.relations,
75                                    x -> Permuted( x, perm ) );
76        lll.transformation := List( lll.transformation,
77                                    x -> Permuted( x, perm ) );
78      fi;
79    fi;
80
81    # 6. Add the components used by the character table functions.
82    lll.irreducibles  := [];
83    lll.remainders    := [];
84    lll.norms         := [];
85    if IsBound( lllrb.transformation ) then
86      lll.irreddecomp := [];
87      lll.reddecomp   := [];
88    fi;
89
90    for i in [ 1 .. Length( lllrb.basis ) ] do
91
92      v:= lllrb.basis[i];
93      if v[1] < 0 then
94        v:= AdditiveInverse( v );
95        if IsBound( lllrb.transformation ) then
96          lll.transformation[i]:= AdditiveInverse( lll.transformation[i] );
97        fi;
98      fi;
99      scpr:= ScalarProduct( tbl, v, v );
100      if scpr = 1 then
101        Add( lll.irreducibles, Character( tbl, v ) );
102        if IsBound( lllrb.transformation ) then
103          Add( lll.irreddecomp, lll.transformation[i] );
104        fi;
105      else
106        Add( lll.remainders, VirtualCharacter( tbl, v ) );
107        Add( lll.norms, scpr );
108        if IsBound( lllrb.transformation ) then
109          Add( lll.reddecomp,   lll.transformation[i] );
110        fi;
111      fi;
112
113    od;
114
115    if not IsEmpty( lll.irreducibles ) then
116      Info( InfoCharacterTable, 2,
117            "LLL: ", Length( lll.irreducibles ), " irreducibles found" );
118    fi;
119
120    # 7. Sort `remainders' and `reddecomp' components if necessary.
121    if "sort" in arg then
122      sorted:= SortedCharacters( tbl, lll.remainders, "degree" );
123      perm:= Sortex( ShallowCopy( lll.remainders ) )
124             / Sortex( ShallowCopy( sorted ) );
125      lll.norms:= Permuted( lll.norms, perm );
126      lll.remainders:= sorted;
127      if "linearcomb" in arg then
128        lll.reddecomp:= Permuted( lll.reddecomp, perm );
129      fi;
130    fi;
131
132    # 7. Unbind components not used for characters.
133    Unbind( lll.transformation );
134
135    # 8. Return the result.
136    return lll;
137end );
138
139
140#############################################################################
141##
142#F  Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing> ] )
143##
144InstallGlobalFunction( Extract, function( arg )
145
146    local
147
148    # indices
149          i, j, k, l, n,
150    # input arrays
151          tbl, y, gram, missing,
152    # booleans
153          deeper, iszero, used, nullbegin, nonmissing,
154          maxnorm, minnorm, normbound, maxsum, solmat,
155          f, squares, sfind, choicecollect, sequence,
156          dependies, solcollect, sum, solcount, max, sumac, kmax,
157          solution,
158    # functions
159          next, zeroset, possiblies, update, correctnorm,
160          maxsquare, square, ident, begin;
161
162    # choosing next vector for combination
163    next := function( lines, solumat, acidx )
164    local i, j, solmat, testvec, idxback;
165
166    while acidx <= n and k + n - acidx >= kmax do
167       solmat := List( solumat, ShallowCopy );
168       if k = 0 then
169          i := acidx;
170          while i <= n and not begin( sequence[i] ) do
171             i := i + 1;
172          od;
173          if i > n then
174             nullbegin := true;
175          else
176             nullbegin := false;
177             if i > acidx then
178                idxback := sequence[i];
179                for j in [acidx + 1..1] do
180                   sequence[j] := sequence[j -1];
181                od;
182                sequence[acidx] := idxback;
183             fi;
184          fi;
185       fi;
186       k := k + 1;
187       f[k] := sequence[acidx];
188       testvec := [];
189       for i in [1..k] do
190          testvec[i] := gram[f[k]][f[i]];
191       od;
192       zeroset( solmat, testvec, lines );
193       acidx := acidx + 1;
194       possiblies( 1, solmat, testvec, acidx, lines );
195       k := k - 1;
196    od;
197    end;
198
199    # filling zero in places that fill already the conditions
200    zeroset := function( solmat, testvec, lines )
201    local i, j;
202
203    for i in [1..k-1] do
204       if testvec[i] = 0 then
205          for j in [1..lines] do
206             if solmat[j][i] <> 0 and not IsBound( solmat[j][k] ) then
207                solmat[j][k] := 0;
208             fi;
209          od;
210       fi;
211    od;
212    end;
213
214    # try and error for the chosen vector
215    possiblies := function( start, solmat, testvect, acidx, lines )
216    local i, j, remainder, toogreat, equal, solmatback, testvec;
217
218    testvec := ShallowCopy( testvect );
219    toogreat := false;
220    equal := true;
221    if k > 1 then
222       for i in [1..k-1] do
223          if testvec[i] < 0 then
224             toogreat := true;
225          fi;
226          if testvec[i] <> 0 then
227             equal := false;
228          fi;
229       od;
230       if testvec[k] < 0 then
231          toogreat := true;
232       fi;
233    else
234       if not nullbegin then
235          while start <= gram[f[k]][f[k]] and start < missing do
236             solmat[start][k] := 1;
237             start := start + 1;
238          od;
239          testvec[k] := 0;
240          if gram[f[k]][f[k]] > lines then
241             lines := gram[f[k]][f[k]];
242          fi;
243       else
244          lines := 0;
245       fi;
246    fi;
247    if not equal and not toogreat then
248       while start < lines and IsBound( solmat[start][k] ) do
249          start := start + 1;
250       od;
251       if start <= lines and not IsBound( solmat[start][k] ) then
252          solmat[start][k] := 0;
253          while not toogreat and not equal do
254             solmat[start][k] := solmat[start][k] + 1;
255             testvec := update( -1, testvec, start, solmat );
256             equal := true;
257             for i in [1..k-1] do
258                if testvec[i] < 0 then
259                   toogreat := true;
260                fi;
261                if testvec[i] <> 0 then
262                   equal := false;
263                fi;
264             od;
265             if testvec[k] < 0 then
266                toogreat := true;
267             fi;
268          od;
269       fi;
270    fi;
271    if equal and not toogreat then
272       solmatback := List( solmat, ShallowCopy );
273       for i in [1..missing] do
274          if not IsBound( solmat[i][k] ) then
275             solmat[i][k] := 0;
276          fi;
277       od;
278       correctnorm( testvec[k], solmat, lines + 1, testvec[k], acidx, lines );
279       solmat := solmatback;
280#T here was a 'Copy' call. WHY?
281    fi;
282    if k > 1 then
283       while start <= lines and solmat[start][k] > 0 do
284          solmat[start][k] := solmat[start][k] - 1;
285          testvec := update( 1, testvec, start, solmat );
286          solmatback := List( solmat, ShallowCopy );
287          zeroset( solmat, testvec, lines );
288          deeper := false;
289          for i in [1..k-1] do
290             if solmat[start][i] <> 0 then
291                deeper := false;
292                if testvec[i] = 0 then
293                   deeper := true;
294                else
295                   for j in [1..missing] do
296                      if solmat[j][i] <> 0 and not IsBound(solmat[j][k]) then
297                        deeper := true;
298                      fi;
299                   od;
300                fi;
301             fi;
302          od;
303          if deeper then
304             possiblies( start + 1, solmat, testvec, acidx, lines );
305          fi;
306          solmat := solmatback;
307#T here was a 'Copy' call. WHY?
308       od;
309    fi;
310    end;
311
312    # update the remaining conditions to fill
313    update := function( x, testvec, start, solmat )
314    local i;
315    for i in [1..k-1] do
316       if solmat[start][i] <> 0 then
317          testvec[i] := testvec[i] + solmat[start][i] * x;
318       fi;
319    od;
320    testvec[k] := testvec[k] - square( solmat[start][k] )
321                             + square( solmat[start][k] + x );
322    return testvec;
323    end;
324
325    # correct the norm if all other conditions are filled
326    correctnorm := function( remainder, solmat, pos, max, acidx, lines )
327    local i, r, newsol, ret;
328    if remainder = 0 and pos <= missing + 1 then
329       newsol := true;
330       for i in [1..solcount[k]] do
331          if ident( solcollect[k][i], solmat ) = missing then
332             newsol := false;
333          fi;
334       od;
335       if newsol then
336          if k > kmax then
337             kmax := k;
338          fi;
339          solcount[k] := solcount[k] + 1;
340          solcollect[k][solcount[k]] := [];
341          choicecollect[k][solcount[k]] := ShallowCopy( f );
342          for i in [1..Length( solmat )] do
343             solcollect[k][solcount[k]][i] := ShallowCopy( solmat[i] );
344          od;
345          if k = n and pos = missing + 1 then
346             ret := 0;
347          else
348             ret := max;
349             if k <> n then
350                next( lines, solmat, acidx );
351             fi;
352          fi;
353       else
354          ret := max;
355       fi;
356    else
357       if pos <= missing then
358          i := maxsquare( remainder, max );
359          while i > 0 do
360             solmat[pos][k] := i;
361             i := correctnorm( remainder-square( i ),
362                               solmat, pos+1, i, acidx, lines + 1);
363             i := i - 1;
364          od;
365          if i < 0 then
366             ret := 0;
367          else
368             ret := max;
369          fi;
370       else
371          ret := 0;
372       fi;
373    fi;
374    return ret;
375    end;
376
377    # compute the maximum squarenumber lower then given integer
378    maxsquare := function( value, max )
379    local i;
380
381    i := 1;
382    while square( i ) <= value and i <= max do
383          i := i + 1;
384    od;
385    return i-1;
386    end;
387
388    square := function( i )
389    if i = 0 then
390       return( 0 );
391    else
392       if not IsBound( squares[i] ) then
393          squares[i] := i * i;
394       fi;
395       return squares[i];
396    fi;
397    end;
398
399    ident := function( a, b )
400    # lists the identities of the two given sequences and counts them
401    local i, j, k, zi, zz, la, lb;
402    la := Length( a );
403    lb := Length( b );
404    zi := [];
405    zz := 0;
406    for i in [1..la] do
407       j := 1;
408       repeat
409          if a[i] = b[j] then
410             k :=1;
411             while k <= zz and j <> zi[k] do
412                k := k + 1;
413             od;
414             if k > zz then
415                zz := k;
416                zi[zz] := j;
417                j := lb;
418             fi;
419          fi;
420          j := j + 1;
421       until j > lb;
422    od;
423    return( zz );
424    end;
425
426    # looking for character that can stand at the beginning
427    begin := function( i )
428    local ind;
429    if y = [] or gram[i][i] < 4 then
430       return true;
431    else
432       if IsBound( ComputedPowerMaps( tbl )[2] ) then
433          if ForAll( ComputedPowerMaps( tbl )[2], IsInt ) then
434#T ??
435             ind := AbsInt( Indicator( tbl, [y[i]], 2 )[1]);
436             if gram[i][i] - 1 <= ind
437             or ( gram[i][i] = 4 and ind = 1 ) then
438                return true;
439             fi;
440          fi;
441       fi;
442    fi;
443    return false;
444    end;
445
446    # check input parameters
447    if IsNearlyCharacterTable( arg[1] ) then
448       tbl := arg[1];
449    else
450       Error( "first argument must be character table\n \
451    usage: Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing>] )" );
452    fi;
453    if IsBound( arg[2] ) and IsList( arg[2] ) and IsList( arg[2][1] ) then
454       y := List( arg[2], ShallowCopy );
455    else
456       Error( "second argument must be list of reducible characters\n \
457    usage: Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing>] )" );
458    fi;
459    if IsBound( arg[2] ) and IsList( arg[3] ) and IsList( arg[3][1] ) then
460       gram := List( arg[3], ShallowCopy );
461    else
462       Error( "third argument must be gram-matrix of reducible characters\n \
463    usage: Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing>] )" );
464    fi;
465    n := Length( gram );
466    if IsBound( arg[4] ) and IsInt( arg[4] ) then
467       missing := arg[4];
468    else
469       missing := n;
470       nonmissing := true;
471    fi;
472
473    # main program
474    maxnorm := 0;
475    minnorm := gram[1][1];
476    normbound := [];
477    maxsum := [];
478    solcollect := [];
479    choicecollect := [];
480    sum := [];
481    solmat := [];
482    used := [];
483    solcount := [];
484    sfind := [];
485    f := [];
486    squares := [];
487    kmax := 0;
488    for i in [1..missing] do
489       solmat[i] := [];
490    od;
491    for i in [1..n] do
492       solcount[i] := 0;
493       used[i] := false;
494       solcollect[i] := [];
495       choicecollect[i] := [];
496    od;
497    for i in [1..n] do
498       if gram[i][i] > maxnorm then
499          maxnorm := gram[i][i];
500       else
501          if gram[i][i] < minnorm then
502             minnorm := gram[i][i];
503          fi;
504       fi;
505    od;
506    j := 0;
507    for i in [minnorm..maxnorm] do
508       k := 1;
509       while k <= n and gram[k][k] <> i do
510          k := k + 1;
511       od;
512       if k <= n then
513           j := j + 1;
514           normbound[j] := rec( norm:=i, first:=k, last:=0 );
515           if k = n then
516              normbound[j].last := k;
517           else
518              k := n;
519              while gram[k][k] <> i and k > 0 do
520                 k := k - 1;
521              od;
522              if k > 0 then
523                 normbound[j].last := k;
524              fi;
525           fi;
526       fi;
527    od;
528    for j in [1..Length( normbound )] do
529       maxsum[j] := 0;
530       for i in [normbound[j].first..normbound[j].last] do
531          if gram[i][i] = normbound[j].norm then
532             sum[i] := 0;
533             for k in [1..n] do
534                sum[i] := sum[i] + gram[i][k];
535             od;
536             if sum[i] > maxsum[j] then
537                maxsum[j] := sum[i];
538             fi;
539          fi;
540       od;
541    od;
542    k := 1;
543    sequence := [];
544    i:= 1;
545    while i <= Length( normbound ) do
546       max := maxsum[i];
547       sumac := 0;
548       for j in [normbound[i].first..normbound[i].last] do
549          if gram[j][j] = normbound[i].norm and sum[j] > sumac
550          and sum[j] <= max and not used[j] then
551             sequence[k] := j;
552             sumac := sum[j];
553          fi;
554       od;
555       if IsBound( sequence[k] ) then
556          max := sumac;
557          used[sequence[k]] := true;
558          k := k + 1;
559       else
560          i := i + 1;
561       fi;
562    od;
563    k := 0;
564    next( 1, solmat, 1 );
565    solution := rec( solution := [], choice := choicecollect[kmax] );
566    for i in [1..solcount[kmax]] do
567       solution.solution[i] := [];
568       l := 0;
569       for j in [1..missing] do
570          iszero := true;
571          for k in [1..kmax] do
572             if solcollect[kmax][i][j][k] <> 0 then
573                iszero := false;
574             fi;
575          od;
576          if not iszero then
577             l := l + 1;
578             solution.solution[i][l] := solcollect[kmax][i][j];
579          fi;
580       od;
581    od;
582    return( solution );
583end );
584
585
586#############################################################################
587##
588#F  Decreased( <tbl>, <chars>, <decompmat>, [ <choice> ] )
589##
590InstallGlobalFunction( Decreased, function( arg )
591    local
592        # indices
593          m, n, m1, n1, i, i1, i2, i3, i4, j, jj, j1, j2, j3,
594        # booleans
595          ende1, ende2, ok, change, delline, delcolumn,
596        # help fields
597          deleted, kgv, l1, l2, l3, dim, ident,
598        # matrices
599          invmat, remmat, remmat2, solmat, nonzero,
600        # double-indices
601          columnidx, lineidx, system, components, compo2,
602        # output-fields
603          sol, red, redcount, irred,
604        # help fields
605          IRS, SFI, lc, nc, char, char1, entries,
606        # input fields
607          tbl, y, choice,
608        # functions
609          Idxset, Identset, Invadd, Invmult, Nonzeroset;
610
611    Idxset := function()
612    # update indices
613    local i1, j1;
614    i1 := 0;
615    for i in [1..m] do
616       if not delline[i] then
617          i1 := i1 + 1;
618          lineidx[i1] := i;
619       fi;
620    od;
621    m1 := i1;
622    j1 := 0;
623    for j in [1..n] do
624       if not delcolumn[j] then
625          j1 := j1 + 1;
626          columnidx[j1] := j;
627       fi;
628    od;
629    n1 := j1;
630    end;
631
632    Identset := function( veca, vecb )
633#T just one place where this is called ...
634    # count identities of veca and vecb and store "non-identities"
635    local la, lb, i, j, n, nonid, nic, r;
636    n := 0;
637    la := Length( veca );
638    lb := Length( vecb );
639    j := 1;
640    nonid := [];
641    nic := 0;
642    for i in [1..la] do
643       while j <= lb and veca[i] > vecb[j] do
644          nic := nic + 1;
645          nonid[nic] := vecb[j];
646          j := j + 1;
647       od;
648       if j <= lb and veca[i] = vecb[j] then
649          n := n + 1;
650          j := j + 1;
651       fi;
652    od;
653    while j <= lb do
654       nic := nic + 1;
655       nonid[nic] := vecb[j];
656       j := j + 1;
657    od;
658    r := rec( nonid := nonid, id := n  );
659    return( r );
660    end;
661
662    Invadd := function( j1, j2, l )
663    # addition of two lines of invmat
664    local i;
665    for i in [1..n] do
666       if invmat[i][j2] <> 0 then
667          invmat[i][j1] := invmat[i][j1] - l * invmat[i][j2];
668       fi;
669    od;
670    end;
671
672    Invmult := function( j1, l )
673    # multiply line of invmat
674    local i;
675    if l <> 1 then
676       for i in [1..n] do
677          if invmat[i][j1] <> 0 then
678             invmat[i][j1] := invmat[i][j1] * l;
679          fi;
680       od;
681    fi;
682    end;
683
684    Nonzeroset := function( j )
685    # entries <> 0 in j-th column of 'solmat'
686
687    local i, j1;
688    nonzero[j] := [];
689    j1 := 0;
690    for i in [1..m] do
691       if solmat[i][j] <> 0 then
692          j1 := j1 + 1;
693          nonzero[j][j1] := i;
694       fi;
695    od;
696    entries[j] := j1;
697    end;
698
699    # check input parameters
700    if Length( arg ) < 3 or Length( arg ) > 4 then
701      Error( "usage: Decreased( <tbl>, <list of char>,\n",
702             "<decomposition matrix>, [<choice>] )" );
703    fi;
704
705    if IsNearlyCharacterTable( arg[1] ) then
706       tbl := arg[1];
707    else
708       Error( "first argument must be a nearly character table\n",
709              "usage: Decreased( <tbl>, <list of char>,\n",
710              "<decomposition matrix>, [<choice>] )" );
711    fi;
712    if IsList( arg[2] ) and IsList( arg[2][1] ) then
713       y := arg[2];
714    else
715       Error( "second argument must be list of characters\n",
716              "usage: Decreased( <tbl>, <list of char>,\n",
717              "<decomposition matrix>, [<choice>] )" );
718    fi;
719    if IsList( arg[3] ) and IsList( arg[3][1] ) then
720       solmat := List( arg[3], ShallowCopy );
721    else
722       Error( "third argument must be decomposition matrix\n",
723              "usage: Decreased( <tbl>, <list of char>,\n",
724              "<decomposition-matrix>, [<choice>] )" );
725    fi;
726    if not IsBound( arg[4] ) then
727       choice := [ 1 .. Length( y ) ];
728    elif IsList( arg[4] ) then
729       choice := arg[4];
730    else
731       Error( "forth argument contains choice of characters\n",
732           "usage: Decreased( <tbl>, <list of char>,\n",
733           "<decomposition-matrix>, [<choice>] )" );
734    fi;
735
736    # initialisations
737    lc := Length( y[1] );
738    nc := [];
739    for i in [1..lc] do
740       nc[i] := 0;
741    od;
742    columnidx := [];
743    lineidx   := [];
744    nonzero   := [];
745    entries   := [];
746    delline   := [];
747    delcolumn := [];
748
749    # number of lines
750    m := Length( solmat );
751
752    # number of columns
753    n := Length( solmat[1] );
754    invmat := IdentityMat( n );
755    for i in [1..m] do
756       delline[i] := false;
757    od;
758    for j in [1..n] do
759       delcolumn[j] := false;
760    od;
761    i := 1;
762
763    # check lines for information
764    while i <= m do
765       if not delline[i] then
766          entries[i] := 0;
767          for j in [1..n] do
768             if solmat[i][j] <> 0 and not delcolumn[j] then
769                entries[i] := entries[i] + 1;
770                if entries[i] = 1 then
771                   nonzero[i] := j;
772                fi;
773             fi;
774          od;
775          if entries[i] = 1 then
776             delcolumn[nonzero[i]] := true;
777             delline[i] := true;
778             j := 1;
779             while j < i and solmat[j][nonzero[i]] = 0 do
780                j := j + 1;
781             od;
782             if j < i then
783                i := j;
784             else
785                i := i + 1;
786             fi;
787          else
788             if entries[i] = 0 then
789                delline[i] := true;
790             fi;
791             i := i + 1;
792          fi;
793       else
794          i := i + 1;
795       fi;
796    od;
797    Idxset();
798
799    deleted := m - Length(lineidx);
800    for j in [1..n] do
801       Nonzeroset( j );
802    od;
803    ende1 := false;
804    while not ende1 and deleted < m do
805       j := 1;
806
807    # check solo-entry-columns
808       while j <= n do
809          if entries[j] = 1 then
810             change := false;
811             for jj in [1..n] do
812                if (delcolumn[j] and delcolumn[jj])
813                or not delcolumn[j] then
814                   if solmat[nonzero[j][1]][jj] <> 0 and jj <> j then
815                      change := true;
816                      kgv := Lcm( solmat[nonzero[j][1]][j],
817                              solmat[nonzero[j][1]][jj] );
818                      l1 := kgv / solmat[nonzero[j][1]][jj];
819                      Invmult( jj, l1 );
820                      for i1 in [1..Length( nonzero[jj] )] do
821                         solmat[nonzero[jj][i1]][jj]
822                               := solmat[nonzero[jj][i1]][jj] * l1;
823                      od;
824                      Invadd( jj, j, kgv/solmat[nonzero[j][1]][j] );
825                      solmat[nonzero[j][1]][jj] := 0;
826                      Nonzeroset( jj );
827                   fi;
828                fi;
829             od;
830             if not delline[nonzero[j][1]] then
831                delline[nonzero[j][1]] := true;
832                delcolumn[j] := true;
833                deleted := deleted + 1;
834                Idxset();
835             fi;
836             if change then
837                j := 1;
838             else
839                j := j + 1;
840             fi;
841          else
842             j := j + 1;
843          fi;
844       od;
845
846    # search for Equality-System
847    # system :     chosen columns
848    # components : entries <> 0 in the chosen columns
849       dim := 2;
850       change := false;
851       ende2 := false;
852       while dim <= n1 and not ende2 do
853          j3 := 1;
854          while j3 <= n1 and not ende2 do
855             j2 := j3;
856             j1 := 0;
857             system := [];
858             components := [];
859             while j2 <= n1 do
860                while j2 <= n1 and entries[columnidx[j2]] > dim do
861                   j2 := j2 + 1;
862                od;
863                if j2 <= n1 then
864                   if j1 = 0 then
865                      j1 := 1;
866                      system[j1] := columnidx[j2];
867                      components := ShallowCopy( nonzero[columnidx[j2]] );
868                   else
869                      ident := Identset( components, nonzero[columnidx[j2]] );
870                      if dim - Length( components ) >= entries[columnidx[j2]]
871                             - ident.id then
872                         j1 := j1 + 1;
873                         system[j1] := columnidx[j2];
874                         if ident.id < entries[columnidx[j2]] then
875                            compo2 := ShallowCopy( components );
876                            components := [];
877                            i1 := 1;
878                            i2 := 1;
879                            i3 := 1;
880
881    # append new entries to "components"
882                            while i1 <= Length( ident.nonid )
883                               or i2 <= Length( compo2 ) do
884                               if i1 <= Length( ident.nonid ) then
885                                  if i2 <= Length( compo2 ) then
886                                     if ident.nonid[i1] < compo2[i2] then
887                                        components[i3] := ident.nonid[i1];
888                                        i1 := i1 + 1;
889                                     else
890                                        components[i3] := compo2[i2];
891                                        i2 := i2 + 1;
892                                     fi;
893                                  else
894                                     components[i3] := ident.nonid[i1];
895                                     i1 := i1 + 1;
896                                  fi;
897                               else
898                                  if i2 <= Length( compo2 ) then
899                                     components[i3] := compo2[i2];
900                                     i2 := i2 + 1;
901                                  fi;
902                               fi;
903                               i3 := i3 + 1;
904                            od;
905                         fi;
906                      fi;
907                   fi;
908                   j2 := j2 + 1;
909                fi;
910             od;
911
912    # try to solve system with Gauss
913             if  Length( system ) > 1 then
914                for i1 in [1..Length( components )] do
915                   i2 := 1;
916                   repeat
917                      ok := true;
918                      if solmat[components[i1]][system[i2]] = 0 then
919                         ok := false;
920                      else
921                         for i3 in [1..i1-1] do
922                            if solmat[components[i3]][system[i2]] <> 0 then
923                               ok := false;
924                            fi;
925                         od;
926                      fi;
927                      if not ok then
928                         i2 := i2 + 1;
929                      fi;
930                   until ok or i2 > Length( system );
931                   if ok then
932                      for i3 in [1..Length( system )] do
933                         if i3 <> i2
934                            and solmat[components[i1]][system[i3]] <> 0 then
935                            change := true;
936                            kgv := Lcm( solmat[components[i1]][system[i3]],
937                                       solmat[components[i1]][system[i2]] );
938                            l2 := kgv / solmat[components[i1]][system[i2]];
939                            l3 := kgv / solmat[components[i1]][system[i3]];
940                            for i4 in [1..Length( nonzero[system[i3]] )] do
941                               solmat[nonzero[system[i3]][i4]][system[i3]]
942                            := solmat[nonzero[system[i3]][i4]][system[i3]]*l3;
943                            od;
944                            Invmult( system[i3], l3 );
945                            for i4 in [1..Length( nonzero[system[i2]] )] do
946                               solmat[nonzero[system[i2]][i4]][system[i3]]
947                               := solmat[nonzero[system[i2]][i4]][system[i3]]
948                            -  solmat[nonzero[system[i2]][i4]][system[i2]]*l2;
949                            od;
950                            Invadd( system[i3], system[i2], l2 );
951                            Nonzeroset( system[i3] );
952                            if entries[system[i3]] = 0 then
953                               delcolumn[system[i3]] := true;
954                               Idxset();
955                            fi;
956                         fi;
957                      od;
958                   fi;
959                od;
960
961   # check for columns with only one entry <> 0
962                for i1 in [1..Length( system )] do
963                   if entries[system[i1]] = 1 then
964                      ende2 := true;
965                   fi;
966                od;
967                if not ende2 then
968                   j3 := j3 + 1;
969                fi;
970             else
971                j3 := j3 + 1;
972             fi;
973          od;
974          dim := dim + 1;
975       od;
976       if dim > n1 and not change and j3 > n1 then
977          ende1 := true;
978       fi;
979    od;
980
981    # check, if
982    #    the transformation of solmat allows computation of new irreducibles
983    remmat := [];
984    for i in [1..m] do
985       remmat[i] := [];
986       delline[i] := true;
987    od;
988    redcount := 0;
989    red := [];
990    irred := [];
991    j := 1;
992    sol := true;
993    while j <= n and sol do
994
995    # computation of character
996       char := ShallowCopy( nc );
997       for i in [1..n] do
998          if invmat[i][j] <> 0 then
999             char := char + invmat[i][j] * y[choice[i]];
1000          fi;
1001       od;
1002
1003    # probably irreducible ==> has to pass tests
1004       if entries[j] = 1 then
1005          if solmat[nonzero[j][1]][j] <> 1 then
1006             char1 := char/solmat[nonzero[j][1]][j];
1007          else
1008             char1 := char;
1009          fi;
1010          if char1[1] < 0 then
1011             char1 := - char1;
1012          fi;
1013
1014    # is 'char1' real?
1015          IRS := ForAll( char1, x -> GaloisCyc(x,-1) = x );
1016
1017   # Frobenius Schur indicator
1018          if IsBound( ComputedPowerMaps( tbl )[2] )
1019             and ForAll( ComputedPowerMaps( tbl )[2], IsInt ) then
1020#T ??
1021            SFI:= Indicator( tbl, [ char1 ], 2 )[1];
1022          else
1023            SFI:= Unknown();
1024            Info( InfoCharacterTable, 2,
1025                  "Decreased: 2nd power map not available or not unique,\n",
1026                  "#I  no test with 'Indicator'" );
1027          fi;
1028
1029   # test if 'char1' can be an irreducible character
1030          if    char1[1] = 0
1031             or ForAny( char1, x -> not IsCycInt(x) )
1032             or ScalarProduct( tbl, char1, char1 ) <> 1
1033             or ( IsCyc( SFI ) and ( ( IRS and AbsInt( SFI ) <> 1 ) or
1034                                     ( not IRS and SFI <> 0 ) ) )   then
1035            Info( InfoCharacterTable, 2,
1036                  "Decreased : computation of ",
1037                  Ordinal( Length( irred ) + 1 ), " character failed" );
1038            return fail;
1039          else
1040
1041    # irreducible character found
1042            Add( irred, Character( tbl, char1 ) );
1043          fi;
1044       else
1045
1046    # what a pity (!), some reducible character remaining
1047          if char[1] < 0 then
1048             char := - char;
1049          fi;
1050          if char <> nc then
1051             redcount := redcount + 1;
1052             red[redcount] := ClassFunction( tbl, char );
1053             for i in [1..m] do
1054                remmat[i][redcount] := solmat[i][j];
1055                if solmat[i][j] <> 0 then
1056                   delline[i] := false;
1057                fi;
1058             od;
1059          fi;
1060       fi;
1061       j := j+1;
1062    od;
1063    i1 := 0;
1064    remmat2 := [];
1065    for i in [1..m] do
1066       if not delline[i] then
1067          i1 := i1 + 1;
1068          remmat2[i1] := remmat[i];
1069       fi;
1070    od;
1071    return rec( irreducibles := irred,
1072                remainders   := red,
1073                matrix       := remmat2 );
1074end );
1075
1076
1077#############################################################################
1078##
1079#F  OrthogonalEmbeddingsSpecialDimension( <tbl>, <reducibles>, <grammat>,
1080#F                                        [, \"positive\"], <dim> )
1081##
1082InstallGlobalFunction( OrthogonalEmbeddingsSpecialDimension, function ( arg )
1083    local  red, dim, reducibles, matrix, tbl, emb, dec, i, s, irred;
1084    # check input
1085    if Length( arg ) < 4 then
1086       Error( "please specify desired dimension\n",
1087              "usage : OrthogonalE...( <tbl>, <reducibles>,\n",
1088              "<gram-matrix>[, \"positive\" ], <dim> )" );
1089    fi;
1090    if IsInt( arg[4] ) then
1091       dim := arg[4];
1092    else
1093       if IsBound( arg[5] ) then
1094          if IsInt( arg[5] ) then
1095             dim := arg[5];
1096          else
1097       Error( "please specify desired dimension\n",
1098              "usage : Orthog...( <tbl>, < reducibles >,\n",
1099              "< gram-matrix >, [, <\"positive\"> ], < integer > )" );
1100          fi;
1101       fi;
1102    fi;
1103    tbl := arg[1];
1104    reducibles := arg[2];
1105    if Length( arg ) = 4 then
1106       emb := OrthogonalEmbeddings( arg[3], arg[4] );
1107    else
1108       emb := OrthogonalEmbeddings( arg[3], arg[4], arg[5] );
1109    fi;
1110    s := [];
1111    for i in [1..Length(emb.solutions)] do
1112       if Length( emb.solutions[i] ) = dim then
1113          Add( s, emb.vectors{ emb.solutions[i] } );
1114       fi;
1115    od;
1116    dec:= List( s, x -> Decreased( tbl, reducibles, x ) );
1117    dec:= Filtered( dec, x -> x <> fail );
1118    if dec = [] then
1119      Info( InfoCharacterTable, 2,
1120            "OrthogonalE...: no embedding corresp. to characters" );
1121      return rec( irreducibles:= [], remainders:= reducibles );
1122    fi;
1123    irred:= Set( dec[1].irreducibles );
1124    for i in [2..Length(dec)] do
1125       IntersectSet( irred, dec[i].irreducibles );
1126    od;
1127    red:= ReducedClassFunctions( tbl, irred, reducibles );
1128    Append( irred, red.irreducibles );
1129    return rec( irreducibles:= irred, remainders:= red.remainders );
1130end );
1131
1132
1133#############################################################################
1134##
1135#F  DnLattice( <tbl>, <g1>, <y1> )
1136##
1137InstallGlobalFunction( DnLattice, function( tbl, g1, y1 )
1138    local
1139    # indices
1140      i, i1, j, j1, k, k1, l, next,
1141    # booleans
1142      empty, change, used, addable, SFIbool,
1143    # dimensions
1144      m, n,
1145    # help fields
1146      found, foundpos,
1147      z, zw, nullcount, nullgenerate,
1148      maxentry, max, ind, irred, irredcount, red,
1149      blockcount, blocks, perm, addtest, preirred,
1150    # Gram matrix
1151      g, gblock,
1152    # characters
1153      y, y2,
1154    # variables for recursion
1155      root, rootcount, solution, ligants, ligantscount, glblock, begin,
1156      depth, choice, ende, sol,
1157    # functions
1158      callreduced, nullset, maxset, Search, Add, DnSearch, test;
1159
1160    # counts zeroes in given line
1161    nullset := function( g, i )
1162    local j;
1163
1164    nullcount[ i ] := 0;
1165    for j in [ 1..n ] do
1166       if g[ j ] = 0 then
1167          nullcount[ i ] := nullcount[ i ] + 1;
1168       fi;
1169    od;
1170    end;
1171
1172    # searches line with most non-zero-entries
1173    maxset := function( )
1174    local i;
1175
1176    maxentry := 1;
1177    max := n;
1178    for i in [ 1..n ] do
1179       if nullcount[ i ] < max then
1180          max := nullcount[ i ];
1181          maxentry := i;
1182       fi;
1183    od;
1184    end;
1185
1186    # searches lines to add in order to produce zeroes
1187    Search := function( j )
1188    local signum;
1189
1190    nullgenerate := 0;
1191    if g[ j ][ maxentry ] > 0 then
1192       signum := -1;
1193       for k in [ 1..n ] do
1194          if k <> maxentry and k <> j then
1195             if g[ maxentry ][ k ] <> 0 then
1196                if g[ j ][ k ] = g[ maxentry ][ k ] then
1197                   nullgenerate := nullgenerate + 1;
1198                else
1199                   nullgenerate := nullgenerate - 1;
1200                fi;
1201             fi;
1202          fi;
1203       od;
1204    else
1205       if g[ j ][ maxentry ] < 0 then
1206          signum := 1;
1207          for k in [ 1..n ] do
1208             if k <> maxentry and k <> j then
1209                if g[ maxentry ][ k ] <> 0 then
1210                   if g[ j ][ k ] = -g[ maxentry ][ k ] then
1211                      nullgenerate := nullgenerate + 1;
1212                   else
1213                      nullgenerate := nullgenerate - 1;
1214                   fi;
1215                fi;
1216             fi;
1217          od;
1218       fi;
1219    fi;
1220    if nullgenerate > 0 then
1221       change := true;
1222       Add( j, maxentry );
1223       j := j + 1;
1224    fi;
1225    end;
1226
1227    # adds two lines/columns
1228    Add := function( i, j )
1229    local k;
1230
1231       y[ i ] := y[ i ] - g[ i ][ j ] * y[ j ];
1232       g[ i ] := g[ i ] - g[ i ][ j ] * g[ j ];
1233       for k in [ 1..i-1 ] do
1234          g[ k ][ i ] := g[ i ][ k ];
1235       od;
1236       g[ i ][ i ] := 2;
1237       for k in [ i+1..n ] do
1238          g[ k ][ i ] := g[ i ][ k ];
1239       od;
1240    end;
1241
1242    # backtrack-search for dn-lattice
1243    DnSearch := function( begin, depth, oldchoice )
1244    local connections, connect, i1, j1, choice, found;
1245
1246    choice := ShallowCopy( oldchoice );
1247    if depth = 3 then
1248       # d4-lattice found !!!
1249       solution := 1;
1250       ende := true;
1251       if n > 4 then
1252          i1 := 0;
1253          found := false;
1254          while not found and i1 < n do
1255             i1 := i1 + 1;
1256             if i1 <> root[ j ] and i1 <> choice[ 1 ]
1257             and i1 <> choice[ 2 ] and i1 <> choice[ 3 ] then
1258                connections := 0;
1259                for j1 in [1..3] do
1260                   if gblock[ i1 ][ choice[ j1 ] ] <> 0 then
1261                      connections := connections + 1;
1262                      connect := choice[ j1 ];
1263                   fi;
1264                od;
1265                if connections = 1 then
1266                   found := true;
1267                   choice[ 4 ] := connect;
1268                   solution := solution + 1;
1269                fi;
1270             fi;
1271             i1 := i1 + 1;
1272          od;
1273       fi;
1274       sol := choice;
1275    else
1276       i1 := begin;
1277       while not ende and i1 <= ligantscount do
1278          found := true;
1279          for j1 in [1..depth] do
1280             if gblock[ ligants[ i1 ] ][ choice[ j1 ] ] <> 0 then
1281                found := false;
1282             fi;
1283          od;
1284          if found then
1285             depth := depth + 1;
1286             choice[ depth ] := ligants[ i1 ];
1287             DnSearch( i1 + 1, depth, choice );
1288             depth := depth - 1;
1289          else
1290             i1 := i1 + 1;
1291          fi;
1292          if ligantscount - i1 + 1 + depth < 3 then
1293             ende := true;
1294          fi;
1295       od;
1296    fi;
1297    end;
1298
1299    test := function(z)
1300    # some tests for the found characters
1301    local result, IRS, SFI, i1, y1, ind, testchar;
1302    testchar := z/2;
1303    result := true;
1304    IRS := ForAll( testchar, x -> GaloisCyc(x,-1) = x );
1305    if IsBound( ComputedPowerMaps( tbl )[2] ) then
1306       if ForAll( ComputedPowerMaps( tbl )[2], IsInt ) then
1307          SFI := Indicator( tbl, [testchar], 2 )[1];
1308          SFIbool := true;
1309       else
1310          Info( InfoCharacterTable, 2,
1311                "DnLattice: 2nd power map not available or not unique,\n",
1312                "#I            cannot test with Indicator" );
1313          SFIbool := false;
1314       fi;
1315    else
1316      Info( InfoCharacterTable, 2,
1317            "DnLattice: 2nd power map not availabe\n",
1318            "#I            cannot test with Indicator" );
1319      SFIbool := false;
1320    fi;
1321    if SFIbool then
1322       if ForAny( testchar, x -> IsRat(x) and not IsInt(x) )
1323          or ScalarProduct( tbl, testchar, testchar ) <> 1
1324          or testchar[1] = 0
1325          or ( IRS and AbsInt( SFI ) <> 1 )
1326          or ( not IRS and SFI <> 0 ) then
1327         result := false;
1328       fi;
1329    else
1330       if ForAny( testchar, x -> IsRat(x) and not IsInt(x) )
1331          or ScalarProduct( tbl, testchar, testchar ) <> 1
1332          or testchar[1] = 0 then
1333         result := false;
1334       fi;
1335    fi;
1336    return result;
1337    end;
1338
1339    # reduce whole lattice with the found irreducible
1340    callreduced := function()
1341    z[ 1 ] := z[ 1 ]/ 2 ;
1342    if ScalarProduct( tbl, z[ 1 ], z[ 1 ] ) = 1 then
1343       irredcount := irredcount + 1;
1344       if z[ 1 ][ 1 ] > 0 then
1345          irred[ irredcount ] := Character( tbl, z[ 1 ] );
1346       else
1347          irred[ irredcount ] := Character( tbl, -z[ 1 ] );
1348       fi;
1349       y1 := y{ [ blocks.begin[i] .. blocks.ende[i] ] };
1350       red := ReducedClassFunctions( tbl, z, y1 );
1351       Append( irred, List( red.irreducibles, x -> Character( tbl, x ) ) );
1352       irredcount := Length( irred );
1353       y2 := Concatenation( y2, red.remainders );
1354    fi;
1355    end;
1356
1357    # check input parameters
1358    if not IsNearlyCharacterTable( tbl ) then
1359       Error( "first argument must be a nearly character table\n",
1360              "usage: DnLattice( <tbl>, <gram-matrix>, <reducibles> )" );
1361    fi;
1362    empty := false;
1363    if not IsEmpty( g1 ) then
1364      if IsList( g1 ) and IsBound( g1[1] ) and IsList( g1[1] ) then
1365        g := List( g1, ShallowCopy );
1366      else
1367        Error( "second argument must be Gram matrix of characters\n",
1368               "usage: DnLattice( <tbl>, <gram-matrix>, <reducibles> )" );
1369      fi;
1370    else
1371      empty := true;
1372    fi;
1373    if not IsEmpty( y1 ) then
1374      if IsList( y1 ) and IsBound( y1[1] ) and IsList( y1[1] ) then
1375        y := List( y1, ShallowCopy );
1376      else
1377        Error( "third argument must be list of reducible characters\n",
1378               "usage: DnLattice( <tbl>, <gram-matrix>, <reducibles> )" );
1379      fi;
1380    else
1381      empty := true;
1382    fi;
1383    y2        := [  ];
1384    irred     := [  ];
1385
1386    if not empty then
1387
1388    n := Length( y );
1389    for i in [1..n] do
1390       if g[i][i] <> 2 then
1391          Error( "reducible characters don't have norm 2\n",
1392                "usage: DnLattice( <tbl>, <gram-matrix>, <reducibles> )" );
1393       fi;
1394    od;
1395    # initialisations
1396    z         := [  ];
1397    used      := [  ];
1398    next      := [  ];
1399    nullcount := [  ];
1400    m := Length( y[ 1 ] );
1401    for i in [1..n] do
1402       used[i] := false;
1403    od;
1404    blocks := rec( begin := [ ], ende := [ ] );
1405    blockcount   := 0;
1406    irredcount   := 0;
1407    change       := true;
1408    while change do
1409       change := false;
1410       for i in [ 1..n ] do
1411          nullset( g[ i ], i );
1412       od;
1413       maxset( );
1414       while max < n-2 and not change do
1415          while maxentry <= n and not change do
1416             if nullcount[ maxentry ] <> max then
1417                maxentry := maxentry + 1;
1418             else
1419                j := 1;
1420                while j < maxentry and not change do
1421                   Search( j );
1422                   j := j + 1;
1423                od;
1424                j := maxentry + 1;
1425                while j <= n and not change do
1426                   Search( j );
1427                   j := j + 1;
1428                od;
1429                if not change then
1430                   maxentry := maxentry + 1;
1431                fi;
1432             fi;
1433          od;
1434          if not change then
1435             max := max + 1;
1436             maxentry := 1;
1437          fi;
1438       od;
1439
1440    # 2 step-search in order to produce zeroes
1441    # 2_0_Box-Method
1442       change := false;
1443       i := 1;
1444       while i <= n and not change do
1445          while i <= n and nullcount[ i ] > n-3 do
1446             i := i + 1;
1447          od;
1448          if i <= n then
1449             j := 1;
1450             while j <= n and not change do
1451                while j <= n and g[ i ][ j ] <> 0 do
1452                   j := j + 1;
1453                od;
1454                if j <= n then
1455                   i1 := 1;
1456                   while i1 <= n and not change do
1457                      while i1 <= n
1458                      and ( i1 = i or i1 = j or g[ i1 ][ j ] = 0 ) do
1459                         i1 := i1 + 1;
1460                      od;
1461                      if i1 <= n then
1462                         addtest := g[ i ] - g[ i ][ i1 ] * g[ i1 ];
1463                         nullgenerate := 0;
1464                         addable := true;
1465                         for k in [ 1..n ] do
1466                            if addtest[ k ] = 0 then
1467                               nullgenerate := nullgenerate + 1;
1468                            else
1469                               if AbsInt( addtest[ k ] ) > 1 then
1470                                  addable := false;
1471                               fi;
1472                            fi;
1473                         od;
1474                         if addable then
1475                            nullgenerate := nullgenerate - nullcount[ i ];
1476                            for k in [ 1..n ] do
1477                               if k <> i and k <> j then
1478                                  if addtest[ k ]
1479                                     = addtest[ j ] * g[ j ][ k ] then
1480                                     if g[ j ][ k ] <> 0 then
1481                                        nullgenerate := nullgenerate + 1;
1482                                     fi;
1483                                  else
1484                                     if addtest[ k ] <> 0 then
1485                                        if g[ j ][ k ] = 0 then
1486                                           nullgenerate := nullgenerate - 1;
1487                                        else
1488                                           addable := false;
1489                                        fi;
1490                                     fi;
1491                                  fi;
1492                               fi;
1493                            od;
1494                            if nullgenerate > 0 and addable then
1495                               Add( i, i1 );
1496                               Add( j, i );
1497                               change := true;
1498                            fi;
1499                         fi;
1500                         i1 := i1 + 1;
1501                      fi;
1502                   od;
1503                   j := j + 1;
1504                fi;
1505             od;
1506             i := i + 1;
1507          fi;
1508       od;
1509    od;
1510    i := 1;
1511    j := 0;
1512    next[ 1	] := 1;
1513    while j < n do
1514       blockcount := blockcount + 1;
1515       blocks.begin[ blockcount ] := i;
1516       l := 0;
1517       used[ next [ i ] ] := true;
1518       j := j + 1;
1519       y2[ j ] := y[ next [ i ] ];
1520       while l >= 0 do
1521          for k in [ 1..n ] do
1522             if g[ next[ i ] ][ k ] <> 0 and not used[ k ] then
1523                l := l + 1;
1524                next[ i + l ] := k;
1525                j := j + 1;
1526                y2[ j ] := y[ k ];
1527                used[ k ] := true;
1528             fi;
1529          od;
1530          i := i + 1;
1531          l := l - 1;
1532       od;
1533       blocks.ende[ blockcount ] := i - 1;
1534       k := 1;
1535       while k <= n and used[ k ] do
1536          k := k + 1;
1537       od;
1538       if k <= n then
1539          next[i] := k;
1540       fi;
1541    od;
1542    perm := PermList( next )^-1;
1543    for i in [1..n] do
1544       g[i] := Permuted( g[i], perm );
1545    od;
1546    g := Permuted( g, perm );
1547    y := y2;
1548    y2 := [  ];
1549
1550    # search for d4/d5 - lattice
1551    for i in [1..blockcount] do
1552       n := blocks.ende[ i ] - blocks.begin[ i ] + 1;
1553       solution := 0;
1554       if n >= 4 then
1555          gblock := [  ];
1556          j1 := 0;
1557          for j in [ blocks.begin[ i ]..blocks.ende[ i ] ] do
1558             j1 := j1 + 1;
1559             gblock[ j1 ] := [  ];
1560             k1 := 0;
1561             for k in [ blocks.begin[ i ]..blocks.ende[ i ] ] do
1562                k1 := k1 + 1;
1563                gblock[ j1 ][ k1 ] := g[ j ][ k ];
1564             od;
1565          od;
1566          root      := [  ];
1567          rootcount := 0;
1568          for j in [1..n] do
1569             nullset( gblock[ j ], j );
1570             if nullcount[ j ] < n - 3 then
1571                rootcount := rootcount + 1;
1572                root[ rootcount ] := j;
1573             fi;
1574          od;
1575          j := 1;
1576          while solution = 0 and j <= rootcount do
1577             ligants := [  ];
1578             ligantscount := 0;
1579             for k in [1..n] do
1580                if k <> root[ j ] and gblock[ root[ j ] ][ k ] <> 0 then
1581                   ligantscount := ligantscount + 1;
1582                   ligants[ ligantscount ] := k;
1583                fi;
1584             od;
1585             begin := 1;
1586             depth := 0;
1587             choice := [  ];
1588             ende := false;
1589             DnSearch( begin, depth, choice );
1590             if solution > 0 then
1591                choice := sol;
1592             fi;
1593             j := j + 1;
1594          od;
1595       fi;
1596
1597    # test of the found irreducibles
1598       if solution = 1 then
1599          # treatment of D4-lattice
1600          found := 0;
1601          preirred := y{ [ blocks.begin[i] .. blocks.ende[i] ] };
1602          z[1] := preirred[choice[1]] + preirred[choice[2]];
1603          if test(z[1]) then
1604             red := ReducedClassFunctions( tbl, preirred, [ z[1] ] );
1605             if ForAll( red.irreducibles, test ) then
1606                found := found + 1;
1607                foundpos := 1;
1608             fi;
1609          fi;
1610          z[2] := preirred[choice[1]] + preirred[choice[3]];
1611          if test(z[2]) then
1612             red := ReducedClassFunctions( tbl, preirred, [ z[2] ] );
1613             if ForAll( red.irreducibles, test ) then
1614                found := found + 1;
1615                foundpos := 2;
1616             fi;
1617          fi;
1618          z[3] := preirred[choice[2]] + preirred[choice[3]];
1619          if test(z[3]) then
1620             red := ReducedClassFunctions( tbl, preirred, [ z[3] ] );
1621             if ForAll( red.irreducibles, test ) then
1622                found := found + 1;
1623                foundpos := 3;
1624             fi;
1625          fi;
1626          if found = 1 then
1627             z := [z[foundpos]];
1628             callreduced();
1629          fi;
1630
1631       else
1632          # treatment of D5-lattice
1633          if solution = 2 then
1634             if choice [ 1 ] <> choice [ 4 ] then
1635                z[ 1 ] := y[ blocks.begin[ i ] + choice[ 1 ] - 1 ];
1636                if choice [ 2 ] <> choice [ 4 ] then
1637                   z[ 1 ]
1638                        := z[ 1 ] + y[ blocks.begin[ i ] + choice[ 2 ] - 1 ];
1639                else
1640                   z[ 1 ]
1641                        := z[ 1 ] + y[ blocks.begin[ i ] + choice[ 3 ] - 1 ];
1642                fi;
1643             else
1644                z[ 1 ] := y[ blocks.begin[ i ] + choice[ 2 ] - 1 ]
1645                        + y[ blocks.begin[ i ] + choice[ 3 ] - 1 ];
1646             fi;
1647             found := 0;
1648             if test(z[1]) then
1649                callreduced();
1650             fi;
1651          else
1652            Append( y2, y{ [ blocks.begin[i] .. blocks.ende[i] ] } );
1653          fi;
1654       fi;
1655    od;
1656
1657    if irredcount > 0 then
1658       g := MatScalarProducts( tbl, y2, y2 );
1659    fi;
1660    else
1661       # input was empty i.e. empty=true
1662       g := [];
1663    fi;
1664    return rec( gram:=g, remainders:=y2, irreducibles:=irred );
1665end );
1666
1667
1668#############################################################################
1669##
1670#F  DnLatticeIterative( <tbl>, <red> )
1671##
1672InstallGlobalFunction( DnLatticeIterative, function( tbl, red )
1673    local dnlat, red1, norms, i, reduc, irred, norm2, g;
1674
1675    # check input parameters
1676    if not IsNearlyCharacterTable( tbl ) then
1677       Error( "first argument must be a nearly character table\n",
1678              "usage: DnLatticeIterative( <tbl>, <record or list> )" );
1679    fi;
1680    if not IsRecord( red ) and not IsList( red ) then
1681       Error( "second argument must be record or list\n",
1682              "usage: DnLatticeIterative( <tbl>, <record or list> )" );
1683    fi;
1684    if IsRecord( red ) and not IsBound( red.remainders ) then
1685       Error( "second record must contain a field 'remainders'\n",
1686              "usage: DnLatticeIterative( <tbl>, <record or list> )" );
1687    fi;
1688    if not IsRecord( red ) then
1689       red := rec( remainders:=red );
1690    fi;
1691    if not IsBound( red.norms ) then
1692       norms := List( red.remainders, x -> ScalarProduct( tbl, x, x ) );
1693    else
1694       norms := ShallowCopy( red.norms );
1695    fi;
1696    reduc := List( red.remainders, ShallowCopy );
1697    irred := [];
1698    repeat
1699       norm2 := [];
1700       for i in [1..Length( reduc )] do
1701          if norms[i] = 2 then
1702             Add( norm2, reduc[i] );
1703          fi;
1704       od;
1705       g := MatScalarProducts( tbl, norm2, norm2 );
1706       dnlat := DnLattice( tbl, g, norm2 );
1707       Append( irred, dnlat.irreducibles );
1708       red1:= ReducedClassFunctions( tbl, dnlat.irreducibles, reduc );
1709       reduc := red1.remainders;
1710       Append( irred, red1.irreducibles );
1711       norms:= List( reduc, x -> ScalarProduct( tbl, x, x ) );
1712    until dnlat.irreducibles=[] and red1.irreducibles=[];
1713    return rec( irreducibles:=irred, remainders:=reduc , norms := norms );
1714end );
1715