SSL zlattice.gi Sprache: unbekannt

#############################################################################
##
## 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>, , <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>, , <depth> ) . . . . . . . . integral solutions
##
## returns the decomposition matrix <X> with `<X> * <A> = ' for <A> and
## 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( " 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>, , <depth> ) . . . . . . . . . . integral solutions
#F Decomposition( <A>, , \"nonnegative\" ) . . . . . . integral solutions
##
## For a matrix <A> of cyclotomics and a list of cyclotomic vectors,
## `Decomposition' tries to find integral solutions of the linear equation
## systems `<x> \* <A> = [i]'.
##
## <A> must have full rank, i.e., there must be a linear independent set of
## columns of same length as <A>.
##
## `Decomposition( <A>, , <depth> )', where <depth> is a nonnegative
## integer, computes for every `[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> = [i]' *is* a solution, we have `<X>[i] = <x>', otherwise
## `<X>[i] = false'.
##
## `Decomposition( <A>, , \"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> = [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>,,<depth> ) for integer <depth>\n",
 " resp. Decomposition( <A>,,\"nonnegative\" )\n",
 " ( solution <X> of <X> * <A> = )" );
 elif not ( IsMatrix(A) and IsMatrix(B)
 and Length(B[1]) = Length(A[1]) ) then
 Error( "<A>, 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 $'[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 );

SSL zlattice.gi Sprache: unbekannt

[ zur Elbe Produktseite wechseln0.56Quellennavigators Analyse erneut starten ]