Ticket #4115: trac_4115-double-cosets.patch

sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd

@@ -8 +8 @@
 
 include '../../../ext/cdefs.pxi'
 include '../../../ext/stdsage.pxi'
-include '../../../misc/bitset_pxd.pxi'
+include 'data_structures_pxd.pxi' # includes bitsets
 
 from sage.rings.integer cimport Integer
-
-cdef struct OrbitPartition:
-    # Disjoint set data structure for representing the orbits of the generators
-    # so far found. Also keeps track of the minimum elements of the cells and
-    # their sizes.
-    int degree
-    int *parent
-    int *rank
-    int *mcr # minimum cell representatives - only valid at the root of a cell
-    int *size # also only valid at the root of a cell
-
-cdef inline OrbitPartition *OP_new(int)
-cdef inline int OP_dealloc(OrbitPartition *)
-cdef inline int OP_find(OrbitPartition *, int)
-cdef inline int OP_join(OrbitPartition *, int, int)
-cdef inline int OP_merge_list_perm(OrbitPartition *, int *)
-
-cdef struct PartitionStack:
-    # Representation of a node of the search tree. A sequence of partitions of
-    # length depth + 1, each of which is finer than the last. Partition k is
-    # represented as PS.entries in order, broken up immediately after each
-    # entry of levels which is at most k.
-    int *entries
-    int *levels
-    int depth
-    int degree
-
-cdef inline PartitionStack *PS_new(int, bint)
-cdef inline PartitionStack *PS_copy(PartitionStack *)
-cdef inline int PS_copy_from_to(PartitionStack *, PartitionStack *)
-cdef inline PartitionStack *PS_from_list(object)
-cdef inline int PS_dealloc(PartitionStack *)
-cdef inline int PS_print(PartitionStack *)
-cdef inline int PS_print_partition(PartitionStack *, int)
-cdef inline bint PS_is_discrete(PartitionStack *)
-cdef inline int PS_num_cells(PartitionStack *)
-cdef inline int PS_move_min_to_front(PartitionStack *, int, int)
-cdef inline bint PS_is_mcr(PartitionStack *, int)
-cdef inline bint PS_is_fixed(PartitionStack *, int)
-cdef inline int PS_clear(PartitionStack *)
-cdef inline int PS_split_point(PartitionStack *, int)
-cdef inline int PS_first_smallest(PartitionStack *, bitset_t)
-cdef inline int PS_get_perm_from(PartitionStack *, PartitionStack *, int *)
-
-cdef inline int split_point_and_refine(PartitionStack *, int, object,
-    int (*)(PartitionStack *, object, int *, int), int *)
 
 cdef struct aut_gp_and_can_lab_return:
     int *generators
@@ -69 +23 @@
 cdef aut_gp_and_can_lab_return *get_aut_gp_and_can_lab( object, int **, int,
     bint (*)(PartitionStack *, object),
     int (*)(PartitionStack *, object, int *, int),
-    int (*)(int *, int *, object), bint, bint, bint)
+    int (*)(int *, int *, object, object), bint, bint, bint)
 
 
 

sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx

@@ -58 +58 @@
     
     Signature:
     
-    \code{int compare_structures(int *gamma_1, int *gamma_2, object S)}
+    \code{int compare_structures(int *gamma_1, int *gamma_2, object S1, object S2)}
     
     This function must implement a total ordering on the set of objects of fixed
-    order. Return -1 if \code{gamma_1(S) < gamma_2(S)}, 0 if
-    \code{gamma_1(S) == gamma_2(S)}, 1 if \code{gamma_1(S) > gamma_2(S)}.
+    order. Return -1 if \code{gamma_1(S1) < gamma_2(S2)}, 0 if
+    \code{gamma_1(S1) == gamma_2(S2)}, 1 if \code{gamma_1(S1) > gamma_2(S2)}.
 
 C. \code{all_children_are_equivalent}:
     
@@ -93 +93 @@
 #                         http://www.gnu.org/licenses/
 #*****************************************************************************
 
-include '../../../misc/bitset.pxi'
-
-cdef inline int min(int a, int b):
-    if a < b: return a
-    else: return b
-
-cdef inline int max(int a, int b):
-    if a > b: return a
-    else: return b
+include 'data_structures_pyx.pxi' # includes bitsets
 
 cdef inline int cmp(int a, int b):
     if a < b: return -1
     elif a == b: return 0
     else: return 1
 
