#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Thomas Breuer.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##
##  This file contains methods for lattices.
##


#############################################################################
##
#M  ScalarProduct( <v>, <w> ) . . . . . . . . . . . . . . for two row vectors
##
InstallMethod( ScalarProduct,
    "method for two row vectors",
    IsIdenticalObj,
    [ IsRowVector, IsRowVector ],
    function( v, w )
    return v * w;
    end );


#############################################################################
##
#F  StandardScalarProduct( <L>, <x>, <y> )
##
InstallGlobalFunction( StandardScalarProduct, function( L, x, y )
    return x * y;
end );


#############################################################################
##
#M  InverseMatMod( <intmat>, <prime> )
##
##  This method is much faster than the generic method in `matrix.gi'
##
InstallMethod( InverseMatMod,
    "method for a matrix, and an integer",
    [ IsMatrix, IsPosInt ],
    function( intmat, p )

    local i, j, k,              # loop variables
          n,                    # dimension
          intmatq, intmatqinv,  # matrix & inverse modulo p
          x,                    # solution of one iteration
          zline;                # help-line for exchange

    n:= Length( intmat );

    # `intmatq': `intmat' reduced mod `p'
    intmatq := intmat mod p;
    intmatqinv := IdentityMat( n );

    for i in [ 1 .. n ] do
      j := i;
      while j <= n and intmatq[j][i] = 0 do
        j := j + 1;
      od;
      if j > n then

        # matrix is singular modulo that `p'
        return fail;
      else

        # exchange lines `i' and `j'
        if j <> i then
           zline := intmatq[j];
           intmatq[j] := intmatq[i];
           intmatq[i] := zline;
           zline := intmatqinv[j];
           intmatqinv[j] := intmatqinv[i];
           intmatqinv[i] := zline;
        fi;

        # normalize line `i'
        zline:= intmatq[i];
        if zline[i] <> 1 then
          x:= (1/zline[i]) mod p;
          zline[i]:= 1;
          for k in [ i+1 .. n ] do
            if zline[k] <> 0 then
              zline[k]:= (x * zline[k]) mod p;
            fi;
          od;
          zline:= intmatqinv[i];
          for k in [1 .. n] do
            if zline[k] <> 0 then
              zline[k]:= (x * zline[k]) mod p;
            fi;
          od;
        fi;

        # elimination in column `i'
        for j in [ 1 .. n ] do
          if j <> i and intmatq[j][i] <> 0 then
            x:= intmatq[j][i];
            for k in [ 1 .. n ] do
              if intmatqinv[i][k] <> 0 then
                intmatqinv[j][k]:=
                    (intmatqinv[j][k] - x * intmatqinv[i][k] ) mod p;
              fi;
              if intmatq[i][k] <> 0 then
                intmatq[j][k]:=
                     (intmatq[j][k] - x * intmatq[i][k] ) mod p;
              fi;
            od;
          fi;
        od;

      fi;

    od;

    return intmatqinv;
    end );


#############################################################################
##
#F  PadicCoefficients( <A>, <Amodpinv>, <b>, <prime>, <depth> )
##
InstallGlobalFunction( PadicCoefficients,
    function( A, Amodpinv, b, prime, depth )
    local i, n, coeff, step, p2, val;

    n:= Length( b );
    coeff:= [];
    step:= 0;
    p2:= ( prime - 1 ) / 2;
    while PositionNonZero( b ) <= n and step < depth do
      step:= step + 1;
      coeff[ step ]:= ShallowCopy( b * Amodpinv );
      for i in [ 1 .. n ] do
        val:= coeff[ step ][i] mod prime;
        if val > p2 then
          coeff[ step ][i]:= val - prime;
        else
          coeff[ step ][i]:= val;
        fi;
      od;
      b:= ( b - coeff[ step ] * A ) / prime;
    od;
    Add( coeff, b );
    return coeff;
end );


#############################################################################
##
#F  LinearIndependentColumns( <mat> )
##
InstallGlobalFunction( LinearIndependentColumns, function( mat )
    local   m, n,       # dimensions of `mat'
            maxrank,    # maximal possible rank of `mat'
            i, j, k, q,
            row,
            zero,
            val,
            choice;     # list of linear independent columns, result

    # Make a copy to avoid changing the original argument.
    m := Length( mat );
    n := Length( mat[1] );
    maxrank:= m;
    if n < m then
      maxrank:= n;
    fi;
    zero := Zero( mat[1][1] );
    mat := List( mat, ShallowCopy );
    choice:= [];

    # run through all columns of the matrix
    i:= 0;
    for k in [1..n]  do

        # find a nonzero entry in this column
        j := i + 1;
        while j <= m and mat[j][k] = zero  do j := j+1;  od;

        # if there is a nonzero entry
        if j <= m  then

            i:= i+1;

            # Choose this column.
            Add( choice, k );
            if Length( choice ) = maxrank then
              return choice;
            fi;

            # Swap rows `j' and `i'.
            row:= mat[j];
            SwapMatrixRows(mat, i, j);

            # Normalize column `k'.
            val:= row[k];
            for q in [ j .. m ] do
              mat[q][k] := mat[q][k] / val;
            od;
            row[k]:= 1;

            # Clear all entries in row `i'.
            for j in [ k+1 .. n ] do
              if mat[i][j] <> zero  then
                val:= mat[i][j];
                for q in [ i .. m ] do
                  mat[q][j] := mat[q][j] - val * mat[q][k];
                od;
              fi;
            od;

        fi;

    od;

    # Return the list of positions of linear independent columns.
    return choice;
end );


#############################################################################
##
#F  DecompositionInt( <A>, <B>, <depth> )  . . . . . . . . integral solutions
##
##  returns the decomposition matrix <X> with `<X> * <A> = <B>' for <A> and
##  <B> integral matrices.
##
##  For an odd prime $p$, each integer $x$ has a unique representation
##  $x = \sum_{i=0}^{n} x_i p^i$ where  $|x_i| \leq \frac{p-1}{2}$ .
##  Let $x$ be a solution of the equation $xA = b$ where $A$ is a square
##  integral matrix and $b$ an integral vector, $\overline{A} = A \bmod p$
##  and $\overline{b} = b \bmod p$;
##  then $\overline{x} \overline{A} \equiv \overline{b} \bmod p$ for
##  $\overline{x} = x \bmod p$.
##  Assume $\overline{A}$ is regular over the field with $p$ elements; then
##  $\overline{x}$ is uniquely determined mod $p$.
##  Define $x^{\prime} = \frac{x - \overline{x}}{p}$ and
##         $b^{\prime} = \frac{b - \overline{x} A }{p}$.
##  If $y$ is a solution of the equation $x^{\prime} A = b^{\prime}$ we
##  have $( \overline{x} + p y ) A = b$ and thus $x = \overline{x} + p y$
##  is the solution of our problem.
##  Note that the process must terminate if an integral solution $x$ exists,
##  since the $p$--adic series for $y$ has one term less than that for $x$.
##
##  If $A$ is not square, it must have full rank,
##  and $'Length( <A> )' \leq `Length( <A>[1] )'$.
##
InstallGlobalFunction( DecompositionInt, function( A, B, depth )
    local i,       # loop variables
          Aqinv,      # inverse of matrix modulo p
          b,          # vector
          sol,        # solution of one step
          result,     # whole solution
          p,          # prime
          nullv,      # zero-vector
          origA,      # store full argument `A' in case of column choice
          n,          # dimension
          choice,
          coeff;

    # check input parameters
    if   Length( A ) > Length( A[1] ) then
      Error( "<A> must have at least `Length(<A>)' columns" );
    elif not IsMatrix( A ) and ForAll( A, x -> ForAll( x, IsInt ) ) then
      Error( "<A> must be integer matrix" );
    elif not ForAll( B, x -> x = fail or ( ForAll( x, IsInt )
                                  and Length( x ) = Length( A[1] ) ) ) then
      Error( "<B> must be list of integer vectors",
             " of same dimension as in <A>" );
    elif not IsInt( depth ) and depth >= 0 then
      Error( "<depth> (of iterations) must be a nonnegative integer" );
    fi;

    # initialisations
    n      := Length( A );
    depth  := depth + 1;
    result := [];
    p      := 83;
    nullv  := List( [ 1 .. n ], x -> 0 );

    # if `A' is not square choose `n' linear independent columns
    origA:= A;
    if Length( A[1] ) > n then

      choice:= LinearIndependentColumns( A );
      if Length( choice ) < Length( A ) then
        Error( "<A> has not full rank" );
      fi;
      A:= List( A, x -> x{ choice } );

    else
      choice:= [ 1 .. n ];
    fi;

    # compute the inverse `Aqinv' of `A' modulo `p';
    Aqinv:= InverseMatMod( A, p );
    while Aqinv = fail do

      # matrix is singular modulo that `p', choose another one
      p := NextPrimeInt( p );
      Info( InfoZLattice, 1,
            "DecompositionInt: choosing new prime ", p );
      Aqinv:= InverseMatMod( A, p );
    od;

    # compute the p-adic coefficients of the solutions,
    # and form the solutions
    for b in B do

      if b = fail then
        Add( result, fail );
      else
        b:= b{ choice };
        coeff:= PadicCoefficients( A, Aqinv, b, p, depth );
        if Last(coeff) = nullv then
          sol := nullv;
          for i in Reversed( [ 1 .. Length( coeff ) - 1 ] ) do
            sol := sol * p + coeff[i];
          od;
          Add( result, ShallowCopy( sol ) );
        else
          Add( result, fail );
        fi;
      fi;

    od;

    # if the argument `A' is not square test if the solutions are correct
    if Length( origA[1] ) > n then
      for i in [ 1 .. Length( result ) ] do
        if result[i] <> fail and result[i] * origA <> B[i] then
          result[i]:= fail;
        fi;
      od;
    fi;

    return result;
end );


