############################################################################# ## #A codemisc.gi GUAVA library Reinald Baart #A Jasper Cramwinckel #A Erik Roijackers #A Eric Minkes ## ## This file contains miscellaneous functions for codes ## #H @(#)$Id: codemisc.gi,v 1.6 2003/02/12 03:49:16 gap Exp $ ## Revision.("guava/lib/codemisc_gi") := "@(#)$Id: codemisc.gi,v 1.6 2003/02/12 03:49:16 gap Exp $"; ######################################################################## ## #F CodeWeightEnumerator( ) ## ## 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( , ) ## ## 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 . 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( ) ## ## 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( ) ## ## 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( [, [, ] ] ) ## 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: 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( ) ## ## Return true if is a complementary code, false otherwise. ## A code is called complementary if for every v \in ## also 1 - v \in (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: 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: is not a binary code" ); else return( AllOneCodeword( WordLength( code ), GF(2) ) in code ); fi; end); ######################################################################## ## #F IsAffineCode( ) ## ## Return true if 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( ) ## ## Return true if 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( ) ## ## Return true if 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: 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( ) ## ## Return the density of , 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 is unknown" ); fi; n := WordLength( code ); q := Size( LeftActingDomain( code ) ); return Size( code ) * SphereContent( n, CoveringRadius( code ), q ) / q^n; end); ######################################################################## ## #F DecreaseMinimumDistanceUpperBound( , , ) ## ## 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: must be a linear code"); fi; end); InstallMethod(DecreaseMinimumDistanceUpperBound, "method for linear code, s, iteration", true, [IsLinearCode, IsInt, IsInt], 0, function ( C, s, iteration ) # is the code to compute the min. dist. for # is the parameter to help find words with # small weight # 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: must lie between 1 and the ", "dimension of ." ); fi; if iteration < 1 then Error( "DecreaseMinimumDistanceLB: must be at least zero." ); fi; # the function PutSemiStandardForm is local ########################################################################### ## #F PutSemiStandardForm( , ) ## ## 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);