There are also coding theory packages in GAP and MAGMA. See, for example, GAP3 and the Higman-Sims group (which uses the GAP coding theory package guava) and linear codes in MAGMA.

> restart;

> with(linalg): with(padic):

Warning, new definition for fibonacci

> read(`d:/maplestuff/codes/codes.mpl`);

Hamming codes

Parity check matrices of the binary Hamming code of length 2^2-1=3:

> hamming_check_binary(2);

> H:=hamming_check(3,5);
hamming_check_nonstandard(3,5);
Parity check matrices (in standard and non-standard form) of the 5-ary Hamming code of length (5^3-1)/(5-1)=31

> H:=hamming_check_binary(3);
Parity check matrices of the binary Hamming code of length 2^3-1=7:

> H:=hamming_check(3,2);

A generator matrix of the binary Hamming code of length 3:

> hamming_generator(2,2);

> G1:=hamming_generator(3,2);
A generator matrix of the binary Hamming code of length 7 (sometimes the transpose of this matrix is called the generator matrix):

The dual code of this Hamming (7,4) code is a simplex (7,3) code.

> dual_code_list(G1,2);

>

> C1:=code_list(G1,2):
C1:=sort(C1,weight_order);
We can print out the list of all codewords, sorted according to weight:

> weight_enumerator_vector(G1,2);
There are 7 code words of weight 3, 7 code words of weight 4, and no code words of weights 1,2, 5 or 6.

> weight_enumeration_polynomial(G1,2,x,y);
This is the weight enumerator polynomial of the Hamming (7,4,3)-code.

> min_distance(G1,2);
The minimum distance of the binary Hamming code of length 3

> v1:=[1,1,1,0,0,0,0];
matrix_times_vector_modp(H,v1,2);
This received word v1 has an error since H*v1 is non-zero.

> decode_hamming_binary(v1);
This command only works for binary Hamming codes.

> decode(v1,G1,2,3);
This command works for any linear code (2 is the prime residue, 3 is the distance of the code).

>

> H:=hamming_check(3,3);
G:=hamming_generator(3,3);
size_of_code:=3^(10);
parity check and generator matrix of Hamming ((p^k-1)/(p-1),(p^k-1)/(p-1)-k,3)-code with k=3,p=3