#############################################################################
##
#F  IntegralizedMat( <A> )
#F  IntegralizedMat( <A>, <inforec> )
##
InstallGlobalFunction( IntegralizedMat, function( arg )
    local i, A, inforec, tr, f,
          stab,       # Galois stabilizer of `f'
          galaut, repr, aut, conj, pos, row, intA,
          col, introw, nofcyc,
          coeffs;     # coefficients of one column basis

    if Length( arg ) = 0 or Length( arg ) > 2 or not IsMatrix( arg[1] )
       or ( Length( arg ) = 2 and not IsRecord( arg[2] ) ) then
      Error( "usage: IntegralizedMat( <A> ) resp. \n",
             "       IntegralizedMat( <A>, <inforec> )" );
    fi;

    A:= arg[1];
    if Length( arg ) = 2 then

      # just use `inforec' to transform `A'
      inforec:= arg[2];

    else

      # compute transformed matrix `intA' and info record `inforec'
      inforec:= rec( intcols  := [],
                     irratcols:= [],
                     fields   := [] );
      tr:= MutableTransposedMat( A );

      for i in [ 1 .. Length( tr ) ] do

        if IsBound( tr[i] ) then

          if ForAll( tr[i], IsInt ) then
            Add( inforec.intcols, i );
          else

            # compute the field and the coefficients of values;
            # if `tr' contains conjugates of the row, delete them
            f:= FieldByGenerators( tr[i] );
            stab:= GaloisStabilizer( f );
            nofcyc:= Conductor( GeneratorsOfField( f ) );
            galaut:= PrimeResidues( nofcyc );
            SubtractSet( galaut, stab );
            repr:= [];
            while galaut <> [] do
              Add( repr, galaut[1] );
              SubtractSet( galaut,
                           List( stab * galaut[1], x -> x mod nofcyc) );
            od;
            for aut in repr do
              conj:= List( tr[i], x-> GaloisCyc( x, aut ) );
              pos:= Position( tr, conj, 0 );
              if pos <> fail then
                Unbind( tr[ pos ] );
              fi;
            od;
            inforec.fields[i]:= f;
            Add( inforec.irratcols, i );

          fi;
        fi;
      od;
    fi;

    intA:= [];
    for row in A do
      introw:= [];
      coeffs:= row{ inforec.intcols };
      if not ForAll( coeffs, IsInt ) then
        coeffs:= fail;
      fi;
      for col in inforec.irratcols do
        if coeffs <> fail then
          Append( introw, coeffs );
          coeffs:= Coefficients( CanonicalBasis( inforec.fields[ col ] ),
                                 row[ col ] );
        fi;
      od;
      if coeffs = fail then
        introw:= fail;
      else
        Append( introw, coeffs );
      fi;
      Add( intA, introw );
    od;

    return rec( mat:= intA, inforec:= inforec );
end );


#############################################################################
##
#F  Decomposition( <A>, <B>, <depth> ) . . . . . . . . . . integral solutions
#F  Decomposition( <A>, <B>, \"nonnegative\" ) . . . . . . integral solutions
##
##  For a matrix <A> of cyclotomics and a list <B> of cyclotomic vectors,
##  `Decomposition' tries to find integral solutions of the linear equation
##  systems `<x> \* <A> = <B>[i]'.
##
##  <A> must have full rank, i.e., there must be a linear independent set of
##  columns of same length as <A>.
##
##  `Decomposition( <A>, <B>, <depth> )', where <depth> is a nonnegative
##  integer, computes for every `<B>[i]' the initial part
##  $\Sum_{k=0}^{<depth>} x_k p^k$ (all $x_k$ integer vectors with entries
##  bounded by $\pm\frac{p-1}{2}$) of the $p$-adic series of a hypothetical
##  solution. The prime $p$ is 83 first; if the reduction of <A>
##  modulo $p$ is singular, the next prime is chosen automatically.
##
##  A list <X> is returned. If the computed initial part for
##  `<x> \* <A> = <B>[i]' *is* a solution, we have `<X>[i] = <x>', otherwise
##  `<X>[i] = false'.
##
##  `Decomposition( <A>, <B>, \"nonnegative\" )' assumes that the solutions
##  have only nonnegative entries.
##  This is e.g.\ satisfied for the decomposition of ordinary characters into
##  Brauer characters.
##  If the first column of <A> consists of positive integers,
##  the necessary number <depth> of iterations can be computed. In that case
##  the `i'-th entry of the returned list is `false' if there *exists* no
##  nonnegative integral solution of the system `<x> \* <A> = <B>[i]', and it
##  is the solution otherwise.
##
##  *Note* that the result is a list of `false' if <A> has not full rank,
##  even if there might be a unique integral solution for some equation
##  system.
##
InstallMethod( Decomposition, "for a matrix of cyclotomics, a vector and a depth",
        [IsMatrix,IsList,IsObject],
        function( A, B, depth_or_nonnegative )
    local i, intA, intB, newintA, newintB, result, choice, inforec;

    # Check the input parameters.
    if not ( IsInt( depth_or_nonnegative ) and depth_or_nonnegative >= 0 )
       and depth_or_nonnegative <> "nonnegative" then
      Error( "usage: Decomposition( <A>,<B>,<depth> ) for integer <depth>\n",
             "        resp. Decomposition( <A>,<B>,\"nonnegative\" )\n",
             "        ( solution <X> of <X> * <A> = <B> )" );
    elif not ( IsMatrix(A) and IsMatrix(B)
               and Length(B[1]) = Length(A[1]) ) then
      Error( "<A>, <B> must be matrices with same number of columns" );
    fi;

    # Transform `A' to an integer matrix `intA'.
    intA:= IntegralizedMat( A );
    inforec:= intA.inforec;
    intA:= intA.mat;

    # Transform `B' to `intB', choose coefficients compatible.
    intB:= IntegralizedMat( B, inforec ).mat;

    # If `intA' is not square then choose linear independent columns.
    if Length( intA ) < Length( intA[1] ) then
      choice:= LinearIndependentColumns( intA );
      newintA:= List( intA, x -> x{ choice } );
      newintB:= [];
      for i in [ 1 .. Length( intB ) ] do
        if intB[i] = fail then
          newintB[i]:= fail;
        else
          newintB[i]:= intB[i]{ choice };
        fi;
      od;
    elif Length( intA ) = Length( intA[1] ) then
      newintA:= intA;
      newintB:= intB;
    else
      Error( "There must be a subset of columns forming a regular matrix" );
    fi;

    # depth of iteration
    if depth_or_nonnegative = "nonnegative" then
      if not ForAll( newintA, x -> IsInt( x[1] ) and x[1] >= 0 ) then
        Error( "option \"nonnegative\" is allowed only if the first column\n",
               "       of <A> consists of positive integers" );
      fi;

      # The smallest value that has length `c' in the p-adic series is
      # p^c + \Sum_{k=0}^{c-1} -\frac{p-1}{2} p^k = \frac{1}{2}(p^c + 1).
      # So if $'<B>[i][1] / Minimum( newintA[1] )' \< \frac{1}{2}(p^c + 1)$
      # we have `depth' at most `c-1'.

      result:= DecompositionInt( newintA, newintB,
                       LogInt( 2*Int( Maximum( List( B, x->x[1] ) )
                               / Minimum( List( A, x -> x[1]) ) ), 83 ) + 2 );
      for i in [ 1 .. Length( result ) ] do
        if IsList( result[i] ) and Minimum( result[i] ) < 0 then
          result[i]:= fail;
        fi;
      od;
    else
      result:= DecompositionInt( newintA, newintB, depth_or_nonnegative );
    fi;

    # If `A' is not integral or `intA' is not square
    # then test if the result is correct.
    if Length( intA ) < Length( intA[1] )
       or not IsEmpty( inforec.irratcols) then
      for i in [ 1 .. Length( result ) ] do
        if result[i] <> fail and result[i] * A <> B[i] then
          result[i]:= fail;
        fi;
      od;
    fi;

    return result;
end );


