Ticket #3112: trac-3112-socodes.patch

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

This is a flat version of the previous three.

  • sage/coding/all.py 

    # HG changeset patch
# User Robert L. Miller <rlm@rlmiller.org>
# Date 1211221535 25200
# Node ID d8cf41c5a6c54e4a954ae8c0f0988e8039407594
# Parent  23d0e6b6a67514295103c8fecd63a9ec0a5c77db
Generate self-orthogonal binary codes.

diff -r 23d0e6b6a675 -r d8cf41c5a6c5 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 23d0e6b6a675 -r d8cf41c5a6c5 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    cdef int put_in_std_form(self)
     27
     28cdef WordPermutation *create_word_perm(object)
     29cdef WordPermutation *create_array_word_perm(int *, int, int)
     30cdef WordPermutation *create_id_word_perm(int)
     31cdef WordPermutation *create_comp_word_perm(WordPermutation *, WordPermutation *)
     32cdef WordPermutation *create_inv_word_perm(WordPermutation *)
     33cdef int dealloc_word_perm(WordPermutation *)
     34cdef codeword permute_word_by_wp(WordPermutation *, codeword)
     35cdef codeword *expand_to_ortho_basis(BinaryCode, int)
    1836
    1937cdef class OrbitPartition:
    2038    cdef int nwords
    … … cdef class PartitionStack: 
    6886    cdef int refine(self, int, int *, int, BinaryCode, int *)
    6987    cdef void clear(self, int)
    7088    cdef int cmp(self, PartitionStack, BinaryCode)
    71     cdef void find_basis(self, int *)
    72     cdef void get_permutation(self, PartitionStack, int *, int *, int *)
     89    cdef int find_basis(self, int *)
     90    cdef void get_permutation(self, PartitionStack, int *, int *)
    7391
    7492cdef class BinaryCodeClassifier:
    7593    cdef int *ham_wts
    … … cdef class BinaryCodeClassifier: 
    84102    cdef int *alpha
    85103    cdef int alpha_size
    86104    cdef int *v, *e
    87     cdef int *aut_gp_gens, *labeling
    88     cdef int aut_gp_index, aut_gens_size
     105    cdef int *aut_gp_gens, *labeling, *base
     106    cdef int aut_gp_index, aut_gens_size, base_size
    89107    cdef object aut_gp_size
    90108
    91109    cdef int Phi_size
    … … cdef class BinaryCodeClassifier: 
    95113
    96114
    97115
    98 
    99 
    100 
    101 
    102 

  • sage/coding/binary_code.pyx 

    diff -r 23d0e6b6a675 -r d8cf41c5a6c5 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
     42
     43cdef enum:
     44    chunk_size = 8
     45
     46cdef inline int min(int a, int b):
     47    if a > b:
     48        return b
     49    else:
     50        return a
    4251
    4352## NOTE - Since most of the functions are used from within the module, cdef'd
    4453## functions come without an underscore, and the def'd equivalents, which are
    4554## essentially only for doctesting and debugging, have underscores.
    4655
    4756cdef int *hamming_weights():
    48     cdef int *ham_wts
     57    cdef int *ham_wts, i
    4958    ham_wts = <int *> sage_malloc( 65536 * sizeof(int) )
    50     if not ham_wts:
     59    if ham_wts is NULL:
    5160        sage_free(ham_wts)
    5261        raise MemoryError("Memory.")
    5362    ham_wts[0] = 0
    … … cdef int *hamming_weights(): 
    6170    for i from 256 <= i < 65536:
    6271        ham_wts[i] = ham_wts[i & 255] + ham_wts[(i>>8) & 255]
    6372    return ham_wts
     73
     74def test_word_perms(t_limit=5.0, n=32):
     75    """
     76    Tests the WordPermutation structs for at least t_limit seconds.
     77   
     78    TESTS:
     79        sage: from sage.coding.binary_code import test_word_perms
     80        sage: test_word_perms()
     81
     82    """
     83    cdef WordPermutation *g, *h, *i
     84    cdef codeword cw1, cw2, cw3
     85    cdef int j
     86    from sage.misc.prandom import randint
     87    from sage.combinat.permutation import Permutations
     88    S = Permutations(range(n))
     89    t = cputime()
     90    while cputime(t) < t_limit:
     91        word = [randint(0,1) for _ in xrange(n)]
     92        cw1 = 0
     93        for j from 0 <= j < n:
     94            cw1 += (<codeword>word[j]) << (<codeword>j)
     95        # test create_word_perm
     96        gg = S.random_element()
     97        g = create_word_perm(gg)
     98        word2 = [0]*n
     99        for j from 0 <= j < n:
     100            word2[gg[j]] = word[j]
     101        cw2 = permute_word_by_wp(g, cw1)
     102        cw3 = 0
     103        for j from 0 <= j < n:
     104            cw3 += (<codeword>word2[j]) << (<codeword>j)
     105        if cw3 != cw2:
     106            print "ERROR1"
     107        dealloc_word_perm(g)
     108        # test create_comp_word_perm
     109        gg = S.random_element()
     110        hh = S.random_element()
     111        g = create_word_perm(gg)
     112        h = create_word_perm(hh)
     113        i = create_comp_word_perm(g, h)
     114        word2 = [0]*n
     115        for j from 0 <= j < n:
     116            word2[gg[hh[j]]] = word[j]
     117        cw2 = permute_word_by_wp(i, cw1)
     118        cw3 = 0
     119        for j from 0 <= j < n:
     120            cw3 += (<codeword>word2[j]) << (<codeword>j)
     121        if cw3 != cw2:
     122            print "ERROR2"
     123        dealloc_word_perm(g)
     124        dealloc_word_perm(h)
     125        dealloc_word_perm(i)
     126        # test create_inv_word_perm
     127        gg = S.random_element()
     128        g = create_word_perm(gg)
     129        h = create_inv_word_perm(g)
     130        cw2 = permute_word_by_wp(g, cw1)
     131        cw2 = permute_word_by_wp(h, cw2)
     132        if cw1 != cw2:
     133            print "ERROR3a"
     134        cw2 = permute_word_by_wp(h, cw1)
     135        cw2 = permute_word_by_wp(g, cw2)
     136        if cw1 != cw2:
     137            print "ERROR3b"
     138        dealloc_word_perm(g)
     139        dealloc_word_perm(h)
     140
     141cdef WordPermutation *create_word_perm(object list_perm):
     142    r"""
     143    Create a word permutation from a Python list permutation L, i.e. such that
     144    $i \mapsto L[i]$.
     145    """
     146    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     147    cdef codeword *images_i, image
     148    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     149    if word_perm is NULL:
     150        raise RuntimeError("Error allocating memory.")
     151    word_perm.degree = len(list_perm)
     152    list_perm = copy(list_perm)
     153    while num_chunks*chunk_size < word_perm.degree:
     154        num_chunks += 1
     155    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     156    if word_perm.images is NULL:
     157        sage_free(word_perm)
     158        raise RuntimeError("Error allocating memory.")
     159    word_perm.chunk_num = num_chunks
     160    words_per_chunk = 1 << chunk_size
     161    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     162    list_perm += range(len(list_perm), chunk_size*num_chunks)
     163    word_perm.chunk_words = words_per_chunk
     164    for i from 0 <= i < num_chunks:
     165        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     166        if images_i is NULL:
     167            for j from 0 <= j < i:
     168                sage_free(word_perm.images[j])
     169            sage_free(word_perm.images)
     170            sage_free(word_perm)
     171            raise RuntimeError("Error allocating memory.")           
     172        word_perm.images[i] = images_i
     173        for j from 0 <= j < chunk_size:
     174            images_i[1 << j] = (<codeword>1) << list_perm[chunk_size*i + j]
     175        image = <codeword> 0
     176        parity = 0
     177        comb = 0
     178        while 1:
     179            images_i[comb] = image
     180            parity ^= 1
     181            j = 0
     182            if not parity:
     183                while not comb & (1 << j): j += 1
     184                j += 1
     185            if j == chunk_size: break
     186            else:
     187                comb ^= (1 << j)
     188                image ^= images_i[1 << j]
     189    return word_perm
     190
     191cdef WordPermutation *create_array_word_perm(int *array, int start, int degree):
     192    """
     193    Create a word permutation of a given degree from a C array, starting at start.
     194    """
     195    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     196    cdef codeword *images_i, image
     197    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     198    if word_perm is NULL:
     199        raise RuntimeError("Error allocating memory.")
     200    word_perm.degree = degree
     201    while num_chunks*chunk_size < word_perm.degree:
     202        num_chunks += 1
     203    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     204    if word_perm.images is NULL:
     205        sage_free(word_perm)
     206        raise RuntimeError("Error allocating memory.")
     207    word_perm.chunk_num = num_chunks
     208    words_per_chunk = 1 << chunk_size
     209    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     210    word_perm.chunk_words = words_per_chunk
     211    for i from 0 <= i < num_chunks:
     212        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     213        if images_i is NULL:
     214            for j from 0 <= j < i:
     215                sage_free(word_perm.images[j])
     216            sage_free(word_perm.images)
     217            sage_free(word_perm)
     218            raise RuntimeError("Error allocating memory.")           
     219        word_perm.images[i] = images_i
     220        for j from 0 <= j < chunk_size:
     221            images_i[1 << j] = (<codeword>1) << array[start + chunk_size*i + j]
     222        image = <codeword> 0
     223        parity = 0
     224        comb = 0
     225        while 1:
     226            images_i[comb] = image
     227            parity ^= 1
     228            j = 0
     229            if not parity:
     230                while not comb & (1 << j): j += 1
     231                j += 1
     232            if j == chunk_size: break
     233            else:
     234                comb ^= (1 << j)
     235                image ^= images_i[1 << j]
     236    return word_perm
     237
     238cdef WordPermutation *create_id_word_perm(int degree):
     239    """
     240    Create the identity word permutation of degree degree.
     241    """
     242    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     243    cdef codeword *images_i, image
     244    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     245    if word_perm is NULL:
     246        raise RuntimeError("Error allocating memory.")
     247    word_perm.degree = degree
     248    while num_chunks*chunk_size < degree:
     249        num_chunks += 1
     250    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     251    if word_perm.images is NULL:
     252        sage_free(word_perm)
     253        raise RuntimeError("Error allocating memory.")
     254    word_perm.chunk_num = num_chunks
     255    words_per_chunk = 1 << chunk_size
     256    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     257    word_perm.chunk_words = words_per_chunk
     258    for i from 0 <= i < num_chunks:
     259        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     260        if images_i is NULL:
     261            for j from 0 <= j < i:
     262                sage_free(word_perm.images[j])
     263            sage_free(word_perm.images)
     264            sage_free(word_perm)
     265            raise RuntimeError("Error allocating memory.")
     266        word_perm.images[i] = images_i
     267        for j from 0 <= j < chunk_size:
     268            images_i[1 << j] = (<codeword>1) << (chunk_size*i + j)
     269        image = <codeword> 0
     270        parity = 0
     271        comb = 0
     272        while 1:
     273            images_i[comb] = image
     274            parity ^= 1
     275            j = 0
     276            if not parity:
     277                while not comb & (1 << j): j += 1
     278                j += 1
     279            if j == chunk_size: break
     280            else:
     281                comb ^= (1 << j)
     282                image ^= images_i[1 << j]
     283    return word_perm
     284
     285cdef WordPermutation *create_comp_word_perm(WordPermutation *g, WordPermutation *h):
     286    r"""
     287    Create the composition of word permutations $g \circ h$.
     288    """
     289    cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
     290    cdef codeword *images_i, image
     291    cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
     292    if word_perm is NULL:
     293        raise RuntimeError("Error allocating memory.")
     294    word_perm.degree = g.degree
     295    while num_chunks*chunk_size < word_perm.degree:
     296        num_chunks += 1
     297    word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
     298    if word_perm.images is NULL:
     299        sage_free(word_perm)
     300        raise RuntimeError("Error allocating memory.")
     301    word_perm.chunk_num = num_chunks
     302    words_per_chunk = 1 << chunk_size
     303    word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
     304    word_perm.chunk_words = words_per_chunk
     305    for i from 0 <= i < num_chunks:
     306        images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
     307        if images_i is NULL:
     308            for j from 0 <= j < i:
     309                sage_free(word_perm.images[j])
     310            sage_free(word_perm.images)
     311            sage_free(word_perm)
     312            raise RuntimeError("Error allocating memory.")           
     313        word_perm.images[i] = images_i
     314        for j from 0 <= j < chunk_size:
     315            image = (<codeword>1) << (chunk_size*i + j)
     316            image = permute_word_by_wp(h, image)
     317            image = permute_word_by_wp(g, image)
     318            images_i[1 << j] = image
     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_inv_word_perm(WordPermutation *g):
     336    r"""
     337    Create the inverse $g^{-1}$ of the word permutation of $g$.
     338    """
     339    cdef int i, j, *array = <int *> sage_malloc( g.degree * sizeof(int) )
     340    cdef codeword temp
     341    cdef WordPermutation *w
     342    for i from 0 <= i < g.degree:
     343        j = 0
     344        temp = permute_word_by_wp(g, (<codeword>1) << i)
     345        while not ((<codeword>1) << j) & temp:
     346            j += 1
     347        array[j] = i
     348    w = create_array_word_perm(array, 0, g.degree)
     349    sage_free(array)
     350    return w
     351
     352cdef int dealloc_word_perm(WordPermutation *wp):
     353    """
     354    Free the memory used by a word permutation.
     355    """
     356    cdef int i
     357    for i from 0 <= i < wp.chunk_num:
     358        sage_free(wp.images[i])
     359    sage_free(wp.images)
     360    sage_free(wp)
     361
     362cdef codeword permute_word_by_wp(WordPermutation *wp, codeword word):
     363    """
     364    Return the codeword obtained by applying the permutation wp to word.
     365    """
     366    cdef int num_chunks = wp.chunk_num
     367    cdef int i
     368    cdef codeword gate = wp.gate
     369    cdef codeword image = 0
     370    cdef codeword **images = wp.images
     371    for i from 0 <= i < num_chunks:
     372        image += images[i][(word >> i*chunk_size) & gate]
     373    return image
     374
     375cdef codeword *expand_to_ortho_basis(BinaryCode B, int n):
     376    r"""
     377    INPUT:
     378    B -- a BinaryCode
     379    n -- the degree
     380   
     381    OUTPUT:
     382    An array of codewords which represent the expansion of a basis for $B$ to a
     383    basis for $(B^\prime)^\perp$, where $B^\prime = B$ if the all-ones vector 1
     384    is in $B$, otherwise $B^\prime = \text{span}(B,1)$ (note that this guarantees
     385    that all the vectors in the span of the output have even weight).
     386    """
     387    # assumes B is already in standard form
     388    cdef codeword *basis, word = 0, temp, new, pivots = 0, combo, parity
     389    cdef codeword n_gate = (~<codeword>0) >> ( (sizeof(codeword)<<3) - n)
     390    cdef int i, j, m, k = B.nrows, dead, d
     391    cdef WordPermutation *wp
     392    basis = <codeword *> sage_malloc( n * sizeof(codeword) )
     393    if basis is NULL:
     394        raise MemoryError()
     395    for i from 0 <= i < k:
     396        basis[i] = B.basis[i]
     397        word ^= basis[i]
     398    # If 11...1 is already a word of the code,
     399    # then being orthogonal to the code guarantees
     400    # being even weight. Otherwise, add this in.
     401    word = (~word) & n_gate
     402    if word:
     403        basis[k] = word
     404        temp = (<codeword>1) << k
     405        i = k
     406        while not word & temp:
     407            temp = temp << 1
     408            i += 1
     409        for j from 0 <= j < k:
     410            if temp & basis[j]:
     411                basis[j] ^= word
     412        temp += (<codeword>1 << k) - 1
     413        i = k
     414        word = <codeword>1 << k
     415        k += 1
     416    else: # NOTE THIS WILL NEVER HAPPEN AS CURRENTLY SET UP!
     417        temp = (<codeword>1 << k) - 1
     418        i = k
     419        word = <codeword>1 << k
     420    # Now:
     421    # k is the length of the basis so far
     422    j = k
     423    while i < n:
     424        # we are now looking at the ith free variable,
     425        # word has a 1 in the ith place, and
     426        # j is the current row we are putting in basis
     427        new = 0
     428        for m from 0 <= m < k:
     429            if basis[m] & word:
     430                new ^= basis[m]
     431        basis[j] = (new & temp) + word
     432        j += ((word ^ temp) >> i) & 1
     433        i += 1
     434        word = word << 1
     435    temp = (<codeword>1 << B.nrows) - 1
     436    for i from k <= i < n:
     437        basis[i-k] = basis[i] ^ B.words[basis[i] & temp]
     438    k = n-k
     439    i = 0
     440    word = (<codeword>1 << B.nrows)
     441    while i < k and (word & n_gate):
     442        m = i
     443        while m < k and not basis[m] & word:
     444            m += 1
     445        if m < k:
     446            pivots += word
     447            if m != i:
     448                new = basis[i]
     449                basis[i] = basis[m]
     450                basis[m] = new
     451            for j from 0 <= j < i:
     452                if basis[j] & word:
     453                    basis[j] ^= basis[i]
     454            for j from i < j < k:
     455                if basis[j] & word:
     456                    basis[j] ^= basis[i]
     457            i += 1
     458        word = word << 1
     459    for j from i <= j < n:
     460        basis[j] = 0
     461    # now basis is length i
     462    perm = range(B.nrows)
     463    perm_c = []
     464    for j from B.nrows <= j < B.ncols:
     465        if (<codeword>1 << j) & pivots:
     466            perm.append(j)
     467        else:
     468            perm_c.append(j)
     469    perm.extend(perm_c)
     470    perm.extend(range(B.ncols, n))
     471    perm_c = [0]*n
     472    for j from 0 <= j < n:
     473        perm_c[perm[j]] = j
     474    wp = create_word_perm(perm_c)
     475    for j from 0 <= j < i:
     476        basis[j] = permute_word_by_wp(wp, basis[j])
     477    for j from 0 <= j < B.nrows:
     478        B.basis[j] = permute_word_by_wp(wp, B.basis[j])
     479    dealloc_word_perm(wp)
     480    word = 0
     481    parity = 0
     482    combo = 0
     483    while 1:
     484        B.words[combo] = word
     485        parity ^= 1
     486        j = 0
     487        if not parity:
     488            while not combo & (1 << j): j += 1
     489            j += 1
     490        if j == B.nrows: break
     491        else:
     492            combo ^= (1 << j)
     493            word ^= B.basis[j]
     494    return basis
    64495
    65496cdef class BinaryCode:
    66497    """
    … … cdef class BinaryCode: 
    101532
    102533    """
    103534    def __new__(self, arg1, arg2=None):
    104         cdef int nrows, i, j
     535        cdef int nrows, i, j, size
    105536        cdef int nwords, other_nwords, parity, combination
    106537        cdef codeword word, glue_word
    107538        cdef BinaryCode other
    … … cdef class BinaryCode: 
    116547            self.nwords = 1 << nrows
    117548            nwords = self.nwords
    118549        elif isinstance(arg1, BinaryCode):
    119             other = arg1
     550            other = <BinaryCode> arg1
    120551            self.nrows = other.nrows + 1
    121552            glue_word = <codeword> arg2
    122             self.ncols = max( other.ncols , floor(log(arg2,2))+1 )
     553            size = 0
     554            while 0 < ((<codeword> 1) << size) <= glue_word:
     555                size += 1
     556            if other.ncols > size:
     557                self.ncols = other.ncols
     558            else:
     559                self.ncols = size
    123560            other_nwords = other.nwords
    124561            self.nwords = 2 * other_nwords
    125562            nrows = self.nrows
    … … cdef class BinaryCode: 
    131568
    132569        self.words = <codeword *> sage_malloc( nwords * sizeof(int) )
    133570        self.basis = <codeword *> sage_malloc( nrows * sizeof(int) )
    134         if not self.words or not self.basis:
    135             if self.words: sage_free(self.words)
    136             if self.basis: sage_free(self.basis)
     571        if self.words is NULL or self.basis is NULL:
     572            if self.words is not NULL: sage_free(self.words)
     573            if self.basis is not NULL: sage_free(self.basis)
    137574            raise MemoryError("Memory.")
    138575        self_words = self.words
    139576        self_basis = self.basis
    … … cdef class BinaryCode: 
    149586            word = <codeword> 0
    150587            parity = 0
    151588            combination = 0
    152             while True:
     589            while 1:
    153590                self_words[combination] = word
    154591                parity ^= 1
    155592                j = 0
    … … cdef class BinaryCode: 
    176613    def __dealloc__(self):
    177614        sage_free(self.words)
    178615        sage_free(self.basis)
     616
     617    def matrix(self):
     618        """
     619        Returns the generator matrix of the BinaryCode, i.e. the code is the
     620        rowspace of B.matrix().
     621       
     622        EXAMPLE:
     623            sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
     624            sage: from sage.coding.binary_code import *
     625            sage: B = BinaryCode(M)
     626            sage: B.matrix()
     627            [1 1 1 1 0 0]
     628            [0 0 1 1 1 1]
     629       
     630        """
     631        cdef int i, j
     632        from sage.matrix.constructor import matrix
     633        from sage.rings.all import GF
     634        rows = []
     635        for i from 0 <= i < self.nrows:
     636            row = [0]*self.ncols
     637            for j from 0 <= j < self.ncols:
     638                if self.basis[i] & ((<codeword>1) << j):
     639                    row[j] = 1
     640            rows.append(row)
     641        return matrix(GF(2), self.nrows, self.ncols, rows)
    179642
    180643    def print_data(self):
    181644        """
    … … cdef class BinaryCode: 
    368831        cdef int *_word_gamma
    369832        _word_gamma = <int *> sage_malloc(self.nwords * sizeof(int))
    370833        _col_gamma = <int *> sage_malloc(self.ncols * sizeof(int))
    371         if not (_col_gamma and _word_gamma):
    372             if _word_gamma: sage_free(_word_gamma)
    373             if _col_gamma: sage_free(_col_gamma)
     834        if _col_gamma is NULL or _word_gamma is NULL:
     835            if _word_gamma is not NULL: sage_free(_word_gamma)
     836            if _col_gamma is not NULL: sage_free(_col_gamma)
    374837            raise MemoryError("Memory.")
    375838        for i from 0 <= i < self.nwords:
    376839            _word_gamma[i] = word_gamma[i]
    … … cdef class BinaryCode: 
    382845        return result
    383846
    384847    cdef int is_automorphism(self, int *col_gamma, int *word_gamma):
    385         # TODO: optimize? check only basis?
    386848        cdef int i, j, self_nwords = self.nwords, self_ncols = self.ncols
    387         for i from 0 <= i < self_nwords:
     849        i = 1
     850        while i < self_nwords:
    388851            for j from 0 <= j < self_ncols:
    389852                if self.is_one(i, j) != self.is_one(word_gamma[i], col_gamma[j]):
    390853                    return 0
     854            i = i << 1
    391855        return 1
     856
     857    def _put_in_std_form(self):
     858        """
     859        Put the code in binary form, which is defined by an identity matrix on
     860        the left, augmented by a matrix of data.
     861       
     862        EXAMPLE:
     863            sage: from sage.coding.binary_code import *
     864            sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
     865            sage: B = BinaryCode(M); B
     866            Binary [6,2] linear code, generator matrix
     867            [111100]
     868            [001111]
     869            sage: B._put_in_std_form(); B
     870            Binary [6,2] linear code, generator matrix
     871            [101011]
     872            [010111]
     873
     874        """
     875        self.put_in_std_form()
     876   
     877    cdef int put_in_std_form(self):
     878        cdef codeword swap, current = 1, pivots = 0, word
     879        cdef int i, j, k, row = 0, parity, combination
     880        cdef object perm
     881        cdef WordPermutation *wp
     882        while row < self.nrows:
     883            i = row
     884            while i < self.nrows and not self.basis[i] & current:
     885                i += 1
     886            if i < self.nrows:
     887                pivots += current
     888                if i != row:
     889                    swap = self.basis[row]
     890                    self.basis[row] = self.basis[i]
     891                    self.basis[i] = swap
     892                for j from 0 <= j < row:
     893                    if self.basis[j] & current:
     894                        self.basis[j] ^= self.basis[row]
     895                for j from row < j < self.nrows:
     896                    if self.basis[j] & current:
     897                        self.basis[j] ^= self.basis[row]
     898                row += 1
     899            current = current << 1
     900        perm = [0]*self.ncols
     901        j = 0
     902        k = self.nrows
     903        for i from 0 <= i < self.ncols:
     904            if ((<codeword>1) << i) & pivots:
     905                perm[i] = j
     906                j += 1
     907            else:
     908                perm[i] = k
     909                k += 1
     910        wp = create_word_perm(perm)
     911        for i from 0 <= i < self.nrows:
     912            self.basis[i] = permute_word_by_wp(wp, self.basis[i])
     913        dealloc_word_perm(wp)
     914        word = 0
     915        parity = 0
     916        combination = 0
     917        while 1:
     918            self.words[combination] = word
     919            parity ^= 1
     920            j = 0
     921            if not parity:
     922                while not combination & (1 << j): j += 1
     923                j += 1
     924            if j == self.nrows: break
     925            else:
     926                combination ^= (1 << j)
     927                word ^= self.basis[j]
     928
    392929
    393930cdef class OrbitPartition:
    394931    """
    … … cdef class OrbitPartition: 
    401938    * http://en.wikipedia.org/wiki/Disjoint-set_data_structure
    402939
    403940    """
    404     def __new__(self, nrows, ncols):
     941    def __new__(self, int nrows, int ncols):
    405942        cdef int col
    406943        cdef int nwords, word
    407944        nwords = (1 << nrows)
    … … cdef class OrbitPartition: 
    415952        self.col_rank =         <int *> sage_malloc( ncols * sizeof(int) )
    416953        self.col_min_cell_rep = <int *> sage_malloc( ncols * sizeof(int) )
    417954        self.col_size =         <int *> sage_malloc( ncols * sizeof(int) )
    418         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):
    419             if self.wd_parent: sage_free(self.wd_parent)
    420             if self.wd_rank: sage_free(self.wd_rank)
    421             if self.wd_min_cell_rep: sage_free(self.wd_min_cell_rep)
    422             if self.wd_size: sage_free(self.wd_size)
    423             if self.col_parent:       sage_free(self.col_parent)
    424             if self.col_rank:         sage_free(self.col_rank)
    425             if self.col_min_cell_rep: sage_free(self.col_min_cell_rep)
    426             if self.col_size:         sage_free(self.col_size)
     955        if self.wd_parent is NULL or self.wd_rank is NULL or self.wd_min_cell_rep is NULL \
     956        or self.wd_size is NULL or self.col_parent is NULL or self.col_rank is NULL \
     957        or self.col_min_cell_rep is NULL or self.col_size is NULL:
     958            if self.wd_parent is not NULL:        sage_free(self.wd_parent)
     959            if self.wd_rank is not NULL:          sage_free(self.wd_rank)
     960            if self.wd_min_cell_rep is not NULL:  sage_free(self.wd_min_cell_rep)
     961            if self.wd_size is not NULL:          sage_free(self.wd_size)
     962            if self.col_parent is not NULL:       sage_free(self.col_parent)
     963            if self.col_rank is not NULL:         sage_free(self.col_rank)
     964            if self.col_min_cell_rep is not NULL: sage_free(self.col_min_cell_rep)
     965            if self.col_size is not NULL:         sage_free(self.col_size)
    427966            raise MemoryError("Memory.")
    428967        for word from 0 <= word < nwords:
    429968            self.wd_parent[word] = word
    … … cdef class OrbitPartition: 
    6611200        cdef int *_wd_gamma
    6621201        _wd_gamma = <int *> sage_malloc(self.nwords * sizeof(int))
    6631202        _col_gamma = <int *> sage_malloc(self.ncols * sizeof(int))
    664         if not (_col_gamma and _wd_gamma):
    665             if _wd_gamma: sage_free(_wd_gamma)
    666             if _col_gamma: sage_free(_col_gamma)
     1203        if _col_gamma is NULL or _wd_gamma is NULL:
     1204            if _wd_gamma is not NULL: sage_free(_wd_gamma)
     1205            if _col_gamma is not NULL: sage_free(_col_gamma)
    6671206            raise MemoryError("Memory.")
    6681207        for i from 0 <= i < self.nwords:
    6691208            _wd_gamma[i] = wd_gamma[i]
    … … cdef class PartitionStack: 
    7081247        sizeof_int = sizeof(int)
    7091248
    7101249        try:
    711             self.nrows = int(arg1)
     1250            self.nrows = <int> arg1
    7121251            self.nwords = 1 << self.nrows
    713             self.ncols = int(arg2)
     1252            self.ncols = <int> arg2
    7141253        except:
    7151254            other = arg1
    7161255            self.nrows = other.nrows
    … … cdef class PartitionStack: 
    7341273        self.wd_counts =  <int *> sage_malloc( (self.ncols+1)  * sizeof_int )
    7351274        self.wd_output =  <int *> sage_malloc( self.nwords * sizeof_int )
    7361275
    737         if not (self.wd_ents  and self.wd_lvls    and self.col_ents   and self.col_lvls  \
    738             and self.col_degs and self.col_counts and self.col_output \
    739             and self.wd_degs  and self.wd_counts  and self.wd_output):
    740             if self.wd_ents:         sage_free(self.wd_ents)
    741             if self.wd_lvls:         sage_free(self.wd_lvls)
    742             if self.col_ents:        sage_free(self.col_ents)
    743             if self.col_lvls:        sage_free(self.col_lvls)
    744             if self.col_degs:        sage_free(self.col_degs)
    745             if self.col_counts:      sage_free(self.col_counts)
    746             if self.col_output:      sage_free(self.col_output)
    747             if self.wd_degs:         sage_free(self.wd_degs)
    748             if self.wd_counts:       sage_free(self.wd_counts)
    749             if self.wd_output:       sage_free(self.wd_output)
     1276        if self.wd_ents is NULL or self.wd_lvls is NULL or self.col_ents is NULL \
     1277        or self.col_lvls is NULL or self.col_degs is NULL or self.col_counts is NULL \
     1278        or self.col_output is NULL or self.wd_degs is NULL or self.wd_counts is NULL \
     1279        or self.wd_output is NULL:
     1280            if self.wd_ents is not NULL:    sage_free(self.wd_ents)
     1281            if self.wd_lvls is not NULL:    sage_free(self.wd_lvls)
     1282            if self.col_ents is not NULL:   sage_free(self.col_ents)
     1283            if self.col_lvls is not NULL:   sage_free(self.col_lvls)
     1284            if self.col_degs is not NULL:   sage_free(self.col_degs)
     1285            if self.col_counts is not NULL: sage_free(self.col_counts)
     1286            if self.col_output is not NULL: sage_free(self.col_output)
     1287            if self.wd_degs is not NULL:    sage_free(self.wd_degs)
     1288            if self.wd_counts is not NULL:  sage_free(self.wd_counts)
     1289            if self.wd_output is not NULL:  sage_free(self.wd_output)
    7501290            raise MemoryError("Memory.")
    7511291
    7521292        nwords = self.nwords
    … … cdef class PartitionStack: 
    12681808        cdef int min_is_col = 1, radix = self.radix
    12691809        cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls
    12701810        cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents
    1271         while True:
     1811        while 1:
    12721812            if self_col_lvls[i] <= k:
    12731813                if i != j and min > i - j + 1:
    12741814                    min = i - j + 1
    … … cdef class PartitionStack: 
    12761816                j = i + 1
    12771817            if self_col_lvls[i] == -1: break
    12781818            i += 1
    1279         i = 0; j = 0
    1280         while True:
    1281             if self_wd_lvls[i] <= k:
    1282                 if i != j and min > i - j + 1:
    1283                     min = i - j + 1
    1284                     min_is_col = 0
    1285                     location = j
    1286                 j = i + 1
    1287             if self_wd_lvls[i] == -1: break
    1288             i += 1
     1819#        i = 0; j = 0
     1820#        while 1:
     1821#            if self_wd_lvls[i] <= k:
     1822#                if i != j and min > i - j + 1:
     1823#                    min = i - j + 1
     1824#                    min_is_col = 0
     1825#                    location = j
     1826#                j = i + 1
     1827#            if self_wd_lvls[i] == -1: break
     1828#            i += 1
    12891829        # location now points to the beginning of the first, smallest,
    12901830        # nontrivial cell
    12911831        j = location
    … … cdef class PartitionStack: 
    12961836        for i from 0 <= i < ell:
    12971837            W[start+i] = 0
    12981838        if min_is_col:
    1299             while True:
     1839            while 1:
    13001840                if self_col_lvls[j] <= k: break
    13011841                j += 1
    13021842            # j now points to the last element of the cell
    … … cdef class PartitionStack: 
    13061846                i += 1
    13071847            return self_col_ents[location]
    13081848        else:
    1309             while True:
     1849            while 1:
    13101850                if self_wd_lvls[j] <= k: break
    13111851                j += 1
    13121852            # j now points to the last element of the cell
    … … cdef class PartitionStack: 
    15682108    cdef int col_degree(self, BinaryCode CG, int col, int wd_ptr, int k):
    15692109        cdef int i = 0
    15702110        cdef int *self_wd_lvls = self.wd_lvls, *self_wd_ents = self.wd_ents
    1571         while True:
     2111        while 1:
    15722112            if CG.is_one(self_wd_ents[wd_ptr], col): i += 1
    15732113            if self_wd_lvls[wd_ptr] > k: wd_ptr += 1
    15742114            else: break
    … … cdef class PartitionStack: 
    18052345        cdef int i, alpha_length = len(alpha)
    18062346        cdef int *_alpha = <int *> sage_malloc( (self.nwords + self.ncols) * sizeof(int) )
    18072347        cdef int *ham_wts = hamming_weights()
    1808         if not _alpha:
    1809             sage_free(_alpha)
     2348        if _alpha is NULL:
    18102349            raise MemoryError("Memory.")
    18112350        for i from 0 <= i < alpha_length:
    18122351            if alpha[i][0]:
    … … cdef class PartitionStack: 
    18362375#                    print self
    18372376                    i = j; s = 0
    18382377                    invariant += 8
    1839                     while True:
     2378                    while 1:
    18402379#                        print 'col_i', self_col_ents[i]
    18412380#                        print 'alpha[m]^flag', alpha[m]^flag
    18422381                        self_col_degs[i-j] = self.col_degree(CG, self_col_ents[i], alpha[m]^flag, k)
    … … cdef class PartitionStack: 
    18562395                                break
    18572396                            q += 1
    18582397                        r = j
    1859                         while True:
     2398                        while 1:
    18602399                            if r == j or self.col_lvls[r-1] == k:
    18612400                                if r != t:
    18622401                                    alpha[alpha_length] = r
    … … cdef class PartitionStack: 
    18732412#                    print self
    18742413                    i = j; s = 0
    18752414                    invariant += 64
    1876                     while True:
     2415                    while 1:
    18772416#                        print 'i', i
    18782417                        self_wd_degs[i-j] = self.wd_degree(CG, self_wd_ents[i], alpha[m], k, ham_wts)
    18792418#                        print 'deg', self_wd_degs[i-j]
    … … cdef class PartitionStack: 
    18932432                            q += 1
    18942433                        j ^= flag
    18952434                        r = j
    1896                         while True:
     2435                        while 1:
    18972436                            if r == j or self.wd_lvls[r-1] == k:
    18982437                                if r != t_w:
    18992438                                    alpha[alpha_length] = r^flag
    … … cdef class PartitionStack: 
    19042443                        invariant += (i-j)
    19052444                    j = i
    19062445            m += 1
    1907         if invariant != -1:
    1908             return invariant
    1909         else:
    1910             return 0
     2446        return invariant
    19112447
    19122448    def _clear(self, k):
    19132449        """
    … … cdef class PartitionStack: 
    19952531        cdef int *self_wd_ents = self.wd_ents
    19962532        cdef codeword *CG_words = CG.words
    19972533        cdef int i, j, l, m, span = 1, ncols = self.ncols, nwords = self.nwords
    1998         for i from 0 <= i < nwords: # TODO: probably don't need to check i == 0 here!
     2534        for i from 0 < i < nwords:
    19992535            for j from 0 <= j < ncols:
    20002536                l = CG.is_one(self.wd_ents[i], self.col_ents[j])
    20012537                m = CG.is_one(other.wd_ents[i], other.col_ents[j])
    … … cdef class PartitionStack: 
    20652601        self.find_basis(ham_wts)
    20662602        sage_free(ham_wts)
    20672603
    2068     cdef void find_basis(self, int *ham_wts):
     2604    cdef int find_basis(self, int *ham_wts):
    20692605        cdef int i = 0, j, k, nwords = self.nwords, weight, basis_elts = 0, nrows = self.nrows
    20702606        cdef int *self_wd_ents = self.wd_ents
    2071         if not self.basis_locations:
    2072             self.basis_locations = <int *> sage_malloc( nrows * sizeof(int) )
    2073         if not self.basis_locations:
     2607        if self.basis_locations is NULL:
     2608            self.basis_locations = <int *> sage_malloc( 2 * nrows * sizeof(int) )
     2609        if self.basis_locations is NULL:
    20742610            raise MemoryError("Memory.")
    20752611        while i < nwords:
    20762612            j = self_wd_ents[i]
    … … cdef class PartitionStack: 
    20832619                self.basis_locations[k] = i
    20842620                if basis_elts == nrows: break
    20852621            i += 1
     2622        for i from 0 <= i < nrows:
     2623            self.basis_locations[nrows + i] = self_wd_ents[1 << i]
    20862624
    20872625    def _get_permutation(self, other):
    20882626        """
    … … cdef class PartitionStack: 
    21252663
    21262664        """
    21272665        cdef int i
    2128         cdef int *ham_wts = hamming_weights()
    21292666        cdef int *word_g = <int *> sage_malloc( self.nwords * sizeof(int) )
    21302667        cdef int *col_g = <int *> sage_malloc( self.ncols * sizeof(int) )
    2131         if not (word_g and col_g):
    2132             if word_g: sage_free(word_g)
    2133             if col_g: sage_free(col_g)
     2668        if word_g is NULL or col_g is NULL:
     2669            if word_g is not NULL: sage_free(word_g)
     2670            if col_g is not NULL: sage_free(col_g)
    21342671            raise MemoryError("Memory.")
    2135         self.get_permutation(other, word_g, col_g, ham_wts)
    2136         sage_free(ham_wts)
     2672        self.get_permutation(other, word_g, col_g)
    21372673        word_l = [word_g[i] for i from 0 <= i < self.nwords]
    21382674        col_l = [col_g[i] for i from 0 <= i < self.ncols]
    21392675        sage_free(word_g)
    21402676        sage_free(col_g)
    21412677        return word_l, col_l
    21422678
    2143     cdef void get_permutation(self, PartitionStack other, int *word_gamma, int *col_gamma, int *ham_wts):
     2679    cdef void get_permutation(self, PartitionStack other, int *word_gamma, int *col_gamma):
    21442680        cdef int i
    21452681        cdef int *self_wd_ents = self.wd_ents, *other_wd_ents = other.wd_ents
    21462682        cdef int *self_col_ents = self.col_ents, *other_col_ents = other.col_ents
    … … cdef class BinaryCodeClassifier: 
    21532689cdef class BinaryCodeClassifier:
    21542690
    21552691    def __new__(self):
    2156         self.radix = sizeof(int) << 3
     2692        self.radix = sizeof(codeword) << 3
    21572693        self.ham_wts = hamming_weights()
    21582694        self.L = 100 # memory limit for Phi and Omega- multiply by 8KB
    21592695        self.aut_gens_size = self.radix * 100
    … … cdef class BinaryCodeClassifier: 
    21682704        self.Omega =   <unsigned int *> sage_malloc( self.Phi_size * self.L         * sizeof(unsigned int) )
    21692705        self.W =       <unsigned int *> sage_malloc( self.Phi_size * self.radix * 2 * sizeof(unsigned int) )
    21702706
    2171         self.aut_gp_gens = <int *> sage_malloc( self.aut_gens_size             * sizeof(int) )
    2172         self.c_gamma =     <int *> sage_malloc( self.radix                     * sizeof(int) )
    2173         self.labeling =    <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2174         self.Lambda1 =     <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2175         self.Lambda2 =     <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2176         self.Lambda3 =     <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2177         self.v =           <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
    2178         self.e =           <int *> sage_malloc( self.radix * 2                 * sizeof(int) )
     2707        self.base =        <int *> sage_malloc( self.radix          * sizeof(int) )
     2708        self.aut_gp_gens = <int *> sage_malloc( self.aut_gens_size  * sizeof(int) )
     2709        self.c_gamma =     <int *> sage_malloc( self.radix          * sizeof(int) )
     2710        self.labeling =    <int *> sage_malloc( self.radix * 3      * sizeof(int) )
     2711        self.Lambda1 =     <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2712        self.Lambda2 =     <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2713        self.Lambda3 =     <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2714        self.v =           <int *> sage_malloc( self.radix * 2      * sizeof(int) )
     2715        self.e =           <int *> sage_malloc( self.radix * 2      * sizeof(int) )
    21792716
    2180         if not (self.Phi and self.Omega and self.W and self.Lambda1 and self.Lambda2 and self.Lambda3 \
    2181             and self.w_gamma and self.c_gamma and self.alpha and self.v and self.e and self.aut_gp_gens \
    2182             and self.labeling):
    2183             if self.Phi:          sage_free(self.Phi)
    2184             if self.Omega:        sage_free(self.Omega)
    2185             if self.W:            sage_free(self.W)
    2186             if self.Lambda1:      sage_free(self.Lambda1)
    2187             if self.Lambda2:      sage_free(self.Lambda2)
    2188             if self.Lambda3:      sage_free(self.Lambda3)
    2189             if self.w_gamma:      sage_free(self.w_gamma)
    2190             if self.c_gamma:      sage_free(self.c_gamma)
    2191             if self.alpha:        sage_free(self.alpha)
    2192             if self.v:            sage_free(self.v)
    2193             if self.e:            sage_free(self.e)
    2194             if self.aut_gp_gens:  sage_free(self.aut_gp_gens)
    2195             if self.labeling:     sage_free(self.labeling)
     2717        if self.Phi is NULL or self.Omega is NULL or self.W is NULL or self.Lambda1 is NULL \
     2718        or self.Lambda2 is NULL or self.Lambda3 is NULL or self.w_gamma is NULL \
     2719        or self.c_gamma is NULL or self.alpha is NULL or self.v is NULL or self.e is NULL \
     2720        or self.aut_gp_gens is NULL or self.labeling is NULL or self.base is NULL:
     2721            if self.Phi is not NULL:          sage_free(self.Phi)
     2722            if self.Omega is not NULL:        sage_free(self.Omega)
     2723            if self.W is not NULL:            sage_free(self.W)
     2724            if self.Lambda1 is not NULL:      sage_free(self.Lambda1)
     2725            if self.Lambda2 is not NULL:      sage_free(self.Lambda2)
     2726            if self.Lambda3 is not NULL:      sage_free(self.Lambda3)
     2727            if self.w_gamma is not NULL:      sage_free(self.w_gamma)
     2728            if self.c_gamma is not NULL:      sage_free(self.c_gamma)
     2729            if self.alpha is not NULL:        sage_free(self.alpha)
     2730            if self.v is not NULL:            sage_free(self.v)
     2731            if self.e is not NULL:            sage_free(self.e)
     2732            if self.aut_gp_gens is not NULL:  sage_free(self.aut_gp_gens)
     2733            if self.labeling is not NULL:     sage_free(self.labeling)
     2734            if self.base is not NULL:         sage_free(self.base)
    21962735            raise MemoryError("Memory.")
    21972736
    21982737    def __dealloc__(self):
    2199         if self.ham_wts: sage_free(self.ham_wts)
    2200         if self.Phi:     sage_free(self.Phi)
    2201         if self.Omega:   sage_free(self.Omega)
    2202         if self.W:       sage_free(self.W)
    2203         if self.Lambda1: sage_free(self.Lambda1)
    2204         if self.Lambda2: sage_free(self.Lambda2)
    2205         if self.Lambda3: sage_free(self.Lambda3)
    2206         if self.c_gamma: sage_free(self.c_gamma)
    2207 
    2208         if self.w_gamma:   sage_free(self.w_gamma)
    2209         if self.alpha:     sage_free(self.alpha)
    2210 
    2211         if self.v:           sage_free(self.v)
    2212         if self.e:           sage_free(self.e)
    2213         if self.aut_gp_gens: sage_free(self.aut_gp_gens)
    2214         if self.labeling:    sage_free(self.labeling)
     2738        sage_free(self.ham_wts)
     2739        sage_free(self.Phi)
     2740        sage_free(self.Omega)
     2741        sage_free(self.W)
     2742        sage_free(self.Lambda1)
     2743        sage_free(self.Lambda2)
     2744        sage_free(self.Lambda3)
     2745        sage_free(self.c_gamma)
     2746        sage_free(self.w_gamma)
     2747        sage_free(self.alpha)
     2748        sage_free(self.v)
     2749        sage_free(self.e)
     2750        sage_free(self.aut_gp_gens)
     2751        sage_free(self.labeling)
     2752        sage_free(self.base)
    22152753
    22162754    cdef void record_automorphism(self, int *gamma, int ncols):
    22172755        cdef int i, j
    22182756        if self.aut_gp_index + ncols > self.aut_gens_size:
    22192757            self.aut_gens_size *= 2
    22202758            self.aut_gp_gens = <int *> sage_realloc( self.aut_gp_gens, self.aut_gens_size * sizeof(int) )
    2221             if not self.aut_gp_gens:
     2759            if self.aut_gp_gens is NULL:
    22222760                raise MemoryError("Memory.")
    22232761        j = self.aut_gp_index
    22242762        for i from 0 <= i < ncols:
    … … cdef class BinaryCodeClassifier: 
    22342772            verbosity - a nonnegative integer
    22352773       
    22362774        OUTPUT:
    2237             a tuple, (gens, labeling, size)
     2775            a tuple, (gens, labeling, size, base)
    22382776            gens -- a list of permutations (in list form) representing generators
    22392777                of the permutation automorphism group of the code CC.
    22402778            labeling -- a permutation representing the canonical labeling of the
    2241                 code. mostly for internal use; if the dimension of the code is k
    2242                 and the degree (number of columns) is n, then the first n entries
    2243                 describe the relabeling on the columns, and the next k describe
    2244                 where the basis is sent.
     2779                code. mostly for internal use; entries describe the relabeling
     2780                on the columns.
    22452781            size -- the order of the automorphism group.
     2782            base -- a set of cols whose action determines the action on all cols
    22462783       
    22472784        EXAMPLES:
    22482785            sage: import sage.coding.binary_code
    … … cdef class BinaryCodeClassifier: 
    22562793            ... [0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],\
    22572794            ... [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]])
    22582795            sage: B = BinaryCode(M)
    2259             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2796            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    22602797            sage: S = SymmetricGroup(M.ncols())
    22612798            sage: L = [S([x+1 for x in g]) for g in gens]
    22622799            sage: PermutationGroup(L).order()
    … … cdef class BinaryCodeClassifier: 
    22702807            ... [0,0,0,0,0,1,0,1,0,0,0,1,1,1,1,1,1],\
    22712808            ... [0,0,0,1,1,0,0,0,0,1,1,0,1,1,0,1,1]])
    22722809            sage: B = BinaryCode(M)
    2273             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2810            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    22742811            sage: S = SymmetricGroup(M.ncols())
    22752812            sage: L = [S([x+1 for x in g]) for g in gens]
    22762813            sage: PermutationGroup(L).order()
    … … cdef class BinaryCodeClassifier: 
    22882825            ... [0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0],\
    22892826            ... [0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1]])
    22902827            sage: B = BinaryCode(M)
    2291             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2828            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    22922829            sage: S = SymmetricGroup(M.ncols())
    22932830            sage: L = [S([x+1 for x in g]) for g in gens]
    22942831            sage: PermutationGroup(L).order()
    … … cdef class BinaryCodeClassifier: 
    23092846            ... [1,0,0,1,0,1,1,1,0,0,0,0,1,0,0,1,0,0,0,1,1,1],\
    23102847            ... [0,0,1,0,1,1,1,0,0,0,1,1,0,0,1,0,0,0,1,1,1,0]])
    23112848            sage: B = BinaryCode(M)
    2312             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2849            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    23132850            sage: S = SymmetricGroup(M.ncols())
    23142851            sage: L = [S([x+1 for x in g]) for g in gens]
    23152852            sage: PermutationGroup(L).order()
    … … cdef class BinaryCodeClassifier: 
    23192856
    23202857            sage: B = BinaryCode(Matrix(GF(2),[[1,0,1],[0,1,1]]))
    23212858            sage: BC._aut_gp_and_can_label(B)
    2322             ([[0, 2, 1], [1, 0, 2]], [0, 1, 2, 2, 1], 6)
     2859            ([[0, 2, 1], [1, 0, 2]], [0, 1, 2], 6, [0, 1])
    23232860
    23242861            sage: B = BinaryCode(Matrix(GF(2),[[1,1,1,1]]))
    23252862            sage: BC._aut_gp_and_can_label(B)
    2326             ([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], [0, 1, 2, 3, 1], 24)
     2863            ([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], [0, 1, 2, 3], 24, [0, 1, 2])
    23272864
    23282865            sage: B = BinaryCode(Matrix(GF(2),[[0,0,0,0,0,0,0,0,0,0,0,0,0,0,1]]))
    2329             sage: gens, labeling, size = BC._aut_gp_and_can_label(B)
     2866            sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
    23302867            sage: size
    23312868            87178291200
    23322869
    … … cdef class BinaryCodeClassifier: 
    23812918            gen = [self.aut_gp_gens[i+j] for j from 0 <= j < C.ncols]
    23822919            py_aut_gp_gens.append(gen)
    23832920            i += C.ncols
    2384         py_labeling = [self.labeling[i] for i from 0 <= i < C.ncols] \
    2385                     + [self.labeling[i+C.ncols] for i from 0 <= i < C.nrows]
     2921        py_labeling = [self.labeling[i] for i from 0 <= i < C.ncols]
     2922        base = []
     2923        for i from 0 <= i < self.radix:
     2924            if self.base[i] == -1:
     2925                break
     2926            base.append(self.base[i])
    23862927        aut_gp_size = self.aut_gp_size
    2387         return py_aut_gp_gens, py_labeling, aut_gp_size
     2928        return py_aut_gp_gens, py_labeling, aut_gp_size, base
    23882929
    23892930    cdef void aut_gp_and_can_label(self, BinaryCode C, int verbosity):
    23902931
    … … cdef class BinaryCodeClassifier: 
    24492990            self.Phi =     <unsigned int *> sage_realloc(self.Phi,   self.Phi_size * self.L         * sizeof(int) )
    24502991            self.Omega =   <unsigned int *> sage_realloc(self.Omega, self.Phi_size * self.L         * sizeof(int) )
    24512992            self.W =       <unsigned int *> sage_realloc(self.W,     self.Phi_size * self.radix * 2 * sizeof(int) )
    2452             if not (self.w_gamma and self.alpha and self.Phi and self.Omega and self.W):
    2453                 if self.w_gamma: sage_free(self.w_gamma)
    2454                 if self.alpha: sage_free(self.alpha)
    2455                 if self.Phi: sage_free(self.Phi)
    2456                 if self.Omega: sage_free(self.Omega)
    2457                 if self.W: sage_free(self.W)
     2993            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:
     2994                if self.w_gamma is not NULL: sage_free(self.w_gamma)
     2995                if self.alpha is not NULL: sage_free(self.alpha)
     2996                if self.Phi is not NULL: sage_free(self.Phi)
     2997                if self.Omega is not NULL: sage_free(self.Omega)
     2998                if self.W is not NULL: sage_free(self.W)
    24582999                raise MemoryError("Memory.")
    24593000        word_gamma = self.w_gamma
    24603001        alpha = self.alpha # think of alpha as of length exactly nwords + ncols
    … … cdef class BinaryCodeClassifier: 
    25503091            elif state == 2: # Move down the search tree one level by refining nu:
    25513092                             # split out a vertex, and refine nu against it
    25523093                k += 1
    2553                 # TODO: Consider removing the following lines
    2554                 if k >= 2*self.radix: raise RuntimeError(\
    2555                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    25563094                nu.clear(k)
    25573095
    25583096                alpha[0] = nu.split_vertex(v[k-1], k)
    … … cdef class BinaryCodeClassifier: 
    26283166                        k = hzb__h_rho
    26293167                    else:
    26303168                        k = hh-1
    2631                 # TODO: Consider removing the following lines
    2632                 if k >= 2*self.radix: raise RuntimeError(\
    2633                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    26343169                # TODO: is the following line necessary?
    26353170                if k == -1: k = 0
    26363171               
    … … cdef class BinaryCodeClassifier: 
    26583193                # equivalent to zeta, i.e. there is no automorphism to discover.
    26593194                if k < hzf__h_zeta: state = 8; continue
    26603195
    2661                 nu.get_permutation(zeta, word_gamma, col_gamma, ham_wts)
     3196                nu.get_permutation(zeta, word_gamma, col_gamma)
    26623197#                print "gamma:", str([[word_gamma[i] for i from 0 <= i < nwords], [col_gamma[i] for i from 0 <= i < ncols]]).replace(' ','')
    26633198#                print Theta
    26643199                # if C^gamma == C, the permutation is an automorphism, goto 10
    … … cdef class BinaryCodeClassifier: 
    26843219                if j < 0: state = 6; continue
    26853220
    26863221                # if C(nu) == C(rho), get the automorphism and goto 10
    2687                 rho.get_permutation(nu, word_gamma, col_gamma, ham_wts)
     3222                rho.get_permutation(nu, word_gamma, col_gamma)
    26883223#                print "gamma:", str([[word_gamma[i] for i from 0 <= i < nwords], [col_gamma[i] for i from 0 <= i < ncols]]).replace(' ','')
    26893224#                print Theta
    26903225                state = 10
    … … cdef class BinaryCodeClassifier: 
    27663301
    27673302                # Otherwise, proceed to where zeta meets nu:
    27683303                k = h
    2769                 # TODO: Consider removing the following lines
    2770                 if k >= 2*self.radix: raise RuntimeError(\
    2771                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    27723304                state = 13
    27733305
    27743306            elif state == 11: # We have just found a new automorphism, and deduced that there may
    27753307                # be a better canonical label below the current branch off of zeta. So go to where
    27763308                # nu meets rho
    27773309                k = hb
    2778                 # TODO: Consider removing the following lines
    2779                 if k >= 2*self.radix: raise RuntimeError(\
    2780                     "A counterexample to an assumption the author made while writing this software has been encountered.")
    27813310                state = 12
    27823311
    27833312            elif state == 12: # Coming here from either state 6 or 11, the algorithm has discovered
    … … cdef class BinaryCodeClassifier: 
    29683497                ht = k # nodes descended from zeta[ht] are all equivalent
    29693498                hzf__h_zeta = k # max such that indicators for zeta and nu agree
    29703499                zeta = PartitionStack(nu)
    2971                 zeta.find_basis(ham_wts)
     3500                for i from 0 <= i < k:
     3501                    self.base[i] = v[i]
     3502                self.base_size = k
     3503                if k != self.radix:
     3504                    self.base[k] = -1
    29723505                # (POINT B)
    29733506                k -= 1
    29743507                rho = PartitionStack(nu)
    … … cdef class BinaryCodeClassifier: 
    29833516        rho.find_basis(ham_wts)
    29843517        for i from 0 <= i < ncols:
    29853518            self.labeling[rho.col_ents[i]] = i
    2986         for i from 0 <= i < nrows:
     3519        for i from 0 <= i < 2*nrows:
    29873520            self.labeling[i+ncols] = rho.basis_locations[i]
    29883521
     3522    def generate_children(self, BinaryCode B, int n, int d=2):
     3523        """
     3524        Use canonical augmentation to generate children of the code B.
     3525       
     3526        INPUT:
     3527        B -- a BinaryCode
     3528        n -- limit on the degree of the code
     3529        d -- test whether new vector has weight divisible by d. If d==4, this
     3530            ensures that all doubly-even canonically augmented children are
     3531            generated.
     3532       
     3533        EXAMPLE:
     3534            sage: from sage.coding.binary_code import *
     3535            sage: BC = BinaryCodeClassifier()
     3536            sage: B = BinaryCode(Matrix(GF(2), [[1,1,1,1]]))
     3537            sage: BC.generate_children(B, 6, 4)
     3538            [[1 1 1 1 0 0]
     3539            [0 1 0 1 1 1]]
     3540       
     3541        """
     3542        cdef BinaryCode m
     3543        cdef codeword *ortho_basis, *B_can_lab, *m_can_lab, current, swap
     3544        cdef codeword word, temp, gate, nonzero_gate, orbit, bwd, k_gate
     3545        cdef codeword *temp_basis, swap, *orbit_checks, orb_chx_size, orb_chx_shift, radix_gate
     3546        cdef WordPermutation *gwp, *hwp, *can_lab, *can_lab_inv
     3547        cdef WordPermutation **parent_generators
     3548        cdef BinaryCode B_aug
     3549        cdef int i, ii, j, jj, ij, k = 0, parity, combo, num_gens
     3550        cdef int base_size, *multimod2_index, row
     3551        cdef int *ham_wts = self.ham_wts
     3552        cdef int *num_inner_gens, *num_outer_gens, *v, log_2_radix
     3553        cdef bint bingo, bingo2, bingo3
    29893554
     3555        B.put_in_std_form()
     3556        ortho_basis = expand_to_ortho_basis(B, n) # modifies B!
    29903557
     3558#        print 'parent:'
     3559#        print B
     3560        aut_gp_gens, labeling, size, base = self._aut_gp_and_can_label(B)
     3561        B_can_lab = <codeword *> sage_malloc(B.nrows * sizeof(codeword))
     3562        can_lab = create_word_perm(labeling[:B.ncols])
     3563        if B_can_lab is NULL or m_can_lab is NULL or can_lab is NULL:
     3564            sage_free(ortho_basis)
     3565            if B_can_lab is not NULL:
     3566                sage_free(B_can_lab)
     3567            if can_lab is not NULL:
     3568                sage_free(can_lab)
     3569            raise MemoryError()
     3570        for i from 0 <= i < B.nrows:
     3571            B_can_lab[i] = permute_word_by_wp(can_lab, B.basis[i])
     3572        dealloc_word_perm(can_lab)
     3573        row = 0
     3574        current = 1
     3575        while row < B.nrows:
     3576            i = row
     3577            while i < B.nrows and not B_can_lab[i] & current:
     3578                i += 1
     3579            if i < B.nrows:
     3580                if i != row:
     3581                    swap = B_can_lab[row]
     3582                    B_can_lab[row] = B_can_lab[i]
     3583                    B_can_lab[i] = swap
     3584                for j from 0 <= j < row:
     3585                    if B_can_lab[j] & current:
     3586                        B_can_lab[j] ^= B_can_lab[row]
     3587                for j from row < j < B.nrows:
     3588                    if B_can_lab[j] & current:
     3589                        B_can_lab[j] ^= B_can_lab[row]
     3590                row += 1
     3591            current = current << 1
     3592        num_gens = len(aut_gp_gens)
     3593        base_size = len(base)
    29913594
     3595#        print 'gens:'
     3596#        for g in aut_gp_gens:
     3597#            print g
     3598
     3599        parent_generators = <WordPermutation **> sage_malloc( len(aut_gp_gens) * sizeof(WordPermutation*) )
     3600        temp_basis = <codeword *> sage_malloc( self.radix * sizeof(codeword) )
     3601
     3602        output = []
     3603
     3604       
     3605        for i from 0 <= i < len(aut_gp_gens):
     3606            parent_generators[i] = create_word_perm(aut_gp_gens[i] + range(B.ncols, n))
     3607
     3608        word = 0
     3609        while ortho_basis[k] & (((<codeword>1) << B.ncols) - 1):
     3610            k += 1
     3611        j = k
     3612        while ortho_basis[j]:
     3613            word ^= ortho_basis[j]
     3614            j += 1
     3615       
     3616#        print "ortho_basis:"
     3617#        for i from 0 <= i < k:
     3618#            print ''.join(reversed(Integer(ortho_basis[i]).binary().zfill(n)))
     3619#        print '-'
     3620#        for i from k <= i < j:
     3621#            print ''.join(reversed(Integer(ortho_basis[i]).binary().zfill(n)))
     3622#        print 'word:'
     3623#        print ''.join(reversed(Integer(word).binary().zfill(n)))
     3624       
     3625        log_2_radix = 0
     3626        while ((<codeword>1) << log_2_radix) < self.radix:
     3627            log_2_radix += 1
     3628        # now we assume (<codeword>1 << log_2_radix) == self.radix
     3629        if k < log_2_radix:
     3630            orb_chx_size = 0
     3631        else:
     3632            orb_chx_size = k - log_2_radix
     3633        orbit_checks = <codeword *> sage_malloc( ((<codeword>1) << orb_chx_size) * sizeof(codeword) )
     3634        if orbit_checks is NULL:
     3635            raise MemoryError()
     3636        for temp from 0 <= temp < ((<codeword>1) << orb_chx_size):
     3637            orbit_checks[temp] = 0
     3638       
     3639       
     3640        combo = 0
     3641        parity = 0
     3642        gate = (<codeword>1 << B.nrows) - 1
     3643        k_gate = (<codeword>1 << k) - 1
     3644        nonzero_gate = ( (<codeword>1 << (n-B.ncols)) - 1 ) << B.ncols
     3645        radix_gate = (((<codeword>1) << log_2_radix) - 1)
     3646#        print 'gate:', ''.join(reversed(Integer(gate).binary().zfill(n)))
     3647#        print 'gate:', ''.join(reversed(Integer(nonzero_gate).binary().zfill(n)))
     3648        while 1:
     3649#            print '    while 1'
     3650#            print '    ' + ''.join(reversed(Integer(word).binary().zfill(n)))
     3651            if nonzero_gate & word == nonzero_gate and \
     3652              (ham_wts[word & 65535] + ham_wts[(word >> 16) & 65535])%d == 0:
     3653#                print ''.join(reversed(Integer(word).binary().zfill(n)))
     3654                temp = (word >> B.nrows) & ((<codeword>1 << k) - 1)
     3655#                print "if not orbit_checks[temp >> log_2_radix] & ((<codeword>1) << (temp & radix_gate)):"
     3656#                print temp >> log_2_radix
     3657#                print temp & radix_gate
     3658                if not orbit_checks[temp >> log_2_radix] & ((<codeword>1) << (temp & radix_gate)):
     3659                    B_aug = BinaryCode(B, word)
     3660#                    print 'child:'
     3661#                    print B_aug
     3662#                    print 'canonically labeling child'
     3663                    aug_aut_gp_gens, aug_labeling, aug_size, aug_base = self._aut_gp_and_can_label(B_aug)
     3664#                    print 'done canonically labeling child'
     3665                    # check if (B, B_aug) ~ (m(B_aug), B_aug)
     3666
     3667                    can_lab = create_word_perm(aug_labeling[:n])
     3668#                    print 'relabeling:'
     3669#                    print [self.labeling[j] for j from 0 <= j < n]
     3670                    can_lab_inv = create_inv_word_perm(can_lab)
     3671                    for j from 0 <= j < B_aug.nrows:
     3672                        temp_basis[j] = permute_word_by_wp(can_lab, B_aug.basis[j])
     3673#                    print 'temp_basis:'
     3674#                    for j from 0 <= j < B_aug.nrows:
     3675#                        print ''.join(reversed(Integer(temp_basis[j]).binary().zfill(n)))
     3676
     3677                    # row reduce to get canonical label
     3678                    i = 0
     3679                    j = 0
     3680                    while j < B_aug.nrows:
     3681                        ii = j
     3682                        while ii < B_aug.nrows and not temp_basis[ii] & (<codeword>1 << i):
     3683                            ii += 1
     3684                        if ii != B_aug.nrows:
     3685                            if ii != j:
     3686                                swap = temp_basis[ii]
     3687                                temp_basis[ii] = temp_basis[j]
     3688                                temp_basis[j] = swap
     3689                            for jj from 0 <= jj < j:
     3690                                if temp_basis[jj] & (<codeword>1 << i):
     3691                                    temp_basis[jj] ^= temp_basis[j]
     3692                            for jj from j < jj < B_aug.nrows:
     3693                                if temp_basis[jj] & (<codeword>1 << i):
     3694                                    temp_basis[jj] ^= temp_basis[j]
     3695                            j += 1
     3696                        i += 1
     3697                    # done row reduction
     3698                           
     3699#                    print 'temp_basis:'
     3700                    for j from 0 <= j < B.nrows:
     3701                        temp_basis[j] = permute_word_by_wp(can_lab_inv, temp_basis[j])
     3702#                        print ''.join(reversed(Integer(temp_basis[j]).binary().zfill(n)))
     3703                    from sage.matrix.constructor import matrix
     3704                    from sage.rings.all import ZZ
     3705                    from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroupElement
     3706                    from sage.interfaces.gap import gap
     3707                    rs = []
     3708                    for i from 0 <= i < B.nrows:
     3709                        r = []
     3710                        for j from 0 <= j < n:
     3711                            r.append((((<codeword>1)<<j)&temp_basis[i])>>j)
     3712                        rs.append(r)
     3713                    m = BinaryCode(matrix(ZZ, rs))
     3714#                    print 'm:'
     3715#                    print m
     3716                    m_aut_gp_gens, m_labeling, m_size, m_base = self._aut_gp_and_can_label(m)
     3717                    from sage.rings.arith import factorial
     3718                    if True:#size*factorial(n-B.ncols) == m_size:
     3719#                        print 'in if'
     3720#                        print 'm_aut_gp_gens:', m_aut_gp_gens
     3721                        if len(m_aut_gp_gens) == 0:
     3722                            aut_m = PermutationGroup([()])
     3723                        else:
     3724                            aut_m = PermutationGroup([PermutationGroupElement([a+1 for a in g]) for g in m_aut_gp_gens])
     3725#                        print 'aut_m:', aut_m
     3726#                        print 'aug_aut_gp_gens:', aug_aut_gp_gens
     3727                        if len(aug_aut_gp_gens) == 0:
     3728                            aut_B_aug = PermutationGroup([()])
     3729                        else:
     3730                            aut_B_aug = PermutationGroup([PermutationGroupElement([a+1 for a in g]) for g in aug_aut_gp_gens])
     3731#                        print 'aut_B_aug:', aut_B_aug
     3732                        H = aut_m.__interface[gap].Intersection2(aut_B_aug.__interface[gap])
     3733#                        print 'H:', H
     3734                        rt_transversal = list(gap('List(RightTransversal( %s,%s ));'\
     3735                          %(str(aut_B_aug.__interface[gap]),str(H))))
     3736#                        print 'rt_transversal:', rt_transversal
     3737                        rt_transversal = [PermutationGroupElement(g) for g in rt_transversal if str(g) != '()']
     3738                        rt_transversal = [[a-1 for a in g.list()] for g in rt_transversal]
     3739                        rt_transversal = [g + range(len(g), n) for g in rt_transversal]
     3740                        rt_transversal.append(range(n))
     3741#                        print 'rt_transversal:', rt_transversal
     3742                        bingo2 = 0
     3743                        for coset_rep in rt_transversal:
     3744#                            print 'coset_rep:'
     3745#                            print coset_rep
     3746                            hwp = create_word_perm(coset_rep)
     3747                            #hwp = create_inv_word_perm(gwp) # since we want a left transversal
     3748                            #dealloc_word_perm(gwp)
     3749                            bingo2 = 1
     3750                            for j from 0 <= j < B.nrows:
     3751                                temp = permute_word_by_wp(hwp, temp_basis[j])
     3752                                if temp != B.words[temp & gate]:
     3753                                    bingo2 = 0
     3754                                    dealloc_word_perm(hwp)
     3755                                    break
     3756                            if bingo2:
     3757                                dealloc_word_perm(hwp)
     3758                                break
     3759                        if bingo2:
     3760                            from sage.matrix.constructor import Matrix
     3761                            from sage.rings.all import GF
     3762                            M = Matrix(GF(2), B_aug.nrows, B_aug.ncols)
     3763                            for i from 0 <= i < B_aug.ncols:
     3764                                for j from 0 <= j < B_aug.nrows:
     3765                                    M[j,i] = B_aug.is_one(1 << j, i)
     3766                            output.append(M)
     3767#                            print "ACCEPT"
     3768                    dealloc_word_perm(can_lab)
     3769                    dealloc_word_perm(can_lab_inv)
     3770                #...
     3771#                    print '    orbit_checks:'
     3772#                    for temp from 0 <= temp < ((<codeword>1) << orb_chx_size):
     3773#                        print '    ' + ''.join(reversed(Integer(orbit_checks[temp]).binary().zfill(n)))
     3774                    orbits = [word]
     3775                    j = 0
     3776                    while j < len(orbits):
     3777                        for i from 0 <= i < len(aut_gp_gens):
     3778#                            print '        i', i
     3779                            temp = <codeword> orbits[j]
     3780                            temp = permute_word_by_wp(parent_generators[i], temp)
     3781#                            print '        temp:', ''.join(reversed(Integer(temp).binary().zfill(n)))
     3782                            temp ^= B.words[temp & gate]
     3783#                            print '        temp:', ''.join(reversed(Integer(temp).binary().zfill(n)))
     3784                            if temp not in orbits:
     3785                                orbits.append(temp)
     3786                        j += 1
     3787                    for temp in orbits:
     3788                        temp = (temp >> B.nrows) & k_gate
     3789#                        print '        temp:', temp
     3790#                        print '        ', temp >> log_2_radix
     3791#                        print '        ', ((<codeword>1) << (temp & radix_gate))                   
     3792                        orbit_checks[temp >> log_2_radix] |= ((<codeword>1) << (temp & radix_gate))
     3793
     3794
     3795            parity ^= 1
     3796            i = 0
     3797            if not parity:
     3798                while not combo & (1 << i): i += 1
     3799                i += 1
     3800            if i == k: break
     3801            else:
     3802                combo ^= (1 << i)
     3803                word ^= ortho_basis[i]
     3804
     3805        for i from 0 <= i < len(aut_gp_gens):
     3806            dealloc_word_perm(parent_generators[i])
     3807        sage_free(B_can_lab)
     3808        sage_free(parent_generators)
     3809        sage_free(orbit_checks)
     3810        return output
     3811
     3812
     3813

  • sage/coding/linear_code.py 

    diff -r 23d0e6b6a675 -r d8cf41c5a6c5 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): 
    13811482        EXAMPLES:
    13821483            sage: C = HammingCode(3,GF(2))
    13831484            sage: G = C.automorphism_group_binary_code(); G
    1384             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)]
     1485            Permutation Group with generators [(3,4)(5,6), (3,5)(4,6), (2,3)(5,7), (1,2)(5,6)]
    13851486            sage: g = G("(2,3)(5,7)")
    13861487            sage: Cg = C.permuted_code(g)
    13871488            sage: Cg