-# OrbitPartitions
-
-cdef inline OrbitPartition *OP_new(int n):
-    """
-    Allocate and return a pointer to a new OrbitPartition of degree n. Returns a
-    null pointer in the case of an allocation failure.
-    """
-    cdef int i
-    cdef OrbitPartition *OP = <OrbitPartition *> \
-                                sage_malloc(sizeof(OrbitPartition))
-    if OP is NULL:
-        return OP
-    OP.degree = n
-    OP.parent = <int *> sage_malloc( n * sizeof(int) )
-    OP.rank = <int *> sage_malloc( n * sizeof(int) )
-    OP.mcr = <int *> sage_malloc( n * sizeof(int) )
-    OP.size = <int *> sage_malloc( n * sizeof(int) )
-    if OP.parent is NULL or OP.rank is NULL or \
-       OP.mcr is NULL or OP.size is NULL:
-        if OP.parent is not NULL: sage_free(OP.parent)
-        if OP.rank is not NULL: sage_free(OP.rank)
-        if OP.mcr is not NULL: sage_free(OP.mcr)
-        if OP.size is not NULL: sage_free(OP.size)
-        return NULL
-    for i from 0 <= i < n:
-        OP.parent[i] = i
-        OP.rank[i] = 0
-        OP.mcr[i] = i
-        OP.size[i] = 1
-    return OP
-
-cdef inline int OP_dealloc(OrbitPartition *OP):
-    sage_free(OP.parent)
-    sage_free(OP.rank)
-    sage_free(OP.mcr)
-    sage_free(OP.size)
-    sage_free(OP)
-    return 0
-
-cdef inline int OP_find(OrbitPartition *OP, int n):
-    """
-    Report the representative ("root") of the cell which contains n.
-    """
-    if OP.parent[n] == n:
-        return n
-    else:
-        OP.parent[n] = OP_find(OP, OP.parent[n])
-        return OP.parent[n]
-
-cdef inline int OP_join(OrbitPartition *OP, int m, int n):
-    """
-    Join the cells containing m and n, if they are different.
-    """
-    cdef int m_root = OP_find(OP, m)
-    cdef int n_root = OP_find(OP, n)
-    if OP.rank[m_root] > OP.rank[n_root]:
-        OP.parent[n_root] = m_root
-        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
-        OP.size[m_root] += OP.size[n_root]
-    elif OP.rank[m_root] < OP.rank[n_root]:
-        OP.parent[m_root] = n_root
-        OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root])
-        OP.size[n_root] += OP.size[m_root]
-    elif m_root != n_root:
-        OP.parent[n_root] = m_root
-        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
-        OP.size[m_root] += OP.size[n_root]
-        OP.rank[m_root] += 1
-
-cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma):
-    """
-    Joins the cells of OP which intersect the same orbit of gamma.
-    
-    INPUT:
-        gamma - an integer array representing i -> gamma[i].
-    
-    OUTPUT:
-        1 - something changed
-        0 - orbits of gamma all contained in cells of OP
-    """
-    cdef int i, i_root, gamma_i_root, changed = 0
-    for i from 0 <= i < OP.degree:
-        if gamma[i] == i: continue
-        i_root = OP_find(OP, i)
-        gamma_i_root = OP_find(OP, gamma[i])
-        if i_root != gamma_i_root:
-            changed = 1
-            OP_join(OP, i_root, gamma_i_root)
-    return changed
-
-def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]):
-    """
-    Demonstration and testing.
-    
-    DOCTEST:
-        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
-        ... import OP_represent
-        sage: OP_represent()
-        Allocating OrbitPartition...
-        Allocation passed.
-        Checking that each element reports itself as its root.
-        Each element reports itself as its root.
-        Merging:
-        Merged 0 and 1.
-        Merged 2 and 3.
-        Merged 3 and 4.
-        Done merging.
-        Finding:
-        0 -> 0, root: size=2, mcr=0, rank=1
-        1 -> 0
-        2 -> 2, root: size=3, mcr=2, rank=1
-        3 -> 2
-        4 -> 2
-        5 -> 5, root: size=1, mcr=5, rank=0
-        6 -> 6, root: size=1, mcr=6, rank=0
-        7 -> 7, root: size=1, mcr=7, rank=0
-        8 -> 8, root: size=1, mcr=8, rank=0
-        Allocating array to test merge_perm.
-        Allocation passed.
-        Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8]
-        Done merging.
-        Finding:
-        0 -> 0, root: size=5, mcr=0, rank=2
-        1 -> 0
-        2 -> 0
-        3 -> 0
-        4 -> 0
-        5 -> 5, root: size=3, mcr=5, rank=1
-        6 -> 5
-        7 -> 5
-        8 -> 8, root: size=1, mcr=8, rank=0
-        Deallocating OrbitPartition.
-        Done.
-
-    """
-    cdef int i
-    print "Allocating OrbitPartition..."
-    cdef OrbitPartition *OP = OP_new(n)
-    if OP is NULL:
-        print "Allocation failed!"
-        return
-    print "Allocation passed."
-    print "Checking that each element reports itself as its root."
-    good = True
-    for i from 0 <= i < n:
-        if not OP_find(OP, i) == i:
-            print "Failed at i = %d!"%i
-            good = False
-    if good: print "Each element reports itself as its root."
-    print "Merging:"
-    for i,j in merges:
-        OP_join(OP, i, j)
-        print "Merged %d and %d."%(i,j)
-    print "Done merging."
-    print "Finding:"
-    for i from 0 <= i < n:
-        j = OP_find(OP, i)
-        s = "%d -> %d"%(i, j)
-        if i == j:
-            s += ", root: size=%d, mcr=%d, rank=%d"%\
-                   (OP.size[i], OP.mcr[i], OP.rank[i])
-        print s
-    print "Allocating array to test merge_perm."
-    cdef int *gamma = <int *> sage_malloc( n * sizeof(int) )
-    if gamma is NULL:
-        print "Allocation failed!"
-        OP_dealloc(OP)
-        return
-    print "Allocation passed."
-    for i from 0 <= i < n:
-        gamma[i] = perm[i]
-    print "Merging permutation: %s"%perm
-    OP_merge_list_perm(OP, gamma)
-    print "Done merging."
-    print "Finding:"
-    for i from 0 <= i < n:
-        j = OP_find(OP, i)
-        s = "%d -> %d"%(i, j)
-        if i == j:
-            s += ", root: size=%d, mcr=%d, rank=%d"%\
-                   (OP.size[i], OP.mcr[i], OP.rank[i])
-        print s    
-    print "Deallocating OrbitPartition."
-    sage_free(gamma)
-    OP_dealloc(OP)
-    print "Done."
-
-# PartitionStacks
-
-cdef inline PartitionStack *PS_new(int n, bint unit_partition):
-    """
-    Allocate and return a pointer to a new PartitionStack of degree n. Returns a
-    null pointer in the case of an allocation failure.
-    """
-    cdef int i
-    cdef PartitionStack *PS = <PartitionStack *> \
-                                sage_malloc(sizeof(PartitionStack))
-    if PS is NULL:
-        return PS
-    PS.entries = <int *> sage_malloc( n * sizeof(int) )
-    PS.levels = <int *> sage_malloc( n * sizeof(int) )
-    if PS.entries is NULL or PS.levels is NULL:
-        if PS.entries is not NULL: sage_free(PS.entries)
-        if PS.levels is not NULL: sage_free(PS.levels)
-        return NULL
-    PS.depth = 0
-    PS.degree = n
-    if unit_partition:
-        for i from 0 <= i < n-1:
-            PS.entries[i] = i
-            PS.levels[i] = n
-        PS.entries[n-1] = n-1
-        PS.levels[n-1] = -1
-    return PS
-
-cdef inline PartitionStack *PS_copy(PartitionStack *PS):
-    """
-    Allocate and return a pointer to a copy of PartitionStack PS. Returns a null
-    pointer in the case of an allocation failure.
-    """
-    cdef int i
-    cdef PartitionStack *PS2 = <PartitionStack *> \
-                                sage_malloc(sizeof(PartitionStack))
-    if PS2 is NULL:
-        return PS2
-    PS2.entries = <int *> sage_malloc( PS.degree * sizeof(int) )
-    PS2.levels = <int *> sage_malloc( PS.degree * sizeof(int) )
-    if PS2.entries is NULL or PS2.levels is NULL:
-        if PS2.entries is not NULL: sage_free(PS2.entries)
-        if PS2.levels is not NULL: sage_free(PS2.levels)
-        return NULL
-    PS_copy_from_to(PS, PS2)
-    return PS2
-
-cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2):
-    """
-    Copy all data from PS to PS2.
-    """
-    PS2.depth = PS.depth
-    PS2.degree = PS.degree
-    for i from 0 <= i < PS.degree:
-        PS2.entries[i] = PS.entries[i]
-        PS2.levels[i] = PS.levels[i]    
-
-cdef inline PartitionStack *PS_from_list(object L):
-    """
-    Allocate and return a pointer to a PartitionStack representing L. Returns a
-    null pointer in the case of an allocation failure.
-    """
-    cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0
-    for cell from 0 <= cell < num_cells:
-        n += len(L[cell])
-    cdef PartitionStack *PS = PS_new(n, 0)
-    cell = 0
-    if PS is NULL:
-        return PS
-    while 1:
-        cur_len = len(L[cell])
-        for i from 0 <= i < cur_len:
-            PS.entries[cur_start + i] = L[cell][i]
-            PS.levels[cur_start + i] = n
-        PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1)
-        cur_start += cur_len
-        cell += 1
-        if cell == num_cells:
-            PS.levels[cur_start-1] = -1
-            break
-        PS.levels[cur_start-1] = 0
-    return PS
-
-cdef inline int PS_dealloc(PartitionStack *PS):
-    sage_free(PS.entries)
-    sage_free(PS.levels)
-    sage_free(PS)
-    return 0
-
-cdef inline int PS_print(PartitionStack *PS):
-    """
-    Print a visual representation of PS.
-    """
-    cdef int i
-    for i from 0 <= i < PS.depth + 1:
-        PS_print_partition(PS, i)
-    print
-
-cdef inline int PS_print_partition(PartitionStack *PS, int k):
-    """
-    Print the partition at depth k.
-    """
-    s = '('
-    for i from 0 <= i < PS.degree:
-        s += str(PS.entries[i])
-        if PS.levels[i] <= k:
-            s += '|'
-        else:
-            s += ','
-    s = s[:-1] + ')'
-    print s
-
-cdef inline bint PS_is_discrete(PartitionStack *PS):
-    """
-    Returns whether the deepest partition consists only of singleton cells.
-    """
-    cdef int i
-    for i from 0 <= i < PS.degree:
-        if PS.levels[i] > PS.depth:
-            return 0
-    return 1
-
-cdef inline int PS_num_cells(PartitionStack *PS):
-    """
-    Returns the number of cells.
-    """
-    cdef int i, ncells = 0
-    for i from 0 <= i < PS.degree:
-        if PS.levels[i] <= PS.depth:
-            ncells += 1
-    return ncells
-
-cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end):
-    """
-    Makes sure that the first element of the segment of entries i with
-    start <= i <= end is minimal.
-    """
-    cdef int i, min_loc = start, minimum = PS.entries[start]
-    for i from start < i <= end:
-        if PS.entries[i] < minimum:
-            min_loc = i
-            minimum = PS.entries[i]
-    if min_loc != start:
-        PS.entries[min_loc] = PS.entries[start]
-        PS.entries[start] = minimum
-
-cdef inline bint PS_is_mcr(PartitionStack *PS, int m):
-    """
-    Returns whether PS.elements[m] (not m!) is the smallest element of its cell.
-    """
-    return m == 0 or PS.levels[m-1] <= PS.depth
-
-cdef inline bint PS_is_fixed(PartitionStack *PS, int m):
-    """
-    Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming
-    PS_is_mcr(PS, m) is already true.
-    """
-    return PS.levels[m] <= PS.depth
-
-cdef inline int PS_clear(PartitionStack *PS):
-    """
-    Sets the current partition to the first shallower one, i.e. forgets about
-    boundaries between cells that are new to the current level.
-    """
-    cdef int i, cur_start = 0
-    for i from 0 <= i < PS.degree:
-        if PS.levels[i] >= PS.depth:
-            PS.levels[i] += 1
-        if PS.levels[i] < PS.depth:
-            PS_move_min_to_front(PS, cur_start, i)
-            cur_start = i+1
-
-cdef inline int PS_split_point(PartitionStack *PS, int v):
-    """
-    Detaches the point v from the cell it is in, putting the singleton cell
-    of just v in front. Returns the position where v is now located.
-    """
-    cdef int i = 0, index_of_v
-    while PS.entries[i] != v:
-        i += 1
-    index_of_v = i
-    while PS.levels[i] > PS.depth:
-        i += 1
-    if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \
-       and PS.levels[index_of_v] > PS.depth:
-        # if v is first (make sure v is not already alone)
-        PS_move_min_to_front(PS, index_of_v+1, i)
-        PS.levels[index_of_v] = PS.depth
-        return index_of_v
-    else:
-        # If v is not at front, i.e. v is not minimal in its cell,
-        # then move_min_to_front is not necessary since v will swap
-        # with the first before making its own cell, leaving it at
-        # the front of the other.
-        i = index_of_v
-        while i != 0 and PS.levels[i-1] > PS.depth:
-            i -= 1
-        PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v
-        PS.entries[i+1] = PS.entries[i] # move the first (min) to second
-        PS.entries[i] = v # place v first
-        PS.levels[i] = PS.depth
-        return i
-
-cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b):
-    """
-    Find the first occurrence of the smallest cell of size greater than one,
-    store its entries to b, and return its minimum element.
-    """
-    cdef int i = 0, j = 0, location = 0, n = PS.degree
-    bitset_zero(b)
-    while 1:
-        if PS.levels[i] <= PS.depth:
-            if i != j and n > i - j + 1:
-                n = i - j + 1
-                location = j
-            j = i + 1
-        if PS.levels[i] == -1: break
-        i += 1
-    # location now points to the beginning of the first, smallest,
-    # nontrivial cell
-    i = location
-    while 1:
-        bitset_flip(b, PS.entries[i])
-        if PS.levels[i] <= PS.depth: break
-        i += 1
-    return PS.entries[location]
-
-cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma):
-    """
-    Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i]
-    denotes the entry of the ith cell of the discrete partition PS.
-    """
-    cdef int i = 0
-    for i from 0 <= i < PS1.degree:
-        gamma[PS2.entries[i]] = PS1.entries[i]
-
-def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]):
-    """
-    Demonstration and testing.
-    
-    DOCTEST:
-        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
-        ... import PS_represent
-        sage: PS_represent()
-        Allocating PartitionStack...
-        Allocation passed:
-        (0,1,2,3,4,5,6,7,8,9)
-        Checking that entries are in order and correct level.
-        Everything seems in order, deallocating.
-        Deallocated.
-        Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]].
-        PartitionStack's data:
-        entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2]
-        levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1]
-        depth = 0, degree = 10
-        (6|3,4,8,7|1,9,5|0,2)
-        Checking PS_is_discrete:
-        False
-        Checking PS_num_cells:
-        4
-        Checking PS_is_mcr, min cell reps are:
-        [6, 3, 1, 0]
-        Checking PS_is_fixed, fixed elements are:
-        [6]
-        Copying PartitionStack:
-        (6|3,4,8,7|1,9,5|0,2)
-        Checking for consistency.
-        Everything is consistent.
-        Clearing copy:
-        (0,3,4,8,7,1,9,5,6,2)
-        Splitting point 6 from original:
-        0
-        (6|3,4,8,7|1,9,5|0,2)
-        Splitting point 1 from original:
-        5
-        (6|3,4,8,7|1|5,9|0,2)
-        Splitting point 8 from original:
-        1
-        (6|8|3,4,7|1|5,9|0,2)
-        Splitting point 2 from original:
-        8
-        (6|8|3,4,7|1|5,9|2|0)
-        Getting permutation from PS2->PS:
-        [6, 1, 0, 8, 3, 9, 2, 7, 4, 5]
-        Finding first smallest:
-        Minimal element is 5, bitset is:
-        0000010001
-        Deallocating PartitionStacks.
-        Done.
-
-    """
-    cdef int i, n = sum([len(cell) for cell in partition]), *gamma
-    cdef bitset_t b
-    print "Allocating PartitionStack..."
-    cdef PartitionStack *PS = PS_new(n, 1), *PS2
-    if PS is NULL:
-        print "Allocation failed!"
-        return
-    print "Allocation passed:"
-    PS_print(PS)
-    print "Checking that entries are in order and correct level."
-    good = True
-    for i from 0 <= i < n-1:
-        if not (PS.entries[i] == i and PS.levels[i] == n):
-            print "Failed at i = %d!"%i
-            print PS.entries[i], PS.levels[i], i, n
-            good = False
-    if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1):
-        print "Failed at i = %d!"%(n-1)
-        good = False
-    if not PS.degree == n or not PS.depth == 0:
-        print "Incorrect degree or depth!"
-        good = False
-    if good: print "Everything seems in order, deallocating."
-    PS_dealloc(PS)
-    print "Deallocated."
-    print "Creating PartitionStack from partition %s."%partition
-    PS = PS_from_list(partition)
-    print "PartitionStack's data:"
-    print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n]
-    print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n]
-    print "depth = %d, degree = %d"%(PS.depth,PS.degree)
-    PS_print(PS)
-    print "Checking PS_is_discrete:"
-    print "True" if PS_is_discrete(PS) else "False"
-    print "Checking PS_num_cells:"
-    print PS_num_cells(PS)
-    print "Checking PS_is_mcr, min cell reps are:"
-    L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)]
-    print L
-    print "Checking PS_is_fixed, fixed elements are:"
-    print [PS.entries[l] for l in L if PS_is_fixed(PS, l)]
-    print "Copying PartitionStack:"
-    PS2 = PS_copy(PS)
-    PS_print(PS2)
-    print "Checking for consistency."
-    good = True
-    for i from 0 <= i < n:
-        if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]:
-            print "Failed at i = %d!"%i
-            good = False
-    if PS.degree != PS2.degree or PS.depth != PS2.depth:
-        print "Failure with degree or depth!"
-        good = False
-    if good:
-        print "Everything is consistent."
-    print "Clearing copy:"
-    PS_clear(PS2)
-    PS_print(PS2)
-    for s in splits:
-        print "Splitting point %d from original:"%s
-        print PS_split_point(PS, s)
-        PS_print(PS)
-    print "Getting permutation from PS2->PS:"
-    gamma = <int *> sage_malloc(n * sizeof(int))
-    PS_get_perm_from(PS, PS2, gamma)
-    print [gamma[i] for i from 0 <= i < n]
-    sage_free(gamma)
-    print "Finding first smallest:"
-    bitset_init(b, n)
-    i = PS_first_smallest(PS, b)
-    print "Minimal element is %d, bitset is:"%i
-    print bitset_string(b)
-    bitset_clear(b)
-    print "Deallocating PartitionStacks."
-    PS_dealloc(PS)
-    PS_dealloc(PS2)
-    print "Done."
-
 # Functions
-
-cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S,
-    int (*refine_and_return_invariant)\
-         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
-    int *cells_to_refine_by):
-    """
-    Make the partition stack one longer by copying the last partition in the
-    stack, split off a given point, and refine. Return the invariant given by
-    the refinement function.
-    
-    INPUT:
-    PS -- the partition stack to refine
-    v -- the point to split
-    S -- the structure
-    refine_and_return_invariant -- the refinement function provided
-    cells_to_refine_by -- an array, contents ignored
-
-    """
-    PS.depth += 1
-    PS_clear(PS)
-    cells_to_refine_by[0] = PS_split_point(PS, v)
-    return refine_and_return_invariant(PS, S, cells_to_refine_by, 1)
 
 cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, object S):
     return 0