#############################################################################
##
#F  LLLReducedBasis( [<L>, ]<vectors>[, <y>][, \"linearcomb\"][, <lllout>] )
##
InstallGlobalFunction( LLLReducedBasis, function( arg )
      local mmue,      # buffer $\mue$
            L,         # the lattice
            y,         # sensitivity $y$ (default $y = \frac{3}{4}$)
            kmax,      # $k_{max}$
            b,         # list $b$ of vectors
            H,         # basechange matrix $H$
            mue,       # matrix $\mue$ of scalar products
            B,         # list $B$ of norms of $b^{\ast}$
            BB,        # buffer $B$
            q,         # buffer $q$ for function `RED'
            i,         # loop variable $i$
            j,         # loop variable $j$
            k,         # loop variable $k$
            l,         # loop variable $l$
            n,         # number of vectors in $b$
            lc,        # boolean: option `linearcomb'?
            lllout,    # record with info about initial part of $b$
            scpr,      # scalar product of lattice `L'
            RED,       # reduction subprocedure; `RED( l )'
                       # means `RED( k, l )' in Cohen's book
            r;         # number of zero vectors found up to now

    RED := function( l )

    # Terminate for $\|\mue_{k,l}\| \leq \frac{1}{2}$.
    if 1 < mue[k][l] * 2 or mue[k][l] * 2 < -1 then

      # Let $q = `Round( mue[k][l] )'$ (is never zero), \ldots
#T Round ?
      q:= Int( mue[k][l] );
      if AbsoluteValue( mue[k][l] - q ) * 2 > 1 then
        q:= q + SignInt( mue[k][l] );
      fi;

      # \ldots and subtract $q b_l$ from $b_k$;
      AddRowVector( b[k], b[l], - q );

      # adjust `mue', \ldots
      mue[k][l]:= mue[k][l] - q;
      for i in [ r+1 .. l-1 ] do
        if mue[l][i] <> 0 then
          mue[k][i]:= mue[k][i] - q * mue[l][i];
        fi;
      od;

      # \ldots and the basechange.
      if lc then
        AddRowVector( H[k], H[l], - q );
      fi;

    fi;
    end;


    # Check the input parameters.
    if   IsLeftModule( arg[1] ) then
      L:= arg[1];
      scpr:= ScalarProduct;
      arg:= arg{ [ 2 .. Length( arg ) ] };
    elif IsList( arg[1] ) then
      # There is no lattice given, take standard scalar product.
      L:= "L";
      scpr:= StandardScalarProduct;
    else
      Error( "usage: LLLReducedBasis( [<L>], <vectors> [,<y>]",
             "[,\"linearcomb\"] )" );
    fi;

    b:= List( arg[1], ShallowCopy );

    # Preset the ``sensitivity'' (value between $\frac{1}{4}$ and $1$).
    if IsBound( arg[2] ) and IsRat( arg[2] ) then
      y:= arg[2];
      if y <= 1/4 or 1 < y then
        Error( "sensitivity `y' must satisfy 1/4 < y <= 1" );
      fi;
    else
      y:= 3/4;
    fi;

    # Get the optional `\"linearcomb\"' parameter
    # and the optional `lllout' record.
    lc     := false;
    lllout := false;

    for i in [ 2 .. Length( arg ) ] do
      if arg[i] = "linearcomb" then
        lc:= true;
      elif IsRecord( arg[i] ) then
        lllout:= arg[i];
      fi;
    od;


    # step 1 (Initialize.)
    n := Length( b );
    r := 0;
    i := 1;
    if lc then
      H:= IdentityMat( n );
    fi;

    if lllout = false or lllout.B = [] then

      k    := 2;
      mue  := [ [] ];
      kmax := 1;

      # Handle the case of leading zero vectors in the input.
      while i <= n and IsZero( b[i] ) do
        i:= i+1;
      od;
      if n < i then

        r:= n;
        k:= n+1;

      elif 1 < i then

        q    := b[i];
        b[i] := b[1];
        b[1] := q;
        if lc then
          q    := H[i];
          H[i] := H[1];
          H[1] := q;
        fi;

      fi;

      if 0 < n then
        B:= [ scpr( L, b[1], b[1] ) ];
      else
        B:= [];
      fi;

      Info( InfoZLattice, 1,
            "LLLReducedBasis called with ", n, " vectors, y = ", y );

    else

      # Note that the first $k_{max}$ vectors are all nonzero.

      mue  := List( lllout.mue, ShallowCopy );
      kmax := Length( mue );
      k    := kmax + 1;
      B    := ShallowCopy( lllout.B );

      Info( InfoZLattice, 1,
            "LLLReducedBasis (incr.) called with ",
            n, " = ", kmax, " + ", n - kmax, " vectors, y = ", y );

    fi;

    while k <= n do

      # step 2 (Incremental Gram-Schmidt)

      # If $k \leq k_{max}$ go to step 3.
      # Otherwise \ldots
      if kmax < k then

        Info( InfoZLattice, 2,
              "LLLReducedBasis: Take ", Ordinal( k ), " vector" );

        # \ldots set $k_{max} \leftarrow k$
        # and for $j = 1, \ldots, k-1$ set
        # $\mue_{k,j} \leftarrow b_k \cdot b_j^{\ast} / B_j$ if
        # $B_j \not= 0$ and $\mue_{k,j} \leftarrow 0$ if $B_j = 0$, \ldots
        kmax:= k;
        mue[k]:= [];
        for j in [ r+1 .. k-1 ] do
          mmue:= scpr( L, b[k], b[j] );
          for i in [ r+1 .. j-1 ] do
            mmue:= mmue - mue[j][i] * mue[k][i];
          od;
          mue[k][j]:= mmue;
        od;

        # (Now `mue[k][j]' contains $\mue_{k,j} B_j$ for all $j$.)
        for j in [ r+1 .. k-1 ] do
          mue[k][j]:= mue[k][j] / B[j];
        od;

        # \ldots then set $b_k^{\ast} \leftarrow
        # b_k - \sum_{j=1}^{k-1} \mue_{k,j} b_j^{\ast}$ and
        # $B_k \leftarrow b_k^{\ast} \cdot b_k^{\ast}$.
        B[k]:= scpr( L, b[k], b[k] );
        for j in [ r+1 .. k-1 ] do
          B[k]:= B[k] - mue[k][j]^2 * B[j];
        od;

      fi;

      # step 3 (Test LLL condition)
      RED( k-1 );
      while B[k] < ( y - mue[k][k-1] * mue[k][k-1] ) * B[k-1] do

        # Execute Sub-algorithm SWAPG$( k )$\:
        # Exchange the vectors $b_k$ and $b_{k-1}$,
        q      := b[k];
        b[k]   := b[k-1];
        b[k-1] := q;

        # $H_k$ and $H_{k-1}$,
        if lc then
          q      := H[k];
          H[k]   := H[k-1];
          H[k-1] := q;
        fi;

        # and if $k > 2$, for all $j$ such that $1 \leq j \leq k-2$
        # exchange $\mue_{k,j}$ with $\mue_{k-1,j}$.
        for j in [ r+1 .. k-2 ] do
          q           := mue[k][j];
          mue[k][j]   := mue[k-1][j];
          mue[k-1][j] := q;
        od;

        # Then set $\mue \leftarrow \mue_{k,k-1}$
        mmue:= mue[k][k-1];

        # and $B \leftarrow B_k + \mue^2 B_{k-1}$.
        BB:= B[k] + mmue^2 * B[k-1];

        # Now, in the case $B = 0$ (i.e. $B_k = \mue = 0$),
        if BB = 0 then

          # exchange $B_k$ and $B_{k-1}$
          B[k]   := B[k-1];
          B[k-1] := 0;

          # and for $i = k+1, k+2, \ldots, k_{max}$
          # exchange $\mue_{i,k}$ and $\mue_{i,k-1}$.
          for i in [ k+1 .. kmax ] do
            q           := mue[i][k];
            mue[i][k]   := mue[i][k-1];
            mue[i][k-1] := q;
          od;

        # In the case $B_k = 0$ and $\mue \not= 0$,
        elif B[k] = 0 and mmue <> 0 then

          # set $B_{k-1} \leftarrow B$,
          B[k-1]:= BB;

          # $\mue_{k,k-1} \leftarrow \frac{1}{\mue}
          mue[k][k-1]:= 1 / mmue;

          # and for $i = k+1, k+2, \ldots, k_{max}$
          # set $\mue_{i,k-1} \leftarrow \mue_{i,k-1} / \mue$.
          for i in [ k+1 .. kmax ] do
            mue[i][k-1]:= mue[i][k-1] / mmue;
          od;

        else

          # Finally, in the case $B_k \not= 0$,
          # set (in this order) $t \leftarrow B_{k-1} / B$,
          q:= B[k-1] / BB;

          # $\mue_{k,k-1} \leftarrow \mue t$,
          mue[k][k-1]:= mmue * q;

          # $B_k \leftarrow B_k t$,
          B[k]:= B[k] * q;

          # $B_{k-1} \leftarrow B$,
          B[k-1]:= BB;

          # then for $i = k+1, k+2, \ldots, k_{max}$ set
          # (in this order) $t \leftarrow \mue_{i,k}$,
          # $\mue_{i,k} \leftarrow \mue_{i,k-1} - \mue t$,
          # $\mue_{i,k-1} \leftarrow t + \mue_{k,k-1} \mue_{i,k}$.
          for i in [ k+1 .. kmax ] do
            q           := mue[i][k];
            mue[i][k]   := mue[i][k-1] - mmue * q;
            mue[i][k-1] := q + mue[k][k-1] * mue[i][k];
          od;

        fi;

        # Terminate the subalgorithm.

        if k > 2 then k:= k-1; fi;

        # Here we have always `k > r' since the loop is entered
        # for `k > r+1' only (because of `B[k-1] \<> 0'),
        # so the only problem might be the case `k = r+1',
        # namely `mue[ r+1 ][r]' is used then; but this is bound
        # provided that the initial list of vectors did not start
        # with zero vectors, and its (perhaps not updated) value
        # does not matter because this would mean just to subtract
        # a multiple of a zero vector.

        RED( k-1 );

      od;

      if B[ r+1 ] = 0 then
        r:= r+1;
        Unbind( b[r] );
      fi;

      for l in [ k-2, k-3 .. r+1 ] do
        RED( l );
      od;
      k:= k+1;

    # step 4 (Finished?)
    # If $k \leq n$ go to step 2.

    od;

    # Otherwise, let $r$ be the number of initial vectors $b_i$
    # which are equal to zero, output $p \leftarrow n - r$,
    # the vectors $b_i$ for $r+1 \leq i \leq n$ (which form an LLL-reduced
    # basis of $L$), the transformation matrix $H \in GL_n(\Z)$
    # and terminate the algorithm.

    # Check whether the last calls of `RED' have produced new zero vectors
    # in `b'; unfortunately this cannot be read off from `B'.
    while r < n and IsZero( b[ r+1 ] ) do
#T if this happens then is `B' outdated???
#T but `B' contains the norms of the orthogonal basis,
#T so this should be impossible!
#T (but if it happens then also `LLLReducedGramMat' should be adjusted!)
Print( "reached special case of increasing r in the last moment\n" );
if B[r+1] <> 0 then
  Print( "strange situation in LLL!\n" );
fi;
      r:= r+1;
    od;

    Info( InfoZLattice, 1,
          "LLLReducedBasis returns basis of length ", n-r );

    mue:= List( [ r+1 .. n ], i -> mue[i]{ [ r+1 .. i-1 ] } );
    MakeImmutable( mue );
    B:= B{ [ r+1 .. n ] };
    MakeImmutable( B );

    if lc then
      return rec( basis          := b{ [ r+1 .. n ] },
                  relations      := H{ [  1  .. r ] },
                  transformation := H{ [ r+1 .. n ] },
                  mue            := mue,
                  B              := B );
    else
      return rec( basis          := b{ [ r+1 .. n ] },
                  mue            := mue,
                  B              := B );
    fi;

end );


