|
#############################################################################
##
#A codemisc.gi GUAVA library Reinald Baart
#A Jasper Cramwinckel
#A Erik Roijackers
#A Eric Minkes
##
## This file contains miscellaneous functions for codes
##
########################################################################
##
#F CodeWeightEnumerator( <code> )
##
## Returns a polynomial over the rationals
## with degree not greater than the length of the code.
## The coefficient of x^i equals
## the number of codewords of weight i.
##
InstallMethod(CodeWeightEnumerator, "unrestricted code", true, [IsCode], 0,
function( code )
return LaurentPolynomialByCoefficients(
ElementsFamily(FamilyObj(Rationals)),
WeightDistribution( code ), 0 );
end);
########################################################################
##
#F CodeDistanceEnumerator( <code>, <word> )
##
## Returns a polynomial over the rationals
## with degree not greater than the length of the code.
## The coefficient of x^i equals
## the number of codewords with distance i to <word>.
InstallMethod(CodeDistanceEnumerator, "unrestricted code, codeword", true,
[IsCode, IsCodeword], 0,
function( code, word )
word := Codeword( word, code );
return LaurentPolynomialByCoefficients(
ElementsFamily(FamilyObj( Rationals )),
DistancesDistribution( code, word ), 0 );
end);
########################################################################
##
#F CodeMacWilliamsTransform( <code> )
##
## Returns a polynomial with the weight
## distribution of the dual code as
## coefficients.
##
InstallMethod(CodeMacWilliamsTransform, "unrestricted code",
true, [IsCode], 0,
function( code )
local weightdist, transform, size, n, x, i, j, tmp;
n := WordLength( code );
# if dimension < n/2, or if non-linear code,
# use weightdistribution of code,
# else use weightdistribution of dual code
if not IsLinearCode( code ) or Dimension( code ) < n / 2 then
weightdist := WeightDistribution( code );
size := Size( code );
transform := List( [ 1 .. n+1 ], x -> 0 );
for j in [ 0 .. n ] do
tmp := 0;
for i in [ 0 .. n ] do
tmp := tmp + weightdist[ i+1 ] * Krawtchouk( j, i, n, 2 );
od;
transform[ j+1 ] := tmp / size;
od;
else
transform := WeightDistribution( DualCode( code ) );
fi;
return LaurentPolynomialByCoefficients(
ElementsFamily(FamilyObj( Rationals )),
transform, 0 );
end);
########################################################################
##
#F WeightVector( <vector> )
##
## Returns the number of non-zeroes in a vector.
InstallMethod(WeightVector, "method for vector", true, [IsVector], 0,
function( vector )
local pos, number, fieldzero;
number := 0;
fieldzero := Zero( Field( vector ) );
for pos in [ 1 .. Length( vector ) ] do
if vector[ pos ] <> fieldzero then
number := number + 1;
fi;
od;
return number;
end);
########################################################################
##
#F RandomVector( <len> [, <weight> [, <field> ] ] )
##
InstallMethod(RandomVector, "length, weight, field", true,
[IsInt, IsInt, IsField], 0,
function(len, wt, field)
local vec, coord, coordlist, elslist, i;
if len <= 0 then
Error( "RandomVector: length must be a positive integer" );
fi;
if wt < -1 or wt > len then
Error( "RandomVector: <weight> must be an integer in the range",
" -1 .. ", len );
fi;
vec := NullVector( len, field );
if wt > 0 then
coordlist := [ 1 .. len ];
elslist := Difference( AsSSortedList( field ), [ Zero(field) ] );
# make wt elements of the vector non-zero,
# choosing uniformly between the other field-elements
for i in [ 1 .. wt ] do
coord := Random( coordlist );
SubtractSet( coordlist, [ coord ] );
vec[ coord ] := Random( elslist );
od;
# do nothing if w = 0
elif wt = -1 then
# for each coordinate, choose uniformly from
# all field elements, including zero
elslist := AsSSortedList( field );
for i in [ 1 .. len ] do
vec[ i ] := Random( elslist );
od;
fi;
return vec;
end);
InstallOtherMethod(RandomVector, "length, weight, fieldsize", true,
[IsInt, IsInt, IsInt], 0,
function(len, wt, q)
return RandomVector(len, wt, GF(q));
end);
InstallOtherMethod(RandomVector, "length, weight", true, [IsInt, IsInt], 0,
function(len, wt)
return RandomVector(len, wt, GF(2));
end);
InstallOtherMethod(RandomVector, "length, field", true, [IsInt, IsField], 0,
function(len, field)
return RandomVector(len, -1, field);
end);
InstallOtherMethod(RandomVector, "length", true, [IsInt], 0,
function(len)
return RandomVector(len, -1, GF(2));
end);
########################################################################
##
#F IsSelfComplementaryCode( <code> )
##
## Return true if <code> is a complementary code, false otherwise.
## A code is called complementary if for every v \in <code>
## also 1 - v \in <code> (where 1 is the all-one word).
##
InstallMethod(IsSelfComplementaryCode, "method for unrestricted code",
true, [IsCode], 0,
function ( code )
local size, els, selfcompl, alloneword, newword;
if LeftActingDomain( code ) <> GF(2) then
Error("IsSelfComplementaryCode: <code> is not a binary code" );
elif IsLinearCode( code ) then
return IsSelfComplementaryCode( code );
else
els := AsSSortedList( code );
selfcompl := true;
alloneword := AllOneCodeword( WordLength( code ), GF(2) );
while Length( els ) > 0 and selfcompl = true do
newword := alloneword - els[ 1 ];
if newword <> els[ 1 ] then
if newword in els then
els := Difference( els, [ newword ] );
else
selfcompl := false;
fi;
els := Difference( els, [ els[ 1 ] ] );
fi;
od;
return selfcompl;
fi;
end);
InstallMethod(IsSelfComplementaryCode, "method for linear code",
true, [IsLinearCode], 0,
function ( code )
if LeftActingDomain( code ) <> GF(2) then
Error("IsSelfComplementaryCode: <code> is not a binary code" );
else
return( AllOneCodeword( WordLength( code ), GF(2) ) in code );
fi;
end);
########################################################################
##
#F IsAffineCode( <code> )
##
## Return true if <code> is affine, i.e. a linear code or
## a coset of a linear code, false otherwise.
##
InstallMethod(IsAffineCode, "method for unrestricted code",
true, [IsCode], 0,
function ( code )
if IsLinearCode( code ) then
return IsAffineCode( code );
elif NullWord( code ) in code then
# code cannot be a coset code of a linear code
return false;
elif not ( Size( code ) in List( [ 0 .. WordLength( code ) ],
x -> Characteristic( LeftActingDomain( code ) ) ^ x ) ) then
# the code must have a "dimension"
return false;
else
# subtract the first codeword from all codewords.
# if the resulting code is linear, then the
# original code is affine.
return IsLinearCode(
CosetCode( code, NullWord( code )
- CodewordNr( code, 1 ) ) );
fi;
end);
InstallTrueMethod(IsAffineCode, IsLinearCode);
########################################################################
##
#F IsAlmostAffineCode( <code> )
##
## Return true if <code> is almost affine, false otherwise.
## A code is called almost affine if the size of any punctured
## code is equal to q^r for some integer r, where q is the
## size of the alphabet of the code.
##
InstallMethod(IsAlmostAffineCode, "method for unrestricted code",
true, [IsCode], 0,
function( code )
local F, n, i, j, subcode, sizelist, coordlist, almostaffine;
if IsAffineCode( code ) then
# every affine code is also almost affine
almostaffine := true;
else
# not affine
almostaffine := true;
F := LeftActingDomain( code );
# however, any code over GF(2) or GF(3) is affine
# if it is almost affine.
# so non-affine codes with q=2,3 are also not almost affine
if Size( F ) = 2
or Size( F ) = 3 then
almostaffine := false;
fi;
n := WordLength( code );
sizelist := List( [ 0 .. n ], x -> Characteristic( F ) ^ x );
# first check whether the code itself is of size q^r
if not ( Size( code ) in sizelist ) then
almostaffine := false;
else
# now check for all possible puncturings
i := 1;
while almostaffine and i < n do
coordlist := List( Tuples( [ 1 .. n ], i ),
x -> Difference( [ 1 .. n ], x ) );
j := 1;
while almostaffine
and j < Length( coordlist ) do
subcode := PuncturedCode( code, coordlist[ j ] );
# one fault is enough !
if not Size( subcode ) in sizelist then
almostaffine := false;
fi;
j := j + 1;
od;
i := i + 1;
od;
fi;
fi;
return almostaffine;
end);
InstallTrueMethod(IsAlmostAffineCode, IsAffineCode);
########################################################################
##
#F IsGriesmerCode( <code> )
##
## Return true if <code> is a Griesmer code, i.e. if
## n = \sum_{i=0}^{k-1} d/(q^i), false otherwise.
##
InstallMethod(IsGriesmerCode, "method for unrestricted code",
true, [IsCode], 0,
function( code )
if IsLinearCode( code ) then
return IsGriesmerCode( code );
else
Error( "IsGriesmerCode: <code> must be a linear code" );
fi;
end);
InstallMethod(IsGriesmerCode, "method for linear code", true,
[IsLinearCode], 0,
function( code )
local n, k, d, q;
n := WordLength( code );
k := Dimension( code );
d := MinimumDistance( code );
q := Size( LeftActingDomain( code ) );
return n = Sum( [ 0 .. k-1 ], x -> IntCeiling( d / q^x ) );
end);
########################################################################
##
#F CodeDensity( <code> )
##
## Return the density of <code>, i.e. M*V_q(n,r)/(q^n).
##
InstallMethod(CodeDensity, "method for unrestricted code", true,
[IsCode], 0,
function ( code )
local n, q, cr;
cr := CoveringRadius( code );
# Linear codes with redundancy >= 20 can return an interval
# for the Covering Radius, so this test is necessary.
if not IsInt( cr ) then
Error( "CodeDensity: the covering radius of <code> is unknown" );
fi;
n := WordLength( code );
q := Size( LeftActingDomain( code ) );
return Size( code )
* SphereContent( n, CoveringRadius( code ), q )
/ q^n;
end);
########################################################################
##
#F DecreaseMinimumDistanceUpperBound( <C>, <s>, <iteration> )
##
## Tries to compute the minimum distance of C.
## The algorithm is Leon's, see for more
## information his article.
InstallMethod(DecreaseMinimumDistanceUpperBound,
"method for unrestricted code, s, iteration",
true, [IsCode, IsInt, IsInt], 0,
function(C, s, iteration)
if IsLinearCode(C) then
return DecreaseMinimumDistanceUpperBound(C, s, iteration);
else
Error("DecreaseMinimumDIstanceUB: <C> must be a linear code");
fi;
end);
InstallMethod(DecreaseMinimumDistanceUpperBound,
"method for linear code, s, iteration",
true, [IsLinearCode, IsInt, IsInt], 0,
function ( C, s, iteration )
# <C> is the code to compute the min. dist. for
# <s> is the parameter to help find words with
# small weight
# <iteration> is number of iterations to perform
local
trials, # the number of trials so far
n, k, # some parameters of the code C
genmat, # the generator matrix of C
d, # the minimum distance so far
cont, # have we computed enough trials ?
N, # the set { 1, ..., n }
S, # a random s-subset of N
h, i, j, # some counters
sigma, # permutation, mapping of N, mapping S on {1,...,s}
tau, # permutation, for eliminating first s columns of Emat
Emat, # genmat ^ sigma
Dmat, # (k-e,n-s) right lower submatrix of Emat
e, # rank of k * s submatrix of Emat
nullrow, # row of zeroes, for appending to Emat
res, # result from PutSemiStandardForm
w, # runs through all words spanned by Dmat
t, # weight of the current codeword
v, # word with current lowest weight
Bmat, # (e, n-s) right upper submatrix of Emat
Bsupp, # supports of differences of rows of Bmat
Bweight, # weights of rows of Bmat
sup1, sup2,# temporary variables holding supports
Znonempty, # true if e < s, false otherwise (indicates whether
# Zmat is a real matrix or not
Zmat, # ( s-e, e ) middle upper submatrix of Emat
Zweight, # weights of differences of rows of Zmat
wsupp, # weight of the current codeword w of D
ij1, # 0: i<>1 and j<>1 1: i=1 xor j=1 2: i=1 and j=1
nullw, # nullword of length s, begin of w
PutSemiStandardForm, # local function for partial Gaussian
# elimination
sups, # the supports of the elements of B
found; # becomes true if a better minimum distance is
# found
# check the arguments
if s < 1 or s > Dimension( C ) then
Error( "DecreaseMinimumDistanceUB: <s> must lie between 1 and the ",
"dimension of <C>." );
fi;
if iteration < 1 then
Error( "DecreaseMinimumDistanceLB: <iteration> must be at least zero." );
fi;
# the function PutSemiStandardForm is local
###########################################################################
##
#F PutSemiStandardForm( <mat>, <s> )
##
## Put first s coordinates of mat in standard form.
## Return e as the rank of the s x s left upper
## matrix. The coordinates s+1, ..., n are not permuted.
##
## This function is based on PutStandardForm.
##
## (maybe it's better to make this function local
## in DecreaseMinimumDistanceUpperBound)
##
PutSemiStandardForm := function ( mat, s )
local k, n, zero,
stop, found,
g, h, i, j,
row, e, tau;
k := Length(mat); # number of rows: dimension
n := Length(mat[1]); # number of columns: wordlength
zero := Zero(GF(2));
stop := false;
e := 0;
tau := ( );
for j in [ 1..s ] do
if not stop then
if mat[j][j] = zero then
# start looking for another pivot
i := j;
found := false;
while ( i <= s ) and not found do
h := j;
while ( h <= k ) and not found do
if mat[h][i] <> zero then
found := true;
else
h := h + 1;
fi; # if mat[h][i] <> zero
od; # while ( h <= k ) and not found
if not found then
i := i + 1;
fi; # if not found
od; # while ( i <= s ) and not found
if not found then
stop := true;
else
# pivot found at position (h,i)
# increase subrank
e := e + 1;
# permutate the matrix so that (h,i) <-> (j,j)
if h <> j then
row := mat[h];
mat[h] := mat[j];
mat[j] := row;
fi; # if h <> j
if i <> j then
tau := tau * (i,j);
for g in [ 1 .. k ] do
mat[g] := Permuted( mat[g], (i,j) );
od; # for g in [ 1..k ]
fi; # if i <> j
fi; # if not found
else
e := e + 1;
fi; # if mat[j][j] = zero
if not stop then
for i in [ 1..k ] do
if i <> j then
if mat[i][j] <> zero then
mat[i] := mat[i] + mat[j];
fi; # if mat[i][j] <> zero
fi; # if i <> j
od; # for i in [ 1..k ]
fi; # if not stop
fi; # if not stop
od; # for j in [ 1..s ] do
return [ e, tau ];
end;
n := WordLength( C );
k := Dimension( C );
genmat := GeneratorMat( C );
# step 1. initialisation
trials := 0;
d := n;
cont := true;
found := false;
while cont do
# step 2.
trials := trials + 1;
InfoMinimumDistance( "Trial nr. ", trials, " distance: ", d, "\n" );
# step 3. choose a random s-elements subset of N
N := [ 1 .. WordLength( C ) ];
S := [ ];
for i in [ 1 .. s ] do
S[ i ] := Random( N ); # pick a random element from N
RemoveSet( N, S[ i ] ); # and remove it from N
od;
Sort( S ); # not really necessary, but
# it doesn't hurt either
# step 4. choose a permutation sigma of N,
# mapping S onto { 1, ..., s }
Append( S, N );
sigma := PermList( S ) ^ (-1);
# step 5. Emat := genmat^sigma (genmat is the generator matrix C)
Emat := [ ];
for i in [ 1 .. k ] do
Emat[ i ] := Permuted( genmat[ i ], sigma );
od;
# step 6. apply elementary row operations to E
# and perhaps a permutation tau so that
# we get the following form:
# [ I | Z | B ]
# [ 0 | 0 | D ]
# where I is the e * e identity matrix,
# e is the rank of the k * s left submatrix of E
# the permutation tau leaves { s+1, ..., n } fixed
InfoMinimumDistance( "Gaussian elimination of E ... \n");
res := PutSemiStandardForm( Emat, s );
e := res[ 1 ]; # rank (in most cases equal to s)
tau := res[ 2 ]; # permutation of { 1, ..., s }
# append null-row to Emat (at front)
nullrow := NullMat( 1, n, GF(2) );
Append( nullrow, Emat );
Emat := nullrow;
InfoMinimumDistance( "Gaussian elimination of E ... done. \n" );
# retrieve Dmat from Emat
Dmat := [ ];
for i in [ e + 1 .. k ] do
Dmat[ i - e ] := List( [ s+1 .. n ], x -> Emat[ i+1 ][ x ] );
od;
# retrieve Bmat from Emat
# we only need the support of the differences of the
# rows of B
Bmat := [ ];
Bmat[ 1 ] := NullVector( n-s, GF(2) );
for j in [ 2 .. e+1 ] do
Bmat[ j ] := List( [ s+1 .. n ], x -> Emat[ j ][ x ] );
od;
InfoMinimumDistance( "Computing supports of B ... \n" );
sups := List( [ 1 .. e+1 ], x -> Support( Codeword( Bmat[ x ] ) ) );
# compute supports of differences of rows of Bmat
# and the weights of these supports
# do this once every trial, instead of for each codeword,
# to save time
Bsupp := List( [ 1 .. e ], x -> [ ] );
Bweight := List( [ 1 .. e ], x -> [ ] );
for i in [ 1 .. e ] do
sup1 := sups[ i ];
# Bsupp[ i ] := List( [ i + 1 - KroneckerDelta( i, 1 ) .. e+1 ],
# x -> Difference( Union( sup1, sups[ x ] ),
# Intersection( sup1, sups[ x ] ) ) );
for j in [ i + 1 - KroneckerDelta( i, 1 ) .. e+1 ] do
sup2 := sups[ j ];
Bsupp[ i ][ j ] := Difference( Union( sup1, sup2 ),
Intersection( sup1, sup2 ) );
Bweight[ i ][ j ] := Length( Bsupp[ i ][ j ] );
od;
od;
InfoMinimumDistance( "Computing supports of B ... done. \n" );
# retrieve Zmat from Emat
# in this case we only need the weights of the supports of
# the differences of the rows of Zmat
# because we don't have to add them to codewords
if e < s then
InfoMinimumDistance( "Computing weights of Z ... \n" );
Znonempty := true;
Zmat := List( [ 1 .. e ], x -> [ ] );
Zmat[ 1 ] := NullVector( s-e, GF(2) );
for i in [ 2 .. e+1 ] do
Zmat[ i ] := List( [ e+1 .. s ], x -> Emat[ i ][ x ] );
od;
Zweight := List( [ 1 ..e ], x -> [ ] );
for i in [ 1 .. e ] do
for j in [ i + 1 - KroneckerDelta( i, 1 ) .. e+1 ] do
Zweight[ i ][ j ] :=
WeightCodeword( Codeword( Zmat[ i ] + Zmat[ j ] ) );
od;
od;
InfoMinimumDistance( "Computing weights of Z ... done. \n" );
else
Znonempty := false;
fi;
# step 7. for each w in (n-s, k-e) code spanned by D
for w in AsSSortedList( GeneratorMatCode( Dmat, GF(2) ) ) do
wsupp := Support( w );
# step 8.
for i in [ 1 .. e ] do
# step 9.
for j in [ i + 1 - KroneckerDelta( i, 1 ) .. e+1 ] do
ij1 := KroneckerDelta( i, 1 ) + KroneckerDelta( j, 1 );
# step 10.
if Znonempty then
t := Zweight[ i ][ j ];
else
t := 0;
fi;
# step 11.
if t <= ij1 then
# step 12.
t := t
+ Bweight[ i ][ j ]
+ Length( wsupp )
- 2 * Length( Intersection(
Bsupp[ i ][ j ], wsupp ) );
t := t + ( 2 - ij1 );
if 0 < t and t < d then
found := true;
# step 13.
d := t;
C!.upperBoundMinimumDistance :=
Minimum( UpperBoundMinimumDistance( C ), t );
# step 14.
nullw := NullVector( s, GF(2) );
Append( nullw, VectorCodeword( w ) );
v := Emat[ i ] + Emat[ j ] + nullw;
v := Permuted( v, tau ^ (-1) );
v := Permuted( v, sigma ^ (-1) );
fi;
fi;
od;
od;
od;
if iteration <= trials then
cont := false;
fi;
od;
if found then
return rec( mindist := d, word := v );
fi;
end);
[ Dauer der Verarbeitung: 0.37 Sekunden
(vorverarbeitet)
]
|