@@ -694 +108 @@
 cdef int refine_and_return_invariant_trivial(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     return 0
 
-cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S):
+cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S1, object S2):
     return 0
 
 def test_get_aut_gp_and_can_lab_trivially(int n=6, partition=[[0,1,2],[3,4],[5]],
@@ -747 +161 @@
     bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
     int (*refine_and_return_invariant)\
          (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
-    int (*compare_structures)(int *gamma_1, int *gamma_2, object S),
+    int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2),
     bint canonical_label, bint base, bint order):
     """
     Traverse the search space for subgroup/canonical label calculation.
@@ -779 +193 @@
         int -- returns an invariant under application of arbitrary permutations
     compare_structures -- pointer to a function
         INPUT:
-        gamma_1, gamma_2 -- (list) permutations of the points of S
-        S -- pointer to the structure
+        gamma_1, gamma_2 -- (list) permutations of the points of S1 and S2
+        S1, S2 -- pointers to the structures
         OUTPUT:
-        int -- 0 if gamma_1(S) = gamma_2(S), otherwise -1 or 1 (see docs for cmp),
+        int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp),
             such that the set of all structures is well-ordered
     OUTPUT:
     pointer to a aut_gp_and_can_lab_return struct
@@ -1020 +434 @@
     # Ignore the invariant this time, since we are
     # creating the root of the search tree.
     refine_and_return_invariant(current_ps, S, cells_to_refine_by, j)
-
+    PS_move_all_mins_to_front(current_ps)
+    
     # Refine down to a discrete partition
     while not PS_is_discrete(current_ps):
         if not all_children_are_equivalent(current_ps, S):
@@ -1029 +444 @@
         bitset_unset(vertices_have_been_reduced, current_ps.depth)
         i = current_ps.depth
         current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
+        PS_move_all_mins_to_front(current_ps)
         first_indicators[i] = current_indicators[i]
         if canonical_label:
             label_indicators[i] = current_indicators[i]
@@ -1074 +490 @@
                 bitset_flip(vertices_have_been_reduced, current_ps.depth)
             # Look for a new point to split.
             i = vertices_determining_current_stack[current_ps.depth] + 1
-            while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
-                i += 1
-            if i < n:
+            i = bitset_next(vertices_to_split[current_ps.depth], i)
+            if i != -1:
                 # There is a new point.
                 vertices_determining_current_stack[current_ps.depth] = i
                 new_vertex = 1
@@ -1106 +521 @@
                     subgroup_primary_orbit_size += 1
                 # Look for a new point to split.
                 i += 1
-                while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
-                    i += 1
-                if i < n:
+                i = bitset_next(vertices_to_split[current_ps.depth], i)
+                if i != -1:
                     # There is a new point.
                     vertices_determining_current_stack[current_ps.depth] = i
                     if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i:
@@ -1139 +553 @@
             current_kids_are_same = current_ps.depth + 1
         if first_and_current_indicator_same > current_ps.depth:
             first_and_current_indicator_same = current_ps.depth
-        if canonical_label and label_and_current_indicator_same > current_ps.depth:
+        if canonical_label and label_and_current_indicator_same >= current_ps.depth:
             label_and_current_indicator_same = current_ps.depth
             compared_current_and_label_indicators = 0
         
@@ -1149 +563 @@
         while 1:
             i = current_ps.depth
             current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
+            PS_move_all_mins_to_front(current_ps)
             if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]:
                 first_and_current_indicator_same += 1
             if canonical_label:
@@ -1177 +592 @@
         if discrete:
             if current_ps.depth == first_and_current_indicator_same:
                 PS_get_perm_from(current_ps, first_ps, permutation)
-                if compare_structures(permutation, id_perm, S) == 0:
+                if compare_structures(permutation, id_perm, S, S) == 0:
                     automorphism = 1
         if not automorphism:
             if (not canonical_label) or (compared_current_and_label_indicators < 0):
@@ -1188 +603 @@
                 if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth):
                     update_label = 1
                 else:
-                    i = compare_structures(current_ps.entries, label_ps.entries, S)
+                    i = compare_structures(current_ps.entries, label_ps.entries, S, S)
                     if i > 0:
                         update_label = 1
                     elif i < 0:

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi

@@ +1 @@
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include '../../../ext/cdefs.pxi'
+include '../../../ext/stdsage.pxi'
+include '../../../misc/bitset_pxd.pxi'
+
+cdef struct OrbitPartition:
+    # Disjoint set data structure for representing the orbits of the generators
+    # so far found. Also keeps track of the minimum elements of the cells and
+    # their sizes.
+    int degree
+    int *parent
+    int *rank
+    int *mcr # minimum cell representatives - only valid at the root of a cell
+    int *size # also only valid at the root of a cell
+
+cdef struct PartitionStack:
+    # Representation of a node of the search tree. A sequence of partitions of
+    # length depth + 1, each of which is finer than the last. Partition k is
+    # represented as PS.entries in order, broken up immediately after each
+    # entry of levels which is at most k.
+    int *entries
+    int *levels
+    int depth
+    int degree
+
+

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi

@@ +1 @@
+r"""
+Data structures
+
+This module implements TODO...
+
+REFERENCE:
+
+    [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
+        Vol. 30 (1981), pp. 45-87.
+
+"""
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include '../../../misc/bitset.pxi'
+
+cdef inline int min(int a, int b):
+    if a < b: return a
+    else: return b
+
+cdef inline int max(int a, int b):
+    if a > b: return a
+    else: return b
+
+# OrbitPartitions
+
+cdef inline OrbitPartition *OP_new(int n):
+    """
+    Allocate and return a pointer to a new OrbitPartition of degree n. Returns a
+    null pointer in the case of an allocation failure.
+    """
+    cdef int i
+    cdef OrbitPartition *OP = <OrbitPartition *> \
+                                sage_malloc(sizeof(OrbitPartition))
+    if OP is NULL:
+        return OP
+    OP.degree = n
+    OP.parent = <int *> sage_malloc( n * sizeof(int) )
+    OP.rank = <int *> sage_malloc( n * sizeof(int) )
+    OP.mcr = <int *> sage_malloc( n * sizeof(int) )
+    OP.size = <int *> sage_malloc( n * sizeof(int) )
+    if OP.parent is NULL or OP.rank is NULL or \
+       OP.mcr is NULL or OP.size is NULL:
+        if OP.parent is not NULL: sage_free(OP.parent)
+        if OP.rank is not NULL: sage_free(OP.rank)
+        if OP.mcr is not NULL: sage_free(OP.mcr)
+        if OP.size is not NULL: sage_free(OP.size)
+        return NULL
+    for i from 0 <= i < n:
+        OP.parent[i] = i
+        OP.rank[i] = 0
+        OP.mcr[i] = i
+        OP.size[i] = 1
+    return OP
+
+cdef inline int OP_dealloc(OrbitPartition *OP):
+    sage_free(OP.parent)
+    sage_free(OP.rank)
+    sage_free(OP.mcr)
+    sage_free(OP.size)
+    sage_free(OP)
+    return 0
+
+cdef inline int OP_find(OrbitPartition *OP, int n):
+    """
+    Report the representative ("root") of the cell which contains n.
+    """
+    if OP.parent[n] == n:
+        return n
+    else:
+        OP.parent[n] = OP_find(OP, OP.parent[n])
+        return OP.parent[n]
+
+cdef inline int OP_join(OrbitPartition *OP, int m, int n):
+    """
+    Join the cells containing m and n, if they are different.
+    """
+    cdef int m_root = OP_find(OP, m)
+    cdef int n_root = OP_find(OP, n)
+    if OP.rank[m_root] > OP.rank[n_root]:
+        OP.parent[n_root] = m_root
+        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
+        OP.size[m_root] += OP.size[n_root]
+    elif OP.rank[m_root] < OP.rank[n_root]:
+        OP.parent[m_root] = n_root
+        OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root])
+        OP.size[n_root] += OP.size[m_root]
+    elif m_root != n_root:
+        OP.parent[n_root] = m_root
+        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
+        OP.size[m_root] += OP.size[n_root]
+        OP.rank[m_root] += 1
+
+cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma):
+    """
+    Joins the cells of OP which intersect the same orbit of gamma.
+    
+    INPUT:
+        gamma - an integer array representing i -> gamma[i].
+    
+    OUTPUT:
+        1 - something changed
+        0 - orbits of gamma all contained in cells of OP
+    """
+    cdef int i, i_root, gamma_i_root, changed = 0
+    for i from 0 <= i < OP.degree:
+        if gamma[i] == i: continue
+        i_root = OP_find(OP, i)
+        gamma_i_root = OP_find(OP, gamma[i])
+        if i_root != gamma_i_root:
+            changed = 1
+            OP_join(OP, i_root, gamma_i_root)
+    return changed
+
+def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]):
+    """
+    Demonstration and testing.
+    
+    DOCTEST:
+        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
+        ... import OP_represent
+        sage: OP_represent()
+        Allocating OrbitPartition...
+        Allocation passed.
+        Checking that each element reports itself as its root.
+        Each element reports itself as its root.
+        Merging:
+        Merged 0 and 1.
+        Merged 2 and 3.
+        Merged 3 and 4.
+        Done merging.
+        Finding:
+        0 -> 0, root: size=2, mcr=0, rank=1
+        1 -> 0
+        2 -> 2, root: size=3, mcr=2, rank=1
+        3 -> 2
+        4 -> 2
+        5 -> 5, root: size=1, mcr=5, rank=0
+        6 -> 6, root: size=1, mcr=6, rank=0
+        7 -> 7, root: size=1, mcr=7, rank=0
+        8 -> 8, root: size=1, mcr=8, rank=0
+        Allocating array to test merge_perm.
+        Allocation passed.
+        Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8]
+        Done merging.
+        Finding:
+        0 -> 0, root: size=5, mcr=0, rank=2
+        1 -> 0
+        2 -> 0
+        3 -> 0
+        4 -> 0
+        5 -> 5, root: size=3, mcr=5, rank=1
+        6 -> 5
+        7 -> 5
+        8 -> 8, root: size=1, mcr=8, rank=0
+        Deallocating OrbitPartition.
+        Done.
+
+    """
+    cdef int i
+    print "Allocating OrbitPartition..."
+    cdef OrbitPartition *OP = OP_new(n)
+    if OP is NULL:
+        print "Allocation failed!"
+        return
+    print "Allocation passed."
+    print "Checking that each element reports itself as its root."
+    good = True
+    for i from 0 <= i < n:
+        if not OP_find(OP, i) == i:
+            print "Failed at i = %d!"%i
+            good = False
+    if good: print "Each element reports itself as its root."
+    print "Merging:"
+    for i,j in merges:
+        OP_join(OP, i, j)
+        print "Merged %d and %d."%(i,j)
+    print "Done merging."
+    print "Finding:"
+    for i from 0 <= i < n:
+        j = OP_find(OP, i)
+        s = "%d -> %d"%(i, j)
+        if i == j:
+            s += ", root: size=%d, mcr=%d, rank=%d"%\
+                   (OP.size[i], OP.mcr[i], OP.rank[i])
+        print s
+    print "Allocating array to test merge_perm."
+    cdef int *gamma = <int *> sage_malloc( n * sizeof(int) )
+    if gamma is NULL:
+        print "Allocation failed!"
+        OP_dealloc(OP)
+        return
+    print "Allocation passed."
+    for i from 0 <= i < n:
+        gamma[i] = perm[i]
+    print "Merging permutation: %s"%perm
+    OP_merge_list_perm(OP, gamma)
+    print "Done merging."
+    print "Finding:"
+    for i from 0 <= i < n:
+        j = OP_find(OP, i)
+        s = "%d -> %d"%(i, j)
+        if i == j:
+            s += ", root: size=%d, mcr=%d, rank=%d"%\
+                   (OP.size[i], OP.mcr[i], OP.rank[i])
+        print s    
+    print "Deallocating OrbitPartition."
+    sage_free(gamma)
+    OP_dealloc(OP)
+    print "Done."
+
+# PartitionStacks
+
+cdef inline PartitionStack *PS_new(int n, bint unit_partition):
+    """
+    Allocate and return a pointer to a new PartitionStack of degree n. Returns a
+    null pointer in the case of an allocation failure.
+    """
+    cdef int i
+    cdef PartitionStack *PS = <PartitionStack *> \
+                                sage_malloc(sizeof(PartitionStack))
+    if PS is NULL:
+        return PS
+    PS.entries = <int *> sage_malloc( n * sizeof(int) )
+    PS.levels = <int *> sage_malloc( n * sizeof(int) )
+    if PS.entries is NULL or PS.levels is NULL:
+        if PS.entries is not NULL: sage_free(PS.entries)
+        if PS.levels is not NULL: sage_free(PS.levels)
+        return NULL
+    PS.depth = 0
+    PS.degree = n
+    if unit_partition:
+        for i from 0 <= i < n-1:
+            PS.entries[i] = i
+            PS.levels[i] = n
+        PS.entries[n-1] = n-1
+        PS.levels[n-1] = -1
+    return PS
+
+cdef inline PartitionStack *PS_copy(PartitionStack *PS):
+    """
+    Allocate and return a pointer to a copy of PartitionStack PS. Returns a null
+    pointer in the case of an allocation failure.
+    """
+    cdef int i
+    cdef PartitionStack *PS2 = <PartitionStack *> \
+                                sage_malloc(sizeof(PartitionStack))
+    if PS2 is NULL:
+        return PS2
+    PS2.entries = <int *> sage_malloc( PS.degree * sizeof(int) )
+    PS2.levels = <int *> sage_malloc( PS.degree * sizeof(int) )
+    if PS2.entries is NULL or PS2.levels is NULL:
+        if PS2.entries is not NULL: sage_free(PS2.entries)
+        if PS2.levels is not NULL: sage_free(PS2.levels)
+        return NULL
+    PS_copy_from_to(PS, PS2)
+    return PS2
+
+cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2):
+    """
+    Copy all data from PS to PS2.
+    """
+    PS2.depth = PS.depth
+    PS2.degree = PS.degree
+    for i from 0 <= i < PS.degree:
+        PS2.entries[i] = PS.entries[i]
+        PS2.levels[i] = PS.levels[i]    
+
+cdef inline PartitionStack *PS_from_list(object L):
+    """
+    Allocate and return a pointer to a PartitionStack representing L. Returns a
+    null pointer in the case of an allocation failure.
+    """
+    cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0
+    for cell from 0 <= cell < num_cells:
+        n += len(L[cell])
+    cdef PartitionStack *PS = PS_new(n, 0)
+    cell = 0
+    if PS is NULL:
+        return PS
+    while 1:
+        cur_len = len(L[cell])
+        for i from 0 <= i < cur_len:
+            PS.entries[cur_start + i] = L[cell][i]
+            PS.levels[cur_start + i] = n
+        PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1)
+        cur_start += cur_len
+        cell += 1
+        if cell == num_cells:
+            PS.levels[cur_start-1] = -1
+            break
+        PS.levels[cur_start-1] = 0
+    return PS
+
+cdef inline int PS_dealloc(PartitionStack *PS):
+    sage_free(PS.entries)
+    sage_free(PS.levels)
+    sage_free(PS)
+    return 0
+
+cdef inline int PS_print(PartitionStack *PS):
+    """
+    Print a visual representation of PS.
+    """
+    cdef int i
+    for i from 0 <= i < PS.depth + 1:
+        PS_print_partition(PS, i)
+    print
+
+cdef inline int PS_print_partition(PartitionStack *PS, int k):
+    """
+    Print the partition at depth k.
+    """
+    s = '('
+    for i from 0 <= i < PS.degree:
+        s += str(PS.entries[i])
+        if PS.levels[i] <= k:
+            s += '|'
+        else:
+            s += ','
+    s = s[:-1] + ')'
+    print s
+
+cdef inline bint PS_is_discrete(PartitionStack *PS):
+    """
+    Returns whether the deepest partition consists only of singleton cells.
+    """
+    cdef int i
+    for i from 0 <= i < PS.degree:
+        if PS.levels[i] > PS.depth:
+            return 0
+    return 1
+
+cdef inline int PS_num_cells(PartitionStack *PS):
+    """
+    Returns the number of cells.
+    """
+    cdef int i, ncells = 0
+    for i from 0 <= i < PS.degree:
+        if PS.levels[i] <= PS.depth:
+            ncells += 1
+    return ncells
+
+cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end):
+    """
+    Makes sure that the first element of the segment of entries i with
+    start <= i <= end is minimal.
+    """
+    cdef int i, min_loc = start, minimum = PS.entries[start]
+    for i from start < i <= end:
+        if PS.entries[i] < minimum:
+            min_loc = i
+            minimum = PS.entries[i]
+    if min_loc != start:
+        PS.entries[min_loc] = PS.entries[start]
+        PS.entries[start] = minimum
+
+cdef inline bint PS_is_mcr(PartitionStack *PS, int m):
+    """
+    Returns whether PS.elements[m] (not m!) is the smallest element of its cell.
+    """
+    return m == 0 or PS.levels[m-1] <= PS.depth
+
+cdef inline bint PS_is_fixed(PartitionStack *PS, int m):
+    """
+    Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming
+    PS_is_mcr(PS, m) is already true.
+    """
+    return PS.levels[m] <= PS.depth
+
+cdef inline int PS_clear(PartitionStack *PS):
+    """
+    Sets the current partition to the first shallower one, i.e. forgets about
+    boundaries between cells that are new to the current level.
+    """
+    cdef int i, cur_start = 0
+    for i from 0 <= i < PS.degree:
+        if PS.levels[i] >= PS.depth:
+            PS.levels[i] += 1
+        if PS.levels[i] < PS.depth:
+            PS_move_min_to_front(PS, cur_start, i)
+            cur_start = i+1
+
+cdef inline int PS_move_all_mins_to_front(PartitionStack *PS):
+    """
+    Move minimal cell elements to the front of each cell.
+    """
+    cdef int i, cur_start = 0
+    for i from 0 <= i < PS.degree:
+        if PS.levels[i] <= PS.depth:
+            PS_move_min_to_front(PS, cur_start, i)
+            cur_start = i+1
+
+cdef inline int PS_split_point(PartitionStack *PS, int v):
+    """
+    Detaches the point v from the cell it is in, putting the singleton cell
+    of just v in front. Returns the position where v is now located.
+    """
+    cdef int i = 0, index_of_v
+    while PS.entries[i] != v:
+        i += 1
+    index_of_v = i
+    while PS.levels[i] > PS.depth:
+        i += 1
+    if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \
+       and PS.levels[index_of_v] > PS.depth:
+        # if v is first (make sure v is not already alone)
+        PS_move_min_to_front(PS, index_of_v+1, i)
+        PS.levels[index_of_v] = PS.depth
+        return index_of_v
+    else:
+        # If v is not at front, i.e. v is not minimal in its cell,
+        # then move_min_to_front is not necessary since v will swap
+        # with the first before making its own cell, leaving it at
+        # the front of the other.
+        i = index_of_v
+        while i != 0 and PS.levels[i-1] > PS.depth:
+            i -= 1
+        PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v
+        PS.entries[i+1] = PS.entries[i] # move the first (min) to second
+        PS.entries[i] = v # place v first
+        PS.levels[i] = PS.depth
+        return i
+
+cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b):
+    """
+    Find the first occurrence of the smallest cell of size greater than one,
+    store its entries to b, and return its minimum element.
+    """
+    cdef int i = 0, j = 0, location = 0, n = PS.degree
+    bitset_zero(b)
+    while 1:
+        if PS.levels[i] <= PS.depth:
+            if i != j and n > i - j + 1:
+                n = i - j + 1
+                location = j
+            j = i + 1
+        if PS.levels[i] == -1: break
+        i += 1
+    # location now points to the beginning of the first, smallest,
+    # nontrivial cell
+    i = location
+    while 1:
+        bitset_flip(b, PS.entries[i])
+        if PS.levels[i] <= PS.depth: break
+        i += 1
+    return PS.entries[location]
+
+cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma):
+    """
+    Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i]
+    denotes the entry of the ith cell of the discrete partition PS.
+    """
+    cdef int i
+    for i from 0 <= i < PS1.degree:
+        gamma[PS2.entries[i]] = PS1.entries[i]
+
+def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]):
+    """
+    Demonstration and testing.
+    
+    DOCTEST:
+        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
+        ... import PS_represent
+        sage: PS_represent()
+        Allocating PartitionStack...
+        Allocation passed:
+        (0,1,2,3,4,5,6,7,8,9)
+        Checking that entries are in order and correct level.
+        Everything seems in order, deallocating.
+        Deallocated.
+        Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]].
+        PartitionStack's data:
+        entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2]
+        levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1]
+        depth = 0, degree = 10
+        (6|3,4,8,7|1,9,5|0,2)
+        Checking PS_is_discrete:
+        False
+        Checking PS_num_cells:
+        4
+        Checking PS_is_mcr, min cell reps are:
+        [6, 3, 1, 0]
+        Checking PS_is_fixed, fixed elements are:
+        [6]
+        Copying PartitionStack:
+        (6|3,4,8,7|1,9,5|0,2)
+        Checking for consistency.
+        Everything is consistent.
+        Clearing copy:
+        (0,3,4,8,7,1,9,5,6,2)
+        Splitting point 6 from original:
+        0
+        (6|3,4,8,7|1,9,5|0,2)
+        Splitting point 1 from original:
+        5
+        (6|3,4,8,7|1|5,9|0,2)
+        Splitting point 8 from original:
+        1
+        (6|8|3,4,7|1|5,9|0,2)
+        Splitting point 2 from original:
+        8
+        (6|8|3,4,7|1|5,9|2|0)
+        Getting permutation from PS2->PS:
+        [6, 1, 0, 8, 3, 9, 2, 7, 4, 5]
+        Finding first smallest:
+        Minimal element is 5, bitset is:
+        0000010001
+        Deallocating PartitionStacks.
+        Done.
+
+    """
+    cdef int i, n = sum([len(cell) for cell in partition]), *gamma
+    cdef bitset_t b
+    print "Allocating PartitionStack..."
+    cdef PartitionStack *PS = PS_new(n, 1), *PS2
+    if PS is NULL:
+        print "Allocation failed!"
+        return
+    print "Allocation passed:"
+    PS_print(PS)
+    print "Checking that entries are in order and correct level."
+    good = True
+    for i from 0 <= i < n-1:
+        if not (PS.entries[i] == i and PS.levels[i] == n):
+            print "Failed at i = %d!"%i
+            print PS.entries[i], PS.levels[i], i, n
+            good = False
+    if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1):
+        print "Failed at i = %d!"%(n-1)
+        good = False
+    if not PS.degree == n or not PS.depth == 0:
+        print "Incorrect degree or depth!"
+        good = False
+    if good: print "Everything seems in order, deallocating."
+    PS_dealloc(PS)
+    print "Deallocated."
+    print "Creating PartitionStack from partition %s."%partition
+    PS = PS_from_list(partition)
+    print "PartitionStack's data:"
+    print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n]
+    print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n]
+    print "depth = %d, degree = %d"%(PS.depth,PS.degree)
+    PS_print(PS)
+    print "Checking PS_is_discrete:"
+    print "True" if PS_is_discrete(PS) else "False"
+    print "Checking PS_num_cells:"
+    print PS_num_cells(PS)
+    print "Checking PS_is_mcr, min cell reps are:"
+    L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)]
+    print L
+    print "Checking PS_is_fixed, fixed elements are:"
+    print [PS.entries[l] for l in L if PS_is_fixed(PS, l)]
+    print "Copying PartitionStack:"
+    PS2 = PS_copy(PS)
+    PS_print(PS2)
+    print "Checking for consistency."
+    good = True
+    for i from 0 <= i < n:
+        if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]:
+            print "Failed at i = %d!"%i
+            good = False
+    if PS.degree != PS2.degree or PS.depth != PS2.depth:
+        print "Failure with degree or depth!"
+        good = False
+    if good:
+        print "Everything is consistent."
+    print "Clearing copy:"
+    PS_clear(PS2)
+    PS_print(PS2)
+    for s in splits:
+        print "Splitting point %d from original:"%s
+        print PS_split_point(PS, s)
+        PS_print(PS)
+    print "Getting permutation from PS2->PS:"
+    gamma = <int *> sage_malloc(n * sizeof(int))
+    PS_get_perm_from(PS, PS2, gamma)
+    print [gamma[i] for i from 0 <= i < n]
+    sage_free(gamma)
+    print "Finding first smallest:"
+    bitset_init(b, n)
+    i = PS_first_smallest(PS, b)
+    print "Minimal element is %d, bitset is:"%i
+    print bitset_string(b)
+    bitset_clear(b)
+    print "Deallocating PartitionStacks."
+    PS_dealloc(PS)
+    PS_dealloc(PS2)
+    print "Done."
+
+# Functions
+
+cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S,
+    int (*refine_and_return_invariant)\
+         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
+    int *cells_to_refine_by):
+    """
+    Make the partition stack one longer by copying the last partition in the
+    stack, split off a given point, and refine. Return the invariant given by
+    the refinement function.
+    
+    INPUT:
+    PS -- the partition stack to refine
+    v -- the point to split
+    S -- the structure
+    refine_and_return_invariant -- the refinement function provided
+    cells_to_refine_by -- an array, contents ignored
+
+    """
+    PS.depth += 1
+    PS_clear(PS)
+    cells_to_refine_by[0] = PS_split_point(PS, v)
+    return refine_and_return_invariant(PS, S, cells_to_refine_by, 1)
+

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/double_coset.pxd

@@ +1 @@
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include '../../../ext/cdefs.pxi'
+include '../../../ext/stdsage.pxi'
+include 'data_structures_pxd.pxi' # includes bitsets
+
+from sage.rings.integer cimport Integer
+
+cdef int *double_coset( object, object, int **, int *, int,
+    bint (*)(PartitionStack *, object),
+    int (*)(PartitionStack *, object, int *, int),
+    int (*)(int *, int *, object, object))
+
+
+

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/double_coset.pyx

@@ +1 @@
+r"""
+Double cosets
+
+This module implements a general algorithm for computing double coset problems
+for pairs of objects. The class of objects in question must be some kind
+of structure for which an isomorphism is a permutation in $S_n$ for some $n$,
+which we call here the order of the object. Given objects $X$ and $Y$,
+the program returns an isomorphism in list permutation form if $X \cong Y$, and
+a NULL pointer otherwise.
+
+In order to take advantage of the algorithms in this module for a specific kind
+of object, one must implement (in Cython) three functions which will be specific
+to the kind of objects in question. Pointers to these functions are passed to
+the main function of the module, which is \code{double_coset}. For specific
+examples of implementations of these functions, see any of the files in
+\code{sage.groups.perm_gps.partn_ref} beginning with "refinement." They are:
+
+A. \code{refine_and_return_invariant}:
+    
+    Signature:
+    
+    \code{int refine_and_return_invariant(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len)}
+    
+    This function should split up cells in the partition at the top of the
+    partition stack in such a way that any automorphism that respects the
+    partition also respects the resulting partition. The array
+    cells_to_refine_by is a list of the beginning positions of some cells which
+    have been changed since the last refinement. It is not necessary to use
+    this in an implementation of this function, but it will affect performance.
+    One should consult \code{refinement_graphs} for more details and ideas for
+    particular implementations.
+    
+    Output:
+    
+    An integer $I$ invariant under the orbits of $S_n$.  That is, if
+    $\gamma \in S_n$, then
+    $$ I(G, PS, cells_to_refine_by) = I( \gamma(G), \gamma(PS), \gamma(cells_to_refine_by) ) .$$
+
+
+B. \code{compare_structures}:
+    
+    Signature:
+    
+    \code{int compare_structures(int *gamma_1, int *gamma_2, object S)}
+    
+    This function must implement a total ordering on the set of objects of fixed
+    order. Return -1 if \code{gamma_1(S) < gamma_2(S)}, 0 if
+    \code{gamma_1(S) == gamma_2(S)}, 1 if \code{gamma_1(S) > gamma_2(S)}.
+
+C. \code{all_children_are_equivalent}:
+    
+    Signature:
+    
+    \code{bint all_children_are_equivalent(PartitionStack *PS, object S)}
+    
+    This function must return False unless it is the case that each discrete
+    partition finer than the top of the partition stack is equivalent to the
+    others under some automorphism of S. The converse need not hold: if this is
+    indeed the case, it still may return False. This function is originally used
+    as a consequence of Lemma 2.25 in [1].
+
+DOCTEST:
+    sage: import sage.groups.perm_gps.partn_ref.double_coset
+
+REFERENCE:
+
+    [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
+        Vol. 30 (1981), pp. 45-87.
+
+    [2] Leon, Jeffrey. Permutation Group Algorithms Based on Partitions, I:
+        Theory and Algorithms. J. Symbolic Computation, Vol. 12 (1991), pp.
+        533-583.
+
+"""
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include 'data_structures_pyx.pxi' # includes bitsets
+
+cdef inline int cmp(int a, int b):
+    if a < b: return -1
+    elif a == b: return 0
+    else: return 1
+
+cdef inline bint stacks_are_equivalent(PartitionStack *PS1, PartitionStack *PS2):
+    cdef int i, j, depth = min(PS1.depth, PS2.depth)
+    for i from 0 <= i < PS1.degree:
+        j = cmp(PS1.levels[i], PS2.levels[i])
+        if j == 0: continue
+        if ( (j < 0 and PS1.levels[i] <= depth and PS2.levels[i] > depth)
+            or (j > 0 and PS2.levels[i] <= depth and PS1.levels[i] > depth) ):
+            return 0
+    return 1
+
+# Functions
+
+cdef int *double_coset(object S1, object S2, int **partition1, int *ordering2, 
+    int n, bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
+    int (*refine_and_return_invariant)\
+         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
+    int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2) ):
+    """
+    Traverse the search space for double coset calculation.
+    
+    INPUT:
+    S1, S2 -- pointers to the structures
+    partition1 -- array representing a partition of the points of S1
+    ordering2 -- an ordering of the points of S2 representing a second partition
+    n -- the number of points (points are assumed to be 0,1,...,n-1)
+    all_children_are_equivalent -- pointer to a function
+        INPUT:
+        PS -- pointer to a partition stack
+        S -- pointer to the structure
+        OUTPUT:
+        bint -- returns True if it can be determined that all refinements below
+            the current one will result in an equivalent discrete partition
+    refine_and_return_invariant -- pointer to a function
+        INPUT:
+        PS -- pointer to a partition stack
+        S -- pointer to the structure
+        alpha -- an array consisting of numbers, which indicate the starting
+            positions of the cells to refine against (will likely be modified)
+        OUTPUT:
+        int -- returns an invariant under application of arbitrary permutations
+    compare_structures -- pointer to a function
+        INPUT:
+        gamma_1, gamma_2 -- (list) permutations of the points of S1 and S2
+        S1, S2 -- pointers to the structures
+        OUTPUT:
+        int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp),
+            such that the set of all structures is well-ordered
+    OUTPUT:
+    If S1 and S2 are isomorphic, a pointer to an integer array representing an
+    isomorphism. Otherwise, a NULL pointer.
+
+    """
+    cdef PartitionStack *current_ps, *first_ps, *left_ps
+    cdef int first_meets_current = -1
+    cdef int current_kids_are_same = 1
+    cdef int first_kids_are_same
+    
+    cdef int *indicators
+    
+    cdef OrbitPartition *orbits_of_subgroup
+    cdef int subgroup_primary_orbit_size = 0
+    cdef int minimal_in_primary_orbit
+    
+    cdef bitset_t *fixed_points_of_generators # i.e. fp
+    cdef bitset_t *minimal_cell_reps_of_generators # i.e. mcr
+    cdef int len_of_fp_and_mcr = 100
+    cdef int index_in_fp_and_mcr = -1
+    
+    cdef bitset_t *vertices_to_split
+    cdef bitset_t vertices_have_been_reduced
+    cdef int *permutation, *id_perm, *cells_to_refine_by
+    cdef int *vertices_determining_current_stack
+    
+    cdef int *output
+    
+    cdef int i, j, k
+    cdef bint discrete, automorphism, update_label
+    cdef bint backtrack, new_vertex, narrow, mem_err = 0
+    
+    if n == 0:
+        return NULL
+    
+    indicators = <int *> sage_malloc(n * sizeof(int))
+    
+    fixed_points_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
+    minimal_cell_reps_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
+    
+    vertices_to_split = <bitset_t *> sage_malloc( n * sizeof(bitset_t) )
+    permutation = <int *> sage_malloc( n * sizeof(int) )
+    id_perm = <int *> sage_malloc( n * sizeof(int) )
+    cells_to_refine_by = <int *> sage_malloc( n * sizeof(int) )
+    vertices_determining_current_stack = <int *> sage_malloc( n * sizeof(int) )
+
+    output = <int *> sage_malloc( n * sizeof(int) )
+    
+    current_ps = PS_new(n, 0)
+    left_ps = PS_new(n, 0)
+    orbits_of_subgroup = OP_new(n)
+    
+    # Check for allocation failures:
+    if indicators                         is NULL or \
+       fixed_points_of_generators         is NULL or \
+       minimal_cell_reps_of_generators    is NULL or \
+       vertices_to_split                  is NULL or \
+       permutation                        is NULL or \
+       id_perm                            is NULL or \
+       cells_to_refine_by                 is NULL or \
+       vertices_determining_current_stack is NULL or \
+       current_ps                         is NULL or \
+       orbits_of_subgroup                 is NULL or \
+       output                             is NULL:
+        if indicators                         is not NULL:
+            sage_free(indicators)
+        if fixed_points_of_generators         is not NULL:
+            sage_free(fixed_points_of_generators)
+        if minimal_cell_reps_of_generators    is not NULL:
+            sage_free(minimal_cell_reps_of_generators)
+        if vertices_to_split                  is not NULL:
+            sage_free(vertices_to_split)
+        if permutation                        is not NULL:
+            sage_free(permutation)
+        if id_perm                            is not NULL:
+            sage_free(id_perm)
+        if cells_to_refine_by                 is not NULL:
+            sage_free(cells_to_refine_by)
+        if vertices_determining_current_stack is not NULL:
+            sage_free(vertices_determining_current_stack)
+        if current_ps is not NULL:
+            PS_dealloc(current_ps)
+        if orbits_of_subgroup is not NULL:
+            OP_dealloc(orbits_of_subgroup)
+        if output is not NULL:
+            sage_free(output)
+        raise MemoryError
+    
+    # Initialize bitsets, checking for allocation failures:
+    cdef bint succeeded = 1
+    for i from 0 <= i < len_of_fp_and_mcr:
+        try:
+            bitset_init(fixed_points_of_generators[i], n)
+        except MemoryError:
+            succeeded = 0
+            for j from 0 <= j < i:
+                bitset_clear(fixed_points_of_generators[j])
+                bitset_clear(minimal_cell_reps_of_generators[j])
+            break
+        try:
+            bitset_init(minimal_cell_reps_of_generators[i], n)
+        except MemoryError:
+            succeeded = 0
+            for j from 0 <= j < i:
+                bitset_clear(fixed_points_of_generators[j])
+                bitset_clear(minimal_cell_reps_of_generators[j])
+            bitset_clear(fixed_points_of_generators[i])
+            break
+    if succeeded:
+        for i from 0 <= i < n:
+            try:
+                bitset_init(vertices_to_split[i], n)
+            except MemoryError:
+                succeeded = 0
+                for j from 0 <= j < i:
+                    bitset_clear(vertices_to_split[j])
+                for j from 0 <= j < len_of_fp_and_mcr:
+                    bitset_clear(fixed_points_of_generators[j])
+                    bitset_clear(minimal_cell_reps_of_generators[j])
+                break
+    if succeeded:
+        try:
+            bitset_init(vertices_have_been_reduced, n)
+        except MemoryError:
+            succeeded = 0
+            for j from 0 <= j < n:
+                bitset_clear(vertices_to_split[j])
+            for j from 0 <= j < len_of_fp_and_mcr:
+                bitset_clear(fixed_points_of_generators[j])
+                bitset_clear(minimal_cell_reps_of_generators[j])
+    if not succeeded:
+        sage_free(indicators)
+        sage_free(fixed_points_of_generators)
+        sage_free(minimal_cell_reps_of_generators)
+        sage_free(vertices_to_split)
+        sage_free(permutation)
+        sage_free(id_perm)
+        sage_free(cells_to_refine_by)
+        sage_free(vertices_determining_current_stack)
+        PS_dealloc(current_ps)
+        PS_dealloc(left_ps)
+        OP_dealloc(orbits_of_subgroup)
+        raise MemoryError
+    
+    bitset_zero(vertices_have_been_reduced)
+    
+    # set up the identity permutation
+    for i from 0 <= i < n:
+        id_perm[i] = i
+    if ordering2 is NULL:
+        ordering2 = id_perm
+    
+    # Copy data of partition to left_ps, and
+    # ordering of that partition to current_ps.
+    i = 0
+    j = 0
+    while partition1[i] is not NULL:
+        k = 0
+        while partition1[i][k] != -1:
+            left_ps.entries[j+k] = partition1[i][k]
+            left_ps.levels[j+k] = n
+            current_ps.entries[j+k] = ordering2[j+k]
+            current_ps.levels[j+k] = n
+            k += 1
+        left_ps.levels[j+k-1] = 0
+        current_ps.levels[j+k-1] = 0
+        PS_move_min_to_front(current_ps, j, j+k-1)
+        j += k
+        i += 1
+    left_ps.levels[j-1] = -1
+    current_ps.levels[j-1] = -1
+    left_ps.depth = 0
+    current_ps.depth = 0
+    left_ps.degree = n
+    current_ps.degree = n
+    
+    # default values of "infinity"
+    for i from 0 <= i < n:
+        indicators[i] = -1
+    
+    cdef bint possible = 1
+    cdef bint unknown = 1
+    
+    # Our first refinement needs to check every cell of the partition,
+    # so cells_to_refine_by needs to be a list of pointers to each cell.
+    j = 1
+    cells_to_refine_by[0] = 0
+    for i from 0 < i < n:
+        if left_ps.levels[i-1] == 0:
+            cells_to_refine_by[j] = i
+            j += 1
+    
+    k = refine_and_return_invariant(left_ps, S1, cells_to_refine_by, j)
+    
+    j = 1
+    cells_to_refine_by[0] = 0
+    for i from 0 < i < n:
+        if current_ps.levels[i-1] == 0:
+            cells_to_refine_by[j] = i
+            j += 1
+    j = refine_and_return_invariant(current_ps, S2, cells_to_refine_by, j)
+    if k != j:
+        possible = 0; unknown = 0
+    elif not stacks_are_equivalent(left_ps, current_ps):
+        possible = 0; unknown = 0
+    else:
+        PS_move_all_mins_to_front(current_ps)
+
+    first_ps = NULL
+    # Refine down to a discrete partition
+    while not PS_is_discrete(left_ps) and possible:
+        k = PS_first_smallest(left_ps, vertices_to_split[left_ps.depth])
+        i = left_ps.depth
+        indicators[i] = split_point_and_refine(left_ps, k, S1, refine_and_return_invariant, cells_to_refine_by)
+        vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
+        bitset_unset(vertices_have_been_reduced, current_ps.depth)
+        while 1:
+            j =  split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by)
+            if indicators[i] != j:
+                possible = 0
+            elif not stacks_are_equivalent(left_ps, current_ps):
+                possible = 0
+            else:
+                PS_move_all_mins_to_front(current_ps)
+            if not possible:
+                j = vertices_determining_current_stack[i] + 1
+                j = bitset_next(vertices_to_split[i], j)
+                if j == -1:
+                    break
+                else:
+                    possible = 1
+                    vertices_determining_current_stack[i] = j
+                    current_ps.depth -= 1 # reset for next refinement
+            else: break
+        if not all_children_are_equivalent(current_ps, S2):
+            current_kids_are_same = current_ps.depth + 1
+    
+    if possible:
+        if compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
+            unknown = 0
+        else:
+            first_meets_current = current_ps.depth
+            first_kids_are_same = current_ps.depth
+            first_ps = PS_copy(current_ps)
+            current_ps.depth -= 1
+   
+    # Main loop:
+    while possible and unknown and current_ps.depth != -1:
+        
+        # I. Search for a new vertex to split, and update subgroup information
+        new_vertex = 0
+        if current_ps.depth > first_meets_current:
+            # If we are not at a node of the first stack, reduce size of
+            # vertices_to_split by using the symmetries we already know.
+            if not bitset_check(vertices_have_been_reduced, current_ps.depth):
+                for i from 0 <= i <= index_in_fp_and_mcr:
+                    j = 0
+                    while j < current_ps.depth and bitset_check(fixed_points_of_generators[i], vertices_determining_current_stack[j]):
+                        j += 1
+                    # If each vertex split so far is fixed by generator i,
+                    # then remove elements of vertices_to_split which are
+                    # not minimal in their orbits under generator i.
+                    if j == current_ps.depth:
+                        for k from 0 <= k < n:
+                            if bitset_check(vertices_to_split[current_ps.depth], k) and not bitset_check(minimal_cell_reps_of_generators[i], k):
+                                bitset_flip(vertices_to_split[current_ps.depth], k)
+                bitset_flip(vertices_have_been_reduced, current_ps.depth)
+            # Look for a new point to split.
+            i = vertices_determining_current_stack[current_ps.depth] + 1
+            i = bitset_next(vertices_to_split[current_ps.depth], i)
+            if i != -1:
+                # There is a new point.
+                vertices_determining_current_stack[current_ps.depth] = i
+                new_vertex = 1
+            else:
+                # No new point: backtrack.
+                current_ps.depth -= 1
+        else:
+            # If we are at a node of the first stack, the above reduction
+            # will not help. Also, we must update information about
+            # primary orbits here.
+            if current_ps.depth < first_meets_current:
+                # If we are done searching under this part of the first
+                # stack, then first_meets_current is one higher, and we
+                # are looking at a new primary orbit (corresponding to a
+                # larger subgroup in the stabilizer chain).
+                first_meets_current = current_ps.depth
+                for i from 0 <= i < n:
+                    if bitset_check(vertices_to_split[current_ps.depth], i):
+                        minimal_in_primary_orbit = i
+                        break
+            while 1:
+                i = vertices_determining_current_stack[current_ps.depth]
+                # This was the last point to be split here.
+                # If it is in the same orbit as minimal_in_primary_orbit,
+                # then count it as an element of the primary orbit.
+                if OP_find(orbits_of_subgroup, i) == OP_find(orbits_of_subgroup, minimal_in_primary_orbit):
+                    subgroup_primary_orbit_size += 1
+                # Look for a new point to split.
+                i += 1
+                i = bitset_next(vertices_to_split[current_ps.depth], i)
+                if i != -1:
+                    # There is a new point.
+                    vertices_determining_current_stack[current_ps.depth] = i
+                    if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i:
+                        new_vertex = 1
+                        break
+                else:
+                    # No new point: backtrack.
+                    # Note that now, we are backtracking up the first stack.
+                    vertices_determining_current_stack[current_ps.depth] = -1
+                    # If every choice of point to split gave something in the
+                    # (same!) primary orbit, then all children of the first
+                    # stack at this point are equivalent.
+                    j = 0
+                    for i from 0 <= i < n:
+                        if bitset_check(vertices_to_split[current_ps.depth], i):
+                            j += 1
+                    if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1:
+                        first_kids_are_same = current_ps.depth
+                    # Backtrack.
+                    subgroup_primary_orbit_size = 0
+                    current_ps.depth -= 1
+                    break
+        if not new_vertex:
+            continue
+    
+        if current_kids_are_same > current_ps.depth + 1:
+            current_kids_are_same = current_ps.depth + 1
+        
+        # II. Refine down to a discrete partition, or until
+        # we leave the part of the tree we are interested in
+        discrete = 0
+        while 1:
+            i = current_ps.depth
+            while 1:
+                k = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by)
+                PS_move_all_mins_to_front(current_ps)
+                if indicators[i] != k:
+                    possible = 0
+                elif not stacks_are_equivalent(left_ps, current_ps):
+                    possible = 0
+                if not possible:
+                    j = vertices_determining_current_stack[i] + 1
+                    j = bitset_next(vertices_to_split[i], j)
+                    if j == -1:
+                        break
+                    else:
+                        possible = 1
+                        vertices_determining_current_stack[i] = j
+                        current_ps.depth -= 1 # reset for next refinement
+                else: break
+            if not possible:
+                break
+            if PS_is_discrete(current_ps):
+                break
+            vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
+            bitset_unset(vertices_have_been_reduced, current_ps.depth)
+            if not all_children_are_equivalent(current_ps, S2):
+                current_kids_are_same = current_ps.depth + 1
+        
+        # III. Check for automorphisms and isomorphisms
+        automorphism = 0
+        if possible:
+            PS_get_perm_from(current_ps, first_ps, permutation)
+            if compare_structures(permutation, id_perm, S2, S2) == 0:
+                automorphism = 1
+        if not automorphism and possible:
+            # if we get here, discrete must be true
+            if current_ps.depth != left_ps.depth:
+                possible = 0
+            elif compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
+                unknown = 0
+                break # main loop
+            else:
+                possible = 0
+        if automorphism:
+            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
+                index_in_fp_and_mcr += 1
+            bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr])
+            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
+            for i from 0 <= i < n:
+                if permutation[i] == i:
+                    bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
+                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+                else:
+                    bitset_unset(fixed_points_of_generators[index_in_fp_and_mcr], i)
+                    k = i
+                    j = permutation[i]
+                    while j != i:
+                        if j < k: k = j
+                        j = permutation[j]
+                    if k == i:
+                        bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+                    else:
+                        bitset_unset(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+            current_ps.depth = first_meets_current
+            if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser
+                if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit:
+                    continue # main loop
+            if bitset_check(vertices_have_been_reduced, current_ps.depth):
+                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
+            continue # main loop
+        if not possible:
+            possible = 1
+            i = current_ps.depth
+            current_ps.depth = min(first_kids_are_same-1, current_kids_are_same-1)
+            if i == current_kids_are_same:
+                continue # main loop
+            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
+                index_in_fp_and_mcr += 1
+            bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr])
+            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
+            j = current_ps.depth
+            current_ps.depth = i # just for mcr and fixed functions...
+            for i from 0 <= i < n:
+                if PS_is_mcr(current_ps, i):
+                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+                    if PS_is_fixed(current_ps, i):
+                        bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
+            current_ps.depth = j
+            if bitset_check(vertices_have_been_reduced, current_ps.depth):
+                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
+    
+    # End of main loop.
+    
+    if possible and not unknown:
+        for i from 0 <= i < n:
+            output[left_ps.entries[i]] = current_ps.entries[i]
+    else:
+        sage_free(output)
+        output = NULL
+    
+    # Deallocate:
+    for i from 0 <= i < len_of_fp_and_mcr:
+        bitset_clear(fixed_points_of_generators[i])
+        bitset_clear(minimal_cell_reps_of_generators[i])
+    for i from 0 <= i < n:
+        bitset_clear(vertices_to_split[i])
+    bitset_clear(vertices_have_been_reduced)
+    sage_free(indicators)
+    sage_free(fixed_points_of_generators)
+    sage_free(minimal_cell_reps_of_generators)
+    sage_free(vertices_to_split)
+    sage_free(permutation)
+    sage_free(id_perm)
+    sage_free(cells_to_refine_by)
+    sage_free(vertices_determining_current_stack)
+    PS_dealloc(current_ps)
+    PS_dealloc(left_ps)
+    OP_dealloc(orbits_of_subgroup)
+    if first_ps is not NULL: PS_dealloc(first_ps)
+    
+    return output
+
+
+
+
+

sage/groups/perm_gps/partn_ref/refinement_binary.pxd

@@ -8 +8 @@
 
 include '../../../ext/cdefs.pxi'
 include '../../../ext/stdsage.pxi'
-include '../../../misc/bitset_pxd.pxi'
+include 'data_structures_pxd.pxi' # includes bitsets
 
 from sage.rings.integer cimport Integer
-from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport PartitionStack, PS_print, PS_new, PS_dealloc, PS_clear, PS_is_discrete, PS_move_min_to_front, PS_num_cells, get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
+from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
+from double_coset cimport double_coset
 
 cdef class BinaryCodeStruct:
     cdef bitset_s *alpha_is_wd # single bitset of length nwords + degree
@@ -36 +37 @@
 cdef int ith_word_nonlinear(BinaryCodeStruct, int, bitset_s *)
 
 cdef int refine_by_bip_degree(PartitionStack *, object, int *, int)
-cdef int compare_linear_codes(int *, int *, object)
-cdef int compare_nonlinear_codes(int *, int *, object)
+cdef int compare_linear_codes(int *, int *, object, object)
+cdef int compare_nonlinear_codes(int *, int *, object, object)
 cdef bint all_children_are_equivalent(PartitionStack *, object)
 cdef inline int word_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *)
 cdef inline int col_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *)

sage/groups/perm_gps/partn_ref/refinement_binary.pyx

@@ -22 +22 @@
 #                         http://www.gnu.org/licenses/
 #*****************************************************************************
 
-include '../../../misc/bitset.pxi'
+include 'data_structures_pyx.pxi' # includes bitsets
+
 from sage.matrix.matrix import is_Matrix
 
 cdef class LinearBinaryCodeStruct(BinaryCodeStruct):
@@ -95 +96 @@
         
         self.output = NULL
         self.ith_word = &ith_word_linear
-        self.first_time = 1
     
     def run(self, partition=None):
         """