#############################################################################
##
#F  LLLReducedGramMat( <G>[, <y>] ) . . . . . . . . . LLL reduced Gram matrix
##
InstallGlobalFunction( LLLReducedGramMat, function( arg )
      local gram,      # the Gram matrix
            mmue,      # buffer $\mue$
            y,         # sensitivity $y$ (default $y = \frac{3}{4}$)
            kmax,      # $k_{max}$
            H,         # basechange matrix $H$
            mue,       # matrix $\mue$ of scalar products
            B,         # list $B$ of norms of $b^{\ast}$
            BB,        # buffer $B$
            q,         # buffer $q$ for function `RED'
            i,         # loop variable $i$
            j,         # loop variable $j$
            k,         # loop variable $k$
            l,         # loop variable $l$
            n,         # length of `gram'
            RED,       # reduction subprocedure; `RED( l )'
                       # means `RED( k, l )' in Cohen's book
            ak,        # buffer vector in Gram-Schmidt procedure
            r;         # number of zero vectors found up to now

    RED := function( l )

    # Terminate for $\|\mue_{k,l}\| \leq \frac{1}{2}$.
    if 1 < mue[k][l] * 2 or mue[k][l] * 2 < -1 then

      # Let $q = `Round( mue[k][l] )'$ (is never zero), \ldots
      q:= Int( mue[k][l] );
      if AbsoluteValue( mue[k][l] - q ) * 2 > 1 then
        q:= q + SignInt( mue[k][l] );
      fi;

      # \ldots adjust the Gram matrix (rows and columns, but only
      # in the lower triangular half), \ldots
      gram[k][k]:= gram[k][k] - q * gram[k][l];
      for i in [ r+1 .. l ] do
        gram[k][i]:= gram[k][i] - q * gram[l][i];
      od;
      for i in [ l+1 .. k ] do
        gram[k][i]:= gram[k][i] - q * gram[i][l];
      od;
      for i in [ k+1 .. n ] do
        gram[i][k]:= gram[i][k] - q * gram[i][l];
#T AddRowVector
      od;

      # \ldots adjust `mue', \ldots
      mue[k][l]:= mue[k][l] - q;
      for i in [ r+1 .. l-1 ] do
        if mue[l][i] <> 0 then
          mue[k][i]:= mue[k][i] - q * mue[l][i];
        fi;
      od;

      # \ldots and the basechange.
      H[k]:= H[k] - q * H[l];

    fi;
    end;


    # Check the input parameters.
    if arg[1] = [] or ( IsList( arg[1] ) and IsList( arg[1][1] ) ) then

      gram:= List( arg[1], ShallowCopy );

      # Preset the ``sensitivity'' (value between $\frac{1}{4}$ and $1$).
      if IsBound( arg[2] ) and IsRat( arg[2] ) then
        y:= arg[2];
        if y <= 1/4 or 1 < y then
          Error( "sensitivity `y' must satisfy 1/4 < y <= 1" );
        fi;
      else
        y:= 3/4;
      fi;

    else
      Error( "usage: LLLReducedGramMat( <gram>[, <y>] )" );
    fi;

    # step 1 (Initialize \ldots
    n    := Length( gram );
    k    := 2;
    kmax := 1;
    mue  := [ [] ];
    r    := 0;
    ak   := [];
    H    := IdentityMat( n );

    Info( InfoZLattice, 1,
          "LLLReducedGramMat called with matrix of length ", n,
          ", y = ", y );

    # \ldots and handle the case of leading zero vectors in the input.)
    i:= 1;
    while i <= n and gram[i][i] = 0 do
      i:= i+1;
    od;
    if i > n then

      r:= n;
      k:= n+1;

    elif i > 1 then

      for j in [ i+1 .. n ] do
        gram[j][1]:= gram[j][i];
        gram[j][i]:= 0;
      od;
      gram[1][1]:= gram[i][i];
      gram[i][i]:= 0;

      q    := H[i];
      H[i] := H[1];
      H[1] := q;

    fi;

    B:= [ gram[1][1] ];

    while k <= n do

      # step 2 (Incremental Gram-Schmidt)

      # If $k \leq k_{max}$ go to step 3.
      if k > kmax then

        Info( InfoZLattice, 2,
              "LLLReducedGramMat: Take ", Ordinal( k ), " vector" );

        # Otherwise \ldots
        kmax:= k;
        B[k]:= gram[k][k];
        mue[k]:= [];
        for j in [ r+1 .. k-1 ] do
          ak[j]:= gram[k][j];
          for i in [ r+1 .. j-1 ] do
            ak[j]:= ak[j] - mue[j][i] * ak[i];
          od;
          mue[k][j]:= ak[j] / B[j];
          B[k]:= B[k] - mue[k][j] * ak[j];
        od;

      fi;

      # step 3 (Test LLL condition)
      RED( k-1 );
      while B[k] < ( y - mue[k][k-1] * mue[k][k-1] ) * B[k-1] do

        # Execute Sub-algorithm SWAPG$( k )$\:
        # Exchange $H_k$ and $H_{k-1}$,
        q      := H[k];
        H[k]   := H[k-1];
        H[k-1] := q;

        # adjust the Gram matrix (rows and columns,
        # but only in the lower triangular half),
        for j in [ r+1 .. k-2 ] do
          q            := gram[k][j];
          gram[k][j]   := gram[k-1][j];
          gram[k-1][j] := q;
        od;
        for j in [ k+1 .. n ] do
          q            := gram[j][k];
          gram[j][k]   := gram[j][k-1];
          gram[j][k-1] := q;
        od;
        q              := gram[k-1][k-1];
        gram[k-1][k-1] := gram[k][k];
        gram[k][k]     := q;

        # and if $k > 2$, for all $j$ such that $1 \leq j \leq k-2$
        # exchange $\mue_{k,j}$ with $\mue_{k-1,j}$.
        for j in [ r+1 .. k-2 ] do
          q           := mue[k][j];
          mue[k][j]   := mue[k-1][j];
          mue[k-1][j] := q;
        od;

        # Then set $\mue \leftarrow \mue_{k,k-1}$
        mmue:= mue[k][k-1];

        # and $B \leftarrow B_k + \mue^2 B_{k-1}$.
        BB:= B[k] + mmue^2 * B[k-1];

        # Now, in the case $B = 0$ (i.e. $B_k = \mue = 0$),
        if BB = 0 then

          # exchange $B_k$ and $B_{k-1}$
          B[k]   := B[k-1];
          B[k-1] := 0;

          # and for $i = k+1, k+2, \ldots, k_{max}$
          # exchange $\mue_{i,k}$ and $\mue_{i,k-1}$.
          for i in [ k+1 .. kmax ] do
            q           := mue[i][k];
            mue[i][k]   := mue[i][k-1];
            mue[i][k-1] := q;
          od;

        # In the case $B_k = 0$ and $\mue \not= 0$,
        elif B[k] = 0 and mmue <> 0 then

          # set $B_{k-1} \leftarrow B$,
          B[k-1]:= BB;

          # $\mue_{k,k-1} \leftarrow \frac{1}{\mue}
          mue[k][k-1]:= 1 / mmue;

          # and for $i = k+1, k+2, \ldots, k_{max}$
          # set $\mue_{i,k-1} \leftarrow \mue_{i,k-1} / \mue$.
          for i in [ k+1 .. kmax ] do
            mue[i][k-1]:= mue[i][k-1] / mmue;
          od;

        else

          # Finally, in the case $B_k \not= 0$,
          # set (in this order) $t \leftarrow B_{k-1} / B$,
          q:= B[k-1] / BB;

          # $\mue_{k,k-1} \leftarrow \mue t$,
          mue[k][k-1]:= mmue * q;

          # $B_k \leftarrow B_k t$,
          B[k]:= B[k] * q;

          # $B_{k-1} \leftarrow B$,
          B[k-1]:= BB;

          # then for $i = k+1, k+2, \ldots, k_{max}$ set
          # (in this order) $t \leftarrow \mue_{i,k}$,
          # $\mue_{i,k} \leftarrow \mue_{i,k-1} - \mue t$,
          # $\mue_{i,k-1} \leftarrow t + \mue_{k,k-1} \mue_{i,k}$.
          for i in [ k+1 .. kmax ] do
            q:= mue[i][k];
            mue[i][k]:= mue[i][k-1] - mmue * q;
            mue[i][k-1]:= q + mue[k][k-1] * mue[i][k];
          od;

        fi;

        # Terminate the subalgorithm.

        if k > 2 then k:= k-1; fi;

        # Here we have always `k > r' since the loop is entered
        # for `k > r+1' only (because of `B[k-1] <> 0'),
        # so the only problem might be the case `k = r+1',
        # namely `mue[ r+1 ][r]' is used then; but this is bound
        # provided that the initial Gram matrix did not start
        # with zero columns, and its (perhaps not updated) value
        # does not matter because this would mean just to subtract
        # a multiple of a zero vector.

        RED( k-1 );

      od;

      if B[ r+1 ] = 0 then
        r:= r+1;
      fi;

      for l in [ k-2, k-3 .. r+1 ] do
        RED( l );
      od;
      k:= k+1;

    # step 4 (Finished?)
    # If $k \leq n$ go to step 2.

    od;

    # Otherwise, let $r$ be the number of initial vectors $b_i$
    # which are equal to zero,
    # take the nonzero rows and columns of the Gram matrix
    # the transformation matrix $H \in GL_n(\Z)$
    # and terminate the algorithm.

    if IsBound( arg[1][1][n] ) then

      # adjust also upper half of the Gram matrix
      gram:= gram{ [ r+1 .. n ] }{ [ r+1 .. n ] };
      for i in [ 2 .. n-r ] do
        for j in [ 1 .. i-1 ] do
          gram[j][i]:= gram[i][j];
        od;
      od;

    else

      # get the triangular matrix
      gram:= gram{ [ r+1 .. n ] };
      for i in [ 1 .. n-r ] do
        gram[i]:= gram[i]{ [ r+1 .. r+i ] };
      od;

    fi;

    Info( InfoZLattice, 1,
          "LLLReducedGramMat returns matrix of length ", n-r );

    mue:= List( [ r+1 .. n ], i -> mue[i]{ [ r+1 .. i-1 ] } );
    MakeImmutable( mue );
    B:= B{ [ r+1 .. n ] };
    MakeImmutable( B );

    return rec( remainder      := gram,
                relations      := H{ [  1  .. r ] },
                transformation := H{ [ r+1 .. n ] },
                mue            := mue,
                B              := B );
end );


