"""
Fast binary code routines.
Some computations with linear binary codes. Fix a basis for $GF(2)^n$.
A linear binary code is a linear subspace of $GF(2)^n$, together with
this choice of basis. A permutation $g \in S_n$ of the fixed basis
gives rise to a permutation of the vectors, or words, in $GF(2)^n$,
sending $(w_i)$ to $(w_{g(i)})$. The permutation automorphism group of
the code $C$ is the set of permutations of the basis that bijectively
map $C$ to itself. Note that if $g$ is such a permutation, then
$$g(a_i) + g(b_i) = (a_{g(i)} + b_{g(i)}) = g((a_i) + (b_i)).$$
Over other fields, it is also required that the map be linear, which
as per above boils down to scalar multiplication. However, over
$GF(2),$ the only scalars are 0 and 1, so the linearity condition has
trivial effect.
AUTHOR:
Robert L Miller (Oct-Nov 2007)
* compiled code data structure
* union-find based orbit partition
* optimized partition stack class
* NICE-based partition refinement algorithm
* canonical generation function
"""
include 'sage/ext/cdefs.pxi'
from cpython.mem cimport *
include 'sage/ext/stdsage.pxi'
include 'sage/ext/interrupt.pxi'
from sage.structure.element import is_Matrix
from sage.misc.misc import cputime
from sage.rings.integer cimport Integer
from copy import copy
WORD_SIZE = sizeof(codeword) << 3
cdef enum:
chunk_size = 8
cdef inline int min(int a, int b):
if a > b:
return b
else:
return a
cdef int *hamming_weights():
cdef int *ham_wts, i
ham_wts = <int *> sage_malloc( 65536 * sizeof(int) )
if ham_wts is NULL:
sage_free(ham_wts)
raise MemoryError("Memory.")
ham_wts[0] = 0
ham_wts[1] = 1
ham_wts[2] = 1
ham_wts[3] = 2
for i from 4 <= i < 16:
ham_wts[i] = ham_wts[i & 3] + ham_wts[(i>>2) & 3]
for i from 16 <= i < 256:
ham_wts[i] = ham_wts[i & 15] + ham_wts[(i>>4) & 15]
for i from 256 <= i < 65536:
ham_wts[i] = ham_wts[i & 255] + ham_wts[(i>>8) & 255]
return ham_wts
include 'sage/misc/bitset.pxi'
def weight_dist(M):
"""
Computes the weight distribution of the row space of M.
EXAMPLES:
sage: from sage.coding.binary_code import weight_dist
sage: M = Matrix(GF(2),[
... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0],
... [0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0],
... [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1],
... [0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],
... [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]])
sage: weight_dist(M)
[1, 0, 0, 0, 0, 0, 0, 0, 30, 0, 0, 0, 0, 0, 0, 0, 1]
sage: M = Matrix(GF(2),[
... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0],
... [0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0],
... [0,0,0,0,0,1,0,1,0,0,0,1,1,1,1,1,1],
... [0,0,0,1,1,0,0,0,0,1,1,0,1,1,0,1,1]])
sage: weight_dist(M)
[1, 0, 0, 0, 0, 0, 0, 0, 11, 0, 0, 0, 4, 0, 0, 0, 0, 0]
sage: M=Matrix(GF(2),[
... [1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0,0],
... [0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0],
... [0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0],
... [0,0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0],
... [0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0,0],
... [0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0],
... [0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0],
... [0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1]])
sage: weight_dist(M)
[1, 0, 0, 0, 0, 0, 68, 0, 85, 0, 68, 0, 34, 0, 0, 0, 0, 0]
"""
cdef bitset_t word
cdef int i,j,k, dim=M.nrows(), deg=M.ncols()
cdef list L
cdef int *LL = <int *> sage_malloc((deg+1) * sizeof(int))
cdef bitset_s *basis = <bitset_s *> sage_malloc(dim * sizeof(bitset_s))
for i from 0 <= i < dim:
bitset_init(&basis[i], deg)
bitset_zero(&basis[i])
for j in M.row(i).nonzero_positions():
bitset_set(&basis[i], j)
for i from 0 <= i < deg+1: LL[i] = 0
bitset_init(word, deg)
bitset_zero(word)
i = 0
j = 0
while 1:
LL[bitset_hamming_weight(word)] += 1
i ^= 1
k = 0
if not i:
while not j & (1 << k): k += 1
k += 1
if k == dim: break
else:
j ^= (1 << k)
bitset_xor(word, word, &basis[k])
bitset_free(word)
L = [int(LL[i]) for i from 0 <= i < deg+1]
for i from 0 <= i < dim:
bitset_free(&basis[i])
sage_free(LL)
sage_free(basis)
return L
def test_word_perms(t_limit=5.0):
"""
Tests the WordPermutation structs for at least t_limit seconds.
These are structures written in pure C for speed, and are tested from this
function, which performs the following tests:
1. Tests create_word_perm, which creates a WordPermutation from a Python
list L representing a permutation i --> L[i]. Takes a random word and
permutes it by a random list permutation, and tests that the result
agrees with doing it the slow way.
1b. Tests create_array_word_perm, which creates a WordPermutation from a
C array. Does the same as above.
2. Tests create_comp_word_perm, which creates a WordPermutation as a
composition of two WordPermutations. Takes a random word and
two random permutations, and tests that the result of permuting by the
composition is correct.
3. Tests create_inv_word_perm and create_id_word_perm, which create a
WordPermutation as the inverse and identity permutations, resp.
Takes a random word and a random permutation, and tests that the result
permuting by the permutation and its inverse in either order, and
permuting by the identity both return the original word.
.. NOTE::
The functions permute_word_by_wp and dealloc_word_perm are implicitly
involved in each of the above tests.
TESTS::
sage: from sage.coding.binary_code import test_word_perms
sage: test_word_perms() # long time (5s on sage.math, 2011)
"""
cdef WordPermutation *g, *h, *i
cdef codeword cw1, cw2, cw3
cdef int n = sizeof(codeword) << 3
cdef int j, *arr = <int*> sage_malloc(n * sizeof(int))
if arr is NULL:
raise MemoryError("Error allocating memory.")
from sage.misc.prandom import randint
from sage.combinat.permutation import Permutations
S = Permutations(range(n))
t = cputime()
while cputime(t) < t_limit:
word = [randint(0,1) for _ in xrange(n)]
cw1 = 0
for j from 0 <= j < n:
cw1 += (<codeword>word[j]) << (<codeword>j)
gg = S.random_element()
g = create_word_perm(gg)
word2 = [0]*n
for j from 0 <= j < n:
word2[gg[j]] = word[j]
cw2 = permute_word_by_wp(g, cw1)
cw3 = 0
for j from 0 <= j < n:
cw3 += (<codeword>word2[j]) << (<codeword>j)
if cw3 != cw2:
print "ERROR1"
dealloc_word_perm(g)
gg = S.random_element()
for j from 0 <= j < n:
arr[j] = gg[j]
g = create_array_word_perm(arr, 0, n)
word2 = [0]*n
for j from 0 <= j < n:
word2[gg[j]] = word[j]
cw2 = permute_word_by_wp(g, cw1)
cw3 = 0
for j from 0 <= j < n:
cw3 += (<codeword>word2[j]) << (<codeword>j)
if cw3 != cw2:
print "ERROR1b"
dealloc_word_perm(g)
gg = S.random_element()
hh = S.random_element()
g = create_word_perm(gg)
h = create_word_perm(hh)
i = create_comp_word_perm(g, h)
word2 = [0]*n
for j from 0 <= j < n:
word2[gg[hh[j]]] = word[j]
cw2 = permute_word_by_wp(i, cw1)
cw3 = 0
for j from 0 <= j < n:
cw3 += (<codeword>word2[j]) << (<codeword>j)
if cw3 != cw2:
print "ERROR2"
dealloc_word_perm(g)
dealloc_word_perm(h)
dealloc_word_perm(i)
gg = S.random_element()
g = create_word_perm(gg)
h = create_inv_word_perm(g)
i = create_id_word_perm(n)
cw2 = permute_word_by_wp(g, cw1)
cw2 = permute_word_by_wp(h, cw2)
if cw1 != cw2:
print "ERROR3a"
cw2 = permute_word_by_wp(h, cw1)
cw2 = permute_word_by_wp(g, cw2)
if cw1 != cw2:
print "ERROR3b"
cw2 = permute_word_by_wp(i, cw1)
if cw1 != cw2:
print "ERROR3c"
dealloc_word_perm(g)
dealloc_word_perm(h)
dealloc_word_perm(i)
sage_free(arr)
cdef WordPermutation *create_word_perm(object list_perm):
r"""
Create a word permutation from a Python list permutation L, i.e. such that
$i \mapsto L[i]$.
"""
cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
cdef codeword *images_i, image
cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
if word_perm is NULL:
raise RuntimeError("Error allocating memory.")
word_perm.degree = len(list_perm)
list_perm = copy(list_perm)
while num_chunks*chunk_size < word_perm.degree:
num_chunks += 1
word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
if word_perm.images is NULL:
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.chunk_num = num_chunks
words_per_chunk = 1 << chunk_size
word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
list_perm += range(len(list_perm), chunk_size*num_chunks)
word_perm.chunk_words = words_per_chunk
for i from 0 <= i < num_chunks:
images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
if images_i is NULL:
for j from 0 <= j < i:
sage_free(word_perm.images[j])
sage_free(word_perm.images)
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.images[i] = images_i
for j from 0 <= j < chunk_size:
images_i[1 << j] = (<codeword>1) << list_perm[chunk_size*i + j]
image = <codeword> 0
parity = 0
comb = 0
while 1:
images_i[comb] = image
parity ^= 1
j = 0
if not parity:
while not comb & (1 << j): j += 1
j += 1
if j == chunk_size: break
else:
comb ^= (1 << j)
image ^= images_i[1 << j]
return word_perm
cdef WordPermutation *create_array_word_perm(int *array, int start, int degree):
"""
Create a word permutation of a given degree from a C array, starting at start.
"""
cdef int i, j, cslim, parity, comb, words_per_chunk, num_chunks = 1
cdef codeword *images_i, image
cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
if word_perm is NULL:
raise RuntimeError("Error allocating memory.")
word_perm.degree = degree
while num_chunks*chunk_size < word_perm.degree:
num_chunks += 1
word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
if word_perm.images is NULL:
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.chunk_num = num_chunks
words_per_chunk = 1 << chunk_size
word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
word_perm.chunk_words = words_per_chunk
for i from 0 <= i < num_chunks:
images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
if images_i is NULL:
for j from 0 <= j < i:
sage_free(word_perm.images[j])
sage_free(word_perm.images)
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.images[i] = images_i
cslim = min(chunk_size, degree - i*chunk_size)
for j from 0 <= j < cslim:
images_i[1 << j] = (<codeword>1) << array[start + chunk_size*i + j]
image = <codeword> 0
parity = 0
comb = 0
while 1:
images_i[comb] = image
parity ^= 1
j = 0
if not parity:
while not comb & (1 << j): j += 1
j += 1
if j == chunk_size: break
else:
comb ^= (1 << j)
image ^= images_i[1 << j]
return word_perm
cdef WordPermutation *create_id_word_perm(int degree):
"""
Create the identity word permutation of degree degree.
"""
cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
cdef codeword *images_i, image
cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
if word_perm is NULL:
raise RuntimeError("Error allocating memory.")
word_perm.degree = degree
while num_chunks*chunk_size < degree:
num_chunks += 1
word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
if word_perm.images is NULL:
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.chunk_num = num_chunks
words_per_chunk = 1 << chunk_size
word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
word_perm.chunk_words = words_per_chunk
for i from 0 <= i < num_chunks:
images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
if images_i is NULL:
for j from 0 <= j < i:
sage_free(word_perm.images[j])
sage_free(word_perm.images)
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.images[i] = images_i
for j from 0 <= j < chunk_size:
images_i[1 << j] = (<codeword>1) << (chunk_size*i + j)
image = <codeword> 0
parity = 0
comb = 0
while 1:
images_i[comb] = image
parity ^= 1
j = 0
if not parity:
while not comb & (1 << j): j += 1
j += 1
if j == chunk_size: break
else:
comb ^= (1 << j)
image ^= images_i[1 << j]
return word_perm
cdef WordPermutation *create_comp_word_perm(WordPermutation *g, WordPermutation *h):
r"""
Create the composition of word permutations $g \circ h$.
"""
cdef int i, j, parity, comb, words_per_chunk, num_chunks = 1
cdef codeword *images_i, image
cdef WordPermutation *word_perm = <WordPermutation *> sage_malloc( sizeof(WordPermutation) )
if word_perm is NULL:
raise RuntimeError("Error allocating memory.")
word_perm.degree = g.degree
while num_chunks*chunk_size < word_perm.degree:
num_chunks += 1
word_perm.images = <codeword **> sage_malloc(num_chunks * sizeof(codeword *))
if word_perm.images is NULL:
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.chunk_num = num_chunks
words_per_chunk = 1 << chunk_size
word_perm.gate = ( (<codeword>1) << chunk_size ) - 1
word_perm.chunk_words = words_per_chunk
for i from 0 <= i < num_chunks:
images_i = <codeword *> sage_malloc(words_per_chunk * sizeof(codeword))
if images_i is NULL:
for j from 0 <= j < i:
sage_free(word_perm.images[j])
sage_free(word_perm.images)
sage_free(word_perm)
raise RuntimeError("Error allocating memory.")
word_perm.images[i] = images_i
for j from 0 <= j < chunk_size:
image = (<codeword>1) << (chunk_size*i + j)
image = permute_word_by_wp(h, image)
image = permute_word_by_wp(g, image)
images_i[1 << j] = image
image = <codeword> 0
parity = 0
comb = 0
while 1:
images_i[comb] = image
parity ^= 1
j = 0
if not parity:
while not comb & (1 << j): j += 1
j += 1
if j == chunk_size: break
else:
comb ^= (1 << j)
image ^= images_i[1 << j]
return word_perm
cdef WordPermutation *create_inv_word_perm(WordPermutation *g):
r"""
Create the inverse $g^{-1}$ of the word permutation of $g$.
"""
cdef int i, j, *array = <int *> sage_malloc( g.degree * sizeof(int) )
cdef codeword temp
cdef WordPermutation *w
for i from 0 <= i < g.degree:
j = 0
temp = permute_word_by_wp(g, (<codeword>1) << i)
while not ((<codeword>1) << j) & temp:
j += 1
array[j] = i
w = create_array_word_perm(array, 0, g.degree)
sage_free(array)
return w
cdef int dealloc_word_perm(WordPermutation *wp):
"""
Free the memory used by a word permutation.
"""
cdef int i
for i from 0 <= i < wp.chunk_num:
sage_free(wp.images[i])
sage_free(wp.images)
sage_free(wp)
cdef codeword permute_word_by_wp(WordPermutation *wp, codeword word):
"""
Return the codeword obtained by applying the permutation wp to word.
"""
cdef int num_chunks = wp.chunk_num
cdef int i
cdef codeword gate = wp.gate
cdef codeword image = 0
cdef codeword **images = wp.images
for i from 0 <= i < num_chunks:
image += images[i][(word >> i*chunk_size) & gate]
return image
def test_expand_to_ortho_basis(B=None):
"""
This function is written in pure C for speed, and is tested from this
function.
INPUT:
B -- a BinaryCode in standard form
OUTPUT:
An array of codewords which represent the expansion of a basis for $B$ to a
basis for $(B^\prime)^\perp$, where $B^\prime = B$ if the all-ones vector 1
is in $B$, otherwise $B^\prime = \text{span}(B,1)$ (note that this guarantees
that all the vectors in the span of the output have even weight).
TESTS:
sage: from sage.coding.binary_code import test_expand_to_ortho_basis, BinaryCode
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]])
sage: B = BinaryCode(M)
sage: B.put_in_std_form()
0
sage: test_expand_to_ortho_basis(B=B)
INPUT CODE:
Binary [10,2] linear code, generator matrix
[1010001111]
[0101111111]
Expanding to the basis of an orthogonal complement...
Basis:
0010000010
0001000010
0000100001
0000010001
0000001001
"""
cdef codeword *output
cdef int k=0, i
cdef BinaryCode C
if not isinstance(B, BinaryCode):
raise TypeError()
C = B
print "INPUT CODE:"
print C
print "Expanding to the basis of an orthogonal complement..."
output = expand_to_ortho_basis(C, C.ncols)
print "Basis:"
while output[k]:
k += 1
for i from 0 <= i < k:
print ''.join(reversed(Integer(output[i]).binary().zfill(C.ncols)))
sage_free(output)
cdef codeword *expand_to_ortho_basis(BinaryCode B, int n):
r"""
INPUT:
B -- a BinaryCode in standard form
n -- the degree
OUTPUT:
An array of codewords which represent the expansion of a basis for $B$ to a
basis for $(B^\prime)^\perp$, where $B^\prime = B$ if the all-ones vector 1
is in $B$, otherwise $B^\prime = \text{span}(B,1)$ (note that this guarantees
that all the vectors in the span of the output have even weight).
"""
cdef codeword *basis, word = 0, temp, new, pivots = 0, combo, parity
cdef codeword n_gate = (~<codeword>0) >> ( (sizeof(codeword)<<3) - n)
cdef int i, j, m, k = B.nrows, dead, d
cdef WordPermutation *wp
basis = <codeword *> sage_malloc( (n+1) * sizeof(codeword) )
if basis is NULL:
raise MemoryError()
for i from 0 <= i < k:
basis[i] = B.basis[i]
word ^= basis[i]
word = (~word) & n_gate
if word:
basis[k] = word
temp = (<codeword>1) << k
i = k
while not word & temp:
temp = temp << 1
i += 1
for j from 0 <= j < k:
if temp & basis[j]:
basis[j] ^= word
temp += (<codeword>1 << k) - 1
i = k
word = <codeword>1 << k
k += 1
else:
temp = (<codeword>1 << k) - 1
i = k
word = <codeword>1 << k
j = k
while i < n:
new = 0
for m from 0 <= m < k:
if basis[m] & word:
new ^= basis[m]
basis[j] = (new & temp) + word
j += ((word ^ temp) >> i) & 1
i += 1
word = word << 1
temp = (<codeword>1 << B.nrows) - 1
for i from k <= i < n:
basis[i-k] = basis[i] ^ B.words[basis[i] & temp]
k = n-k
i = 0
word = (<codeword>1 << B.nrows)
while i < k and (word & n_gate):
m = i
while m < k and not basis[m] & word:
m += 1
if m < k:
pivots += word
if m != i:
new = basis[i]
basis[i] = basis[m]
basis[m] = new
for j from 0 <= j < i:
if basis[j] & word:
basis[j] ^= basis[i]
for j from i < j < k:
if basis[j] & word:
basis[j] ^= basis[i]
i += 1
word = word << 1
for j from i <= j < n:
basis[j] = 0
perm = range(B.nrows)
perm_c = []
for j from B.nrows <= j < B.ncols:
if (<codeword>1 << j) & pivots:
perm.append(j)
else:
perm_c.append(j)
perm.extend(perm_c)
perm.extend(range(B.ncols, n))
perm_c = [0]*n
for j from 0 <= j < n:
perm_c[perm[j]] = j
wp = create_word_perm(perm_c)
for j from 0 <= j < i:
basis[j] = permute_word_by_wp(wp, basis[j])
for j from 0 <= j < B.nrows:
B.basis[j] = permute_word_by_wp(wp, B.basis[j])
dealloc_word_perm(wp)
word = 0
parity = 0
combo = 0
while 1:
B.words[combo] = word
parity ^= 1
j = 0
if not parity:
while not combo & (1 << j): j += 1
j += 1
if j == B.nrows: break
else:
combo ^= (1 << j)
word ^= B.basis[j]
return basis
cdef class BinaryCode:
"""
Minimal, but optimized, binary code object.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1]])
sage: B = BinaryCode(M) # create from matrix
sage: C = BinaryCode(B, 60) # create using glue
sage: D = BinaryCode(C, 240)
sage: E = BinaryCode(D, 85)
sage: B
Binary [4,1] linear code, generator matrix
[1111]
sage: C
Binary [6,2] linear code, generator matrix
[111100]
[001111]
sage: D
Binary [8,3] linear code, generator matrix
[11110000]
[00111100]
[00001111]
sage: E
Binary [8,4] linear code, generator matrix
[11110000]
[00111100]
[00001111]
[10101010]
sage: M = Matrix(GF(2), [[1]*32])
sage: B = BinaryCode(M)
sage: B
Binary [32,1] linear code, generator matrix
[11111111111111111111111111111111]
"""
def __cinit__(self, arg1, arg2=None):
cdef int nrows, i, j, size
cdef int nwords, other_nwords, parity, combination
cdef codeword word, glue_word
cdef BinaryCode other
cdef codeword *self_words, *self_basis, *other_basis
self.radix = sizeof(int) << 3
if is_Matrix(arg1):
self.ncols = arg1.ncols()
self.nrows = arg1.nrows()
nrows = self.nrows
self.nwords = 1 << nrows
nwords = self.nwords
elif isinstance(arg1, BinaryCode):
other = <BinaryCode> arg1
self.nrows = other.nrows + 1
glue_word = <codeword> arg2
size = 0
while 0 < ((<codeword> 1) << size) <= glue_word:
size += 1
if other.ncols > size:
self.ncols = other.ncols
else:
self.ncols = size
other_nwords = other.nwords
self.nwords = 2 * other_nwords
nrows = self.nrows
nwords = self.nwords
else: raise NotImplementedError("!")
if self.nrows >= self.radix or self.ncols > self.radix:
raise NotImplementedError("Columns and rows are stored as ints. This code is too big.")
self.words = <codeword *> sage_malloc( nwords * sizeof(int) )
self.basis = <codeword *> sage_malloc( nrows * sizeof(int) )
if self.words is NULL or self.basis is NULL:
if self.words is not NULL: sage_free(self.words)
if self.basis is not NULL: sage_free(self.basis)
raise MemoryError("Memory.")
self_words = self.words
self_basis = self.basis
if is_Matrix(arg1):
rows = arg1.rows()
for i from 0 <= i < nrows:
word = <codeword> 0
for j in rows[i].nonzero_positions():
word += (1<<j)
self_basis[i] = word
word = <codeword> 0
parity = 0
combination = 0
while 1:
self_words[combination] = word
parity ^= 1
j = 0
if not parity:
while not combination & (1 << j): j += 1
j += 1
if j == nrows: break
else:
combination ^= (1 << j)
word ^= self_basis[j]
else:
other_basis = other.basis
for i from 0 <= i < nrows-1:
self_basis[i] = other_basis[i]
i = nrows - 1
self_basis[i] = glue_word
memcpy(self_words, other.words, other_nwords*(self.radix>>3))
for combination from 0 <= combination < other_nwords:
self_words[combination+other_nwords] = self_words[combination] ^ glue_word
def __dealloc__(self):
sage_free(self.words)
sage_free(self.basis)
def __reduce__(self):
"""
Method for pickling and unpickling BinaryCodes.
TESTS:
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1]])
sage: B = BinaryCode(M)
sage: loads(dumps(B)) == B
True
"""
return BinaryCode, (self.matrix(),)
def __cmp__(self, other):
"""
Comparison of BinaryCodes.
TESTS:
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1]])
sage: B = BinaryCode(M)
sage: C = BinaryCode(B.matrix())
sage: B == C
True
"""
return cmp(self.matrix(), other.matrix())
def matrix(self):
"""
Returns the generator matrix of the BinaryCode, i.e. the code is the
rowspace of B.matrix().
EXAMPLE:
sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
sage: from sage.coding.binary_code import *
sage: B = BinaryCode(M)
sage: B.matrix()
[1 1 1 1 0 0]
[0 0 1 1 1 1]
"""
cdef int i, j
from sage.matrix.constructor import matrix
from sage.rings.all import GF
rows = []
for i from 0 <= i < self.nrows:
row = [0]*self.ncols
for j from 0 <= j < self.ncols:
if self.basis[i] & ((<codeword>1) << j):
row[j] = 1
rows.append(row)
return matrix(GF(2), self.nrows, self.ncols, rows)
def print_data(self):
"""
Print all data for self.
EXAMPLES:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1]])
sage: B = BinaryCode(M)
sage: C = BinaryCode(B, 60)
sage: D = BinaryCode(C, 240)
sage: E = BinaryCode(D, 85)
sage: B.print_data() # random - actually "print P.print_data()"
ncols: 4
nrows: 1
nwords: 2
radix: 32
basis:
1111
words:
0000
1111
sage: C.print_data() # random - actually "print P.print_data()"
ncols: 6
nrows: 2
nwords: 4
radix: 32
basis:
111100
001111
words:
000000
111100
001111
110011
sage: D.print_data() # random - actually "print P.print_data()"
ncols: 8
nrows: 3
nwords: 8
radix: 32
basis:
11110000
00111100
00001111
words:
00000000
11110000
00111100
11001100
00001111
11111111
00110011
11000011
sage: E.print_data() # random - actually "print P.print_data()"
ncols: 8
nrows: 4
nwords: 16
radix: 32
basis:
11110000
00111100
00001111
10101010
words:
00000000
11110000
00111100
11001100
00001111
11111111
00110011
11000011
10101010
01011010
10010110
01100110
10100101
01010101
10011001
01101001
"""
from sage.graphs.generic_graph_pyx import binary
cdef int ui
cdef int i
s = ''
s += "ncols:" + str(self.ncols)
s += "\nnrows:" + str(self.nrows)
s += "\nnwords:" + str(self.nwords)
s += "\nradix:" + str(self.radix)
s += "\nbasis:\n"
for i from 0 <= i < self.nrows:
b = list(binary(self.basis[i]).zfill(self.ncols))
b.reverse()
b.append('\n')
s += ''.join(b)
s += "\nwords:\n"
for ui from 0 <= ui < self.nwords:
b = list(binary(self.words[ui]).zfill(self.ncols))
b.reverse()
b.append('\n')
s += ''.join(b)
def __repr__(self):
"""
String representation of self.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]])
sage: B = BinaryCode(M)
sage: B
Binary [8,4] linear code, generator matrix
[11110000]
[00111100]
[00001111]
[10101010]
"""
cdef int i, j
s = 'Binary [%d,%d] linear code, generator matrix\n'%(self.ncols, self.nrows)
for i from 0 <= i < self.nrows:
s += '[' + self._word((<codeword> 1)<<i) + ']\n'
return s
def _word(self, coords):
"""
Considering coords as an integer in binary, think of the 0's and 1's as
coefficients of the basis given by self.matrix(). This function returns
a string representation of that word.
EXAMPLE:
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1]])
sage: B = BinaryCode(M)
sage: B._word(0)
'0000'
sage: B._word(1)
'1111'
Note that behavior under input which does not represent a word in
the code is unspecified (gives nonsense).
"""
s = ''
for j from 0 <= j < self.ncols:
s += '%d'%self.is_one(coords,j)
return s
def _is_one(self, word, col):
"""
Returns the col-th letter of word, i.e. 0 or 1. Words are expressed
as integers, which represent linear combinations of the rows of the
generator matrix of the code.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]])
sage: B = BinaryCode(M)
sage: B
Binary [8,4] linear code, generator matrix
[11110000]
[00111100]
[00001111]
[10101010]
sage: B._is_one(7, 4)
0
sage: B._is_one(8, 4)
1
sage: B._is_automorphism([1,0,3,2,4,5,6,7], [0, 1, 2, 3, 4, 5, 6, 7, 9, 8, 11, 10, 13, 12, 15, 14])
1
"""
return self.is_one(word, col) != 0
cdef int is_one(self, int word, int column):
return (self.words[word] & (<codeword> 1 << column)) >> column
def _is_automorphism(self, col_gamma, word_gamma):
"""
Check whether a given permutation is an automorphism of the code.
INPUT:
col_gamma -- permutation sending i |--> col_gamma[i] acting
on the columns.
word_gamma -- permutation sending i |--> word_gamma[i] acting
on the words.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]])
sage: B = BinaryCode(M)
sage: B
Binary [8,4] linear code, generator matrix
[11110000]
[00111100]
[00001111]
[10101010]
sage: B._is_automorphism([1,0,3,2,4,5,6,7], [0, 1, 2, 3, 4, 5, 6, 7, 9, 8, 11, 10, 13, 12, 15, 14])
1
"""
cdef int i
cdef int *_col_gamma
cdef int *_word_gamma
_word_gamma = <int *> sage_malloc(self.nwords * sizeof(int))
_col_gamma = <int *> sage_malloc(self.ncols * sizeof(int))
if _col_gamma is NULL or _word_gamma is NULL:
if _word_gamma is not NULL: sage_free(_word_gamma)
if _col_gamma is not NULL: sage_free(_col_gamma)
raise MemoryError("Memory.")
for i from 0 <= i < self.nwords:
_word_gamma[i] = word_gamma[i]
for i from 0 <= i < self.ncols:
_col_gamma[i] = col_gamma[i]
result = self.is_automorphism(_col_gamma, _word_gamma)
sage_free(_col_gamma)
sage_free(_word_gamma)
return result
cdef int is_automorphism(self, int *col_gamma, int *word_gamma):
cdef int i, j, self_nwords = self.nwords, self_ncols = self.ncols
i = 1
while i < self_nwords:
for j from 0 <= j < self_ncols:
if self.is_one(i, j) != self.is_one(word_gamma[i], col_gamma[j]):
return 0
i = i << 1
return 1
def apply_permutation(self, labeling):
"""
Apply a column permutation to the code.
INPUT:
labeling -- a list permutation of the columns
EXAMPLE:
sage: from sage.coding.binary_code import *
sage: B = BinaryCode(codes.ExtendedBinaryGolayCode().gen_mat())
sage: B
Binary [24,12] linear code, generator matrix
[100000000000101011100011]
[010000000000111110010010]
[001000000000110100101011]
[000100000000110001110110]
[000010000000110011011001]
[000001000000011001101101]
[000000100000001100110111]
[000000010000101101111000]
[000000001000010110111100]
[000000000100001011011110]
[000000000010101110001101]
[000000000001010111000111]
sage: B.apply_permutation(range(11,-1,-1) + range(12, 24))
sage: B
Binary [24,12] linear code, generator matrix
[000000000001101011100011]
[000000000010111110010010]
[000000000100110100101011]
[000000001000110001110110]
[000000010000110011011001]
[000000100000011001101101]
[000001000000001100110111]
[000010000000101101111000]
[000100000000010110111100]
[001000000000001011011110]
[010000000000101110001101]
[100000000000010111000111]
"""
self._apply_permutation_to_basis(labeling)
self._update_words_from_basis()
cdef void _apply_permutation_to_basis(self, object labeling):
cdef WordPermutation *wp
cdef int i
wp = create_word_perm(labeling)
for i from 0 <= i < self.nrows:
self.basis[i] = permute_word_by_wp(wp, self.basis[i])
dealloc_word_perm(wp)
cdef void _update_words_from_basis(self):
cdef codeword word
cdef int j, parity, combination
word = 0
parity = 0
combination = 0
while 1:
self.words[combination] = word
parity ^= 1
j = 0
if not parity:
while not combination & (1 << j): j += 1
j += 1
if j == self.nrows: break
else:
combination ^= (1 << j)
word ^= self.basis[j]
cpdef int put_in_std_form(self):
"""
Put the code in binary form, which is defined by an identity matrix on
the left, augmented by a matrix of data.
EXAMPLE:
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
sage: B = BinaryCode(M); B
Binary [6,2] linear code, generator matrix
[111100]
[001111]
sage: B.put_in_std_form(); B
0
Binary [6,2] linear code, generator matrix
[101011]
[010111]
"""
cdef codeword swap, current = 1, pivots = 0
cdef int i, j, k, row = 0
cdef object perm
while row < self.nrows:
i = row
while i < self.nrows and not self.basis[i] & current:
i += 1
if i < self.nrows:
pivots += current
if i != row:
swap = self.basis[row]
self.basis[row] = self.basis[i]
self.basis[i] = swap
for j from 0 <= j < row:
if self.basis[j] & current:
self.basis[j] ^= self.basis[row]
for j from row < j < self.nrows:
if self.basis[j] & current:
self.basis[j] ^= self.basis[row]
row += 1
current = current << 1
perm = [0]*self.ncols
j = 0
k = self.nrows
for i from 0 <= i < self.ncols:
if ((<codeword>1) << i) & pivots:
perm[i] = j
j += 1
else:
perm[i] = k
k += 1
self._apply_permutation_to_basis(perm)
self._update_words_from_basis()
cdef class OrbitPartition:
"""
Structure which keeps track of which vertices are equivalent
under the part of the automorphism group that has already been
seen, during search. Essentially a disjoint-set data structure*,
which also keeps track of the minimum element and size of each
cell of the partition, and the size of the partition.
* http://en.wikipedia.org/wiki/Disjoint-set_data_structure
"""
def __cinit__(self, int nrows, int ncols):
cdef int col
cdef int nwords, word
nwords = (1 << nrows)
self.nwords = nwords
self.ncols = ncols
self.wd_parent = <int *> sage_malloc( nwords * sizeof(int) )
self.wd_rank = <int *> sage_malloc( nwords * sizeof(int) )
self.wd_min_cell_rep = <int *> sage_malloc( nwords * sizeof(int) )
self.wd_size = <int *> sage_malloc( nwords * sizeof(int) )
self.col_parent = <int *> sage_malloc( ncols * sizeof(int) )
self.col_rank = <int *> sage_malloc( ncols * sizeof(int) )
self.col_min_cell_rep = <int *> sage_malloc( ncols * sizeof(int) )
self.col_size = <int *> sage_malloc( ncols * sizeof(int) )
if self.wd_parent is NULL or self.wd_rank is NULL or self.wd_min_cell_rep is NULL \
or self.wd_size is NULL or self.col_parent is NULL or self.col_rank is NULL \
or self.col_min_cell_rep is NULL or self.col_size is NULL:
if self.wd_parent is not NULL: sage_free(self.wd_parent)
if self.wd_rank is not NULL: sage_free(self.wd_rank)
if self.wd_min_cell_rep is not NULL: sage_free(self.wd_min_cell_rep)
if self.wd_size is not NULL: sage_free(self.wd_size)
if self.col_parent is not NULL: sage_free(self.col_parent)
if self.col_rank is not NULL: sage_free(self.col_rank)
if self.col_min_cell_rep is not NULL: sage_free(self.col_min_cell_rep)
if self.col_size is not NULL: sage_free(self.col_size)
raise MemoryError("Memory.")
for word from 0 <= word < nwords:
self.wd_parent[word] = word
self.wd_rank[word] = 0
self.wd_min_cell_rep[word] = word
self.wd_size[word] = 1
for col from 0 <= col < ncols:
self.col_parent[col] = col
self.col_rank[col] = 0
self.col_min_cell_rep[col] = col
self.col_size[col] = 1
def __dealloc__(self):
sage_free(self.wd_parent)
sage_free(self.wd_rank)
sage_free(self.wd_min_cell_rep)
sage_free(self.wd_size)
sage_free(self.col_parent)
sage_free(self.col_rank)
sage_free(self.col_min_cell_rep)
sage_free(self.col_size)
def __repr__(self):
"""
Return a string representation of the orbit partition.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: O = OrbitPartition(4, 8)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15
Columns:
0,1,2,3,4,5,6,7
"""
cdef int i
cdef int j
s = 'OrbitPartition on %d words and %d columns. Data:\n'%(self.nwords, self.ncols)
s += 'Words:\n'
for i from 0 <= i < self.nwords:
s += '%d,'%self.wd_parent[i]
s = s[:-1] + '\nColumns:\n'
for j from 0 <= j < self.ncols:
s += '%d,'%self.col_parent[j]
return s[:-1]
def _wd_find(self, word):
"""
Returns the root of word.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: O = OrbitPartition(4, 8)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15
Columns:
0,1,2,3,4,5,6,7
sage: O._wd_find(12)
12
"""
return self.wd_find(word)
cdef int wd_find(self, int word):
if self.wd_parent[word] == word:
return word
else:
self.wd_parent[word] = self.wd_find(self.wd_parent[word])
return self.wd_parent[word]
def _wd_union(self, x, y):
"""
Join the cells containing x and y.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: O = OrbitPartition(4, 8)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15
Columns:
0,1,2,3,4,5,6,7
sage: O._wd_union(1,10)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,1,11,12,13,14,15
Columns:
0,1,2,3,4,5,6,7
sage: O._wd_find(10)
1
"""
self.wd_union(x, y)
cdef void wd_union(self, int x, int y):
cdef int x_root, y_root
x_root = self.wd_find(x)
y_root = self.wd_find(y)
if self.wd_rank[x_root] > self.wd_rank[y_root]:
self.wd_parent[y_root] = x_root
self.wd_min_cell_rep[x_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root])
self.wd_size[x_root] += self.wd_size[y_root]
elif self.wd_rank[x_root] < self.wd_rank[y_root]:
self.wd_parent[x_root] = y_root
self.wd_min_cell_rep[y_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root])
self.wd_size[y_root] += self.wd_size[x_root]
elif x_root != y_root:
self.wd_parent[y_root] = x_root
self.wd_min_cell_rep[x_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root])
self.wd_size[x_root] += self.wd_size[y_root]
self.wd_rank[x_root] += 1
def _col_find(self, col):
"""
Returns the root of col.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: O = OrbitPartition(4, 8)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15
Columns:
0,1,2,3,4,5,6,7
sage: O._col_find(6)
6
"""
return self.col_find(col)
cdef int col_find(self, int col):
if self.col_parent[col] == col:
return col
else:
self.col_parent[col] = self.col_find(self.col_parent[col])
return self.col_parent[col]
def _col_union(self, x, y):
"""
Join the cells containing x and y.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: O = OrbitPartition(4, 8)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15
Columns:
0,1,2,3,4,5,6,7
sage: O._col_union(1,4)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15
Columns:
0,1,2,3,1,5,6,7
sage: O._col_find(4)
1
"""
self.col_union(x, y)
cdef void col_union(self, int x, int y):
cdef int x_root, y_root
x_root = self.col_find(x)
y_root = self.col_find(y)
if self.col_rank[x_root] > self.col_rank[y_root]:
self.col_parent[y_root] = x_root
self.col_min_cell_rep[x_root] = min(self.col_min_cell_rep[x_root],self.col_min_cell_rep[y_root])
self.col_size[x_root] += self.col_size[y_root]
elif self.col_rank[x_root] < self.col_rank[y_root]:
self.col_parent[x_root] = y_root
self.col_min_cell_rep[y_root] = min(self.col_min_cell_rep[x_root],self.col_min_cell_rep[y_root])
self.col_size[y_root] += self.col_size[x_root]
elif x_root != y_root:
self.col_parent[y_root] = x_root
self.col_min_cell_rep[x_root] = min(self.col_min_cell_rep[x_root],self.col_min_cell_rep[y_root])
self.col_size[x_root] += self.col_size[y_root]
self.col_rank[x_root] += 1
def _merge_perm(self, col_gamma, wd_gamma):
"""
Merges the cells of self under the given permutation. If gamma[a] = b,
then after merge_perm, a and b will be in the same cell. Returns 0 if
nothing was done, otherwise returns 1.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: O = OrbitPartition(4, 8)
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15
Columns:
0,1,2,3,4,5,6,7
sage: O._merge_perm([1,0,3,2,4,5,6,7], [0,1,2,3,4,5,6,7,9,8,11,10,13,12,15,14])
1
sage: O
OrbitPartition on 16 words and 8 columns. Data:
Words:
0,1,2,3,4,5,6,7,8,8,10,10,12,12,14,14
Columns:
0,0,2,2,4,5,6,7
"""
cdef int i
cdef int *_col_gamma
cdef int *_wd_gamma
_wd_gamma = <int *> sage_malloc(self.nwords * sizeof(int))
_col_gamma = <int *> sage_malloc(self.ncols * sizeof(int))
if _col_gamma is NULL or _wd_gamma is NULL:
if _wd_gamma is not NULL: sage_free(_wd_gamma)
if _col_gamma is not NULL: sage_free(_col_gamma)
raise MemoryError("Memory.")
for i from 0 <= i < self.nwords:
_wd_gamma[i] = wd_gamma[i]
for i from 0 <= i < self.ncols:
_col_gamma[i] = col_gamma[i]
result = self.merge_perm(_col_gamma, _wd_gamma)
sage_free(_col_gamma)
sage_free(_wd_gamma)
return result
cdef int merge_perm(self, int *col_gamma, int *wd_gamma):
cdef int i, gamma_i_root
cdef int j, gamma_j_root, return_value = 0
cdef int *self_wd_parent = self.wd_parent
cdef int *self_col_parent = self.col_parent
for i from 0 <= i < self.nwords:
gamma_i_root = self.wd_find(wd_gamma[i])
if gamma_i_root != i:
return_value = 1
self.wd_union(i, gamma_i_root)
for j from 0 <= j < self.ncols:
gamma_j_root = self.col_find(col_gamma[j])
if gamma_j_root != j:
return_value = 1
self.col_union(j, gamma_j_root)
return return_value
cdef class PartitionStack:
"""
Partition stack structure for traversing the search tree during automorphism
group computation.
"""
def __cinit__(self, arg1, arg2=None):
cdef int k, nwords, ncols, sizeof_int
cdef PartitionStack other = None
cdef int *wd_ents, *wd_lvls, *col_ents, *col_lvls
cdef int *col_degs, *col_counts, *col_output, *wd_degs, *wd_counts, *wd_output
sizeof_int = sizeof(int)
try:
self.nrows = <int> arg1
self.nwords = 1 << self.nrows
self.ncols = <int> arg2
except Exception:
other = arg1
self.nrows = other.nrows
self.nwords = other.nwords
self.ncols = other.ncols
self.radix = sizeof_int << 3
self.flag = (1 << (self.radix-1))
self.wd_ents = <int *> sage_malloc( self.nwords * sizeof_int )
self.wd_lvls = <int *> sage_malloc( self.nwords * sizeof_int )
self.col_ents = <int *> sage_malloc( self.ncols * sizeof_int )
self.col_lvls = <int *> sage_malloc( self.ncols * sizeof_int )
self.col_degs = <int *> sage_malloc( self.ncols * sizeof_int )
self.col_counts = <int *> sage_malloc( self.nwords * sizeof_int )
self.col_output = <int *> sage_malloc( self.ncols * sizeof_int )
self.wd_degs = <int *> sage_malloc( self.nwords * sizeof_int )
self.wd_counts = <int *> sage_malloc( (self.ncols+1) * sizeof_int )
self.wd_output = <int *> sage_malloc( self.nwords * sizeof_int )
if self.wd_ents is NULL or self.wd_lvls is NULL or self.col_ents is NULL \
or self.col_lvls is NULL or self.col_degs is NULL or self.col_counts is NULL \
or self.col_output is NULL or self.wd_degs is NULL or self.wd_counts is NULL \
or self.wd_output is NULL:
if self.wd_ents is not NULL: sage_free(self.wd_ents)
if self.wd_lvls is not NULL: sage_free(self.wd_lvls)
if self.col_ents is not NULL: sage_free(self.col_ents)
if self.col_lvls is not NULL: sage_free(self.col_lvls)
if self.col_degs is not NULL: sage_free(self.col_degs)
if self.col_counts is not NULL: sage_free(self.col_counts)
if self.col_output is not NULL: sage_free(self.col_output)
if self.wd_degs is not NULL: sage_free(self.wd_degs)
if self.wd_counts is not NULL: sage_free(self.wd_counts)
if self.wd_output is not NULL: sage_free(self.wd_output)
raise MemoryError("Memory.")
nwords = self.nwords
ncols = self.ncols
if other:
memcpy(self.wd_ents, other.wd_ents, self.nwords * sizeof_int)
memcpy(self.wd_lvls, other.wd_lvls, self.nwords * sizeof_int)
memcpy(self.col_ents, other.col_ents, self.ncols * sizeof_int)
memcpy(self.col_lvls, other.col_lvls, self.ncols * sizeof_int)
else:
wd_ents = self.wd_ents
wd_lvls = self.wd_lvls
col_ents = self.col_ents
col_lvls = self.col_lvls
for k from 0 <= k < nwords-1:
wd_ents[k] = k
wd_lvls[k] = 2*ncols
for k from 0 <= k < ncols-1:
col_ents[k] = k
col_lvls[k] = 2*ncols
wd_ents[nwords-1] = nwords-1
wd_lvls[nwords-1] = -1
col_ents[ncols-1] = ncols-1
col_lvls[ncols-1] = -1
col_degs = self.col_degs
col_counts = self.col_counts
col_output = self.col_output
wd_degs = self.wd_degs
wd_counts = self.wd_counts
wd_output = self.wd_output
for k from 0 <= k < ncols:
col_degs[k]=0
col_output[k]=0
wd_counts[k]=0
wd_counts[ncols]=0
for k from 0 <= k < nwords:
col_counts[k]=0
wd_degs[k]=0
wd_output[k]=0
def __dealloc__(self):
if self.basis_locations: sage_free(self.basis_locations)
sage_free(self.wd_ents)
sage_free(self.wd_lvls)
sage_free(self.col_ents)
sage_free(self.col_lvls)
sage_free(self.col_degs)
sage_free(self.col_counts)
sage_free(self.col_output)
sage_free(self.wd_degs)
sage_free(self.wd_counts)
sage_free(self.wd_output)
def print_data(self):
"""
Prints all data for self.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: print P.print_data()
nwords:4
nrows:2
ncols:6
radix:32
wd_ents:
0
1
2
3
wd_lvls:
12
12
12
-1
col_ents:
0
1
2
3
4
5
col_lvls:
12
12
12
12
12
-1
col_degs:
0
0
0
0
0
0
col_counts:
0
0
0
0
col_output:
0
0
0
0
0
0
wd_degs:
0
0
0
0
wd_counts:
0
0
0
0
0
0
0
wd_output:
0
0
0
0
"""
cdef int i, j
s = ''
s += "nwords:" + str(self.nwords) + '\n'
s += "nrows:" + str(self.nrows) + '\n'
s += "ncols:" + str(self.ncols) + '\n'
s += "radix:" + str(self.radix) + '\n'
s += "wd_ents:" + '\n'
for i from 0 <= i < self.nwords:
s += str(self.wd_ents[i]) + '\n'
s += "wd_lvls:" + '\n'
for i from 0 <= i < self.nwords:
s += str(self.wd_lvls[i]) + '\n'
s += "col_ents:" + '\n'
for i from 0 <= i < self.ncols:
s += str(self.col_ents[i]) + '\n'
s += "col_lvls:" + '\n'
for i from 0 <= i < self.ncols:
s += str(self.col_lvls[i]) + '\n'
s += "col_degs:" + '\n'
for i from 0 <= i < self.ncols:
s += str(self.col_degs[i]) + '\n'
s += "col_counts:" + '\n'
for i from 0 <= i < self.nwords:
s += str(self.col_counts[i]) + '\n'
s += "col_output:" + '\n'
for i from 0 <= i < self.ncols:
s += str(self.col_output[i]) + '\n'
s += "wd_degs:" + '\n'
for i from 0 <= i < self.nwords:
s += str(self.wd_degs[i]) + '\n'
s += "wd_counts:" + '\n'
for i from 0 <= i < self.ncols + 1:
s += str(self.wd_counts[i]) + '\n'
s += "wd_output:" + '\n'
for i from 0 <= i < self.nwords:
s += str(self.wd_output[i]) + '\n'
if self.basis_locations:
s += "basis_locations:" + '\n'
j = 1
while (1 << j) < self.nwords:
j += 1
for i from 0 <= i < j:
s += str(self.basis_locations[i]) + '\n'
return s
def __repr__(self):
"""
Return a string representation of self.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
"""
cdef int i, j, k
s = ''
last = ''
current = ''
for k from 0 <= k < 2*self.ncols:
current = self._repr_at_k(k)
if current == last: break
s += current
last = current
return s
def _repr_at_k(self, k):
"""
Gives a string representing the partition at level k:
EXAMPLE:
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6); P
({0,1,2,3}) ({0,1,2,3,4,5})
sage: P._repr_at_k(0)
'({0,1,2,3}) ({0,1,2,3,4,5})\n'
"""
s = '({'
for j from 0 <= j < self.nwords:
s += str(self.wd_ents[j])
if self.wd_lvls[j] <= k:
s += '},{'
else:
s += ','
s = s[:-2] + ') '
s += '({'
for j from 0 <= j < self.ncols:
s += str(self.col_ents[j])
if self.col_lvls[j] <= k:
s += '},{'
else:
s += ','
s = s[:-2] + ')\n'
return s
def _is_discrete(self, k):
"""
Returns whether the partition at level k is discrete.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P._sort_wds(0, [0,2,3,1], 5)
0
sage: P
({0,3,1,2}) ({0,1,2,3,4,5})
({0,3,1,2}) ({0},{1,2,3,4,5})
({0,3,1,2}) ({0},{1},{2,3,4,5})
({0,3,1,2}) ({0},{1},{2},{3,4,5})
({0,3,1,2}) ({0},{1},{2},{3},{4,5})
({0},{3},{1},{2}) ({0},{1},{2},{3},{4},{5})
sage: P._is_discrete(4)
0
sage: P._is_discrete(5)
1
"""
return self.is_discrete(k)
cdef int is_discrete(self, int k):
cdef int i, self_ncols = self.ncols, self_nwords = self.nwords
cdef int *self_col_lvls = self.col_lvls
cdef int *self_wd_lvls = self.wd_lvls
for i from 0 <= i < self_ncols:
if self_col_lvls[i] > k:
return 0
for i from 0 <= i < self_nwords:
if self_wd_lvls[i] > k:
return 0
return 1
def _num_cells(self, k):
"""
Returns the number of cells in the partition at level k.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
({0,1,2,3}) ({0},{1,2,3,4,5})
({0,1,2,3}) ({0},{1},{2,3,4,5})
({0,1,2,3}) ({0},{1},{2},{3,4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4},{5})
sage: P._num_cells(3)
5
"""
return self.num_cells(k)
cdef int num_cells(self, int k):
cdef int i, j = 0
cdef int *self_wd_lvls = self.wd_lvls
cdef int *self_col_lvls = self.col_lvls
for i from 0 <= i < self.nwords:
if self_wd_lvls[i] <= k:
j += 1
for i from 0 <= i < self.ncols:
if self_col_lvls[i] <= k:
j += 1
return j
def _sat_225(self, k):
"""
Returns whether the partition at level k satisfies the hypotheses of
Lemma 2.25 in Brendan McKay's Practical Graph Isomorphism paper (see
sage/graphs/graph_isom.pyx.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P._sat_225(3)
0
sage: P._sat_225(4)
1
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
({0,1,2,3}) ({0},{1,2,3,4,5})
({0,1,2,3}) ({0},{1},{2,3,4,5})
({0,1,2,3}) ({0},{1},{2},{3,4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4},{5})
"""
return self.sat_225(k)
cdef int sat_225(self, int k):
cdef int i, n = self.nwords + self.ncols, in_cell = 0
cdef int nontrivial_cells = 0, total_cells = self.num_cells(k)
cdef int *self_wd_lvls = self.wd_lvls
cdef int *self_col_lvls = self.col_lvls
if n <= total_cells + 4:
return 1
for i from 0 <= i < self.nwords:
if self_wd_lvls[i] <= k:
if in_cell:
nontrivial_cells += 1
in_cell = 0
else:
in_cell = 1
in_cell = 0
for i from 0 <= i < self.ncols:
if self_col_lvls[i] <= k:
if in_cell:
nontrivial_cells += 1
in_cell = 0
else:
in_cell = 1
if n == total_cells + nontrivial_cells:
return 1
if n == total_cells + nontrivial_cells + 1:
return 1
return 0
cdef void new_min_cell_reps(self, int k, unsigned int *Omega, int start):
cdef int i, j
cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls
cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents
cdef int reps = (1 << self_col_ents[0]), length, word
cdef int radix = self.radix, nwords = self.nwords, ncols = self.ncols
length = 1 + nwords/radix
if nwords%radix:
length += 1
for i from 0 <= i < length:
Omega[start+i] = 0
for i from 0 < i < ncols:
Omega[start] += ((self_col_lvls[i-1] <= k) << self_col_ents[i])
Omega[start+1] = (1 << self_wd_ents[0])
for i from 0 < i < nwords:
if self_wd_lvls[i-1] <= k:
word = self_wd_lvls[i-1]
Omega[start+1+word/radix] += (1 << word%radix)
cdef void fixed_vertices(self, int k, unsigned int *Phi, unsigned int *Omega, int start):
cdef int i, j, length, ell, fixed = 0
cdef int radix = self.radix, nwords = self.nwords, ncols = self.ncols
cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls
cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents
for i from 0 <= i < ncols:
fixed += ((self_col_lvls[i] <= k) << self_col_ents[i])
Phi[start] = fixed & Omega[start]
length = 1 + nwords/self.radix
if nwords%self.radix:
length += 1
for i from 0 < i < length:
Phi[start+i] = 0
for i from 0 <= i < nwords:
ell = self_wd_ents[i]
Phi[start+1+ell/radix] = ((self_wd_lvls[i] <= k) << ell%radix)
for i from 0 < i < length:
Phi[i] &= Omega[i]
cdef int new_first_smallest_nontrivial(self, int k, unsigned int *W, int start):
cdef int ell
cdef int i = 0, j = 0, location = 0, min = self.ncols, nwords = self.nwords
cdef int min_is_col = 1, radix = self.radix
cdef int *self_col_lvls = self.col_lvls, *self_wd_lvls = self.wd_lvls
cdef int *self_col_ents = self.col_ents, *self_wd_ents = self.wd_ents
while 1:
if self_col_lvls[i] <= k:
if i != j and min > i - j + 1:
min = i - j + 1
location = j
j = i + 1
if self_col_lvls[i] == -1: break
i += 1
j = location
ell = 1 + nwords/radix
if nwords%radix:
ell += 1
for i from 0 <= i < ell:
W[start+i] = 0
if min_is_col:
while 1:
if self_col_lvls[j] <= k: break
j += 1
i = location
while i <= j:
W[start] ^= (1 << self_col_ents[i])
i += 1
return self_col_ents[location]
else:
while 1:
if self_wd_lvls[j] <= k: break
j += 1
i = location
while i <= j:
ell = self_wd_ents[i]
W[start+1+ell/radix] ^= (1 << ell%radix)
i += 1
return self_wd_ents[location] ^ self.flag
def _dangerous_dont_use_set_ents_lvls(self, col_ents, col_lvls, wd_ents, wd_lvls):
"""
For debugging only.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
sage: P._dangerous_dont_use_set_ents_lvls([99]*6, [0,3,2,3,5,-1], [4,3,5,6], [3,2,1,-1])
sage: P
({4,3,5,6}) ({99},{99,99,99,99,99})
({4,3,5},{6}) ({99},{99,99,99,99,99})
({4,3},{5},{6}) ({99},{99,99},{99,99,99})
({4},{3},{5},{6}) ({99},{99},{99},{99},{99,99})
"""
cdef int i
for i from 0 <= i < len(col_ents):
self.col_ents[i] = col_ents[i]
for i from 0 <= i < len(col_lvls):
self.col_lvls[i] = col_lvls[i]
for i from 0 <= i < len(wd_ents):
self.wd_ents[i] = wd_ents[i]
for i from 0 <= i < len(wd_lvls):
self.wd_lvls[i] = wd_lvls[i]
def _col_percolate(self, start, end):
"""
Do one round of bubble sort on ents.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: P._dangerous_dont_use_set_ents_lvls(range(5,-1,-1), [1,2,2,3,3,-1], range(3,-1,-1), [1,1,2,-1])
sage: P
({3,2,1,0}) ({5,4,3,2,1,0})
({3},{2},{1,0}) ({5},{4,3,2,1,0})
({3},{2},{1},{0}) ({5},{4},{3},{2,1,0})
({3},{2},{1},{0}) ({5},{4},{3},{2},{1},{0})
sage: P._wd_percolate(0,3)
sage: P._col_percolate(0,5)
sage: P
({0,3,2,1}) ({0,5,4,3,2,1})
({0},{3},{2,1}) ({0},{5,4,3,2,1})
({0},{3},{2},{1}) ({0},{5},{4},{3,2,1})
({0},{3},{2},{1}) ({0},{5},{4},{3},{2},{1})
"""
self.col_percolate(start, end)
cdef void col_percolate(self, int start, int end):
cdef int i, temp
cdef int *self_col_ents = self.col_ents
for i from end >= i > start:
if self_col_ents[i] < self_col_ents[i-1]:
temp = self.col_ents[i]
self_col_ents[i] = self_col_ents[i-1]
self_col_ents[i-1] = temp
def _wd_percolate(self, start, end):
"""
Do one round of bubble sort on ents.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: P._dangerous_dont_use_set_ents_lvls(range(5,-1,-1), [1,2,2,3,3,-1], range(3,-1,-1), [1,1,2,-1])
sage: P
({3,2,1,0}) ({5,4,3,2,1,0})
({3},{2},{1,0}) ({5},{4,3,2,1,0})
({3},{2},{1},{0}) ({5},{4},{3},{2,1,0})
({3},{2},{1},{0}) ({5},{4},{3},{2},{1},{0})
sage: P._wd_percolate(0,3)
sage: P._col_percolate(0,5)
sage: P
({0,3,2,1}) ({0,5,4,3,2,1})
({0},{3},{2,1}) ({0},{5,4,3,2,1})
({0},{3},{2},{1}) ({0},{5},{4},{3,2,1})
({0},{3},{2},{1}) ({0},{5},{4},{3},{2},{1})
"""
self.wd_percolate(start, end)
cdef void wd_percolate(self, int start, int end):
cdef int i, temp
cdef int *self_wd_ents = self.wd_ents
for i from end >= i > start:
if self_wd_ents[i] < self_wd_ents[i-1]:
temp = self.wd_ents[i]
self_wd_ents[i] = self_wd_ents[i-1]
self_wd_ents[i-1] = temp
def _split_vertex(self, v, k):
"""
Split vertex v out, placing it before the rest of the cell it was in.
Returns the location of the split vertex.
NOTE:
There is a convention regarding whether a vertex is a word or a
column. See the 'flag' attribute of the PartitionStack object:
If vertex&flag is not zero, it is a word.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
({0,1,2,3}) ({0},{1,2,3,4,5})
({0,1,2,3}) ({0},{1},{2,3,4,5})
({0,1,2,3}) ({0},{1},{2},{3,4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4},{5})
"""
return self.split_vertex(v, k)
cdef int split_vertex(self, int v, int k):
cdef int i = 0, j, flag = self.flag
cdef int *ents
cdef int *lvls
if v & flag:
ents = self.wd_ents
lvls = self.wd_lvls
v = v ^ flag
while ents[i] != v: i += 1
v = v ^ flag
else:
ents = self.col_ents
lvls = self.col_lvls
while ents[i] != v: i += 1
j = i
while lvls[i] > k: i += 1
if j == 0 or lvls[j-1] <= k:
if v & self.flag:
self.wd_percolate(j+1, i)
else:
self.col_percolate(j+1, i)
else:
while j != 0 and lvls[j-1] > k:
ents[j] = ents[j-1]
j -= 1
if v & flag:
ents[j] = v ^ flag
else:
ents[j] = v
lvls[j] = k
if v & flag:
return j ^ flag
else:
return j
def _col_degree(self, C, col, wd_ptr, k):
"""
Returns the number of words in the cell specified by wd_ptr that have a
1 in the col-th column.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
sage: B = BinaryCode(M)
sage: B
Binary [6,2] linear code, generator matrix
[111100]
[001111]
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
({0,1,2,3}) ({0},{1,2,3,4,5})
({0,1,2,3}) ({0},{1},{2,3,4,5})
({0,1,2,3}) ({0},{1},{2},{3,4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4},{5})
sage: P._col_degree(B, 2, 0, 2)
2
"""
return self.col_degree(C, col, wd_ptr, k)
cdef int col_degree(self, BinaryCode CG, int col, int wd_ptr, int k):
cdef int i = 0
cdef int *self_wd_lvls = self.wd_lvls, *self_wd_ents = self.wd_ents
while 1:
if CG.is_one(self_wd_ents[wd_ptr], col): i += 1
if self_wd_lvls[wd_ptr] > k: wd_ptr += 1
else: break
return i
def _wd_degree(self, C, wd, col_ptr, k):
"""
Returns the number of columns in the cell specified by col_ptr that are
1 in wd.
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])
sage: B = BinaryCode(M)
sage: B
Binary [6,2] linear code, generator matrix
[111100]
[001111]
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
({0,1,2,3}) ({0},{1,2,3,4,5})
({0,1,2,3}) ({0},{1},{2,3,4,5})
({0,1,2,3}) ({0},{1},{2},{3,4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4},{5})
sage: P._wd_degree(B, 1, 1, 1)
3
"""
cdef int *ham_wts = hamming_weights()
result = self.wd_degree(C, wd, col_ptr, k, ham_wts)
sage_free(ham_wts)
return result
cdef int wd_degree(self, BinaryCode CG, int wd, int col_ptr, int k, int *ham_wts):
cdef int *self_col_lvls = self.col_lvls, *self_col_ents = self.col_ents
cdef int mask = (1 << self_col_ents[col_ptr])
while self_col_lvls[col_ptr] > k:
col_ptr += 1
mask += (1 << self_col_ents[col_ptr])
mask &= CG.words[wd]
return ham_wts[mask & 65535] + ham_wts[(mask >> 16) & 65535]
def _sort_cols(self, start, degrees, k):
"""
Essentially a counting sort, but on only one cell of the partition.
INPUT:
start -- location of the beginning of the cell
k -- at what level of refinement the partition of interest lies
degrees -- the counts to sort by
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P._sort_cols(1, [0,2,2,1,1], 1)
2
sage: P
({0,1,2,3}) ({0,1,4,5,2,3})
({0,1,2,3}) ({0},{1},{4,5},{2,3})
"""
cdef int i
for i from 0 <= i < len(degrees):
self.col_degs[i] = degrees[i]
return self.sort_cols(start, k)
cdef int sort_cols(self, int start, int k):
cdef int i, j, max, max_location, self_ncols = self.ncols
cdef int self_nwords = self.nwords, ii
cdef int *self_col_counts = self.col_counts
cdef int *self_col_lvls = self.col_lvls
cdef int *self_col_degs = self.col_degs
cdef int *self_col_ents = self.col_ents
cdef int *self_col_output = self.col_output
for ii from 0 <= ii < self_nwords:
self_col_counts[ii] = 0
i = 0
while self_col_lvls[i+start] > k:
self_col_counts[self_col_degs[i]] += 1
i += 1
self_col_counts[self_col_degs[i]] += 1
max = self_col_counts[0]
max_location = 0
for ii from 0 < ii < self_nwords:
if self_col_counts[ii] > max:
max = self_col_counts[ii]
max_location = ii
self_col_counts[ii] += self_col_counts[ii-1]
for j from i >= j >= 0:
self_col_counts[self_col_degs[j]] -= 1
self_col_output[self_col_counts[self_col_degs[j]]] = self_col_ents[start+j]
max_location = self_col_counts[max_location] + start
for j from 0 <= j <= i:
self_col_ents[start+j] = self_col_output[j]
ii = 1
while ii < self_nwords and self_col_counts[ii] <= i:
if self_col_counts[ii] > 0:
self_col_lvls[start + self_col_counts[ii] - 1] = k
self.col_percolate(start + self_col_counts[ii-1], start + self_col_counts[ii] - 1)
ii += 1
return max_location
def _sort_wds(self, start, degrees, k):
"""
Essentially a counting sort, but on only one cell of the partition.
INPUT:
start -- location of the beginning of the cell
k -- at what level of refinement the partition of interest lies
degrees -- the counts to sort by
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(3, 6)
sage: P._sort_wds(0, [0,0,3,3,3,3,2,2], 1)
4
sage: P
({0,1,6,7,2,3,4,5}) ({0,1,2,3,4,5})
({0,1},{6,7},{2,3,4,5}) ({0,1,2,3,4,5})
"""
cdef int i
for i from 0 <= i < len(degrees):
self.wd_degs[i] = degrees[i]
return self.sort_wds(start, k)
cdef int sort_wds(self, int start, int k):
cdef int i, j, max, max_location, self_nwords = self.nwords
cdef int ii, self_ncols = self.ncols
cdef int *self_wd_counts = self.wd_counts
cdef int *self_wd_lvls = self.wd_lvls
cdef int *self_wd_degs = self.wd_degs
cdef int *self_wd_ents = self.wd_ents
cdef int *self_wd_output = self.wd_output
for ii from 0 <= ii < self_ncols+1:
self_wd_counts[ii] = 0
i = 0
while self_wd_lvls[i+start] > k:
self_wd_counts[self_wd_degs[i]] += 1
i += 1
self_wd_counts[self_wd_degs[i]] += 1
max = self_wd_counts[0]
max_location = 0
for ii from 0 < ii < self_ncols+1:
if self_wd_counts[ii] > max:
max = self_wd_counts[ii]
max_location = ii
self_wd_counts[ii] += self_wd_counts[ii-1]
for j from i >= j >= 0:
if j > i: break
self_wd_counts[self_wd_degs[j]] -= 1
self_wd_output[self_wd_counts[self_wd_degs[j]]] = self_wd_ents[start+j]
max_location = self_wd_counts[max_location] + start
for j from 0 <= j <= i:
self_wd_ents[start+j] = self_wd_output[j]
ii = 1
while ii < self_ncols+1 and self_wd_counts[ii] <= i:
if self_wd_counts[ii] > 0:
self_wd_lvls[start + self_wd_counts[ii] - 1] = k
self.wd_percolate(start + self_wd_counts[ii-1], start + self_wd_counts[ii] - 1)
ii += 1
return max_location
def _refine(self, k, alpha, CG):
"""
EXAMPLE::
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]])
sage: B = BinaryCode(M)
sage: P = PartitionStack(4, 8)
sage: P._refine(1, [[0,0],[1,0]], B)
181
sage: P._split_vertex(0, 2)
0
sage: P._refine(2, [[0,0]], B)
290
sage: P._split_vertex(1, 3)
1
sage: P._refine(3, [[0,1]], B)
463
sage: P._split_vertex(2, 4)
2
sage: P._refine(4, [[0,2]], B)
1500
sage: P._split_vertex(3, 5)
3
sage: P._refine(5, [[0,3]], B)
641
sage: P._split_vertex(4, 6)
4
sage: P._refine(6, [[0,4]], B)
1224
sage: P._is_discrete(5)
0
sage: P._is_discrete(6)
1
sage: P
({0,4,6,2,13,9,11,15,10,14,12,8,7,3,1,5}) ({0,1,2,3,4,7,6,5})
({0},{4,6,2,13,9,11,15,10,14,12,8,7,3,1},{5}) ({0,1,2,3,4,7,6,5})
({0},{4,6,2,13,9,11,15},{10,14,12,8,7,3,1},{5}) ({0},{1,2,3,4,7,6,5})
({0},{4,6,2},{13,9,11,15},{10,14,12,8},{7,3,1},{5}) ({0},{1},{2,3,4,7,6,5})
({0},{4},{6,2},{13,9},{11,15},{10,14},{12,8},{7,3},{1},{5}) ({0},{1},{2},{3,4,7,6,5})
({0},{4},{6,2},{13,9},{11,15},{10,14},{12,8},{7,3},{1},{5}) ({0},{1},{2},{3},{4,7,6,5})
({0},{4},{6},{2},{13},{9},{11},{15},{10},{14},{12},{8},{7},{3},{1},{5}) ({0},{1},{2},{3},{4},{7},{6},{5})
"""
cdef int i, alpha_length = len(alpha)
cdef int *_alpha = <int *> sage_malloc( (self.nwords + self.ncols) * sizeof(int) )
cdef int *ham_wts = hamming_weights()
if _alpha is NULL:
raise MemoryError("Memory.")
for i from 0 <= i < alpha_length:
if alpha[i][0]:
_alpha[i] = alpha[i][1] ^ self.flag
else:
_alpha[i] = alpha[i][1]
result = self.refine(k, _alpha, alpha_length, CG, ham_wts)
sage_free(_alpha)
sage_free(ham_wts)
return result
cdef int refine(self, int k, int *alpha, int alpha_length, BinaryCode CG, int *ham_wts):
cdef int q, r, s, t, flag = self.flag, self_ncols = self.ncols
cdef int t_w, self_nwords = self.nwords, invariant = 0, i, j, m = 0
cdef int *self_wd_degs = self.wd_degs, *self_wd_lvls = self.wd_lvls, *self_wd_ents = self.wd_ents
cdef int *self_col_degs = self.col_degs, *self_col_lvls = self.col_lvls, *self_col_ents = self.col_ents
while not self.is_discrete(k) and m < alpha_length:
invariant += 1
j = 0
if alpha[m] & flag:
while j < self_ncols:
i = j; s = 0
invariant += 8
while 1:
self_col_degs[i-j] = self.col_degree(CG, self_col_ents[i], alpha[m]^flag, k)
if s == 0 and self_col_degs[i-j] != self_col_degs[0]: s = 1
i += 1
if self_col_lvls[i-1] <= k: break
if s:
invariant += 8
t = self.sort_cols(j, k)
invariant += t
q = m
while q < alpha_length:
if alpha[q] == j:
alpha[q] = t
break
q += 1
r = j
while 1:
if r == j or self.col_lvls[r-1] == k:
if r != t:
alpha[alpha_length] = r
alpha_length += 1
r += 1
if r >= i: break
invariant += self.col_degree(CG, self_col_ents[i-1], alpha[m]^flag, k)
invariant += (i-j)
j = i
else:
while j < self.nwords:
i = j; s = 0
invariant += 64
while 1:
self_wd_degs[i-j] = self.wd_degree(CG, self_wd_ents[i], alpha[m], k, ham_wts)
if s == 0 and self_wd_degs[i-j] != self_wd_degs[0]: s = 1
i += 1
if self_wd_lvls[i-1] <= k: break
if s:
invariant += 64
t_w = self.sort_wds(j, k)
invariant += t_w
q = m
j ^= flag
while q < alpha_length:
if alpha[q] == j:
alpha[q] = t_w ^ flag
break
q += 1
j ^= flag
r = j
while 1:
if r == j or self.wd_lvls[r-1] == k:
if r != t_w:
alpha[alpha_length] = r^flag
alpha_length += 1
r += 1
if r >= i: break
invariant += self.wd_degree(CG, self_wd_ents[i-1], alpha[m], k, ham_wts)
invariant += (i-j)
j = i
m += 1
return invariant
def _clear(self, k):
"""
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(2, 6)
sage: [P._split_vertex(i,i+1) for i in range(5)]
[0, 1, 2, 3, 4]
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
({0,1,2,3}) ({0},{1,2,3,4,5})
({0,1,2,3}) ({0},{1},{2,3,4,5})
({0,1,2,3}) ({0},{1},{2},{3,4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4,5})
({0,1,2,3}) ({0},{1},{2},{3},{4},{5})
sage: P._clear(2)
sage: P
({0,1,2,3}) ({0,1,2,3,4,5})
({0,1,2,3}) ({0},{1,2,3,4,5})
"""
self.clear(k)
cdef void clear(self, int k):
cdef int i, j = 0, nwords = self.nwords, ncols = self.ncols
cdef int *wd_lvls = self.wd_lvls, *col_lvls = self.col_lvls
for i from 0 <= i < nwords:
if wd_lvls[i] >= k:
wd_lvls[i] += 1
if wd_lvls[i] < k:
self.wd_percolate(j, i)
j = i + 1
j = 0
for i from 0 <= i < ncols:
if col_lvls[i] >= k:
col_lvls[i] += 1
if col_lvls[i] < k:
self.col_percolate(j, i)
j = i + 1
def _cmp(self, other, C):
"""
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]])
sage: B = BinaryCode(M)
sage: P = PartitionStack(4, 8)
sage: P._refine(0, [[0,0],[1,0]], B)
181
sage: P._split_vertex(0, 1)
0
sage: P._refine(1, [[0,0]], B)
290
sage: P._split_vertex(1, 2)
1
sage: P._refine(2, [[0,1]], B)
463
sage: P._split_vertex(2, 3)
2
sage: P._refine(3, [[0,2]], B)
1500
sage: P._split_vertex(4, 4)
4
sage: P._refine(4, [[0,4]], B)
1224
sage: P._is_discrete(4)
1
sage: Q = PartitionStack(P)
sage: Q._clear(4)
sage: Q._split_vertex(5, 4)
4
sage: Q._refine(4, [[0,4]], B)
1224
sage: Q._is_discrete(4)
1
sage: Q._cmp(P, B)
0
"""
return self.cmp(other, C)
cdef int cmp(self, PartitionStack other, BinaryCode CG):
cdef int *self_wd_ents = self.wd_ents
cdef codeword *CG_words = CG.words
cdef int i, j, l, m, span = 1, ncols = self.ncols, nwords = self.nwords
for i from 0 < i < nwords:
for j from 0 <= j < ncols:
l = CG.is_one(self.wd_ents[i], self.col_ents[j])
m = CG.is_one(other.wd_ents[i], other.col_ents[j])
if l != m:
return l - m
return 0
def print_basis(self):
"""
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(4, 8)
sage: P._dangerous_dont_use_set_ents_lvls(range(8), range(7)+[-1], [4,7,12,11,1,9,3,0,2,5,6,8,10,13,14,15], [0]*16)
sage: P
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1,2,3,4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2,3,4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3,4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4},{5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4},{5},{6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4},{5},{6},{7})
sage: P._find_basis()
sage: P.print_basis()
basis_locations:
4
8
0
11
"""
cdef int i, j
if self.basis_locations:
print "basis_locations:"
j = 1
while (1 << j) < self.nwords:
j += 1
for i from 0 <= i < j:
print self.basis_locations[i]
def _find_basis(self):
"""
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: P = PartitionStack(4, 8)
sage: P._dangerous_dont_use_set_ents_lvls(range(8), range(7)+[-1], [4,7,12,11,1,9,3,0,2,5,6,8,10,13,14,15], [0]*16)
sage: P
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1,2,3,4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2,3,4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3,4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4,5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4},{5,6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4},{5},{6,7})
({4},{7},{12},{11},{1},{9},{3},{0},{2},{5},{6},{8},{10},{13},{14},{15}) ({0},{1},{2},{3},{4},{5},{6},{7})
sage: P._find_basis()
sage: P.print_basis()
basis_locations:
4
8
0
11
"""
cdef int i
cdef int *ham_wts = hamming_weights()
self.find_basis(ham_wts)
sage_free(ham_wts)
cdef int find_basis(self, int *ham_wts):
cdef int i = 0, j, k, nwords = self.nwords, weight, basis_elts = 0, nrows = self.nrows
cdef int *self_wd_ents = self.wd_ents
if self.basis_locations is NULL:
self.basis_locations = <int *> sage_malloc( 2 * nrows * sizeof(int) )
if self.basis_locations is NULL:
raise MemoryError("Memory.")
while i < nwords:
j = self_wd_ents[i]
weight = ham_wts[j & 65535] + ham_wts[(j>>16) & 65535]
if weight == 1:
basis_elts += 1
k = 0
while not (1<<k) & j:
k += 1
self.basis_locations[k] = i
if basis_elts == nrows: break
i += 1
for i from 0 <= i < nrows:
self.basis_locations[nrows + i] = self_wd_ents[1 << i]
def _get_permutation(self, other):
"""
EXAMPLE:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]])
sage: B = BinaryCode(M)
sage: P = PartitionStack(4, 8)
sage: P._refine(0, [[0,0],[1,0]], B)
181
sage: P._split_vertex(0, 1)
0
sage: P._refine(1, [[0,0]], B)
290
sage: P._split_vertex(1, 2)
1
sage: P._refine(2, [[0,1]], B)
463
sage: P._split_vertex(2, 3)
2
sage: P._refine(3, [[0,2]], B)
1500
sage: P._split_vertex(4, 4)
4
sage: P._refine(4, [[0,4]], B)
1224
sage: P._is_discrete(4)
1
sage: Q = PartitionStack(P)
sage: Q._clear(4)
sage: Q._split_vertex(5, 4)
4
sage: Q._refine(4, [[0,4]], B)
1224
sage: Q._is_discrete(4)
1
sage: P._get_permutation(Q)
([0, 1, 2, 3, 4, 5, 6, 7, 12, 13, 14, 15, 8, 9, 10, 11], [0, 1, 2, 3, 5, 4, 7, 6])
"""
cdef int i
cdef int *word_g = <int *> sage_malloc( self.nwords * sizeof(int) )
cdef int *col_g = <int *> sage_malloc( self.ncols * sizeof(int) )
if word_g is NULL or col_g is NULL:
if word_g is not NULL: sage_free(word_g)
if col_g is not NULL: sage_free(col_g)
raise MemoryError("Memory.")
self.get_permutation(other, word_g, col_g)
word_l = [word_g[i] for i from 0 <= i < self.nwords]
col_l = [col_g[i] for i from 0 <= i < self.ncols]
sage_free(word_g)
sage_free(col_g)
return word_l, col_l
cdef void get_permutation(self, PartitionStack other, int *word_gamma, int *col_gamma):
cdef int i
cdef int *self_wd_ents = self.wd_ents, *other_wd_ents = other.wd_ents
cdef int *self_col_ents = self.col_ents, *other_col_ents = other.col_ents
for i from 0 <= i < self.nwords:
word_gamma[other_wd_ents[i]] = self_wd_ents[i]
for i from 0 <= i < self.ncols:
col_gamma[other_col_ents[i]] = self_col_ents[i]
cdef class BinaryCodeClassifier:
def __cinit__(self):
self.radix = sizeof(codeword) << 3
self.ham_wts = hamming_weights()
self.L = 100
self.aut_gens_size = self.radix * 100
self.w_gamma_size = 1 << (self.radix/2)
self.alpha_size = self.w_gamma_size + self.radix
self.Phi_size = self.w_gamma_size/self.radix + 1
self.w_gamma = <int *> sage_malloc( self.w_gamma_size * sizeof(int) )
self.alpha = <int *> sage_malloc( self.alpha_size * sizeof(int) )
self.Phi = <unsigned int *> sage_malloc( self.Phi_size * (self.L+1) * sizeof(unsigned int) )
self.Omega = <unsigned int *> sage_malloc( self.Phi_size * self.L * sizeof(unsigned int) )
self.W = <unsigned int *> sage_malloc( self.Phi_size * self.radix * 2 * sizeof(unsigned int) )
self.base = <int *> sage_malloc( self.radix * sizeof(int) )
self.aut_gp_gens = <int *> sage_malloc( self.aut_gens_size * sizeof(int) )
self.c_gamma = <int *> sage_malloc( self.radix * sizeof(int) )
self.labeling = <int *> sage_malloc( self.radix * 3 * sizeof(int) )
self.Lambda1 = <int *> sage_malloc( self.radix * 2 * sizeof(int) )
self.Lambda2 = <int *> sage_malloc( self.radix * 2 * sizeof(int) )
self.Lambda3 = <int *> sage_malloc( self.radix * 2 * sizeof(int) )
self.v = <int *> sage_malloc( self.radix * 2 * sizeof(int) )
self.e = <int *> sage_malloc( self.radix * 2 * sizeof(int) )
if self.Phi is NULL or self.Omega is NULL or self.W is NULL or self.Lambda1 is NULL \
or self.Lambda2 is NULL or self.Lambda3 is NULL or self.w_gamma is NULL \
or self.c_gamma is NULL or self.alpha is NULL or self.v is NULL or self.e is NULL \
or self.aut_gp_gens is NULL or self.labeling is NULL or self.base is NULL:
if self.Phi is not NULL: sage_free(self.Phi)
if self.Omega is not NULL: sage_free(self.Omega)
if self.W is not NULL: sage_free(self.W)
if self.Lambda1 is not NULL: sage_free(self.Lambda1)
if self.Lambda2 is not NULL: sage_free(self.Lambda2)
if self.Lambda3 is not NULL: sage_free(self.Lambda3)
if self.w_gamma is not NULL: sage_free(self.w_gamma)
if self.c_gamma is not NULL: sage_free(self.c_gamma)
if self.alpha is not NULL: sage_free(self.alpha)
if self.v is not NULL: sage_free(self.v)
if self.e is not NULL: sage_free(self.e)
if self.aut_gp_gens is not NULL: sage_free(self.aut_gp_gens)
if self.labeling is not NULL: sage_free(self.labeling)
if self.base is not NULL: sage_free(self.base)
raise MemoryError("Memory.")
def __dealloc__(self):
sage_free(self.ham_wts)
sage_free(self.Phi)
sage_free(self.Omega)
sage_free(self.W)
sage_free(self.Lambda1)
sage_free(self.Lambda2)
sage_free(self.Lambda3)
sage_free(self.c_gamma)
sage_free(self.w_gamma)
sage_free(self.alpha)
sage_free(self.v)
sage_free(self.e)
sage_free(self.aut_gp_gens)
sage_free(self.labeling)
sage_free(self.base)
cdef void record_automorphism(self, int *gamma, int ncols):
cdef int i, j
if self.aut_gp_index + ncols > self.aut_gens_size:
self.aut_gens_size *= 2
self.aut_gp_gens = <int *> sage_realloc( self.aut_gp_gens, self.aut_gens_size * sizeof(int) )
if self.aut_gp_gens is NULL:
raise MemoryError("Memory.")
j = self.aut_gp_index
for i from 0 <= i < ncols:
self.aut_gp_gens[i+j] = gamma[i]
self.aut_gp_index += ncols
def _aut_gp_and_can_label(self, CC, verbosity=0):
"""
Compute the automorphism group and canonical label of the code CC.
INPUT:
CC - a BinaryCode object
verbosity - a nonnegative integer
OUTPUT:
a tuple, (gens, labeling, size, base)
gens -- a list of permutations (in list form) representing generators
of the permutation automorphism group of the code CC.
labeling -- a permutation representing the canonical labeling of the
code. mostly for internal use; entries describe the relabeling
on the columns.
size -- the order of the automorphism group.
base -- a set of cols whose action determines the action on all cols
EXAMPLES:
sage: import sage.coding.binary_code
sage: from sage.coding.binary_code import *
sage: BC = BinaryCodeClassifier()
sage: M = Matrix(GF(2),[
... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0],
... [0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0],
... [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1],
... [0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],
... [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]])
sage: B = BinaryCode(M)
sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
sage: S = SymmetricGroup(M.ncols())
sage: L = [S([x+1 for x in g]) for g in gens]
sage: PermutationGroup(L).order()
322560
sage: size
322560
sage: M = Matrix(GF(2),[
... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0],
... [0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0],
... [0,0,0,0,0,1,0,1,0,0,0,1,1,1,1,1,1],
... [0,0,0,1,1,0,0,0,0,1,1,0,1,1,0,1,1]])
sage: B = BinaryCode(M)
sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
sage: S = SymmetricGroup(M.ncols())
sage: L = [S([x+1 for x in g]) for g in gens]
sage: PermutationGroup(L).order()
2304
sage: size
2304
sage: M=Matrix(GF(2),[
... [1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0,0],
... [0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0],
... [0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0],
... [0,0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0],
... [0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0,0],
... [0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0],
... [0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0],
... [0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1]])
sage: B = BinaryCode(M)
sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
sage: S = SymmetricGroup(M.ncols())
sage: L = [S([x+1 for x in g]) for g in gens]
sage: PermutationGroup(L).order()
136
sage: size
136
sage: M=Matrix(GF(2),[
... [0,1,0,1,1,1,0,0,0,1,0,0,0,1,0,0,0,1,1,1,0,1],
... [1,0,1,1,1,0,0,0,1,0,0,0,1,0,0,0,1,1,1,0,1,0],
... [0,1,1,1,0,0,0,1,0,0,1,1,0,0,0,1,1,1,0,1,0,0],
... [1,1,1,0,0,0,1,0,0,1,0,0,0,0,1,1,1,0,1,0,0,1],
... [1,1,0,0,0,1,0,0,1,0,1,0,0,1,1,1,0,1,0,0,1,0],
... [1,0,0,0,1,0,0,1,0,1,1,0,1,1,1,0,1,0,0,1,0,0],
... [0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,0,1,0,0,0],
... [0,0,1,0,0,1,0,1,1,1,0,1,1,0,1,0,0,1,0,0,0,1],
... [0,1,0,0,1,0,1,1,1,0,0,1,0,1,0,0,1,0,0,0,1,1],
... [1,0,0,1,0,1,1,1,0,0,0,0,1,0,0,1,0,0,0,1,1,1],
... [0,0,1,0,1,1,1,0,0,0,1,1,0,0,1,0,0,0,1,1,1,0]])
sage: B = BinaryCode(M)
sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
sage: S = SymmetricGroup(M.ncols())
sage: L = [S([x+1 for x in g]) for g in gens]
sage: PermutationGroup(L).order()
887040
sage: size
887040
sage: B = BinaryCode(Matrix(GF(2),[[1,0,1],[0,1,1]]))
sage: BC._aut_gp_and_can_label(B)
([[0, 2, 1], [1, 0, 2]], [0, 1, 2], 6, [0, 1])
sage: B = BinaryCode(Matrix(GF(2),[[1,1,1,1]]))
sage: BC._aut_gp_and_can_label(B)
([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], [0, 1, 2, 3], 24, [0, 1, 2])
sage: B = BinaryCode(Matrix(GF(2),[[0,0,0,0,0,0,0,0,0,0,0,0,0,0,1]]))
sage: gens, labeling, size, base = BC._aut_gp_and_can_label(B)
sage: size
87178291200
sage: M = Matrix(GF(2),[
... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0],
... [0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0],
... [0,0,0,0,1,1,0,0,0,0,0,0,1,1,1,1,1,1],
... [0,0,1,1,0,0,0,0,0,0,1,1,1,1,0,0,1,1],
... [0,0,0,1,0,0,0,1,0,1,0,1,0,1,1,1,0,1],
... [0,1,0,0,0,1,0,0,0,1,1,1,0,1,0,1,1,0]])
sage: B = BinaryCode(M)
sage: BC._aut_gp_and_can_label(B)[2]
2160
sage: M = Matrix(GF(2),[
... [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
... [0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
... [0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0],
... [0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0],
... [0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0],
... [0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0],
... [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1],
... [1,0,1,0,1,0,1,0,1,1,0,0,0,0,0,0,1,1,0,0],
... [1,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,0,0],
... [1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,0,1,0]])
sage: B = BinaryCode(M)
sage: BC._aut_gp_and_can_label(B)[2]
294912
sage: M = Matrix(GF(2), [
... [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
... [0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
... [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0],
... [0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0],
... [0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1,0],
... [0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0],
... [0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1],
... [0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,1,0,0,1,1,1,0,1],
... [0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,1,1,1,1,0,0,0,1]])
sage: B = BinaryCode(M)
sage: BC = BinaryCodeClassifier()
sage: BC._aut_gp_and_can_label(B)[2]
442368
"""
cdef int i, j
cdef BinaryCode C = CC
self.aut_gp_and_can_label(C, verbosity)
i = 0
py_aut_gp_gens = []
while i < self.aut_gp_index:
gen = [self.aut_gp_gens[i+j] for j from 0 <= j < C.ncols]
py_aut_gp_gens.append(gen)
i += C.ncols
py_labeling = [self.labeling[i] for i from 0 <= i < C.ncols]
base = []
for i from 0 <= i < self.radix:
if self.base[i] == -1:
break
base.append(self.base[i])
aut_gp_size = self.aut_gp_size
return py_aut_gp_gens, py_labeling, aut_gp_size, base
cdef void aut_gp_and_can_label(self, BinaryCode C, int verbosity):
cdef int i, j, ii, jj, iii, jjj, iiii
cdef PartitionStack nu, zeta, rho
cdef int k = 0
cdef int k_rho
cdef int *v = self.v
cdef int h = -1
cdef int hb
cdef int hh = 1
cdef int ht
cdef int *alpha
cdef int tvc
cdef OrbitPartition Theta
cdef unsigned int *Phi
cdef unsigned int *Omega
cdef int l = -1
cdef unsigned int *W
cdef int *e
cdef int index = 0
cdef int *Lambda = self.Lambda1
cdef int *zf__Lambda_zeta = self.Lambda2
cdef int *zb__Lambda_rho = self.Lambda3
cdef int qzb
cdef int hzf__h_zeta
cdef int hzb__h_rho = -1
cdef int *word_gamma, *col_gamma = self.c_gamma
cdef int nwords = C.nwords, ncols = C.ncols, nrows = C.nrows
cdef int *ham_wts = self.ham_wts
cdef int state
self.aut_gp_index = 0
self.aut_gp_size = Integer(1)
if self.w_gamma_size < nwords:
while self.w_gamma_size < nwords:
self.w_gamma_size *= 2
self.alpha_size = self.w_gamma_size + self.radix
self.Phi_size = self.w_gamma_size/self.radix + 1
self.w_gamma = <int *> sage_realloc(self.w_gamma, self.w_gamma_size * sizeof(int) )
self.alpha = <int *> sage_realloc(self.alpha, self.alpha_size * sizeof(int) )
self.Phi = <unsigned int *> sage_realloc(self.Phi, self.Phi_size * self.L * sizeof(int) )
self.Omega = <unsigned int *> sage_realloc(self.Omega, self.Phi_size * self.L * sizeof(int) )
self.W = <unsigned int *> sage_realloc(self.W, self.Phi_size * self.radix * 2 * sizeof(int) )
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:
if self.w_gamma is not NULL: sage_free(self.w_gamma)
if self.alpha is not NULL: sage_free(self.alpha)
if self.Phi is not NULL: sage_free(self.Phi)
if self.Omega is not NULL: sage_free(self.Omega)
if self.W is not NULL: sage_free(self.W)
raise MemoryError("Memory.")
for i from 0 <= i < self.Phi_size * self.L:
self.Omega[i] = 0
word_gamma = self.w_gamma
alpha = self.alpha
Phi = self.Phi
Omega = self.Omega
W = self.W
e = self.e
nu = PartitionStack(nrows, ncols)
Theta = OrbitPartition(nrows, ncols)
if ncols == 0 or nrows == 0:
raise NotImplementedError("Must supply a nontrivial code.")
state = 1
while state != -1:
if False:
print '-----'
print "k:", k
print "h:", h
if False:
if k != -1:
if v[k]&nu.flag:
print "v[k]: word ", v[k]^nu.flag
else:
print "v[k]: col ", v[k]
if tvc&nu.flag:
print "tvc- wd", tvc^nu.flag
else:
print "tvc- col", tvc
if W[self.Phi_size * k]:
print "W[k]: cols", Integer(W[self.Phi_size * k]).binary()
else:
j = nwords/self.radix
if nwords%self.radix:
j += 1
L = ''
for i from 0 <= i < j:
if i == j - 1:
jj = nwords%self.radix
if jj == 0:
jj = self.radix
else:
jj = self.radix
for ii from 0 <= ii < jj:
if W[self.Phi_size * k + 1 + i] & (1 << ii):
L += '1'
else:
L += '0'
print "W[k]: words", L
if False:
print 'nu'
print nu
if tvc&nu.flag:
print 'tvc is word', tvc^nu.flag
else:
print 'tvc is col', tvc
if v[k]&nu.flag:
print 'v[k] is word', v[k]^nu.flag
else:
print 'v[k] is col', v[k]
if False:
if h != -1:
print 'zeta'
print zeta
print 'rho'
print rho
print "hzf:", hzf__h_zeta
print "aut_gp_index", self.aut_gp_index
print 'hh', hh
print 'ht', ht
print 'hzf__h_zeta', hzf__h_zeta
print 'qzb', qzb
if False:
print '-----'
print "state:", state
if state == 1:
alpha[0] = 0
alpha[1] = nu.flag
nu.refine(k, alpha, 2, C, ham_wts)
if nu.sat_225(k): hh = k
if nu.is_discrete(k): state = 18; continue
v[k] = nu.new_first_smallest_nontrivial(k, W, self.Phi_size * k)
Lambda[k] = 0
e[k] = 0
state = 2
elif state == 2:
k += 1
nu.clear(k)
alpha[0] = nu.split_vertex(v[k-1], k)
Lambda[k] = nu.refine(k, alpha, 1, C, ham_wts)
if h == -1: state = 5; continue
if hzf__h_zeta == k-1 and Lambda[k] == zf__Lambda_zeta[k]: hzf__h_zeta = k
if qzb == 0:
if zb__Lambda_rho[k] == -1 or Lambda[k] < zb__Lambda_rho[k]:
qzb = -1
elif Lambda[k] > zb__Lambda_rho[k]:
qzb = 1
else:
qzb = 0
if hzb__h_rho == k-1 and qzb == 0: hzb__h_rho = k
if qzb > 0: zb__Lambda_rho[k] = Lambda[k]
state = 3
elif state == 3:
if hzf__h_zeta <= k or qzb >= 0: state = 4
else: state = 6
elif state == 4:
if nu.is_discrete(k): state = 7; continue
v[k] = nu.new_first_smallest_nontrivial(k, W, self.Phi_size * k)
if not nu.sat_225(k): hh = k + 1
e[k] = 0
state = 2
elif state == 5:
zf__Lambda_zeta[k] = Lambda[k]
zb__Lambda_rho[k] = Lambda[k]
state = 4
elif state == 6:
j = k
if ht-1 > hzb__h_rho:
if ht-1 < hh-1:
k = ht-1
else:
k = hh-1
else:
if hzb__h_rho < hh-1:
k = hzb__h_rho
else:
k = hh-1
if k == -1: k = 0
if hb > k:
hb = k
if j == hh: state = 13; continue
if l < self.L-1: l += 1
nu.new_min_cell_reps(hh, Omega, self.Phi_size*l)
nu.fixed_vertices(hh, Phi, Omega, self.Phi_size*l)
state = 12
elif state == 7:
if h == -1: state = 18; continue
if k < hzf__h_zeta: state = 8; continue
nu.get_permutation(zeta, word_gamma, col_gamma)
if C.is_automorphism(col_gamma, word_gamma):
state = 10
else:
state = 8
elif state == 8:
if qzb < 0: state = 6; continue
if (qzb > 0) or (k < k_rho): state = 9; continue
j = nu.cmp(rho, C)
if j > 0: state = 9; continue
if j < 0: state = 6; continue
rho.get_permutation(nu, word_gamma, col_gamma)
state = 10
elif state == 9:
rho = PartitionStack(nu)
k_rho = k
qzb = 0
hb = k
hzb__h_rho = k
zb__Lambda_rho[k+1] = -1
state = 6
elif state == 10:
if l < self.L-1: l += 1
ii = self.Phi_size*l
jj = 1 + nwords/self.radix
for i from 0 <= i < jj:
Omega[ii+i] = ~0
Phi[ii+i] = 0
if nwords%self.radix:
jj += 1
i = 0
while i < ncols:
j = col_gamma[i]
while j != i:
Omega[ii] ^= (1<<j)
j = col_gamma[j]
i += 1
while i < ncols and not Omega[ii]&(1<<i):
i += 1
i = 0
jj = self.radix
while i < nwords:
j = word_gamma[i]
while j != i:
Omega[ii+1+j/jj] ^= (1<<(j%jj))
j = word_gamma[j]
i += 1
while i < nwords and not Omega[ii+1+i/jj]&(1<<(i%jj)):
i += 1
for i from 0 <= i < ncols:
if col_gamma[i] == i:
Phi[ii] ^= (1 << i)
for i from 0 <= i < nwords:
if word_gamma[i] == i:
Phi[ii+1+i/jj] ^= (1<<(i%jj))
j = Theta.merge_perm(col_gamma, word_gamma)
if not j: state = 11; continue
self.record_automorphism(col_gamma, ncols)
if tvc & nu.flag:
i = tvc^nu.flag
if Theta.wd_min_cell_rep[Theta.wd_find(i)] == i:
state = 11; continue
else:
if Theta.col_min_cell_rep[Theta.col_find(tvc)] == tvc:
state = 11; continue
k = h
state = 13
elif state == 11:
k = hb
state = 12
elif state == 12:
if e[k] == 1:
ii = self.Phi_size*l
jj = self.Phi_size*k
j = 1 + nwords/self.radix
if nwords%self.radix:
j += 1
W[jj] &= Omega[ii]
for i from 0 < i < j:
W[jj+i] &= Omega[ii+i]
state = 13
elif state == 13:
if k == -1: state = -1; continue
if k > h: state = 17; continue
if k == h: state = 14; continue
h = k
tvc = 0
jj = self.Phi_size*k
if W[jj]:
while not (1 << tvc) & W[jj]:
tvc += 1
else:
ii = 0
while not W[jj+1+ii]:
ii += 1
while not W[jj+1+ii] & (1 << tvc):
tvc += 1
tvc = (ii*self.radix + tvc) ^ nu.flag
state = 14
elif state == 14:
if v[k]&nu.flag == tvc&nu.flag:
if tvc&nu.flag:
if Theta.wd_find(v[k]^nu.flag) == Theta.wd_find(tvc^nu.flag):
index += 1
else:
if Theta.col_find(v[k]) == Theta.col_find(tvc):
index += 1
jj = self.Phi_size*k
if v[k]&nu.flag:
ii = self.radix
i = (v[k]^nu.flag) + 1
while i < nwords and not (1 << i%ii) & W[jj+1+i/ii]:
i += 1
if i < nwords:
v[k] = i^nu.flag
else:
state = 16; continue
if Theta.wd_min_cell_rep[Theta.wd_find(i)] == i:
state = 15
else:
state = 14
else:
i = v[k] + 1
while i < ncols and not (1 << i) & W[jj]:
i += 1
if i < ncols:
v[k] = i
else:
state = 16; continue
if Theta.col_min_cell_rep[Theta.col_find(v[k])] == v[k]:
state = 15
else:
state = 14
elif state == 15:
if k + 1 < hh:
hh = k + 1
if k < hzf__h_zeta:
hzf__h_zeta = k
if hzb__h_rho >= k:
hzb__h_rho = k
qzb = 0
state = 2
elif state == 16:
jj = self.Phi_size*k
if W[jj]:
i = W[jj]
j = ham_wts[i & 65535] + ham_wts[(i >> 16) & 65535]
else:
i = 0; j = 0
ii = self.radix
while i*ii < nwords:
iii = W[jj+1+i]
j += ham_wts[iii & 65535] + ham_wts[(iii >> 16) & 65535]
i += 1
if j == index and ht == k + 1: ht = k
self.aut_gp_size *= index
index = 0
k -= 1
if hb > k:
hb = k
state = 13
elif state == 17:
jjj = self.Phi_size*k
if e[k] == 0:
jj = self.Phi_size*self.L
iii = nwords/self.radix
if nwords%self.radix:
iii += 1
for ii from 0 <= ii < iii:
Phi[jj+ii] = 0
for ii from 0 <= ii < k:
if v[ii]&nu.flag:
i = v[ii]^nu.flag
Phi[jj+1+i/self.radix] ^= (1 << i%self.radix)
else:
Phi[jj] ^= (1 << v[ii])
for i from 0 <= i <= l:
ii = self.Phi_size*i
iiii = 1
for j from 0 <= j < iii:
if Phi[ii + j] & Phi[jj + j] != Phi[jj + j]:
iiii = 0
break
if iiii:
for j from 0 <= j < iii:
W[jjj + j] &= Omega[ii + j]
e[k] = 1
if nu.flag&v[k]:
i = (v[k]^nu.flag)
while i < nwords:
i += 1
if (1 << i%self.radix) & W[jjj+1+i/self.radix]: break
if i < nwords:
v[k] = i^nu.flag
state = 15; continue
else:
i = v[k]
while i < ncols:
i += 1
if (1 << i) & W[jjj]: break
if i < ncols:
v[k] = i
state = 15; continue
k -= 1
state = 13
elif state == 18:
h = k
ht = k
hzf__h_zeta = k
zeta = PartitionStack(nu)
for i from 0 <= i < k:
self.base[i] = v[i]
self.base_size = k
if k != self.radix:
self.base[k] = -1
k -= 1
rho = PartitionStack(nu)
k_rho = k+1
hzb__h_rho = k
hb = k
qzb = 0
state = 13
rho.find_basis(ham_wts)
for i from 0 <= i < ncols:
self.labeling[rho.col_ents[i]] = i
for i from 0 <= i < 2*nrows:
self.labeling[i+ncols] = rho.basis_locations[i]
def put_in_canonical_form(self, BinaryCode B):
"""
Puts the code into canonical form.
Canonical form is obtained by performing row reduction, permuting the
pivots to the front so that the generator matrix is of the form: the
identity matrix augmented to the right by arbitrary data.
EXAMPLE:
sage: from sage.coding.binary_code import *
sage: BC = BinaryCodeClassifier()
sage: B = BinaryCode(codes.ExtendedBinaryGolayCode().gen_mat())
sage: B.apply_permutation(range(24,-1,-1))
sage: B
Binary [24,12] linear code, generator matrix
[011000111010100000000000]
[001001001111100000000001]
[011010100101100000000010]
[001101110001100000000100]
[010011011001100000001000]
[010110110011000000010000]
[011101100110000000100000]
[000011110110100001000000]
[000111101101000010000000]
[001111011010000100000000]
[010110001110101000000000]
[011100011101010000000000]
sage: BC.put_in_canonical_form(B)
sage: B
Binary [24,12] linear code, generator matrix
[100000000000001100111001]
[010000000000001010001111]
[001000000000001111010010]
[000100000000010110101010]
[000010000000010110010101]
[000001000000010001101101]
[000000100000011000110110]
[000000010000011111001001]
[000000001000010101110011]
[000000000100010011011110]
[000000000010001011110101]
[000000000001001101101110]
"""
aut_gp_gens, labeling, size, base = self._aut_gp_and_can_label(B)
B._apply_permutation_to_basis(labeling)
B.put_in_std_form()
def generate_children(self, BinaryCode B, int n, int d=2):
"""
Use canonical augmentation to generate children of the code B.
INPUT:
B -- a BinaryCode
n -- limit on the degree of the code
d -- test whether new vector has weight divisible by d. If d==4, this
ensures that all doubly-even canonically augmented children are
generated.
EXAMPLE:
sage: from sage.coding.binary_code import *
sage: BC = BinaryCodeClassifier()
sage: B = BinaryCode(Matrix(GF(2), [[1,1,1,1]]))
sage: BC.generate_children(B, 6, 4)
[
[1 1 1 1 0 0]
[0 1 0 1 1 1]
]
NOTE:
The function self_orthogonal_binary_codes makes heavy use of this
function.
MORE EXAMPLES:
sage: soc_iter = self_orthogonal_binary_codes(12, 6, 4)
sage: L = list(soc_iter)
sage: for n in range(0, 13):
... s = 'n=%2d : '%n
... for k in range(1,7):
... s += '%3d '%len([C for C in L if C.length() == n and C.dimension() == k])
... print s
n= 0 : 0 0 0 0 0 0
n= 1 : 0 0 0 0 0 0
n= 2 : 0 0 0 0 0 0
n= 3 : 0 0 0 0 0 0
n= 4 : 1 0 0 0 0 0
n= 5 : 0 0 0 0 0 0
n= 6 : 0 1 0 0 0 0
n= 7 : 0 0 1 0 0 0
n= 8 : 1 1 1 1 0 0
n= 9 : 0 0 0 0 0 0
n=10 : 0 1 1 1 0 0
n=11 : 0 0 1 1 0 0
n=12 : 1 2 3 4 2 0
"""
cdef BinaryCode m
cdef codeword *ortho_basis, *B_can_lab, current, swap
cdef codeword word, temp, gate, nonzero_gate, orbit, bwd, k_gate
cdef codeword *temp_basis, *orbit_checks, orb_chx_size, orb_chx_shift, radix_gate
cdef WordPermutation *gwp, *hwp, *can_lab, *can_lab_inv
cdef WordPermutation **parent_generators
cdef BinaryCode B_aug
cdef int i, ii, j, jj, ij, k = 0, parity, combo, num_gens
cdef int base_size, *multimod2_index, row
cdef int *ham_wts = self.ham_wts
cdef int *num_inner_gens, *num_outer_gens, *v, log_2_radix
cdef bint bingo, bingo2, bingo3
B.put_in_std_form()
ortho_basis = expand_to_ortho_basis(B, n)
aut_gp_gens, labeling, size, base = self._aut_gp_and_can_label(B)
B_can_lab = <codeword *> sage_malloc(B.nrows * sizeof(codeword))
can_lab = create_word_perm(labeling[:B.ncols])
if B_can_lab is NULL or can_lab is NULL:
sage_free(ortho_basis)
if B_can_lab is not NULL:
sage_free(B_can_lab)
if can_lab is not NULL:
sage_free(can_lab)
raise MemoryError()
for i from 0 <= i < B.nrows:
B_can_lab[i] = permute_word_by_wp(can_lab, B.basis[i])
dealloc_word_perm(can_lab)
row = 0
current = 1
while row < B.nrows:
i = row
while i < B.nrows and not B_can_lab[i] & current:
i += 1
if i < B.nrows:
if i != row:
swap = B_can_lab[row]
B_can_lab[row] = B_can_lab[i]
B_can_lab[i] = swap
for j from 0 <= j < row:
if B_can_lab[j] & current:
B_can_lab[j] ^= B_can_lab[row]
for j from row < j < B.nrows:
if B_can_lab[j] & current:
B_can_lab[j] ^= B_can_lab[row]
row += 1
current = current << 1
num_gens = len(aut_gp_gens)
base_size = len(base)
parent_generators = <WordPermutation **> sage_malloc( len(aut_gp_gens) * sizeof(WordPermutation*) )
temp_basis = <codeword *> sage_malloc( self.radix * sizeof(codeword) )
output = []
for i from 0 <= i < len(aut_gp_gens):
parent_generators[i] = create_word_perm(aut_gp_gens[i] + range(B.ncols, n))
word = 0
while ortho_basis[k] & (((<codeword>1) << B.ncols) - 1):
k += 1
j = k
while ortho_basis[j]:
word ^= ortho_basis[j]
j += 1
log_2_radix = 0
while ((<codeword>1) << log_2_radix) < self.radix:
log_2_radix += 1
if k < log_2_radix:
orb_chx_size = 0
else:
orb_chx_size = k - log_2_radix
orbit_checks = <codeword *> sage_malloc( ((<codeword>1) << orb_chx_size) * sizeof(codeword) )
if orbit_checks is NULL:
raise MemoryError()
for temp from 0 <= temp < ((<codeword>1) << orb_chx_size):
orbit_checks[temp] = 0
combo = 0
parity = 0
gate = (<codeword>1 << B.nrows) - 1
k_gate = (<codeword>1 << k) - 1
nonzero_gate = ( (<codeword>1 << (n-B.ncols)) - 1 ) << B.ncols
radix_gate = (((<codeword>1) << log_2_radix) - 1)
while 1:
if nonzero_gate & word == nonzero_gate and \
(ham_wts[word & 65535] + ham_wts[(word >> 16) & 65535])%d == 0:
temp = (word >> B.nrows) & ((<codeword>1 << k) - 1)
if not orbit_checks[temp >> log_2_radix] & ((<codeword>1) << (temp & radix_gate)):
B_aug = BinaryCode(B, word)
aug_aut_gp_gens, aug_labeling, aug_size, aug_base = self._aut_gp_and_can_label(B_aug)
can_lab = create_word_perm(aug_labeling[:n])
can_lab_inv = create_inv_word_perm(can_lab)
for j from 0 <= j < B_aug.nrows:
temp_basis[j] = permute_word_by_wp(can_lab, B_aug.basis[j])
i = 0
j = 0
while j < B_aug.nrows:
ii = j
while ii < B_aug.nrows and not temp_basis[ii] & (<codeword>1 << i):
ii += 1
if ii != B_aug.nrows:
if ii != j:
swap = temp_basis[ii]
temp_basis[ii] = temp_basis[j]
temp_basis[j] = swap
for jj from 0 <= jj < j:
if temp_basis[jj] & (<codeword>1 << i):
temp_basis[jj] ^= temp_basis[j]
for jj from j < jj < B_aug.nrows:
if temp_basis[jj] & (<codeword>1 << i):
temp_basis[jj] ^= temp_basis[j]
j += 1
i += 1
for j from 0 <= j < B.nrows:
temp_basis[j] = permute_word_by_wp(can_lab_inv, temp_basis[j])
from sage.matrix.constructor import matrix
from sage.rings.all import ZZ
from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroupElement
from sage.interfaces.gap import gap
rs = []
for i from 0 <= i < B.nrows:
r = []
for j from 0 <= j < n:
r.append((((<codeword>1)<<j)&temp_basis[i])>>j)
rs.append(r)
m = BinaryCode(matrix(ZZ, rs))
m_aut_gp_gens, m_labeling, m_size, m_base = self._aut_gp_and_can_label(m)
from sage.rings.arith import factorial
if True:
if len(m_aut_gp_gens) == 0:
aut_m = PermutationGroup([()])
else:
aut_m = PermutationGroup([PermutationGroupElement([a+1 for a in g]) for g in m_aut_gp_gens])
if len(aug_aut_gp_gens) == 0:
aut_B_aug = PermutationGroup([()])
else:
aut_B_aug = PermutationGroup([PermutationGroupElement([a+1 for a in g]) for g in aug_aut_gp_gens])
H = aut_m._gap_(gap).Intersection2(aut_B_aug._gap_(gap))
rt_transversal = list(gap('List(RightTransversal( %s,%s ));'\
%(str(aut_B_aug.__interface[gap]),str(H))))
rt_transversal = [PermutationGroupElement(g) for g in rt_transversal if str(g) != '()']
rt_transversal = [[a-1 for a in g.domain()] for g in rt_transversal]
rt_transversal = [g + range(len(g), n) for g in rt_transversal]
rt_transversal.append(range(n))
bingo2 = 0
for coset_rep in rt_transversal:
hwp = create_word_perm(coset_rep)
bingo2 = 1
for j from 0 <= j < B.nrows:
temp = permute_word_by_wp(hwp, temp_basis[j])
if temp != B.words[temp & gate]:
bingo2 = 0
dealloc_word_perm(hwp)
break
if bingo2:
dealloc_word_perm(hwp)
break
if bingo2:
from sage.matrix.constructor import Matrix
from sage.rings.all import GF
M = Matrix(GF(2), B_aug.nrows, B_aug.ncols)
for i from 0 <= i < B_aug.ncols:
for j from 0 <= j < B_aug.nrows:
M[j,i] = B_aug.is_one(1 << j, i)
output.append(M)
dealloc_word_perm(can_lab)
dealloc_word_perm(can_lab_inv)
orbits = [word]
j = 0
while j < len(orbits):
for i from 0 <= i < len(aut_gp_gens):
temp = <codeword> orbits[j]
temp = permute_word_by_wp(parent_generators[i], temp)
temp ^= B.words[temp & gate]
if temp not in orbits:
orbits.append(temp)
j += 1
for temp in orbits:
temp = (temp >> B.nrows) & k_gate
orbit_checks[temp >> log_2_radix] |= ((<codeword>1) << (temp & radix_gate))
parity ^= 1
i = 0
if not parity:
while not combo & (1 << i): i += 1
i += 1
if i == k: break
else:
combo ^= (1 << i)
word ^= ortho_basis[i]
for i from 0 <= i < len(aut_gp_gens):
dealloc_word_perm(parent_generators[i])
sage_free(B_can_lab)
sage_free(parent_generators)
sage_free(orbit_checks)
sage_free(ortho_basis)
sage_free(temp_basis)
return output