@@ -238 +238 @@
                 part[i][j] = partition[i][j]
             part[i][len(partition[i])] = -1
         part[len(partition)] = NULL
+        self.first_time = 1
         
         self.output = get_aut_gp_and_can_lab(self, part, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes, 1, 1, 1)
         
@@ -295 +296 @@
         if self.output is NULL:
             self.run()
         return [self.output.relabeling[i] for i from 0 <= i < self.degree]
+    
+    def is_isomorphic(self, LinearBinaryCodeStruct other):
+        """
+        Calculate whether self is isomorphic to other.
+        
+        EXAMPLES:
+            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import LinearBinaryCodeStruct
+
+            sage: B = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]))
+            sage: C = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0,0,1],[1,1,0,1,1,0]]))
+            sage: B.is_isomorphic(C)
+            [0, 1, 2, 5, 3, 4]
+            
+        """
+        cdef int **part, i, j
+        cdef int *output, *ordering
+        partition = [range(self.degree)]
+        part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
+        ordering = <int *> sage_malloc(self.degree * sizeof(int))
+        if part is NULL or ordering is NULL:
+            if part is not NULL: sage_free(part)
+            if ordering is not NULL: sage_free(ordering)
+            raise MemoryError
+        for i from 0 <= i < len(partition):
+            part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
+            if part[i] is NULL:
+                for j from 0 <= j < i:
+                    sage_free(part[j])
+                sage_free(part)
+                raise MemoryError
+            for j from 0 <= j < len(partition[i]):
+                part[i][j] = partition[i][j]
+            part[i][len(partition[i])] = -1
+        part[len(partition)] = NULL
+        for i from 0 <= i < self.degree:
+            ordering[i] = i
+        self.first_time = 1
+        other.first_time = 1
+        
+        output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes)
+        
+        for i from 0 <= i < len(partition):
+            sage_free(part[i])
+        sage_free(part)
+        sage_free(ordering)
+
+        if output is NULL:
+            return False
+        else:
+            output_py = [output[i] for i from 0 <= i < self.degree]
+            sage_free(output)
+            return output_py
     
     def __dealloc__(self):
         cdef int j