#############################################################################
##
#F  ShortestVectors( <mat>, <bound> [, \"positive\" ] )
##
InstallGlobalFunction( ShortestVectors, function( arg )
    local
    # variables
          n,  checkpositiv, a, llg, nullv, m, c, anz, con, b, v,
    # procedures
          srt, vschr;

    # search for shortest vectors
    srt := function( d, dam )
    local i, j, x, k, k1, q;
    if d = 0 then
       if v = nullv then
          con := false;
       else
          anz := anz + 1;
          vschr( dam );
       fi;
    else
       x := 0;
       for j in [d+1..n] do
          x := x + v[j] * llg.mue[j][d];
       od;
       i := - Int( x );
       if AbsInt( -x-i ) * 2 > 1 then
          i := i - SignInt( x );
       fi;
       k := i + x;
       q := ( m + 1/1000 - dam ) / llg.B[d];
       if k * k < q then
          repeat
             i := i + 1;
             k := k + 1;
          until k * k >= q and k > 0;
          i := i - 1;
          k := k - 1;
          while k * k < q and con do
             v[d] := i;
             k1 := llg.B[d] * k * k + dam;
             srt( d-1, k1 );
             i := i - 1;
             k := k - 1;
          od;
       fi;
    fi;
    end;

    # output of vector
    vschr := function( dam )
    local i, j, w, neg;
    c.vectors[anz] := [];
    neg := false;
    for i in [1..n] do
       w := 0;
       for j in [1..n] do
          w := w + v[j] * llg.transformation[j][i];
       od;
       if w < 0 then
          neg := true;
#T better here check testpositiv and return!
       fi;
       c.vectors[anz][i] := w;
    od;
    if checkpositiv and neg then
       Unbind(c.vectors[anz]);
       anz := anz - 1;
    else
       c.norms[anz] := dam;
    fi;
    end;

    # main program
    # check input
    if    not IsBound( arg[1] )
       or not IsList( arg[1] ) or not IsList( arg[1][1] ) then
       Error ( "first argument must be Gram matrix\n",
          "usage: ShortestVectors( <mat>, <integer> [,<\"positive\">] )" );
    elif not IsBound( arg[2] ) or not IsInt( arg[2] ) then
       Error ( "second argument must be integer\n",
          "usage: ShortestVectors( <mat>, <integer> [,<\"positive\">] )");
    elif IsBound( arg[3] ) then
       if IsString( arg[3] ) then
          if arg[3] = "positive" then
             checkpositiv := true;
          else
             checkpositiv := false;
          fi;
       else
          Error ( "third argument must be string\n",
          "usage: ShortestVectors( <mat>, <integer> [,<\"positive\">] )");
       fi;
    else
       checkpositiv := false;
    fi;

    a := arg[1];
    m := arg[2];
    n := Length( a );
    b := List( a, ShallowCopy );
    c     := rec( vectors:= [], norms:= [] );
    v     := ListWithIdenticalEntries( n, 0 );
    nullv := ListWithIdenticalEntries( n, 0 );

    llg:= LLLReducedGramMat( b );
#T here check that the matrix is really regular
#T (empty relations component)

    anz := 0;
    con := true;
    srt( n, 0 );

    Info( InfoZLattice, 2,
          "ShortestVectors: ", Length( c.vectors ), " vectors found" );
    return c;
end );


