source:sage/coding/linear_code.py@7880:a41a364cde98

Revision 7880:a41a364cde98, 70.0 KB checked in by William Stein <wstein@…>, 5 years ago (diff)

2.9.1

Line
1r"""
2Linear Codes
3
4VERSION: 0.9
5
6Let $F$ be a finite field (we denote the finite field with $q$ elements
7$GF(q)$ by $\FF_q$). A subspace of $F^n$ (with the standard basis)
8is called a {\it linear code} of length $n$. If its
9dimension is denoted $k$ then we typically store a basis
10of $C$ as a $k\times n$ matrix (the rows are the basis vectors)
11called the {\it generator matrix} of $C$.
12The rows of the {\it parity check matrix} of $C$ are a basis
13for the code,
14
15$C^* = \{ v \in GF(q)^n\ |\ v\cdot c = 0,\ for \ all\ c \in C \}, 16$
17called the {\it dual space} of  $C$.
18
19If $F=\FF_2$ then $C$ is called a {\it binary code}.
20If $F = \FF_q$ then $C$ is called a {\it $q$-ary code}.
21The elements of a code $C$ are called {\it codewords}.
22
23Let $F$ be a finite field with $q$ elements. Here's a constructive
24definition of a cyclic code of length $n$.
25
26\begin{enumerate}
27\item
28      Pick a monic polynomial $g(x)\in F[x]$ dividing $x^n-1$.
29      This is called the {\it generating polynomial} of the code.
30\item
31      For each polynomial $p(x)\in F[x]$, compute
32      $p(x)g(x)\ ({\rm mod}\ x^n-1).$
33      Denote the answer by $c_0+c_1x+...+c_{n-1}x^{n-1}$.
34\item
35      ${\bf c} =(c_0,c_1,...,c_{n-1})$ is a codeword in $C$. Every
36      codeword in $C$ arises in this way (from some $p(x)$).
37\end{enumerate}
38The {\it polynomial notation} for the code is to call
39$c_0+c_1x+...+c_{n-1}x^{n-1}$ the codeword (instead of
40$(c_0,c_1,...,c_{n-1})$). The polynomial $h(x)=(x^n-1)/g(x)$
41is called the {\it check polynomial} of $C$.
42
43Let $n$ be a positive integer relatively prime to $q$ and
44let $\alpha$ be a primitive $n$-th root of unity. Each generator
45polynomial $g$ of a cyclic code $C$ of length $n$ has a factorization
46of the form
47
48$49g(x) = (x - \alpha^{k_1})...(x - \alpha^{k_r}), 50$
51where $\{k_1,...,k_r\} \subset \{0,...,n-1\}$. The numbers
52$\alpha^{k_i}$, $1 \leq i \leq r$, are called the {\it zeros}
53of the code $C$. Many families of cyclic codes (such as
54the quadratic residue codes) are defined using properties of the
55zeros of $C$.
56
57The symmetric group $S_n$ acts on $F^n$ by permuting coordinates.
58If an element $p\in S_n$ sends a code $C$ of length $n$ to itself
59(in other words, every codeword of $C$ is sent to some other codeword
60of $C$) then $p$ is called a {\it permutation automorphism} of $C$.
61The (permutation) automorphism group is denoted $Aut(C)$.
62
63This file contains
64\begin{enumerate}
65\item
66LinearCode, Codeword class definitions; LinearCode_from_vectorspace
67conversion function
68\item
69The spectrum (weight distribution), minimum distance
70programs (calling Steve Linton's C programs), characteristic function,
71and several implementations of the Duursma zeta function.
72\item
73interface with best_known_linear_code (interface with A. Brouwer's online tables has
74been disabled), bounds_minimum_distance which call tables
75in GUAVA (updated May 2006) created by Cen Tjhai instead of the online
76internet tables.
77\item
78ToricCode, TrivialCode (in a separate "guava.py" module, you will find the constructions
79Hamming codes, "random" linear codes, Golay codes, binary Reed-Muller codes,
81wrapped form the corresponding GUAVA codes),
82\item
83gen_mat, list, check_mat, decode, dual_code, extended_code, genus, binomial_moment,
84and divisor methods for LinearCode,
85\item
86permutation methods:
87is_permutation_automorphism, permutation_automorphism_group,
88permuted_code, standard_form, module_composition_factors.
89\item
90design-theoretic methods:
91assmus_mattson_designs (implementing Assmus-Mattson Theorem).
92\end{enumerate}
93
94
95    EXAMPLES:
96        sage: MS = MatrixSpace(GF(2),4,7)
97        sage: G = MS([[1,1,1,0,0,0,0], [1,0,0,1,1,0,0], [0,1,0,1,0,1,0], [1,1,0,1,0,0,1]])
98        sage: C = LinearCode(G)
99        sage: C.basis()
100        [(1, 1, 1, 0, 0, 0, 0),
101         (1, 0, 0, 1, 1, 0, 0),
102         (0, 1, 0, 1, 0, 1, 0),
103         (1, 1, 0, 1, 0, 0, 1)]
104        sage: c = C.basis()[1]
105        sage: c in C
106        True
107        sage: c.nonzero_positions()
108        [0, 3, 4]
109        sage: c.support()
110        [0, 3, 4]
111        sage: c.parent()
112        Vector space of dimension 7 over Finite Field of size 2
113
115\begin{enumerate}
116\item More wrappers
117\item GRS codes and special decoders.
118\item $P^1$ Goppa codes and group actions on $P^1$ RR space codes.
119\end{enumerate}
120
121REFERENCES:
122   [HP] W. C. Huffman and V. Pless, {\bf Fundamentals of error-correcting codes},
123        Cambridge Univ. Press, 2003.
124
125   [Gu] GUAVA manual, http://www.gap-system.org/Packages/guava.html
126
127AUTHOR:
128    -- David Joyner (2005-11-22, 2006-12-03): initial version
129    -- William Stein (2006-01-23) -- Inclusion in SAGE
130    -- David Joyner (2006-01-30, 2006-04): small fixes
131    -- DJ (2006-07): added documentation, group-theoretical methods,
133    -- DJ (2006-08): hopeful latex fixes to documention, added list
134and __iter__ methods to LinearCode and examples, added hamming_weight
135function, fixed random method to return a vector, TrivialCode,
136fixed subtle bug in dual_code, added galois_closure method,
137fixed mysterious bug in permutation_automorphism_group (GAP
138was over-using "G" somehow?)
139    -- DJ (2006-08): hopeful latex fixes to documention,
141assmus_mattson_designs (implementing Assmus-Mattson Theorem).
142    -- DJ (2006-09): modified decode syntax, fixed bug in is_galois_closed,
144    -- Nick Alexander (2006-12-10): factor GUAVA code to guava.py
145    -- DJ (2007-05): added methods punctured, shortened, divisor,
146characteristic_polynomial, binomial_moment, support for LinearCode.
147Completely rewritten zeta_function (old version is now zeta_function2)
148and a new function, LinearCodeFromVectorSpace.
149    -- DJ (2007-11): added zeta_polynomial, weight_enumerator,
150                     chinen_polynomial;  improved best_known_code;
152                     added is_equivalent (for binary codes)
153
154TESTS:
155   sage: MS = MatrixSpace(GF(2),4,7)
156   sage: G  = MS([[1,1,1,0,0,0,0], [1,0,0,1,1,0,0], [0,1,0,1,0,1,0], [1,1,0,1,0,0,1]])
157   sage: C  = LinearCode(G)
159   True
160"""
161
162#*****************************************************************************
163#       Copyright (C) 2005 David Joyner <wdjoyner@gmail.com>
164#                     2006 William Stein <wstein@gmail.com>
165#
167#  version 2 or later (at your preference).
168#
170#*****************************************************************************
171
172import copy
173import sage.modules.free_module as fm
174import sage.modules.module as module
175import sage.modules.free_module_element as fme
176from sage.interfaces.all import gap
177from sage.rings.finite_field import GF
178from sage.groups.perm_gps.permgroup import PermutationGroup
179from sage.matrix.matrix_space import MatrixSpace
180from sage.rings.arith import GCD, rising_factorial, binomial
181from sage.groups.all import SymmetricGroup
182from sage.misc.sage_eval import sage_eval
184from sage.misc.functional import log, is_even, is_odd
185from sage.rings.rational_field import QQ
186from sage.structure.parent_gens import ParentWithGens
187from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
188from sage.rings.fraction_field import FractionField
189from sage.rings.integer_ring import IntegerRing
190from sage.combinat.set_partition import SetPartitions
191from sage.rings.real_mpfr import RR      ## RealField
192
193ZZ = IntegerRing()
194VectorSpace = fm.VectorSpace
195
196###################### coding theory functions ##############################
197
198def hamming_weight(v):
199    return len(v.nonzero_positions())
200
201def wtdist(Gmat, F):
202    r"""
203    INPUT:
204        Gmat -- a string representing a GAP generator matrix G of a
205                linear code.
206        F -- a (SAGE) finite field - the base field of the code.
207
208    OUTPUT:
209        Returns the spectrum of the associated code.
210
211    EXAMPLES:
212        sage: Gstr = 'Z(2)*[[1,1,1,0,0,0,0], [1,0,0,1,1,0,0], [0,1,0,1,0,1,0], [1,1,0,1,0,0,1]]'
213        sage: F = GF(2)
214        sage: sage.coding.linear_code.wtdist(Gstr, F)
215        [1, 0, 0, 7, 7, 0, 0, 1]
216
217    Here Gstr is a generator matrix of the Hamming [7,4,3] binary code.
218
219    ALGORITHM: Uses C programs written by Steve Linton in the kernel
220        of GAP, so is fairly fast.
221
222    AUTHOR: David Joyner (2005-11)
223    """
224    q = F.order()
225    G = gap(Gmat)
226    k = gap(F)
227    C = G.GeneratorMatCode(k)
228    n = int(C.WordLength())
229    z = 'Z(%s)*%s'%(q, [0]*n)     # GAP zero vector as a string
230    dist = gap.eval("w:=DistancesDistributionMatFFEVecFFE("+Gmat+", GF("+str(q)+"),"+z+")")
231    ## for some reason, this commented code doesn't work:
232    #dist0 = gap("DistancesDistributionMatFFEVecFFE("+Gmat+", GF("+str(q)+"),"+z+")")
233    #v0 = dist0._matrix_(F)
234    #print dist0,v0
235    #d = G.DistancesDistributionMatFFEVecFFE(k, z)
236    v = [eval(gap.eval("w["+str(i)+"]")) for i in range(1,n+2)] ## because GAP returns vectors in compressed form
237    return v
238
239def min_wt_vec(Gmat,F):
240    r"""
241    Uses C programs written by Steve Linton in the kernel of GAP, so
242    is fairly fast.
243
244    INPUT:
245        Same as wtdist.
246
247    OUTPUT:
248        Returns a minimum weight vector v of the code generated by Gmat
249         ## , the"message" vector m such that m*G = v, and the (minimum) distance, as a triple.
250
251    EXAMPLES:
252        sage: Gstr = "Z(2)*[[1,1,1,0,0,0,0], [1,0,0,1,1,0,0], [0,1,0,1,0,1,0], [1,1,0,1,0,0,1]]"
253        sage: sage.coding.linear_code.min_wt_vec(Gstr,GF(2))
254        (0, 0, 1, 0, 1, 1, 0)
255
256    Here Gstr is a generator matrix of the Hamming [7,4,3] binary code.
257
258    AUTHOR: David Joyner (11-2005)
259    """
260    from sage.rings.finite_field import gap_to_sage
261    gap.eval("G:="+Gmat)
262    k = int(gap.eval("Length(G)"))
263    q = F.order()
264    qstr = str(q)
265    C = gap(Gmat).GeneratorMatCode(F)
266    n = int(C.WordLength())
267    cg = C.MinimumDistanceCodeword()
268    c = [gap_to_sage(cg[j],F) for j in range(1,n+1)]
269    V = VectorSpace(F,n)
270    return V(c)
271    ## this older code returns more info but may be slower:
272    #zerovec = [0 for i in range(n)]
273    #zerovecstr = "Z("+qstr+")*"+str(zerovec)
274    #all = []
275    #for i in range(1,k+1):
276    #    P = gap.eval("P:=AClosestVectorCombinationsMatFFEVecFFECoords("+Gmat+", GF("+qstr+"),"+zerovecstr+","+str(i)+","+str(0)+"); d:=WeightVecFFE(P[1])")
277    #    v = gap("[List(P[1], i->i)]")
278    #    m = gap("[List(P[2], i->i)]")
279    #    dist = gap.eval("d")
280    #    #print v,m,dist
281    #    #print [gap.eval("v["+str(i+1)+"]") for i in range(n)]
282    #    all.append([v._matrix_(F),m._matrix_(F),int(dist)])
283    #ans = all[0]
284    #for x in all:
285    #    if x[2]<ans[2] and x[2]>0:
286    #        ans = x
287    #return ans
288
289## def minimum_distance_lower_bound(n,k,F):
290##     r"""
291##     Connects to \verb+http://www.win.tue.nl/~aeb/voorlincod.html+
292##     Tables of A. E. Brouwer,   Techn. Univ. Eindhoven,
293##     via Steven Sivek's linear_code_bound.
294
295##     EXAMPLES:
296##     sage: sage.coding.linear_code.minimum_distance_upper_bound(7,4,GF(2))     # optional (net connection)
297##         3
298
299##     Obviously requires an internet connection.
300##     """
301##     q = F.order()
302##     bounds = linear_code_bound(q,n,k)
303##     return bounds[0]
304
305## def minimum_distance_upper_bound(n,k,F):
306##     r"""
307##     Connects to \verb+http://www.win.tue.nl/~aeb/voorlincod.html+
308##     Tables of A. E. Brouwer,   Techn. Univ. Eindhoven
309##     via Steven Sivek's linear_code_bound.
310
311##     EXAMPLES:
312##         sage: sage.coding.linear_code.minimum_distance_upper_bound(7,4,GF(2))  # optional (net connection)
313##         3
314
315##     Obviously requires an internet connection.
316##     """
317##     q = F.order()
318##     bounds = linear_code_bound(q,n,k)
319##     return bounds[1]
320
321## def minimum_distance_why(n,k,F):
322##     r"""
323##     Connects to http://www.win.tue.nl/~aeb/voorlincod.html
324##     Tables of A. E. Brouwer,   Techn. Univ. Eindhoven
325##     via Steven Sivek's linear_code_bound.
326
327##     EXAMPLES:
328##         sage: sage.coding.linear_code.minimum_distance_why(7,4,GF(2))  # optional (net connection)
329##         Lb(7,4) = 3 is found by truncation of:
330##         Lb(8,4) = 4 is found by the (u|u+v) construction
331##         applied to [4,3,2] and [4,1,4]-codes
332##         Ub(7,4) = 3 follows by the Griesmer bound.
333
334##     Obviously requires an internet connection.
335##     """
336##     q = F.order()
337##     bounds = linear_code_bound(q,n,k)
338##     print bounds[2]
339
340def best_known_linear_code(n,k,F):
341    r"""
342    best_known_linear_code returns the best known (as of 11 May 2006)
343    linear code of length n, dimension k over field F. The function
344    uses the tables described in bounds_minimum_distance to construct
345    this code.
346
347    This does not require an internet connection.
348
349    EXAMPLES:
350        sage: best_known_linear_code(10,5,GF(2))    # long time
351        Linear code of length 10, dimension 5 over Finite Field of size 2
352        sage: gap.eval("C:=BestKnownLinearCode(10,5,GF(2))")     # long time
353        'a linear [10,5,4]2..4 shortened code'
354
355    This means that best possible binary linear code of length 10 and dimension 5
356    is a code with minimum distance 4 and covering radius somewhere between 2 and 4.
357    Use "minimum_distance_why(10,5,GF(2))" or "print bounds_minimum_distance(10,5,GF(2))"
358    for further details.
359    """
360    q = F.order()
361    C = gap("BestKnownLinearCode(%s,%s,GF(%s))"%(n,k,q))
362    G = C.GeneratorMat()
363    k = G.Length()
364    n = G[1].Length()
365    Gs = G._matrix_(F)
366    MS = MatrixSpace(F,k,n)
367    return LinearCode(MS(Gs))
368    #return gap.eval("BestKnownLinearCode(%s,%s,GF(%s))"%(n,k,q))
369
370def bounds_minimum_distance(n,k,F):
371    r"""
372    The function bounds_minimum_distance calculates a lower and upper
373    bound for the minimum distance of an optimal linear code with word
374    length n, dimension k over field F. The function returns a record
375    with the two bounds and an explanation for each bound. The
376    function Display can be used to show the explanations.
377
378    The values for the lower and upper bound are obtained from a table
379    constructed by Cen Tjhai for GUAVA, derived from the table of
380    Brouwer. (See http://www.win.tue.nl/~aeb/voorlincod.html or use
381    the SAGE function minimum_distance_why for the most recent data.)
382    These tables contain lower and upper bounds for q=2 (n <= 257), 3
383    (n <= 243), 4 (n <= 256). (Current as of 11 May 2006.)  For codes
384    over other fields and for larger word lengths, trivial bounds are
385    used.
386
387    This does not require an internet connection. The format of the
388    output is a little non-intuitive. Try print
389    bounds_minimum_distance(10,5,GF(2)) for example.
390    """
391    q = F.order()
392    gap.eval("data := BoundsMinimumDistance(%s,%s,GF(%s))"%(n,k,q))
393    Ldata = gap.eval("Display(data)")
394    return Ldata
395
396def LinearCode_from_vectorspace(self):
397    r"""
398    Converts a VectorSpace over GF(q) into a LinearCode
399
400    EXAMPLES:
401        sage: V = VectorSpace(GF(2),7)
402        sage: W = V.subspace([[1,0,0,1,1,1,1], [1,1,0,0,0,1,1]])
403        sage: C = LinearCode_from_vectorspace(W)
404        sage: C
405        Linear code of length 7, dimension 2 over Finite Field of size 2
406        sage: C.gen_mat()
407        [1 0 0 1 1 1 1]
408        [0 1 0 1 1 0 0]
409    """
410    F = self.base_ring()
411    B = self.basis()
412    n = len(B[0].list())
413    k = len(B)
414    MS = MatrixSpace(F,k,n)
415    G = MS([B[i].list() for i in range(k)])
416    return LinearCode(G)
417
418
419########################### linear codes python class #######################
420
421class LinearCode(module.Module):
422    r"""
423    A class for linear codes over a finite field or finite ring.  Each
424    instance is a linear code determined by a generator matrix $G$
425    (i.e., a k x n matrix of (full) rank $k$, $k\leq n$ over a finite
426    field $F$.
427
428    INPUT:
429        G -- a generator matrix over $F$. (G can be defined over a
430             finite ring but the matrices over that ring must have
431             certain attributes, such as"rank".)
432
433    OUTPUT:
434        The linear code of length $n$ over $F$ having $G$ as a
435        generator matrix.
436
437    EXAMPLES:
438        sage: MS = MatrixSpace(GF(2),4,7)
439        sage: G  = MS([[1,1,1,0,0,0,0], [1,0,0,1,1,0,0], [0,1,0,1,0,1,0], [1,1,0,1,0,0,1]])
440        sage: C  = LinearCode(G)
441        sage: C
442        Linear code of length 7, dimension 4 over Finite Field of size 2
443        sage: C.base_ring()
444        Finite Field of size 2
445        sage: C.dimension()
446        4
447        sage: C.length()
448        7
449        sage: C.minimum_distance()
450        3
451        sage: C.spectrum()
452        [1, 0, 0, 7, 7, 0, 0, 1]
453        sage: C.weight_distribution()
454        [1, 0, 0, 7, 7, 0, 0, 1]
455        sage: MS = MatrixSpace(GF(5),4,7)
456        sage: G  = MS([[1,1,1,0,0,0,0], [1,0,0,1,1,0,0], [0,1,0,1,0,1,0], [1,1,0,1,0,0,1]])
457        sage: C  = LinearCode(G)
458        sage: C
459        Linear code of length 7, dimension 4 over Finite Field of size 5
460
461    AUTHOR: David Joyner (11-2005)
462    """
463    #    sage: C.minimum_distance_upper_bound()   # optional (net connection)
464    #    3
465    #    sage: C.minimum_distance_why()     # optional (net connection)
466    #    Ub(7,4) = 3 follows by the Griesmer bound.
467    def __init__(self, gen_mat):
468        base_ring = gen_mat[0,0].parent()
469        ParentWithGens.__init__(self, base_ring)
470        self.__gens = gen_mat.rows()
471        self.__gen_mat = gen_mat
472        self.__length = len(gen_mat.row(0))
473        self.__dim = gen_mat.rank()
474
475    def length(self):
476        return self.__length
477
478    def dimension(self):
479        return self.__dim
480
481    def _repr_(self):
482        return "Linear code of length %s, dimension %s over %s"%(self.length(), self.dimension(), self.base_ring())
483
484    def random(self):
485        """
486        EXAMPLES:
487            sage: C = HammingCode(3,GF(4,'a'))
488            sage: Cc = C.galois_closure(GF(2))
489            sage: c = C.random()
490            sage: V = VectorSpace(GF(4,'a'),21)
491            sage: c2 = V([x^2 for x in c.list()])
492            sage: c2 in C
493            False
494
495        """
496        from sage.rings.finite_field import gap_to_sage
497        F = self.base_ring()
498        q = F.order()
499        G = self.gen_mat()
500        n = len(G.columns())
501        Cg = gap(G).GeneratorMatCode(F)
502        c = Cg.Random()
503        ans = [gap_to_sage(c[i],F) for i in range(1,n+1)]
504        V = VectorSpace(F,n)
505        return V(ans)
506
507    def gen_mat(self):
508        return self.__gen_mat
509
510    def gens(self):
511        return self.__gens
512
513    def basis(self):
514        return self.__gens
515
516    def list(self):
517        r"""
518        Return list of all elements of this linear code.
519
520        EXAMPLES:
521            sage: C = HammingCode(3,GF(2))
522            sage: Clist = C.list()
523            sage: Clist[5]; Clist[5] in C
524            (1, 0, 1, 0, 1, 0, 1)
525            True
526        """
527        n = self.length()
528        k = self.dimension()
529        F = self.base_ring()
530        Cs,p = self.standard_form()
531        Gs = Cs.gen_mat()
532        V = VectorSpace(F,k)
533        MS = MatrixSpace(F,n,n)
534        ans = []
535        perm_mat = MS(p.matrix().rows())**(-1)
536        for v in V:
537            ans.append((v*Gs)*perm_mat)
538        return ans
539
540    def __iter__(self):
541        """
542        Return an iterator over the elements of this linear code.
543
544        EXAMPLES:
545            sage: C = HammingCode(3,GF(2))
546            sage: [list(c) for c in C if hamming_weight(c) < 4]
547            [[0, 0, 0, 0, 0, 0, 0],
548             [1, 0, 0, 0, 0, 1, 1],
549             [0, 1, 0, 0, 1, 0, 1],
550             [0, 0, 1, 0, 1, 1, 0],
551             [1, 1, 1, 0, 0, 0, 0],
552             [1, 0, 0, 1, 1, 0, 0],
553             [0, 1, 0, 1, 0, 1, 0],
554             [0, 0, 1, 1, 0, 0, 1]]
555        """
556        n = self.length()
557        k = self.dimension()
558        F = self.base_ring()
559        Cs,p = self.standard_form()
560        Gs = Cs.gen_mat()
561        V = VectorSpace(F,k)
562        MS = MatrixSpace(F,n,n)
563        perm_mat = MS(p.matrix().rows())**(-1)
564        for v in V:
565            yield (v*Gs)*perm_mat
566
567    def ambient_space(self):
568        return VectorSpace(self.base_ring(),self.__length)
569
570    def __contains__(self,v):
571        A = self.ambient_space()
572        C = A.subspace(self.gens())
573        return C.__contains__(v)
574
575    def characteristic(self):
576        return (self.base_ring()).characteristic()
577
578##     def minimum_distance_lower_bound(self):
579##         r"""
580##         Connects to \verb+http://www.win.tue.nl/~aeb/voorlincod.html+
581##         Tables of A. E. Brouwer,   Techn. Univ. Eindhoven
582
583##         Obviously requires an internet connection
584##         """
585##         q = (self.base_ring()).order()
586##         n = self.length()
587##         k = self.dimension()
588##         bounds = linear_code_bound(q,n,k)
589##         return bounds[0]
590
591##     def minimum_distance_upper_bound(self):
592##         r"""
593##         Connects to http://www.win.tue.nl/~aeb/voorlincod.html
594##         Tables of A. E. Brouwer,   Techn. Univ. Eindhoven
595
596##         Obviously requires an internet connection
597##         """
598##         q = (self.base_ring()).order()
599##         n = self.length()
600##         k = self.dimension()
601##         bounds = linear_code_bound(q,n,k)
602##         return bounds[1]
603
604##     def minimum_distance_why(self):
605##         r"""
606##         Connects to http://www.win.tue.nl/~aeb/voorlincod.html
607##         Tables of A. E. Brouwer,   Techn. Univ. Eindhoven
608
609##         Obviously requires an internet connection.
610##         """
611##         q = (self.base_ring()).order()
612##         n = self.length()
613##         k = self.dimension()
614##         bounds = linear_code_bound(q,n,k)
615##         lines = bounds[2].split("\n")
616##         for line in lines:
617##             if len(line)>0:
618##                 if line[0] == "U":
619##                     print line
620
621    def minimum_distance(self):
622        r"""
623        Uses a GAP kernel function (in C) written by Steve Linton.
624
625        EXAMPLES:
626            sage: MS = MatrixSpace(GF(3),4,7)
627            sage: G = MS([[1,1,1,0,0,0,0], [1,0,0,1,1,0,0], [0,1,0,1,0,1,0], [1,1,0,1,0,0,1]])
628            sage: C = LinearCode(G)
629            sage: C.minimum_distance()
630            3
631        """
632        #sage: C.minimum_distance_upper_bound()  # optional (net connection)
633        #5
634        #    sage: C.minimum_distance_why()          # optional (net connection)
635        #    Ub(10,5) = 5 follows by the Griesmer bound.
636
637        F = self.base_ring()
638        q = F.order()
639        G = self.gen_mat()
640        gapG = gap(G)
641        Gstr = "%s*Z(%s)^0"%(gapG, q)
642        return hamming_weight(min_wt_vec(Gstr,F))
643
644    def genus(self):
645        r"""
646        Returns the "Duursma genus" of the code, gamma_C = n+1-k-d.
647        """
648        d = self.minimum_distance()
649        n = self.length()
650        k = self.dimension()
651        gammaC = n+1-k-d
652        return gammaC
653
654    def spectrum(self):
655        r"""
656        Uses a GAP kernel function (in C) written by Steve Linton.
657
658        EXAMPLES:
659            sage: MS = MatrixSpace(GF(2),4,7)
660            sage: G = MS([[1,1,1,0,0,0,0], [ 1, 0, 0, 1, 1, 0, 0], [ 0, 1, 0,1, 0, 1, 0], [1, 1, 0, 1, 0, 0, 1]])
661            sage: C = LinearCode(G)
662            sage: C.spectrum()
663            [1, 0, 0, 7, 7, 0, 0, 1]
664
665        """
666        F = self.base_ring()
667        q = F.order()
668        G = self.gen_mat()
669        Glist = [list(x) for x in G]
670        Gstr = "Z("+str(q)+")*"+str(Glist)
671        spec = wtdist(Gstr,F)
672        return spec
673
674    def weight_distribution(self):
675        #same as spectrum
676        return self.spectrum()
677
678    def __cmp__(self, right):
679        if not isinstance(right, LinearCode):
680            return cmp(type(self), type(right))
681        return cmp(self.__gen_mat, right.__gen_mat)
682
683    def decode(self, right):
684        r"""
685        Wraps GUAVA's Decodeword. Hamming codes have a special
686        decoding algorithm. Otherwise, syndrome decoding is used.
687
688        INPUT:
689            right must be a vector of length = length(self)
690
691        OUTPUT:
692            The codeword c in C closest to r.
693
694        EXAMPLES:
695            sage: C = HammingCode(3,GF(2))
696            sage: MS = MatrixSpace(GF(2),1,7)
697            sage: F = GF(2); a = F.gen()
698            sage: v = [a,a,F(0),a,a,F(0),a]
699            sage: C.decode(v)
700            (1, 1, 0, 1, 0, 0, 1)
701
702        Does not work for very long codes since the syndrome table grows too large.
703        """
704        from sage.rings.finite_field import gap_to_sage
705        F = self.base_ring()
706        q = F.order()
707        G = self.gen_mat()
708        n = len(G.columns())
709        k = len(G.rows())
710        Gstr = str(gap(G))
711        vstr = str(gap(right))
712        if vstr[:3] == '[ [':
713            vstr = vstr[1:-1]     # added by William Stein so const.tex works 2006-10-01
714        gap.eval("C:=GeneratorMatCode("+Gstr+",GF("+str(q)+"))")
715        gap.eval("c:=VectorCodeword(Decodeword( C, Codeword( "+vstr+" )))") # v->vstr, 8-27-2006
716        ans = [gap_to_sage(gap.eval("c["+str(i)+"]"),F) for i in range(1,n+1)]
717        V = VectorSpace(F,n)
718        return V(ans)
719
720    def dual_code(self):
721        r"""
722        Wraps GUAVA's DualCode.
723
724        OUTPUT:
725            The dual code.
726
727        EXAMPLES:
728            sage: C = HammingCode(3,GF(2))
729            sage: C.dual_code()
730            Linear code of length 7, dimension 3 over Finite Field of size 2
731            sage: C = HammingCode(3,GF(4,'a'))
732            sage: C.dual_code()
733            Linear code of length 21, dimension 3 over Finite Field in a of size 2^2
734        """
735        F = self.base_ring()
736        q = F.order()
737        G = self.gen_mat()
738        n = len(G.columns())
739        k = len(G.rows())
740        if n==k:
741            return TrivialCode(F,n)
742        Gstr = str(gap(G))
743        #H = C.CheckMat()
744        #A = H._matrix_(GF(q))
745        #return LinearCode(A)       ## This does not work when k = n-1 for a mysterious reason.
746        ##  less pythonic way :
747        C = gap("DualCode(GeneratorMatCode("+Gstr+",GF("+str(q)+")))")
748        G = C.GeneratorMat()
749        Gs = G._matrix_(F)
750        MS = MatrixSpace(F,n-k,n)
751        return LinearCode(MS(Gs))
752
753    def extended_code(self):
754        r"""
755        If self is a linear code of length n defined over F then this
756        returns the code of length n+1 where the last digit $c_n$
757        satisfies the check condition $c_0+...+c_n=0$. If self is an
758        $[n,k,d]$ binary code then the extended code $C^{\vee}$ is an
759        $[n+1,k,d^{\vee}]$ code, where $d^=d$ (if d is even) and
760        $d^{\vee}=d+1$ (if $d$ is odd).
761
762        EXAMPLES:
763            sage: C = HammingCode(3,GF(4,'a'))
764            sage: C
765            Linear code of length 21, dimension 18 over Finite Field in a of size 2^2
766            sage: Cx = C.extended_code()
767            sage: Cx
768            Linear code of length 22, dimension 18 over Finite Field in a of size 2^2
769        """
770        G = self.gen_mat()
771        F = self.base_ring()
772        q = F.order()
773        n = len(G.columns())
774        k = len(G.rows())
775        g = gap(G)
776        gap.eval( "G:="+g.name())
777        C = gap("GeneratorMatCode(G,GF("+str(q)+"))")
778        Cx = C.ExtendedCode()
779        Gx = Cx.GeneratorMat()
780        Gxs = Gx._matrix_(F)           # this is the killer
781        MS = MatrixSpace(F,k,n+1)
782        return LinearCode(MS(Gxs))
783
784    def check_mat(self):
785        r"""
786        Returns the check matrix of self.
787
788        EXAMPLES:
789            sage: C = HammingCode(3,GF(2))
790            sage: Cperp = C.dual_code()
791            sage: C; Cperp
792            Linear code of length 7, dimension 4 over Finite Field of size 2
793            Linear code of length 7, dimension 3 over Finite Field of size 2
794            sage: C.gen_mat()
795            [1 1 1 0 0 0 0]
796            [1 0 0 1 1 0 0]
797            [0 1 0 1 0 1 0]
798            [1 1 0 1 0 0 1]
799            sage: C.check_mat()
800            [0 1 1 1 1 0 0]
801            [1 0 1 1 0 1 0]
802            [1 1 0 1 0 0 1]
803            sage: Cperp.check_mat()
804            [1 1 1 0 0 0 0]
805            [1 0 0 1 1 0 0]
806            [0 1 0 1 0 1 0]
807            [1 1 0 1 0 0 1]
808            sage: Cperp.gen_mat()
809            [0 1 1 1 1 0 0]
810            [1 0 1 1 0 1 0]
811            [1 1 0 1 0 0 1]
812        """
813        Cperp = self.dual_code()
814        return Cperp.gen_mat()
815
816    def __eq__(self, right):
817        """
818        Checks if self == right.
819
820        EXAMPLES:
821        """
822        slength = self.length()
823        rlength = right.length()
824        sdim = self.dimension()
825        rdim = right.dimension()
826        sF = self.base_ring()
827        rF = right.base_ring()
828        if slength != rlength:
829            return False
830        if sdim != rdim:
831            return False
832        if sF != rF:
833            return False
834        V = VectorSpace(sF,sdim)
835        sbasis = self.gens()
836        rbasis = right.gens()
837        scheck = self.check_mat()
838        rcheck = right.check_mat()
839        for c in sbasis:
840            if rcheck*c:
841                return False
842        for c in rbasis:
843            if scheck*c:
844                return False
845        return True
846
847    def is_permutation_automorphism(self,g):
848        r"""
849        Returns $1$ if $g$ is an element of $S_n$ ($n$ = length of
850        self) and if $g$ is an automorphism of self.
851
852        EXAMPLES:
853            sage: C = HammingCode(3,GF(3))
854            sage: g = SymmetricGroup(13).random_element()
855            sage: C.is_permutation_automorphism(g)
856            0
857            sage: MS = MatrixSpace(GF(2),4,8)
858            sage: G  = MS([[1,0,0,0,1,1,1,0],[0,1,1,1,0,0,0,0],[0,0,0,0,0,0,0,1],[0,0,0,0,0,1,0,0]])
859            sage: C  = LinearCode(G)
860            sage: S8 = SymmetricGroup(8)
861            sage: g = S8("(2,3)")
862            sage: C.is_permutation_automorphism(g)
863            1
864            sage: g = S8("(1,2,3,4)")
865            sage: C.is_permutation_automorphism(g)
866            0
867
868        """
869        basis = self.gen_mat().rows()
870        H = self.check_mat()
871        V = H.column_space()
872        HGm = H*g.matrix()
873        # raise TypeError, (type(H), type(V), type(basis[0]), type(Gmc))
874        for c in basis:
875            if HGm*c != V(0):
876                return 0
877        return 1
878
879    def permutation_automorphism_group(self,mode=None):
880        r"""
881        If $C$ is an $[n,k,d]$ code over $F$, this function computes
882        the subgroup $Aut(C) \subset S_n$ of all permutation
883        automorphisms of $C$.  If mode="verbose" then
884        code-theoretic data is printed out at several stages
885        of the computation.
886
887        Combines an idea of mine with an improvement suggested by Cary
888        Huffman.
889
890        EXAMPLES:
891            sage: MS = MatrixSpace(GF(2),4,8)
892            sage: G  = MS([[1,0,0,0,1,1,1,0],[0,1,1,1,0,0,0,0],[0,0,0,0,0,0,0,1],[0,0,0,0,0,1,0,0]])
893            sage: C  = LinearCode(G)
894            sage: C
895            Linear code of length 8, dimension 4 over Finite Field of size 2
896            sage: G = C.permutation_automorphism_group()   # long time
897            sage: G.order()                                # long time
898            144
899
900        A less easy example involves showing that the permutation automorphism
901        group of the extended ternary Golay code is the Mathieu group $M_{11}$.
902
903            sage: C = ExtendedTernaryGolayCode()
904            sage: M11 = MathieuGroup(11)
905            sage: G = C.permutation_automorphism_group()  ## this should take < 15 seconds  # long time
906            sage: G.is_isomorphic(M11)        # long time
907            True
908
909        WARNING: - *Ugly code, which should be replaced by a call to Robert Miller's nice program.*
910                 - Known to mysteriously crash in one example.
911        """
912        F = self.base_ring()
913        q = F.order()
914        G = self.gen_mat()
915        n = len(G.columns())
916        k = len(G.rows())                                 ## G is always full rank
917        gap.eval("Gp:=SymmetricGroup(%s)"%n)               ## initializing G in gap
918        Sn = SymmetricGroup(n)
919        wts = self.spectrum()                                            ## bottleneck 1
920        Gstr = str(gap(G))
921        gap.eval("C:=GeneratorMatCode("+Gstr+",GF("+str(q)+"))")
922        gap.eval("eltsC:=Elements(C)")
923        nonzerowts = [i for i in range(len(wts)) if wts[i]!=0]
924        if mode=="verbose":
925            print "\n Minimum distance: %s \n Weight distribution: \n %s"%(nonzerowts[1],wts)
926        stop = 0                                          ## only stop if all gens are autos
927        for i in range(1,len(nonzerowts)):
928          if stop == 1:
929              break
930          wt = nonzerowts[i]
931          if mode=="verbose":
932              size = eval(gap.eval("Size(Gp)"))
933              print "\n Using the %s codewords of weight %s \n Supergroup size: \n %s\n "%(wts[wt],wt,size)
934          gap.eval("Cwt:=Filtered(eltsC,c->WeightCodeword(c)=%s)"%wt)   ## bottleneck 2 (repeated
935          gap.eval("matCwt:=List(Cwt,c->VectorCodeword(c))")            ##        for each i until stop = 1)
936          gap.eval("A:=MatrixAutomorphisms(matCwt); GG:=Intersection(Gp,A)")    ## bottleneck 3
937          #print i," = i \n",gap.eval("matCwt")," = matCwt\n"
938          #print gap.eval("A")," = A \n",gap.eval("GG")," = GG\n\n"
939          if eval(gap.eval("Size(GG)"))==0:
940              #print gap.eval("GG; Size(GG)")
941              #print "GG=0 ", gap.eval("A; Gp")
942              return PermutationGroup([()])
943          gap.eval("autgp_gens:=GeneratorsOfGroup(GG); Gp:=GG")
944          #gap.eval("autgp_gens:=GeneratorsOfGroup(G)")
945          N = eval(gap.eval("Length(autgp_gens)"))
946          gens = [Sn(gap.eval("autgp_gens[%s]"%i).replace("\n","")) for i in range(1,N+1)]
947          stop = 1                                        ## get ready to stop
948          for x in gens:                                  ## if one of these gens is not an auto then don't stop
949              if not(self.is_permutation_automorphism(x)):
950                  stop = 0
951                  break
952        G = PermutationGroup(gens)
953        return G
954
955    def permuted_code(self,p):
956        r"""
957        Returns the permuted code -- the code $C$ which is equivalent
958        to self via the column permutation $p$.
959
960        EXAMPLES:
962            sage: G = C.permutation_automorphism_group()     # long time
963            sage: p = G("(1,6,3,5)(2,7,4,8)")                # long time
964            sage: Cp = C.permuted_code(p)                    # long time
965            sage: C.gen_mat()
966            [1 1 0 1 0 0 0 1]
967            [0 1 1 0 1 0 0 1]
968            [0 0 1 1 0 1 0 1]
969            [0 0 0 1 1 0 1 1]
970            sage: Cp.gen_mat()                               # long time
971            [0 1 0 0 0 1 1 1]
972            [1 1 0 0 1 0 1 0]
973            [0 1 1 0 1 0 0 1]
974            [1 1 0 1 0 0 0 1]
975            sage: Cs1,p1 = C.standard_form(mode="verbose"); p1
976            1 . . . 1 1 . 1
977            . 1 . . . 1 1 1
978             . . 1 . 1 1 1 .
979             . . . 1 1 . 1 1
980            ()
981            sage: Cs2,p2 = Cp.standard_form(mode="verbose"); p2   # long time
982             1 . . . 1 1 . 1
983             . 1 . . . 1 1 1
984             . . 1 . 1 1 1 .
985             . . . 1 1 . 1 1
986            ()
987
988        Therefore you can see that Cs1 and Cs2 are the same, so C and Cp are
989        equivalent.
990
991        """
992        F = self.base_ring()
993        G = self.gen_mat()
994        n = len(G.columns())
995        MS = MatrixSpace(F,n,n)
996        Gp = G*MS(p.matrix().rows())
997        return LinearCode(Gp)
998
999    def standard_form(self,mode=None):
1000        r"""
1001        An $[n,k]$ linear code with generator matrix $G$ is in
1002        standard form is the row-reduced echelon form of $G$ is
1003        $(I,A)$, where $I$ denotes the $k \times k$ identity matrix
1004        and $A$ is a $k \times (n-k)$ block.  This method returns a
1005        pair $(C,p)$ where $C$ is a code permutation equivalent to
1006        self and $p$ in $S_n$ ($n$ = length of $C$) is the permutation
1007        sending self to $C$.  When mode ="verbose" the new generator
1008        matrix in the standard form $(I,A)$ is "Display"'d.
1009
1010        EXAMPLES:
1012            sage: C.gen_mat()
1013            [1 1 0 1 0 0 0 1]
1014            [0 1 1 0 1 0 0 1]
1015            [0 0 1 1 0 1 0 1]
1016            [0 0 0 1 1 0 1 1]
1017            sage: Cs,p = C.standard_form()
1018            sage: Cs.gen_mat()
1019            [1 0 0 0 1 1 0 1]
1020            [0 1 0 0 0 1 1 1]
1021            [0 0 1 0 1 1 1 0]
1022            [0 0 0 1 1 0 1 1]
1023            sage: p
1024            ()
1025            sage: C.standard_form(mode="verbose")
1026            <BLANKLINE>
1027            1 . . . 1 1 . 1
1028            . 1 . . . 1 1 1
1029            . . 1 . 1 1 1 .
1030            . . . 1 1 . 1 1
1031            <BLANKLINE>
1032            (Linear code of length 8, dimension 4 over Finite Field of size 2, ())
1033
1034        """
1035        F = self.base_ring()
1036        q = F.order()
1037        G = self.gen_mat()
1038        n = len(G.columns())
1039        k = len(G.rows())                                 ## G is always full rank
1040        Sn = SymmetricGroup(n)
1041        Gstr = str(gap(G))
1042        gap.eval( "G:="+Gstr )
1043        p = Sn(gap.eval("p:=PutStandardForm(G)"))
1044        if mode=="verbose":
1045            print "\n",gap.eval("Display(G)"),"\n"
1046        C = gap("GeneratorMatCode(G,GF("+str(q)+"))")
1047        Gp = C.GeneratorMat()
1048        Gs = Gp._matrix_(F)
1049        MS = MatrixSpace(F,k,n)
1050        return LinearCode(MS(Gs)),p
1051
1052    def module_composition_factors(self,gp):
1053        r"""
1054        Prints the GAP record of the Meataxe composition factors module in
1055        Meataxe notation.
1056
1057        EXAMPLES:
1058            sage: MS = MatrixSpace(GF(2),4,8)
1059            sage: G  = MS([[1,0,0,0,1,1,1,0],[0,1,1,1,0,0,0,0],[0,0,0,0,0,0,0,1],[0,0,0,0,0,1,0,0]])
1060            sage: C  = LinearCode(G)
1061            sage: gp = C.permutation_automorphism_group()    # long time
1062
1063        Now type "C.module_composition_factors(gp)" to get the record printed.
1064
1065        """
1066        F = self.base_ring()
1067        q = F.order()
1068        gens = gp.gens()
1069        G = self.gen_mat()
1070        n = len(G.columns())
1071        k = len(G.rows())
1072        MS = MatrixSpace(F,n,n)
1073        mats = [] # initializing list of mats by which the gens act on self
1074        Fn = VectorSpace(F, n)
1075        W = Fn.subspace_with_basis(G.rows()) # this is self
1076        for g in gens:
1077            p = MS(g.matrix())
1078            m = [W.coordinate_vector(r*p) for r in G.rows()]
1079            mats.append(m)
1080        mats_str = str(gap([[list(r) for r in m] for m in mats]))
1081        gap.eval("M:=GModuleByMats("+mats_str+", GF("+str(q)+"))")
1082        print gap("MTX.CompositionFactors( M )")
1083
1084    def galois_closure(self, F0):
1085        r"""
1086        If self is a linear code defined over $F$ and $F_0$ is a subfield
1087        with Galois group $G = Gal(F/F_0)$ then this returns the $G$-module
1088        $C^-$ containing $C$.
1089
1090        EXAMPLES:
1091            sage: C = HammingCode(3,GF(4,'a'))
1092            sage: Cc = C.galois_closure(GF(2))
1093            sage: C; Cc
1094            Linear code of length 21, dimension 18 over Finite Field in a of size 2^2
1095            Linear code of length 21, dimension 20 over Finite Field in a of size 2^2
1096            sage: c = C.random()
1097            sage: V = VectorSpace(GF(4,'a'),21)
1098            sage: c2 = V([x^2 for x in c.list()])
1099            sage: c2 in C
1100            False
1101            sage: c2 in Cc
1102            True
1103
1104        """
1105        G = self.gen_mat()
1106        F = self.base_ring()
1107        q = F.order()
1108        q0 = F0.order()
1109        a = log(q,q0)  ## test if F/F0 is a field extension
1110        if not(a.integer_part() == a):
1111            raise ValueError,"Base field must be an extension of given field %s"%F0
1112        n = len(G.columns())
1113        k = len(G.rows())
1114        G0 = [[x**q0 for x in g.list()] for g in G.rows()]
1115        G1 = [[x for x in g.list()] for g in G.rows()]
1116        G2 = G0+G1
1117        MS = MatrixSpace(F,2*k,n)
1118        G3 = MS(G2)
1119        r = G3.rank()
1120        MS = MatrixSpace(F,r,n)
1121        Grref = G3.echelon_form()
1122        G = MS([Grref.row(i) for i in range(r)])
1123        return LinearCode(G)
1124
1125    def is_galois_closed(self):
1126        r"""
1127        Checks if $C$ is equal to its Galois closure.
1128        """
1129        F = self.base_ring()
1130        p = F.characteristic()
1131        return self == self.galois_closure(GF(p))
1132
1133    def assmus_mattson_designs(self,t,mode=None):
1134        r"""
1135        Assmus and Mattson Theorem (section 8.4, page 303 of [HP]):
1136        Let $A_0, A_1, ..., A_n$ be the weights of the codewords in a
1137        binary linear $[n , k, d]$ code $C$, and let $A_0^*, A_1^*, 1138 ..., A_n^*$ be the weights of the codewords in its dual $[n, 1139 n-k, d^*]$ code $C^*$.  Fix a $t$, $0<t<d$, and let
1140        $$1141 s = |{ i | A_i^* \not= 0, 0<i\leq n-t}|. 1142$$
1143        Assume $s\leq d-t$.
1144
1145        (1) If $Ai\not= 0$ and $d\leq i\leq n then Ci = { c in C | wt(c) = i}$
1146            holds a simple t-design.
1147
1148        (2) If $Ai*\not= 0$ and $d*\leq i\leq n-t then Ci* = { c in C* 1149 | wt(c) = i}$ holds a simple t-design.
1150
1151        A {\bf block design} is a pair $(X,B)$, where $X$ is a
1152        non-empty finite set of $v>0$ elements called {\bf points},
1153        and $B$ is a non-empty finite multiset of size b whose
1154        elements are called {\bf blocks}, such that each block is a
1155        non-empty finite multiset of $k$ points. $A$ design without
1156        repeated blocks is called a {\bf simple} block design. If
1157        every subset of points of size $t$ is contained in exactly
1158        lambda blocks the the block design is called a {\bf
1159        $t-(v,k,lambda)$ design} (or simply a $t$-design when the
1160        parameters are not specfied). When $\lambda=1$ then the block
1161        design is called a {\bf $S(t,k,v)$ Steiner system}.
1162
1163        In the Assmus and Mattson Theorem (1), $X$ is the set
1164        $\{1,2,...,n\}$ of coordinate locations and
1165        $B = \{supp(c) | c in C_i\}$ is the set of
1166        supports of the codewords of $C$ of weight $i$.
1167        Therefore, the parameters of the $t$-design for $C_i$ are
1168        \begin{verbatim}
1169             t =       given
1170             v =       n
1171             k =       i   (k not to be confused with dim(C))
1172             b =       Ai
1173             lambda = b*binomial(k,t)/binomial(v,t) (by Theorem 8.1.6,
1174                      p 294, in [HP])
1175        \end{verbatim}
1176
1177        Setting the mode="verbose" option prints out the values of the parameters.
1178
1179        The first example below means that the binary [24,12,8]-code C
1180        has the property that the (support of the) codewords of weight
1181        8 (resp, 12, 16) form a 5-design.  Similarly for its dual code
1182        C* (of course C=C* in this case, so this info is extraneous).
1183        The test fails to produce 6-designs (ie, the hypotheses of the
1184        theorem fail to hold, not that the 6-designs definitely don't
1185        exist). The command assmus_mattson_designs(C,5,mode="verbose")
1186        returns the same value but prints out more detailed information.
1187
1188        The second example below illustrates the blocks of the 5-(24,
1189        8, 1) design (ie, the S(5,8,24) Steiner system).
1190
1191        EXAMPLES:
1192            sage: C = ExtendedBinaryGolayCode()             #  example 1
1193            sage: C.assmus_mattson_designs(5)
1194            ['weights from C: ',
1195            [8, 12, 16, 24],
1196            'designs from C: ',
1197            [[5, (24, 8, 1)], [5, (24, 12, 48)], [5, (24, 16, 78)], [5, (24, 24, 1)]],
1198            'weights from C*: ',
1199            [8, 12, 16],
1200            'designs from C*: ',
1201            [[5, (24, 8, 1)], [5, (24, 12, 48)], [5, (24, 16, 78)]]]
1202            sage: C.assmus_mattson_designs(6)
1203            0
1204            sage: X = range(24)                           #  example 2
1205            sage: blocks = [c.support() for c in C if hamming_weight(c)==8]; len(blocks)  ## long time computation
1206            759
1207
1208        REFERENCE:
1209            [HP] W. C. Huffman and V. Pless, {\bf Fundamentals of ECC}, Cambridge Univ. Press, 2003.
1210        """
1211        from sage.rings.arith import binomial
1212        C = self
1213        ans = []
1214        F = C.base_ring()
1215        q = F.order()
1216        G = C.gen_mat()
1217        n = len(G.columns())
1218        Cp = C.dual_code()
1219        k = len(G.rows())  ## G is always full rank
1220        wts = C.spectrum()
1221        d = min([i for i in range(1,len(wts)) if wts[i]!=0])
1222        if t>=d:
1223            return 0
1224        nonzerowts = [i for i in range(len(wts)) if wts[i]!=0 and i<=n and i>=d]
1225        #print d,t,len(nonzerowts)
1226        if mode=="verbose":
1227            #print "\n"
1228            for w in nonzerowts:
1229                print "The weight ",w," codewords of C form a t-(v,k,lambda) design, where"
1230                print "      t = ",t," , v = ",n," , k = ",w," , lambda = ",wts[w]*binomial(w,t)/binomial(n,t)#,"\n"
1231                print "      There are ",wts[w]," blocks of this design."
1232        wtsp = Cp.spectrum()
1233        dp = min([i for i in range(1,len(wtsp)) if wtsp[i]!=0])
1234        nonzerowtsp = [i for i in range(len(wtsp)) if wtsp[i]!=0 and i<=n-t and i>=dp]
1235        s = len([i for i in range(1,n) if wtsp[i]!=0 and i<=n-t and i>0])
1236        if mode=="verbose":
1237            #print "\n"
1238            for w in nonzerowtsp:
1239                print "The weight ",w," codewords of C* form a t-(v,k,lambda) design, where"
1240                print "      t = ",t," , v = ",n," , k = ",w," , lambda = ",wtsp[w]*binomial(w,t)/binomial(n,t)#,"\n"
1241                print "      There are ",wts[w]," blocks of this design."
1242        if s<=d-t:
1243            des = [[t,(n,w,wts[w]*binomial(w,t)/binomial(n,t))] for w in nonzerowts]
1244            ans = ans + ["weights from C: ",nonzerowts,"designs from C: ",des]
1245            desp = [[t,(n,w,wtsp[w]*binomial(w,t)/binomial(n,t))] for w in nonzerowtsp]
1246            ans = ans + ["weights from C*: ",nonzerowtsp,"designs from C*: ",desp]
1247            return ans
1248        return 0
1249
1250    def weight_enumerator(self,names="xy"):
1251        """
1252        Returns the weight enumerator of the code.
1253
1254        EXAMPLES:
1255            sage: C = HammingCode(3,GF(2))
1256            sage: C.weight_enumerator()
1257            x^7 + 7*x^4*y^3 + 7*x^3*y^4 + y^7
1258            sage: C.weight_enumerator(names="st")
1259            s^7 + 7*s^4*t^3 + 7*s^3*t^4 + t^7
1260        """
1261        spec = self.spectrum()
1262        n = self.length()
1263        from sage.rings.polynomial.polynomial_ring import PolynomialRing
1264        R = PolynomialRing(QQ,2,names)
1265        x,y = R.gens()
1266        we = sum([spec[i]*x**(n-i)*y**i for i in range(n+1)])
1267        return we
1268
1269    def zeta_polynomial(self,name = "T"):
1270        r"""
1271        Returns the Duursma zeta polynomial of the code.
1272
1273        Assumes minimum_distance(C) > 1 and minimum_distance$(C^\perp) > 1$.
1274
1275        EXAMPLES:
1276            sage: C = HammingCode(3,GF(2))
1277            sage: C.zeta_polynomial()
1278            2/5*T^2 + 2/5*T + 1/5
1279            sage: C = best_known_linear_code(6,3,GF(2))  # long time
1280            sage: C.minimum_distance()                   # long time (because of above)
1281            3
1282            sage: C.zeta_polynomial()                    # long time (because of above)
1283            2/5*T^2 + 2/5*T + 1/5
1284            sage: C = HammingCode(4,GF(2))
1285            sage: C.zeta_polynomial()
1286            16/429*T^6 + 16/143*T^5 + 80/429*T^4 + 32/143*T^3 + 30/143*T^2 + 2/13*T + 1/13
1287
1288        REFERENCES:
1289
1290             I. Duursma, "From weight enumerators to zeta functions",
1291             in {\bf Discrete Applied Mathematics}, vol. 111, no. 1-2,
1292             pp. 55-73, 2001.
1293        """
1294        n = self.length()
1295        q = (self.base_ring()).characteristic()
1296        d = self.minimum_distance()
1297        dperp = (self.dual_code()).minimum_distance()
1298        if d == 1 or dperp == 1:
1299            print "\n WARNING: There is no guarantee this function works when the minimum distance"
1300            print "            of the code or of the dual code equals 1.\n"
1301        from sage.rings.polynomial.polynomial_ring import PolynomialRing
1302        from sage.rings.fraction_field import FractionField
1303        from sage.rings.power_series_ring import PowerSeriesRing
1304        RT = PolynomialRing(QQ,"%s"%name)
1305        R = PolynomialRing(QQ,3,"xy%s"%name)
1306        x,y,T = R.gens()
1307        we = self.weight_enumerator()
1308        A = R(we)
1309        B = A(x+y,y,T)-(x+y)**n
1310        Bs = B.coefficients()
1311        b = [Bs[i]/binomial(n,i+d) for i in range(len(Bs))]
1312        r = n-d-dperp+2
1313        #print B,Bs,b,r
1314        P_coeffs = []
1315        for i in range(len(b)):
1316           if i == 0:
1317              P_coeffs.append(b[0])
1318           if i == 1:
1319              P_coeffs.append(b[1] - (q+1)*b[0])
1320           if i>1:
1321              P_coeffs.append(b[i] - (q+1)*b[i-1] + q*b[i-2])
1322        #print P_coeffs
1323        P = sum([P_coeffs[i]*T**i for i in range(r+1)])
1324        return RT(P)
1325
1326    def chinen_polynomial(self):
1327        """
1328        Returns the Chinen zeta polynomial of the code.
1329
1330        EXAMPLES:
1331            sage: C = HammingCode(3,GF(2))
1332            sage: C.chinen_polynomial()       # long time
1333            (2*sqrt(2)*t^3/5 + 2*sqrt(2)*t^2/5 + 2*t^2/5 + sqrt(2)*t/5 + 2*t/5 + 1/5)/(sqrt(2) + 1)
1334            sage: C = TernaryGolayCode()
1335            sage: C.chinen_polynomial()       # long time
1336            (6*sqrt(3)*t^3/7 + 6*sqrt(3)*t^2/7 + 6*t^2/7 + 2*sqrt(3)*t/7 + 6*t/7 + 2/7)/(2*sqrt(3) + 2)
1337
1338        This last output agrees with the corresponding example given in Chinen's paper below.
1339
1340        REFERENCES:
1341            Chinen, K. "An abundance of invariant polynomials
1342            satisfying the Riemann hypothesis", April 2007 preprint.
1343
1344        """
1345        from sage.rings.polynomial.polynomial_ring import PolynomialRing, polygen
1346        #from sage.calculus.functional import expand
1347        from sage.calculus.calculus import sqrt, SymbolicExpressionRing, var
1348        C = self
1349        n = C.length()
1350        RT = PolynomialRing(QQ,2,"Ts")
1351        T,s = FractionField(RT).gens()
1352        t = PolynomialRing(QQ,"t").gen()
1353        Cd = C.dual_code()
1354        k = C.dimension()
1355        q = (C.base_ring()).characteristic()
1356        d = C.minimum_distance()
1357        dperp = Cd.minimum_distance()
1358        if dperp > d:
1359            P = RT(C.zeta_polynomial())
1360            ## SAGE does not find dealing with sqrt(int) *as an algebraic object*
1361            ## an easy thing to do. Some tricky gymnastics are used to
1362            ## make SAGE deal with objects over QQ(sqrt(q)) nicely.
1363            if is_even(n):
1364                Pd = q**(k-n/2)*RT(Cd.zeta_polynomial())*T**(dperp - d)
1365            if not(is_even(n)):
1366                Pd = s*q**(k-(n+1)/2)*RT(Cd.zeta_polynomial())*T**(dperp - d)
1367            CP = P+Pd
1368            f = CP/CP(1,s)
1369            return f(t,sqrt(q))
1370        if dperp < d:
1371            P = RT(C.zeta_polynomial())*T**(d - dperp)
1372            if is_even(n):
1373                Pd = q**(k-n/2)*RT(Cd.zeta_polynomial())
1374            if not(is_even(n)):
1375                Pd = s*q**(k-(n+1)/2)*RT(Cd.zeta_polynomial())
1376            CP = P+Pd
1377            f = CP/CP(1,s)
1378            return f(t,sqrt(q))
1379        if dperp == d:
1380            P = RT(C.zeta_polynomial())
1381            if is_even(n):
1382                Pd = q**(k-n/2)*RT(Cd.zeta_polynomial())
1383            if not(is_even(n)):
1384                Pd = s*q**(k-(n+1)/2)*RT(Cd.zeta_polynomial())
1385            CP = P+Pd
1386            f = CP/CP(1,s)
1387            return f(t,sqrt(q))
1388
1389    def zeta_function(self,name = "T"):
1390        r"""
1391        Returns the Duursma zeta function of the code.
1392
1393        EXAMPLES:
1394            sage: C = HammingCode(3,GF(2))
1395            sage: C.zeta_function()
1396            (2/5*T^2 + 2/5*T + 1/5)/(2*T^2 - 3*T + 1)
1397        """
1398        P =  self.zeta_polynomial()
1399        q = (self.base_ring()).characteristic()
1400        RT = PolynomialRing(QQ,"%s"%name)
1401        T = RT.gen()
1402        return P/((1-T)*(1-q*T))
1403
1404    def zeta_function2(self,mode=None):
1405        r"""
1406        Returns the Duursma zeta function of the code.
1407
1408        NOTE: This is somewhat experimental code. It sometimes
1409        returns "fail" for a reason I don't fully understand.
1410        However, when it does return a polynomial, the answer is
1411        (as far as I know) correct.  *Experimental code* included to
1412        study a particular implementation.
1413
1414        INPUT:
1415
1416           mode -- string;
1417                -- "dual" computes both the zeta function
1418                   of $C$ and that of $C*$,
1419                -- "normalized" computes the normalized zeta function
1420                of $C$, $zeta_C(T)*T(1-genus(C))$ (NOTE: if xi(T,C)
1421                denotes the normalized zeta function then $xi(T,C*) = 1422 xi(1/(qT),C)$ is equivalent to $zeta_{C*}(T) = 1423 zeta_C(1/(qT))*q^(gamma-1)T^(gamma+gamma*-2)$, where
1424                $gamma = gamma(C)$ and $gamma* = gamma(C*)$.)
1425
1426        EXAMPLES:
1427            sage: C = HammingCode(3,GF(2))
1428            sage: C.zeta_function2()
1429            (2/5*T^2 + 2/5*T + 1/5)/(2*T^2 - 3*T + 1)
1430            sage: C = ExtendedTernaryGolayCode()
1431            sage: C.zeta_function2()
1432            (3/7*T^2 + 3/7*T + 1/7)/(3*T^2 - 4*T + 1)
1433
1434        Both these examples occur in Duursma's paper below.
1435
1436        REFERENCES:
1437
1438             I. Duursma, "Weight distributions of geometric Goppa codes,"
1439             Trans. A.M.S., 351 (1999)3609-3639.
1440
1441        """
1442        from sage.rings.polynomial.polynomial_ring import PolynomialRing
1443        from sage.rings.fraction_field import FractionField
1444        from sage.rings.power_series_ring import PowerSeriesRing
1445        R = PolynomialRing(QQ,3,"xyT")
1446        #F = FractionField(R)
1447        x,y,T = R.gens()
1448        d = self.minimum_distance()
1449        n = self.length()
1450        k = self.dimension()
1451        dperp = (self.dual_code()).minimum_distance()
1452        if d == 1 or dperp == 1:
1453            print "\n WARNING: There is no guarantee this function works when the minimum distance\n"
1454            print "            of the code or of the dual code equals 1.\n"
1455        gammaC = n+1-d #  C.genus()
1456        if mode=="dual":
1457            Cp = self.dual_code()
1458            #dp = Cp.minimum_distance()
1459            #kp = Cp.dimension()
1460            gammaCp = Cp.genus() # = n+1-kp-dp
1461        q = self.characteristic()
1462        W = self.weight_distribution()
1463        A = add([W[i]*y**i*x**(n-i) for i in range(n+1)])
1464        f = (x*T+y*(1-T))**n*add([T**i for i in range(n+3)])*add([(q*T)**i for i in range(n+3)])
1465        Mn = [f.coefficient(T**(n-i)) for i in range(d,n)]
1466        coeffs = [[(Mn[j]+x**(n)).monomial_coefficient(x**i*y**(n-i)) for i in range(n+1)] for j in range(n-d)]
1467        MS = MatrixSpace(QQ,n-d,n+1)
1468        M = MS(coeffs)
1469        F = PowerSeriesRing(PolynomialRing(QQ,2,"xy"),"T")
1470        for m in range(2,n-d+1):
1471            #x,y,T = F.gens()
1472            V = VectorSpace(QQ,m)
1473            v = V([W[i] for i in range(m)])
1474            Mmm = M.matrix_from_rows_and_columns(range(m),range(m))
1475            if Mmm.determinant()!=0:
1476                Pcoeffs = v*Mmm**(-1)
1477                P = add([Pcoeffs[i]*T**i for i in range(len(Pcoeffs))])
1478                #x,y,T = F.gens()
1479                Z = P*(1-T)**(-1)*(1-q*T)**(-1)
1480                lhs = P*f
1481                r = lhs.coefficient(T**(n-d))/((A-x**n)/(q-1))
1482                if mode=="verbose": print m, P, r
1483                #if r(x,y,T)==r(1,y,T) and mode=="dual":         ## check if r is a constant in x ...
1484                #    Z = F(Z*r(1,1,1)**(-1))
1485                #    x,y,T = F.gens()
1486                #    Zp = Z(x,y,1/(q*T))*q**(gammaC-1)*T**(gammaC+gammaCp-2)
1487                #    return Z,Zp
1488                if r(x,y,T)==r(1,y,T) and mode=="normalized":   ## check if r is a constant in x ...
1489                    Z = Z*r(1,1,1)**(-1)
1490                    #x,y,T = F.gens()
1491                    return Z*T**(1-gammaC)
1492                if r(x,y,T)==r(1,y,T):                         ## check if r is a constant in x ...
1493                    Z = Z*r(1,1,1)**(-1)
1494                    return Z
1495                if mode=="verbose":
1496                    print "det = ",Mmm.determinant(),"   m = ",m,"  Z = ",Z
1497            #if Mmm.determinant()==0:
1498            #    print "det = ",Mmm.determinant(),"   m = ",m
1499        return "fails"
1500
1501    def zeta_function3(self,mode=None):
1502        r"""
1503        Returns the Duursma zeta function of the code.
1504
1505        NOTE: This sometimes returns "fail" for a reason I don't fully
1506        understand. However, when it does return a polynomial, the answer is
1507        (as far as I know) correct. *Experimental code* included to
1508        study a particular implementation.
1509
1510        INPUT:
1511            mode -- string
1512               mode = "dual" computes both the zeta function of C and that of C*
1513               mode = "normalized" computes the normalized zeta function of C, zeta_C(T)*T(1-genus(C))
1514            (NOTE: if xi(T,C) denotes the normalized zeta function
1515            then xi(T,C*) = xi(1/(qT),C) is equivalent to zeta_{C*}(T)
1516            = zeta_C(1/(qT))*q^(gamma-1)T^(gamma+gamma*-2), where
1517            gamma = gamma(C) and gamma* = gamma(C*).)
1518
1519        EXAMPLES:
1520            sage.: C = HammingCode(3,GF(2))
1521            sage.: C.zeta_function3()
1522            (2/5*T^2 + 2/5*T + 1/5)/(2*T^2 - 3*T + 1)
1523            sage.: C = ExtendedTernaryGolayCode()
1524            sage.: C.zeta_function3()
1525            (3/7*T^2 + 3/7*T + 1/7)/(3*T^2 - 4*T + 1)
1526            sage.: C.zeta_function3(mode="dual")
1527            ((3/7*T^2 + 3/7*T + 1/7)/(3*T^2 - 4*T + 1),
1528             (729/7*T^2 + 729/7*T + 243/7)/(729*T^2 - 972*T + 243))
1529
1530        REFERENCES:
1531            I. Duursma, "Weight distributions of geometric Goppa codes," Trans. A.M.S., 351 (1999)3609-3639.
1532        """
1533        R = PolynomialRing(QQ,3,"xyT")
1534        F = FractionField(R)
1535        x,y,T = R.gens()
1536        d = self.minimum_distance()
1537        n = self.length()
1538        k = self.dimension()
1539        dperp = (self.dual_code()).minimum_distance()
1540        if d == 1 or dperp == 1:
1541            print "\n WARNING: There is no guarantee this function works when the minimum distance\n"
1542            print "            of the code or of the dual code equals 1.\n"
1543        gammaC = n+1-k-d
1544        if mode=="dual":
1545            Cp = self.dual_code()
1546            dp = Cp.minimum_distance()
1547            kp = Cp.dimension()
1548            gammaCp = n+1-kp-dp
1549        q = self.characteristic()
1550        W = self.weight_distribution()
1551        A = add([W[i]*y**i*x**(n-i) for i in range(n+1)])
1552        f = (x*T+y*(1-T))**n*add([T**i for i in range(n+3)])*add([(q*T)**i for i in range(n+3)])
1553        Mn = [f.coefficient(T**(n-i)) for i in range(d,n)]
1554        coeffs = [[(Mn[j]+x**(n)).coefficient(x**i*y**(n-i))(1,1,1) for i in range(n+1)] for j in range(n-d)]
1555        MS = MatrixSpace(QQ,n-d,n+1)
1556        #print [type(x[0]) for x in coeffs], MS
1557        M = MS(coeffs)
1558        for m in range(2,n-d+1):
1559            V = VectorSpace(QQ,m)
1560            v = V([W[i] for i in range(m)])
1561            Mmm = M.matrix_from_rows_and_columns(range(m),range(m))
1562            if Mmm.determinant()!=0:
1563                Pcoeffs = v*Mmm**(-1)
1564                P = add([Pcoeffs[i]*T**i for i in range(len(Pcoeffs))])
1565                Z = P*(1-T)**(-1)*(1-q*T)**(-1)
1566                lhs = P*f
1567                r = lhs.coefficient(T**(n-d))/((A-x**n)/(q-1))
1568                if mode=="verbose": print m, P, r
1569                if r(x,y,T)==r(x+T,y,T) and mode=="dual":
1570                    Z = Z*r(1,1,1)**(-1)
1571                    x,y,T = F.gens()
1572                    Zp = Z(x,y,1/(q*T))*q**(gammaC-1)*T**(gammaC+gammaCp-2)
1573                    return Z,Zp
1574                if r(x,y,T)==r(x+T,y,T) and mode=="normalized":
1575                    Z = Z*r(1,1,1)**(-1)
1576                    x,y,T = F.gens()
1577                    Zp = Z(x,y,1/(q*T))*q**(gammaC-1)*T**(gammaC+gammaCp-2)
1578                    return Z*T**(1-gammaC)
1579                if r(x,y,T)==r(1,y,T): # check if r is a constant in x ...
1580                    Z = Z*r(1,1,1)**(-1)
1581                    return Z
1582                if mode=="verbose":
1583                    print "det = ",Mmm.determinant(),"   m = ",m,"  Z = ",Z
1584            if Mmm.determinant()==0:
1585                pass
1586                #print "det = ",Mmm.determinant(),"   m = ",m
1587        return "fails"
1588
1589    def punctured(self,L):
1590        r"""
1591        Returns the code punctured at the positions L, $L \subset \{1,2,...,n\}$.
1592        If C is a code of length n in GF(q) then the code $C^L$ obtained from C
1593        by puncturing at the positions in L is the code of length n-|L|
1594        consisting of codewords of $C$ which have their $i-th$ coordinate
1595        deleted if $i \in L$ and left alone if $i\notin L$:
1596        $$1597 C^L = \{(c_{i_1},...,c_{i_N})\ |\ (c_1,...,c_n)\in C\}, 1598$$
1599        where $\{1,2,...,n\}-T = \{i_1,...,i_N\}$. In particular, if $L=\{j\}$ then
1600        $C^L$ is simply the code obtained from $C$ by deleting the $j-th$ coordinate
1601        of each codeword. The code $C^L$ is called the {\it punctured code at $L$}.
1602        The dimension of $C^L$ can decrease if $|L|>d-1$.
1603
1604        EXAMPLES:
1605            sage: C = HammingCode(3,GF(2))
1606            sage: C.punctured([1,2])
1607            Linear code of length 5, dimension 4 over Finite Field of size 2
1608
1609        """
1610        G = self.gen_mat()
1611        F = self.base_ring()
1612        q = F.order()
1613        n = len(G.columns())
1614        k = len(G.rows())
1615        nL = n-len(L)
1616        Gstr = str(gap(G))
1617        gap.eval( "G:="+Gstr )
1618        C = gap("GeneratorMatCode(G,GF("+str(q)+"))")
1619        CP = C.PuncturedCode(L)
1620        Gp = CP.GeneratorMat()
1621        kL = Gp.Length()  ## this is = dim(CL), usually = k
1622        G2 = Gp._matrix_(F)
1623        MS = MatrixSpace(F,kL,nL)
1624        return LinearCode(MS(G2))
1625
1626    def shortened(self,L):
1627        r"""
1628        Returns the code shortened at the positions L, $L \subset \{1,2,...,n\}$.
1629        Consider the subcode $C(L)$ consisting of all codewords $c\in C$ which
1630        satisfy $c_i=0$ for all $i\in L$. The punctured code $C(L)^L$ is called
1631        the {\it shortened code on $L$} and is denoted $C_L$. The code
1632        constructed is actually only isomorphic to the shortened code defined
1633        in this way.
1634
1635        EXAMPLES:
1636            sage: C = HammingCode(3,GF(2))
1637            sage: C.shortened([1,2])
1638            Linear code of length 5, dimension 2 over Finite Field of size 2
1639
1640        """
1641        G = self.gen_mat()
1642        F = self.base_ring()
1643        q = F.order()
1644        n = len(G.columns())
1645        k = len(G.rows())
1646        nL = n-len(L)
1647        Gstr = str(gap(G))
1648        gap.eval("G:="+Gstr )
1649        C = gap("GeneratorMatCode(G,GF("+str(q)+"))")
1650        Cd = C.DualCode()
1651        Cdp = Cd.PuncturedCode(L)
1652        Cdpd = Cdp.DualCode()
1653        Gs = Cdpd.GeneratorMat()
1654        kL = Gs.Length()  ## this is = dim(CL), usually = k
1655        Gss = Gs._matrix_(F)
1656        MS = MatrixSpace(F,kL,nL)
1657        return LinearCode(MS(Gss))
1658
1659    def binomial_moment(self,i):
1660        r"""
1661        Returns the i-th binomial moment of the $[n,k,d]_q$-code $C$:
1662        $1663 B_i(C) = \sum_{S, |S|=i} \frac{q^{k_S}-1}{q-1} 1664$
1665        where $k_S$ is the dimension of the shortened code $C_{J-S}$, $J=[1,2,...,n]$.
1666        (The normalized binomial moment is $b_i(C) = binomial(n,d+i)^{-1}B_{d+i}(C)$.)
1667        In other words, $C_{J-S}$ is the isomorphic to the subcode of C of codewords
1668        supported on S.
1669
1670        EXAMPLES:
1671            sage: C = HammingCode(3,GF(2))
1672            sage: C.binomial_moment(2)
1673            0
1674            sage: C.binomial_moment(3)    # long time
1675            0
1676            sage: C.binomial_moment(4)    # long time
1677            35
1678
1679        WARNING: This is slow.
1680
1681        REFERENCE:
1682            I. Duursma, "Combinatorics of the two-variable zeta function",
1683            Finite fields and applications, 109--136, Lecture Notes in Comput. Sci.,
1684            2948, Springer, Berlin, 2004.
1685        """
1686        n = self.length()
1687        #ii = n-i
1688        k = self.dimension()
1689        d = self.minimum_distance()
1690        F = self.base_ring()
1691        q = F.order()
1692        J = range(1,n+1)
1693        Cp = self.dual_code()
1694        dp = Cp.minimum_distance()
1695        if i<d:
1696            return 0
1697        if i>n-dp and i<=n:
1698            return binomial(n,i)*(q**(i+k-n) -1)/(q-1)
1699        P = SetPartitions(J,2).list()
1700        b = QQ(0)
1701        for p in P:
1702            p = list(p)
1703            S = p[0]
1704            if len(S)==n-i:
1705                C_S = self.shortened(S)
1706                k_S = C_S.dimension()
1707                b = b + (q**(k_S) -1)/(q-1)
1708        return b
1709
1710    def divisor(self):
1711        r"""
1712        Returns the divisor of a code (the divisor is the smallest
1713        integer $d_0>0$ such that each $A_i>0$ iff $i$ is divisible by $d_0$).
1714
1715        EXAMPLES:
1716            sage: C = ExtendedBinaryGolayCode()
1717            sage: C.divisor()   ## type 2
1718            4
1720            sage: C.divisor()   ## type 1
1721            2
1722
1723        """
1724        C = self
1725        A = C.spectrum()
1726        n = C.length()
1727        V = VectorSpace(QQ,n+1)
1728        S = V(A).nonzero_positions()
1729        S0 = [S[i] for i in range(1,len(S))]
1730        #print S0
1731        if len(S)>1: return GCD(S0)
1732        return 1
1733
1734    def support(self):
1735        r"""
1736        Returns the set of indices $j$ where $A_j$ is nonzero,
1737        where spectrum(self) = $[A_0,A_1,...,A_n]$.
1738
1739        EXAMPLES:
1740            sage: C = HammingCode(3,GF(2))
1741            sage: C.spectrum()
1742            [1, 0, 0, 7, 7, 0, 0, 1]
1743            sage: C.support()
1744            [0, 3, 4, 7]
1745
1746        """
1747        n = self.length()
1748        F = self.base_ring()
1749        V = VectorSpace(F,n+1)
1750        return V(self.spectrum()).support()
1751
1752    def characteristic_polynomial(self):
1753        r"""
1754        Returns the characteristic polynomial of a linear code, as
1755        defined in van Lint's text [vL].
1756
1757        EXAMPLES:
1759            sage: C.characteristic_polynomial()
1760            -2*x + 16
1761
1762        REFERENCES:
1763            van Lint, {\it Introduction to coding theory, 3rd ed.}, Springer-Verlag GTM, 86, 1999.
1764
1765        """
1766        R = PolynomialRing(QQ,"x")
1767        x = R.gen()
1768        C = self
1769        Cd = C.dual_code()
1770        Sd = Cd.support()
1771        k = C.dimension()
1772        n = C.length()
1773        q = (C.base_ring()).order()
1774        return q**(n-k)*prod([1-x/j for j in Sd if j>0])
1775
1776    def sd_duursma_data(C,i):
1777        r"""
1778        INPUT:
1779           The formally s.d. code C and the type number (1,2,3,4)
1780           (does not check if C is actually sd)
1781
1782        RETURN:
1783           The data v,m as in Duursama [D]
1784
1785        EXAMPLES:
1786
1787        REFERENCES:
1788            [D] I. Duursma, "Extremal weight enumerators and ultraspherical polynomials"
1789
1790        """
1791        n = C.length()
1792        d = C.minimum_distance()
1793        if i == 1:
1794            v = (n-4*d)/2 + 4
1795            m = d-3
1796        if i == 2:
1797            v = (n-6*d)/8 + 3
1798            m = d-5
1799        if i == 3:
1800            v = (n-4*d)/4 + 3
1801            m = d-4
1802        if i == 4:
1803            v = (n-3*d)/2 + 3
1804            m = d-3
1805        return [v,m]
1806
1807    def sd_duursma_q(C,i,d0):
1808        r"""
1809
1810        INPUT:
1811           C -- an sd code (does not check if C is actually an sd code),
1812           i -- the type number, one of 1,2,3,4,
1813           d0 -- and the divisor d0 (the smallest integer d0>0 such that each
1814                 A_i>0 iff i is divisible by d0).
1815
1816        RETURN:
1817           The coefficients $q_0, q_1, ...,$  of $q(T)$ as in Duursama [D].
1818
1819        EXAMPLES:
1820
1821        REFERENCES:
1822            [D] I. Duursma, "Extremal weight enumerators and ultraspherical polynomials"
1823
1824        """
1825        q = (C.base_ring()).order()
1826        n = C.length()
1827        d = C.minimum_distance()
1828        d0 = C.divisor()
1829        if i==1 or i==2:
1830            if d>d0:
1831                c0 = QQ((n-d)*rising_factorial(d-d0,d0+1)*C.spectrum()[d])/rising_factorial(n-d0-1,d0+2)
1832            else:
1833                c0 = QQ((n-d)*C.spectrum()[d])/rising_factorial(n-d0-1,d0+2)
1834        if i==3 or i==4:
1835            if d>d0:
1836                c0 = rising_factorial(d-d0,d0+1)*C.spectrum()[d]/((q-1)*rising_factorial(n-d0,d0+1))
1837            else:
1838                c0 = C.spectrum()[d]/((q-1)*rising_factorial(n-d0,d0+1))
1839        v = ZZ(C.sd_duursma_data(i)[0])
1840        m = ZZ(C.sd_duursma_data(i)[1])
1841        #print m,v,d,d0,c0
1842        if m<0 or v<0:
1843            raise ValueError("This case not implemented.")
1844        PR = PolynomialRing(QQ,"T")
1845        T = PR.gen()
1846        if i == 1:
1847            coefs = PR(c0*(1+3*T+2*T**2)**m*(2*T**2+2*T+1)**v).list()
1848            #print coefs, len(coefs)
1849            qc = [coefs[j]/binomial(4*m+2*v,m+j) for j in range(2*m+2*v+1)]
1850            q = PR(qc)
1851        if i == 2:
1852            F = ((T+1)**8+14*T**4*(T+1)**4+T**8)**v
1853            coefs = (c0*(1+T)**m*(1+4*T+6*T**2+4*T**3)**m*F).coeffs()
1854            qc = [coefs[j]/binomial(6*m+8*v,m+j) for j in range(4*m+8*v+1)]
1855            q = PR(qc)
1856        if i == 3:
1857            F = (3*T^2+4*T+1)**v*(1+3*T^2)**v
1858            ## Note that: (3*T^2+4*T+1)(1+3*T^2)=(T+1)**4+8*T**3*(T+1)
1859            coefs = (c0*(1+3*T+3*T**2)**m*F).coeffs()
1860            qc = [coefs[j]/binomial(4*m+4*v,m+j) for j in range(2*m+4*v+1)]
1861            q = PR(qc)
1862        if i == 4:
1863            coefs = (c0*(1+2*T)**m*(4*T**2+2*T+1)**v).coeffs()
1864            qc = [coefs[j]/binomial(3*m+2*v,m+j) for j in range(m+2*v+1)]
1865            q = PR(qc)
1866        return q/q(1)
1867
1868    def sd_zeta_polynomial(C,type=1):
1869        r"""
1870        Returns the Duursma zeta function of a self-dual code using
1871        the construction in [D].
1872
1873        INPUT:
1874            type -- type of the s.d. code; one of 1,2,3, or 4.
1875
1876        EXAMPLES:
1878            sage: C.sd_zeta_polynomial()
1879            2/5*T^2 + 2/5*T + 1/5
1880            sage: C.zeta_function()
1881            (2/5*T^2 + 2/5*T + 1/5)/(2*T^2 - 3*T + 1)
1883            sage: P = C.sd_zeta_polynomial(); P
1884            8/91*T^8 + 16/91*T^7 + 212/1001*T^6 + 28/143*T^5 + 43/286*T^4 + 14/143*T^3 + 53/1001*T^2 + 2/91*T + 1/182
1885            sage: P(1)
1886            1
1887            sage: C.sd_zeta_polynomial()
1888            8/91*T^8 + 16/91*T^7 + 212/1001*T^6 + 28/143*T^5 + 43/286*T^4 + 14/143*T^3 + 53/1001*T^2 + 2/91*T + 1/182
1889
1890        It is a general fact about Duursma zeta polynomials that P(1) = 1.
1891
1892        REFERENCES:
1893            [D] I. Duursma, "Extremal weight enumerators and ultraspherical polynomials"
1894        """
1895        d0 = C.divisor()
1896        q = (C.base_ring()).order()
1897        P = C.sd_duursma_q(type,d0)
1898        PR = P.parent()
1899        T = FractionField(PR).gen()
1900        if type == 1: return P
1901        if type == 2: return P/(1-2*T+2*T**2)
1902        if type == 3: return P/(1+3*T**2)
1903        if type == 4: return P/(1+2*T)
1904
1905def LinearCodeFromVectorSpace(self):
1906    """
1907    """
1908    F = self.base_ring()
1909    B = self.basis()
1910    n = len(B[0].list())
1911    k = len(B)
1912    MS = MatrixSpace(F,k,n)
1913    G = MS([B[i].list() for i in range(k)])
1914    return LinearCode(G)
Note: See TracBrowser for help on using the repository browser.