home *** CD-ROM | disk | FTP | other *** search
/ Math Solutions 1995 October / Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso / pc / mac / discrete / lib / ctlattic.g < prev    next >
Encoding:
Text File  |  1993-05-05  |  75.1 KB  |  2,665 lines
  1. #############################################################################
  2. ##
  3. #A  ctlattic.g                  GAP library                       Ansgar Kaup
  4. ##
  5. #A  @(#)$Id: ctlattic.g,v 3.10 1993/02/09 14:25:55 martin Rel $
  6. ##
  7. #Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
  8. ##
  9. ##  This file contains those functions which mainly deal with lattices in the
  10. ##  context of character tables.
  11. ##
  12. #H  $Log: ctlattic.g,v $
  13. #H  Revision 3.10  1993/02/09  14:25:55  martin
  14. #H  made undefined globals local
  15. #H
  16. #H  Revision 3.9  1992/12/08  09:08:38  sam
  17. #H  fixed bug in 'OrthogonalEmbeddings'
  18. #H
  19. #H  Revision 3.8  1992/10/20  15:24:46  sam
  20. #H  fixed little bug in 'LLL'
  21. #H
  22. #H  Revision 3.7  1992/09/23  08:40:08  sam
  23. #H  added some comments
  24. #H
  25. #H  Revision 3.6  1992/09/14  13:24:04  sam
  26. #H  fixed bug in 'Decreased'
  27. #H
  28. #H  Revision 3.5  1992/08/07  12:04:51  sam
  29. #H  bug in 'DnLattice' fixed
  30. #H
  31. #H  Revision 3.4  1992/07/10  10:09:36  sam
  32. #H  corrected calls of 'Reduced' in 'DnLattice'
  33. #H
  34. #H  Revision 3.3  1992/07/02  13:01:12  sam
  35. #H  changed test for solution of D4 lattice in 'DnLattice'
  36. #H
  37. #H  Revision 3.2  1992/06/25  10:51:45  sam
  38. #H  fixed bug in 'Decreased':  No return of 'false' if the indicators
  39. #H  cannot be computed, just warning
  40. #H
  41. #H  Revision 3.1  1992/06/12  15:24:56  sam
  42. #H  bug in OESD fixed
  43. #H
  44. #H  Revision 3.0  1992/06/04  07:57:05  sam
  45. #H  initial revision under RCS
  46. #H
  47. ##
  48.  
  49.  
  50. #############################################################################
  51. ##
  52. #F  InfoCharTable1
  53. #F  InfoCharTable2
  54. ##
  55.     if not IsBound( InfoCharTable1 ) then InfoCharTable1:= Ignore; fi;
  56.     if not IsBound( InfoCharTable2 ) then InfoCharTable2:= Ignore; fi;
  57.  
  58.  
  59. #############################################################################
  60. ##
  61. #F  'LLL( <tbl>, <reducibles>, [<y>], [\"sort\"], [\"linearcomb\"] )'
  62. ##
  63. ##
  64. LLL := function( arg )
  65.       local                               
  66. # help fields
  67.             perm, muezw, mmue, Bz, Bzw, sczw, 
  68. # group
  69.             tbl, 
  70. # sensibility    
  71.             y, 
  72. # matrices
  73.             a, b, lcmat, bs, mue, 
  74. # records
  75.             c, 
  76. # vectors
  77.             bsn, bz, B, nullv, scvec, 
  78. # computational fields
  79.             q, ggt, kgv, 
  80. # indices
  81.             i, j, k, ka, l, m, n,
  82. # booleans
  83.             sort, lc,
  84. # sub-procedures
  85.             so, reduce;
  86.  
  87.       reduce := function( k, l ) 
  88.       local r, j;
  89.          if AbsInt( mue[k][l] ) * 2 > 1 then
  90.             r := Int( mue[k][l] );
  91.             if AbsInt( mue[k][l] - r ) * 2 > 1 then
  92.                r := r + SignInt( mue[k][l] );
  93.             fi;
  94.             if r <> 0 then
  95.                if lc then
  96.                   lcmat[k] := lcmat[k] - r * lcmat[l];
  97.                fi;
  98.                b[k] := b[k] -r * b[l];
  99.                for j in [1..l-1] do  
  100.                   if mue[l][j] <> 0 then
  101.                      mue[k][j] := mue[k][j] - r * mue[l][j];
  102.                   fi;
  103.                od;
  104.                mue[k][l] := mue[k][l] - r;
  105.             fi;
  106.          fi;
  107.       end;
  108.  
  109.    # check input parameters
  110.    if not IsRec( arg[1] ) then
  111.       Error( "first argument must be character table\n",
  112.              "usage: LLL( <tbl>, <reducibles>, [<y>], [\"sort\"],",
  113.              " [\"linearcomb\"] )" );
  114.    fi;
  115.    if not( IsList( arg[2] ) and IsList( arg[2][1] ) ) then   
  116.       Error( "second argument must be list of characters\n",
  117.              "usage: LLL( <tbl>, <seq of seq>, [<y>], [\"sort\"],",
  118.              " [\"linearcomb\"] )" );
  119.    fi;  
  120.    y := 3/4;
  121.    if IsBound( arg[3] ) and IsRat( arg[3] ) then
  122.       y := arg[3];
  123.    fi;
  124.    tbl := arg[1];
  125.    a := Filtered( arg[2], x -> ForAny( x, y -> y <> 0 ) );
  126.    sort  := false;
  127.    lc    := false;
  128.    for i in [2..5] do  
  129.       if IsBound( arg[i] ) then
  130.          if arg[i] = "sort" then
  131.             sort := true;
  132.          fi; 
  133.          if arg[i] = "linearcomb" then
  134.             lc   := true;
  135.          fi;
  136.       fi;
  137.    od;
  138.  
  139.       n     := Length( a );
  140.       while a[n] = [] do
  141.          n  := n - 1;
  142.       od;
  143.       if lc then
  144.          lcmat := IdentityMat( n );
  145.       fi;
  146.       if sort then
  147.          perm := SortCharactersCharTable( tbl, a, "degree" );
  148.          if lc then
  149.             lcmat := Permuted( lcmat, perm );
  150.          fi;
  151.       fi;
  152.       m      := Length( a[1] );
  153.       bs     := [];  
  154.       bsn    := [];
  155.       B      := [];
  156.       mue    := [];
  157.       nullv  := [];
  158.       ka     := 1;
  159.       b := [];
  160.       for j in [1..m] do  
  161.          nullv[j] := 0;
  162.       od;
  163.       for i in [1..n] do  
  164.          b[i] := Copy( a[i] );
  165.          mue[i] := [];
  166.          bs[i]   := b[i];
  167.          bsn[i]  := 1;
  168.       od;
  169.  
  170.       # calculate orthogonal base 'bs' (start with 'bs' = 'b' = 'a')
  171.       B[1] := tbl.operations.ScalarProduct( tbl, a[1], a[1] );
  172.       for i in [ 2 .. n ] do
  173.          for j in [ 1 .. i-1 ] do
  174.             Bzw := tbl.operations.ScalarProduct( tbl, a[i], bs[j] );
  175.             if Bzw = 0 then
  176.  
  177.                # 'bs[j]' is already orthogonal to 'a[i]'
  178.                mue[i][j] := 0;
  179.             else
  180.                mue[i][j] := ( Bzw / B[j] ) / bsn[j];
  181.                q := Denominator( mue[i][j] ) * bsn[j];
  182.                ggt := GcdInt( Numerator( mue[i][j] ), bsn[j] );
  183.                kgv    := LcmInt( q / ggt, bsn[i] );
  184.                bs[i]  :=  bs[i] * ( kgv/bsn[i] ) 
  185.                         - bs[j] * ( ( kgv*mue[i][j] )/bsn[j] );
  186.                bsn[i] := kgv;
  187.             fi;
  188.          od;
  189.          B[i] := tbl.operations.ScalarProduct( tbl, bs[i], bs[i] )/
  190.                  ( bsn[i]*bsn[i] );
  191.       od;
  192.       InfoCharTable2( "#I LLL: orthogonal base calculated, ", n,
  193.                       " actual characters\n" );
  194.  
  195.       # calculate vectors
  196.       k := 1;
  197.       InfoCharTable2( "#I LLL: calculating vector 1\n" );
  198.       k := 2;
  199.       repeat 
  200.          if k > ka then
  201.             ka := k;
  202.             InfoCharTable2( "#I LLL: calculating vector", k, "\n" );
  203.          fi;
  204.          reduce( k, k-1 );
  205.          q  := mue[k][k-1] * mue[k][k-1] * B[k-1];
  206.          if y * B[k-1] - q - B[k] > 0 then
  207.             if b[k] = nullv then
  208.  
  209.                # delete!
  210.                if lc then
  211.                   for i in [k..n-1] do
  212.                      lcmat[i] := lcmat[i+1];
  213.                   od;
  214.                fi;
  215.                for i in [k..n-1] do   
  216.                   b[i] := b[i+1];
  217.                   B[i] := B[i+1];
  218.                   for j in [1..k-1] do  
  219.                      mue[i][j] := mue[i+1][j];
  220.                   od;
  221.                   for j in [k..i-1] do  
  222.                      mue[i][j] := mue[i+1][j+1];
  223.                   od;
  224.                od;
  225.                b[n] := [];
  226.                mue[n] := [];
  227.                n := n - 1;
  228.                InfoCharTable2( "#I LLL: ", n, " actual characters\n" );
  229.  
  230.             else
  231.                mmue := mue[k][k-1];
  232.                Bz := B[k] + q;
  233.                if ( B[k] <> 0 and mmue <> 0 ) then
  234.                   mue[k][k-1] := mmue * B[k-1] / Bz;
  235.                   B[k] := B[k] * B[k-1] / Bz;
  236.                   for i in [k+1..n] do  
  237.                      muezw := mue[i][k]; 
  238.                      mue[i][k] := mue[i][k-1] - mmue * mue[i][k];
  239.                      mue[i][k-1] := muezw + mue[k][k-1] * mue[i][k];
  240.                   od;
  241.                else
  242.                   if mmue <> 0 then
  243.                      for i in [k+1..n] do  
  244.                         mue[i][k-1] := mue[i][k-1] / mmue;
  245.                         mue[i][k] := ( mue[i][k-1] - mue[i][k] ) * mmue;
  246.                      od;
  247.                      mue[k][k-1] := 1/mmue;
  248.                   else 
  249.                      for i in [k+1..n] do  
  250.                         muezw       := mue[i][k-1];
  251.                         mue[i][k-1] := mue[i][k];
  252.                         mue[i][k]   := muezw;
  253.                      od;
  254.                      B[k]     := B[k-1];
  255.                   fi;
  256.                fi;
  257.                B[k-1]         := Bz;
  258.                bz := b[k-1];
  259.                b[k-1] := b[k];
  260.                b[k] := bz;
  261.                if lc then
  262.                   bz           := lcmat[k-1];
  263.                   lcmat[k-1]   := lcmat[k];
  264.                   lcmat[k]     := bz;
  265.                fi;
  266.                for j in [1..k-2] do  
  267.                   muezw       := mue[k-1][j];
  268.                   mue[k-1][j] := mue[k][j];
  269.                   mue[k][j]   := muezw;
  270.                od;
  271.                if k > 2 then
  272.                   k := k - 1;
  273.                fi;
  274.             fi;
  275.          else
  276.             for l in [0..k-3] do  
  277.                reduce( k, k-2-l );
  278.             od;
  279.             k := k + 1;
  280.          fi;
  281.       until k > n;
  282.  
  283.       for i in [1..n] do  
  284.          if b[i][1] < 0 then
  285.             b[i] := - b[i];
  286.             if lc then
  287.                lcmat[i] := - lcmat[i];
  288.             fi;
  289.          fi;
  290.       od;
  291.  
  292.    # look for irreducibles
  293.    scvec := [];
  294.    for j in [1..n] do
  295.       scvec[j] := tbl.operations.ScalarProduct( tbl, b[j], b[j] );
  296.    od;
  297.    if ForAny( scvec, x -> x = 1 ) then
  298.      InfoCharTable2( "#I LLL: ", Length( Filtered( scvec, x -> x = 1 ) ),
  299.                      " irreducible characters found\n" );
  300.    fi; 
  301.  
  302.    # prepare the output
  303.    c := rec( irreducibles := [], remainders := [], norms := [] );
  304.    if lc then
  305.       c.irreddecomp:= [];
  306.       c.reddecomp:= [];
  307.       for i in [1..n] do  
  308.          if scvec[i] = 1 then
  309.             Add( c.irreducibles, b[i] );
  310.             Add( c.irreddecomp, lcmat[i] );
  311.          else
  312.             Add( c.remainders, b[i] );
  313.             Add( c.reddecomp, lcmat[i] );
  314.             Add( c.norms, scvec[i] );
  315.          fi; 
  316.       od;
  317.       if sort then
  318.          perm := SortCharactersCharTable( tbl, c.irreducibles, "degree" );
  319.          c.irreddecomp := Permuted( c.irreddecomp, perm );
  320.          perm := SortCharactersCharTable( tbl, c.remainders, "degree" );
  321.          c.reddecomp := Permuted( c.reddecomp, perm );
  322.          c.norms := Permuted( c.norms, perm );
  323.       fi;
  324.    else
  325.       for i in [1..n] do  
  326.          if scvec[i] = 1 then
  327.             Add( c.irreducibles, b[i] );
  328.          else
  329.             Add( c.remainders, b[i] );
  330.             Add( c.norms, scvec[i] );
  331.          fi;
  332.       od;
  333.       if sort then
  334.          SortCharactersCharTable( tbl, c.irreducibles, "degree" );
  335.          perm := SortCharactersCharTable( tbl, c.remainders, "degree" );
  336.          c.norms := Permuted( c.norms, perm );
  337.       fi;
  338.    fi;
  339.    return c;
  340.    end;
  341.  
  342.  
  343. #############################################################################
  344. ##
  345. #F  LLLReducedGramMat( <grammat> )
  346. ##
  347. ##  LLL, working with Gram matrix <grammat>
  348. ##
  349. LLLReducedGramMat := function( grammat )
  350.     local 
  351.     # variables
  352.              y, i, j, k, n, IdMat,
  353.              m, q, bz, l, ba, b, mue, B, v, 
  354.     # procedures
  355.              red, chg, ort, scp, model; 
  356.     
  357.     # shortening of vectors
  358.     red := function( l )
  359.     local i, j, q; 
  360.     if mue[k][l] - 501/1000 > 0 then
  361.        q := Int( mue[k][l] );
  362.        if AbsInt( mue[k][l] - q ) * 2 - 1 > 0 then
  363.           q := q + SignInt( mue[k][l] );
  364.        fi;
  365.        if q <> 0 then
  366.           mue[k] := mue[k] - q * mue[l];
  367.           b[k] := b[k] - q * b[l];
  368.           ba[k] := ba[k] - q * ba[l];
  369.           for i in [1..n] do
  370.              b[i][k] := b[i][k] - q * b[i][l];
  371.           od;
  372.        fi;
  373.     fi;
  374.     end;
  375.     
  376.     # exchange of vectors
  377.     chg := function(  )
  378.     local i, l;
  379.     bz := ba[k];
  380.     ba[k] := ba[k-1];
  381.     ba[k-1] := bz;
  382.     for l in [1..n] do  
  383.        q := b[k][l];
  384.        b[k][l] := b[k-1][l];
  385.        b[k-1][l] := q;
  386.     od;
  387.     for l in [1..n] do  
  388.        q := b[l][k];
  389.        b[l][k] := b[l][k-1];
  390.        b[l][k-1] := q;
  391.     od;
  392.     if k > 2 then
  393.        k := k - 1;
  394.     fi;
  395.     end;
  396.     
  397.     # Norms of b_i^*
  398.     ort := function( i ) 
  399.     local l, k;
  400.     k := b[i][i];
  401.     for l in [1..i-1] do
  402.        k := k - mue[i][l] * mue[i][l] * B[l];
  403.     od;
  404.     B[i] := k;
  405.     end; 
  406.     
  407.     # Scalar Product
  408.     scp := function( i, j )
  409.     local l, k;
  410.     k := b[i][j];
  411.     for l in [1..j-1] do  
  412.        k := k - B[l] * mue[i][l] * mue[j][l];
  413.     od;
  414.     mue[i][j] := k/B[j];
  415.     end; 
  416.     
  417.     # making a model
  418.     model := function(  )
  419.     if k = 2 then
  420.        B[1] := b[1][1];
  421.        i := 2;
  422.     else
  423.        i := k-1;
  424.     fi;
  425.     while i <= k do
  426.        j := 1;
  427.        while j < i do
  428.           scp( i, j );
  429.           j := j + 1;
  430.        od;
  431.        ort( i );
  432.        i := i + 1;
  433.     od;
  434.     end;
  435.     
  436.     # main program
  437.     # check the input
  438.     if not IsList( grammat ) or not IsList( grammat[1] ) then
  439.        Error( "usage: LLLReducedGramMat( <grammat> )" );
  440.     fi;
  441.     n := Length( grammat );
  442.     y := 3/4;
  443.     b := [];
  444.     for i in [1..n] do  
  445.        b[i] := Copy( grammat[i] );
  446.     od;
  447.     B   := [];
  448.     IdMat := IdentityMat( n );
  449.     ba  := Copy( IdMat );
  450.     mue := Copy( IdMat );
  451.     k   := 2;
  452.     while k <= n do
  453.        model(  );
  454.        red( k-1 );
  455.        if B[k] - ( y-mue[k][k-1]*mue[k][k-1] )*B[k-1] < 0 then
  456.           chg(  );
  457.           model(  );
  458.           red( k-1 );
  459.        fi;
  460.        if k > 2 then
  461.           for l in [2..k-1] do  
  462.              red( k-l );
  463.           od;
  464.        fi;
  465.        k := k + 1;
  466.     od;   
  467.     return rec( remainder:= b, transformation:=ba, scalarproducts:=mue,
  468.               bsnorms := B);
  469.     end;  
  470.  
  471.  
  472. #############################################################################
  473. ##
  474. #F  ShortestVectors( <mat>, <bound> [, \"positive\" ] )
  475. ##
  476. ##  kurvec
  477. ##
  478. ShortestVectors := function( arg ) 
  479.     local
  480.     # variables
  481.           n, i, checkpositiv, a, llg, nullv, m, c, q, anz, con, b, v, 
  482.     # procedures
  483.           kur, srt, vschr;
  484.     
  485.     # sub-procedures
  486.     kur := function(  )
  487.     local l;
  488.     for l in [1..n] do  
  489.        v[l] := 0;
  490.     od;
  491.     anz := 0;
  492.     con := true;
  493.     srt( n, 0 );
  494.     end;
  495.     
  496.     # search for shortest vectors
  497.     srt := function( d, dam )
  498.     local i, j, x, k, k1, q;
  499.     if d = 0 then
  500.        if v = nullv then
  501.           con := false;
  502.        else
  503.           anz := anz + 1;
  504.           vschr( dam );
  505.        fi;
  506.     else
  507.        x := 0;
  508.        for j in [d+1..n] do  
  509.           x := x + v[j] * llg.scalarproducts[j][d];
  510.        od;
  511.        i := - Int( x );
  512.        if AbsInt( -x-i ) * 2 - 1 > 0 then             
  513.           i := i - SignInt( x );
  514.        fi;
  515.        k := i + x;
  516.        q := ( m + 1/1000 - dam ) / llg.bsnorms[d];
  517.        if k * k - q < 0 then
  518.           repeat
  519.              i := i + 1;
  520.              k := k + 1;
  521.           until k * k - q > 0 and k > 0;
  522.           i := i - 1;
  523.           k := k - 1;
  524.           while k * k - q  < 0 and con do
  525.              v[d] := i;
  526.              k1 := llg.bsnorms[d] * k * k + dam;
  527.              srt( d-1, k1 );
  528.              i := i - 1;
  529.              k := k - 1;
  530.           od;
  531.        fi;
  532.     fi;
  533.     end;
  534.     
  535.     # output of vector
  536.     vschr := function( dam )
  537.     local i, j, w, neg;
  538.     c.vectors[anz] := [];
  539.     neg := false;
  540.     for i in [1..n] do  
  541.        w := 0;
  542.        for j in [1..n] do   
  543.           w := w + v[j] * llg.transformation[j][i];
  544.        od;
  545.        if w < 0 then 
  546.           neg := true;
  547.        fi;
  548.        c.vectors[anz][i] := w;
  549.     od;
  550.     if checkpositiv and neg then
  551.        c.vectors[anz] := [];
  552.        anz := anz - 1;
  553.     else
  554.        c.norms[anz] := dam;
  555.     fi;
  556.     end;
  557.     
  558.     # main program
  559.     # check input
  560.     if not IsBound( arg[1] ) 
  561.     or not IsList( arg[1] ) or not IsList( arg[1][1] ) then
  562.        Error ( "first argument must be Gram matrix\n",
  563.           "usage: ShortestVectors( <mat>, <integer> [,<\"positive\">] )" );
  564.     fi;
  565.     if not IsBound( arg[2] ) or not IsInt( arg[2] ) then
  566.        Error ( "second argument must be integer\n",
  567.           "usage: ShortestVectors( <mat>, <integer> [,<\"positive\">] )");
  568.     fi;
  569.     if IsBound( arg[3] ) then
  570.        if IsString( arg[3] ) then
  571.           if arg[3] = "positive" then
  572.              checkpositiv := true;
  573.           else
  574.              checkpositiv := false;
  575.           fi;
  576.        else
  577.           Error ( "third argument must be string\n",
  578.           "usage: ShortestVectors( <mat>, <integer> [,<\"positive\">] )");
  579.        fi;
  580.     else
  581.        checkpositiv := false;
  582.     fi;
  583.     a := arg[1];
  584.     m := arg[2];
  585.     n := Length( a );
  586.     b := [];
  587.     for i in [1..n] do  
  588.        b[i] := Copy( a[i] );
  589.     od;
  590.     c    := rec( vectors:=[],norms:=[]);
  591.     v    := [];
  592.     nullv := [];
  593.     for i in [1..n] do  
  594.        nullv[i] := 0;
  595.     od;
  596.     llg:=LLLReducedGramMat(b);
  597.     kur();
  598.     InfoCharTable2( "#I ShortestVectors: ",
  599.                     Length( c.vectors ), " vectors found\n" );
  600.     return( c );
  601.     end;
  602.  
  603.  
  604. #############################################################################
  605. ##
  606. #F  Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing> ] )
  607. ##
  608. Extract := function( arg )
  609.     local  
  610.     
  611.     # indices
  612.           i, j, k, l, n,
  613.     # input arrays
  614.           tbl, y, gram, missing,
  615.     # booleans
  616.           deeper, iszero, used, nullbegin, nonmissing,
  617.           maxnorm, minnorm, normbound, maxsum, solmat, 
  618.           f, squares, sfind, choicecollect, sequence,
  619.           dependies, solcollect, sum, solcount, max, sumac, kmax, 
  620.           solution, 
  621.     # functions
  622.           next, zeroset, possiblies, update, correctnorm,  
  623.           maxsquare, square, ident, begin;
  624.     
  625.     # choosing next vector for combination
  626.     next := function( lines, solumat, acidx )
  627.     local i, j, solmat, testvec, idxback;
  628.     
  629.     while acidx <= n and k + n - acidx >= kmax do
  630.        solmat := Copy( solumat );
  631.        if k = 0 then
  632.           i := acidx;
  633.           while i <= n and not begin( sequence[i] ) do
  634.              i := i + 1;
  635.           od;
  636.           if i > n then
  637.              nullbegin := true;
  638.           else
  639.              nullbegin := false;
  640.              if i > acidx then
  641.                 idxback := sequence[i];
  642.                 for j in [acidx + 1..1] do
  643.                    sequence[j] := sequence[j -1];
  644.                 od;
  645.                 sequence[acidx] := idxback;
  646.              fi;
  647.           fi;
  648.        fi;
  649.        k := k + 1;
  650.        f[k] := sequence[acidx];
  651.        testvec := [];
  652.        for i in [1..k] do 
  653.           testvec[i] := gram[f[k]][f[i]];
  654.        od;
  655.        zeroset( solmat, testvec, lines );
  656.        acidx := acidx + 1; 
  657.        possiblies( 1, solmat, testvec, acidx, lines );
  658.        k := k - 1;
  659.     od;
  660.     end;
  661.     
  662.     # filling zero in places that fill already the conditions
  663.     zeroset := function( solmat, testvec, lines )
  664.     local i, j;
  665.     
  666.     for i in [1..k-1] do  
  667.        if testvec[i] = 0 then
  668.           for j in [1..lines] do  
  669.              if solmat[j][i] <> 0 and not IsBound( solmat[j][k] ) then 
  670.                 solmat[j][k] := 0;
  671.              fi;
  672.           od;
  673.        fi;
  674.     od;
  675.     end;
  676.     
  677.     # try and error for the chosen vector
  678.     possiblies := function( start, solmat, testvect, acidx, lines )
  679.     local i, j, remainder, toogreat, equal, solmatback, testvec;
  680.     
  681.     testvec := Copy( testvect );
  682.     toogreat := false;
  683.     equal := true; 
  684.     if k > 1 then
  685.        for i in [1..k-1] do    
  686.           if testvec[i] < 0 then
  687.              toogreat := true;
  688.           fi;
  689.           if testvec[i] <> 0 then
  690.              equal := false;
  691.           fi;
  692.        od;
  693.        if testvec[k] < 0 then
  694.           toogreat := true;
  695.        fi;
  696.     else
  697.        if not nullbegin then
  698.           while start <= gram[f[k]][f[k]] and start < missing do
  699.              solmat[start][k] := 1;
  700.              start := start + 1;
  701.           od;
  702.           testvec[k] := 0;
  703.           if gram[f[k]][f[k]] > lines then
  704.              lines := gram[f[k]][f[k]];
  705.           fi;
  706.        else
  707.           lines := 0;
  708.        fi;
  709.     fi;
  710.     if not equal and not toogreat then
  711.        while start < lines and IsBound( solmat[start][k] ) do
  712.           start := start + 1;
  713.        od;
  714.        if start <= lines and not IsBound( solmat[start][k] ) then
  715.           solmat[start][k] := 0;
  716.           while not toogreat and not equal do
  717.              solmat[start][k] := solmat[start][k] + 1;
  718.              testvec := update( -1, testvec, start, solmat );
  719.              equal := true;
  720.              for i in [1..k-1] do    
  721.                 if testvec[i] < 0 then
  722.                    toogreat := true;
  723.                 fi;
  724.                 if testvec[i] <> 0 then
  725.                    equal := false;
  726.                 fi;
  727.              od;
  728.              if testvec[k] < 0 then
  729.                 toogreat := true;
  730.              fi;
  731.           od;
  732.        fi;
  733.     fi;
  734.     if equal and not toogreat then 
  735.        solmatback := Copy( solmat );
  736.        for i in [1..missing] do          
  737.           if not IsBound( solmat[i][k] ) then        
  738.              solmat[i][k] := 0;
  739.           fi;
  740.        od;
  741.        correctnorm( testvec[k], solmat, lines + 1, testvec[k], acidx, lines );
  742.        solmat := Copy( solmatback );
  743.     fi;
  744.     if k > 1 then
  745.        while start <= lines and solmat[start][k] > 0 do
  746.           solmat[start][k] := solmat[start][k] - 1;
  747.           testvec := update( 1, testvec, start, solmat );
  748.           solmatback := Copy( solmat );
  749.           zeroset( solmat, testvec, lines );
  750.           deeper := false;
  751.           for i in [1..k-1] do
  752.              if solmat[start][i] <> 0 then
  753.                 deeper := false;
  754.                 if testvec[i] = 0 then
  755.                    deeper := true;
  756.                 else 
  757.                    for j in [1..missing] do
  758.                       if solmat[j][i] <> 0 and not IsBound(solmat[j][k]) then 
  759.                         deeper := true;
  760.                       fi;
  761.                    od;
  762.                 fi;
  763.              fi;
  764.           od;
  765.           if deeper then
  766.              possiblies( start + 1, solmat, testvec, acidx, lines );
  767.           fi;
  768.           solmat := Copy( solmatback );
  769.        od;
  770.     fi;
  771.     end;
  772.     
  773.     # update the remaining conditions to fill
  774.     update := function( x, testvec, start, solmat )
  775.     local i;
  776.     for i in [1..k-1] do    
  777.        if solmat[start][i] <> 0 then
  778.           testvec[i] := testvec[i] + solmat[start][i] * x;
  779.        fi;
  780.     od;
  781.     testvec[k] := testvec[k] - square( solmat[start][k] ) 
  782.                              + square( solmat[start][k] + x );
  783.     return( testvec );
  784.     end;
  785.     
  786.     # correct the norm if all other conditions are filled
  787.     correctnorm := function( remainder, solmat, pos, max, acidx, lines )           
  788.     local i, r, newsol, ret;
  789.     if remainder = 0 and pos <= missing + 1 then
  790.        newsol := true;
  791.        for i in [1..solcount[k]] do
  792.           if ident( solcollect[k][i], solmat ) = missing then
  793.              newsol := false;
  794.           fi;
  795.        od;
  796.        if newsol then
  797.           if k > kmax then
  798.              kmax := k;
  799.           fi;
  800.           solcount[k] := solcount[k] + 1;
  801.           solcollect[k][solcount[k]] := [];
  802.           choicecollect[k][solcount[k]] := Copy( f );
  803.           for i in [1..Length( solmat )] do
  804.              solcollect[k][solcount[k]][i] := Copy( solmat[i] );
  805.           od;
  806.           if k = n and pos = missing + 1 then
  807.              ret := 0;
  808.           else
  809.              ret := max;
  810.              if k <> n then
  811.                 next( lines, solmat, acidx );
  812.              fi;
  813.           fi;
  814.        else 
  815.           ret := max;
  816.        fi;
  817.     else
  818.        if pos <= missing then
  819.           i := maxsquare( remainder, max );
  820.           while i > 0 do
  821.              solmat[pos][k] := i;
  822.              i := correctnorm( remainder-square( i ), 
  823.                                solmat, pos+1, i, acidx, lines + 1);
  824.              i := i - 1;
  825.           od;
  826.           if i < 0 then
  827.              ret := 0;
  828.           else
  829.              ret := max;
  830.           fi;
  831.        else
  832.           ret := 0;
  833.        fi;
  834.     fi;
  835.     return( ret );
  836.     end;
  837.     
  838.     # compute the maximum squarenumber lower then given integer
  839.     maxsquare := function( value, max ) 
  840.     local i;
  841.     
  842.     i := 1;
  843.     while square( i ) <= value and i <= max do
  844.           i := i + 1;
  845.     od;
  846.     return( i-1 );
  847.     end;
  848.     
  849.     square := function( i ) 
  850.     if i = 0 then
  851.        return( 0 );
  852.     else
  853.        if not IsBound( squares[i] ) then
  854.           squares[i] := i * i;
  855.        fi;
  856.        return( squares[i] );
  857.     fi;
  858.     end;
  859.     
  860.     ident := function( a, b ) 
  861.     # lists the identities of the two given sequences and counts them
  862.     local i, j, k, zi, zz, la, lb;
  863.     la := Length( a );
  864.     lb := Length( b );
  865.     zi := [];
  866.     zz := 0;
  867.     for i in [1..la] do  
  868.        j := 1;
  869.        repeat 
  870.           if a[i] = b[j] then
  871.              k :=1;
  872.              while k <= zz and j <> zi[k] do
  873.                 k := k + 1;
  874.              od;
  875.              if k > zz then
  876.                 zz := k;
  877.                 zi[zz] := j;
  878.                 j := lb;
  879.              fi;
  880.           fi;
  881.           j := j + 1;
  882.        until j > lb;
  883.     od;
  884.     return( zz );
  885.     end;
  886.     
  887.     # looking for character that can stand at the beginning
  888.     begin := function( i )
  889.     local ind;
  890.     if y = [] or gram[i][i] < 4 then
  891.        return true;
  892.     else
  893.        if IsBound( tbl.powermap ) and IsBound( tbl.powermap[2] ) then
  894.           if IsList( tbl.powermap[2] ) and ForAll( tbl.powermap[2], IsInt ) then
  895.              ind := AbsInt( Indicator( tbl, [y[i]], 2 )[1]);
  896.              if gram[i][i] - 1 <= ind 
  897.              or ( gram[i][i] = 4 and ind = 1 ) then
  898.                 return true;
  899.              fi;
  900.           fi;
  901.        fi;
  902.     fi;
  903.     return false;
  904.     end;
  905.     
  906.     # check input parameters
  907.     if IsCharTable( arg[1] ) then
  908.        tbl := arg[1];
  909.     else
  910.        Error( "first argument must be character-table\n \
  911.     usage: Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing>] )" );
  912.     fi;
  913.     if IsBound( arg[2] ) and IsList( arg[2] ) and IsList( arg[2][1] ) then
  914.        y := Copy( arg[2] );
  915.     else
  916.        Error( "second argument must be list of reducible characters\n \
  917.     usage: Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing>] )" );
  918.     fi;
  919.     if IsBound( arg[2] ) and IsList( arg[3] ) and IsList( arg[3][1] ) then
  920.        gram := Copy( arg[3] );
  921.     else
  922.        Error( "third argument must be gram-matrix of reducible characters\n \
  923.     usage: Extract( <tbl>, <reducibles>, <gram-matrix> [, <missing>] )" );
  924.     fi;
  925.     n := Length( gram );
  926.     if IsBound( arg[4] ) and IsInt( arg[4] ) then
  927.        missing := arg[4];
  928.     else
  929.        missing := n;
  930.        nonmissing := true;
  931.     fi;
  932.     
  933.     # main program
  934.     maxnorm := 0;
  935.     minnorm := gram[1][1];
  936.     normbound := [];
  937.     maxsum := [];
  938.     solcollect := [];
  939.     choicecollect := [];
  940.     sum := [];
  941.     solmat := [];
  942.     used := [];
  943.     solcount := [];
  944.     sfind := [];
  945.     f := [];
  946.     squares := [];
  947.     kmax := 0;
  948.     for i in [1..missing] do
  949.        solmat[i] := [];
  950.     od;
  951.     for i in [1..n] do  
  952.        solcount[i] := 0;
  953.        used[i] := false;
  954.        solcollect[i] := [];
  955.        choicecollect[i] := [];
  956.     od;
  957.     for i in [1..n] do  
  958.        if gram[i][i] > maxnorm then
  959.           maxnorm := gram[i][i];
  960.        else
  961.           if gram[i][i] < minnorm then
  962.              minnorm := gram[i][i];
  963.           fi;
  964.        fi; 
  965.     od;
  966.     j := 0;
  967.     for i in [minnorm..maxnorm] do
  968.        k := 1;
  969.        while k <= n and gram[k][k] <> i do
  970.           k := k + 1;
  971.        od;
  972.        if k <= n then
  973.            j := j + 1;
  974.            normbound[j] := rec( norm:=i, first:=k, last:=0 );
  975.            if k = n then
  976.               normbound[j].last := k;
  977.            else
  978.               k := n;
  979.               while gram[k][k] <> i and k > 0 do
  980.                  k := k - 1;
  981.               od;
  982.               if k > 0 then
  983.                  normbound[j].last := k;
  984.               fi;
  985.            fi;
  986.        fi;
  987.     od;
  988.     for j in [1..Length( normbound )] do
  989.        maxsum[j] := 0;
  990.        for i in [normbound[j].first..normbound[j].last] do
  991.           if gram[i][i] = normbound[j].norm then
  992.              sum[i] := 0;
  993.              for k in [1..n] do  
  994.                 sum[i] := sum[i] + gram[i][k];
  995.              od;
  996.              if sum[i] > maxsum[j] then
  997.                 maxsum[j] := sum[i];
  998.              fi;
  999.           fi;
  1000.        od;
  1001.     od; 
  1002.     k := 1;
  1003.     sequence := [];
  1004.     i:= 1;
  1005.     while i <= Length( normbound ) do
  1006.        max := maxsum[i];
  1007.        sumac := 0;
  1008.        for j in [normbound[i].first..normbound[i].last] do
  1009.           if gram[j][j] = normbound[i].norm and sum[j] > sumac 
  1010.           and sum[j] <= max and not used[j] then
  1011.              sequence[k] := j;
  1012.              sumac := sum[j];
  1013.           fi;
  1014.        od;
  1015.        if IsBound( sequence[k] ) then
  1016.           max := sumac;
  1017.           used[sequence[k]] := true;
  1018.           k := k + 1;
  1019.        else
  1020.           i := i + 1;
  1021.        fi;
  1022.     od;
  1023.     k := 0;
  1024.     next( 1, solmat, 1 );
  1025.     solution := rec( solution := [], choice := choicecollect[kmax] );
  1026.     for i in [1..solcount[kmax]] do
  1027.        solution.solution[i] := [];
  1028.        l := 0;
  1029.        for j in [1..missing] do
  1030.           iszero := true;
  1031.           for k in [1..kmax] do
  1032.              if solcollect[kmax][i][j][k] <> 0 then
  1033.                 iszero := false;
  1034.              fi;
  1035.           od;
  1036.           if not iszero then
  1037.              l := l + 1;
  1038.              solution.solution[i][l] := solcollect[kmax][i][j];  
  1039.           fi;
  1040.        od;
  1041.     od;
  1042.     return( solution );
  1043.     end;
  1044.  
  1045.  
  1046. #############################################################################
  1047. ##
  1048. #F  Decreased( <tbl>, <chars>, <decompmat>, [ <choice> ] )
  1049. ##
  1050. Decreased := function( arg ) 
  1051.     local 
  1052.         # indices
  1053.           m, n, m1, n1, i, i1, i2, i3, i4, j, jj, j1, j2, j3, 
  1054.         # booleans
  1055.           ende1, ende2, ok, change, delline, delcolumn,
  1056.         # help fields
  1057.           deleted, kgv, l1, l2, l3, dim, ident, 
  1058.         # matrices
  1059.           invmat, remmat, remmat2, solmat, nonzero, 
  1060.         # double-indices
  1061.           columnidx, lineidx, system, components, compo2, 
  1062.         # output-fields
  1063.           sol, red, redcount, irred,
  1064.         # help fields
  1065.           IRS, SFI, lc, nc, char, char1, entries, 
  1066.         # input fields
  1067.           tbl, y, solmat, choice, 
  1068.         # functions
  1069.           Idxset, Identset, Invadd, Invmult, Nonzeroset;
  1070.     
  1071.     Idxset := function()
  1072.     # update indices
  1073.     local i1, j1;
  1074.     i1 := 0;
  1075.     for i in [1..m] do
  1076.        if not delline[i] then
  1077.           i1 := i1 + 1;
  1078.           lineidx[i1] := i;
  1079.        fi;
  1080.     od;
  1081.     m1 := i1;
  1082.     j1 := 0;
  1083.     for j in [1..n] do
  1084.        if not delcolumn[j] then
  1085.           j1 := j1 + 1;
  1086.           columnidx[j1] := j;
  1087.        fi;
  1088.     od;
  1089.     n1 := j1;
  1090.     end;
  1091.     
  1092.     Identset := function( veca, vecb )
  1093.     # count identities of veca and vecb and store "non-identities"
  1094.     local la, lb, i, j, n, nonid, nic, r;
  1095.     n := 0;
  1096.     la := Length( veca );
  1097.     lb := Length( vecb );
  1098.     j := 1;
  1099.     nonid := [];
  1100.     nic := 0;
  1101.     for i in [1..la] do  
  1102.        while j <= lb and veca[i] > vecb[j] do
  1103.           nic := nic + 1;
  1104.           nonid[nic] := vecb[j];
  1105.           j := j + 1;
  1106.        od;
  1107.        if j <= lb and veca[i] = vecb[j] then
  1108.           n := n + 1;
  1109.           j := j + 1;
  1110.        fi;
  1111.     od;
  1112.     while j <= lb do
  1113.        nic := nic + 1;
  1114.        nonid[nic] := vecb[j];
  1115.        j := j + 1;
  1116.     od;
  1117.     r := rec( nonid := nonid, id := n  );
  1118.     return( r );
  1119.     end;
  1120.     
  1121.     Invadd := function( j1, j2, l )
  1122.     # addition of two lines of invmat
  1123.     local i;
  1124.     for i in [1..n] do  
  1125.        if invmat[i][j2] <> 0 then
  1126.           invmat[i][j1] := invmat[i][j1] - l * invmat[i][j2];
  1127.        fi;
  1128.     od;
  1129.     end;
  1130.     
  1131.     Invmult := function( j1, l )
  1132.     # multiply line of invmat
  1133.     local i;
  1134.     if l <> 1 then
  1135.        for i in [1..n] do  
  1136.           if invmat[i][j1] <> 0 then
  1137.              invmat[i][j1] := invmat[i][j1] * l;
  1138.           fi;
  1139.        od;
  1140.     fi;
  1141.     end;
  1142.     
  1143.     Nonzeroset := function( j ) 
  1144.     # entries <> 0 in j-th column
  1145.     
  1146.     local i, j1;
  1147.     nonzero[j] := [];
  1148.     j1 := 0;
  1149.     for i in [1..m] do  
  1150.        if solmat[i][j] <> 0 then
  1151.           j1 := j1 + 1;
  1152.           nonzero[j][j1] := i;
  1153.        fi;
  1154.     od;
  1155.     entries[j] := j1;
  1156.     end;
  1157.     
  1158.     # check input parameters
  1159.     if IsCharTable( arg[1] ) then
  1160.        tbl := arg[1];
  1161.     else
  1162.        Error( "first argument must be character-table\n",
  1163.               "usage: Decreased( <tbl>, <list of char>,\n", 
  1164.               "<decomposition-matrix>, [<choice>] )" );
  1165.     fi;
  1166.     if IsList( arg[2] ) and IsList( arg[2][1] ) then
  1167.        y := arg[2];
  1168.     else
  1169.        Error( "second argument must be list of characters\n",
  1170.               "usage: Decreased( <tbl>, <list of char>,\n", 
  1171.               "<decomposition-matrix>, [<choice>] )" );
  1172.     fi;
  1173.     if IsList( arg[3] ) and IsList( arg[3][1] ) then
  1174.        solmat := Copy( arg[3] );
  1175.     else
  1176.        Error( "third argument must be decomposition matrix\n",
  1177.               "usage: Decreased( <tbl>, <list of char>,\n", 
  1178.               "<decomposition-matrix>, [<choice>] )" );
  1179.     fi;
  1180.     if not IsBound( arg[4] ) then
  1181.        choice := [];
  1182.        for i in [1..Length( y )] do
  1183.           choice[i] := i;
  1184.        od;
  1185.     else
  1186.        if IsList( arg[4] ) then
  1187.           choice := arg[4];
  1188.        else
  1189.           Error( "forth argument contains choice of characters\n",
  1190.               "usage: Decreased( <tbl>, <list of char>,\n", 
  1191.               "<decomposition-matrix>, [<choice>] )" );
  1192.        fi;
  1193.     fi;
  1194.      
  1195.     # initialisations
  1196.     lc := Length( y[1] );
  1197.     nc := [];
  1198.     for i in [1..lc] do  
  1199.        nc[i] := 0;
  1200.     od;
  1201.     columnidx := [];
  1202.     lineidx   := [];
  1203.     nonzero   := [];
  1204.     entries   := [];
  1205.     delline   := [];
  1206.     delcolumn := [];
  1207.  
  1208.     # number of lines
  1209.     m := Length( solmat );
  1210.  
  1211.     # number of columns
  1212.     n := Length( solmat[1] );
  1213.     invmat := [];
  1214.     for i in [1..n] do
  1215.        invmat[i] := [];
  1216.        for j in [1..n] do
  1217.           invmat[i][j] := 0;
  1218.        od;
  1219.     od;
  1220.     invmat := invmat^0;
  1221.     for i in [1..m] do  
  1222.        delline[i] := false;
  1223.     od;
  1224.     for j in [1..n] do   
  1225.        delcolumn[j] := false;
  1226.     od;
  1227.     i := 1;
  1228.  
  1229.     # check lines for information
  1230.     while i <= m do
  1231.        if not delline[i] then
  1232.           entries[i] := 0;
  1233.           for j in [1..n] do  
  1234.              if solmat[i][j] <> 0 and not delcolumn[j] then
  1235.                 entries[i] := entries[i] + 1;
  1236.                 if entries[i] = 1 then
  1237.                    nonzero[i] := j;
  1238.                 fi;
  1239.              fi;
  1240.           od;
  1241.           if entries[i] = 1 then
  1242.              delcolumn[nonzero[i]] := true;
  1243.              delline[i] := true;
  1244.              j := 1; 
  1245.              while j < i and solmat[j][nonzero[i]] = 0 do
  1246.                 j := j + 1;
  1247.              od;
  1248.              if j < i then
  1249.                 i := j;
  1250.              else
  1251.                 i := i + 1;
  1252.              fi;
  1253.           else
  1254.              if entries[i] = 0 then
  1255.                 delline[i] := true;
  1256.              fi;
  1257.              i := i + 1;
  1258.           fi;
  1259.        else
  1260.           i := i + 1;
  1261.        fi;
  1262.     od;
  1263.     Idxset();
  1264.     
  1265.     deleted := m - Length(lineidx);
  1266.     for j in [1..n] do  
  1267.        Nonzeroset( j );
  1268.     od;
  1269.     ende1 := false;
  1270.     while not ende1 and deleted < m do
  1271.        j := 1;
  1272.  
  1273.     # check solo-entry-columns
  1274.        while j <= n do
  1275.           if entries[j] = 1 then
  1276.              change := false;
  1277.              for jj in [1..n] do  
  1278.                 if (delcolumn[j] and delcolumn[jj])
  1279.                 or not delcolumn[j] then
  1280.                    if solmat[nonzero[j][1]][jj] <> 0 and jj <> j then 
  1281.                       change := true;
  1282.                       kgv := Lcm( solmat[nonzero[j][1]][j], 
  1283.                               solmat[nonzero[j][1]][jj] );
  1284.                       l1 := kgv / solmat[nonzero[j][1]][jj];
  1285.                       Invmult( jj, l1 );
  1286.                       for i1 in [1..Length( nonzero[jj] )] do
  1287.                          solmat[nonzero[jj][i1]][jj] 
  1288.                                := solmat[nonzero[jj][i1]][jj] * l1;
  1289.                       od;
  1290.                       Invadd( jj, j, kgv/solmat[nonzero[j][1]][j] );
  1291.                       solmat[nonzero[j][1]][jj] := 0;
  1292.                       Nonzeroset( jj );
  1293.                    fi;
  1294.                 fi;
  1295.              od;
  1296.              if not delline[nonzero[j][1]] then
  1297.                 delline[nonzero[j][1]] := true;
  1298.                 delcolumn[j] := true;
  1299.                 deleted := deleted + 1;
  1300.                 Idxset();
  1301.              fi;
  1302.              if change then
  1303.                 j := 1;
  1304.              else
  1305.                 j := j + 1;
  1306.              fi;
  1307.           else
  1308.              j := j + 1;
  1309.           fi;
  1310.        od;
  1311.  
  1312.     # search for Equality-System
  1313.     # system :     chosen columns
  1314.     # components : entries <> 0 in the chosen columns
  1315.        dim := 2;
  1316.        change := false;
  1317.        ende2 := false;
  1318.        while dim <= n1 and not ende2 do
  1319.           j3 := 1;
  1320.           while j3 <= n1 and not ende2 do
  1321.              j2 := j3;
  1322.              j1 := 0;
  1323.              system := [];
  1324.              components := [];
  1325.              while j2 <= n1 do
  1326.                 while j2 <= n1 and entries[columnidx[j2]] > dim do
  1327.                    j2 := j2 + 1;
  1328.                 od;
  1329.                 if j2 <= n1 then
  1330.                    if j1 = 0 then
  1331.                       j1 := 1;
  1332.                       system[j1] := columnidx[j2];
  1333.                       components := Copy( nonzero[columnidx[j2]] );
  1334.                    else
  1335.                       ident := Identset( components, nonzero[columnidx[j2]] );
  1336.                       if dim - Length( components ) >= entries[columnidx[j2]]
  1337.                              - ident.id then
  1338.                          j1 := j1 + 1;
  1339.                          system[j1] := columnidx[j2];
  1340.                          if ident.id < entries[columnidx[j2]] then
  1341.                             compo2 := Copy( components );
  1342.                             components := [];
  1343.                             i1 := 1;
  1344.                             i2 := 1;
  1345.                             i3 := 1;
  1346.  
  1347.     # append new entries to "components"
  1348.                             while i1 <= Length( ident.nonid )
  1349.                                or i2 <= Length( compo2 ) do
  1350.                                if i1 <= Length( ident.nonid ) then
  1351.                                   if i2 <= Length( compo2 ) then
  1352.                                      if ident.nonid[i1] < compo2[i2] then
  1353.                                         components[i3] := ident.nonid[i1];
  1354.                                         i1 := i1 + 1;
  1355.                                      else
  1356.                                         components[i3] := compo2[i2];
  1357.                                         i2 := i2 + 1;
  1358.                                      fi;
  1359.                                   else
  1360.                                      components[i3] := ident.nonid[i1];
  1361.                                      i1 := i1 + 1;
  1362.                                   fi;
  1363.                                else
  1364.                                   if i2 <= Length( compo2 ) then
  1365.                                      components[i3] := compo2[i2];
  1366.                                      i2 := i2 + 1;
  1367.                                   fi;
  1368.                                fi;
  1369.                                i3 := i3 + 1;
  1370.                             od;
  1371.                          fi;
  1372.                       fi;
  1373.                    fi;
  1374.                    j2 := j2 + 1;
  1375.                 fi;
  1376.              od;
  1377.  
  1378.     # try to solve system with Gauss
  1379.              if  Length( system ) > 1 then
  1380.                 for i1 in [1..Length( components )] do
  1381.                    i2 := 1;
  1382.                    repeat
  1383.                       ok := true;
  1384.                       if solmat[components[i1]][system[i2]] = 0 then
  1385.                          ok := false;
  1386.                       else
  1387.                          for i3 in [1..i1-1] do    
  1388.                             if solmat[components[i3]][system[i2]] <> 0 then
  1389.                                ok := false;
  1390.                             fi;
  1391.                          od;
  1392.                       fi;
  1393.                       if not ok then
  1394.                          i2 := i2 + 1;
  1395.                       fi;
  1396.                    until ok or i2 > Length( system );
  1397.                    if ok then
  1398.                       for i3 in [1..Length( system )] do
  1399.                          if i3 <> i2 
  1400.                             and solmat[components[i1]][system[i3]] <> 0 then
  1401.                             change := true;
  1402.                             kgv := Lcm( solmat[components[i1]][system[i3]], 
  1403.                                        solmat[components[i1]][system[i2]] );
  1404.                             l2 := kgv / solmat[components[i1]][system[i2]];
  1405.                             l3 := kgv / solmat[components[i1]][system[i3]];
  1406.                             for i4 in [1..Length( nonzero[system[i3]] )] do
  1407.                                solmat[nonzero[system[i3]][i4]][system[i3]]
  1408.                             := solmat[nonzero[system[i3]][i4]][system[i3]]*l3;
  1409.                             od;
  1410.                             Invmult( system[i3], l3 );
  1411.                             for i4 in [1..Length( nonzero[system[i2]] )] do
  1412.                                solmat[nonzero[system[i2]][i4]][system[i3]]
  1413.                                := solmat[nonzero[system[i2]][i4]][system[i3]]
  1414.                             -  solmat[nonzero[system[i2]][i4]][system[i2]]*l2;
  1415.                             od;
  1416.                             Invadd( system[i3], system[i2], l2 );
  1417.                             Nonzeroset( system[i3] );
  1418.                             if entries[system[i3]] = 0 then
  1419.                                delcolumn[system[i3]] := true;
  1420.                                Idxset();
  1421.                             fi;
  1422.                          fi;
  1423.                       od;
  1424.                    fi;
  1425.                 od;
  1426.  
  1427.    # check for columns with only one entry <> 0
  1428.                 for i1 in [1..Length( system )] do
  1429.                    if entries[system[i1]] = 1 then
  1430.                       ende2 := true;
  1431.                    fi;
  1432.                 od;
  1433.                 if not ende2 then
  1434.                    j3 := j3 + 1;
  1435.                 fi;
  1436.              else
  1437.                 j3 := j3 + 1;
  1438.              fi;
  1439.           od;
  1440.           dim := dim + 1;
  1441.        od;
  1442.        if dim > n1 and not change and j3 > n1 then
  1443.           ende1 := true;
  1444.        fi;
  1445.     od;
  1446.     
  1447.     # check, if
  1448.     #    the transformation of solmat allows computation of new irreducibles
  1449.     remmat := [];
  1450.     for i in [1..m] do  
  1451.        remmat[i] := [];
  1452.        delline[i] := true;
  1453.     od;
  1454.     redcount := 0;
  1455.     red := [];
  1456.     irred := [];
  1457.     j := 1;
  1458.     sol := true;
  1459.     while j <= n and sol do
  1460.  
  1461.     # computation of character
  1462.        char := Copy( nc );
  1463.        for i in [1..n] do  
  1464.           if invmat[i][j] <> 0 then
  1465.              char := char + invmat[i][j] * y[choice[i]];
  1466.           fi;
  1467.        od;
  1468.  
  1469.     # probably irreducible ==> has to pass tests
  1470.        if entries[j] = 1 then
  1471.           if solmat[nonzero[j][1]][j] <> 1 then
  1472.              char1 := char/solmat[nonzero[j][1]][j];
  1473.           else
  1474.              char1 := char;
  1475.           fi;
  1476.           if char1[1] < 0 then
  1477.              char1 := - char1;
  1478.           fi;
  1479.  
  1480.     # is 'char1' real?
  1481.           IRS := ForAll( char1, x -> GaloisCyc(x,-1) = x );
  1482.  
  1483.    # Frobenius Schur indicator
  1484.           if IsBound( tbl.powermap ) and IsBound( tbl.powermap[2] ) then
  1485.              SFI:= Indicator( tbl, [ char1 ], 2 )[1];
  1486.           else
  1487.              SFI:= Unknown();
  1488.           fi;
  1489.           if IsUnknown( SFI ) then
  1490.              InfoCharTable2( "#I Decreased: 2nd power map not available ",
  1491.                              "or not unique,\n",
  1492.                              "#I  no test with 'Indicator'\n" );
  1493.           fi;
  1494.  
  1495.    # test if 'char1' can be an irreducible character
  1496.           if ForAny( char1, x -> IsRat(x) and not IsInt(x) ) 
  1497.              or ScalarProduct( tbl, char1, char1 ) <> 1 
  1498.              or char1[1] = 0
  1499.              or ( IsInt( SFI ) and ( ( IRS and AbsInt( SFI ) <> 1 ) or
  1500.                                      ( not IRS and SFI <> 0 ) ) )   then
  1501.             InfoCharTable2( "#E Decreased : computation of ",
  1502.                   Ordinal( Length( irred ) ), " character failed\n" );
  1503.             return false;   
  1504.           else
  1505.  
  1506.     # irreducible character found
  1507.             Add( irred, Copy( char1 ) );
  1508.           fi;
  1509.        else
  1510.  
  1511.     # what a pity (!), some reducible character remaining
  1512.           if char[1] < 0 then
  1513.              char := - char;
  1514.           fi;
  1515.           if char <> nc then
  1516.              redcount := redcount + 1;
  1517.              red[redcount] := Copy( char );
  1518.              for i in [1..m] do  
  1519.                 remmat[i][redcount] := solmat[i][j];
  1520.                 if solmat[i][j] <> 0 then
  1521.                    delline[i] := false;
  1522.                 fi;
  1523.              od;
  1524.           fi;
  1525.        fi;
  1526.        j := j+1;
  1527.     od;
  1528.     i1 := 0;
  1529.     remmat2 := [];
  1530.     for i in [1..m] do  
  1531.        if not delline[i] then
  1532.           i1 := i1 + 1;
  1533.           remmat2[i1] := remmat[i];
  1534.        fi;
  1535.     od;
  1536.     return rec( irreducibles := irred,
  1537.                 remainders := red, matrix := remmat2 );
  1538.     end;
  1539.  
  1540.  
  1541. #############################################################################
  1542. ##
  1543. #F  OrthogonalEmbeddings( <grammat> [, \"positive\" ] [, <integer> ] )
  1544. ##
  1545. OrthogonalEmbeddings := function( arg )
  1546.     local 
  1547.     # sonstige prozeduren
  1548.           Symmatinv,
  1549.     # variablen fuer Embed
  1550.           maxdim, M, D, s, phi, mult, m, x, t, x2, sumg, sumh, 
  1551.           f, invg, sol, solcount, out, 
  1552.           l, g, nullv, i, j, k, n, kgv, a, IdMat, chpo,
  1553.     # booleans
  1554.           positiv, checkpositiv, checkdim,
  1555.     # prozeduren fuer Embed
  1556.           comp1, comp2, scp2, multiples, solvevDMtr, 
  1557.           Dextend, Mextend, inca, rnew,
  1558.           deca, algorithm;
  1559.     
  1560.     Symmatinv := function( b )
  1561.     # inverts symmetric matrices
  1562.     
  1563.     local n, i, j, l, k, c, d, ba, B, kgv, kgv1;
  1564.     n := Length( b );
  1565.     c := Copy( IdMat );
  1566.     d := [];
  1567.     B := [];
  1568.     kgv1 := 1;
  1569.     ba := Copy( IdMat );
  1570.     for i in [1..n] do  
  1571.        k := b[i][i];
  1572.        for j in [1..i-1] do  
  1573.           k := k - c[i][j] * c[i][j] * B[j];
  1574.        od;
  1575.        B[i] := k;
  1576.        for j in [i+1..n] do  
  1577.           k := b[j][i];
  1578.           for l in [1..i-1] do  
  1579.              k := k - c[i][l] * c[j][l] * B[l];
  1580.           od;
  1581.           if B[i] <> 0 then
  1582.              c[j][i] := k / B[i];
  1583.           else
  1584.              Error ( "matrix not invertable, ", Ordinal( i ), 
  1585.                      " column is linearly dependent" );
  1586.           fi;
  1587.        od;
  1588.     od;
  1589.     if B[n] = 0 then
  1590.        Error ( "matrix not invertable, ", Ordinal( i ), 
  1591.                " column is linearly dependent" );
  1592.     fi;
  1593.     for i in [1..n-1] do
  1594.        for j in [i+1..n] do  
  1595.           if c[j][i] <> 0 then
  1596.              for l in [1..i] do  
  1597.                 ba[j][l] := ba[j][l] - c[j][i] * ba[i][l];
  1598.              od;
  1599.           fi;
  1600.        od;
  1601.     od;
  1602.     for i in [1..n] do  
  1603.        for j in [1..i-1] do  
  1604.           ba[j][i] := ba[i][j];
  1605.           ba[i][j] := ba[i][j] / B[i];
  1606.        od;
  1607.        ba[i][i] := 1/B[i];
  1608.     od;
  1609.     for i in [1..n] do  
  1610.        d[i] := [];
  1611.        for j in [1..n] do  
  1612.           if i >= j then
  1613.              k := ba[i][j];
  1614.              l := i + 1;
  1615.           else
  1616.              l := j;
  1617.              k := 0;
  1618.           fi;
  1619.           while l <= n do
  1620.              k := k + ba[i][l] * ba[l][j];
  1621.              l := l + 1;
  1622.           od;
  1623.           d[i][j] := k;
  1624.           kgv1 := Lcm( kgv1, Denominator( k ) );
  1625.        od;
  1626.     od;
  1627.     for i in [1..n] do  
  1628.        for j in [1..n] do  
  1629.           d[i][j] := kgv1 * d[i][j];
  1630.        od;
  1631.     od;
  1632.     return rec( inverse := d, enuminator := kgv1 );
  1633.     end;
  1634.     
  1635.     # program embed
  1636.      
  1637.     comp1 := function( a, b ) 
  1638.     local i;
  1639.     if ( a[n+1] < b[n+1] ) then
  1640.       return false;
  1641.     elif ( a[n+1] > b[n+1] ) then
  1642.       return true;
  1643.     else
  1644.       for i in [ 1 .. n ] do
  1645.         if AbsInt( a[i] ) > AbsInt( b[i] ) then  
  1646.           return true;
  1647.         elif AbsInt( a[i] ) < AbsInt( b[i] ) then
  1648.           return false;
  1649.         fi;
  1650.       od;
  1651.     fi;
  1652.     return false;
  1653.     end;
  1654.     
  1655.     comp2 := function( a, b )
  1656.     local i, t;
  1657.     t := Length(a)-1;
  1658.     if a[t+1] < b[t+1] then
  1659.       return true;
  1660.     elif a[t+1] > b[t+1] then
  1661.       return false;
  1662.     else
  1663.       for i in [ 1 .. t ] do
  1664.         if a[i] < b[i] then
  1665.           return false;
  1666.         elif a[i] > b[i] then
  1667.           return true;
  1668.         fi;
  1669.       od;
  1670.     fi;
  1671.     return false;
  1672.     end;
  1673.     
  1674.     scp2 := function( k, l )
  1675.     # uses    x, invg, 
  1676.     # changes 
  1677.     local   i, j, sum, sum1;
  1678.     
  1679.     sum := 0;
  1680.     for i in [1..n] do
  1681.        sum1 := 0;
  1682.        for j in [1..n] do  
  1683.           sum1 := sum1 + x[k][j] * invg[j][i];
  1684.        od;
  1685.        sum := sum + sum1 * x[l][i];
  1686.     od;
  1687.     return sum;
  1688.     end;
  1689.     
  1690.     multiples := function( l )
  1691.     # uses    m, phi, 
  1692.     # changes mult, s, k, a, sumh, sumg, 
  1693.     local   v, r, i, j, break;
  1694.     
  1695.     for j in [1..n] do  
  1696.        sumh[j] := 0;
  1697.     od;
  1698.     i := l;
  1699.     while i <= t and ( not checkdim or s <= maxdim ) do
  1700.        if mult[i] * phi[i][i] < m then
  1701.           break := false;
  1702.           repeat
  1703.              v := solvevDMtr( i );
  1704.              if not IsBound( v[1] ) or not IsList( v[1] ) then
  1705.                 r := rnew( v, i );
  1706.                 if r >= 0 then
  1707.                    if r > 0 then
  1708.                       Dextend( r );
  1709.                    fi;
  1710.                    Mextend( v, i );
  1711.                    a[i] := a[i] + 1;
  1712.                 else
  1713.                    break := true;
  1714.                 fi;
  1715.              else
  1716.                 break := true;
  1717.              fi;
  1718.           until a[i] * phi[i][i] >= m or ( checkdim and s > maxdim )
  1719.                 or break;
  1720.           mult[i] := a[i];
  1721.           while a[i] > 0 do
  1722.              s := s - 1;
  1723.              if M[s][Length( M[s] )] = 1 then
  1724.                 k := k -1;
  1725.              fi;
  1726.              a[i] := a[i] - 1;
  1727.           od;
  1728.        fi;  
  1729.        if mult[i] <> 0 then
  1730.           for j in [1..n] do  
  1731.              sumh[j] := sumh[j] + mult[i] * x2[i][j];
  1732.           od;
  1733.        fi;
  1734.        i := i + 1;
  1735.     od;
  1736.     end;
  1737.     
  1738.     solvevDMtr := function( l )
  1739.     #  uses    M, D, phi, f, 
  1740.     #  changes 
  1741.     local   M1, M2, i, j, k1, v, sum;
  1742.     
  1743.     k1 := 1;
  1744.     v := [];
  1745.     i := 1;
  1746.     while i < s do
  1747.        sum := 0;
  1748.        M1 := Length( M[i] );
  1749.        M2 := M[i][M1];
  1750.        for j in [1..M1-1] do  
  1751.           sum := sum + v[j] * M[i][j];
  1752.        od;
  1753.        if M2 = 1 then
  1754.           v[k1] := -( phi[l][f[i]] + sum ) / D[k1];
  1755.           k1 := k1 + 1;
  1756.        else
  1757.           if Denominator( sum ) <> 1
  1758.           or Numerator( sum ) <> -phi[l][f[i]] then
  1759.              v[1] := [];
  1760.              i := s;
  1761.           fi;
  1762.        fi;
  1763.        i := i + 1;
  1764.     od;
  1765.     return( v );
  1766.     end;
  1767.     
  1768.     inca := function( l )
  1769.     #  uses    x2, 
  1770.     #  changes l, a, sumg, sumh, 
  1771.     local   v, r, break, i;
  1772.     
  1773.     while l <= t and ( not checkdim or s <= maxdim ) do
  1774.        break := false;
  1775.        repeat
  1776.           v := solvevDMtr( l );
  1777.           if not IsBound( v[1] ) or not IsList( v[1] ) then
  1778.              r := rnew( v, l );
  1779.              if r >= 0 then
  1780.                 if r > 0 then
  1781.                    Dextend( r );
  1782.                 fi;
  1783.                 Mextend( v, l );
  1784.                 a[l] := a[l] + 1;
  1785.                 for i in [1..n] do  
  1786.                    sumg[i] := sumg[i] + x2[l][i];
  1787.                 od;
  1788.              else
  1789.                 break := true;
  1790.              fi;
  1791.           else
  1792.              break := true;
  1793.           fi;
  1794.        until a[l] >= mult[l] or ( checkdim and s > maxdim ) or break;
  1795.        mult[l] := 0;
  1796.        l := l + 1;
  1797.     od;
  1798.     return l; 
  1799.     end; 
  1800.     
  1801.     rnew := function( v, l )
  1802.     #  uses    phi, m, k, D, 
  1803.     #  changes v, 
  1804.     local   sum, i;
  1805.     sum := m - phi[l][l];
  1806.     for i in [1..k-1] do  
  1807.        sum := sum - v[i] * D[i] * v[i];
  1808.     od;
  1809.     if sum >= 0 then
  1810.      if sum > 0 then
  1811.        v[k] := 1;
  1812.      else
  1813.        v[k] := 0;
  1814.      fi;
  1815.     fi;
  1816.     return sum;
  1817.     end;
  1818.     
  1819.     Mextend := function( line, l ) 
  1820.     #  uses    D,  
  1821.     #  changes M, s, f, 
  1822.     local   i;
  1823.     for i in [1..Length( line )-1] do
  1824.        line[i] := line[i] * D[i];
  1825.     od;
  1826.     M[s] := line;
  1827.     f[s] := l;
  1828.     s := s + 1;
  1829.     end;
  1830.     
  1831.     Dextend := function( r )
  1832.     #  uses    a, 
  1833.     #  changes k, D, 
  1834.     D[k] := r;
  1835.     k := k + 1;
  1836.     end;
  1837.     
  1838.     deca := function( l )
  1839.     #  uses    x2, t, M, 
  1840.     #  changes l, k, s, a, sumg, 
  1841.     local   i;
  1842.     if l <> 1 then
  1843.        l := l - 1;
  1844.        if l = t - 1 then
  1845.           while a[l] > 0 do
  1846.              s := s -1;
  1847.              if M[s][Length( M[s] )] = 1 then
  1848.                 k := k - 1;
  1849.              fi;
  1850.              a[l] := a[l] - 1;
  1851.              for i in [1..n] do  
  1852.                 sumg[i] := sumg[i] - x2[l][i];
  1853.              od;
  1854.           od;
  1855.           l := deca( l );
  1856.        else
  1857.           if a[l] <> 0 then
  1858.              s := s - 1;
  1859.              if M[s][Length( M[s] )] = 1 then
  1860.                 k := k - 1;
  1861.              fi;
  1862.              a[l] := a[l] - 1;
  1863.              for i in [1..n] do  
  1864.                 sumg[i] := sumg[i] - x2[l][i];
  1865.              od;
  1866.              l := l + 1;
  1867.           else
  1868.              l := deca( l );
  1869.           fi;
  1870.        fi;
  1871.     fi;
  1872.     return l;
  1873.     end;
  1874.     
  1875.     # check input
  1876.     if not IsList( arg[1] ) or not IsList( arg[1][1] ) then 
  1877.        Error( "first argument must be symmetric Gram matrix\n",
  1878.               "usage : Orthog... ( < gram-matrix > \n",
  1879.               " [, <\"positive\"> ] [, < integer > ] )" );
  1880.     fi;
  1881.     if Length( arg[1] ) <> Length( arg[1][1] ) then
  1882.        Error( "number of lines and columns not identic\n",
  1883.               "usage : Orthog... ( < gram-matrix >\n",
  1884.               " [, <\"positive\"> ] [, < integer > ] )" );
  1885.     fi;
  1886.     g := Copy ( arg[1] );
  1887.     checkpositiv := false;
  1888.     checkdim := false;
  1889.     chpo := "xxx";
  1890.     if IsBound( arg[2] ) then
  1891.        if IsString( arg[2] ) then
  1892.           chpo := arg[2];
  1893.           if arg[2] = "positive" then
  1894.              checkpositiv := true;
  1895.           fi;
  1896.        else
  1897.           if IsInt( arg[2] ) then
  1898.              maxdim := arg[2];
  1899.              checkdim := true;
  1900.           else
  1901.              Error( "second argument must be string or integer\n",
  1902.                     "usage : Orthog... ( < gram-matrix >\n",
  1903.                     " [, <\"positive\"> ] [, < integer > ] )" );
  1904.           fi;
  1905.        fi;
  1906.     fi;
  1907.     if IsBound( arg[3] ) then
  1908.        if IsString( arg[3] ) then
  1909.           chpo := arg[3];
  1910.           if arg[3] = "positive" then
  1911.              checkpositiv := true;
  1912.           fi;
  1913.        else
  1914.           if IsInt( arg[3] ) then
  1915.              maxdim := arg[3];
  1916.              checkdim := true;
  1917.           else
  1918.              Error( "third argument must be string or integer\n",
  1919.                     "usage : Orthog... ( < gram-matrix >\n",
  1920.                     " [, <\"positive\"> ] [, < integer > ] )" );
  1921.           fi;
  1922.        fi;
  1923.     fi;
  1924.     n := Length( g );
  1925.     for i in [1..n] do  
  1926.        for j in [1..i-1] do  
  1927.           if g[i][j] <> g[j][i] then
  1928.              Error( "matrix not symmetric \n",
  1929.                     "usage : Orthog... ( < gram-matrix >\n",
  1930.                     " [, <\"positive\"> ] [, < integer > ] )" );
  1931.           fi;
  1932.        od;
  1933.     od;
  1934.     
  1935.     # main program
  1936.     IdMat := IdentityMat( n );
  1937.     invg  := Symmatinv( g );
  1938.     m     := invg.enuminator;
  1939.     invg  := invg.inverse;
  1940.     x     := ShortestVectors( invg, m, chpo );
  1941.     t     := Length(x.vectors); 
  1942.     for i in [1..t] do
  1943.        x.vectors[i][n+1] := x.norms[i];
  1944.     od;
  1945.     x    := x.vectors;
  1946.     M    := [];
  1947.     M[1] := [];
  1948.     D    := [];
  1949.     mult := [];
  1950.     sol  := [];
  1951.     f    := [];
  1952.     solcount := 0;
  1953.     s        := 1;
  1954.     k        := 1;
  1955.     l        := 1; 
  1956.     a        := [];
  1957.     for i in [1..t] do  
  1958.        a[i]      := 0;
  1959.        x[i][n+2] := 0;
  1960.        mult[i]   := 0;
  1961.     od;
  1962.     sumg := [];
  1963.     sumh := [];
  1964.     for i in [1..n] do  
  1965.        sumg[i] := 0;
  1966.        sumh[i] := 0;
  1967.     od;
  1968.     Sort(x,comp1);
  1969.     x2 := [];
  1970.     for i in [1..t] do  
  1971.        x2[i] := [];
  1972.        for j in [1..n] do  
  1973.           x2[i][j]  := x[i][j] * x[i][j];
  1974.           x[i][n+2] := x[i][n+2] + x2[i][j];
  1975.        od;
  1976.     od;
  1977.     phi := [];
  1978.     for i in [1..t] do  
  1979.        phi[i] := [];
  1980.        for j in [1..i-1] do  
  1981.           phi[i][j] := scp2( i, j );
  1982.        od;
  1983.        phi[i][i] := x[i][n+1];
  1984.     od;
  1985.     repeat
  1986.        multiples( l );
  1987.  
  1988.        # (former call of 'tracecond')    
  1989.        if ForAll( [ 1 .. n ], i -> g[i][i] - sumg[i] <= sumh[i] ) then
  1990.           l := inca( l );
  1991.           if s-k = n then
  1992.              solcount := solcount + 1;
  1993.              InfoCharTable2("#I OrthogonalEmbeddings: ",
  1994.                                 solcount," solutions found\n");
  1995.              sol[solcount] := [];
  1996.              for i in [1..t] do  
  1997.                 sol[solcount][i] := a[i];
  1998.              od;
  1999.              sol[solcount][t+1]  := s - 1;
  2000.           fi;
  2001.        fi;
  2002.        l := deca( l );
  2003.     until l <= 1; 
  2004.     out := rec( vectors := [], norms := [], solutions := [] );
  2005.     for i in [1..t] do  
  2006.        out.vectors[i] := [];
  2007.        out.norms[i]   := x[i][n+1]/m;
  2008.        for j in [1..n] do  
  2009.           out.vectors[i][j] := x[i][j];
  2010.        od;
  2011.     od;
  2012.     Sort( sol, comp2 );
  2013.     for i in [1..solcount] do
  2014.        out.solutions[i]  := [];
  2015.        for j in [1..t] do  
  2016.           for k in [1..sol[i][j]] do
  2017.             Add( out.solutions[i], j );
  2018.           od;
  2019.        od;
  2020.     od;
  2021.     return out;
  2022.     end;
  2023.  
  2024.  
  2025. #############################################################################
  2026. ##
  2027. #F  OrthogonalEmbeddingsSpecialDimension( <tbl>, <reducibles>, <grammat>,
  2028. #F                                        [, \"positive\" ], <integer> )
  2029. ##
  2030. OrthogonalEmbeddingsSpecialDimension := function ( arg )  
  2031.     local  red, dim, reducibles, matrix, tbl, emb, dec, i, s, irred;
  2032.     # check input
  2033.     if Length( arg ) < 4 then
  2034.        Error( "please specify desired dimension\n",
  2035.               "usage : Orthog...( <tbl>, <reducibles>,\n",
  2036.               "<gram-matrix>, [, \"positive\" ], <integer> )" );
  2037.     fi;
  2038.     if IsInt( arg[4] ) then
  2039.        dim := arg[4];
  2040.     else
  2041.        if IsBound( arg[5] ) then
  2042.           if IsInt( arg[5] ) then
  2043.              dim := arg[5];
  2044.           else
  2045.        Error( "please specify desired dimension\n",
  2046.               "usage : Orthog...( <tbl>, < reducibles >,\n",
  2047.               "< gram-matrix >, [, <\"positive\"> ], < integer > )" );
  2048.           fi;
  2049.        fi;
  2050.     fi;
  2051.     tbl := arg[1];
  2052.     reducibles := arg[2];
  2053.     if Length( arg ) = 4 then
  2054.        emb := OrthogonalEmbeddings( arg[3], arg[4] );
  2055.     else
  2056.        emb := OrthogonalEmbeddings( arg[3], arg[4], arg[5] );
  2057.     fi;
  2058.     s := [];
  2059.     for i in [1..Length(emb.solutions)] do
  2060.        if Length( emb.solutions[i] ) = dim then
  2061.           Add( s, Sublist( emb.vectors, emb.solutions[i] ) );
  2062.        fi;
  2063.     od;
  2064.     dec:= List( s, x -> Decreased( tbl, reducibles, x ) );
  2065.     dec:= Filtered( dec, x -> x <> false );
  2066.     if dec = [] then
  2067.       InfoCharTable2( "#I OrthogonalE...: no embedding",
  2068.                       " corresp. to characters\n" );
  2069.       return rec( irreducibles:= [], remainders:= reducibles );
  2070.     fi;
  2071.     irred:= Set( dec[1].irreducibles );
  2072.     for i in [2..Length(dec)] do
  2073.        IntersectSet( irred, dec[i].irreducibles );
  2074.     od;
  2075.     red:= Reduced( tbl, irred, reducibles );
  2076.     Append( irred, red.irreducibles );
  2077.     return rec( irreducibles:= irred, remainders:= red.remainders );
  2078.     end;
  2079.  
  2080.  
  2081. #############################################################################
  2082. ##
  2083. #F  DnLattice( <tbl>, <g1>, <y1> )
  2084. ##
  2085. DnLattice := function( tbl, g1, y1 )
  2086.     local
  2087.     # indices
  2088.       i, i1, j, j1, k, k1, l, next,
  2089.     # booleans
  2090.       empty, change, used, addable, SFIbool,
  2091.     # dimensions
  2092.       m, n,
  2093.     # help fields
  2094.       found, foundpos,
  2095.       z, zw, nullcount, nullgenerate, 
  2096.       maxentry, max, ind, irred, irredcount, red,
  2097.       blockcount, blocks, perm, addtest, preirred,
  2098.     # Gram matrix
  2099.       g, gblock,
  2100.     # characters
  2101.       y, y2,
  2102.     # variables for recursion
  2103.       root, rootcount, solution, ligants, ligantscount, glblock, begin,
  2104.       depth, choice, ende, sol,
  2105.     # functions    
  2106.       callreduced, nullset, maxset, Search, Add, DnSearch, test;
  2107.     
  2108.     # counts zeroes in given line
  2109.     nullset := function( g, i )
  2110.     local j;
  2111.     
  2112.     nullcount[ i ] := 0;
  2113.     for j in [ 1..n ] do
  2114.        if g[ j ] = 0 then
  2115.           nullcount[ i ] := nullcount[ i ] + 1;
  2116.        fi;
  2117.     od;
  2118.     end;
  2119.     
  2120.     # searches line with most non-zero-entries
  2121.     maxset := function( )
  2122.     local i;
  2123.     
  2124.     maxentry := 1;
  2125.     max := n;
  2126.     for i in [ 1..n ] do
  2127.        if nullcount[ i ] < max then
  2128.           max := nullcount[ i ];
  2129.           maxentry := i;
  2130.        fi;
  2131.     od;
  2132.     end;
  2133.     
  2134.     # searches lines to add in order to produce zeroes
  2135.     Search := function( j )
  2136.     local signum;
  2137.     
  2138.     nullgenerate := 0;
  2139.     if g[ j ][ maxentry ] > 0 then
  2140.        signum := -1;
  2141.        for k in [ 1..n ] do
  2142.           if k <> maxentry and k <> j then
  2143.              if g[ maxentry ][ k ] <> 0 then
  2144.                 if g[ j ][ k ] = g[ maxentry ][ k ] then
  2145.                    nullgenerate := nullgenerate + 1;
  2146.                 else
  2147.                    nullgenerate := nullgenerate - 1;
  2148.                 fi;
  2149.              fi;
  2150.           fi;
  2151.        od;
  2152.     else
  2153.        if g[ j ][ maxentry ] < 0 then
  2154.           signum := 1;
  2155.           for k in [ 1..n ] do
  2156.              if k <> maxentry and k <> j then 
  2157.                 if g[ maxentry ][ k ] <> 0 then
  2158.                    if g[ j ][ k ] = -g[ maxentry ][ k ] then
  2159.                       nullgenerate := nullgenerate + 1;
  2160.                    else
  2161.                       nullgenerate := nullgenerate - 1;
  2162.                    fi;
  2163.                 fi;
  2164.              fi;
  2165.           od;
  2166.        fi;
  2167.     fi;
  2168.     if nullgenerate > 0 then
  2169.        change := true;
  2170.        Add( j, maxentry );
  2171.        j := j + 1;
  2172.     fi;
  2173.     end;
  2174.     
  2175.     # adds two lines/columns
  2176.     Add := function( i, j )
  2177.     local k;
  2178.     
  2179.        y[ i ] := y[ i ] - g[ i ][ j ] * y[ j ];  
  2180.        g[ i ] := g[ i ] - g[ i ][ j ] * g[ j ];
  2181.        for k in [ 1..i-1 ] do
  2182.           g[ k ][ i ] := g[ i ][ k ];
  2183.        od;
  2184.        g[ i ][ i ] := 2;
  2185.        for k in [ i+1..n ] do
  2186.           g[ k ][ i ] := g[ i ][ k ];
  2187.        od;
  2188.     end;
  2189.     
  2190.     # backtrack-search for dn-lattice
  2191.     DnSearch := function( begin, depth, oldchoice )
  2192.     local connections, connect, i1, j1, choice, found;
  2193.     
  2194.     choice := Copy( oldchoice );
  2195.     if depth = 3 then
  2196.        # d4-lattice found !!!
  2197.        solution := 1;
  2198.        ende := true;
  2199.        if n > 4 then
  2200.           i1 := 0;
  2201.           found := false;
  2202.           while not found and i1 < n do
  2203.              i1 := i1 + 1;
  2204.              if i1 <> root[ j ] and i1 <> choice[ 1 ] 
  2205.              and i1 <> choice[ 2 ] and i1 <> choice[ 3 ] then
  2206.                 connections := 0;
  2207.                 for j1 in [1..3] do
  2208.                    if gblock[ i1 ][ choice[ j1 ] ] <> 0 then
  2209.                       connections := connections + 1;
  2210.                       connect := choice[ j1 ];
  2211.                    fi;
  2212.                 od;
  2213.                 if connections = 1 then
  2214.                    found := true;
  2215.                    choice[ 4 ] := connect;
  2216.                    solution := solution + 1;
  2217.                 fi;
  2218.              fi;
  2219.              i1 := i1 + 1;
  2220.           od;
  2221.        fi;
  2222.        sol := choice;
  2223.     else
  2224.        i1 := begin;
  2225.        while not ende and i1 <= ligantscount do
  2226.           found := true;
  2227.           for j1 in [1..depth] do
  2228.              if gblock[ ligants[ i1 ] ][ choice[ j1 ] ] <> 0 then
  2229.                 found := false;
  2230.              fi;
  2231.           od;
  2232.           if found then
  2233.              depth := depth + 1;
  2234.              choice[ depth ] := ligants[ i1 ];
  2235.              DnSearch( i1 + 1, depth, choice );
  2236.              depth := depth - 1;
  2237.           else
  2238.              i1 := i1 + 1;
  2239.           fi;
  2240.           if ligantscount - i1 + 1 + depth < 3 then
  2241.              ende := true;
  2242.           fi;
  2243.        od;
  2244.     fi;
  2245.     end;
  2246.     
  2247.     test := function(z)
  2248.     # some tests for the found characters
  2249.     local result, IRS, SFI, i1, y1, ind, testchar;
  2250.     testchar := Copy( z )/2;
  2251.     result := true;
  2252.     IRS := ForAll( testchar, x -> GaloisCyc(x,-1) = x );
  2253.     if IsBound( tbl.powermap ) and IsBound( tbl.powermap[2] ) then
  2254.        if IsList( tbl.powermap[2] ) and
  2255.           ForAll( tbl.powermap[2], IsInt ) then
  2256.           SFI := Indicator( tbl, [testchar], 2 )[1];
  2257.           SFIbool := true;
  2258.        else
  2259.           InfoCharTable2
  2260.              ("#I DnLattice: 2nd power map not available or not unique,\n",
  2261.               "#I            cannot test with Indicator\n");
  2262.           SFIbool := false;
  2263.        fi;
  2264.     else
  2265.       InfoCharTable2
  2266.          ("#I DnLattice: 2nd power map not availabe\n",
  2267.           "#I            cannot test with Indicator\n");
  2268.       SFIbool := false;
  2269.     fi;
  2270.     if SFIbool then
  2271.        if ForAny( testchar, x -> IsRat(x) and not IsInt(x) ) 
  2272.           or ScalarProduct( tbl, testchar, testchar ) <> 1 
  2273.           or testchar[1] = 0  
  2274.           or ( IRS and AbsInt( SFI ) <> 1 ) 
  2275.           or ( not IRS and SFI <> 0 ) then
  2276.          result := false;
  2277.        fi;
  2278.     else
  2279.        if ForAny( testchar, x -> IsRat(x) and not IsInt(x) ) 
  2280.           or ScalarProduct( tbl, testchar, testchar ) <> 1
  2281.           or testchar[1] = 0 then
  2282.          result := false;
  2283.        fi;
  2284.     fi;
  2285.     return result;
  2286.     end;
  2287.     
  2288.     # reduce whole lattice with the found irreducible
  2289.     callreduced := function() 
  2290.     z[ 1 ] := z[ 1 ]/ 2 ;
  2291.     if ScalarProduct( tbl, z[ 1 ], z[ 1 ] ) = 1 then
  2292.        irredcount := irredcount + 1;
  2293.        if z[ 1 ][ 1 ] > 0 then
  2294.           irred[ irredcount ] := z[ 1 ];
  2295.        else
  2296.           irred[ irredcount ] := -z[ 1 ];
  2297.        fi;
  2298.        y1 := Sublist( y, [ blocks.begin[i] .. blocks.ende[i] ] );
  2299.        red := Reduced( tbl, z, y1 );
  2300.        irred := Concatenation( irred, red.irreducibles );
  2301.        irredcount := Length( irred );
  2302.        y2 := Concatenation( y2, red.remainders );
  2303.     fi;
  2304.     end;
  2305.     
  2306.     # check input parameters
  2307.     if not IsCharTable( tbl ) then
  2308.        Error( "first argument must be character-table\n",
  2309.               "usage: DnLattice( <tbl>, <gram-matrix>, <reducibles> )" );
  2310.     fi;
  2311.     empty := false;
  2312.     if g1 <> [] then
  2313.     if IsList( g1 ) and IsBound( g1[1] ) and IsList( g1[1] ) then
  2314.        g := Copy( g1 );
  2315.     else
  2316.        Error( "second argument must be Gram matrix of reducible characters\n",
  2317.               "usage: DnLattice( <tbl>, <grammat>, <reducibles> )" );
  2318.     fi;
  2319.     else
  2320.     empty := true;
  2321.     fi;
  2322.     if y1 <> [] then
  2323.     if IsList( y1 ) and IsBound( y1[1] ) and IsList( y1[1] ) then
  2324.        y := Copy( y1 );
  2325.     else
  2326.        Error( "third argument must be list of reducible characters\n",
  2327.               "usage: DnLattice( <tbl>, <gram-matrix>, <reducibles> )" );
  2328.     fi;
  2329.     else
  2330.     empty := true;
  2331.     fi;
  2332.     y2        := [  ];
  2333.     irred     := [  ];
  2334.     
  2335.     if not empty then
  2336.     
  2337.     n := Length( y );
  2338.     for i in [1..n] do
  2339.        if g[i][i] <> 2 then
  2340.           Error( "reducible characters don't have norm 2\n",
  2341.                 "usage: DnLattice( <tbl>, <gram-matrix>, <reducibles> )" );
  2342.        fi;
  2343.     od;
  2344.     # initialisations
  2345.     z         := [  ];
  2346.     used      := [  ];
  2347.     next      := [  ];
  2348.     nullcount := [  ];
  2349.     m := Length( y[ 1 ] );
  2350.     for i in [1..n] do
  2351.        used[i] := false;
  2352.     od;
  2353.     blocks := rec( begin := [ ], ende := [ ] );
  2354.     blockcount   := 0;
  2355.     irredcount   := 0;
  2356.     change       := true;
  2357.     while change do
  2358.        change := false;
  2359.        for i in [ 1..n ] do
  2360.           nullset( g[ i ], i );
  2361.        od;
  2362.        maxset( );
  2363.        while max < n-2 and not change do
  2364.           while maxentry <= n and not change do
  2365.              if nullcount[ maxentry ] <> max then
  2366.                 maxentry := maxentry + 1;
  2367.              else
  2368.                 j := 1;
  2369.                 while j < maxentry and not change do
  2370.                    Search( j );
  2371.                    j := j + 1;
  2372.                 od;
  2373.                 j := maxentry + 1;
  2374.                 while j <= n and not change do
  2375.                    Search( j );
  2376.                    j := j + 1;
  2377.                 od;
  2378.                 if not change then
  2379.                    maxentry := maxentry + 1;
  2380.                 fi;
  2381.              fi;
  2382.           od;
  2383.           if not change then
  2384.              max := max + 1;
  2385.              maxentry := 1;
  2386.           fi;
  2387.        od;
  2388.     
  2389.     # 2 step-search in order to produce zeroes
  2390.     # 2_0_Box-Method
  2391.        change := false;
  2392.        i := 1;
  2393.        while i <= n and not change do
  2394.           while i <= n and nullcount[ i ] > n-3 do
  2395.              i := i + 1;
  2396.           od;
  2397.           if i <= n then
  2398.              j := 1;
  2399.              while j <= n and not change do
  2400.                 while j <= n and g[ i ][ j ] <> 0 do
  2401.                    j := j + 1;
  2402.                 od;
  2403.                 if j <= n then
  2404.                    i1 := 1;
  2405.                    while i1 <= n and not change do
  2406.                       while i1 <= n
  2407.                       and ( i1 = i or i1 = j or g[ i1 ][ j ] = 0 ) do
  2408.                          i1 := i1 + 1;
  2409.                       od;
  2410.                       if i1 <= n then
  2411.                          addtest := g[ i ] - g[ i ][ i1 ] * g[ i1 ];
  2412.                          nullgenerate := 0;
  2413.                          addable := true;
  2414.                          for k in [ 1..n ] do
  2415.                             if addtest[ k ] = 0 then
  2416.                                nullgenerate := nullgenerate + 1;
  2417.                             else 
  2418.                                if AbsInt( addtest[ k ] ) > 1 then
  2419.                                   addable := false;
  2420.                                fi;
  2421.                             fi;
  2422.                          od;
  2423.                          if addable then
  2424.                             nullgenerate := nullgenerate - nullcount[ i ];
  2425.                             for k in [ 1..n ] do
  2426.                                if k <> i and k <> j then
  2427.                                   if addtest[ k ] 
  2428.                                      = addtest[ j ] * g[ j ][ k ] then
  2429.                                      if g[ j ][ k ] <> 0 then
  2430.                                         nullgenerate := nullgenerate + 1;
  2431.                                      fi;
  2432.                                   else
  2433.                                      if addtest[ k ] <> 0 then
  2434.                                         if g[ j ][ k ] = 0 then
  2435.                                            nullgenerate := nullgenerate - 1;
  2436.                                         else
  2437.                                            addable := false;
  2438.                                         fi;
  2439.                                      fi;
  2440.                                   fi;
  2441.                                fi;
  2442.                             od;
  2443.                             if nullgenerate > 0 and addable then
  2444.                                Add( i, i1 );
  2445.                                Add( j, i );
  2446.                                change := true;
  2447.                             fi;
  2448.                          fi;
  2449.                          i1 := i1 + 1;
  2450.                       fi;
  2451.                    od;
  2452.                    j := j + 1;
  2453.                 fi;
  2454.              od;
  2455.              i := i + 1;
  2456.           fi;
  2457.        od;
  2458.     od;
  2459.     i := 1;
  2460.     j := 0;
  2461.     next[ 1    ] := 1;
  2462.     while j < n do
  2463.        blockcount := blockcount + 1;
  2464.        blocks.begin[ blockcount ] := i;
  2465.        l := 0;
  2466.        used[ next [ i ] ] := true;
  2467.        j := j + 1;
  2468.        y2[ j ] := y[ next [ i ] ];
  2469.        while l >= 0 do
  2470.           for k in [ 1..n ] do
  2471.              if g[ next[ i ] ][ k ] <> 0 and not used[ k ] then
  2472.                 l := l + 1;
  2473.                 next[ i + l ] := k;
  2474.                 j := j + 1;
  2475.                 y2[ j ] := y[ k ];
  2476.                 used[ k ] := true;
  2477.              fi;
  2478.           od;
  2479.           i := i + 1;
  2480.           l := l - 1;
  2481.        od;
  2482.        blocks.ende[ blockcount ] := i - 1;
  2483.        k := 1;
  2484.        while k <= n and used[ k ] do
  2485.           k := k + 1;
  2486.        od;
  2487.        if k <= n then
  2488.           next[i] := k;
  2489.        fi;
  2490.     od;
  2491.     perm := PermList( next )^-1;
  2492.     for i in [1..n] do
  2493.        g[i] := Permuted( g[i], perm );
  2494.     od;
  2495.     g := Permuted( g, perm );
  2496.     y := y2;
  2497.     y2 := [  ];
  2498.     
  2499.     # search for d4/d5 - lattice
  2500.     for i in [1..blockcount] do
  2501.        n := blocks.ende[ i ] - blocks.begin[ i ] + 1;
  2502.        solution := 0;
  2503.        if n >= 4 then
  2504.           gblock := [  ];
  2505.           j1 := 0;
  2506.           for j in [ blocks.begin[ i ]..blocks.ende[ i ] ] do
  2507.              j1 := j1 + 1;
  2508.              gblock[ j1 ] := [  ];
  2509.              k1 := 0;
  2510.              for k in [ blocks.begin[ i ]..blocks.ende[ i ] ] do
  2511.                 k1 := k1 + 1;
  2512.                 gblock[ j1 ][ k1 ] := g[ j ][ k ];
  2513.              od;
  2514.           od;
  2515.           root      := [  ];
  2516.           rootcount := 0;
  2517.           for j in [1..n] do
  2518.              nullset( gblock[ j ], j );
  2519.              if nullcount[ j ] < n - 3 then
  2520.                 rootcount := rootcount + 1;
  2521.                 root[ rootcount ] := j;
  2522.              fi;
  2523.           od;
  2524.           j := 1;
  2525.           while solution = 0 and j <= rootcount do
  2526.              ligants := [  ];
  2527.              ligantscount := 0;
  2528.              for k in [1..n] do
  2529.                 if k <> root[ j ] and gblock[ root[ j ] ][ k ] <> 0 then
  2530.                    ligantscount := ligantscount + 1;
  2531.                    ligants[ ligantscount ] := k;
  2532.                 fi;
  2533.              od;
  2534.              begin := 1;
  2535.              depth := 0;
  2536.              choice := [  ];
  2537.              ende := false;
  2538.              DnSearch( begin, depth, choice );
  2539.              if solution > 0 then
  2540.                 choice := sol;
  2541.              fi;
  2542.              j := j + 1;
  2543.           od;
  2544.        fi;
  2545.  
  2546.     # test of the found irreducibles  
  2547.        if solution = 1 then
  2548.           # treatment of D4-lattice
  2549.           found := 0;
  2550.           preirred := Sublist( y, [ blocks.begin[i] .. blocks.ende[i] ] );
  2551.           z[1] := preirred[choice[1]] + preirred[choice[2]];
  2552.           if test(z[1]) then
  2553.              red := Reduced( tbl, preirred, [ z[1] ] );
  2554.              if ForAll( red.irreducibles, test ) then
  2555.                 found := found + 1;
  2556.                 foundpos := 1;
  2557.              fi;
  2558.           fi;
  2559.           z[2] := preirred[choice[1]] + preirred[choice[3]];
  2560.           if test(z[2]) then
  2561.              red := Reduced( tbl, preirred, [ z[2] ] );
  2562.              if ForAll( red.irreducibles, test ) then
  2563.                 found := found + 1;
  2564.                 foundpos := 2;
  2565.              fi;
  2566.           fi;
  2567.           z[3] := preirred[choice[2]] + preirred[choice[3]];
  2568.           if test(z[3]) then
  2569.              red := Reduced( tbl, preirred, [ z[3] ] );
  2570.              if ForAll( red.irreducibles, test ) then
  2571.                 found := found + 1;
  2572.                 foundpos := 3;
  2573.              fi;
  2574.           fi;
  2575.           if found = 1 then
  2576.              z := [z[foundpos]];
  2577.              callreduced();
  2578.           fi;  
  2579.         
  2580.        else
  2581.           # treatment of D5-lattice
  2582.           if solution = 2 then
  2583.              if choice [ 1 ] <> choice [ 4 ] then
  2584.                 z[ 1 ] := y[ blocks.begin[ i ] + choice[ 1 ] - 1 ];
  2585.                 if choice [ 2 ] <> choice [ 4 ] then
  2586.                    z[ 1 ] 
  2587.                         := z[ 1 ] + y[ blocks.begin[ i ] + choice[ 2 ] - 1 ];
  2588.                 else
  2589.                    z[ 1 ] 
  2590.                         := z[ 1 ] + y[ blocks.begin[ i ] + choice[ 3 ] - 1 ];
  2591.                 fi;
  2592.              else
  2593.                 z[ 1 ] := y[ blocks.begin[ i ] + choice[ 2 ] - 1 ]
  2594.                         + y[ blocks.begin[ i ] + choice[ 3 ] - 1 ];
  2595.              fi;
  2596.              found := 0;
  2597.              if test(z[1]) then
  2598.                 callreduced();
  2599.              fi;
  2600.           else
  2601.             Append( y2, Sublist( y, [ blocks.begin[i] .. blocks.ende[i] ] ) );
  2602.           fi;
  2603.        fi;
  2604.     od;            
  2605.     
  2606.     if irredcount > 0 then
  2607.        g := MatScalarProducts( tbl, y2, y2 );
  2608.     fi;
  2609.     else
  2610.        # input was empty i.e. empty=true
  2611.        g := [];
  2612.     fi;
  2613.     return rec( gram:=g, remainders:=y2, irreducibles:=irred );
  2614.     end;
  2615.  
  2616.  
  2617. #############################################################################
  2618. ##
  2619. #F  DnLatticeIterative( <tbl>, <red> )
  2620. ##
  2621. DnLatticeIterative := function( tbl, red )
  2622.     local dnlat, red1, norms, i, reduc, irred, norm2, g;
  2623.     
  2624.     # check input parameters
  2625.     if not IsCharTable( tbl ) then
  2626.        Error( "first argument must be character-table\n",
  2627.               "usage: DnLatticeIterative( <tbl>, <record or list> )" );
  2628.     fi;
  2629.     if not IsRec( red ) and not IsList( red ) then
  2630.        Error( "second argument must be record or list\n",
  2631.               "usage: DnLatticeIterative( <tbl>, <record or list> )" );
  2632.     fi;
  2633.     if IsRec( red ) and not IsBound( red.remainders ) then
  2634.        Error( "second record must contain a field 'remainders'\n",
  2635.               "usage: DnLatticeIterative( <tbl>, <record or list> )" );
  2636.     fi;
  2637.     if not IsRec( red ) then
  2638.        red := rec( remainders:=red );
  2639.     fi;
  2640.     if not IsBound( red.norms ) then
  2641.        norms := List( red.remainders, x -> ScalarProduct( tbl, x, x ) );
  2642.     else
  2643.        norms := Copy( red.norms );
  2644.     fi;
  2645.     reduc := Copy( red.remainders );
  2646.     irred := [];
  2647.     repeat
  2648.        norm2 := [];
  2649.        for i in [1..Length( reduc )] do
  2650.           if norms[i] = 2 then
  2651.              Add( norm2, reduc[i] );
  2652.           fi;
  2653.        od;
  2654.        g := MatScalarProducts( tbl, norm2, norm2 ); 
  2655.        dnlat := DnLattice( tbl, g, norm2 );
  2656.        Append( irred, dnlat.irreducibles );
  2657.        red1:= Reduced( tbl, dnlat.irreducibles, reduc );
  2658.        reduc := red1.remainders;
  2659.        Append( irred, red1.irreducibles );
  2660.        norms:= List( reduc, x -> ScalarProduct( tbl, x, x ) );
  2661.     until dnlat.irreducibles=[] and red1.irreducibles=[];
  2662.     return rec( irreducibles:=irred, remainders:=reduc , norms := norms );
  2663.     end;
  2664.  
  2665.