#############################################################################
##
#F  OrthogonalEmbeddings( <grammat> [, \"positive\" ] [, <integer> ] )
##
InstallGlobalFunction( OrthogonalEmbeddings, function( arg )
    local
    # sonstige prozeduren
          Symmatinv,
    # variablen fuer Embed
          maxdim, M, D, s, phi, mult, m, x, t, x2, sumg, sumh,
          f, invg, sol, solcount, out,
          l, g, i, j, k, n, a, IdMat, chpo,
    # booleans
          checkdim,
    # prozeduren fuer Embed
          comp1, comp2, scp2, multiples, solvevDMtr,
          Dextend, Mextend, inca, rnew,
          deca;

    Symmatinv := function( b )
    # inverts symmetric matrices

    local n, i, j, l, k, c, d, ba, B, kgv1;
    n := Length( b );
    c := List( IdMat, ShallowCopy );
    d := [];
    B := [];
    kgv1 := 1;
    ba := List( IdMat, ShallowCopy );
    for i in [1..n] do
       k := b[i][i];
       for j in [1..i-1] do
          k := k - c[i][j] * c[i][j] * B[j];
       od;
       B[i] := k;
       for j in [i+1..n] do
          k := b[j][i];
          for l in [1..i-1] do
             k := k - c[i][l] * c[j][l] * B[l];
          od;
          if B[i] <> 0 then
             c[j][i] := k / B[i];
          else
             Error ( "matrix not invertible, ", Ordinal( i ),
                     " column is linearly dependent" );
          fi;
       od;
    od;
    if B[n] = 0 then
       Error ( "matrix not invertible, ", Ordinal( i ),
               " column is linearly dependent" );
    fi;
    for i in [1..n-1] do
       for j in [i+1..n] do
          if c[j][i] <> 0 then
             for l in [1..i] do
                ba[j][l] := ba[j][l] - c[j][i] * ba[i][l];
             od;
          fi;
       od;
    od;
    for i in [1..n] do
       for j in [1..i-1] do
          ba[j][i] := ba[i][j];
          ba[i][j] := ba[i][j] / B[i];
       od;
       ba[i][i] := 1/B[i];
    od;
    for i in [1..n] do
       d[i] := [];
       for j in [1..n] do
          if i >= j then
             k := ba[i][j];
             l := i + 1;
          else
             l := j;
             k := 0;
          fi;
          while l <= n do
             k := k + ba[i][l] * ba[l][j];
             l := l + 1;
          od;
          d[i][j] := k;
          kgv1 := Lcm( kgv1, DenominatorRat( k ) );
       od;
    od;
    for i in [1..n] do
       for j in [1..n] do
          d[i][j] := kgv1 * d[i][j];
       od;
    od;
    return rec( inverse := d, enuminator := kgv1 );
    end;

    # program embed

    comp1 := function( a, b )
    local i;
    if ( a[n+1] < b[n+1] ) then
      return false;
    elif ( a[n+1] > b[n+1] ) then
      return true;
    else
      for i in [ 1 .. n ] do
        if AbsInt( a[i] ) > AbsInt( b[i] ) then
          return true;
        elif AbsInt( a[i] ) < AbsInt( b[i] ) then
          return false;
        fi;
      od;
    fi;
    return false;
    end;

    comp2 := function( a, b )
    local i, t;
    t := Length(a)-1;
    if a[t+1] < b[t+1] then
      return true;
    elif a[t+1] > b[t+1] then
      return false;
    else
      for i in [ 1 .. t ] do
        if a[i] < b[i] then
          return false;
        elif a[i] > b[i] then
          return true;
        fi;
      od;
    fi;
    return false;
    end;

    scp2 := function( k, l )
    # uses    x, invg,
    # changes
    local   i, j, sum, sum1;

    sum := 0;
    for i in [1..n] do
       sum1 := 0;
       for j in [1..n] do
          sum1 := sum1 + x[k][j] * invg[j][i];
       od;
       sum := sum + sum1 * x[l][i];
    od;
    return sum;
    end;

    multiples := function( l )
    # uses    m, phi,
    # changes mult, s, k, a, sumh, sumg,
    local   v, r, i, j, brk;

    for j in [1..n] do
       sumh[j] := 0;
    od;
    i := l;
    while i <= t and ( not checkdim or s <= maxdim ) do
       if mult[i] * phi[i][i] < m then
          brk := false;
          repeat
             v := solvevDMtr( i );
             if not IsBound( v[1] ) or not IsList( v[1] ) then
                r := rnew( v, i );
                if r >= 0 then
                   if r > 0 then
                      Dextend( r );
                   fi;
                   Mextend( v, i );
                   a[i] := a[i] + 1;
                else
                   brk := true;
                fi;
             else
                brk := true;
             fi;
          until a[i] * phi[i][i] >= m or ( checkdim and s > maxdim )
                or brk;
          mult[i] := a[i];
          while a[i] > 0 do
             s := s - 1;
             if Last(M[s]) = 1 then
                k := k -1;
             fi;
             a[i] := a[i] - 1;
          od;
       fi;
       if mult[i] <> 0 then
          for j in [1..n] do
             sumh[j] := sumh[j] + mult[i] * x2[i][j];
          od;
       fi;
       i := i + 1;
    od;
    end;

    solvevDMtr := function( l )
    #  uses    M, D, phi, f,
    #  changes
    local   M1, M2, i, j, k1, v, sum;

    k1 := 1;
    v := [];
    i := 1;
    while i < s do
       sum := 0;
       M1 := Length( M[i] );
       M2 := M[i][M1];
       for j in [1..M1-1] do
          sum := sum + v[j] * M[i][j];
       od;
       if M2 = 1 then
          v[k1] := -( phi[l][f[i]] + sum ) / D[k1];
          k1 := k1 + 1;
       else
          if DenominatorRat( sum ) <> 1
          or NumeratorRat( sum ) <> -phi[l][f[i]] then
             v[1] := [];
             i := s;
          fi;
       fi;
       i := i + 1;
    od;
    return( v );
    end;

    inca := function( l )
    #  uses    x2,
    #  changes l, a, sumg, sumh,
    local   v, r, brk, i;

    while l <= t and ( not checkdim or s <= maxdim ) do
       brk := false;
       repeat
          v := solvevDMtr( l );
          if not IsBound( v[1] ) or not IsList( v[1] ) then
             r := rnew( v, l );
             if r >= 0 then
                if r > 0 then
                   Dextend( r );
                fi;
                Mextend( v, l );
                a[l] := a[l] + 1;
                for i in [1..n] do
                   sumg[i] := sumg[i] + x2[l][i];
                od;
             else
                brk := true;
             fi;
          else
             brk := true;
          fi;
       until a[l] >= mult[l] or ( checkdim and s > maxdim ) or brk;
       mult[l] := 0;
       l := l + 1;
    od;
    return l;
    end;

    rnew := function( v, l )
    #  uses    phi, m, k, D,
    #  changes v,
    local   sum, i;
    sum := m - phi[l][l];
    for i in [1..k-1] do
       sum := sum - v[i] * D[i] * v[i];
    od;
    if sum >= 0 then
     if sum > 0 then
       v[k] := 1;
     else
       v[k] := 0;
     fi;
    fi;
    return sum;
    end;

    Mextend := function( line, l )
    #  uses    D,
    #  changes M, s, f,
    local   i;
    for i in [1..Length( line )-1] do
       line[i] := line[i] * D[i];
    od;
    M[s] := line;
    f[s] := l;
    s := s + 1;
    end;

    Dextend := function( r )
    #  uses    a,
    #  changes k, D,
    D[k] := r;
    k := k + 1;
    end;

    deca := function( l )
    #  uses    x2, t, M,
    #  changes l, k, s, a, sumg,
    local   i;
if l = 0 then return l; fi;
#   if l <> 1 then
#      l := l - 1;
#      if l = t - 1 then
if l = t then
          while a[l] > 0 do
             s := s -1;
             if Last(M[s]) = 1 then
                k := k - 1;
             fi;
             a[l] := a[l] - 1;
             for i in [1..n] do
                sumg[i] := sumg[i] - x2[l][i];
             od;
          od;
#         l := deca( l );
l:= deca( t-1 );
       else
          if a[l] <> 0 then
             s := s - 1;
             if Last(M[s]) = 1 then
                k := k - 1;
             fi;
             a[l] := a[l] - 1;
             for i in [1..n] do
                sumg[i] := sumg[i] - x2[l][i];
             od;
             l := l + 1;
          else
#            l := deca( l );
l := deca( l-1 );
          fi;
       fi;
#   fi;
    return l;
    end;

    # check input
    if not IsList( arg[1] ) or not IsList( arg[1][1] ) then
       Error( "first argument must be symmetric Gram matrix\n",
              "usage : Orthog... ( < gram-matrix > \n",
              " [, <\"positive\"> ] [, < integer > ] )" );
    elif Length( arg[1] ) <> Length( arg[1][1] ) then
       Error( "Gram matrix must be quadratic\n",
              "usage : Orthog... ( < gram-matrix >\n",
              " [, <\"positive\"> ] [, < integer > ] )" );
    fi;
    g := List( arg[1], ShallowCopy );
    checkdim := false;
    chpo := "xxx";
    if IsBound( arg[2] ) then
       if IsString( arg[2] ) then
          chpo := arg[2];
       else
          if IsInt( arg[2] ) then
             maxdim := arg[2];
             checkdim := true;
          else
             Error( "second argument must be string or integer\n",
                    "usage : Orthog... ( < gram-matrix >\n",
                    " [, <\"positive\"> ] [, < integer > ] )" );
          fi;
       fi;
    fi;
    if IsBound( arg[3] ) then
       if IsString( arg[3] ) then
          chpo := arg[3];
       else
          if IsInt( arg[3] ) then
             maxdim := arg[3];
             checkdim := true;
          else
             Error( "third argument must be string or integer\n",
                    "usage : Orthog... ( < gram-matrix >\n",
                    " [, <\"positive\"> ] [, < integer > ] )" );
          fi;
       fi;
    fi;
    n := Length( g );
    for i in [1..n] do
       for j in [1..i-1] do
          if g[i][j] <> g[j][i] then
             Error( "matrix not symmetric \n",
                    "usage : Orthog... ( < gram-matrix >\n",
                    " [, <\"positive\"> ] [, < integer > ] )" );
          fi;
       od;
    od;

    # main program
    IdMat := IdentityMat( n );
    invg  := Symmatinv( g );
    m     := invg.enuminator;
    invg  := invg.inverse;
    x     := ShortestVectors( invg, m, chpo );
    t     := Length(x.vectors);
    for i in [1..t] do
       x.vectors[i][n+1] := x.norms[i];
    od;
    x    := x.vectors;
    M    := [];
    M[1] := [];
    D    := [];
    mult := [];
    sol  := [];
    f    := [];
    solcount := 0;
    s        := 1;
    k        := 1;
    l        := 1;
    a        := [];
    for i in [1..t] do
       a[i]      := 0;
       x[i][n+2] := 0;
       mult[i]   := 0;
    od;
    sumg := [];
    sumh := [];
    for i in [1..n] do
       sumg[i] := 0;
       sumh[i] := 0;
    od;
    Sort(x,comp1);
    x2 := [];
    for i in [1..t] do
       x2[i] := [];
       for j in [1..n] do
          x2[i][j]  := x[i][j] * x[i][j];
          x[i][n+2] := x[i][n+2] + x2[i][j];
       od;
    od;
    phi := [];
    for i in [1..t] do
       phi[i] := [];
       for j in [1..i-1] do
          phi[i][j] := scp2( i, j );
       od;
       phi[i][i] := x[i][n+1];
    od;
    repeat
       multiples( l );

       # (former call of `tracecond')
       if ForAll( [ 1 .. n ], i -> g[i][i] - sumg[i] <= sumh[i] ) then
          l := inca( l );
# Here we have l = t+1.
          if s-k = n then
             solcount := solcount + 1;
             Info( InfoZLattice, 2,
                   "OrthogonalEmbeddings: ", solcount, " solutions found" );
             sol[solcount] := [];
             for i in [1..t] do
                sol[solcount][i] := a[i];
             od;
             sol[solcount][t+1]  := s - 1;
          fi;
       fi;
l:= t;
       l := deca( l );
    until l <= 1;
    out := rec( vectors := [], norms := [], solutions := [] );
    for i in [1..t] do
       out.vectors[i] := [];
       out.norms[i]   := x[i][n+1]/m;
       for j in [1..n] do
          out.vectors[i][j] := x[i][j];
       od;
    od;
    Sort( sol, comp2 );
    for i in [1..solcount] do
       out.solutions[i]  := [];
       for j in [1..t] do
          for k in [1..sol[i][j]] do
            Add( out.solutions[i], j );
          od;
       od;
    od;
    return out;
end );


