Ticket #3112: trac-3112-rebased.patch

File trac-3112-rebased.patch, 84.4 KB (added by rlm, 13 years ago)

This should replace the previous two.

  • sage/coding/all.py

    # HG changeset patch
    # User Robert Miller <rlm@rlmiller.org>
    # Date 1211495827 25200
    # Node ID 8bcb24714f113ab57bdbbca78c924b73884fcacc
    # Parent  f2c1d993c672d751f09130924124fe368626e4bd
    Flattened patch for 3112.
    
    diff -r f2c1d993c672 -r 8bcb24714f11 sage/coding/all.py
    a b from linear_code import (LinearCode, Lin 
    5050                         hamming_weight,
    5151                         best_known_linear_code,
    5252                         best_known_linear_code_www,
    53                          bounds_minimum_distance)
     53                         bounds_minimum_distance,
     54                         self_orthogonal_binary_codes)
    5455
    5556from sd_codes import self_dual_codes_binary
    5657
  • sage/coding/binary_code.pxd

    diff -r f2c1d993c672 -r 8bcb24714f11 sage/coding/binary_code.pxd
    a b cdef int *hamming_weights() 
    44cdef int *hamming_weights()
    55
    66ctypedef unsigned int codeword
     7
     8cdef struct WordPermutation:
     9    # A word permutation is a permutation of the bits in a codeword.
     10    int chunk_num
     11    int chunk_words
     12    int degree
     13    codeword **images
     14    codeword gate
    715
    816cdef class BinaryCode:
    917    cdef codeword *basis
    cdef class BinaryCode: 
    1523
    1624    cdef int is_one(self, int, int)
    1725    cdef int is_automorphism(self, int *, int *)
     26    cpdef int put_in_std_form(self)
     27    cdef void _apply_permutation_to_basis(self, object labeling)
     28    cdef void _update_words_from_basis(self)
     29
     30cdef WordPermutation *create_word_perm(object)
     31cdef WordPermutation *create_array_word_perm(int *, int, int)
     32cdef WordPermutation *create_id_word_perm(int)
     33cdef WordPermutation *create_comp_word_perm(WordPermutation *, WordPermutation *)
     34cdef WordPermutation *create_inv_word_perm(WordPermutation *)
     35cdef int dealloc_word_perm(WordPermutation *)
     36cdef codeword permute_word_by_wp(WordPermutation *, codeword)
     37cdef codeword *expand_to_ortho_basis(BinaryCode, int)
    1838
    1939cdef class OrbitPartition:
    2040    cdef int nwords
    cdef class PartitionStack: 
    6888    cdef int refine(self, int, int *, int, BinaryCode, int *)
    6989    cdef void clear(self, int)
    7090    cdef int cmp(self, PartitionStack, BinaryCode)
    71     cdef void find_basis(self, int *)
    72     cdef void get_permutation(self, PartitionStack, int *, int *, int *)
     91    cdef int find_basis(self, int *)
     92    cdef void get_permutation(self, PartitionStack, int *, int *)
    7393
    7494cdef class BinaryCodeClassifier:
    7595    cdef int *ham_wts
    cdef class BinaryCodeClassifier: 
    84104    cdef int *alpha
    85105    cdef int alpha_size
    86106    cdef int *v, *e
    87     cdef int *aut_gp_gens, *labeling
    88     cdef int aut_gp_index, aut_gens_size
     107    cdef int *aut_gp_gens, *labeling, *base
     108    cdef int aut_gp_index, aut_gens_size, base_size
    89109    cdef object aut_gp_size
    90110
    91111    cdef int Phi_size
    cdef class BinaryCodeClassifier: 
    95115
    96116
    97117
    98 
    99 
    100 
    101 
    102 
  • sage/coding/binary_code.pyx

    diff -r f2c1d993c672 -r 8bcb24714f11 sage/coding/binary_code.pyx
    a b AUTHOR: 
    2020        * union-find based orbit partition
    2121        * optimized partition stack class
    2222        * NICE-based partition refinement algorithm
    23 
    2423        * canonical generation function
    2524
    2625"""
    include '../ext/cdefs.pxi' 
    3534include '../ext/cdefs.pxi'
    3635include '../ext/python_mem.pxi'
    3736include '../ext/stdsage.pxi'
     37include '../ext/interrupt.pxi'
    3838from sage.structure.element import is_Matrix
    3939from sage.misc.misc import cputime
    40 from math import log, floor
    41 from sage.rings.integer import Integer
     40from sage.rings.integer cimport Integer
     41from copy import copy
    4242
    4343WORD_SIZE = sizeof(codeword) << 3
     44
     45cdef enum:
     46    chunk_size = 8
     47
     48cdef inline int min(int a, int b):
     49    if a > b:
     50        return b
     51    else:
     52        return a
    4453
    4554## NOTE - Since most of the functions are used from within the module, cdef'd
    4655## functions come without an underscore, and the def'd equivalents, which are
    4756## essentially only for doctesting and debugging, have underscores.
    4857
    4958cdef int *hamming_weights():
    50     cdef int *ham_wts
     59    cdef int *ham_wts, i
    5160    ham_wts = <int *> sage_malloc( 65536 * sizeof(int) )
    52     if not ham_wts:
     61    if ham_wts is NULL:
    5362        sage_free(ham_wts)
    5463        raise MemoryError("Memory.")
    5564    ham_wts[0] = 0
    cdef int *hamming_weights(): 
    6372    for i from 256 <= i < 65536:
    6473        ham_wts[i] = ham_wts[i & 255] + ham_wts[(i>>8) & 255]
    6574    return ham_wts
     75
     76def test_word_perms(t_limit=5.0):
     77    """
     78    Tests the WordPermutation structs for at least t_limit seconds.
     79   
     80    These are structures written in pure C for speed, and are tested from this
     81    function, which performs the following tests:
     82   
     83    1. Tests create_word_perm, which creates a WordPermutation from a Python
     84        list L representing a permutation i --> L[i]. Takes a random word and
     85        permutes it by a random list permutation, and tests that the result
     86        agrees with doing it the slow way.
     87   
     88    1b. Tests create_array_word_perm, which creates a WordPermutation from a
     89        C array. Does the same as above.
     90   
     91    2. Tests create_comp_word_perm, which creates a WordPermutation as a
     92        composition of two WordPermutations. Takes a random word and
     93        two random permutations, and tests that the result of permuting by the
     94        composition is correct.
     95   
     96    3. Tests create_inv_word_perm and create_id_word_perm, which create a
     97        WordPermutation as the inverse and identity permutations, resp.
     98        Takes a random word and a random permutation, and tests that the result
     99        permuting by the permutation and its inverse in either order, and
     100        permuting by the identity both return the original word.
     101   
     102    NOTE:
     103        The functions permute_word_by_wp and dealloc_word_perm are implicitly
     104        involved in each of the above tests.
     105   
     106    TESTS:
     107        sage: from sage.coding.binary_code import test_word_perms
     108        sage: test_word_perms()
     109
     110    """
     111    cdef WordPermutation *g, *h, *i
     112    cdef codeword cw1, cw2, cw3
     113    cdef int n = sizeof(codeword) << 3
     114    cdef int j, *arr = <int*> sage_malloc(n * sizeof(int))
     115    if arr is NULL:
     116        raise MemoryError("Error allocating memory.")
     117    from sage.misc.prandom import randint
     118    from sage.combinat.permutation import Permutations
     119    S = Permutations(range(n))
     120    t = cputime()
     121    while cputime(t) < t_limit:
     122        word = [randint(0,1) for _ in xrange(n)]
     123        cw1 = 0
     124        for j from 0 <= j < n:
     125            cw1 += (<codeword>word[j]) << (<codeword>j)
     126        # 1. test create_word_perm
     127        gg = S.random_element()
     128        g = create_word_perm(gg)
     129        word2 = [0]*n
     130        for j from 0 <= j < n:
     131            word2[gg[j]] = word[j]
     132        cw2 = permute_word_by_wp(g, cw1)
     133        cw3 = 0
     134        for j from 0 <= j < n:
     135            cw3 += (<codeword>word2[j]) << (<codeword>j)
     136        if cw3 != cw2:
     137            print "ERROR1"
     138        dealloc_word_perm(g)
     139        # 1b. test create_array_word_perm
     140        gg = S.random_element()
     141        for j from 0 <= j < n:
     142            arr[j] = gg[j]
     143        g = create_array_word_perm(arr, 0, n)
     144        word2 = [0]*n
     145        for j from 0 <= j < n:
     146            word2[gg[j]] = word[j]
     147        cw2 = permute_word_by_wp(g, cw1)
     148        cw3 = 0
     149        for j from 0 <= j < n:
     150            cw3 += (<codeword>word2[j]) << (<codeword>j)
     151        if cw3 != cw2:
     152            print "ERROR1b"
     153        dealloc_word_perm(g)
     154        # 2. test create_comp_word_perm
     155        gg = S.random_element()
     156        hh = S.random_element()
     157        g = create_word_perm(gg)
     158        h = create_word_perm(hh)
     159        i = create_comp_word_perm(g, h)
     160        word2 = [0]*n
     161        for j from 0 <= j < n:
     162            word2[gg[hh[j]]] = word[j]
     163        cw2 = permute_word_by_wp(i, cw1)
     164        cw3 = 0
     165        for j from 0 <= j < n:
     166            cw3 += (<codeword>word2[j]) << (<codeword>j)
     167        if cw3 != cw2:
     168            print "ERROR2"
     169        dealloc_word_perm(g)
     170        dealloc_word_perm(h)
     171        dealloc_word_perm(i)
     172        # 3. test create_inv_word_perm and create_id_word_perm
     173        gg = S.random_element()
     174        g = create_word_perm(gg)
     175        h = create_inv_word_perm(g)
     176        i = create_id_word_perm(n)
     177        cw2 = permute_word_by_wp(g, cw1)
     178        cw2 = permute_word_by_wp(h, cw2)
     179        if cw1 != cw2:
     180            print "ERROR3a"
     181        cw2 = permute_word_by_wp(h, cw1)
     182        cw2 = permute_word_by_wp(g, cw2)
     183        if cw1 != cw2:
     184            print "ERROR3b"
     185        cw2 = permute_word_by_wp(i, cw1)
     186        if cw1 != cw2:
     187            print "ERROR3c"
     188        dealloc_word_perm(g)
     189        dealloc_word_perm(h)
     190
     191cdef WordPermutation *create_word_perm(object list_perm):
     192    r"""
     193    Create a word permutation from a Python list permutation L, i.e. such that
     194    $i \mapsto L[i]$.
     195    """
     196    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     197    cdef codeword *images_i, image
     198    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     199    if word_perm is NULL:
     200        raise RuntimeError("Error allocating memory.")
     201    word_perm.degree = len(list_perm)
     202    list_perm = copy(list_perm)
     203    while num_chunks*chunk_size < word_perm.degree:
     204        num_chunks += 1
     205    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     206    if word_perm.images is NULL:
     207        sage_free(word_perm)
     208        raise RuntimeError("Error allocating memory.")
     209    word_perm.chunk_num = num_chunks
     210    words_per_chunk = 1 << chunk_size
     211    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     212    list_perm += range(len(list_perm), chunk_size*num_chunks)
     213    word_perm.chunk_words = words_per_chunk
     214    for i from 0 <= i < num_chunks:
     215        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     216        if images_i is NULL:
     217            for j from 0 <= j < i:
     218                sage_free(word_perm.images[j])
     219            sage_free(word_perm.images)
     220            sage_free(word_perm)
     221            raise RuntimeError("Error allocating memory.")           
     222        word_perm.images[i] = images_i
     223        for j from 0 <= j < chunk_size:
     224            images_i[1 << j] = (<codeword>1) << list_perm[chunk_size*i + j]
     225        image = <codeword> 0
     226        parity = 0
     227        comb = 0
     228        while 1:
     229            images_i[comb] = image
     230            parity ^= 1
     231            j = 0
     232            if not parity:
     233                while not comb & (1 << j): j += 1
     234                j += 1
     235            if j == chunk_size: break
     236            else:
     237                comb ^= (1 << j)
     238                image ^= images_i[1 << j]
     239    return word_perm
     240
     241cdef WordPermutation *create_array_word_perm(int *array, int start, int degree):
     242    """
     243    Create a word permutation of a given degree from a C array, starting at start.
     244    """
     245    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     246    cdef codeword *images_i, image
     247    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     248    if word_perm is NULL:
     249        raise RuntimeError("Error allocating memory.")
     250    word_perm.degree = degree
     251    while num_chunks*chunk_size < word_perm.degree:
     252        num_chunks += 1
     253    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     254    if word_perm.images is NULL:
     255        sage_free(word_perm)
     256        raise RuntimeError("Error allocating memory.")
     257    word_perm.chunk_num = num_chunks
     258    words_per_chunk = 1 << chunk_size
     259    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     260    word_perm.chunk_words = words_per_chunk
     261    for i from 0 <= i < num_chunks:
     262        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     263        if images_i is NULL:
     264            for j from 0 <= j < i:
     265                sage_free(word_perm.images[j])
     266            sage_free(word_perm.images)
     267            sage_free(word_perm)
     268            raise RuntimeError("Error allocating memory.")           
     269        word_perm.images[i] = images_i
     270        for j from 0 <= j < chunk_size:
     271            images_i[1 << j] = (<codeword>1) << array[start + chunk_size*i + j]
     272        image = <codeword> 0
     273        parity = 0
     274        comb = 0
     275        while 1:
     276            images_i[comb] = image
     277            parity ^= 1
     278            j = 0
     279            if not parity:
     280                while not comb & (1 << j): j += 1
     281                j += 1
     282            if j == chunk_size: break
     283            else:
     284                comb ^= (1 << j)
     285                image ^= images_i[1 << j]
     286    return word_perm
     287
     288cdef WordPermutation *create_id_word_perm(int degree):
     289    """
     290    Create the identity word permutation of degree degree.
     291    """
     292    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     293    cdef codeword *images_i, image
     294    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     295    if word_perm is NULL:
     296        raise RuntimeError("Error allocating memory.")
     297    word_perm.degree = degree
     298    while num_chunks*chunk_size < degree:
     299        num_chunks += 1
     300    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     301    if word_perm.images is NULL:
     302        sage_free(word_perm)
     303        raise RuntimeError("Error allocating memory.")
     304    word_perm.chunk_num = num_chunks
     305    words_per_chunk = 1 << chunk_size
     306    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     307    word_perm.chunk_words = words_per_chunk
     308    for i from 0 <= i < num_chunks:
     309        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     310        if images_i is NULL:
     311            for j from 0 <= j < i:
     312                sage_free(word_perm.images[j])
     313            sage_free(word_perm.images)
     314            sage_free(word_perm)
     315            raise RuntimeError("Error allocating memory.")
     316        word_perm.images[i] = images_i
     317        for j from 0 <= j < chunk_size:
     318            images_i[1 << j] = (<codeword>1) << (chunk_size*i + j)
     319        image = <codeword> 0
     320        parity = 0
     321        comb = 0
     322        while 1:
     323            images_i[comb] = image
     324            parity ^= 1
     325            j = 0
     326            if not parity:
     327                while not comb & (1 << j): j += 1
     328                j += 1
     329            if j == chunk_size: break
     330            else:
     331                comb ^= (1 << j)
     332                image ^= images_i[1 << j]
     333    return word_perm
     334
     335cdef WordPermutation *create_comp_word_perm(WordPermutation *g, WordPermutation *h):
     336    r"""
     337    Create the composition of word permutations $g \circ h$.
     338    """
     339    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     340    cdef codeword *images_i, image
     341    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     342    if word_perm is NULL:
     343        raise RuntimeError("Error allocating memory.")
     344    word_perm.degree = g.degree
     345    while num_chunks*chunk_size < word_perm.degree:
     346        num_chunks += 1
     347    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     348    if word_perm.images is NULL:
     349        sage_free(word_perm)
     350        raise RuntimeError("Error allocating memory.")
     351    word_perm.chunk_num = num_chunks
     352    words_per_chunk = 1 << chunk_size
     353    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     354    word_perm.chunk_words = words_per_chunk
     355    for i from 0 <= i < num_chunks:
     356        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     357        if images_i is NULL:
     358            for j from 0 <= j < i:
     359                sage_free(word_perm.images[j])
     360            sage_free(word_perm.images)
     361            sage_free(word_perm)
     362            raise RuntimeError("Error allocating memory.")           
     363        word_perm.images[i] = images_i
     364        for j from 0 <= j < chunk_size:
     365            image = (<codeword>1) << (chunk_size*i + j)
     366            image = permute_word_by_wp(h, image)
     367            image = permute_word_by_wp(g, image)
     368            images_i[1 << j] = image
     369        image = <codeword> 0
     370        parity = 0
     371        comb = 0
     372        while 1:
     373            images_i[comb] = image
     374            parity ^= 1
     375            j = 0
     376            if not parity:
     377                while not comb & (1 << j): j += 1
     378                j += 1
     379            if j == chunk_size: break
     380            else:
     381                comb ^= (1 << j)
     382                image ^= images_i[1 << j]
     383    return word_perm
     384
     385cdef WordPermutation *create_inv_word_perm(WordPermutation *g):
     386    r"""
     387    Create the inverse $g^{-1}$ of the word permutation of $g$.
     388    """
     389    cdef int i, j, *array = <int *> sage_malloc( g.degree * sizeof(int) )
     390    cdef codeword temp
     391    cdef WordPermutation *w
     392    for i from 0 <= i < g.degree:
     393        j = 0
     394        temp = permute_word_by_wp(g, (<codeword>1) << i)
     395        while not ((<codeword>1) << j) & temp:
     396            j += 1
     397        array[j] = i
     398    w = create_array_word_perm(array, 0, g.degree)
     399    sage_free(array)
     400    return w
     401
     402cdef int dealloc_word_perm(WordPermutation *wp):
     403    """
     404    Free the memory used by a word permutation.
     405    """
     406    cdef int i
     407    for i from 0 <= i < wp.chunk_num:
     408        sage_free(wp.images[i])
     409    sage_free(wp.images)
     410    sage_free(wp)
     411
     412cdef codeword permute_word_by_wp(WordPermutation *wp, codeword word):
     413    """
     414    Return the codeword obtained by applying the permutation wp to word.
     415    """
     416    cdef int num_chunks = wp.chunk_num
     417    cdef int i
     418    cdef codeword gate = wp.gate
     419    cdef codeword image = 0
     420    cdef codeword **images = wp.images
     421    for i from 0 <= i < num_chunks:
     422        image += images[i][(word >> i*chunk_size) & gate]
     423    return image
     424
     425def test_expand_to_ortho_basis(B=None):
     426    """
     427    This function is written in pure C for speed, and is tested from this
     428    function.
     429   
     430    INPUT:
     431    B -- a BinaryCode in standard form
     432
     433    OUTPUT:
     434    An array of codewords which represent the expansion of a basis for $B$ to a
     435    basis for $(B^\prime)^\perp$, where $B^\prime = B$ if the all-ones vector 1
     436    is in $B$, otherwise $B^\prime = \text{span}(B,1)$ (note that this guarantees
     437    that all the vectors in the span of the output have even weight).
     438   
     439    TESTS:
     440        sage: from sage.coding.binary_code import test_expand_to_ortho_basis, BinaryCode
     441        sage: M = Matrix(GF(2), [[1,1,1,1,1,1,0,0,0,0],[0,0,1,1,1,1,1,1,1,1]])
     442        sage: B = BinaryCode(M)
     443        sage: B.put_in_std_form()
     444        0
     445        sage: test_expand_to_ortho_basis(B=B)
     446        INPUT CODE:
     447        Binary [10,2] linear code, generator matrix
     448        [1010001111]
     449        [0101111111]
     450        Expanding to the basis of an orthogonal complement...
     451        Basis:
     452        0010000010
     453        0001000010
     454        0000100001
     455        0000010001
     456        0000001001
     457
     458    """
     459    cdef codeword *output
     460    cdef int k=0, i
     461    cdef BinaryCode C
     462    if not isinstance(B, BinaryCode):
     463        raise TypeError()
     464    C = B
     465    print "INPUT CODE:"
     466    print C
     467    print "Expanding to the basis of an orthogonal complement..."
     468    output = expand_to_ortho_basis(C, C.ncols)
     469    print "Basis:"
     470    while output[k]:
     471        k += 1
     472    for i from 0 <= i < k:
     473        print ''.join(reversed(Integer(output[i]).binary().zfill(C.ncols)))
     474
     475cdef codeword *expand_to_ortho_basis(BinaryCode B, int n):
     476    r"""
     477    INPUT:
     478    B -- a BinaryCode in standard form
     479    n -- the degree
     480   
     481    OUTPUT:
     482    An array of codewords which represent the expansion of a basis for $B$ to a
     483    basis for $(B^\prime)^\perp$, where $B^\prime = B$ if the all-ones vector 1
     484    is in $B$, otherwise $B^\prime = \text{span}(B,1)$ (note that this guarantees
     485    that all the vectors in the span of the output have even weight).
     486    """
     487    # assumes B is already in standard form
     488    cdef codeword *basis, word = 0, temp, new, pivots = 0, combo, parity
     489    cdef codeword n_gate = (~<codeword>0) >> ( (sizeof(codeword)<<3) - n)
     490    cdef int i, j, m, k = B.nrows, dead, d
     491    cdef WordPermutation *wp
     492    basis = <codeword *> sage_malloc( n * sizeof(codeword) )
     493    if basis is NULL:
     494        raise MemoryError()
     495    for i from 0 <= i < k:
     496        basis[i] = B.basis[i]
     497        word ^= basis[i]
     498    # If 11...1 is already a word of the code,
     499    # then being orthogonal to the code guarantees
     500    # being even weight. Otherwise, add this in.
     501    word = (~word) & n_gate
     502    if word:
     503        basis[k] = word
     504        temp = (<codeword>1) << k
     505        i = k
     506        while not word & temp:
     507            temp = temp << 1
     508            i += 1
     509        for j from 0 <= j < k:
     510            if temp & basis[j]:
     511                basis[j] ^= word
     512        temp += (<codeword>1 << k) - 1
     513        i = k
     514        word = <codeword>1 << k
     515        k += 1
     516    else: # NOTE THIS WILL NEVER HAPPEN AS CURRENTLY SET UP!
     517        temp = (<codeword>1 << k) - 1
     518        i = k
     519        word = <codeword>1 << k
     520    # Now:
     521    # k is the length of the basis so far
     522    j = k
     523    while i < n:
     524        # we are now looking at the ith free variable,
     525        # word has a 1 in the ith place, and
     526        # j is the current row we are putting in basis
     527        new = 0
     528        for m from 0 <= m < k:
     529            if basis[m] & word:
     530                new ^= basis[m]
     531        basis[j] = (new & temp) + word
     532        j += ((word ^ temp) >> i) & 1
     533        i += 1
     534        word = word << 1
     535    temp = (<codeword>1 << B.nrows) - 1
     536    for i from k <= i < n:
     537        basis[i-k] = basis[i] ^ B.words[basis[i] & temp]
     538    k = n-k
     539    i = 0
     540    word = (<codeword>1 << B.nrows)
     541    while i < k and (word & n_gate):
     542        m = i
     543        while m < k and not basis[m] & word:
     544            m += 1
     545        if m < k:
     546            pivots += word
     547            if m != i:
     548                new = basis[i]
     549                basis[i] = basis[m]
     550                basis[m] = new
     551            for j from 0 <= j < i:
     552                if basis[j] & word:
     553                    basis[j] ^= basis[i]
     554            for j from i < j < k:
     555                if basis[j] & word:
     556                    basis[j] ^= basis[i]
     557            i += 1
     558        word = word << 1
     559    for j from i <= j < n:
     560        basis[j] = 0
     561    # now basis is length i
     562    perm = range(B.nrows)
     563    perm_c = []
     564    for j from B.nrows <= j < B.ncols:
     565        if (<codeword>1 << j) & pivots:
     566            perm.append(j)
     567        else:
     568            perm_c.append(j)
     569    perm.extend(perm_c)
     570    perm.extend(range(B.ncols, n))
     571    perm_c = [0]*n
     572    for j from 0 <= j < n:
     573        perm_c[perm[j]] = j
     574    wp = create_word_perm(perm_c)
     575    for j from 0 <= j < i:
     576        basis[j] = permute_word_by_wp(wp, basis[j])
     577    for j from 0 <= j < B.nrows:
     578        B.basis[j] = permute_word_by_wp(wp, B.basis[j])
     579    dealloc_word_perm(wp)
     580    word = 0
     581    parity = 0
     582    combo = 0
     583    while 1:
     584        B.words[combo] = word
     585        parity ^= 1
     586        j = 0
     587        if not parity:
     588            while not combo & (1 << j): j += 1
     589            j += 1
     590        if j == B.nrows: break
     591        else:
     592            combo ^= (1 << j)
     593            word ^= B.basis[j]
     594    return basis
    66595
    67596cdef class BinaryCode:
    68597    """
    cdef class BinaryCode: 
    103632
    104633    """
    105634    def __new__(self, arg1, arg2=None):
    106         cdef int nrows, i, j
     635        cdef int nrows, i, j, size
    107636        cdef int nwords, other_nwords, parity, combination
    108637        cdef codeword word, glue_word
    109638        cdef BinaryCode other
    cdef class BinaryCode: 
    118647            self.nwords = 1 << nrows
    119648            nwords = self.nwords
    120649        elif isinstance(arg1, BinaryCode):
    121             other = arg1
     650            other = <BinaryCode> arg1
    122651            self.nrows = other.nrows + 1
    123652            glue_word = <codeword> arg2
    124             self.ncols = max( other.ncols , floor(log(arg2,2))+1 )
     653            size = 0
     654            while 0 < ((<codeword> 1) << size) <= glue_word:
     655                size += 1
     656            if other.ncols > size:
     657                self.ncols = other.ncols
     658            else:
     659                self.ncols = size
    125660            other_nwords = other.nwords
    126661            self.nwords = 2 * other_nwords
    127662            nrows = self.nrows
    cdef class BinaryCode: 
    133668
    134669        self.words = <codeword *> sage_malloc( nwords * sizeof(int) )
    135670        self.basis = <codeword *> sage_malloc( nrows * sizeof(int) )
    136         if not self.words or not self.basis:
    137             if self.words: sage_free(self.words)
    138             if self.basis: sage_free(self.basis)
     671        if self.words is NULL or self.basis is NULL:
     672            if self.words is not NULL: sage_free(self.words)
     673            if self.basis is not NULL: sage_free(self.basis)
    139674            raise MemoryError("Memory.")
    140675        self_words = self.words
    141676        self_basis = self.basis
    cdef class BinaryCode: 
    151686            word = <codeword> 0
    152687            parity = 0
    153688            combination = 0
    154             while True:
     689            while 1:
    155690                self_words[combination] = word
    156691                parity ^= 1
    157692                j = 0
    cdef class BinaryCode: 
    178713    def __dealloc__(self):
    179714        sage_free(self.words)
    180715        sage_free(self.basis)
     716
     717    def matrix(self):
     718        """
     719        Returns the generator matrix of the BinaryCode, i.e. the code is the
     720        rowspace of B.matrix().
     721       
     722        EXAMPLE:
     723            sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
     724            sage: from sage.coding.binary_code import *
     725            sage: B = BinaryCode(M)
     726            sage: B.matrix()
     727            [1 1 1 1 0 0]
     728            [0 0 1 1 1 1]
     729       
     730        """
     731        cdef int i, j
     732        from sage.matrix.constructor import matrix
     733        from sage.rings.all import GF
     734        rows = []
     735        for i from 0 <= i < self.nrows:
     736            row = [0]*self.ncols
     737            for j from 0 <= j < self.ncols:
     738                if self.basis[i] & ((<codeword>1) << j):
     739                    row[j] = 1
     740            rows.append(row)
     741        return matrix(GF(2), self.nrows, self.ncols, rows)
    181742
    182743    def print_data(self):
    183744        """
    cdef class BinaryCode: 
    370931        cdef int *_word_gamma
    371932        _word_gamma = <int *> sage_malloc(self.nwords * sizeof(int))
    372933        _col_gamma = <int *> sage_malloc(self.ncols * sizeof(int))
    373         if not (_col_gamma and _word_gamma):
    374             if _word_gamma: sage_free(_word_gamma)
    375             if _col_gamma: sage_free(_col_gamma)
     934        if _col_gamma is NULL or _word_gamma is NULL:
     935            if _word_gamma is not NULL: sage_free(_word_gamma)
     936            if _col_gamma is not NULL: sage_free(_col_gamma)
    376937            raise MemoryError("Memory.")
    377938        for i from 0 <= i < self.nwords:
    378939            _word_gamma[i] = word_gamma[i]
    cdef class BinaryCode: 
    384945        return result
    385946
    386947    cdef int is_automorphism(self, int *col_gamma, int *word_gamma):
    387         # TODO: optimize? check only basis?
    388948        cdef int i, j, self_nwords = self.nwords, self_ncols = self.ncols
    389         for i from 0 <= i < self_nwords:
     949        i = 1
     950        while i < self_nwords:
    390951            for j from 0 <= j < self_ncols:
    391952                if self.is_one(i, j) != self.is_one(word_gamma[i], col_gamma[j]):
    392953                    return 0
     954            i = i << 1
    393955        return 1
     956
     957    def apply_permutation(self, labeling):
     958        """
     959        Apply a column permutation to the code.
     960       
     961        INPUT:
     962        labeling -- a list permutation of the columns
     963       
     964        EXAMPLE:
     965            sage: from sage.coding.binary_code import *
     966            sage: B = BinaryCode(ExtendedBinaryGolayCode().gen_mat())
     967            sage: B
     968            Binary [24,12] linear code, generator matrix
     969            [100000000000101011100011]
     970            [010000000000111110010010]
     971            [001000000000110100101011]
     972            [000100000000110001110110]
     973            [000010000000110011011001]
     974            [000001000000011001101101]
     975            [000000100000001100110111]
     976            [000000010000101101111000]
     977            [000000001000010110111100]
     978            [000000000100001011011110]
     979            [000000000010101110001101]
     980            [000000000001010111000111]
     981            sage: B.apply_permutation(range(11,-1,-1) + range(12, 24))
     982            sage: B
     983            Binary [24,12] linear code, generator matrix
     984            [000000000001101011100011]
     985            [000000000010111110010010]
     986            [000000000100110100101011]
     987            [000000001000110001110110]
     988            [000000010000110011011001]
     989            [000000100000011001101101]
     990            [000001000000001100110111]
     991            [000010000000101101111000]
     992            [000100000000010110111100]
     993            [001000000000001011011110]
     994            [010000000000101110001101]
     995            [100000000000010111000111]
     996
     997        """
     998        # Tests for this function implicitly test _apply_permutation_to_basis
     999        # and _update_words_from_basis. These functions should not be used
     1000        # individually by the user, so they remain cdef'd.
     1001        self._apply_permutation_to_basis(labeling)
     1002        self._update_words_from_basis()
     1003
     1004    cdef void _apply_permutation_to_basis(self, object labeling):
     1005        cdef WordPermutation *wp
     1006        cdef int i
     1007        wp = create_word_perm(labeling)
     1008        for i from 0 <= i < self.nrows:
     1009            self.basis[i] = permute_word_by_wp(wp, self.basis[i])
     1010        dealloc_word_perm(wp)
     1011
     1012    cdef void _update_words_from_basis(self):
     1013        cdef codeword word
     1014        cdef int j, parity, combination
     1015        word = 0
     1016        parity = 0
     1017        combination = 0
     1018        while 1:
     1019            self.words[combination] = word
     1020            parity ^= 1
     1021            j = 0
     1022            if not parity:
     1023                while not combination & (1 << j): j += 1
     1024                j += 1
     1025            if j == self.nrows: break
     1026            else:
     1027                combination ^= (1 << j)
     1028                word ^= self.basis[j]
     1029       
     1030
     1031    cpdef int put_in_std_form(self):
     1032        """
     1033        Put the code in binary form, which is defined by an identity matrix on
     1034        the left, augmented by a matrix of data.
     1035       
     1036        EXAMPLE:
     1037            sage: from sage.coding.binary_code import *
     1038            sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
     1039            sage: B = BinaryCode(M); B
     1040            Binary [6,2] linear code, generator matrix
     1041            [111100]
     1042            [001111]
     1043            sage: B.put_in_std_form(); B
     1044            0
     1045            Binary [6,2] linear code, generator matrix
     1046            [101011]
     1047            [010111]
     1048
     1049        """
     1050        cdef codeword swap, current = 1, pivots = 0
     1051        cdef int i, j, k, row = 0
     1052        cdef object perm
     1053        while row < self.nrows:
     1054            i = row
     1055            while i < self.nrows and not self.basis[i] & current:
     1056                i += 1
     1057            if i < self.nrows:
     1058                pivots += current
     1059                if i != row:
     1060                    swap = self.basis[row]
     1061                    self.basis[row] = self.basis[i]
     1062                    self.basis[i] = swap
     1063                for j from 0 <= j < row:
     1064                    if self.basis[j] & current:
     1065                        self.basis[j] ^= self.basis[row]
     1066                for j from row < j < self.nrows:
     1067                    if self.basis[j] & current:
     1068                        self.basis[j] ^= self.basis[row]
     1069                row += 1
     1070            current = current << 1
     1071        perm = [0]*self.ncols
     1072        j = 0
     1073        k = self.nrows
     1074        for i from 0 <= i < self.ncols:
     1075            if ((<codeword>1) << i) & pivots:
     1076                perm[i] = j
     1077                j += 1
     1078            else:
     1079                perm[i] = k
     1080                k += 1
     1081        self._apply_permutation_to_basis(perm)
     1082        self._update_words_from_basis()
    3941083
    3951084cdef class OrbitPartition:
    3961085    """
    cdef class OrbitPartition: 
    4031092    * http://en.wikipedia.org/wiki/Disjoint-set_data_structure
    4041093
    4051094    """
    406     def __new__(self, nrows, ncols):
     1095    def __new__(self, int nrows, int ncols):
    4071096        cdef int col
    4081097        cdef int nwords, word
    4091098        nwords = (1 << nrows)
    cdef class OrbitPartition: 
    4171106        self.col_rank =         <int *> sage_malloc( ncols * sizeof(int) )
    4181107        self.col_min_cell_rep = <int *> sage_malloc( ncols * sizeof(int) )
    4191108        self.col_size =         <int *> sage_malloc( ncols * sizeof(int) )
    420         if not (self.wd_parent and self.wd_rank and self.wd_min_cell_rep and self.wd_size and self.col_parent and self.col_rank and self.col_min_cell_rep and self.col_size):
    421             if self.wd_parent: sage_free(self.wd_parent)
    422             if self.wd_rank: sage_free(self.wd_rank)
    423             if self.wd_min_cell_rep: sage_free(self.wd_min_cell_rep)
    424             if self.wd_size: sage_free(self.wd_size)
    425             if self.col_parent:       sage_free(self.col_parent)
    426             if self.col_rank:         sage_free(self.col_rank)
    427             if self.col_min_cell_rep: sage_free(self.col_min_cell_rep)
    428             if self.col_size:         sage_free(self.col_size)
     1109        if self.wd_parent is NULL or self.wd_rank is NULL or self.wd_min_cell_rep is NULL \
     1110        or self.wd_size is NULL or self.col_parent is NULL or self.col_rank is NULL \
     1111        or self.col_min_cell_rep is NULL or self.col_size is NULL:
     1112            if self.wd_parent is not NULL:        sage_free(self.wd_parent)
     1113            if self.wd_rank is not NULL:          sage_free(self.wd_rank)
     1114            if self.wd_min_cell_rep is not NULL:  sage_free(self.wd_min_cell_rep)
     1115            if self.wd_size is not NULL:          sage_free(self.wd_size)
     1116            if self.col_parent is not NULL:       sage_free(self.col_parent)
     1117            if self.col_rank is not NULL:         sage_free(self.col_rank)
     1118            if self.col_min_cell_rep is not NULL: sage_free(self.col_min_cell_rep)
     1119            if self.col_size is not NULL:         sage_free(self.col_size)
    4291120            raise MemoryError("Memory.")
    4301121        for word from 0 <= word < nwords:
    4311122            self.wd_parent[word] = word
    cdef class OrbitPartition: 
    6631354        cdef int *_wd_gamma
    6641355        _wd_gamma = <int *> sage_malloc(self.nwords * sizeof(int))
    6651356        _col_gamma = <int *> sage_malloc(self.ncols * sizeof(int))
    666         if not (_col_gamma and _wd_gamma):
    667             if _wd_gamma: sage_free(_wd_gamma)
    668             if _col_gamma: sage_free(_col_gamma)
     1357        if _col_gamma is NULL or _wd_gamma is NULL:
     1358            if _wd_gamma is not NULL: sage_free(_wd_gamma)
     1359            if _col_gamma is not NULL: sage_free(_col_gamma)
    6691360            raise MemoryError("Memory.")
    6701361        for i from 0 <= i < self.nwords:
    6711362            _wd_gamma[i] = wd_gamma[i]
    cdef class PartitionStack: 
    7101401        sizeof_int = sizeof(int)
    7111402
    7121403        try:
    713             self.nrows = int(arg1)
     1404            self.nrows = <int> arg1
    7141405            self.nwords = 1 << self.nrows
    715             self.ncols = int(arg2)
     1406            self.ncols = <int> arg2
    7161407        except:
    7171408            other = arg1
    7181409            self.nrows = other.nrows
    cdef class PartitionStack: 
    7361427        self.wd_counts =  <int *> sage_malloc( (self.ncols+1)  * sizeof_int )
    7371428        self.wd_output =  <int *> sage_malloc( self.nwords * sizeof_int )
    7381429
    739         if not (self.wd_ents  and self.wd_lvls    and self.col_ents   and self.col_lvls  \
    740             and self.col_degs and self.col_counts and self.col_output \
    741             and self.wd_degs  and self.wd_counts  and self.wd_output):
    742             if self.wd_ents:         sage_free(self.wd_ents)
    743             if self.wd_lvls:         sage_free(self.wd_lvls)
    744             if self.col_ents:        sage_free(self.col_ents)
    745             if self.col_lvls:        sage_free(self.col_lvls)
    746             if self.col_degs:        sage_free(self.col_degs)
    747             if self.col_counts:      sage_free(self.col_counts)
    748             if self.col_output:      sage_free(self.col_output)
    749             if self.wd_degs:         sage_free(self.wd_degs)
    750             if self.wd_counts:       sage_free(self.wd_counts)
    751             if self.wd_output:       sage_free(self.wd_output)
     1430        if self.wd_ents is NULL or self.wd_lvls is NULL or self.col_ents is NULL \
     1431        or self.col_lvls is NULL or self.col_degs is NULL or self.col_counts is NULL \
     1432        or self.col_output is NULL or self.wd_degs is NULL or self.wd_counts is NULL \
     1433        or self.wd_output is NULL:
     1434            if self.wd_ents is not NULL:    sage_free(self.wd_ents)
     1435            if self.wd_lvls is not NULL:    sage_free(self.wd_lvls)
     1436            if self.col_ents is not NULL:   sage_free(self.col_ents)
     1437            if self.col_lvls is not NULL:   sage_free(self.col_lvls)
     1438            if self.col_degs is not NULL:   sage_free(self.col_degs)
     1439            if self.col_counts is not NULL: sage_free(self.col_counts)
     1440            if self.col_output is not NULL: sage_free(self.col_output)
     1441            if self.wd_degs is not NULL:    sage_free(self.wd_degs)
     1442            if self.wd_counts is not NULL:  sage_free(self.wd_counts)
     1443            if self.wd_output is not NULL:  sage_free(self.wd_output)
    7521444            raise MemoryError("Memory.")
    7531445
    7541446        nwords = self.nwords
    cdef class PartitionStack: 
    12701962        cdef int min_is_col = 1, radix = self.radix
    12711963        cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls
    12721964        cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents
    1273         while True:
     1965        while 1:
    12741966            if self_col_lvls[i] <= k:
    12751967                if i != j and min > i - j + 1:
    12761968                    min = i - j + 1
    cdef class PartitionStack: 
    12781970                j = i + 1
    12791971            if self_col_lvls[i] == -1: break
    12801972            i += 1
    1281         i = 0; j = 0
    1282         while True:
    1283             if self_wd_lvls[i] <= k:
    1284                 if i != j and min > i - j + 1:
    1285                     min = i - j + 1
    1286                     min_is_col = 0
    1287                     location = j
    1288                 j = i + 1
    1289             if self_wd_lvls[i] == -1: break
    1290             i += 1
     1973#        i = 0; j = 0
     1974#        while 1:
     1975#            if self_wd_lvls[i] <= k:
     1976#                if i != j and min > i - j + 1:
     1977#                    min = i - j + 1
     1978#                    min_is_col = 0
     1979#                    location = j
     1980#                j = i + 1
     1981#            if self_wd_lvls[i] == -1: break
     1982#            i += 1
    12911983        # location now points to the beginning of the first, smallest,
    12921984        # nontrivial cell
    12931985        j = location
    cdef class PartitionStack: 
    12981990        for i from 0 <= i < ell:
    12991991            W[start+i] = 0
    13001992        if min_is_col:
    1301             while True:
     1993            while 1:
    13021994                if self_col_lvls[j] <= k: break
    13031995                j += 1
    13041996            # j now points to the last element of the cell
    cdef class PartitionStack: 
    13082000                i += 1
    13092001            return self_col_ents[location]
    13102002        else:
    1311             while True:
     2003            while 1:
    13122004                if self_wd_lvls[j] <= k: break
    13132005                j += 1
    13142006            # j now points to the last element of the cell
    cdef class PartitionStack: 
    15702262    cdef int col_degree(self, BinaryCode CG, int col, int wd_ptr, int k):
    15712263        cdef int i = 0
    15722264        cdef int *self_wd_lvls = self.wd_lvls, *self_wd_ents = self.wd_ents
    1573         while True:
     2265        while 1:
    15742266            if CG.is_one(self_wd_ents[wd_ptr], col): i += 1
    15752267            if self_wd_lvls[wd_ptr] > k: wd_ptr += 1
    15762268            else: break
    cdef class PartitionStack: 
    18072499        cdef int i, alpha_length = len(alpha)
    18082500        cdef int *_alpha = <int *> sage_malloc( (self.nwords + self.ncols) * sizeof(int) )
    18092501        cdef int *ham_wts = hamming_weights()
    1810         if not _alpha:
    1811             sage_free(_alpha)
     2502        if _alpha is NULL:
    18122503            raise MemoryError("Memory.")
    18132504        for i from 0 <= i < alpha_length:
    18142505            if alpha[i][0]:
    cdef class PartitionStack: 
    18382529#                    print self
    18392530                    i = j; s = 0
    18402531                    invariant += 8
    1841                     while True:
     2532                    while 1:
    18422533#                        print 'col_i', self_col_ents[i]
    18432534#                        print 'alpha[m]^flag', alpha[m]^flag
    18442535                        self_col_degs[i-j] = self.col_degree(CG, self_col_ents[i], alpha[m]^flag, k)
    cdef class PartitionStack: 
    18582549                                break
    18592550                            q += 1
    18602551                        r = j
    1861                         while True:
     2552                        while 1:
    18622553                            if r == j or self.col_lvls[r-1] == k:
    18632554                                if r != t:
    18642555                                    alpha[alpha_length] = r
    cdef class PartitionStack: 
    18752566#                    print self
    18762567                    i = j; s = 0
    18772568                    invariant += 64
    1878                     while True:
     2569                    while 1:
    18792570#                        print 'i', i
    18802571                        self_wd_degs[i-j] = self.wd_degree(CG, self_wd_ents[i], alpha[m], k, ham_wts)
    18812572#                        print 'deg', self_wd_degs[i-j]
    cdef class PartitionStack: 
    18952586                            q += 1
    18962587                        j ^= flag
    18972588                        r = j
    1898                         while True:
     2589                        while 1:
    18992590                            if r == j or self.wd_lvls[r-1] == k:
    19002591                                if r != t_w:
    19012592                                    alpha[alpha_length] = r^flag
    cdef class PartitionStack: 
    19062597                        invariant += (i-j)
    19072598                    j = i
    19082599            m += 1
    1909         if invariant != -1:
    1910             return invariant
    1911         else:
    1912             return 0
     2600        return invariant
    19132601
    19142602    def _clear(self, k):
    19152603        """
    cdef class PartitionStack: 
    19972685        cdef int *self_wd_ents = self.wd_ents
    19982686        cdef codeword *CG_words = CG.words
    19992687        cdef int i, j, l, m, span = 1, ncols = self.ncols, nwords = self.nwords
    2000         for i from 0 <= i < nwords: # TODO: probably don't need to check i == 0 here!
     2688        for i from 0 < i < nwords:
    20012689            for j from 0 <= j < ncols:
    20022690                l = CG.is_one(self.wd_ents[i], self.col_ents[j])
    20032691                m = CG.is_one(other.wd_ents[i], other.col_ents[j])
    cdef class PartitionStack: 
    20672755        self.find_basis(ham_wts)
    20682756        sage_free(ham_wts)
    20692757
    2070     cdef void find_basis(self, int *ham_wts):
     2758    cdef int find_basis(self, int *ham_wts):
    20712759        cdef int i = 0, j, k, nwords = self.nwords, weight, basis_elts = 0, nrows = self.nrows
    20722760        cdef int *self_wd_ents = self.wd_ents
    2073         if not self.basis_locations:
    2074             self.basis_locations = <int *> sage_malloc( nrows * sizeof(int) )
    2075         if not self.basis_locations:
     2761        if self.basis_locations is NULL:
     2762            self.basis_locations = <int *> sage_malloc( 2 * nrows * sizeof(int) )
     2763        if self.basis_locations is NULL:
    20762764            raise MemoryError("Memory.")
    20772765        while i < nwords:
    20782766            j = self_wd_ents[i]
    cdef class PartitionStack: 
    20852773                self.basis_locations[k] = i
    20862774                if basis_elts == nrows: break
    20872775            i += 1
     2776        for i from 0 <= i < nrows:
     2777            self.basis_locations[nrows + i] = self_wd_ents[1 << i]
    20882778
    20892779    def _get_permutation(self, other):
    20902780        """
    cdef class PartitionStack: 
    21272817
    21282818        """
    21292819        cdef int i
    2130         cdef int *ham_wts = hamming_weights()
    21312820        cdef int *word_g = <int *> sage_malloc( self.nwords * sizeof(int) )
    21322821        cdef int *col_g = <int *> sage_malloc( self.ncols * sizeof(int) )
    2133         if not (word_g and col_g):
    2134             if word_g: sage_free(word_g)
    2135             if col_g: sage_free(col_g)
    2136             raise MemoryError("Memory.")
    2137         self.get_permutation(other, word_g, col_g, ham_wts)
    2138         sage_free(ham_wts)
     2822        if word_g is NULL or col_g is NULL:
     2823            if word_g is not NULL: sage_free(word_g)
     2824            if col_g is not NULL: sage_free(col_g)
     2825            raise MemoryError("Memory.")
     2826        self.get_permutation(other, word_g, col_g)
    21392827        word_l = [word_g[i] for i from 0 <= i < self.nwords]
    21402828        col_l = [col_g[i] for i from 0 <= i < self.ncols]
    21412829        sage_free(word_g)
    21422830        sage_free(col_g)
    21432831        return word_l, col_l
    21442832
    2145     cdef void get_permutation(self, PartitionStack other, int *word_gamma, int *col_gamma, int *ham_wts):
     2833    cdef void get_permutation(self, PartitionStack other, int *word_gamma, int *col_gamma):
    21462834        cdef int i
    21472835        cdef int *self_wd_ents = self.wd_ents, *other_wd_ents = other.wd_ents
    21482836        cdef int *self_col_ents = self.col_ents, *other_col_ents = other.col_ents
    cdef class BinaryCodeClassifier: 
    21552843cdef class BinaryCodeClassifier:
    21562844
    21572845    def __new__(self):
    2158         self.radix = sizeof(int) << 3
     2846        self.radix = sizeof(codeword) << 3
    21592847        self.ham_wts = hamming_weights()
    21602848        self.L = 100 # memory limit for Phi and Omega- multiply by 8KB
    21612849        self.aut_gens_size = self.radix * 100
    cdef class BinaryCodeClassifier: 
    21702858        self.Omega =   <unsigned int *> sage_malloc( self.Phi_size * self.L         * sizeof(unsigned int) )
    21712859        self.W =       <unsigned int *> sage_malloc( self.Phi_size * self.radix * 2 * sizeof(unsigned int) )
    21722860
    2173         self.aut_gp_gens = <int *> sage_malloc( self.aut_gens_size             * sizeof(int) )
    2174         self.c_gamma =     <int *> sage_malloc( self.radix                     * sizeof(int) )
    2175         self.labeling =    <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2176         self.Lambda1 =     <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2177         self.Lambda2 =     <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2178         self.Lambda3 =     <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2179         self.v =           <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2180         self.e =           <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2181 
    2182         if not (self.Phi and self.Omega and self.W and self.Lambda1 and self.Lambda2 and self.Lambda3 \
    2183             and self.w_gamma and self.c_gamma and self.alpha and self.v and self.e and self.aut_gp_gens \
    2184             and self.labeling):
    2185             if self.Phi:          sage_free(self.Phi)
    2186             if self.Omega:        sage_free(self.Omega)
    2187             if self.W:            sage_free(self.W)
    2188             if self.Lambda1:      sage_free(self.Lambda1)
    2189             if self.Lambda2:      sage_free(self.Lambda2)
    2190             if self.Lambda3:      sage_free(self.Lambda3)
    2191             if self.w_gamma:      sage_free(self.w_gamma)
    2192             if self.c_gamma:      sage_free(self.c_gamma)
    2193             if self.alpha:        sage_free(self.alpha)
    2194             if self.v:            sage_free(self.v)
    2195             if self.e:            sage_free(self.e)
    2196             if self.aut_gp_gens:  sage_free(self.aut_gp_gens)
    2197             if self.labeling:     sage_free(self.labeling)
     2861        self.base =        <int *> sage_malloc( self.radix          * sizeof(int) )
     2862        self.aut_gp_gens = <int *> sage_malloc( self.aut_gens_size  * sizeof(int) )
     2863        self.c_gamma =     <int *> sage_malloc( self.radix          * sizeof(int) )
     2864        self.labeling =    <int *> sage_malloc( self.radix * 3      * sizeof(int) )
     2865        self.Lambda1 =     <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2866        self.Lambda2 =     <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2867        self.Lambda3 =     <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2868        self.v =           <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2869        self.e =           <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2870
     2871        if self.Phi is NULL or self.Omega is NULL or self.W is NULL or self.Lambda1 is NULL \
     2872        or self.Lambda2 is NULL or self.Lambda3 is NULL or self.w_gamma is NULL \
     2873        or self.c_gamma is NULL or self.alpha is NULL or self.v is NULL or self.e is NULL \
     2874        or self.aut_gp_gens is NULL or self.labeling is NULL or self.base is NULL:
     2875            if self.Phi is not NULL:          sage_free(self.Phi)
     2876            if self.Omega is not NULL:        sage_free(self.Omega)
     2877            if self.W is not NULL:            sage_free(self.W)
     2878            if self.Lambda1 is not NULL:      sage_free(self.Lambda1)
     2879            if self.Lambda2 is not NULL:      sage_free(self.Lambda2)
     2880            if self.Lambda3 is not NULL:      sage_free(self.Lambda3)
     2881            if self.w_gamma is not NULL:      sage_free(self.w_gamma)
     2882            if self.c_gamma is not NULL:      sage_free(self.c_gamma)
     2883            if self.alpha is not NULL:        sage_free(self.alpha)
     2884            if self.v is not NULL:            sage_free(self.v)
     2885            if self.e is not NULL:            sage_free(self.e)
     2886            if self.aut_gp_gens is not NULL:  sage_free(self.aut_gp_gens)
     2887            if self.labeling is not NULL:     sage_free(self.labeling)
     2888            if self.base is not NULL:         sage_free(self.base)
    21982889            raise MemoryError("Memory.")
    21992890
    22002891    def __dealloc__(self):
    2201         if self.ham_wts: sage_free(self.ham_wts)
    2202         if self.Phi:     sage_free(self.Phi)
    2203         if self.Omega:   sage_free(self.Omega)
    2204         if self.W:       sage_free(self.W)
    2205         if self.Lambda1: sage_free(self.Lambda1)
    2206         if self.Lambda2: sage_free(self.Lambda2)
    2207         if self.Lambda3: sage_free(self.Lambda3)
    2208         if self.c_gamma: sage_free(self.c_gamma)
    2209 
    2210         if self.w_gamma:   sage_free(self.w_gamma)
    2211         if self.alpha:     sage_free(self.alpha)
    2212 
    2213         if self.v:           sage_free(self.v)
    2214         if self.e:           sage_free(self.e)
    2215         if self.aut_gp_gens: sage_free(self.aut_gp_gens)
    2216         if self.labeling:    sage_free(self.labeling)
     2892        sage_free(self.ham_wts)
     2893        sage_free(self.Phi)
     2894        sage_free(self.Omega)
     2895        sage_free(self.W)
     2896        sage_free(self.Lambda1)
     2897        sage_free(self.Lambda2)
     2898        sage_free(self.Lambda3)
     2899        sage_free(self.c_gamma)
     2900        sage_free(self.w_gamma)
     2901        sage_free(self.alpha)
     2902        sage_free(self.v)
     2903        sage_free(self.e)
     2904        sage_free(self.aut_gp_gens)
     2905        sage_free(self.labeling)
     2906        sage_free(self.base)
    22172907
    22182908    cdef void record_automorphism(self, int *gamma, int ncols):
    22192909        cdef int i, j
    22202910        if self.aut_gp_index + ncols > self.aut_gens_size:
    22212911            self.aut_gens_size *= 2
    22222912            self.aut_gp_gens = <int *> sage_realloc( self.aut_gp_gens, self.aut_gens_size * sizeof(int) )
    2223             if not self.aut_gp_gens:
     2913            if self.aut_gp_gens is NULL:
    22242914                raise MemoryError("Memory.")
    22252915        j = self.aut_gp_index
    22262916        for i from 0 <= i < ncols:
    cdef class BinaryCodeClassifier: 
    22362926            verbosity - a nonnegative integer
    22372927       
    22382928        OUTPUT:
    2239             a tuple, (gens, labeling, size)
     2929            a tuple, (gens, labeling, size, base)
    22402930            gens -- a list of permutations (in list form) representing generators
    22412931                of the permutation automorphism group of the code CC.
    22422932            labeling -- a permutation representing the canonical labeling of the
    2243                 code. mostly for internal use; if the dimension of the code is k
    2244                 and the degree (number of columns) is n, then the first n entries
    2245                 describe the relabeling on the columns, and the next k describe
    2246                 where the basis is sent.
     2933                code. mostly for internal use; entries describe the relabeling
     2934                on the columns.
    22472935            size -- the order of the automorphism group.
     2936            base -- a set of cols whose action determines the action on all cols
    22482937       
    22492938        EXAMPLES:
    22502939            sage: import sage.coding.binary_code
    cdef class BinaryCodeClassifier: 
    22582947            ... [0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],\
    22592948            ... [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]])
    22602949            sage: B = BinaryCode(M)
    2261             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2950            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    22622951            sage: S = SymmetricGroup(M.ncols())
    22632952            sage: L = [S([x+1 for x in g]) for g in gens]
    22642953            sage: PermutationGroup(L).order()
    cdef class BinaryCodeClassifier: 
    22722961            ... [0,0,0,0,0,1,0,1,0,0,0,1,1,1,1,1,1],\
    22732962            ... [0,0,0,1,1,0,0,0,0,1,1,0,1,1,0,1,1]])
    22742963            sage: B = BinaryCode(M)
    2275             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2964            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    22762965            sage: S = SymmetricGroup(M.ncols())
    22772966            sage: L = [S([x+1 for x in g]) for g in gens]
    22782967            sage: PermutationGroup(L).order()
    cdef class BinaryCodeClassifier: 
    22902979            ... [0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0],\
    22912980            ... [0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1]])
    22922981            sage: B = BinaryCode(M)
    2293             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2982            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    22942983            sage: S = SymmetricGroup(M.ncols())
    22952984            sage: L = [S([x+1 for x in g]) for g in gens]
    22962985            sage: PermutationGroup(L).order()
    cdef class BinaryCodeClassifier: 
    23113000            ... [1,0,0,1,0,1,1,1,0,0,0,0,1,0,0,1,0,0,0,1,1,1],\
    23123001            ... [0,0,1,0,1,1,1,0,0,0,1,1,0,0,1,0,0,0,1,1,1,0]])
    23133002            sage: B = BinaryCode(M)
    2314             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     3003            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    23153004            sage: S = SymmetricGroup(M.ncols())
    23163005            sage: L = [S([x+1 for x in g]) for g in gens]
    23173006            sage: PermutationGroup(L).order()
    cdef class BinaryCodeClassifier: 
    23213010
    23223011            sage: B = BinaryCode(Matrix(GF(2),[[1,0,1],[0,1,1]]))
    23233012            sage: BC._aut_gp_and_can_label(B)
    2324             ([[0, 2, 1], [1, 0, 2]], [0, 1, 2, 2, 1], 6)
     3013            ([[0, 2, 1], [1, 0, 2]], [0, 1, 2], 6, [0, 1])
    23253014
    23263015            sage: B = BinaryCode(Matrix(GF(2),[[1,1,1,1]]))
    23273016            sage: BC._aut_gp_and_can_label(B)
    2328             ([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], [0, 1, 2, 3, 1], 24)
     3017            ([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], [0, 1, 2, 3], 24, [0, 1, 2])
    23293018
    23303019            sage: B = BinaryCode(Matrix(GF(2),[[0,0,0,0,0,0,0,0,0,0,0,0,0,0,1]]))
    2331             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     3020            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    23323021            sage: size
    23333022            87178291200
    23343023
    cdef class BinaryCodeClassifier: 
    23833072            gen = [self.aut_gp_gens[i+j] for j from 0 <= j < C.ncols]
    23843073            py_aut_gp_gens.append(gen)
    23853074            i += C.ncols
    2386         py_labeling = [self.labeling[i] for i from 0 <= i < C.ncols] \
    2387                     + [self.labeling[i+C.ncols] for i from 0 <= i < C.nrows]
     3075        py_labeling = [self.labeling[i] for i from 0 <= i < C.ncols]
     3076        base = []
     3077        for i from 0 <= i < self.radix:
     3078            if self.base[i] == -1:
     3079                break
     3080            base.append(self.base[i])
    23883081        aut_gp_size = self.aut_gp_size
    2389         return py_aut_gp_gens, py_labeling, aut_gp_size
     3082        return py_aut_gp_gens, py_labeling, aut_gp_size, base
    23903083
    23913084    cdef void aut_gp_and_can_label(self, BinaryCode C, int verbosity):
    23923085
    cdef class BinaryCodeClassifier: 
    24513144            self.Phi =     <unsigned int *> sage_realloc(self.Phi,   self.Phi_size * self.L         * sizeof(int) )
    24523145            self.Omega =   <unsigned int *> sage_realloc(self.Omega, self.Phi_size * self.L         * sizeof(int) )
    24533146            self.W =       <unsigned int *> sage_realloc(self.W,     self.Phi_size * self.radix * 2 * sizeof(int) )
    2454             if not (self.w_gamma and self.alpha and self.Phi and self.Omega and self.W):
    2455                 if self.w_gamma: sage_free(self.w_gamma)
    2456                 if self.alpha: sage_free(self.alpha)
    2457                 if self.Phi: sage_free(self.Phi)
    2458                 if self.Omega: sage_free(self.Omega)
    2459                 if self.W: sage_free(self.W)
     3147            if self.w_gamma is NULL or self.alpha is NULL or self.Phi is NULL or self.Omega is NULL or self.W is NULL:
     3148                if self.w_gamma is not NULL: sage_free(self.w_gamma)
     3149                if self.alpha is not NULL: sage_free(self.alpha)
     3150                if self.Phi is not NULL: sage_free(self.Phi)
     3151                if self.Omega is not NULL: sage_free(self.Omega)
     3152                if self.W is not NULL: sage_free(self.W)
    24603153                raise MemoryError("Memory.")
    24613154        word_gamma = self.w_gamma
    24623155        alpha = self.alpha # think of alpha as of length exactly nwords + ncols
    cdef class BinaryCodeClassifier: 
    25523245            elif state == 2: # Move down the search tree one level by refining nu:
    25533246                             # split out a vertex, and refine nu against it
    25543247                k += 1
    2555                 # TODO: Consider removing the following lines
    2556                 if k >= 2*self.radix: raise RuntimeError(\
    2557                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    25583248                nu.clear(k)
    25593249
    25603250                alpha[0] = nu.split_vertex(v[k-1], k)
    cdef class BinaryCodeClassifier: 
    26303320                        k = hzb__h_rho
    26313321                    else:
    26323322                        k = hh-1
    2633                 # TODO: Consider removing the following lines
    2634                 if k >= 2*self.radix: raise RuntimeError(\
    2635                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    26363323                # TODO: is the following line necessary?
    26373324                if k == -1: k = 0
    26383325               
    cdef class BinaryCodeClassifier: 
    26603347                # equivalent to zeta, i.e. there is no automorphism to discover.
    26613348                if k < hzf__h_zeta: state = 8; continue
    26623349
    2663                 nu.get_permutation(zeta, word_gamma, col_gamma, ham_wts)
     3350                nu.get_permutation(zeta, word_gamma, col_gamma)
    26643351#                print "gamma:", str([[word_gamma[i] for i from 0 <= i < nwords], [col_gamma[i] for i from 0 <= i < ncols]]).replace(' ','')
    26653352#                print Theta
    26663353                # if C^gamma == C, the permutation is an automorphism, goto 10
    cdef class BinaryCodeClassifier: 
    26863373                if j < 0: state = 6; continue
    26873374
    26883375                # if C(nu) == C(rho), get the automorphism and goto 10
    2689                 rho.get_permutation(nu, word_gamma, col_gamma, ham_wts)
     3376                rho.get_permutation(nu, word_gamma, col_gamma)
    26903377#                print "gamma:", str([[word_gamma[i] for i from 0 <= i < nwords], [col_gamma[i] for i from 0 <= i < ncols]]).replace(' ','')
    26913378#                print Theta
    26923379                state = 10
    cdef class BinaryCodeClassifier: 
    27683455
    27693456                # Otherwise, proceed to where zeta meets nu:
    27703457                k = h
    2771                 # TODO: Consider removing the following lines
    2772                 if k >= 2*self.radix: raise RuntimeError(\
    2773                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    27743458                state = 13
    27753459
    27763460            elif state == 11: # We have just found a new automorphism, and deduced that there may
    27773461                # be a better canonical label below the current branch off of zeta. So go to where
    27783462                # nu meets rho
    27793463                k = hb
    2780                 # TODO: Consider removing the following lines
    2781                 if k >= 2*self.radix: raise RuntimeError(\
    2782                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    27833464                state = 12
    27843465
    27853466            elif state == 12: # Coming here from either state 6 or 11, the algorithm has discovered
    cdef class BinaryCodeClassifier: 
    29703651                ht = k # nodes descended from zeta[ht] are all equivalent
    29713652                hzf__h_zeta = k # max such that indicators for zeta and nu agree
    29723653                zeta = PartitionStack(nu)
    2973                 zeta.find_basis(ham_wts)
     3654                for i from 0 <= i < k:
     3655                    self.base[i] = v[i]
     3656                self.base_size = k
     3657                if k != self.radix:
     3658                    self.base[k] = -1
    29743659                # (POINT B)
    29753660                k -= 1
    29763661                rho = PartitionStack(nu)
    cdef class BinaryCodeClassifier: 
    29853670        rho.find_basis(ham_wts)
    29863671        for i from 0 <= i < ncols:
    29873672            self.labeling[rho.col_ents[i]] = i
    2988         for i from 0 <= i < nrows:
     3673        for i from 0 <= i < 2*nrows:
    29893674            self.labeling[i+ncols] = rho.basis_locations[i]
    29903675
    2991 
    2992 
    2993 
     3676    def put_in_canonical_form(self, BinaryCode B):
     3677        """
     3678        Puts the code into canonical form.
     3679       
     3680        Canonical form is obtained by performing row reduction, permuting the
     3681        pivots to the front so that the generator matrix is of the form: the
     3682        identity matrix augmented to the right by arbitrary data.
     3683       
     3684        EXAMPLE:
     3685            sage: from sage.coding.binary_code import *
     3686            sage: BC = BinaryCodeClassifier()
     3687            sage: B = BinaryCode(ExtendedBinaryGolayCode().gen_mat())
     3688            sage: B.apply_permutation(range(24,-1,-1))
     3689            sage: B
     3690            Binary [24,12] linear code, generator matrix
     3691            [011000111010100000000000]
     3692            [001001001111100000000001]
     3693            [011010100101100000000010]
     3694            [001101110001100000000100]
     3695            [010011011001100000001000]
     3696            [010110110011000000010000]
     3697            [011101100110000000100000]
     3698            [000011110110100001000000]
     3699            [000111101101000010000000]
     3700            [001111011010000100000000]
     3701            [010110001110101000000000]
     3702            [011100011101010000000000]
     3703            sage: BC.put_in_canonical_form(B)
     3704            sage: B
     3705            Binary [24,12] linear code, generator matrix
     3706            [100000000000001100111001]
     3707            [010000000000001010001111]
     3708            [001000000000001111010010]
     3709            [000100000000010110101010]
     3710            [000010000000010110010101]
     3711            [000001000000010001101101]
     3712            [000000100000011000110110]
     3713            [000000010000011111001001]
     3714            [000000001000010101110011]
     3715            [000000000100010011011110]
     3716            [000000000010001011110101]
     3717            [000000000001001101101110]
     3718
     3719        """
     3720        aut_gp_gens, labeling, size, base = self._aut_gp_and_can_label(B)
     3721        B._apply_permutation_to_basis(labeling)
     3722        B.put_in_std_form()
     3723
     3724    def generate_children(self, BinaryCode B, int n, int d=2):
     3725        """
     3726        Use canonical augmentation to generate children of the code B.
     3727       
     3728        INPUT:
     3729        B -- a BinaryCode
     3730        n -- limit on the degree of the code
     3731        d -- test whether new vector has weight divisible by d. If d==4, this
     3732            ensures that all doubly-even canonically augmented children are
     3733            generated.
     3734       
     3735        EXAMPLE:
     3736            sage: from sage.coding.binary_code import *
     3737            sage: BC = BinaryCodeClassifier()
     3738            sage: B = BinaryCode(Matrix(GF(2), [[1,1,1,1]]))
     3739            sage: BC.generate_children(B, 6, 4)
     3740            [[1 1 1 1 0 0]
     3741            [0 1 0 1 1 1]]
     3742       
     3743        NOTE:
     3744        The function self_orthogonal_binary_codes makes heavy use of this
     3745        function.
     3746       
     3747        MORE EXAMPLES:
     3748            sage: soc_iter = self_orthogonal_binary_codes(12, 6, 4)
     3749            sage: L = list(soc_iter)
     3750            sage: for n in range(0, 13):
     3751            ...     s = 'n=%2d : '%n
     3752            ...     for k in range(1,7):
     3753            ...         s += '%3d '%len([C for C in L if C.length() == n and C.dimension() == k])
     3754            ...     print s
     3755            n= 0 :   0   0   0   0   0   0
     3756            n= 1 :   0   0   0   0   0   0
     3757            n= 2 :   0   0   0   0   0   0
     3758            n= 3 :   0   0   0   0   0   0
     3759            n= 4 :   1   0   0   0   0   0
     3760            n= 5 :   0   0   0   0   0   0
     3761            n= 6 :   0   1   0   0   0   0
     3762            n= 7 :   0   0   1   0   0   0
     3763            n= 8 :   1   1   1   1   0   0
     3764            n= 9 :   0   0   0   0   0   0
     3765            n=10 :   0   1   1   1   0   0
     3766            n=11 :   0   0   1   1   0   0
     3767            n=12 :   1   2   3   4   2   0
     3768       
     3769        """
     3770        cdef BinaryCode m
     3771        cdef codeword *ortho_basis, *B_can_lab, *m_can_lab, current, swap
     3772        cdef codeword word, temp, gate, nonzero_gate, orbit, bwd, k_gate
     3773        cdef codeword *temp_basis, swap, *orbit_checks, orb_chx_size, orb_chx_shift, radix_gate
     3774        cdef WordPermutation *gwp, *hwp, *can_lab, *can_lab_inv
     3775        cdef WordPermutation **parent_generators
     3776        cdef BinaryCode B_aug
     3777        cdef int i, ii, j, jj, ij, k = 0, parity, combo, num_gens
     3778        cdef int base_size, *multimod2_index, row
     3779        cdef int *ham_wts = self.ham_wts
     3780        cdef int *num_inner_gens, *num_outer_gens, *v, log_2_radix
     3781        cdef bint bingo, bingo2, bingo3
     3782
     3783        B.put_in_std_form()
     3784        ortho_basis = expand_to_ortho_basis(B, n) # modifies B!
     3785
     3786#        print 'parent:'
     3787#        print B
     3788        aut_gp_gens, labeling, size, base = self._aut_gp_and_can_label(B)
     3789        B_can_lab = <codeword *> sage_malloc(B.nrows * sizeof(codeword))
     3790        can_lab = create_word_perm(labeling[:B.ncols])
     3791        if B_can_lab is NULL or m_can_lab is NULL or can_lab is NULL:
     3792            sage_free(ortho_basis)
     3793            if B_can_lab is not NULL:
     3794                sage_free(B_can_lab)
     3795            if can_lab is not NULL:
     3796                sage_free(can_lab)
     3797            raise MemoryError()
     3798        for i from 0 <= i < B.nrows:
     3799            B_can_lab[i] = permute_word_by_wp(can_lab, B.basis[i])
     3800        dealloc_word_perm(can_lab)
     3801        row = 0
     3802        current = 1
     3803        while row < B.nrows:
     3804            i = row
     3805            while i < B.nrows and not B_can_lab[i] & current:
     3806                i += 1
     3807            if i < B.nrows:
     3808                if i != row:
     3809                    swap = B_can_lab[row]
     3810                    B_can_lab[row] = B_can_lab[i]
     3811                    B_can_lab[i] = swap
     3812                for j from 0 <= j < row:
     3813                    if B_can_lab[j] & current:
     3814                        B_can_lab[j] ^= B_can_lab[row]
     3815                for j from row < j < B.nrows:
     3816                    if B_can_lab[j] & current:
     3817                        B_can_lab[j] ^= B_can_lab[row]
     3818                row += 1
     3819            current = current << 1
     3820        num_gens = len(aut_gp_gens)
     3821        base_size = len(base)
     3822
     3823#        print 'gens:'
     3824#        for g in aut_gp_gens:
     3825#            print g
     3826
     3827        parent_generators = <WordPermutation **> sage_malloc( len(aut_gp_gens) * sizeof(WordPermutation*) )
     3828        temp_basis = <codeword *> sage_malloc( self.radix * sizeof(codeword) )
     3829
     3830        output = []
     3831
     3832       
     3833        for i from 0 <= i < len(aut_gp_gens):
     3834            parent_generators[i] = create_word_perm(aut_gp_gens[i] + range(B.ncols, n))
     3835
     3836        word = 0
     3837        while ortho_basis[k] & (((<codeword>1) << B.ncols) - 1):
     3838            k += 1
     3839        j = k
     3840        while ortho_basis[j]:
     3841            word ^= ortho_basis[j]
     3842            j += 1
     3843       
     3844#        print "ortho_basis:"
     3845#        for i from 0 <= i < k:
     3846#            print ''.join(reversed(Integer(ortho_basis[i]).binary().zfill(n)))
     3847#        print '-'
     3848#        for i from k <= i < j:
     3849#            print ''.join(reversed(Integer(ortho_basis[i]).binary().zfill(n)))
     3850#        print 'word:'
     3851#        print ''.join(reversed(Integer(word).binary().zfill(n)))
     3852       
     3853        log_2_radix = 0
     3854        while ((<codeword>1) << log_2_radix) < self.radix:
     3855            log_2_radix += 1
     3856        # now we assume (<codeword>1 << log_2_radix) == self.radix
     3857        if k < log_2_radix:
     3858            orb_chx_size = 0
     3859        else:
     3860            orb_chx_size = k - log_2_radix
     3861        orbit_checks = <codeword *> sage_malloc( ((<codeword>1) << orb_chx_size) * sizeof(codeword) )
     3862        if orbit_checks is NULL:
     3863            raise MemoryError()
     3864        for temp from 0 <= temp < ((<codeword>1) << orb_chx_size):
     3865            orbit_checks[temp] = 0
     3866       
     3867       
     3868        combo = 0
     3869        parity = 0
     3870        gate = (<codeword>1 << B.nrows) - 1
     3871        k_gate = (<codeword>1 << k) - 1
     3872        nonzero_gate = ( (<codeword>1 << (n-B.ncols)) - 1 ) << B.ncols
     3873        radix_gate = (((<codeword>1) << log_2_radix) - 1)
     3874#        print 'gate:', ''.join(reversed(Integer(gate).binary().zfill(n)))
     3875#        print 'gate:', ''.join(reversed(Integer(nonzero_gate).binary().zfill(n)))
     3876        while 1:
     3877#            print '    while 1'
     3878#            print '    ' + ''.join(reversed(Integer(word).binary().zfill(n)))
     3879            if nonzero_gate & word == nonzero_gate and \
     3880              (ham_wts[word & 65535] + ham_wts[(word >> 16) & 65535])%d == 0:
     3881#                print ''.join(reversed(Integer(word).binary().zfill(n)))
     3882                temp = (word >> B.nrows) & ((<codeword>1 << k) - 1)
     3883#                print "if not orbit_checks[temp >> log_2_radix] & ((<codeword>1) << (temp & radix_gate)):"
     3884#                print temp >> log_2_radix
     3885#                print temp & radix_gate
     3886                if not orbit_checks[temp >> log_2_radix] & ((<codeword>1) << (temp & radix_gate)):
     3887                    B_aug = BinaryCode(B, word)
     3888#                    print 'child:'
     3889#                    print B_aug
     3890#                    print 'canonically labeling child'
     3891                    aug_aut_gp_gens, aug_labeling, aug_size, aug_base = self._aut_gp_and_can_label(B_aug)
     3892#                    print 'done canonically labeling child'
     3893                    # check if (B, B_aug) ~ (m(B_aug), B_aug)
     3894
     3895                    can_lab = create_word_perm(aug_labeling[:n])
     3896#                    print 'relabeling:'
     3897#                    print [self.labeling[j] for j from 0 <= j < n]
     3898                    can_lab_inv = create_inv_word_perm(can_lab)
     3899                    for j from 0 <= j < B_aug.nrows:
     3900                        temp_basis[j] = permute_word_by_wp(can_lab, B_aug.basis[j])
     3901#                    print 'temp_basis:'
     3902#                    for j from 0 <= j < B_aug.nrows:
     3903#                        print ''.join(reversed(Integer(temp_basis[j]).binary().zfill(n)))
     3904
     3905                    # row reduce to get canonical label
     3906                    i = 0
     3907                    j = 0
     3908                    while j < B_aug.nrows:
     3909                        ii = j
     3910                        while ii < B_aug.nrows and not temp_basis[ii] & (<codeword>1 << i):
     3911                            ii += 1
     3912                        if ii != B_aug.nrows:
     3913                            if ii != j:
     3914                                swap = temp_basis[ii]
     3915                                temp_basis[ii] = temp_basis[j]
     3916                                temp_basis[j] = swap
     3917                            for jj from 0 <= jj < j:
     3918                                if temp_basis[jj] & (<codeword>1 << i):
     3919                                    temp_basis[jj] ^= temp_basis[j]
     3920                            for jj from j < jj < B_aug.nrows:
     3921                                if temp_basis[jj] & (<codeword>1 << i):
     3922                                    temp_basis[jj] ^= temp_basis[j]
     3923                            j += 1
     3924                        i += 1
     3925                    # done row reduction
     3926                           
     3927#                    print 'temp_basis:'
     3928                    for j from 0 <= j < B.nrows:
     3929                        temp_basis[j] = permute_word_by_wp(can_lab_inv, temp_basis[j])
     3930#                        print ''.join(reversed(Integer(temp_basis[j]).binary().zfill(n)))
     3931                    from sage.matrix.constructor import matrix
     3932                    from sage.rings.all import ZZ
     3933                    from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroupElement
     3934                    from sage.interfaces.gap import gap
     3935                    rs = []
     3936                    for i from 0 <= i < B.nrows:
     3937                        r = []
     3938                        for j from 0 <= j < n:
     3939                            r.append((((<codeword>1)<<j)&temp_basis[i])>>j)
     3940                        rs.append(r)
     3941                    m = BinaryCode(matrix(ZZ, rs))
     3942#                    print 'm:'
     3943#                    print m
     3944                    m_aut_gp_gens, m_labeling, m_size, m_base = self._aut_gp_and_can_label(m)
     3945                    from sage.rings.arith import factorial
     3946                    if True:#size*factorial(n-B.ncols) == m_size:
     3947#                        print 'in if'
     3948#                        print 'm_aut_gp_gens:', m_aut_gp_gens
     3949                        if len(m_aut_gp_gens) == 0:
     3950                            aut_m = PermutationGroup([()])
     3951                        else:
     3952                            aut_m = PermutationGroup([PermutationGroupElement([a+1 for a in g]) for g in m_aut_gp_gens])
     3953#                        print 'aut_m:', aut_m
     3954#                        print 'aug_aut_gp_gens:', aug_aut_gp_gens
     3955                        if len(aug_aut_gp_gens) == 0:
     3956                            aut_B_aug = PermutationGroup([()])
     3957                        else:
     3958                            aut_B_aug = PermutationGroup([PermutationGroupElement([a+1 for a in g]) for g in aug_aut_gp_gens])
     3959#                        print 'aut_B_aug:', aut_B_aug
     3960                        H = aut_m.__interface[gap].Intersection2(aut_B_aug.__interface[gap])
     3961#                        print 'H:', H
     3962                        rt_transversal = list(gap('List(RightTransversal( %s,%s ));'\
     3963                          %(str(aut_B_aug.__interface[gap]),str(H))))
     3964#                        print 'rt_transversal:', rt_transversal
     3965                        rt_transversal = [PermutationGroupElement(g) for g in rt_transversal if str(g) != '()']
     3966                        rt_transversal = [[a-1 for a in g.list()] for g in rt_transversal]
     3967                        rt_transversal = [g + range(len(g), n) for g in rt_transversal]
     3968                        rt_transversal.append(range(n))
     3969#                        print 'rt_transversal:', rt_transversal
     3970                        bingo2 = 0
     3971                        for coset_rep in rt_transversal:
     3972#                            print 'coset_rep:'
     3973#                            print coset_rep
     3974                            hwp = create_word_perm(coset_rep)
     3975                            #hwp = create_inv_word_perm(gwp) # since we want a left transversal
     3976                            #dealloc_word_perm(gwp)
     3977                            bingo2 = 1
     3978                            for j from 0 <= j < B.nrows:
     3979                                temp = permute_word_by_wp(hwp, temp_basis[j])
     3980                                if temp != B.words[temp & gate]:
     3981                                    bingo2 = 0
     3982                                    dealloc_word_perm(hwp)
     3983                                    break
     3984                            if bingo2:
     3985                                dealloc_word_perm(hwp)
     3986                                break
     3987                        if bingo2:
     3988                            from sage.matrix.constructor import Matrix
     3989                            from sage.rings.all import GF
     3990                            M = Matrix(GF(2), B_aug.nrows, B_aug.ncols)
     3991                            for i from 0 <= i < B_aug.ncols:
     3992                                for j from 0 <= j < B_aug.nrows:
     3993                                    M[j,i] = B_aug.is_one(1 << j, i)
     3994                            output.append(M)
     3995#                            print "ACCEPT"
     3996                    dealloc_word_perm(can_lab)
     3997                    dealloc_word_perm(can_lab_inv)
     3998                #...
     3999#                    print '    orbit_checks:'
     4000#                    for temp from 0 <= temp < ((<codeword>1) << orb_chx_size):
     4001#                        print '    ' + ''.join(reversed(Integer(orbit_checks[temp]).binary().zfill(n)))
     4002                    orbits = [word]
     4003                    j = 0
     4004                    while j < len(orbits):
     4005                        for i from 0 <= i < len(aut_gp_gens):
     4006#                            print '        i', i
     4007                            temp = <codeword> orbits[j]
     4008                            temp = permute_word_by_wp(parent_generators[i], temp)
     4009#                            print '        temp:', ''.join(reversed(Integer(temp).binary().zfill(n)))
     4010                            temp ^= B.words[temp & gate]
     4011#                            print '        temp:', ''.join(reversed(Integer(temp).binary().zfill(n)))
     4012                            if temp not in orbits:
     4013                                orbits.append(temp)
     4014                        j += 1
     4015                    for temp in orbits:
     4016                        temp = (temp >> B.nrows) & k_gate
     4017#                        print '        temp:', temp
     4018#                        print '        ', temp >> log_2_radix
     4019#                        print '        ', ((<codeword>1) << (temp & radix_gate))                   
     4020                        orbit_checks[temp >> log_2_radix] |= ((<codeword>1) << (temp & radix_gate))
     4021
     4022
     4023            parity ^= 1
     4024            i = 0
     4025            if not parity:
     4026                while not combo & (1 << i): i += 1
     4027                i += 1
     4028            if i == k: break
     4029            else:
     4030                combo ^= (1 << i)
     4031                word ^= ortho_basis[i]
     4032
     4033        for i from 0 <= i < len(aut_gp_gens):
     4034            dealloc_word_perm(parent_generators[i])
     4035        sage_free(B_can_lab)
     4036        sage_free(parent_generators)
     4037        sage_free(orbit_checks)
     4038        return output
     4039
     4040
     4041
  • sage/coding/linear_code.py

    diff -r f2c1d993c672 -r 8bcb24714f11 sage/coding/linear_code.py
    a b from sage.rings.finite_field import Fini 
    157157from sage.rings.finite_field import FiniteField as GF
    158158from sage.groups.perm_gps.permgroup import PermutationGroup
    159159from sage.matrix.matrix_space import MatrixSpace
     160from sage.matrix.constructor import Matrix
    160161from sage.rings.arith import GCD, rising_factorial, binomial
    161162from sage.groups.all import SymmetricGroup
    162163from sage.misc.misc import prod
    def bounds_minimum_distance(n,k,F): 
    378379    Ldata = gap.eval("Display(data)")
    379380    return Ldata
    380381
     382def self_orthogonal_binary_codes(n, k, b=2, parent=None, BC=None, equal=False,
     383    in_test=None):
     384    """
     385    Returns a Python iterator which generates a complete set of representatives
     386    of all permutation equivalence classes of self-orthogonal binary linear codes
     387    of length in [1..n] and dimension in [1..k].
     388       
     389    INPUT:
     390    n -- maximal length
     391    k -- maximal dimension
     392    b -- require that the generators all have weight divisible by b (if b=2, all
     393        self-orthogonal codes are generated, and if b=4, all doubly even codes
     394        are generated). Must be an even positive integer.
     395    parent -- dafault None, used in recursion
     396    BC -- dafault None, used in recursion
     397    equal -- default False, if True generates only [n, k] codes
     398    in_test -- default None, used in recursion
     399   
     400    EXAMPLES:
     401    Generate all self-orthogonal codes of length up to 7 and dimension up to 3:
     402        sage: for B in self_orthogonal_binary_codes(7,3):
     403        ...    print B
     404        ...
     405        Linear code of length 2, dimension 1 over Finite Field of size 2
     406        Linear code of length 4, dimension 2 over Finite Field of size 2
     407        Linear code of length 6, dimension 3 over Finite Field of size 2
     408        Linear code of length 4, dimension 1 over Finite Field of size 2
     409        Linear code of length 6, dimension 2 over Finite Field of size 2
     410        Linear code of length 6, dimension 2 over Finite Field of size 2
     411        Linear code of length 7, dimension 3 over Finite Field of size 2
     412        Linear code of length 6, dimension 1 over Finite Field of size 2
     413   
     414    Generate all doubly-even codes of length up to 7 and dimension up to 3:
     415        sage: for B in self_orthogonal_binary_codes(7,3,4):
     416        ...    print B; print B.gen_mat()
     417        ...
     418        Linear code of length 4, dimension 1 over Finite Field of size 2
     419        [1 1 1 1]
     420        Linear code of length 6, dimension 2 over Finite Field of size 2
     421        [1 1 1 1 0 0]
     422        [0 1 0 1 1 1]
     423        Linear code of length 7, dimension 3 over Finite Field of size 2
     424        [1 0 1 1 0 1 0]
     425        [0 1 0 1 1 1 0]
     426        [0 0 1 0 1 1 1]
     427
     428    Generate all doubly-even codes of length up to 7 and dimension up to 2:
     429        sage: for B in self_orthogonal_binary_codes(7,2,4):
     430        ...    print B; print B.gen_mat()
     431        Linear code of length 4, dimension 1 over Finite Field of size 2
     432        [1 1 1 1]
     433        Linear code of length 6, dimension 2 over Finite Field of size 2
     434        [1 1 1 1 0 0]
     435        [0 1 0 1 1 1]
     436   
     437    Generate all self-orthogonal codes of length equal to 8 and dimension
     438    equal to 4:
     439        sage: for B in self_orthogonal_binary_codes(8, 4, equal=True):
     440        ...     print B; print B.gen_mat()
     441        Linear code of length 8, dimension 4 over Finite Field of size 2
     442        [1 0 0 1 0 0 0 0]
     443        [0 1 0 0 1 0 0 0]
     444        [0 0 1 0 0 1 0 0]
     445        [0 0 0 0 0 0 1 1]
     446        Linear code of length 8, dimension 4 over Finite Field of size 2
     447        [1 0 0 1 1 0 1 0]
     448        [0 1 0 1 1 1 0 0]
     449        [0 0 1 0 1 1 1 0]
     450        [0 0 0 1 0 1 1 1]
     451
     452    Since all the codes will be self-orthogonal, b must be divisible by 2:
     453        sage: list(self_orthogonal_binary_codes(8, 4, 1, equal=True))
     454        Traceback (most recent call last):
     455        ...
     456        ValueError: b (1) must be a positive even integer.
     457
     458    """
     459    d=int(b)
     460    if d!=b or d%2==1 or d <= 0:
     461        raise ValueError("b (%s) must be a positive even integer."%b)
     462    from binary_code import BinaryCode, BinaryCodeClassifier
     463    if k < 1 or n < 2:
     464        return
     465    if equal:
     466        in_test = lambda M : (M.ncols() - M.nrows()) <= (n-k)
     467        out_test = lambda C : (C.dimension() == k) and (C.length() == n)
     468    else:
     469        in_test = lambda M : True
     470        out_test = lambda C : True
     471    if BC is None:
     472        BC = BinaryCodeClassifier()
     473    if parent is None:
     474        for j in xrange(d, n+1, d):
     475            M = Matrix(GF(2), [[1]*j])
     476            if in_test(M):
     477                for N in self_orthogonal_binary_codes(n, k, d, M, BC, in_test=in_test):
     478                    if out_test(N): yield N
     479    else:
     480        C = LinearCode(parent)
     481        if out_test(C): yield C
     482        if k == parent.nrows():
     483            return
     484        for nn in xrange(parent.ncols()+1, n+1):
     485            if in_test(parent):
     486                for child in BC.generate_children(BinaryCode(parent), nn, d):
     487                    for N in self_orthogonal_binary_codes(n, k, d, child, BC, in_test=in_test):
     488                        if out_test(N): yield N
    381489
    382490########################### linear codes python class #######################
    383491
    class LinearCode(module.Module): 
    440548
    441549    def automorphism_group_binary_code(self):
    442550        r"""
    443         ** Experimental Interface with Robert Miller's program \code{aut_gp_and_can_label}.**
    444 
    445551        This only applies to linear binary codes and returns its
    446552        (permutation) automorphism group. In other words, if
    447553        the code C has length $n$ then it returns the subgroup of the
    class LinearCode(module.Module): 
    454560        EXAMPLES:
    455561            sage: C = HammingCode(3,GF(2))
    456562            sage: G = C.automorphism_group_binary_code(); G
    457             Permutation Group with generators [(2,3)(5,7), (2,5)(3,7), (2,3,7,5)(4,6), (2,4)(6,7), (1,2)(3,4)]
     563            Permutation Group with generators [(3,4)(5,6), (3,5)(4,6), (2,3)(5,7), (1,2)(5,6)]
    458564            sage: G.order()
    459565            168
    460 
    461         WARNING: This is preliminary - can lock up or return RuntimeError.
    462         (This is a bug in the interface, not in \code{aut_gp_and_can_label}.)
    463         For example, "C = QuasiQuadraticResidueCode(11); C.automorphism_group_binary_code()"
    464         causes a lock-up. Please report other bugs to wdjoyner@gmail.com or sage-support.
    465566
    466567        """
    467568        C = self
    class LinearCode(module.Module): 
    13901491        EXAMPLES:
    13911492            sage: C = HammingCode(3,GF(2))
    13921493            sage: G = C.automorphism_group_binary_code(); G
    1393             Permutation Group with generators [(2,3)(5,7), (2,5)(3,7), (2,3,7,5)(4,6), (2,4)(6,7), (1,2)(3,4)]
     1494            Permutation Group with generators [(3,4)(5,6), (3,5)(4,6), (2,3)(5,7), (1,2)(5,6)]
    13941495            sage: g = G("(2,3)(5,7)")
    13951496            sage: Cg = C.permuted_code(g)
    13961497            sage: Cg