> v1:=convert(row(H,1),list)+[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
matrix_times_vector_modp(H,v1,3);
This received word v1 has an error since H*v1 is non-zero

> decode_hamming(v1,3);

Reed-Muller codes

> G:=RM_generator(3);

This is a generator matrix of the binary 1st order Reed-Muller code of length 8. The matrix for the RM code of length 16 is:

> G_RM:=RM_generator(4);

> C:=code_list(G_RM,2):
C_RM:=sort(C,weight_order);

> nops(C);

Cyclic codes

One can test if g(x) is the generating polynomial of a cyclic code of length n over F_p using iscyclic(g(x),n,p)

> iscyclic(x+1,7,2);

> iscyclic(x^2+1,7,2);

>

> iscyclic(x^3+x+1,7,2);

Enter the generator polynomial of cyclic code and generator_matrix will return the corresponding generator matrix.

> G:=generator_matrix(x+1,7,2);

> min_distance(G,2);

The following generator matrix is the same size as that of the above Hamming code. Do they yield the same code?

> G2:=generator_matrix(x^3+x^2+1,7,2);

> C2:=code_list(G2,2):
C2:=sort(C2,weight_order);

> evalb(C1=C2);
This cyclic code generated by x^3+x^2+1 is not the Hamming code of length 7 and distance 3, over F_2, given above.

> min_distance(G2,2);

There is some built-in error checking in generator_matrix :

> generator_matrix(x^2+1,7,2);

Error, (in generator_matrix) the polynomial doesn't generate a cyclic code

Enter the generator polynomial of cyclic code and check_matrix will return the corresponding parity check matrix.

> H:=check_matrix(x^3+x^2+1,7,2);

> v1:=[1,1,1,1,0,0,0];
matrix_times_vector_modp(H,v1,2);
The word v1 is not in the cyclic code generated by x^3+x^2+1 since H*v1 is non-zero.

> decode(v1,G2,2,2);

> check_matrix(x^3+y^2+1,7,2);
check_matrix(x^3+1,7,2);
check_matrix has some built-in error checking

Error, (in check_matrix) g should be a polynomial in x

Error, (in check_matrix) the polynomial doesn't generate a cyclic code

>

Golay codes

The generator matrix for the ternary (12,6) code:

> G:=golay12();

> min_distance(G,3);
This is very time-consuming on some machines (about 15 seconds on a 300Mhz pentium).

The parity check matrix for the ternary (12,6) code:

> H:=golay_check12();

> v1:=[1, 2, 2, 2, 2, 1, 0, 1, 0, 0, 1, 2];
matrix_times_vector_modp(H,v1,3);
v1 is not a codeword in the (12,6)-Golay code. It has only 1 error (in the first position).

> decode(v1,G,3,5);
The decoding:

> v2:=[1, 1, 2, 2, 2, 1, 0, 1, 0, 0, 1, 2];
decode(v2,G,3,5);
Another example, with 2 errors.

>

The generator matrix for the perfect ternary (11,6) code:

> G:=golay11();

> min_distance(G,3);
This is very time-consuming.

The parity check matrix for the ternary (11,6) code:

> H:=golay_check11();

The generator matrix for the perfect binary (23,12) code:

> G:=golay23();

The parity check matrix for the perfect binary (23,12) code:

> H:=golay_check23();

> v1:=[0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1];
matrix_times_vector_modp(H,v1,2);
v1 is not a codeword in the (23,12)-Golay code. It has only 1 error (in the first position).

> decode(v1,G,2,7);
This is time-consuming to decode than the above Hamming code since the Golay code is larger.

> v2:=[0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0];
matrix_times_vector_modp(H,v2,2);
decode(v2,G,2,7);
This "string" has 3 errors

> G:=golay24();
The generator matrix for the (24,12)-Golay code

> H:=golay_check24();
The parity check matrix for the (24,12)-Golay code.

> coldim(G);coldim(H);

BCH codes

We shall construct by hand some BCH codes. A BCH code is a cyclic code which included the binary Hamming codes as a special case.

To construct a (binary, narrow) BCH code, you must have a primitive element alpha of degree k over F=GF(2). You must know the minimal polynomials of a sequence of powers of alpha: alpha, alpha^2, ..., alpha^(2t). Take the lcm of these. This new polynomial is the generating polynomial of the BCH code. It has minimum distance 2t+1.

> ### WARNING: persistent store makes one-argument readlib obsolete
readlib(lattice):
minpoly_lattice:=minpoly;
#linalg has a minpoly command too
with(numtheory):
### WARNING: persistent store makes one-argument readlib obsolete
readlib(GF):

Warning, new definition for order

>

> f:=x^4+x+1;
alias(beta=RootOf(f)):
Factor(f) mod 2;
Factor(f,beta) mod 2;
Roots(f,beta) mod 2;
for i from 1 to 19 do
pp[i]:=unapply(minpoly_lattice(beta^i,4),_X):
print(beta^i,pp[i](x));
od:

> GenPoly:=x->Lcm(seq(pp[j](x),j=1..4)) mod 2:
GenPoly(x);
Quo(x^(15)-1,GenPoly(x),x) mod 2;
iscyclic(GenPoly(x),15,2);

In other words, GenPoly(x) is indeed the generating polynomial for some cyclic code over F_2.

> G:=generator_matrix(GenPoly(x),15,2);
linalg[rowdim](G);
linalg[coldim](G);
H:=check_matrix(GenPoly(x),15,2);
linalg[rowdim](H);
linalg[coldim](H);
This code is 2-error correcting.

>

> f:=x^6+x^3+1;
alias(beta=RootOf(f)):
Irreduc(f) mod 2;
Factor(f,beta) mod 2;
Roots(f,beta) mod 2;
for i from 1 to 12 do
pp[i]:=unapply(minpoly_lattice(beta^i,6),_X):
print(beta^i,pp[i](x));
od:

> GenPoly:=x->Lcm(seq(pp[j](x),j=1..8)) mod 2:
GenPoly(x);
iscyclic(GenPoly(x),63,2);

In other words, GenPoly(x) is indeed the generating polynomial for some cyclic code over F_2. This code is 4-error correcting.

> G:=generator_matrix(GenPoly(x),63,2):
linalg[rowdim](G);
linalg[coldim](G);
H:=check_matrix(GenPoly(x),63,2);

>

When are linear codes equivalent?

We must load the group21_v5.mpl package

(see http://web.usna.navy.mil/~wdj/symm_gp.htm)

> with(combinat):
with(group):
read(`d:/maplestuff/group/group21_v5.mpl`);

Warning, new definition for fibonacci

>

> G1:=hamming_generator(3,2);
C1:=code_list(G1,2):
C1:=sort(C1,weight_lexorder);
Gausselim(G1) mod 2;
G2:=generator_matrix(x^3+x^2+1,7,2);
C2:=code_list(G2,2):
C2:=sort(C2,weight_lexorder);

The following program (on a 350Mz pentium II) returns true , if the binary linear code C1 with generator matrix G1 is equivalent to the binary linear code C2 with gen matrix G2 plus an associated column permutation for the second matrix, and returns false otherwise.

> t0:=time():
isequivalent(G1,G2);
time()-t0;

The following program does the same thing as the program above but uses a different algorithm. (Usually this second algorithm seems to be faster but obviously not in this case.)

> t0:=time():
isequivalent2(G1,G2);
time()-t0;

Next, we compute the automorphism group of a binary linear code.

> G1:=matrix([[1,1,0,1,0],[1,1,1,0,1],[1,1,0,1,0]]);

> t0:=time():
automorphismgroup(G1);
time()-t0;

>