#############################################################################
##
#F  LLLint(<lat>) . . . . . . . . . . . . . . . . . . . .. . integer only LLL
##
InstallGlobalFunction( LLLint, function( arg )
  local lat,b,mu,i,j,k,dim,l,d,dkp,n,r,za,ne,nne,dkm,dkma,mue,muea,muk,mum,
          ca1,ca2,cb1,cb2,tw,sel,s,dkpv,cn,cd;

  lat:=arg[1];
  if Length(arg)>1 then
    cn:=arg[2];
  else
    cn:=3/4;
  fi;
  cd:=DenominatorRat(cn);
  cn:=NumeratorRat(cn);
  b:= List( lat, ShallowCopy );
  mu:=[];
  d:=[1,b[1]*b[1]];
  n:=Length(lat);
  Info( InfoZLattice, 1, "integer LLL in dimension ", n );
  dim:=1;
  k:=2;

  while dim<n do

    # mu[k][j] berechnen (Gram-Schmidt ohne Berechnung der Vektoren selbst)
    mu[k]:=[];
    for j in [1..k-1] do
      za:=d[j]*(b[k]*b[j]);
      ne:=1;
      for i in [1..j-1] do
        nne:=d[i]*d[i+1];
        za:=za*nne-d[j]*mu[k][i]*mu[j][i]*ne;
        ne:=ne*nne;
      od;
      mu[k][j]:=za/ne;
    od;

    # berechne d[k]=\prod B_j

    za:=d[k]*(b[k]*b[k]);
    ne:=1;
    for i in [1..k-1] do
      nne:=d[i]*d[i+1];
      za:=za*nne-d[k]*mu[k][i]*mu[k][i]*ne;
      ne:=ne*nne;
    od;
    d[k+1]:=za/ne;
    if d[k+1]=0 then Error("notbas");fi;

    dim:=dim+1;

    Info( InfoZLattice, 1, "computing vector ", dim );

    while k<=dim do

      #reduce(k-1);

      ne:=d[k];
      za:=mu[k][k-1];
      if za<0 then
        za:=-za;
        s:=-1;
      else
        s:=1;
      fi;
      nne:=ne/2;
      if IsInt(nne) then
        za:=za+nne;
      else
        za:=2*za+ne;
        ne:=ne*2;
      fi;
      r:=s*QuoInt(za,ne);
      if r<>0 then
        b[k]:=b[k]-r*b[k-1];
        for j in [1..k-2] do
          mu[k][j]:=mu[k][j]-r*mu[k-1][j];
        od;
        mu[k][k-1]:=mu[k][k-1]-r*d[k];
      fi;

      mue:=mu[k][k-1];
      dkp:=d[k+1]*d[k-1];
      dkpv:=dkp*4;

      if d[k]*d[k]*cn-mue*mue*cd>dkpv then

        #(2)
        Info( InfoZLattice, 2, "swap ", k-1, " <-> ", k );

        muea:=mue;
        dkm:=d[k];
        dkma:=dkm;

        ca1:=1;
        ca2:=0;
        cb1:=0;
        cb2:=1;

        # iterierter vektor-ggT
        repeat
          dkm:=(dkp+mue*mue)/dkm;
          tw:=ca1;
          ca1:=cb1;
          cb1:=tw;
          tw:=ca2;
          ca2:=cb2;
          cb2:=tw;

          ne:=dkm;
          za:=mue;
          if za<0 then
            za:=-za;
            s:=-1;
          else
            s:=1;
          fi;
          nne:=ne/2;
          if IsInt(nne) then
            za:=za+nne;
          else
            za:=2*za+ne;
            ne:=ne*2;
          fi;
          r:=s*QuoInt(za,ne);

          if r<>0 then
            cb1:=cb1-r*ca1;
            cb2:=cb2-r*ca2;
            mue:=mue-r*dkm;
          fi;
        until dkm*dkm*cn-mue*mue*cd<=dkpv;

        d[k]:=dkm;
        mu[k][k-1]:=mue;

        tw:=ca1*b[k-1]+ca2*b[k];
        b[k]:=cb1*b[k-1]+cb2*b[k];
        b[k-1]:=tw;

        if k>2 then
          sel:=[1..k-2];
          muk:=mu[k]{sel};
          mum:=mu[k-1]{sel};
          tw:=ca1*mum+ca2*muk;
          mu[k]{sel}:=cb1*mum+cb2*muk;
          mu[k-1]{sel}:=tw;
        fi;

        for j in [k+1..dim] do
          za:=ca1*dkma+ca2*muea;
          tw:=(za*mu[j][k-1]+ca2*mu[j][k]*d[k-1])/dkma;
          mu[j][k]:=(((cb1*dkma+cb2*muea)*dkm-mue*za)*mu[j][k-1]+
                     (cb2*dkm-ca2*mue)*d[k-1]*mu[j][k])/dkma/d[k-1];
          mu[j][k-1]:=tw;
        od;

        if k>2 then
          k:=k-1;
        fi;
      else
        for l in [2..k-1] do
          #reduce(k-l);

          ne:=d[k-l+1];
          za:=mu[k][k-l];
          if za<0 then
            za:=-za;
            s:=-1;
          else
            s:=1;
          fi;
          nne:=ne/2;
          if IsInt(nne) then
            za:=za+nne;
          else
            za:=2*za+ne;
            ne:=ne*2;
          fi;
          r:=s*QuoInt(za,ne);
          if r<>0 then
            b[k]:=b[k]-r*b[k-l];
            for j in [1..k-l-1] do
              mu[k][j]:=mu[k][j]-r*mu[k-l][j];
            od;
            mu[k][k-l]:=mu[k][k-l]-r*d[k-l+1];
          fi;
        od;
        k:=k+1;
      fi;
    od;

  od;
  return b;
end );