@@ -402 +455 @@
         
         self.output = NULL
         self.ith_word = &ith_word_nonlinear
-        self.first_time = 1
 
     def __dealloc__(self):
         cdef int j
@@ -477 +529 @@
                 part[i][j] = partition[i][j]
             part[i][len(partition[i])] = -1
         part[len(partition)] = NULL
+        self.first_time = 1
         
         self.output = get_aut_gp_and_can_lab(self, part, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes, 1, 1, 1)
         
@@ -535 +588 @@
             self.run()
         return [self.output.relabeling[i] for i from 0 <= i < self.degree]
 
+    def is_isomorphic(self, NonlinearBinaryCodeStruct other):
+        """
+        Calculate whether self is isomorphic to other.
+        
+        EXAMPLES:
+            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
+
+            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]))
+            sage: C = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,0,0,1,1],[1,1,1,1,0,0]]))
+            sage: B.is_isomorphic(C)
+            [2, 3, 0, 1, 4, 5]
+            
+        """
+        cdef int **part, i, j
+        cdef int *output, *ordering
+        partition = [range(self.degree)]
+        part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
+        ordering = <int *> sage_malloc(self.degree * sizeof(int))
+        if part is NULL or ordering is NULL:
+            if part is not NULL: sage_free(part)
+            if ordering is not NULL: sage_free(ordering)
+            raise MemoryError
+        for i from 0 <= i < len(partition):
+            part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
+            if part[i] is NULL:
+                for j from 0 <= j < i:
+                    sage_free(part[j])
+                sage_free(part)
+                raise MemoryError
+            for j from 0 <= j < len(partition[i]):
+                part[i][j] = partition[i][j]
+            part[i][len(partition[i])] = -1
+        part[len(partition)] = NULL
+        for i from 0 <= i < self.degree:
+            ordering[i] = i
+        self.first_time = 1
+        other.first_time = 1
+        
+        output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes)
+        
+        for i from 0 <= i < len(partition):
+            sage_free(part[i])
+        sage_free(part)
+        sage_free(ordering)
+
+        if output is NULL:
+            return False
+        else:
+            output_py = [output[i] for i from 0 <= i < self.degree]
+            sage_free(output)
+            return output_py
+    
 cdef int ith_word_nonlinear(BinaryCodeStruct self, int i, bitset_s *word):
     cdef NonlinearBinaryCodeStruct NBCS = <NonlinearBinaryCodeStruct> self
     bitset_copy(word, &NBCS.words[i])
@@ -606 +711 @@
                 if necessary_to_split_cell:
                     invariant += 8
                     first_largest_subcell = sort_by_function(col_ps, current_cell, col_degrees, col_counts, col_output, BCS.nwords+1)
+                    invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps)
                     invariant += first_largest_subcell
                     against_index = current_cell_against
                     while against_index < ctrb_len:
@@ -622 +728 @@
                         r += 1
                         if r >= i:
                             break
-                    invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps)
                     invariant += (i - current_cell)
                 current_cell = i
         else:
@@ -641 +746 @@
                 if necessary_to_split_cell:
                     invariant += 64
                     first_largest_subcell = sort_by_function(word_ps, current_cell, word_degrees, word_counts, word_output, BCS.degree+1)
+                    invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps)
                     invariant += first_largest_subcell
                     against_index = current_cell_against
                     while against_index < ctrb_len:
@@ -658 +764 @@
                         r += 1
                         if r >= i:
                             break
-                    invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps)
                     invariant += (i - current_cell)
                 current_cell = i
         current_cell_against += 1
     return invariant
 
-cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S):
+cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S1, object S2):
     """
-    Compare gamma_1(B) and gamma_2(B).
+    Compare gamma_1(S1) and gamma_2(S2).
     
-    Return return -1 if gamma_1(B) < gamma_2(B), 0 if gamma_1(B) == gamma_2(B),
-    1 if gamma_1(B) > gamma_2(B).  (Just like the python \code{cmp}) function.
+    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
+    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
     
     Abstractly, what this function does is relabel the basis of B by gamma_1 and
     gamma_2, run a row reduction on each, and verify that the matrices are the
@@ -680 +785 @@
     
     INPUT:
     gamma_1, gamma_2 -- list permutations (inverse)
-    S -- a binary code struct object
+    S1, S2 -- binary code struct objects
     
     """
     cdef int i, piv_loc_1, piv_loc_2, cur_col, cur_row=0
     cdef bint is_pivot_1, is_pivot_2
-    cdef LinearBinaryCodeStruct BCS = <LinearBinaryCodeStruct> S
-    cdef bitset_s *basis_1 = BCS.scratch_bitsets                   # len = dim
-    cdef bitset_s *basis_2 = &BCS.scratch_bitsets[BCS.dimension]   # len = dim
-    cdef bitset_s *pivots  = &BCS.scratch_bitsets[2*BCS.dimension] # len 1
-    cdef bitset_s *temp  = &BCS.scratch_bitsets[2*BCS.dimension+1] # len 1
-    for i from 0 <= i < BCS.dimension:
-        bitset_copy(&basis_1[i], &BCS.basis[i])
-        bitset_copy(&basis_2[i], &BCS.basis[i])
+    cdef LinearBinaryCodeStruct BCS1 = <LinearBinaryCodeStruct> S1
+    cdef LinearBinaryCodeStruct BCS2 = <LinearBinaryCodeStruct> S2
+    cdef bitset_s *basis_1 = BCS1.scratch_bitsets                   # len = dim
+    cdef bitset_s *basis_2 = &BCS1.scratch_bitsets[BCS1.dimension]   # len = dim
+    cdef bitset_s *pivots  = &BCS1.scratch_bitsets[2*BCS1.dimension] # len 1
+    cdef bitset_s *temp  = &BCS1.scratch_bitsets[2*BCS1.dimension+1] # len 1
+    for i from 0 <= i < BCS1.dimension:
+        bitset_copy(&basis_1[i], &BCS1.basis[i])
+        bitset_copy(&basis_2[i], &BCS2.basis[i])
     bitset_zero(pivots)
-    for cur_col from 0 <= cur_col < BCS.degree:
+    for cur_col from 0 <= cur_col < BCS1.degree:
         is_pivot_1 = 0
         is_pivot_2 = 0
-        for i from cur_row <= i < BCS.dimension:
+        for i from cur_row <= i < BCS1.dimension:
             if bitset_check(&basis_1[i], gamma_1[cur_col]):
                 is_pivot_1 = 1
                 piv_loc_1 = i
                 break
-        for i from cur_row <= i < BCS.dimension:
+        for i from cur_row <= i < BCS1.dimension:
             if bitset_check(&basis_2[i], gamma_2[cur_col]):
                 is_pivot_2 = 1
                 piv_loc_2 = i
@@ -724 +830 @@
                     bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
                 if bitset_check(&basis_2[i], gamma_2[cur_col]):
                     bitset_xor(&basis_2[i], &basis_2[i], &basis_2[cur_row])
-            for i from cur_row < i < BCS.dimension:
+            for i from cur_row < i < BCS1.dimension:
                 if bitset_check(&basis_1[i], gamma_1[cur_col]):
                     bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
                 if bitset_check(&basis_2[i], gamma_2[cur_col]):
@@ -736 +842 @@
                     return <int>bitset_check(&basis_2[i], gamma_2[cur_col]) - <int>bitset_check(&basis_1[i], gamma_1[cur_col])
     return 0
 
-cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S):
+cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S1, object S2):
     """
-    Compare gamma_1(B) and gamma_2(B).
+    Compare gamma_1(S1) and gamma_2(S2).
     
-    Return return -1 if gamma_1(B) < gamma_2(B), 0 if gamma_1(B) == gamma_2(B),
-    1 if gamma_1(B) > gamma_2(B).  (Just like the python \code{cmp}) function.
+    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
+    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
     
     INPUT:
     gamma_1, gamma_2 -- list permutations (inverse)
-    S -- a binary code struct object
+    S1, S2 -- a binary code struct object
     
     """
     cdef int side=0, i, start, end, n_one_1, n_one_2, cur_col
     cdef int where_0, where_1
-    cdef NonlinearBinaryCodeStruct BCS = <NonlinearBinaryCodeStruct> S
-    cdef bitset_s *B_1_0 = BCS.scratch_bitsets                   # nwords of len degree
-    cdef bitset_s *B_1_1 = &BCS.scratch_bitsets[BCS.nwords]      # nwords of len degree
-    cdef bitset_s *B_2_0 = &BCS.scratch_bitsets[2*BCS.nwords]    # nwords of len degree
-    cdef bitset_s *B_2_1 = &BCS.scratch_bitsets[3*BCS.nwords]    # nwords of len degree
-    cdef bitset_s *dividers = &BCS.scratch_bitsets[4*BCS.nwords] # 1 of len nwords
+    cdef NonlinearBinaryCodeStruct BCS1 = <NonlinearBinaryCodeStruct> S1
+    cdef NonlinearBinaryCodeStruct BCS2 = <NonlinearBinaryCodeStruct> S2
+    cdef bitset_s *B_1_0 = BCS1.scratch_bitsets                   # nwords of len degree
+    cdef bitset_s *B_1_1 = &BCS1.scratch_bitsets[BCS1.nwords]      # nwords of len degree
+    cdef bitset_s *B_2_0 = &BCS1.scratch_bitsets[2*BCS1.nwords]    # nwords of len degree
+    cdef bitset_s *B_2_1 = &BCS1.scratch_bitsets[3*BCS1.nwords]    # nwords of len degree
+    cdef bitset_s *dividers = &BCS1.scratch_bitsets[4*BCS1.nwords] # 1 of len nwords
     cdef bitset_s *B_1_this, *B_1_other, *B_2_this, *B_2_other
-    for i from 0 <= i < BCS.nwords:
-        bitset_copy(&B_1_0[i], &BCS.words[i])
-        bitset_copy(&B_2_0[i], &BCS.words[i])
+    for i from 0 <= i < BCS1.nwords:
+        bitset_copy(&B_1_0[i], &BCS1.words[i])
+        bitset_copy(&B_2_0[i], &BCS2.words[i])
     bitset_zero(dividers)
-    bitset_set(dividers, BCS.nwords-1)
+    bitset_set(dividers, BCS1.nwords-1)
     
-    for cur_col from 0 <= cur_col < BCS.degree:
+    for cur_col from 0 <= cur_col < BCS1.degree:
         if side == 0:
             B_1_this  = B_1_0
             B_1_other = B_1_1
@@ -776 +883 @@
             B_2_other = B_2_0
         side ^= 1
         start = 0
-        while start < BCS.nwords:
+        while start < BCS1.nwords:
             end = start
             while not bitset_check(dividers, end):
                 end += 1
@@ -960 +1067 @@
 def random_tests(t=10.0, n_max=50, k_max=6, nwords_max=200, perms_per_code=10, density_range=(.1,.9)):
     """
     Tests to make sure that C(gamma(B)) == C(B) for random permutations gamma
-    and random codes B.
+    and random codes B, and that is_isomorphic returns an isomorphism.
 
     INPUT: 
     t -- run tests for approximately this many seconds
@@ -977 +1084 @@
 
     For each code B (B2) generated, we uniformly generate perms_per_code random 
     permutations and verify that the canonical labels of B and the image of B
-    under the generated permutation are equal.
+    under the generated permutation are equal, and check that the double coset
+    function returns an isomorphism.
     
     DOCTEST:
         sage: import sage.groups.perm_gps.partn_ref.refinement_binary
@@ -992 +1100 @@
     from sage.combinat.permutation import Permutations
     from sage.matrix.constructor import random_matrix, matrix
     from sage.rings.finite_field import FiniteField as GF
-    cdef int h, i, j, n, k, num_tests = 0, num_codes = 0, passed = 1
+    cdef int h, i, j, n, k, num_tests = 0, num_codes = 0
     cdef LinearBinaryCodeStruct B, C
     cdef NonlinearBinaryCodeStruct B_n, C_n
     t_0 = walltime()
@@ -1040 +1148 @@
                     B_n_M[j,B_n_relab[h]] = bitset_check(&B_n.words[j], h)
                     C_n_M[j,C_n_relab[h]] = bitset_check(&C_n.words[j], h)
             if B_M.row_space() != C_M.row_space():
-                print "B:"
+                print "can_lab error -- B:"
                 for j from 0 <= j < B.dimension:
                     print bitset_string(&B.basis[j])
                 print perm
                 return
             if sorted(B_n_M.rows()) != sorted(C_n_M.rows()):
-                print "B_n:"
+                print "can_lab error -- B_n:"
                 for j from 0 <= j < B_n.nwords:
                     print bitset_string(&B_n.words[j])
                 print perm
                 return
+            isom = B.is_isomorphic(C)
+            if not isom:
+                print "isom -- B:"
+                for j from 0 <= j < B.dimension:
+                    print bitset_string(&B.basis[j])
+                print perm
+                print isom
+                return
+            isom = B_n.is_isomorphic(C_n)
+            if not isom:
+                print "isom -- B_n:"
+                for j from 0 <= j < B_n.nwords:
+                    print bitset_string(&B_n.words[j])
+                print perm
+                print isom
+                return
                 
-        num_tests += perms_per_code
+        num_tests += 4*perms_per_code
         num_codes += 2
-    if passed:
-        print "All passed: %d random tests on %d codes."%(num_tests, num_codes)
+        
+    print "All passed: %d random tests on %d codes."%(num_tests, num_codes)
 
 

sage/groups/perm_gps/partn_ref/refinement_graphs.pxd

@@ -8 +8 @@
 
 include '../../../ext/cdefs.pxi'
 include '../../../ext/stdsage.pxi'
+include 'data_structures_pxd.pxi' # includes bitsets
 
 from sage.graphs.base.c_graph cimport CGraph
 from sage.graphs.base.sparse_graph cimport SparseGraph
 from sage.graphs.base.dense_graph cimport DenseGraph
 from sage.rings.integer cimport Integer
-from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport PartitionStack, PS_is_discrete, PS_move_min_to_front, PS_num_cells, get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
+from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
+from double_coset cimport double_coset
 
 cdef class GraphStruct:
     cdef CGraph G
@@ -22 +24 @@
     cdef int *scratch # length 3n+1
 
 cdef int refine_by_degree(PartitionStack *, object, int *, int)
-cdef int compare_graphs(int *, int *, object)
+cdef int compare_graphs(int *, int *, object, object)
 cdef bint all_children_are_equivalent(PartitionStack *, object)
 cdef inline int degree(PartitionStack *, CGraph, int, int, bint)
 cdef inline int sort_by_function(PartitionStack *, int, int *)

sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

@@ -17 +17 @@
 # Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
 #                         http://www.gnu.org/licenses/
 #*****************************************************************************
+
+include 'data_structures_pyx.pxi' # includes bitsets
+
+def isomorphic(G1, G2, partition, ordering2, dig, use_indicator_function, sparse=False):
+    """
+    Tests whether two graphs are isomorphic.
+    
+    sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import isomorphic
+
+    sage: G = Graph(2)
+    sage: H = Graph(2)
+    sage: isomorphic(G, H, [[0,1]], [0,1], 0, 1)
+    [0, 1]
+    sage: isomorphic(G, H, [[0,1]], [0,1], 0, 1)
+    [0, 1]
+    sage: isomorphic(G, H, [[0],[1]], [0,1], 0, 1)
+    [0, 1]
+    sage: isomorphic(G, H, [[0],[1]], [1,0], 0, 1)
+    [1, 0]
+
+    sage: G = Graph(3)
+    sage: H = Graph(3)
+    sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1)
+    [0, 1, 2]
+    sage: G.add_edge(0,1)
+    sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1)
+    False
+    sage: H.add_edge(1,2)
+    sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1)
+    [1, 2, 0]
+
+    """
+    cdef int **part
+    cdef int *output, *ordering
+    cdef CGraph G
+    cdef GraphStruct GS1 = GraphStruct()
+    cdef GraphStruct GS2 = GraphStruct()
+    cdef GraphStruct GS
+    cdef int i, j, n = -1
+
+    from sage.graphs.graph import GenericGraph, Graph, DiGraph
+    for G_in in [G1, G2]:
+        if G_in is G1:
+            GS = GS1
+        else:
+            GS = GS2
+        if isinstance(G_in, GenericGraph):
+            if n == -1:
+                n = G_in.num_verts()
+            elif n != G_in.num_verts():
+                return False
+            G_in = G_in.copy()
+            if G_in.vertices() != range(n):
+                to = G_in.relabel(return_map=True)
+                frm = {}
+                for v in to.iterkeys():
+                    frm[to[v]] = v
+                partition = [[to[v] for v in cell] for cell in partition]
+            else:
+                to = range(n)
+                frm = to
+            if sparse:
+                G = SparseGraph(n)
+            else:
+                G = DenseGraph(n)
+            if G_in.is_directed():
+                for i from 0 <= i < n:
+                    for _,j,_ in G_in.outgoing_edge_iterator(i):
+                        G.add_arc(i,j)
+            else:
+                for i from 0 <= i < n:
+                    for _,j,_ in G_in.edge_iterator(i):
+                        if j <= i:
+                            G.add_arc(i,j)
+                            G.add_arc(j,i)
+        elif isinstance(G_in, CGraph):
+            G = <CGraph> G_in
+            if n == -1:
+                n = G.num_verts
+            elif n != G.num_verts:
+                return False
+            to = range(n)
+            frm = to
+        else:
+            raise TypeError("G must be a Sage graph.")
+        GS.G = G
+        GS.directed = 1 if dig else 0
+        GS.use_indicator = 1 if use_indicator_function else 0
+
+    if n == 0:
+        return {}
+    
+    part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
+    ordering = <int *> sage_malloc(n * sizeof(int))
+    if part is NULL or ordering is NULL:
+        if part is not NULL: sage_free(part)
+        if ordering is not NULL: sage_free(ordering)
+        raise MemoryError
+    for i from 0 <= i < len(partition):
+        part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
+        if part[i] is NULL:
+            for j from 0 <= j < i:
+                sage_free(part[j])
+            sage_free(part)
+            raise MemoryError
+        for j from 0 <= j < len(partition[i]):
+            part[i][j] = partition[i][j]
+        part[i][len(partition[i])] = -1
+    part[len(partition)] = NULL
+    for i from 0 <= i < n:
+        ordering[i] = ordering2[i]
+
+    GS1.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
+    GS2.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
+    if GS1.scratch is NULL or GS2.scratch is NULL:
+        if GS1.scratch is not NULL: sage_free(GS1.scratch)
+        if GS2.scratch is not NULL: sage_free(GS2.scratch)
+        for j from 0 <= j < len(partition):
+            sage_free(part[j])
+        sage_free(part)
+        raise MemoryError
+
+    output = double_coset(GS1, GS2, part, ordering, n, &all_children_are_equivalent, &refine_by_degree, &compare_graphs)
+
+    for i from 0 <= i < len(partition):
+        sage_free(part[i])
+    sage_free(part)
+    sage_free(ordering)
+    sage_free(GS1.scratch)
+    sage_free(GS2.scratch)
+    
+    if output is NULL:
+        return False
+    else:
+        output_py = [output[i] for i from 0 <= i < n]
+        sage_free(output)
+        # TODO: figure out frm, to stuff to relabel for consistency with input
+        return output_py
 
 def search_tree(G_in, partition, lab=True, dig=False, dict_rep=False, certify=False,
                     verbosity=0, use_indicator_function=True, sparse=False,
@@ -200 +338 @@
         sage: st(G, [range(G.num_verts())])[1] == st(H, [range(H.num_verts())])[1]
         True
 
+        sage: from sage.graphs.graph import graph_isom_equivalent_non_multi_graph
+        sage: G = Graph(multiedges=True, implementation='networkx')
+        sage: G.add_edge(('a', 'b'))
+        sage: G.add_edge(('a', 'b'))
+        sage: G.add_edge(('a', 'b'))
+        sage: G, Pi = graph_isom_equivalent_non_multi_graph(G, [['a','b']])
+        sage: s,b = st(G, Pi, lab=False, dict_rep=True)
+        sage: sorted(b.items())
+        [(('o', 'a'), 5), (('o', 'b'), 1), (('x', 0), 2), (('x', 1), 3), (('x', 2), 4)]
+
+        sage: st(Graph(':Dkw'), [range(5)], lab=False, dig=True)
+        [[4, 1, 2, 3, 0], [0, 2, 1, 3, 4]]
+
     """
     cdef CGraph G
     cdef int i, j, n
@@ -277 +428 @@
     if dict_rep:
         ddd = {}
         for v in frm.iterkeys():
-            if frm[v] != 0:
-                ddd[v] = frm[v]
-            else:
-                ddd[v] = n
+            ddd[frm[v]] = v if v != 0 else n
         return_tuple.append(ddd)
     if lab:
         if isinstance(G_in, GenericGraph):
@@ -348 +496 @@
     cdef int current_cell, i, r
     cdef int first_largest_subcell
     cdef int invariant = 1
+    cdef int max_degree
     cdef int *degrees = GS.scratch # length 3n+1
     cdef bint necessary_to_split_cell
     cdef int against_index
@@ -358 +507 @@
             invariant += 50
             i = current_cell
             necessary_to_split_cell = 0
+            max_degree = 0
             while 1:
                 degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 0)
                 if degrees[i-current_cell] != degrees[0]:
                     necessary_to_split_cell = 1
+                if degrees[i-current_cell] > max_degree:
+                    max_degree = degrees[i-current_cell]
                 i += 1
                 if PS.levels[i-1] <= PS.depth:
                     break
@@ -369 +521 @@
             if necessary_to_split_cell:
                 invariant += 10
                 first_largest_subcell = sort_by_function(PS, current_cell, degrees)
-                invariant += first_largest_subcell
+                invariant += first_largest_subcell + max_degree
                 against_index = current_cell_against
                 while against_index < ctrb_len:
                     if cells_to_refine_by[against_index] == current_cell:
@@ -385 +537 @@
                     r += 1
                     if r >= i:
                         break
-                invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 0)
                 invariant += (i - current_cell)
             current_cell = i
         if GS.directed:
@@ -396 +547 @@
                 invariant += 20
                 i = current_cell
                 necessary_to_split_cell = 0
+                max_degree = 0
                 while 1:
                     degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 1)
                     if degrees[i-current_cell] != degrees[0]:
                         necessary_to_split_cell = 1
+                    if degrees[i-current_cell] > max_degree:
+                        max_degree = degrees[i-current_cell]
                     i += 1
                     if PS.levels[i-1] <= PS.depth:
                         break
@@ -407 +561 @@
                 if necessary_to_split_cell:
                     invariant += 7
                     first_largest_subcell = sort_by_function(PS, current_cell, degrees)
-                    invariant += first_largest_subcell
+                    invariant += first_largest_subcell + max_degree
                     against_index = current_cell_against
                     while against_index < ctrb_len:
                         if cells_to_refine_by[against_index] == current_cell:
@@ -425 +579 @@
                         r += 1
                         if r >= i:
                             break
-                    invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 1)
                     invariant += (i - current_cell)
                 current_cell = i
         current_cell_against += 1
@@ -434 +587 @@
     else:
         return 0
     
-cdef int compare_graphs(int *gamma_1, int *gamma_2, object S):
+cdef int compare_graphs(int *gamma_1, int *gamma_2, object S1, object S2):
     r"""
-    Compare gamma_1(G) and gamma_2(G).
+    Compare gamma_1(S1) and gamma_2(S2).
 
-    Return return -1 if gamma_1(G) < gamma_2(G), 0 if gamma_1(G) ==
-    gamma_2(G), 1 if gamma_1(G) > gamma_2(G).  (Just like the python
+    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) ==
+    gamma_2(S2), 1 if gamma_1(S1) > gamma_2(S2).  (Just like the python
     \code{cmp}) function.
     
     INPUT:
     gamma_1, gamma_2 -- list permutations (inverse)
-    S -- a graph struct object
+    S1, S2 -- graph struct objects
     
     """
     cdef int i, j
-    cdef GraphStruct GS = <GraphStruct> S
-    cdef CGraph G = GS.G
-    for i from 0 <= i < G.num_verts:
-        for j from 0 <= j < G.num_verts:
-            if G.has_arc_unsafe(gamma_1[i], gamma_1[j]):
-                if not G.has_arc_unsafe(gamma_2[i], gamma_2[j]):
+    cdef GraphStruct GS1 = <GraphStruct> S1
+    cdef GraphStruct GS2 = <GraphStruct> S2
+    cdef CGraph G1 = GS1.G
+    cdef CGraph G2 = GS2.G
+    for i from 0 <= i < G1.num_verts:
+        for j from 0 <= j < G1.num_verts:
+            if G1.has_arc_unsafe(gamma_1[i], gamma_1[j]):
+                if not G2.has_arc_unsafe(gamma_2[i], gamma_2[j]):
                     return 1
-            elif G.has_arc_unsafe(gamma_2[i], gamma_2[j]):
+            elif G2.has_arc_unsafe(gamma_2[i], gamma_2[j]):
                 return -1
     return 0
 
@@ -472 +627 @@
     S -- a graph struct object
     """
     cdef GraphStruct GS = <GraphStruct> S
+    if GS.directed:
+        return 0
     cdef CGraph G = GS.G
     cdef int i, n = PS.degree
     cdef bint in_cell = 0
@@ -622 +779 @@
 def random_tests(t=10.0, n_max=60, perms_per_graph=10):
     """
     Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma
-    and random graphs G.
+    and random graphs G, and that isomorphic returns an isomorphism.
 
     INPUT:
     t -- run tests for approximately this many seconds
@@ -637 +794 @@
 
     For each graph G generated, we uniformly generate perms_per_graph random
     permutations and verify that the canonical labels of G and the image of G
-    under the generated permutation are equal.
+    under the generated permutation are equal, and that the isomorphic function
+    returns an isomorphism.
         
     TESTS:
         sage: import sage.groups.perm_gps.partn_ref.refinement_graphs
@@ -651 +809 @@
     from sage.misc.prandom import random, randint
     from sage.graphs.graph_generators import GraphGenerators, DiGraphGenerators
     from sage.combinat.permutation import Permutations
-    cdef int i, j, num_tests = 0, num_graphs = 0, passed = 1
+    cdef int i, j, num_tests = 0, num_graphs = 0
     GG = GraphGenerators()
     DGG = DiGraphGenerators()
     t_0 = walltime()
@@ -663 +821 @@
         G = GG.RandomGNP(n, p)
         H = G.copy()
         for i from 0 <= i < perms_per_graph:
+            G = H.copy()
             G1 = search_tree(G, [G.vertices()])[1]
             perm = list(S.random_element())
             perm = [perm[j]-1 for j from 0 <= j < n]
             G.relabel(perm)
             G2 = search_tree(G, [G.vertices()])[1]
             if G1 != G2:
-                print "FAILURE: graph6-"
+                print "search_tree FAILURE: graph6-"
                 print H.graph6_string()
                 print perm
-                passed = 0
+                return
+            isom = isomorphic(G, H, [range(n)], range(n), 0, 1)
+            if not isom or G.relabel(isom, inplace=False) != H:
+                print "isom FAILURE: graph6-"
+                print H.graph6_string()
+                print perm
+                return
 
         D = DGG.RandomDirectedGNP(n, p)
         D.loops(True)
@@ -681 +846 @@
                 D.add_edge(i,i)
         E = D.copy()
         for i from 0 <= i < perms_per_graph:
+            D = E.copy()
             D1 = search_tree(D, [D.vertices()], dig=True)[1]
             perm = list(S.random_element())
             perm = [perm[j]-1 for j from 0 <= j < n]
             D.relabel(perm)
             D2 = search_tree(D, [D.vertices()], dig=True)[1]
             if D1 != D2:
-                print "FAILURE: dig6-"
+                print "search_tree FAILURE: dig6-"
                 print E.dig6_string()
                 print perm
-                passed = 0
-        num_tests += 2*perms_per_graph
+                return
+            isom = isomorphic(D, E, [range(n)], range(n), 1, 1)
+            if not isom or D.relabel(isom, inplace=False) != E:
+                print "isom FAILURE: dig6-"
+                print E.dig6_string()
+                print perm
+                print isom
+                return
+        num_tests += 4*perms_per_graph
         num_graphs += 2
-    if passed:
-        print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs)
+    print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs)
 
 

sage/groups/perm_gps/partn_ref/refinement_matrices.pxd

@@ -8 +8 @@
 
 include '../../../ext/cdefs.pxi'
 include '../../../ext/stdsage.pxi'
-include '../../../misc/bitset_pxd.pxi'
+include 'data_structures_pxd.pxi' # includes bitsets
 
 from sage.rings.integer cimport Integer
-from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport PartitionStack, PS_new, PS_dealloc, PS_copy_from_to, get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
-from sage.groups.perm_gps.partn_ref.refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree
-from sage.groups.perm_gps.partn_ref.refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent
+from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
+from refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree
+from refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent
+from double_coset cimport double_coset
 
 cdef class MatrixStruct:
     cdef list symbol_structs
@@ -26 +27 @@
     cdef aut_gp_and_can_lab_return *output
 
 cdef int refine_matrix(PartitionStack *, object, int *, int)
-cdef int compare_matrices(int *, int *, object)
+cdef int compare_matrices(int *, int *, object, object)