Ticket #10804: trac_10804.patch

File trac_10804.patch, 191.3 KB (added by rlm, 11 years ago)
  • module_list.py

    # HG changeset patch
    # User Robert L. Miller <rlm@rlmiller.org>
    # Date 1279700962 -7200
    # Node ID 8ddd7727b63998b1237c9429fdc1d9b95d85bb3d
    # Parent  872cd84f620c8419cc1dd3e4ad4e8fef92f29442
    #10804 - Stabilizer chains and other improvements to partn_ref
    
    diff -r 872cd84f620c -r 8ddd7727b639 module_list.py
    a b  
    435435   
    436436    Extension('sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label',
    437437              sources = ['sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx'],
    438               libraries = ['gmp']),
     438              libraries = ['gmp', 'flint'],
     439              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     440              extra_compile_args = ['-std=c99'],
     441              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    439442
    440443    Extension('sage.groups.perm_gps.partn_ref.double_coset',
    441               sources = ['sage/groups/perm_gps/partn_ref/double_coset.pyx']),
     444              sources = ['sage/groups/perm_gps/partn_ref/double_coset.pyx'],
     445              libraries = ['gmp', 'flint'],
     446              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     447              extra_compile_args = ['-std=c99'],
     448              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    442449
    443450    Extension('sage.groups.perm_gps.partn_ref.refinement_binary',
    444451              sources = ['sage/groups/perm_gps/partn_ref/refinement_binary.pyx'],
    445               libraries = ['gmp']),
     452              libraries = ['gmp', 'flint'],
     453              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     454              extra_compile_args = ['-std=c99'],
     455              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    446456
    447457    Extension('sage.groups.perm_gps.partn_ref.refinement_graphs',
    448458              sources = ['sage/groups/perm_gps/partn_ref/refinement_graphs.pyx'],
    449               libraries = ['gmp']),
     459              libraries = ['gmp', 'flint'],
     460              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     461              extra_compile_args = ['-std=c99'],
     462              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    450463
    451464    Extension('sage.groups.perm_gps.partn_ref.refinement_lists',
    452465              sources = ['sage/groups/perm_gps/partn_ref/refinement_lists.pyx'],
    453               libraries = ['gmp']),
     466              libraries = ['gmp', 'flint'],
     467              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     468              extra_compile_args = ['-std=c99'],
     469              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    454470
    455471    Extension('sage.groups.perm_gps.partn_ref.refinement_matrices',
    456472              sources = ['sage/groups/perm_gps/partn_ref/refinement_matrices.pyx'],
    457               libraries = ['gmp']),
     473              libraries = ['gmp', 'flint'],
     474              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     475              extra_compile_args = ['-std=c99'],
     476              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    458477
    459478    Extension('sage.groups.perm_gps.partn_ref.refinement_python',
    460479              sources = ['sage/groups/perm_gps/partn_ref/refinement_python.pyx'],
    461               libraries = ['gmp']),
     480              libraries = ['gmp', 'flint'],
     481              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     482              extra_compile_args = ['-std=c99'],
     483              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    462484
    463485    ################################
    464486    ##
     
    15701592    ################################
    15711593
    15721594    Extension('sage.sets.disjoint_set',
    1573               sources = ['sage/sets/disjoint_set.pyx']),
     1595              sources = ['sage/sets/disjoint_set.pyx'],
     1596              libraries = ['gmp', 'flint'],
     1597              include_dirs = [SAGE_ROOT + '/local/include/FLINT/'],
     1598              extra_compile_args = ['-std=c99'],
     1599              depends = [SAGE_ROOT + "/local/include/FLINT/flint.h"]),
    15741600
    15751601    ################################
    15761602    ##
  • sage/graphs/matchpoly.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/graphs/matchpoly.pyx
    a b  
    266266   
    267267    # Computing the signless matching polynomial
    268268
    269     _sig_on
     269    sig_on()
    270270    delete_and_add(edges, nverts, nedges, nverts, 0, pol)
    271     _sig_off
     271    sig_off()
    272272
    273273    # Building the actual matching polynomial
    274274   
  • sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#
    55# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    66#                         http://www.gnu.org/licenses/
     
    1212
    1313from sage.rings.integer cimport Integer
    1414
    15 cdef struct aut_gp_and_can_lab_return:
     15cdef struct aut_gp_and_can_lab:
    1616    int *generators
    1717    int num_gens
     18    StabilizerChain *group
    1819    int *relabeling
    19     int *base
    20     int base_size
    21     mpz_t order
    2220
    23 cdef aut_gp_and_can_lab_return *get_aut_gp_and_can_lab( object, int **, int,
    24     bint (*)(PartitionStack *, object),
    25     int (*)(PartitionStack *, object, int *, int),
    26     int (*)(int *, int *, object, object), bint, bint, bint)
     21cdef aut_gp_and_can_lab *get_aut_gp_and_can_lab( void *,
     22    PartitionStack *, int,
     23    bint (*)(PartitionStack *, void *),
     24    int (*)(PartitionStack *, void *, int *, int),
     25    int (*)(int *, int *, void *, void *), bint, StabilizerChain *) except NULL
    2726
    2827
    2928
  • sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx
    a b  
    3636   
    3737    Signature:
    3838   
    39     \code{int refine_and_return_invariant(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len)}
     39    \code{int refine_and_return_invariant(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len)}
    4040   
    4141    This function should split up cells in the partition at the top of the
    4242    partition stack in such a way that any automorphism that respects the
     
    5858   
    5959    Signature:
    6060   
    61     \code{int compare_structures(int *gamma_1, int *gamma_2, object S1, object S2)}
     61    \code{int compare_structures(int *gamma_1, int *gamma_2, void *S1, void *S2)}
    6262   
    6363    This function must implement a total ordering on the set of objects of fixed
    6464    order. Return:
     
    7777   
    7878    Signature:
    7979   
    80     \code{bint all_children_are_equivalent(PartitionStack *PS, object S)}
     80    \code{bint all_children_are_equivalent(PartitionStack *PS, void *S)}
    8181   
    8282    This function must return False unless it is the case that each discrete
    8383    partition finer than the top of the partition stack is equivalent to the
     
    9696"""
    9797
    9898#*****************************************************************************
    99 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     99#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    100100#
    101101# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    102102#                         http://www.gnu.org/licenses/
    103103#*****************************************************************************
    104104
    105105include 'data_structures_pyx.pxi' # includes bitsets
     106include 'sage/ext/interrupt.pxi'
    106107
    107108cdef inline int cmp(int a, int b):
    108109    if a < b: return -1
     
    111112
    112113# Functions
    113114
    114 cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, object S):
     115cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, void *S):
    115116    return 0
    116117
    117 cdef int refine_and_return_invariant_trivial(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     118cdef int refine_and_return_invariant_trivial(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len):
    118119    return 0
    119120
    120 cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S1, object S2):
     121cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, void *S1, void *S2):
    121122    return 0
    122123
    123 def test_get_aut_gp_and_can_lab_trivially(int n=6, partition=[[0,1,2],[3,4],[5]],
    124     canonical_label=True, base=False):
     124def test_get_aut_gp_and_can_lab_trivially(int n=6,
     125    list partition=[[0,1,2],[3,4],[5]], canonical_label=True, base=False):
    125126    """
    126127    sage: tttt = sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label.test_get_aut_gp_and_can_lab_trivially
    127128    sage: tttt()
     
    142143    1
    143144   
    144145    """
    145     cdef int i, j
    146     cdef aut_gp_and_can_lab_return *output
     146    cdef aut_gp_and_can_lab *output
    147147    cdef Integer I = Integer(0)
    148     cdef int **part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
    149     for i from 0 <= i < len(partition):
    150         part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    151         for j from 0 <= j < len(partition[i]):
    152             part[i][j] = partition[i][j]
    153         part[i][len(partition[i])] = -1
    154     part[len(partition)] = NULL
    155     output = get_aut_gp_and_can_lab([], part, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_structures_trivial, canonical_label, base, True)
    156     mpz_set(I.value, output.order)
     148    cdef PartitionStack *part
     149    part = PS_from_list(partition)
     150    if part is NULL:
     151        raise MemoryError
     152    cdef object empty_list = []
     153    output = get_aut_gp_and_can_lab(<void *> empty_list, part, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_structures_trivial, canonical_label, NULL)
     154    SC_order(output.group, 0, I.value)
    157155    print I
    158     for i from 0 <= i < len(partition):
    159         sage_free(part[i])
    160     sage_free(part)
    161     mpz_clear(output.order)
     156    PS_dealloc(part)
     157    SC_dealloc(output.group)
    162158    sage_free(output.generators)
    163     if base:
    164         sage_free(output.base)
    165159    if canonical_label:
    166160        sage_free(output.relabeling)
    167161    sage_free(output)
    168162
    169 cdef aut_gp_and_can_lab_return *get_aut_gp_and_can_lab(object S, int **partition, int n,
    170     bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
     163def test_intersect_parabolic_with_alternating(int n=9, list partition=[[0,1,2],[3,4],[5,6,7,8]]):
     164    """
     165    A test for nontrivial input group in computing automorphism groups.
     166   
     167    TESTS::
     168
     169        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label import test_intersect_parabolic_with_alternating as tipwa
     170        sage: tipwa()
     171        144
     172        sage: tipwa(5, [[0],[1],[2],[3,4]])
     173        1
     174        sage: tipwa(5, [[0],[1],[2,3,4]])
     175        3
     176        sage: tipwa(5, [[0,1],[2,3,4]])
     177        6
     178        sage: tipwa(7, [[0,1,2,3,4,5,6]])
     179        2520
     180        sage: factorial(7)/2
     181        2520
     182        sage: tipwa(9, [[0,1,2,3],[4,5,6,7,8]])
     183        1440
     184
     185    """
     186    cdef aut_gp_and_can_lab *output
     187    cdef Integer I = Integer(0)
     188    cdef PartitionStack *part
     189    part = PS_from_list(partition)
     190    if part is NULL:
     191        raise MemoryError
     192    cdef StabilizerChain *group = SC_alternating_group(part.degree)
     193    if group is NULL:
     194        PS_dealloc(part)
     195        raise MemoryError
     196    cdef object empty_list = []
     197    output = get_aut_gp_and_can_lab(<void *> empty_list, part, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_structures_trivial, 0, group)
     198    SC_order(output.group, 0, I.value)
     199    print I
     200    PS_dealloc(part)
     201    SC_dealloc(output.group)
     202    SC_dealloc(group)
     203    sage_free(output.generators)
     204    sage_free(output)
     205
     206cdef int compare_perms(int *gamma_1, int *gamma_2, void *S1, void *S2):
     207    cdef list MS1 = <list> S1
     208    cdef list MS2 = <list> S2
     209    cdef int i, j
     210    for i from 0 <= i < len(MS1):
     211        j = cmp(MS1[gamma_1[i]], MS2[gamma_2[i]])
     212        if j != 0: return j
     213    return 0
     214
     215def coset_rep(list perm=[0,1,2,3,4,5], list gens=[[1,2,3,4,5,0]]):
     216    """
     217    Given a group G generated by the given generators, defines a map from the
     218    Symmetric group to G so that any two elements from the same right coset go
     219    to the same element. Tests nontrivial input group when computing canonical
     220    labels.
     221   
     222    TESTS::
     223
     224        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label import coset_rep
     225        sage: coset_rep()
     226        [5, 0, 1, 2, 3, 4]
     227        sage: gens = [[1,2,3,0]]
     228        sage: reps = []
     229        sage: for p in SymmetricGroup(4):
     230        ...     p = [a-1 for a in p.list()]
     231        ...     r = coset_rep(p, gens)
     232        ...     if r not in reps:
     233        ...         reps.append(r)
     234        sage: len(reps)
     235        6
     236        sage: gens = [[1,0,2,3],[0,1,3,2]]
     237        sage: reps = []
     238        sage: for p in SymmetricGroup(4):
     239        ...     p = [a-1 for a in p.list()]
     240        ...     r = coset_rep(p, gens)
     241        ...     if r not in reps:
     242        ...         reps.append(r)
     243        sage: len(reps)
     244        6
     245        sage: gens = [[1,2,0,3]]
     246        sage: reps = []
     247        sage: for p in SymmetricGroup(4):
     248        ...     p = [a-1 for a in p.list()]
     249        ...     r = coset_rep(p, gens)
     250        ...     if r not in reps:
     251        ...         reps.append(r)
     252        sage: len(reps)
     253        8
     254
     255    """
     256    cdef int i, n = len(perm)
     257    assert all(len(g) == n for g in gens)
     258    cdef aut_gp_and_can_lab *output
     259    cdef Integer I = Integer(0)
     260    cdef PartitionStack *part
     261    part = PS_new(n, 1)
     262    cdef int *c_perm = <int *> sage_malloc(n * sizeof(int))
     263    cdef StabilizerChain *group = SC_new(n, 1)
     264    if part is NULL or c_perm is NULL or group is NULL:
     265        sage_free(c_perm)
     266        PS_dealloc(part)
     267        SC_dealloc(group)
     268        raise MemoryError
     269    for g in gens:
     270        for i from 0 <= i < n:
     271            c_perm[i] = g[i]
     272        SC_insert(group, 0, c_perm, 1)
     273    output = get_aut_gp_and_can_lab(<void *> perm, part, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_perms, 1, group)
     274    SC_order(output.group, 0, I.value)
     275    assert I == 1
     276    r_inv = range(n)
     277    for i from 0 <= i < n:
     278        r_inv[output.relabeling[i]] = i
     279    label = [perm[r_inv[i]] for i in range(n)]
     280    PS_dealloc(part)
     281    SC_dealloc(output.group)
     282    SC_dealloc(group)
     283    sage_free(output.generators)
     284    sage_free(output.relabeling)
     285    sage_free(output)
     286    sage_free(c_perm)
     287    return label
     288
     289cdef aut_gp_and_can_lab *get_aut_gp_and_can_lab(void *S,
     290    PartitionStack *partition, int n,
     291    bint (*all_children_are_equivalent)(PartitionStack *PS, void *S),
    171292    int (*refine_and_return_invariant)\
    172          (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
    173     int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2),
    174     bint canonical_label, bint base, bint order):
     293         (PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len),
     294    int (*compare_structures)(int *gamma_1, int *gamma_2, void *S1, void *S2),
     295    bint canonical_label, StabilizerChain *input_group) except NULL:
    175296    """
    176297    Traverse the search space for subgroup/canonical label calculation.
    177298   
    178299    INPUT:
    179300    S -- pointer to the structure
    180     partition -- array representing a partition of the points
     301    partition -- PartitionStack representing a partition of the points
    181302    len_partition -- length of the partition
    182303    n -- the number of points (points are assumed to be 0,1,...,n-1)
    183304    canonical_label -- whether to search for canonical label - if True, return
    184305        the permutation taking S to its canonical label
    185     base -- if True, return the base for which the given generating set is a
    186         strong generating set
    187     order -- if True, return the order of the automorphism group
    188306    all_children_are_equivalent -- pointer to a function
    189307        INPUT:
    190308        PS -- pointer to a partition stack
     
    207325        OUTPUT:
    208326        int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp),
    209327            such that the set of all structures is well-ordered
     328
     329    NOTE:
     330    The partition ``partition1`` *must* satisfy the property that in each cell,
     331    the smallest element occurs first!
     332
    210333    OUTPUT:
    211     pointer to a aut_gp_and_can_lab_return struct
     334    pointer to a aut_gp_and_can_lab struct
    212335
    213336    """
    214337    cdef PartitionStack *current_ps, *first_ps, *label_ps
     
    222345    cdef int label_and_current_indicator_same = -1
    223346    cdef int compared_current_and_label_indicators
    224347   
    225     cdef OrbitPartition *orbits_of_subgroup
    226     cdef mpz_t subgroup_size
     348    cdef OrbitPartition *orbits_of_subgroup, *orbits_of_permutation
    227349    cdef int subgroup_primary_orbit_size = 0
    228350    cdef int minimal_in_primary_orbit
    229351    cdef int size_of_generator_array = 2*n*n
     
    235357   
    236358    cdef bitset_t *vertices_to_split
    237359    cdef bitset_t vertices_have_been_reduced
    238     cdef int *permutation, *id_perm, *cells_to_refine_by
     360    cdef int *permutation, *label_perm, *id_perm, *cells_to_refine_by
    239361    cdef int *vertices_determining_current_stack
     362    cdef int *perm_stack
     363    cdef StabilizerChain *group = NULL, *old_group, *tmp_gp
    240364   
    241365    cdef int i, j, k
    242366    cdef bint discrete, automorphism, update_label
    243367    cdef bint backtrack, new_vertex, narrow, mem_err = 0
    244368   
    245     cdef aut_gp_and_can_lab_return *output = <aut_gp_and_can_lab_return *> sage_malloc(sizeof(aut_gp_and_can_lab_return))
     369    cdef aut_gp_and_can_lab *output = <aut_gp_and_can_lab *> sage_malloc(sizeof(aut_gp_and_can_lab))
    246370    if output is NULL:
    247371        raise MemoryError
    248    
     372    output.group = SC_new(n)
     373    if output.group is NULL:
     374        sage_free(output)
     375        raise MemoryError
    249376    if n == 0:
    250377        output.generators = NULL
    251378        output.num_gens = 0
    252379        output.relabeling = NULL
    253         output.base = NULL
    254         output.base_size = 0
    255         mpz_init_set_si(output.order, 1)
    256380        return output
    257381   
    258382    # Allocate:
    259     mpz_init_set_si(subgroup_size, 1)
    260    
    261383    output.generators = <int *> sage_malloc( size_of_generator_array * sizeof(int) )
    262384    output.num_gens = 0
    263     if base:
    264         output.base = <int *> sage_malloc( n * sizeof(int) )
    265         output.base_size = 0
     385    cdef int *int_array = <int *> sage_malloc( 7*n * sizeof(int) )
     386    if input_group is not NULL:
     387        perm_stack = <int *> sage_malloc( n*n * sizeof(int) )
     388        group = SC_copy(input_group, n)
     389        old_group = SC_new(n)
     390        if perm_stack is NULL or group is NULL or old_group is NULL:
     391            mem_err = 1
     392        else:
     393            SC_identify(perm_stack, n)
     394    if canonical_label:
     395        label_indicators  = <int *> sage_malloc( n * sizeof(int) )
     396        output.relabeling = <int *> sage_malloc( n * sizeof(int) )
     397    else:
     398        output.relabeling = NULL
     399    cdef bitset_t *bitset_array = <bitset_t *> sage_malloc( (n+2*len_of_fp_and_mcr) * sizeof(bitset_t) )
     400    orbits_of_subgroup    = OP_new(n)
     401    orbits_of_permutation = OP_new(n)
    266402   
    267     current_indicators = <int *> sage_malloc( n * sizeof(int) )
    268     first_indicators = <int *> sage_malloc( n * sizeof(int) )
     403    if output.generators is NULL or int_array          is NULL or \
     404       bitset_array      is NULL or orbits_of_subgroup is NULL or \
     405       orbits_of_permutation is NULL:
     406        mem_err = 1
     407    elif canonical_label and (label_indicators is NULL or output.relabeling is NULL):
     408        mem_err = 1
    269409   
    270     if canonical_label:
    271         label_indicators = <int *> sage_malloc( n * sizeof(int) )
    272         output.relabeling = <int *> sage_malloc( n * sizeof(int) )
    273 
    274     fixed_points_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
    275     minimal_cell_reps_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
    276    
    277     vertices_to_split = <bitset_t *> sage_malloc( n * sizeof(bitset_t) )
    278     permutation = <int *> sage_malloc( n * sizeof(int) )
    279     id_perm = <int *> sage_malloc( n * sizeof(int) )
    280     cells_to_refine_by = <int *> sage_malloc( n * sizeof(int) )
    281     vertices_determining_current_stack = <int *> sage_malloc( n * sizeof(int) )
    282    
    283     current_ps = PS_new(n, 0)
    284     orbits_of_subgroup = OP_new(n)
    285    
    286     # Check for allocation failures:
    287     cdef bint succeeded = 1
    288     if canonical_label:
    289         if label_indicators is NULL or output.relabeling is NULL:
    290             succeeded = 0
    291             if label_indicators is not NULL:
    292                 sage_free(label_indicators)
    293             if output.relabeling is not NULL:
    294                 sage_free(output.relabeling)
    295     if base:
    296         if output.base is NULL:
    297             succeeded = 0
    298     if not succeeded                              or \
    299        output.generators                  is NULL or \
    300        current_indicators                 is NULL or \
    301        first_indicators                   is NULL or \
    302        fixed_points_of_generators         is NULL or \
    303        minimal_cell_reps_of_generators    is NULL or \
    304        vertices_to_split                  is NULL or \
    305        permutation                        is NULL or \
    306        id_perm                            is NULL or \
    307        cells_to_refine_by                 is NULL or \
    308        vertices_determining_current_stack is NULL or \
    309        current_ps                         is NULL or \
    310        orbits_of_subgroup                 is NULL:
    311         if output.generators                  is not NULL:
    312             sage_free(output.generators)
    313         if base:
    314             if output.base                        is not NULL:
    315                 sage_free(output.base)
    316         if current_indicators                 is not NULL:
    317             sage_free(current_indicators)
    318         if first_indicators                   is not NULL:
    319             sage_free(first_indicators)
    320         if fixed_points_of_generators         is not NULL:
    321             sage_free(fixed_points_of_generators)
    322         if minimal_cell_reps_of_generators    is not NULL:
    323             sage_free(minimal_cell_reps_of_generators)
    324         if vertices_to_split                  is not NULL:
    325             sage_free(vertices_to_split)
    326         if permutation                        is not NULL:
    327             sage_free(permutation)
    328         if id_perm                            is not NULL:
    329             sage_free(id_perm)
    330         if cells_to_refine_by                 is not NULL:
    331             sage_free(cells_to_refine_by)
    332         if vertices_determining_current_stack is not NULL:
    333             sage_free(vertices_determining_current_stack)
    334         if current_ps is not NULL:
    335             PS_dealloc(current_ps)
    336         if orbits_of_subgroup is not NULL:
    337             OP_dealloc(orbits_of_subgroup)
    338         sage_free(output)
    339         mpz_clear(subgroup_size)
    340         raise MemoryError
    341    
    342     # Initialize bitsets, checking for allocation failures:
    343     succeeded = 1
    344     for i from 0 <= i < len_of_fp_and_mcr:
     410    if not mem_err:
     411        current_indicators                 = int_array
     412        first_indicators                   = int_array +   n
     413        permutation                        = int_array + 2*n
     414        id_perm                            = int_array + 3*n
     415        cells_to_refine_by                 = int_array + 4*n
     416        vertices_determining_current_stack = int_array + 5*n
     417        label_perm                         = int_array + 6*n
     418       
     419        fixed_points_of_generators      = bitset_array
     420        minimal_cell_reps_of_generators = bitset_array + len_of_fp_and_mcr
     421        vertices_to_split               = bitset_array + 2*len_of_fp_and_mcr
    345422        try:
    346             bitset_init(fixed_points_of_generators[i], n)
    347         except MemoryError:
    348             succeeded = 0
    349             for j from 0 <= j < i:
    350                 bitset_free(fixed_points_of_generators[j])
    351                 bitset_free(minimal_cell_reps_of_generators[j])
    352             break
    353         try:
    354             bitset_init(minimal_cell_reps_of_generators[i], n)
    355         except MemoryError:
    356             succeeded = 0
    357             for j from 0 <= j < i:
    358                 bitset_free(fixed_points_of_generators[j])
    359                 bitset_free(minimal_cell_reps_of_generators[j])
    360             bitset_free(fixed_points_of_generators[i])
    361             break
    362     if succeeded:
    363         for i from 0 <= i < n:
    364             try:
    365                 bitset_init(vertices_to_split[i], n)
    366             except MemoryError:
    367                 succeeded = 0
    368                 for j from 0 <= j < i:
    369                     bitset_free(vertices_to_split[j])
    370                 for j from 0 <= j < len_of_fp_and_mcr:
    371                     bitset_free(fixed_points_of_generators[j])
    372                     bitset_free(minimal_cell_reps_of_generators[j])
    373                 break
    374     if succeeded:
    375         try:
     423            for i from 0 <= i < n+2*len_of_fp_and_mcr:
     424                bitset_init(bitset_array[i], n)
    376425            bitset_init(vertices_have_been_reduced, n)
    377426        except MemoryError:
    378             succeeded = 0
    379             for j from 0 <= j < n:
    380                 bitset_free(vertices_to_split[j])
    381             for j from 0 <= j < len_of_fp_and_mcr:
    382                 bitset_free(fixed_points_of_generators[j])
    383                 bitset_free(minimal_cell_reps_of_generators[j])
    384     if not succeeded:
    385         sage_free(output.generators)
    386         if base:
    387             sage_free(output.base)
    388         sage_free(current_indicators)
    389         sage_free(first_indicators)
     427            mem_err = 1
     428   
     429    if not mem_err:
     430        bitset_zero(vertices_have_been_reduced)
     431        current_ps = partition
     432        current_ps.depth = 0
     433       
     434        # default values of "infinity"
     435        for i from 0 <= i < n:
     436            first_indicators[i] = -1
    390437        if canonical_label:
    391             sage_free(output.relabeling)
    392             sage_free(label_indicators)
    393         sage_free(fixed_points_of_generators)
    394         sage_free(minimal_cell_reps_of_generators)
    395         sage_free(vertices_to_split)
    396         sage_free(permutation)
    397         sage_free(id_perm)
    398         sage_free(cells_to_refine_by)
    399         sage_free(vertices_determining_current_stack)
    400         PS_dealloc(current_ps)
    401         sage_free(output)
    402         mpz_clear(subgroup_size)
    403         raise MemoryError
    404    
    405     bitset_zero(vertices_have_been_reduced)
    406    
    407     # Copy data of partition to current_ps.
    408     i = 0
    409     j = 0
    410     while partition[i] is not NULL:
    411         k = 0
    412         while partition[i][k] != -1:
    413             current_ps.entries[j+k] = partition[i][k]
    414             current_ps.levels[j+k] = n
    415             k += 1
    416         current_ps.levels[j+k-1] = 0
    417         PS_move_min_to_front(current_ps, j, j+k-1)
    418         j += k
    419         i += 1
    420     current_ps.levels[j-1] = -1
    421     current_ps.depth = 0
    422     current_ps.degree = n
    423    
    424     # default values of "infinity"
    425     for i from 0 <= i < n:
    426         first_indicators[i] = -1
    427     if canonical_label:
     438            for i from 0 <= i < n:
     439                label_indicators[i] = -1
     440       
     441        # set up the identity permutation
    428442        for i from 0 <= i < n:
    429             label_indicators[i] = -1
    430    
    431     # set up the identity permutation
    432     for i from 0 <= i < n:
    433         id_perm[i] = i
    434    
    435     # Our first refinement needs to check every cell of the partition,
    436     # so cells_to_refine_by needs to be a list of pointers to each cell.
    437     j = 1
    438     cells_to_refine_by[0] = 0
    439     for i from 0 < i < n:
    440         if current_ps.levels[i-1] == 0:
    441             cells_to_refine_by[j] = i
    442             j += 1
    443     # Ignore the invariant this time, since we are
    444     # creating the root of the search tree.
    445     refine_and_return_invariant(current_ps, S, cells_to_refine_by, j)
    446     PS_move_all_mins_to_front(current_ps)
    447    
    448     # Refine down to a discrete partition
    449     while not PS_is_discrete(current_ps):
    450         if not all_children_are_equivalent(current_ps, S):
    451             current_kids_are_same = current_ps.depth + 1
    452         vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
    453         bitset_unset(vertices_have_been_reduced, current_ps.depth)
    454         i = current_ps.depth
    455         current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
     443            id_perm[i] = i
     444       
     445        # Our first refinement needs to check every cell of the partition,
     446        # so cells_to_refine_by needs to be a list of pointers to each cell.
     447        j = 1
     448        cells_to_refine_by[0] = 0
     449        for i from 0 < i < n:
     450            if current_ps.levels[i-1] == 0:
     451                cells_to_refine_by[j] = i
     452                j += 1
     453        # Ignore the invariant this time, since we are
     454        # creating the root of the search tree.
     455        if input_group is NULL:
     456            refine_and_return_invariant(current_ps, S, cells_to_refine_by, j)
     457        else:
     458            refine_also_by_orbits(current_ps, S, refine_and_return_invariant,
     459                cells_to_refine_by, j, group, perm_stack)
    456460        PS_move_all_mins_to_front(current_ps)
    457         first_indicators[i] = current_indicators[i]
     461
     462        # Refine down to a discrete partition
     463        while not PS_is_discrete(current_ps):
     464            i = current_ps.depth
     465            if not all_children_are_equivalent(current_ps, S):
     466                current_kids_are_same = i + 1
     467            vertices_determining_current_stack[i] = PS_first_smallest(current_ps, vertices_to_split[i])
     468            bitset_unset(vertices_have_been_reduced, i)
     469            if input_group is not NULL:
     470                # split the point
     471                current_ps.depth += 1
     472                PS_clear(current_ps)
     473                cells_to_refine_by[0] = PS_split_point(current_ps, vertices_determining_current_stack[i])
     474                # update the base
     475                tmp_gp = group
     476                group = old_group
     477                old_group = tmp_gp
     478                if SC_insert_base_point_nomalloc(group, old_group, i, vertices_determining_current_stack[i]):
     479                    mem_err = 1
     480                    break
     481                # update perm_stack
     482                SC_identify(perm_stack + n*current_ps.depth, n)
     483                # do the refinements
     484                current_indicators[i] = refine_also_by_orbits(current_ps, S, refine_and_return_invariant, cells_to_refine_by, 1, group, perm_stack)
     485            else:
     486                current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
     487            PS_move_all_mins_to_front(current_ps)
     488            first_indicators[i] = current_indicators[i]
     489            if canonical_label:
     490                label_indicators[i] = current_indicators[i]
     491            SC_add_base_point(output.group, vertices_determining_current_stack[i])
     492       
     493    if not mem_err:
     494        first_meets_current = current_ps.depth
     495        first_kids_are_same = current_ps.depth
     496        first_and_current_indicator_same = current_ps.depth
     497        first_ps = PS_copy(current_ps)
     498        if first_ps is NULL:
     499            mem_err = 1
    458500        if canonical_label:
    459             label_indicators[i] = current_indicators[i]
    460         if base:
    461             output.base[i] = vertices_determining_current_stack[i]
    462             output.base_size += 1
     501            label_meets_current = current_ps.depth
     502            label_and_current_indicator_same = current_ps.depth
     503            compared_current_and_label_indicators = 0
     504            label_ps = PS_copy(current_ps)
     505            if label_ps is NULL:
     506                mem_err = 1
     507            if input_group is not NULL:
     508                if compute_relabeling(group, old_group, label_ps.entries, label_perm):
     509                    mem_err = 1
     510        current_ps.depth -= 1
    463511
    464     first_meets_current = current_ps.depth
    465     first_kids_are_same = current_ps.depth
    466     first_and_current_indicator_same = current_ps.depth
    467     first_ps = PS_copy(current_ps)
    468     if canonical_label:
    469         label_meets_current = current_ps.depth
    470         label_and_current_indicator_same = current_ps.depth
    471         compared_current_and_label_indicators = 0
    472         label_ps = PS_copy(current_ps)
    473     current_ps.depth -= 1
    474    
    475512    # Main loop:
    476513    while current_ps.depth != -1:
    477514       
     515        if mem_err:
     516            sage_free(int_array)
     517            OP_dealloc(orbits_of_subgroup)
     518            OP_dealloc(orbits_of_permutation)
     519            if canonical_label:
     520                sage_free(label_indicators)
     521                sage_free(output.relabeling)
     522                PS_dealloc(label_ps)
     523            sage_free(output.generators)
     524            SC_dealloc(output.group)
     525            if input_group is not NULL:
     526                sage_free(perm_stack)
     527                SC_dealloc(old_group)
     528                SC_dealloc(group)
     529            sage_free(output)
     530            if bitset_array is not NULL:
     531                for i from 0 <= i < n+2*len_of_fp_and_mcr:
     532                    bitset_free(bitset_array[i])
     533            bitset_free(vertices_have_been_reduced)
     534            sage_free(bitset_array)
     535            PS_dealloc(first_ps)
     536            raise MemoryError
     537
    478538        # If necessary, update label_meets_current
    479539        if canonical_label and label_meets_current > current_ps.depth:
    480540            label_meets_current = current_ps.depth
     
    517577                # are looking at a new primary orbit (corresponding to a
    518578                # larger subgroup in the stabilizer chain).
    519579                first_meets_current = current_ps.depth
    520                 for i from 0 <= i < n:
    521                     if bitset_check(vertices_to_split[current_ps.depth], i):
    522                         minimal_in_primary_orbit = i
    523                         break
     580                minimal_in_primary_orbit = output.group.base_orbits[current_ps.depth][0]
    524581            while 1:
    525582                i = vertices_determining_current_stack[current_ps.depth]
    526583                # This was the last point to be split here.
    527                 # If it is in the same orbit as minimal_in_primary_orbit,
    528                 # then count it as an element of the primary orbit.
    529                 if OP_find(orbits_of_subgroup, i) == OP_find(orbits_of_subgroup, minimal_in_primary_orbit):
     584                # If it has been added to the primary orbit, update size info.
     585                if output.group.parents[current_ps.depth][i] != -1:
    530586                    subgroup_primary_orbit_size += 1
    531587                # Look for a new point to split.
    532588                i += 1
     
    550606                            j += 1
    551607                    if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1:
    552608                        first_kids_are_same = current_ps.depth
    553                     # Update group size and backtrack.
    554                     mpz_mul_si(subgroup_size, subgroup_size, subgroup_primary_orbit_size)
    555609                    subgroup_primary_orbit_size = 0
    556610                    current_ps.depth -= 1
    557611                    break
     
    569623        # II. Refine down to a discrete partition, or until
    570624        # we leave the part of the tree we are interested in
    571625        discrete = 0
     626        backtrack = 0
    572627        while 1:
    573628            i = current_ps.depth
    574             current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
     629            if input_group is not NULL:
     630                current_indicators[i] = split_point_and_refine_by_orbits(current_ps,
     631                    vertices_determining_current_stack[i], S,
     632                    refine_and_return_invariant, cells_to_refine_by,
     633                    group, perm_stack)
     634            else:
     635                current_indicators[i] = split_point_and_refine(current_ps,
     636                    vertices_determining_current_stack[i], S,
     637                    refine_and_return_invariant, cells_to_refine_by)
    575638            PS_move_all_mins_to_front(current_ps)
    576639            if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]:
    577640                first_and_current_indicator_same += 1
     
    586649                if compared_current_and_label_indicators > 0:
    587650                    label_indicators[i] = current_indicators[i]
    588651            if first_and_current_indicator_same < current_ps.depth and (not canonical_label or compared_current_and_label_indicators < 0):
     652                backtrack = 1
    589653                break
    590654            if PS_is_discrete(current_ps):
    591655                discrete = 1
     
    597661
    598662        # III. Check for automorphisms and labels
    599663        automorphism = 0
    600         backtrack = 1 - discrete
    601664        if discrete:
    602665            if current_ps.depth == first_and_current_indicator_same:
    603666                PS_get_perm_from(current_ps, first_ps, permutation)
    604667                if compare_structures(permutation, id_perm, S, S) == 0:
    605                     automorphism = 1
     668                    if input_group is NULL or SC_contains(group, 0, permutation, 0):
     669                        # TODO: might be slight optimization for containment using perm_stack
     670                        automorphism = 1
    606671        if not automorphism:
    607672            if (not canonical_label) or (compared_current_and_label_indicators < 0):
    608673                backtrack = 1
     
    612677                if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth):
    613678                    update_label = 1
    614679                else:
    615                     i = compare_structures(current_ps.entries, label_ps.entries, S, S)
     680                    if input_group is NULL:
     681                        i = compare_structures(current_ps.entries, label_ps.entries, S, S)
     682                    else:
     683                        if compute_relabeling(group, old_group, current_ps.entries, permutation):
     684                            mem_err = 1
     685                            break
     686                        i = compare_structures(permutation, label_perm, S, S)
    616687                    if i > 0:
    617688                        update_label = 1
    618689                    elif i < 0:
    619690                        backtrack = 1
    620691                    else:
    621                         PS_get_perm_from(current_ps, label_ps, permutation)
     692                        if input_group is NULL:
     693                            PS_get_perm_from(current_ps, label_ps, permutation)
     694                        else:
     695                            SC_invert_perm(group.perm_scratch, permutation, n)
     696                            SC_mult_perms(permutation, group.perm_scratch, label_perm, n)
    622697                        automorphism = 1
    623698                if update_label:
    624699                    PS_copy_from_to(current_ps, label_ps)
     700                    if input_group is not NULL:
     701                        memcpy(label_perm, permutation, n*sizeof(int))
    625702                    compared_current_and_label_indicators = 0
    626703                    label_meets_current = current_ps.depth
    627704                    label_and_current_indicator_same = current_ps.depth
     
    632709                index_in_fp_and_mcr += 1
    633710            bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr])
    634711            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     712            OP_clear(orbits_of_permutation)
     713            for i from 0 <= i < n:
     714                OP_join(orbits_of_permutation, i, permutation[i])
     715            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
    635716            for i from 0 <= i < n:
    636717                if permutation[i] == i:
    637718                    bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
    638719                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
    639720                else:
    640721                    bitset_unset(fixed_points_of_generators[index_in_fp_and_mcr], i)
    641                     k = i
    642                     j = permutation[i]
    643                     while j != i:
    644                         if j < k: k = j
    645                         j = permutation[j]
    646                     if k == i:
    647                         bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
    648                     else:
    649                         bitset_unset(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     722                    if orbits_of_permutation.parent[i] == i:
     723                        bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], orbits_of_permutation.mcr[i])
    650724            if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser
     725                # add permutation as a generator of the automorphism group
    651726                if n*output.num_gens == size_of_generator_array:
    652727                    # must double its size
    653728                    size_of_generator_array *= 2
    654                     output.generators = <int *> sage_realloc( output.generators, size_of_generator_array )
     729                    output.generators = <int *> sage_realloc( output.generators, size_of_generator_array * sizeof(int) )
    655730                    if output.generators is NULL:
    656731                        mem_err = True
    657                         break # main loop
     732                        continue # main loop
    658733                j = n*output.num_gens
    659734                for i from 0 <= i < n:
    660735                    output.generators[j + i] = permutation[i]
    661736                output.num_gens += 1
     737                if SC_update_tree(output.group, first_meets_current, permutation, 1):
     738                    mem_err = True
     739                    continue # main loop
    662740                if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit:
    663741                    current_ps.depth = first_meets_current
    664742                    continue # main loop
     
    691769
    692770    # End of main loop.
    693771
    694     mpz_init_set(output.order, subgroup_size)
    695772    if canonical_label:
    696773        for i from 0 <= i < n:
    697774            output.relabeling[label_ps.entries[i]] = i
    698775   
    699776    # Deallocate:
    700     for i from 0 <= i < len_of_fp_and_mcr:
    701         bitset_free(fixed_points_of_generators[i])
    702         bitset_free(minimal_cell_reps_of_generators[i])
    703     for i from 0 <= i < n:
    704         bitset_free(vertices_to_split[i])
     777    sage_free(int_array)
     778    OP_dealloc(orbits_of_subgroup)
     779    OP_dealloc(orbits_of_permutation)
     780    if input_group is not NULL:
     781        sage_free(perm_stack)
     782        SC_dealloc(old_group)
     783        SC_dealloc(group)
     784    if canonical_label:
     785        sage_free(label_indicators)
     786        PS_dealloc(label_ps)
     787    for i from 0 <= i < n+2*len_of_fp_and_mcr:
     788        bitset_free(bitset_array[i])
    705789    bitset_free(vertices_have_been_reduced)
    706     sage_free(current_indicators)
    707     sage_free(first_indicators)
    708     if canonical_label:
    709         PS_dealloc(label_ps)
    710         sage_free(label_indicators)
    711     sage_free(fixed_points_of_generators)
    712     sage_free(minimal_cell_reps_of_generators)
    713     sage_free(vertices_to_split)
    714     sage_free(permutation)
    715     sage_free(id_perm)
    716     sage_free(cells_to_refine_by)
    717     sage_free(vertices_determining_current_stack)
    718     PS_dealloc(current_ps)
     790    sage_free(bitset_array)
    719791    PS_dealloc(first_ps)
    720     OP_dealloc(orbits_of_subgroup)
    721     mpz_clear(subgroup_size)
    722     if not mem_err:
    723         return output
    724     else:
    725         mpz_clear(output.order)
    726         sage_free(output.generators)
    727         if base:
    728             sage_free(output.base)
    729         if canonical_label:
    730             sage_free(output.relabeling)
    731         sage_free(output)
    732         output = NULL
    733         raise MemoryError
     792   
     793    return output
    734794
    735795
    736796
  • sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#
    55# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    66#                         http://www.gnu.org/licenses/
     
    1010include '../../../ext/stdsage.pxi'
    1111include '../../../misc/bitset_pxd.pxi'
    1212
     13cdef extern from "stdlib.h":
     14    int rand()
     15
     16cdef extern from "FLINT/long_extras.h":
     17    int z_isprime(unsigned long n)
     18
    1319cdef struct OrbitPartition:
    1420    # Disjoint set data structure for representing the orbits of the generators
    1521    # so far found. Also keeps track of the minimum elements of the cells and
     
    3137    int depth
    3238    int degree
    3339
     40cdef struct StabilizerChain:
     41    # A representation of a permutation group acting on 0, 1, ..., degree-1.
     42    int degree
     43    int base_size
    3444
     45    int *orbit_sizes
     46    int *num_gens      # dimension of generator cube on each level
     47    int *array_size    # size of space to hold generators on each level (number of permutations)
     48
     49    int **base_orbits #
     50    int **parents     # three n*n squares, orbits and tree structures
     51    int **labels      #
     52
     53    int **generators   # generators for each level,
     54    int **gen_inverses # and their inverses
     55   
     56    bitset_s gen_used
     57    bitset_s gen_is_id
     58    int *perm_scratch
     59    OrbitPartition *OP_scratch
  • sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi
    a b  
    11r"""
    22Data structures
    33
    4 This module implements TODO...
     4This module implements basic data structures essential to the rest of the
     5partn_ref module.
    56
    6 REFERENCE:
     7REFERENCES:
    78
    89    [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
    910        Vol. 30 (1981), pp. 45-87.
     11    [2] Fredman, M. and Saks, M. The cell probe complexity of dynamic data
     12        structures. Proceedings of the Twenty-First Annual ACM Symposium on
     13        Theory of Computing, pp. 345–354. May 1989.
     14    [3] Seress, Akos. Permutation Group Algorithms. Cambridge University Press,
     15        2003.
    1016
    1117"""
    1218
    1319#*****************************************************************************
    14 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     20#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    1521#
    1622# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    1723#                         http://www.gnu.org/licenses/
     
    1925
    2026include '../../../misc/bitset.pxi'
    2127
    22 cdef inline int min(int a, int b):
    23     if a < b: return a
    24     else: return b
     28cdef extern from "math.h":
     29    float log(float x)
     30    float ceil(float x)
    2531
    26 cdef inline int max(int a, int b):
    27     if a > b: return a
    28     else: return b
     32from sage.groups.perm_gps.permgroup import PermutationGroup
     33from sage.rings.integer cimport Integer
    2934
    3035# OrbitPartitions
    3136
     
    3742    cdef int i
    3843    cdef OrbitPartition *OP = <OrbitPartition *> \
    3944                                sage_malloc(sizeof(OrbitPartition))
    40     if OP is NULL:
    41         return OP
     45    cdef int *int_array = <int *> sage_malloc( 4*n * sizeof(int) )
     46    if OP is NULL or int_array is NULL:
     47        sage_free(OP)
     48        sage_free(int_array)
     49        return NULL
    4250    OP.degree = n
    4351    OP.num_cells = n
    44     OP.parent = <int *> sage_malloc( n * sizeof(int) )
    45     OP.rank = <int *> sage_malloc( n * sizeof(int) )
    46     OP.mcr = <int *> sage_malloc( n * sizeof(int) )
    47     OP.size = <int *> sage_malloc( n * sizeof(int) )
    48     if OP.parent is NULL or OP.rank is NULL or \
    49        OP.mcr is NULL or OP.size is NULL:
    50         if OP.parent is not NULL: sage_free(OP.parent)
    51         if OP.rank is not NULL: sage_free(OP.rank)
    52         if OP.mcr is not NULL: sage_free(OP.mcr)
    53         if OP.size is not NULL: sage_free(OP.size)
    54         return NULL
     52    OP.parent = int_array
     53    OP.rank   = int_array +   n
     54    OP.mcr    = int_array + 2*n
     55    OP.size   = int_array + 3*n
     56    OP_clear(OP)
     57    return OP
     58
     59cdef inline OP_clear(OrbitPartition *OP):
     60    cdef int i, n = OP.degree
    5561    for i from 0 <= i < n:
    5662        OP.parent[i] = i
    5763        OP.rank[i] = 0
    5864        OP.mcr[i] = i
    5965        OP.size[i] = 1
    60     return OP
    6166
    6267cdef inline int OP_dealloc(OrbitPartition *OP):
    63     sage_free(OP.parent)
    64     sage_free(OP.rank)
    65     sage_free(OP.mcr)
    66     sage_free(OP.size)
     68    if OP is not NULL:
     69        sage_free(OP.parent)
    6770    sage_free(OP)
    68     return 0
    6971
    7072cdef inline int OP_find(OrbitPartition *OP, int n):
    7173    """
     
    227229    cdef int i
    228230    cdef PartitionStack *PS = <PartitionStack *> \
    229231                                sage_malloc(sizeof(PartitionStack))
    230     if PS is NULL:
    231         return PS
    232     PS.entries = <int *> sage_malloc( n * sizeof(int) )
    233     PS.levels = <int *> sage_malloc( n * sizeof(int) )
    234     if PS.entries is NULL or PS.levels is NULL:
    235         if PS.entries is not NULL: sage_free(PS.entries)
    236         if PS.levels is not NULL: sage_free(PS.levels)
     232    cdef int *int_array = <int *> sage_malloc( 2*n * sizeof(int) )
     233    if PS is NULL or int_array is NULL:
     234        sage_free(PS)
     235        sage_free(int_array)
    237236        return NULL
    238     PS.depth = 0
     237    PS.entries = int_array
     238    PS.levels  = int_array + n
     239    PS.depth  = 0
    239240    PS.degree = n
    240241    if unit_partition:
    241242        for i from 0 <= i < n-1:
     
    250251    Allocate and return a pointer to a copy of PartitionStack PS. Returns a null
    251252    pointer in the case of an allocation failure.
    252253    """
    253     cdef int i
     254    cdef int i, n = PS.degree
     255   
    254256    cdef PartitionStack *PS2 = <PartitionStack *> \
    255257                                sage_malloc(sizeof(PartitionStack))
    256     if PS2 is NULL:
    257         return PS2
    258     PS2.entries = <int *> sage_malloc( PS.degree * sizeof(int) )
    259     PS2.levels = <int *> sage_malloc( PS.degree * sizeof(int) )
    260     if PS2.entries is NULL or PS2.levels is NULL:
    261         if PS2.entries is not NULL: sage_free(PS2.entries)
    262         if PS2.levels is not NULL: sage_free(PS2.levels)
     258    cdef int *int_array = <int *> sage_malloc( 2*n * sizeof(int) )
     259    if PS2 is NULL or int_array is NULL:
     260        sage_free(PS2)
     261        sage_free(int_array)
    263262        return NULL
     263    PS2.entries = int_array
     264    PS2.levels  = int_array + n
    264265    PS_copy_from_to(PS, PS2)
    265266    return PS2
    266267
     
    270271    """
    271272    PS2.depth = PS.depth
    272273    PS2.degree = PS.degree
    273     for i from 0 <= i < PS.degree:
    274         PS2.entries[i] = PS.entries[i]
    275         PS2.levels[i] = PS.levels[i]   
     274    memcpy(PS2.entries, PS.entries, 2*PS.degree * sizeof(int) )
    276275
    277 cdef inline PartitionStack *PS_from_list(object L):
     276cdef PartitionStack *PS_from_list(list L):
    278277    """
    279278    Allocate and return a pointer to a PartitionStack representing L. Returns a
    280279    null pointer in the case of an allocation failure.
     
    283282    for cell from 0 <= cell < num_cells:
    284283        n += len(L[cell])
    285284    cdef PartitionStack *PS = PS_new(n, 0)
    286     cell = 0
    287285    if PS is NULL:
    288         return PS
    289     while 1:
     286        return NULL
     287    for cell from 0 <= cell < num_cells:
    290288        cur_len = len(L[cell])
    291289        for i from 0 <= i < cur_len:
    292290            PS.entries[cur_start + i] = L[cell][i]
    293291            PS.levels[cur_start + i] = n
    294292        PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1)
    295293        cur_start += cur_len
    296         cell += 1
    297         if cell == num_cells:
    298             PS.levels[cur_start-1] = -1
    299             break
    300294        PS.levels[cur_start-1] = 0
     295    if n > 0:
     296        PS.levels[n-1] = -1
     297    PS.depth = 0
     298    PS.degree = n
    301299    return PS
    302300
    303 cdef inline int PS_dealloc(PartitionStack *PS):
    304     sage_free(PS.entries)
    305     sage_free(PS.levels)
     301cdef inline PS_dealloc(PartitionStack *PS):
     302    if PS is not NULL:
     303        sage_free(PS.entries)
    306304    sage_free(PS)
    307     return 0
    308305
    309 cdef inline int PS_print(PartitionStack *PS):
     306cdef PS_print(PartitionStack *PS):
    310307    """
    311308    Print a visual representation of PS.
    312309    """
     
    315312        PS_print_partition(PS, i)
    316313    print
    317314
    318 cdef inline int PS_print_partition(PartitionStack *PS, int k):
     315cdef PS_print_partition(PartitionStack *PS, int k):
    319316    """
    320317    Print the partition at depth k.
    321318    """
     
    349346            ncells += 1
    350347    return ncells
    351348
    352 cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end):
     349cdef inline PS_move_min_to_front(PartitionStack *PS, int start, int end):
    353350    """
    354351    Makes sure that the first element of the segment of entries i with
    355352    start <= i <= end is minimal.
     
    383380    """
    384381    cdef int i, cur_start = 0
    385382    for i from 0 <= i < PS.degree:
    386         if PS.levels[i] >= PS.depth:
     383        if PS.levels[i] == PS.depth:
    387384            PS.levels[i] += 1
    388385        if PS.levels[i] < PS.depth:
    389386            PS_move_min_to_front(PS, cur_start, i)
     
    399396            PS_move_min_to_front(PS, cur_start, i)
    400397            cur_start = i+1
    401398
    402 cdef inline int PS_split_point(PartitionStack *PS, int v):
     399cdef int PS_split_point(PartitionStack *PS, int v):
    403400    """
    404401    Detaches the point v from the cell it is in, putting the singleton cell
    405402    of just v in front. Returns the position where v is now located.
     
    430427        PS.levels[i] = PS.depth
    431428        return i
    432429
    433 cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b):
     430cdef int PS_first_smallest(PartitionStack *PS, bitset_t b):
    434431    """
    435432    Find the first occurrence of the smallest cell of size greater than one,
    436433    store its entries to b, and return its minimum element.
     
    454451        i += 1
    455452    return PS.entries[location]
    456453
     454cdef int PS_find_element(PartitionStack *PS, bitset_t b, int x):
     455    """
     456    Find the cell containing x, store its entries to b and return the location
     457    of the beginning of the cell.
     458    """
     459    cdef int i, location, n = PS.degree
     460    bitset_zero(b)
     461    for i from 0 <= i < n:
     462        if PS.entries[i] == x:
     463            location = i
     464            break
     465    while location > 0 and PS.levels[location-1] > PS.depth:
     466        location -= 1
     467    i = 0
     468    while 1:
     469        bitset_set(b, PS.entries[location+i])
     470        if PS.levels[location+i] <= PS.depth:
     471            break
     472        i += 1
     473    return location
     474
    457475cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma):
    458476    """
    459477    Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i]
     
    463481    for i from 0 <= i < PS1.degree:
    464482        gamma[PS2.entries[i]] = PS1.entries[i]
    465483
     484cdef inline bint stacks_are_equivalent(PartitionStack *PS1, PartitionStack *PS2):
     485    cdef int i, j, depth = min(PS1.depth, PS2.depth)
     486    for i from 0 <= i < PS1.degree:
     487        j = cmp(PS1.levels[i], PS2.levels[i])
     488        if j == 0: continue
     489        if ( (j < 0 and PS1.levels[i] <= depth and PS2.levels[i] > depth)
     490            or (j > 0 and PS2.levels[i] <= depth and PS1.levels[i] > depth) ):
     491            return 0
     492    return 1
     493
    466494def PS_represent(partition, splits):
    467495    """
    468496    Demonstration and testing.
     
    514542        Finding first smallest:
    515543        Minimal element is 5, bitset is:
    516544        0000010001
     545        Finding element 1:
     546        Location is: 5
     547        Bitset is:
     548        0100000000
    517549        Deallocating PartitionStacks.
    518550        Done.
    519551
     
    591623    print "Minimal element is %d, bitset is:"%i
    592624    print bitset_string(b)
    593625    bitset_free(b)
     626    print "Finding element 1:"
     627    bitset_init(b, n)
     628    print "Location is:", PS_find_element(PS, b, 1)
     629    print "Bitset is:"
     630    print bitset_string(b)
     631    bitset_free(b)
    594632    print "Deallocating PartitionStacks."
    595633    PS_dealloc(PS)
    596634    PS_dealloc(PS2)
    597635    print "Done."
    598636
     637# StabilizerChains
     638
     639cdef enum:
     640    default_num_gens = 8
     641    default_num_bits = 64
     642
     643cdef StabilizerChain *SC_new(int n, bint init_gens=True):
     644    """
     645    Allocate and return a pointer to a new StabilizerChain of degree n. Returns
     646    a null pointer in the case of an allocation failure.
     647    """
     648    cdef int i
     649    cdef StabilizerChain *SC = <StabilizerChain *> \
     650                                sage_malloc(sizeof(StabilizerChain))
     651    cdef int *array1, *array2, *array3
     652    cdef bint mem_err = 0
     653    if SC is NULL:
     654        return NULL
     655    SC.degree = n
     656    SC.base_size = 0
     657    if n == 0:
     658        SC.orbit_sizes    = NULL
     659        SC.num_gens       = NULL
     660        SC.array_size     = NULL
     661        SC.base_orbits    = NULL
     662        SC.parents        = NULL
     663        SC.labels         = NULL
     664        SC.generators     = NULL
     665        SC.gen_inverses   = NULL
     666        SC.gen_used.bits  = NULL
     667        SC.gen_is_id.bits = NULL
     668        SC.perm_scratch   = NULL
     669        SC.OP_scratch     = NULL
     670        return SC
     671
     672    # first level allocations
     673    cdef int *int_array = <int *>  sage_malloc( (3*n*n + 6*n + 1) * sizeof(int) )
     674    cdef int **int_ptrs = <int **> sage_malloc( 5*n * sizeof(int *) )
     675    SC.OP_scratch = OP_new(n)
     676    # bitset_init without the MemoryError:
     677    cdef long limbs = (default_num_bits - 1)/(8*sizeof(unsigned long)) + 1
     678    SC.gen_used.size   = default_num_bits
     679    SC.gen_is_id.size  = default_num_bits
     680    SC.gen_used.limbs  = limbs
     681    SC.gen_is_id.limbs = limbs
     682    SC.gen_used.bits   = <unsigned long*>sage_malloc(limbs * sizeof(unsigned long))
     683    SC.gen_is_id.bits  = <unsigned long*>sage_malloc(limbs * sizeof(unsigned long))
     684    SC.gen_used.bits[limbs-1] = 0
     685    SC.gen_is_id.bits[limbs-1] = 0
     686   
     687    # check for allocation failures
     688    if int_array        is NULL or int_ptrs          is NULL or \
     689       SC.gen_used.bits is NULL or SC.gen_is_id.bits is NULL or \
     690       SC.OP_scratch    is NULL:
     691        SC_dealloc(SC)
     692        return NULL
     693
     694    SC.orbit_sizes  = int_array
     695    SC.num_gens     = int_array +   n
     696    SC.array_size   = int_array + 2*n
     697    SC.perm_scratch = int_array + 3*n # perm_scratch is length 3*n+1 for sorting
     698    int_array += 6*n + 1
     699   
     700    SC.generators   = int_ptrs
     701    SC.gen_inverses = int_ptrs +   n
     702    SC.base_orbits  = int_ptrs + 2*n
     703    SC.parents      = int_ptrs + 3*n
     704    SC.labels       = int_ptrs + 4*n
     705    for i from 0 <= i < n:
     706        SC.base_orbits[i] = int_array
     707        SC.parents[i]     = int_array +   n
     708        SC.labels[i]      = int_array + 2*n
     709        int_array += 3*n
     710   
     711    # second level allocations
     712    if init_gens:
     713        for i from 0 <= i < n:
     714            SC.array_size[i] = default_num_gens
     715            SC.generators[i]   = <int *> sage_malloc( default_num_gens*n * sizeof(int) )
     716            SC.gen_inverses[i] = <int *> sage_malloc( default_num_gens*n * sizeof(int) )
     717            if SC.generators[i] is NULL or SC.gen_inverses[i] is NULL:
     718                SC_dealloc(SC)
     719                return NULL
     720   
     721    return SC
     722
     723cdef StabilizerChain *SC_symmetric_group(int n):
     724    """
     725    Returns a stabilizer chain for the symmetric group on {0, 1, ..., n-1}.
     726   
     727    Returns NULL in the case of an allocation failure.
     728    """
     729    cdef int i, j, b
     730    cdef StabilizerChain *SC = SC_new(n, False)
     731    if SC is NULL:
     732        return NULL
     733    SC.base_size = n-1
     734    for i from 0 <= i < n-1:
     735        SC.array_size[i] = n-i-1
     736    SC.array_size[n-1] = default_num_gens
     737    for i from 0 <= i < n:
     738        SC.generators[i]   = <int *> sage_malloc( SC.array_size[i]*n * sizeof(int) )
     739        SC.gen_inverses[i] = <int *> sage_malloc( SC.array_size[i]*n * sizeof(int) )
     740        if SC.generators[i] is NULL or SC.gen_inverses[i] is NULL:
     741            SC_dealloc(SC)
     742            return NULL
     743    cdef int *id_perm = SC.perm_scratch
     744    for i from 0 <= i < n:
     745        id_perm[i] = i
     746    for i from 0 <= i < n-1:
     747        b = i
     748        SC.orbit_sizes[i] = n-i
     749        SC.num_gens[i] = n-i-1
     750        for j from 0 <= j < i:
     751            SC.parents[i][j] = -1
     752        for j from 0 <= j < n-i:
     753            SC.base_orbits[i][j] = i+j
     754            SC.parents[i][i+j] = b
     755            SC.labels[i][i+j] = j
     756        for j from 0 <= j < n-i-1:
     757            #j-th generator sends i+j+1 to b
     758            memcpy(SC.generators[i] + n*j, id_perm, n * sizeof(int) )
     759            SC.generators[i][n*j + i+j+1] = b
     760            SC.generators[i][n*j + b] = i+j+1
     761            memcpy(SC.gen_inverses[i] + n*j, SC.generators[i] + n*j, n * sizeof(int) )
     762    return SC
     763
     764cdef StabilizerChain *SC_alternating_group(int n):
     765    """
     766    Returns a stabilizer chain for the alternating group on {0, 1, ..., n-1}.
     767   
     768    Returns NULL in the case of an allocation failure.
     769    """
     770    cdef int i, j, b
     771    cdef StabilizerChain *SC = SC_new(n, False)
     772    if SC is NULL:
     773        return NULL
     774    SC.base_size = n-2
     775    for i from 0 <= i < n-2:
     776        SC.array_size[i] = n-i-1
     777    SC.array_size[n-2] = default_num_gens
     778    SC.array_size[n-1] = default_num_gens
     779    for i from 0 <= i < n:
     780        SC.generators[i]   = <int *> sage_malloc( SC.array_size[i]*n * sizeof(int) )
     781        SC.gen_inverses[i] = <int *> sage_malloc( SC.array_size[i]*n * sizeof(int) )
     782        if SC.generators[i] is NULL or SC.gen_inverses[i] is NULL:
     783            SC_dealloc(SC)
     784            return NULL
     785    cdef int *id_perm = SC.perm_scratch
     786    for i from 0 <= i < n:
     787        id_perm[i] = i
     788    for i from 0 <= i < n-2:
     789        b = i
     790        SC.orbit_sizes[i] = n-i
     791        SC.num_gens[i] = n-i-2
     792        for j from 0 <= j < i:
     793            SC.parents[i][j] = -1
     794        for j from 0 <= j < n-i:
     795            SC.base_orbits[i][j] = i+j
     796            SC.parents[i][i+j] = b
     797            SC.labels[i][i+j] = j
     798        SC.labels[i][n-1] = -(n-i-2)
     799        for j from 0 <= j < n-i-2:
     800            #j-th generator sends i+j+1 to b, i+j+2 to i+j+1, and b to i+j+2
     801            memcpy(SC.generators[i] + n*j, id_perm, n * sizeof(int) )
     802            SC.generators[i][n*j + i+j+1] = b
     803            SC.generators[i][n*j + b    ] = i+j+2
     804            SC.generators[i][n*j + i+j+2] = i+j+1
     805            SC_invert_perm(SC.gen_inverses[i] + n*j, SC.generators[i] + n*j, n)
     806    return SC
     807
     808cdef inline int SC_realloc_gens(StabilizerChain *SC, int level, int size):
     809    """
     810    Reallocate generator array at level `level` to size `size`.
     811   
     812    Returns 1 in case of an allocation failure.
     813    """
     814    cdef int *temp, n = SC.degree
     815
     816    temp = <int *> sage_realloc( SC.generators[level],   n * size * sizeof(int) )
     817    if temp is NULL: return 1
     818    SC.generators[level] = temp
     819   
     820    temp = <int *> sage_realloc( SC.gen_inverses[level], n * size * sizeof(int) )
     821    if temp is NULL: return 1
     822    SC.gen_inverses[level] = temp
     823   
     824    SC.array_size[level] = size
     825    return 0
     826
     827cdef int SC_realloc_bitsets(StabilizerChain *SC, long size):
     828    """
     829    If size is larger than current allocation, double the size of the bitsets
     830    until it is not.
     831   
     832    Returns 1 in case of an allocation failure.
     833    """
     834    cdef unsigned long size_old = SC.gen_used.size
     835    if size <= size_old: return 0
     836    cdef unsigned long new_size = size_old
     837    while new_size < size:
     838        new_size *= 2
     839    cdef unsigned long limbs_old = SC.gen_used.limbs
     840    cdef long limbs = (new_size - 1)/(8*sizeof(unsigned long)) + 1
     841    cdef unsigned long *tmp = <unsigned long *> sage_realloc(SC.gen_used.bits, limbs * sizeof(unsigned long))
     842    if tmp is not NULL:
     843        SC.gen_used.bits = tmp
     844    else:
     845        return 1
     846    tmp = <unsigned long *> sage_realloc(SC.gen_is_id.bits, limbs * sizeof(unsigned long))
     847    if tmp is not NULL:
     848        SC.gen_is_id.bits = tmp
     849    else:
     850        return 1
     851    SC.gen_used.limbs = limbs
     852    SC.gen_is_id.limbs = limbs
     853    SC.gen_used.size = new_size
     854    SC.gen_is_id.size = new_size
     855    SC.gen_used.bits[size_old >> index_shift] &= ((<unsigned long>1 << (size_old & offset_mask)) - 1)
     856    memset(SC.gen_used.bits + (size_old >> index_shift) + 1, 0, (limbs - (size_old >> index_shift) - 1) * sizeof(unsigned long))
     857    SC.gen_is_id.bits[size_old >> index_shift] &= ((<unsigned long>1 << (size_old & offset_mask)) - 1)
     858    memset(SC.gen_is_id.bits + (size_old >> index_shift) + 1, 0, (limbs - (size_old >> index_shift) - 1) * sizeof(unsigned long))
     859    return 0
     860
     861cdef inline SC_dealloc(StabilizerChain *SC):
     862    cdef int i, n = SC.degree
     863    if SC is not NULL:
     864        if SC.generators is not NULL:
     865            for i from 0 <= i < n:
     866                sage_free(SC.generators[i])
     867                sage_free(SC.gen_inverses[i])
     868        sage_free(SC.generators) # frees int_ptrs
     869        sage_free(SC.orbit_sizes) # frees int_array
     870        sage_free(SC.gen_used.bits)
     871        sage_free(SC.gen_is_id.bits)
     872        OP_dealloc(SC.OP_scratch)
     873    sage_free(SC)
     874
     875cdef StabilizerChain *SC_copy(StabilizerChain *SC, int level):
     876    """
     877    Creates a copy of the first `level` levels of SC. Must have 0 < level <= SC.base_size.
     878
     879    Returns a null pointer in case of allocation failure.
     880    """
     881    cdef int i, n = SC.degree
     882    cdef StabilizerChain *SCC = SC_new(n, False)
     883    if SCC is NULL:
     884        return NULL
     885    level = min(level, SC.base_size)
     886    for i from 0 <= i < level:
     887        SCC.generators[i]   = <int *> sage_malloc( SC.array_size[i]*n * sizeof(int) )
     888        SCC.gen_inverses[i] = <int *> sage_malloc( SC.array_size[i]*n * sizeof(int) )
     889        if SCC.generators[i] is NULL or SCC.gen_inverses[i] is NULL:
     890            SC_dealloc(SCC)
     891            return NULL
     892        SCC.array_size[i] = SC.array_size[i]
     893    for i from level <= i < n:
     894        SCC.generators[i]   = <int *> sage_malloc( default_num_gens*n * sizeof(int) )
     895        SCC.gen_inverses[i] = <int *> sage_malloc( default_num_gens*n * sizeof(int) )
     896        if SCC.generators[i] is NULL or SCC.gen_inverses[i] is NULL:
     897            SC_dealloc(SCC)
     898            return NULL
     899        SCC.array_size[i] = default_num_gens
     900    SC_copy_nomalloc(SCC, SC, level) # no chance for memory error here...
     901    return SCC
     902
     903cdef int SC_copy_nomalloc(StabilizerChain *SC_dest, StabilizerChain *SC, int level):
     904    cdef int i, n = SC.degree
     905    level = min(level, SC.base_size)
     906    SC_dest.base_size = level
     907    memcpy(SC_dest.orbit_sizes, SC.orbit_sizes, 2*n * sizeof(int) ) # copies orbit_sizes, num_gens
     908    memcpy(SC_dest.base_orbits[0], SC.base_orbits[0], 3*n*n * sizeof(int) ) # copies base_orbits, parents, labels
     909    for i from 0 <= i < level:
     910        if SC.num_gens[i] > SC_dest.array_size[i]:
     911            if SC_realloc_gens(SC_dest, i, max(SC.num_gens[i], 2*SC_dest.array_size[i])):
     912                return 1
     913        memcpy(SC_dest.generators[i],   SC.generators[i],   SC.num_gens[i]*n * sizeof(int) )
     914        memcpy(SC_dest.gen_inverses[i], SC.gen_inverses[i], SC.num_gens[i]*n * sizeof(int) )
     915    return 0
     916
     917cdef inline int SC_perm_is_identity(int *perm, int degree):
     918    for i from 0 <= i < degree:
     919        if perm[i] != i:
     920            break
     921    else:
     922        return 1
     923    return 0
     924
     925cdef inline SC_mult_perms(int *out, int *first, int *second, int degree):
     926    """
     927    DON'T DO THIS WITH out == second!
     928    """
     929    cdef int i
     930    for i from 0 <= i < degree:
     931        out[i] = second[first[i]]
     932
     933cdef inline SC_invert_perm(int *out, int *input, int degree):
     934    """
     935    DON'T DO THIS WITH out == in!
     936    """
     937    cdef int i
     938    for i from 0 <= i < degree:
     939        out[input[i]] = i
     940
     941cdef inline SC_identify(int *perm, int degree):
     942    cdef int i
     943    for i from 0 <= i < degree:
     944        perm[i] = i
     945
     946cdef SC_print_level(StabilizerChain *SC, int level):
     947    cdef int i, j, n = SC.degree
     948    if level < SC.base_size:
     949        print '/ level ', level
     950        print '| orbit   ', [SC.base_orbits[level][i] for i from 0 <= i < SC.orbit_sizes[level]]
     951        print '| parents ', [SC.parents       [level][i] for i from 0 <= i < n]
     952        print '| labels  ', [SC.labels        [level][i] for i from 0 <= i < n]
     953        print '|'
     954        print '| generators  ', [[SC.generators  [level][n*i + j] for j from 0 <= j < n] for i from 0 <= i < SC.num_gens[level]]
     955        print '\ inverses    ', [[SC.gen_inverses[level][n*i + j] for j from 0 <= j < n] for i from 0 <= i < SC.num_gens[level]]
     956    else:
     957        print '/ level ', level
     958        print '|'
     959        print '\ base_size ', SC.base_size
     960
     961cdef inline SC_add_base_point(StabilizerChain *SC, int b):
     962    """
     963    Adds base point b to the end of SC. Assumes b is not already in the base.
     964    """
     965    cdef int i, n = SC.degree
     966    SC.orbit_sizes[SC.base_size] = 1
     967    SC.num_gens[SC.base_size] = 0
     968    SC.base_orbits[SC.base_size][0] = b
     969    for i from 0 <= i < n:
     970        SC.parents[SC.base_size][i] = -1
     971    SC.parents[SC.base_size][b] = b
     972    SC.labels[SC.base_size][b] = 0
     973    SC.base_size += 1
     974
     975cdef int SC_update(StabilizerChain *dest, StabilizerChain *source, int level):
     976    cdef mpz_t src_order, dst_order
     977    cdef int *perm = dest.perm_scratch
     978    mpz_init(src_order)
     979    mpz_init(dst_order)
     980    SC_order(source, level, src_order)
     981    SC_order(dest, level, dst_order)
     982    while mpz_cmp(dst_order, src_order):
     983        SC_random_element(source, level, perm)
     984        if SC_insert(dest, level, perm, 1):
     985            mpz_clear(dst_order)
     986            mpz_clear(src_order)
     987            return 1
     988        SC_order(dest, level, dst_order)
     989    mpz_clear(src_order)
     990    mpz_clear(dst_order)
     991    return 0
     992
     993cdef StabilizerChain *SC_insert_base_point(StabilizerChain *SC, int level, int p):
     994    """
     995    Insert the point ``p`` as a base point on level ``level``. Return a new
     996    StabilizerChain with this new base. Original StabilizerChain is unmodified.
     997   
     998    Use SC_cleanup to remove redundant base points.
     999   
     1000    Returns a null pointer in case of an allocation failure.
     1001    """
     1002    cdef int i, b, n = SC.degree
     1003    cdef StabilizerChain *NEW = SC_copy(SC, level)
     1004    if NEW is NULL:
     1005        return NULL
     1006    SC_add_base_point(NEW, p)
     1007    for i from level <= i < SC.base_size:
     1008        b = SC.base_orbits[i][0]
     1009        if b != p:
     1010            SC_add_base_point(NEW, b)
     1011    if SC_update(NEW, SC, level):
     1012        SC_dealloc(NEW)
     1013        return NULL
     1014    return NEW
     1015
     1016cdef int SC_insert_base_point_nomalloc(StabilizerChain *SC_dest, StabilizerChain *SC, int level, int p):
     1017    cdef int i, b, n = SC.degree
     1018    SC_copy_nomalloc(SC_dest, SC, level)
     1019    SC_add_base_point(SC_dest, p)
     1020    for i from level <= i < SC.base_size:
     1021        b = SC.base_orbits[i][0]
     1022        if b != p:
     1023            SC_add_base_point(SC_dest, b)
     1024    if SC_update(SC_dest, SC, level):
     1025        return 1
     1026    return 0
     1027
     1028cdef inline StabilizerChain *SC_new_base(StabilizerChain *SC, int *base, int base_len):
     1029    """
     1030    Create a new stabilizer chain whose base starts with the given base, and
     1031    which represents the same permutation group. Original StabilizerChain is
     1032    unmodified.
     1033   
     1034    Use SC_cleanup to remove redundant base points.
     1035   
     1036    Returns a null pointer in case of an allocation failure.
     1037    """
     1038    cdef StabilizerChain *NEW = SC_new(SC.degree)
     1039    if NEW is NULL:
     1040        return NULL
     1041    if SC_new_base_nomalloc(NEW, SC, base, base_len):
     1042        SC_dealloc(NEW)
     1043        return NULL
     1044    return NEW
     1045
     1046cdef inline int SC_new_base_nomalloc(StabilizerChain *SC_dest, StabilizerChain *SC, int *base, int base_len):
     1047    cdef int i, n = SC.degree
     1048    SC_dest.base_size = 0
     1049    for i from 0 <= i < base_len:
     1050        SC_add_base_point(SC_dest, base[i])
     1051    if SC_update(SC_dest, SC, 0):
     1052        SC_dealloc(SC_dest)
     1053        return 1
     1054    return 0
     1055
     1056cdef inline int SC_cleanup(StabilizerChain *SC):
     1057    """
     1058    Remove redundant base elements from SC.
     1059
     1060    Returns 1 if nothing changed, and 2 in case of an allocation failure.
     1061    """
     1062    cdef int old, new = 0, i, n = SC.degree
     1063    for old from 0 <= old < SC.base_size:
     1064        if SC.orbit_sizes[old] != 1:
     1065            if old != new:
     1066                # copy row old to row new
     1067                SC.orbit_sizes[new] = SC.orbit_sizes[old]
     1068                SC.num_gens[new] = SC.num_gens[old]
     1069                if SC.array_size[new] < SC.array_size[old]:
     1070                    if SC_realloc_gens(SC, new, max(SC.array_size[old], 2*SC.array_size[new])):
     1071                        return 2
     1072                memcpy(SC.base_orbits[new],  SC.base_orbits[old],    n * sizeof(int))
     1073                memcpy(SC.parents[new],      SC.parents[old],        n * sizeof(int))
     1074                memcpy(SC.labels[new],       SC.labels[old],         n * sizeof(int))
     1075                memcpy(SC.generators[new],   SC.generators[old],     n*SC.num_gens[old] * sizeof(int))
     1076                memcpy(SC.gen_inverses[new], SC.gen_inverses[old],   n*SC.num_gens[old] * sizeof(int))
     1077            new += 1
     1078    old = SC.base_size
     1079    SC.base_size = new
     1080    return (old == new)
     1081
     1082cdef inline SC_compose_up_to_base(StabilizerChain *SC, int level, int x, int *perm):
     1083    """
     1084    Repeatedly compose the given perm by labels on the Schreier tree, starting
     1085    with x, until the base is reached. The composition is stored to perm.
     1086    """
     1087    cdef int b = SC.base_orbits[level][0], n = SC.degree
     1088    cdef int *label, label_no
     1089    while x != b:
     1090        label_no = SC.labels[level][x]
     1091        if label_no < 0:
     1092            label_no = -label_no - 1
     1093            label = SC.gen_inverses[level] + n*label_no
     1094        else:
     1095            label_no = label_no - 1
     1096            label = SC.generators[level] + n*label_no
     1097        x = SC.parents[level][x]
     1098        SC_mult_perms(perm, perm, label, n)
     1099
     1100cdef inline SC_scan(StabilizerChain *SC, int level, int x, int gen_index, int *gen, int sign):
     1101    """
     1102    See whether the point x is moved to a point outside the
     1103    tree by gen, and if so add it to the tree (arc label is gen_inv).
     1104   
     1105    gen_index - where in the generator array the generator is located
     1106    gen - points to the generator
     1107    gen_inv - points to the inverse
     1108    sign - whether to take SC.generators or SC.gen_inverses to go *up* the tree
     1109    """
     1110    cdef int y = gen[x], n = SC.degree
     1111    if SC.parents[level][y] == -1:
     1112        SC.base_orbits[level][SC.orbit_sizes[level]] = y
     1113        SC.orbit_sizes[level] += 1
     1114        SC.parents[level][y] = x
     1115        SC.labels[level][y] = sign*(gen_index+1)
     1116
     1117cdef int SC_re_tree(StabilizerChain *SC, int level, int *perm, int x):
     1118    """
     1119    Return values:
     1120    0 - No errors.
     1121    1 - Allocation failure.
     1122    """
     1123    cdef int *gen, *gen_inv
     1124    cdef int i, b, gen_index, error, n = SC.degree
     1125
     1126    # make sure we have room for the new generator:
     1127    if SC.array_size[level] == SC.num_gens[level]:
     1128        if SC_realloc_gens(SC, level, 2*SC.array_size[level]):
     1129            return 1
     1130    cdef int *new_gen     = SC.generators  [level] + n*SC.num_gens[level]
     1131    cdef int *new_gen_inv = SC.gen_inverses[level] + n*SC.num_gens[level]
     1132
     1133    # new generator is perm^(-1) * (path from x to base) (left to right composition)
     1134    SC_invert_perm(new_gen, perm, n)
     1135    SC_compose_up_to_base(SC, level, x, new_gen)
     1136    SC_invert_perm(new_gen_inv, new_gen, n)
     1137    SC.num_gens[level] += 1
     1138
     1139    # now that we have our generators, regenerate the tree, breadth-first
     1140    b = SC.base_orbits[level][0]
     1141    for i from 0 <= i < n:
     1142        SC.parents[level][i] = -1
     1143    SC.parents[level][b] = b
     1144    i = 0
     1145    SC.orbit_sizes[level] = 1
     1146    while i < SC.orbit_sizes[level]:
     1147        x = SC.base_orbits[level][i]
     1148        for gen_index from SC.num_gens[level] > gen_index >= 0:
     1149            gen_inv = SC.gen_inverses[level] + n*gen_index
     1150            SC_scan(SC, level, x, gen_index, gen_inv, 1)
     1151        for gen_index from 0 <= gen_index < SC.num_gens[level]:
     1152            gen     = SC.generators  [level] + n*gen_index
     1153            SC_scan(SC, level, x, gen_index, gen, -1)
     1154        i += 1
     1155    return 0
     1156
     1157cdef int SC_sift(StabilizerChain *SC, int level, int x, int *gens, int num_gens, int *new_gens):
     1158    """
     1159    Apply Schreier's subgroup lemma[1] as follows. Given a level, a point x, and
     1160    a generator s, find the coset traversal element r coming from x.
     1161    Try inserting rsR(rs)^(-1) (composing left to right) on level+1, where
     1162    R(g) is the coset representative of g.
     1163
     1164    level - which level we are currently working on
     1165    x - the object representing r
     1166    gens - points to the generators to sift by
     1167    num_gens - how many of these there are
     1168    new_gens - space of size at least num_gens*n for the sifted perms to go
     1169   
     1170    Returns 1 in case of an allocation failure.
     1171    """
     1172    cdef int n = SC.degree
     1173    if num_gens == 0:
     1174        return 0
     1175
     1176    # copy a representative taking base to the point x to each of these
     1177    cdef int i
     1178    cdef int *temp = SC.gen_inverses[level] + n*SC.num_gens[level] # one more scratch space
     1179                                                                   # (available since num_gens > 0)
     1180    cdef int *rep_inv = temp
     1181    SC_identify(rep_inv, n)
     1182    SC_compose_up_to_base(SC, level, x, rep_inv)
     1183    SC_invert_perm(new_gens, rep_inv, n)
     1184    for i from 0 < i < num_gens:
     1185        memcpy(new_gens + n*i, new_gens, n*sizeof(int))
     1186
     1187    # post-compose each one with a generator
     1188    for i from 0 <= i < num_gens:
     1189        SC_mult_perms(new_gens + n*i, new_gens + n*i, gens + n*i, n)
     1190
     1191    # for each one, look up the representative it now gives, and post-compose with that inverse
     1192    cdef int y, b = SC.base_orbits[level][0]
     1193    cdef int *perm
     1194    cdef int *perm_rep_inv = temp
     1195    cdef int j
     1196    for i from 0 <= i < num_gens:
     1197        perm = new_gens + n*i # this is now rs
     1198        y = perm[b]
     1199        SC_identify(perm_rep_inv, n)
     1200        SC_compose_up_to_base(SC, level, y, perm_rep_inv)
     1201        SC_mult_perms(perm, perm, perm_rep_inv, n)
     1202    return SC_insert(SC, level+1, new_gens, num_gens)
     1203
     1204cdef int SC_insert_and_sift(StabilizerChain *SC, int level, int *pi, int num_perms, bint sift):
     1205    cdef int i, j, b, n = SC.degree
     1206    cdef int perm_gen_index
     1207    cdef int max_orbit_size, max_orbit_place
     1208    if sift:
     1209        if SC_realloc_bitsets(SC, num_perms):
     1210            return 1
     1211        bitset_clear(&SC.gen_used)
     1212        bitset_clear(&SC.gen_is_id)
     1213        b = -1
     1214        for perm_gen_index from 0 <= perm_gen_index < num_perms:
     1215            for i from 0 <= i < n:
     1216                if pi[n*perm_gen_index + i] != i:
     1217                    b = i
     1218                    break
     1219            else:
     1220                bitset_set(&SC.gen_is_id, perm_gen_index)
     1221            if b != -1: break
     1222        if b == -1: return 0
     1223    if sift and level == SC.base_size:
     1224        SC_add_base_point(SC, b)
     1225    else:
     1226        b = SC.base_orbits[level][0]
     1227    # Now b is the base element at level `level`
     1228
     1229    # Record the old orbit elements and the old generators (see sifting phase)
     1230    cdef int old_num_gens = SC.num_gens[level]
     1231    cdef int old_num_points = SC.orbit_sizes[level]
     1232
     1233    # Add new points to the tree:
     1234    cdef int x
     1235    cdef int *perm
     1236    cdef int start_over = 1
     1237    cdef int error
     1238    cdef int re_treed = 0
     1239    while start_over:
     1240        start_over = 0
     1241        for i from 0 <= i < SC.orbit_sizes[level]:
     1242            x = SC.base_orbits[level][i]
     1243            for perm_gen_index from 0 <= perm_gen_index < num_perms:
     1244                if sift and bitset_check(&SC.gen_is_id, perm_gen_index): continue
     1245                perm = pi + n*perm_gen_index
     1246                if SC.parents[level][perm[x]] == -1:
     1247                    # now we have an x which maps to a new point under perm,
     1248                    re_treed = 1
     1249                    if sift: bitset_set(&SC.gen_used, perm_gen_index)
     1250                    if SC_re_tree(SC, level, perm, x):
     1251                        return 1
     1252                    start_over = 1 # we must look anew
     1253                    break
     1254            if start_over: break
     1255            if not re_treed: continue
     1256            for perm_gen_index from 0 <= perm_gen_index < old_num_gens:
     1257                perm = SC.generators[level] + n*perm_gen_index
     1258                if SC.parents[level][perm[x]] == -1:
     1259                    # now we have an x which maps to a new point under perm,
     1260                    if SC_re_tree(SC, level, perm, x):
     1261                        return 1
     1262                    start_over = 1 # we must look anew
     1263                    break
     1264            if start_over: break
     1265            for j from level < j < SC.base_size:
     1266                for perm_gen_index from 0 <= perm_gen_index < SC.num_gens[j]:
     1267                    perm = SC.generators[j] + n*perm_gen_index
     1268                    if SC.parents[level][perm[x]] == -1:
     1269                        # now we have an x which maps to a new point under perm,
     1270                        if SC_re_tree(SC, level, perm, x):
     1271                            return 1
     1272                        start_over = 1 # we must look anew
     1273                        break
     1274    if not sift:
     1275        return 0
     1276
     1277    # store the rest of the unused generators in the generator array, after the added ones
     1278    cdef int new_size
     1279    cdef int unused_gens = 0
     1280    for perm_gen_index from 0 <= perm_gen_index < num_perms:
     1281        if not bitset_check(&SC.gen_used, perm_gen_index) and not bitset_check(&SC.gen_is_id, perm_gen_index):
     1282            unused_gens += 1
     1283    if 2*(SC.num_gens[level] + unused_gens) > SC.array_size[level]:
     1284        new_size = max(2*(SC.num_gens[level] + unused_gens), 2*SC.array_size[level])
     1285        if SC_realloc_gens(SC, level, new_size):
     1286            return 1
     1287    i = 0
     1288    for perm_gen_index from 0 <= perm_gen_index < num_perms:
     1289        if not bitset_check(&SC.gen_used, perm_gen_index) and not bitset_check(&SC.gen_is_id, perm_gen_index):
     1290            memcpy(SC.generators[level] + n*(SC.num_gens[level]+i), pi + n*perm_gen_index, n*sizeof(int))
     1291            i += 1
     1292
     1293    # Sift:
     1294    cdef int *gens = SC.generators[level]
     1295    cdef int total_gens = SC.num_gens[level] + unused_gens
     1296    cdef int section, gens_in_section
     1297    for i from 0 <= i < SC.orbit_sizes[level]:
     1298        x = SC.base_orbits[level][i]
     1299        section = 0
     1300        while section*n < total_gens:
     1301            gens_in_section = min(n, total_gens - section*n)
     1302            if SC_sift(SC, level, x, gens + n*n*section, gens_in_section, gens + n*total_gens):
     1303                return 1
     1304            section += 1
     1305    return 0
     1306
     1307cdef inline int SC_insert(StabilizerChain *SC, int level, int *pi, int num_perms):
     1308    """
     1309    Add permutations in pi to the stabilizer chain. The array pi is a sequence
     1310    of num_perms permutations, each in list representation, hence pi should be
     1311    at least length SC.degree*num_perms. There must be at most SC.degree perms.
     1312    (Simply call the function again if you want to add more.)
     1313
     1314    The variable ``level`` is used for recursion. From the outside, should be
     1315    set to zero. On the inside, used to bring the data structure up to date of
     1316    level ``level``, given that it is up to date on ``level + 1``.
     1317
     1318    Return values:
     1319    0 - No errors.
     1320    1 - Allocation failure.
     1321    """
     1322    return SC_insert_and_sift(SC, level, pi, num_perms, 1)
     1323
     1324cdef inline int SC_update_tree(StabilizerChain *SC, int level, int *pi, int num_perms):
     1325    return SC_insert_and_sift(SC, level, pi, num_perms, 0)
     1326
     1327cdef inline SC_order(StabilizerChain *SC, int i, mpz_t order):
     1328    """
     1329    Gives the order of the stabilizer of base points up to but not including the
     1330    i-th, storing it to ``order``, which must be already initialized.
     1331
     1332    To get the order of the full group, let ``i = 0``.
     1333    """
     1334    cdef int k
     1335    mpz_set_si(order, 1)
     1336    for k from i <= k < SC.base_size:
     1337        mpz_mul_si(order, order, SC.orbit_sizes[k])
     1338
     1339cdef inline bint SC_contains(StabilizerChain *SC, int level, int *pi, bint modify):
     1340    """
     1341    Test whether pi is in the level-th stabilizer.
     1342   
     1343    Assumes that pi stabilizes the first level base points.
     1344    """
     1345    cdef int b, i, j, x, n = SC.degree
     1346    cdef int *perm
     1347    if modify:
     1348        perm = pi
     1349    else:
     1350        perm = SC.perm_scratch
     1351        memcpy(perm, pi, n*sizeof(int))
     1352    for i from level <= i < SC.base_size:
     1353        b = SC.base_orbits[i][0]
     1354        x = perm[b]
     1355        if x == b: continue
     1356        if SC.parents[i][x] == -1: return 0
     1357        SC_compose_up_to_base(SC, i, x, perm)
     1358    return SC_perm_is_identity(perm, n)
     1359
     1360cdef inline SC_random_element(StabilizerChain *SC, int level, int *perm):
     1361    """
     1362    Gives a random element of the level-th stabilizer. For a random element of the
     1363    whole group, set level to 0. Must have level < SC.base_size.
     1364    """
     1365    cdef int i, x, n = SC.degree
     1366    SC_identify(perm, n)
     1367    for i from level <= i < SC.base_size:
     1368        x = SC.base_orbits[i][rand()%SC.orbit_sizes[i]]
     1369        SC_compose_up_to_base(SC, i, x, perm)
     1370
     1371cdef bint SC_is_giant(int n, int num_perms, int *perms, float p, bitset_t support):
     1372    """
     1373    Test whether the group generated by the input permutations is a giant, i.e.,
     1374    the alternating or symmetric group.
     1375   
     1376    If the group is not a giant, this routine will return False. This could also
     1377    indicate an allocation failure.
     1378   
     1379    If the group is a giant, this routine will return True with approximate
     1380    probability p. It will set `support' to the support of the group in this
     1381    case. Use bitset_len to get the size of support.
     1382   
     1383    The bitset `support' must be initialized. Must have 0 <= p < 1.
     1384   
     1385    Running time is roughly O(-ln(1-p)*n*ln(m)) where m <= n is the size of the
     1386    support of the group.
     1387    """
     1388    cdef int i, j, num_steps, m = 1, support_root
     1389    cdef unsigned long q
     1390    cdef int *gen, *perm = <int *> sage_malloc(n*sizeof(int))
     1391    cdef OrbitPartition *OP = OP_new(n)
     1392    if OP is NULL or perm is NULL:
     1393        OP_dealloc(OP)
     1394        sage_free(perm)
     1395        return False
     1396   
     1397    # giants are transitive
     1398    for j from 0 <= j < num_perms:
     1399        gen = perms + n*j
     1400        for i from 0 <= i < n:
     1401            OP_join(OP, i, gen[i])
     1402    for i from 0 <= i < n:
     1403        if OP.parent[i] == i:
     1404            if OP.size[i] != 1:
     1405                if m != 1:
     1406                    m = 1
     1407                    break
     1408                else:
     1409                    m = OP.size[i]
     1410                    support_root = i
     1411    if m == 1:
     1412        OP_dealloc(OP)
     1413        sage_free(perm)
     1414        return False
     1415    bitset_zero(support)
     1416    for i from 0 <= i < n:
     1417        if OP_find(OP, i) == support_root:
     1418            bitset_set(support, i)
     1419   
     1420    # get a bit lost in the group, so our random elements are more random:
     1421    SC_identify(perm, n)
     1422    for i from 0 <= i < 10:
     1423        SC_mult_perms(perm, perm, perms + n*(rand()%num_perms), n)
     1424   
     1425    # look for elements with cycles of prime length q, m/2 < q < m-2
     1426    num_steps = <int> ceil(-log(1-p)*log(m)/log(2))
     1427    for j from 0 <= j < num_steps:
     1428        OP_clear(OP)
     1429        for i from 0 <= i < n:
     1430            OP_join(OP, i, perm[i])
     1431        for i from 0 <= i < n:
     1432            if OP.parent[i] == i:
     1433                q = OP.size[i]
     1434                if m < 2*q and q < m-2:
     1435                    if z_isprime(q):
     1436                        sage_free(perm)
     1437                        OP_dealloc(OP)
     1438                        return True
     1439        SC_mult_perms(perm, perm, perms + n*(rand()%num_perms), n)
     1440    OP_dealloc(OP)
     1441    sage_free(perm)
     1442    return False
     1443
     1444def SC_test_list_perms(list L, int n, int limit, bint gap, bint limit_complain, bint test_contains):
     1445    """
     1446    Test that the permutation group generated by list perms in L of degree n
     1447    is of the correct order, by comparing with GAP. Don't test if the group is
     1448    of size greater than limit.
     1449   
     1450    TESTS::
     1451
     1452        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label import SC_test_list_perms
     1453        sage: limit = 10^7
     1454        sage: def test_Sn_on_m_points(n, m, gap, contains):
     1455        ...     perm1 = [1,0] + range(m)[2:]
     1456        ...     perm2 = [(i+1)%n for i in range( n )] + range(m)[n:]
     1457        ...     SC_test_list_perms([perm1, perm2], m, limit, gap, 0, contains)
     1458        sage: for i in range(2,9):
     1459        ...     test_Sn_on_m_points(i,i,1,0)
     1460        sage: for i in range(2,9):
     1461        ...     test_Sn_on_m_points(i,i,0,1)
     1462        sage: for i in range(2,9):           # long time
     1463        ...     test_Sn_on_m_points(i,i,1,1) # long time
     1464        sage: test_Sn_on_m_points(8,8,1,1)
     1465        sage: def test_stab_chain_fns_1(n, gap, contains):
     1466        ...     perm1 = sum([[2*i+1,2*i] for i in range(n)], [])
     1467        ...     perm2 = [(i+1)%(2*n) for i in range( 2*n)]
     1468        ...     SC_test_list_perms([perm1, perm2], 2*n, limit, gap, 0, contains)
     1469        sage: for n in range(1,11):
     1470        ...     test_stab_chain_fns_1(n, 1, 0)
     1471        sage: for n in range(1,11):
     1472        ...     test_stab_chain_fns_1(n, 0, 1)
     1473        sage: for n in range(1,9):              # long time
     1474        ...     test_stab_chain_fns_1(n, 1, 1)  # long time
     1475        sage: test_stab_chain_fns_1(11, 1, 1)
     1476        sage: def test_stab_chain_fns_2(n, gap, contains):
     1477        ...     perms = []
     1478        ...     for p,e in factor(n):
     1479        ...         perm1 = [(p*(i//p)) + ((i+1)%p) for i in range(n)]
     1480        ...         perms.append(perm1)
     1481        ...     SC_test_list_perms(perms, n, limit, gap, 0, contains)
     1482        sage: for n in range(2,11):
     1483        ...     test_stab_chain_fns_2(n, 1, 0)
     1484        sage: for n in range(2,11):
     1485        ...     test_stab_chain_fns_2(n, 0, 1)
     1486        sage: for n in range(2,11):            # long time
     1487        ...     test_stab_chain_fns_2(n, 1, 1) # long time
     1488        sage: test_stab_chain_fns_2(11, 1, 1)
     1489        sage: def test_stab_chain_fns_3(n, gap, contains):
     1490        ...     perm1 = [(-i)%n for i in range( n )]
     1491        ...     perm2 = [(i+1)%n for i in range( n )]
     1492        ...     SC_test_list_perms([perm1, perm2], n, limit, gap, 0, contains)
     1493        sage: for n in range(2,20):
     1494        ...     test_stab_chain_fns_3(n, 1, 0)
     1495        sage: for n in range(2,20):
     1496        ...     test_stab_chain_fns_3(n, 0, 1)
     1497        sage: for n in range(2,14):            # long time
     1498        ...     test_stab_chain_fns_3(n, 1, 1) # long time
     1499        sage: test_stab_chain_fns_3(20, 1, 1)
     1500        sage: def test_stab_chain_fns_4(n, g, gap, contains):
     1501        ...     perms = []
     1502        ...     for _ in range(g):
     1503        ...         perm = range(n)
     1504        ...         shuffle(perm)
     1505        ...         perms.append(perm)
     1506        ...     SC_test_list_perms(perms, n, limit, gap, 0, contains)
     1507        sage: for n in range(4,9):                # long time
     1508        ...     test_stab_chain_fns_4(n, 1, 1, 0) # long time
     1509        ...     test_stab_chain_fns_4(n, 2, 1, 0) # long time
     1510        ...     test_stab_chain_fns_4(n, 2, 1, 0) # long time
     1511        ...     test_stab_chain_fns_4(n, 2, 1, 0) # long time
     1512        ...     test_stab_chain_fns_4(n, 2, 1, 0) # long time
     1513        ...     test_stab_chain_fns_4(n, 3, 1, 0) # long time
     1514        sage: for n in range(4,9):
     1515        ...     test_stab_chain_fns_4(n, 1, 0, 1)
     1516        ...     for j in range(6):
     1517        ...         test_stab_chain_fns_4(n, 2, 0, 1)
     1518        ...     test_stab_chain_fns_4(n, 3, 0, 1)
     1519        sage: for n in range(4,8):                # long time
     1520        ...     test_stab_chain_fns_4(n, 1, 1, 1) # long time
     1521        ...     test_stab_chain_fns_4(n, 2, 1, 1) # long time
     1522        ...     test_stab_chain_fns_4(n, 2, 1, 1) # long time
     1523        ...     test_stab_chain_fns_4(n, 3, 1, 1) # long time
     1524        sage: test_stab_chain_fns_4(8, 2, 1, 1)
     1525        sage: def test_stab_chain_fns_5(n, gap, contains):
     1526        ...     perms = []
     1527        ...     m = n//3
     1528        ...     perm1 = range(2*m)
     1529        ...     shuffle(perm1)
     1530        ...     perm1 += range(2*m,n)
     1531        ...     perm2 = range(m,n)
     1532        ...     shuffle(perm2)
     1533        ...     perm2 = range(m) + perm2
     1534        ...     SC_test_list_perms([perm1, perm2], n, limit, gap, 0, contains)
     1535        sage: for n in [4..9]:                     # long time
     1536        ...     for _ in range(2):                 # long time
     1537        ...         test_stab_chain_fns_5(n, 1, 0) # long time
     1538        sage: for n in [4..8]:                     # long time
     1539        ...     test_stab_chain_fns_5(n, 0, 1)     # long time
     1540        sage: for n in [4..9]:                     # long time
     1541        ...     test_stab_chain_fns_5(n, 1, 1)     # long time
     1542        sage: def random_perm(x):
     1543        ...     shuffle(x)
     1544        ...     return x
     1545        sage: def test_stab_chain_fns_6(m,n,k, gap, contains):
     1546        ...     perms = []
     1547        ...     for i in range(k):
     1548        ...         perm = sum([random_perm(range(i*(n//m),min(n,(i+1)*(n//m)))) for i in range(m)], [])
     1549        ...         perms.append(perm)
     1550        ...     SC_test_list_perms(perms, m*(n//m), limit, gap, 0, contains)
     1551        sage: for m in range(2,9):                         # long time
     1552        ...     for n in range(m,3*m):                     # long time
     1553        ...         for k in range(1,3):                   # long time
     1554        ...             test_stab_chain_fns_6(m,n,k, 1, 0) # long time
     1555        sage: for m in range(2,10):
     1556        ...     for n in range(m,4*m):
     1557        ...         for k in range(1,3):
     1558        ...             test_stab_chain_fns_6(m,n,k, 0, 1)
     1559        sage: test_stab_chain_fns_6(10,20,2, 1, 1)
     1560        sage: test_stab_chain_fns_6(8,16,2, 1, 1)
     1561        sage: test_stab_chain_fns_6(6,36,2, 1, 1)
     1562        sage: test_stab_chain_fns_6(4,40,3, 1, 1)
     1563        sage: test_stab_chain_fns_6(4,40,2, 1, 1)
     1564        sage: def test_stab_chain_fns_7(n, cop, gap, contains):
     1565        ...     perms = []
     1566        ...     for i in range(0,n//2,2):
     1567        ...         p = range(n)
     1568        ...         p[i] = i+1
     1569        ...         p[i+1] = i
     1570        ...     if cop:
     1571        ...         perms.append([c for c in p])
     1572        ...     else:
     1573        ...         perms.append(p)
     1574        ...     SC_test_list_perms(perms, n, limit, gap, 0, contains)
     1575        sage: for n in [6..14]:
     1576        ...     test_stab_chain_fns_7(n, 1, 1, 0)
     1577        ...     test_stab_chain_fns_7(n, 0, 1, 0)
     1578        sage: for n in [6..30]:
     1579        ...     test_stab_chain_fns_7(n, 1, 0, 1)
     1580        ...     test_stab_chain_fns_7(n, 0, 0, 1)
     1581        sage: for n in [6..14]:                   # long time
     1582        ...     test_stab_chain_fns_7(n, 1, 1, 1) # long time
     1583        ...     test_stab_chain_fns_7(n, 0, 1, 1) # long time
     1584        sage: test_stab_chain_fns_7(20, 1, 1, 1)
     1585        sage: test_stab_chain_fns_7(20, 0, 1, 1)
     1586
     1587    """
     1588    if gap:
     1589        from sage.all import PermutationGroup, PermutationGroupElement, shuffle
     1590    cdef StabilizerChain *SC, *SCC, *SCCC, *SC_nb
     1591    cdef Integer order, order2
     1592    cdef int i, j, m, SC_says
     1593    cdef bitset_t giant_support
     1594    if gap:
     1595        G = PermutationGroup([[i+1 for i in p] for p in L])
     1596        if G.order() > limit:
     1597            if limit_complain: print 'TOO BIG'
     1598            return
     1599    SC = SC_new(n)
     1600    cdef int *perm = <int *>sage_malloc(n * (len(L)+3) * sizeof(int))
     1601    try:
     1602        bitset_init(giant_support, n)
     1603    except MemoryError:
     1604        sage_free(perm)
     1605        SC_dealloc(SC)
     1606        raise MemoryError
     1607    if perm is NULL or SC is NULL:
     1608        bitset_free(giant_support)
     1609        sage_free(perm)
     1610        SC_dealloc(SC)
     1611        raise MemoryError
     1612    cdef int *perm2 = perm +   n
     1613    cdef int *perm3 = perm + 2*n
     1614    for Lperm in L:
     1615        for i from 0 <= i < n:
     1616            perm[i] = Lperm[i]
     1617        if SC_insert(SC, 0, perm, 1):
     1618            bitset_free(giant_support)
     1619            sage_free(perm)
     1620            SC_dealloc(SC)
     1621            raise MemoryError
     1622    SCC = SC_copy(SC, n)
     1623    SCCC = SC_insert_base_point(SC, 0, n-1)
     1624    for i from 0 <= i < n:
     1625        perm[i] = n-i-1
     1626    SC_nb = SC_new_base(SC, perm, n)
     1627    if SCC is NULL or SCCC is NULL or SC_nb is NULL:
     1628        bitset_free(giant_support)
     1629        sage_free(perm)
     1630        SC_dealloc(SC)
     1631        SC_dealloc(SCC)
     1632        SC_dealloc(SCCC)
     1633        SC_dealloc(SC_nb)
     1634        raise MemoryError
     1635    giant = False
     1636    try:
     1637        order = Integer(0)
     1638        SC_order(SC,0,order.value)
     1639        j = 0
     1640        for Lperm in L:
     1641            for i from 0 <= i < n:
     1642                perm[n*j + i] = Lperm[i]
     1643            j += 1
     1644        if SC_is_giant(n, len(L), perm, 0.9, giant_support):
     1645            giant = True
     1646            m = bitset_len(giant_support)
     1647            from sage.rings.arith import factorial
     1648            if not (order == factorial(m) or order == factorial(m)/2):
     1649                print "SC_is_giant failed: %s %s"%(str(L), order)
     1650                raise AssertionError
     1651            if order == factorial(n):
     1652                SC_dealloc(SC)
     1653                SC = SC_symmetric_group(n)
     1654                SC_order(SC,0,order.value)
     1655                if not order == factorial(n):
     1656                    print "SC_symmetric_group failed: %s %s"%(str(L), order)
     1657                    raise AssertionError
     1658            elif order == factorial(n)/2:
     1659                SC_dealloc(SC)
     1660                SC = SC_alternating_group(n)
     1661                SC_order(SC,0,order.value)
     1662                if not order == factorial(n)/2:
     1663                    print "SC_alternating_group failed: %s %s"%(str(L), order)
     1664                    raise AssertionError
     1665        order2 = Integer(0)
     1666        SC_order(SCC,0,order2.value)
     1667        if order != order2:
     1668            print "FAIL", L
     1669            print 'SC_copy(n) does not agree with order of original', order, order2
     1670            raise AssertionError
     1671        SC_order(SCCC,0,order2.value)
     1672        if order != order2:
     1673            print "FAIL", L
     1674            print 'does not agree with order of inserted base point', order, order2
     1675            raise AssertionError
     1676        SC_order(SC_nb,0,order2.value)
     1677        if order != order2:
     1678            print "FAIL", L
     1679            print 'does not agree with order of new base', order, order2
     1680            raise AssertionError
     1681        if test_contains:
     1682            for i from 0 <= i < 3:
     1683                SC_random_element(SC, 0, perm)
     1684                if not SC_contains(SC, 0, perm, 0):
     1685                    print "FAIL", L
     1686                    print 'element', [perm[ii] for ii in range(n)]
     1687                    print 'SC_random_element says it is an element, SC_contains(modify=0) does not'
     1688                    raise AssertionError
     1689                if not SC_contains(SC, 0, perm, 1):
     1690                    print "FAIL", L
     1691                    print 'element', [perm[ii] for ii in range(n)]
     1692                    print 'SC_random_element says it is an element, SC_contains(modify=1) does not'
     1693                    raise AssertionError
     1694                if not SC_contains(SCC, 0, perm, 0):
     1695                    print "FAIL", L
     1696                    print 'element', [perm[ii] for ii in range(n)]
     1697                    print 'SC_random_element says it is an element, SC_contains(modify=0) does not on copy'
     1698                    raise AssertionError
     1699                if not SC_contains(SCC, 0, perm, 1):
     1700                    print "FAIL", L
     1701                    print 'element', [perm[ii] for ii in range(n)]
     1702                    print 'SC_random_element says it is an element, SC_contains(modify=1) does not on copy'
     1703                    raise AssertionError
     1704                if not SC_contains(SCCC, 0, perm, 0):
     1705                    print "FAIL", L
     1706                    print 'element', [perm[ii] for ii in range(n)]
     1707                    print 'SC_random_element says it is an element, SC_contains(modify=0) does not on inserted base point'
     1708                    raise AssertionError
     1709                if not SC_contains(SCCC, 0, perm, 1):
     1710                    print "FAIL", L
     1711                    print 'element', [perm[ii] for ii in range(n)]
     1712                    print 'SC_random_element says it is an element, SC_contains(modify=1) does not on inserted base point'
     1713                    raise AssertionError
     1714                if not SC_contains(SC_nb, 0, perm, 0):
     1715                    print "FAIL", L
     1716                    print 'element', [perm[ii] for ii in range(n)]
     1717                    print 'SC_random_element says it is an element, SC_contains(modify=0) does not on new base'
     1718                    raise AssertionError
     1719                if not SC_contains(SC_nb, 0, perm, 1):
     1720                    print "FAIL", L
     1721                    print 'element', [perm[ii] for ii in range(n)]
     1722                    print 'SC_random_element says it is an element, SC_contains(modify=1) does not on new base'
     1723                    raise AssertionError
     1724                SC_random_element(SCC, 0, perm2)
     1725                if not SC_contains(SC, 0, perm2, 0):
     1726                    print "FAIL", L
     1727                    print 'element', [perm[ii] for ii in range(n)]
     1728                    print 'SC_random_element says it is an element of copy, SC_contains(modify=0) does not'
     1729                    raise AssertionError
     1730                if not SC_contains(SC, 0, perm2, 1):
     1731                    print "FAIL", L
     1732                    print 'element', [perm[ii] for ii in range(n)]
     1733                    print 'SC_random_element says it is an element of copy, SC_contains(modify=1) does not'
     1734                    raise AssertionError
     1735                SC_random_element(SCCC, 0, perm3)
     1736                if not SC_contains(SC, 0, perm3, 0):
     1737                    print "FAIL", L
     1738                    print 'element', [perm[ii] for ii in range(n)]
     1739                    print 'SC_random_element says it is an element of inserted base point, SC_contains(modify=0) does not'
     1740                    raise AssertionError
     1741                if not SC_contains(SC, 0, perm3, 1):
     1742                    print "FAIL", L
     1743                    print 'element', [perm[ii] for ii in range(n)]
     1744                    print 'SC_random_element says it is an element of inserted base point, SC_contains(modify=1) does not'
     1745                    raise AssertionError
     1746                SC_random_element(SC_nb, 0, perm3)
     1747                if not SC_contains(SC, 0, perm3, 0):
     1748                    print "FAIL", L
     1749                    print 'element', [perm[ii] for ii in range(n)]
     1750                    print 'SC_random_element says it is an element of new base, SC_contains(modify=0) does not'
     1751                    raise AssertionError
     1752                if not SC_contains(SC, 0, perm3, 1):
     1753                    print "FAIL", L
     1754                    print 'element', [perm[ii] for ii in range(n)]
     1755                    print 'SC_random_element says it is an element of new base, SC_contains(modify=1) does not'
     1756                    raise AssertionError
     1757        if gap:
     1758            order = Integer(0)
     1759            SC_order(SC,0,order.value)
     1760            j = 0
     1761            for Lperm in L:
     1762                for i from 0 <= i < n:
     1763                    perm[n*j + i] = Lperm[i]
     1764                j += 1
     1765            if SC_is_giant(n, len(L), perm, 0.9, giant_support):
     1766                from sage.rings.arith import factorial
     1767                m = bitset_len(giant_support)
     1768                if order != factorial(m) and order != factorial(m)/2:
     1769                    print "SC_is_giant failed: %s %s"%(str(L), order)
     1770                    raise AssertionError
     1771            if order != G.order():
     1772                print "FAIL", L
     1773                print 'order was computed to be', order
     1774                print 'GAP says it is', G.order()
     1775                raise AssertionError
     1776            if test_contains:
     1777                for i from 0 <= i < 3:
     1778                    permy = G.random_element().list()
     1779                    for j from 0 <= j < n:
     1780                        perm[j] = permy[j]-1
     1781                    if not SC_contains(SC, 0, perm, 0):
     1782                        print "FAIL", L
     1783                        print 'element', permy
     1784                        print 'GAP says it is an element, SC_contains(modify=0) does not'
     1785                        raise AssertionError
     1786                    if not SC_contains(SC, 0, perm, 1):
     1787                        print "FAIL", L
     1788                        print 'element', permy
     1789                        print 'GAP says it is an element, SC_contains(modify=1) does not'
     1790                        raise AssertionError
     1791                    permy = range(1,n+1)
     1792                    shuffle(permy)
     1793                    gap_says = (PermutationGroupElement(permy) in G)
     1794                    for j from 0 <= j < n:
     1795                        perm[j] = permy[j]-1
     1796                    SC_says = SC_contains(SC, 0, perm, 0)
     1797                    if bool(SC_says) != bool(gap_says):
     1798                        print "FAIL", L
     1799                        print 'element', permy
     1800                        print 'GAP says %d, SC_contains(modify=0) says %d'%(gap_says, SC_says)
     1801                        raise AssertionError
     1802                    SC_says = SC_contains(SC, 0, perm, 1)
     1803                    if bool(SC_says) != bool(gap_says):
     1804                        print "FAIL", L
     1805                        print 'element', permy
     1806                        print 'GAP says %d, SC_contains(modify=0) says %d'%(gap_says, SC_says)
     1807                        raise AssertionError
     1808                    SC_random_element(SC, 0, perm)
     1809                    for j from 0 <= j < n:
     1810                        permy[j] = perm[j]+1
     1811                    gap_says = (PermutationGroupElement(permy) in G)
     1812                    if not SC_contains(SC, 0, perm, 0):
     1813                        print "FAIL", L
     1814                        print 'element', permy
     1815                        print 'random element not contained in group, modify=false'
     1816                        raise AssertionError
     1817                    if not SC_contains(SC, 0, perm, 1):
     1818                        print "FAIL", L
     1819                        print 'element', permy
     1820                        print 'random element not contained in group, modify=true'
     1821                        raise AssertionError
     1822                    if not gap_says:
     1823                        print "FAIL", L
     1824                        print 'element', permy
     1825                        print 'random element not contained in group, according to GAP'
     1826                        raise AssertionError
     1827    except AssertionError:
     1828        bitset_free(giant_support)
     1829        sage_free(perm)
     1830        SC_dealloc(SC)
     1831        SC_dealloc(SCC)
     1832        SC_dealloc(SCCC)
     1833        SC_dealloc(SC_nb)
     1834        if giant:
     1835            print "detected giant!"
     1836        raise AssertionError
     1837    bitset_free(giant_support)
     1838    sage_free(perm)
     1839    SC_dealloc(SC)
     1840    SC_dealloc(SCC)
     1841    SC_dealloc(SCCC)
     1842    SC_dealloc(SC_nb)
     1843
    5991844# Functions
    6001845
    601 cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S,
     1846cdef int sort_by_function(PartitionStack *PS, int start, int *degrees):
     1847    """
     1848    A simple counting sort, given the degrees of vertices to a certain cell.
     1849   
     1850    INPUT:
     1851    PS -- the partition stack to be checked
     1852    start -- beginning index of the cell to be sorted
     1853    degrees -- the values to be sorted by, must have extra scratch space for a
     1854        total of 3*n+1
     1855
     1856    """
     1857    cdef int n = PS.degree
     1858    cdef int i, j, max, max_location
     1859    cdef int *counts = degrees + n, *output = degrees + 2*n + 1
     1860    for i from 0 <= i <= n:
     1861        counts[i] = 0
     1862    i = 0
     1863    while PS.levels[i+start] > PS.depth:
     1864        counts[degrees[i]] += 1
     1865        i += 1
     1866    counts[degrees[i]] += 1
     1867    # i+start is the right endpoint of the cell now
     1868    max = counts[0]
     1869    max_location = 0
     1870    for j from 0 < j <= n:
     1871        if counts[j] > max:
     1872            max = counts[j]
     1873            max_location = j
     1874        counts[j] += counts[j - 1]
     1875    for j from i >= j >= 0:
     1876        counts[degrees[j]] -= 1
     1877        output[counts[degrees[j]]] = PS.entries[start+j]
     1878    max_location = counts[max_location]+start
     1879    for j from 0 <= j <= i:
     1880        PS.entries[start+j] = output[j]
     1881    j = 1
     1882    while j <= n and counts[j] <= i:
     1883        if counts[j] > 0:
     1884            PS.levels[start + counts[j] - 1] = PS.depth
     1885        PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1)
     1886        j += 1
     1887    return max_location
     1888
     1889cdef inline int split_point_and_refine(PartitionStack *PS, int v, void *S,
    6021890    int (*refine_and_return_invariant)\
    603          (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
     1891         (PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len),
    6041892    int *cells_to_refine_by):
    6051893    """
    6061894    Make the partition stack one longer by copying the last partition in the
     
    6131901    S -- the structure
    6141902    refine_and_return_invariant -- the refinement function provided
    6151903    cells_to_refine_by -- an array, contents ignored
     1904    group -- the containing group, NULL for full S_n
     1905    perm_stack -- represents a partial traversal decomposition for group
    6161906
    6171907    """
    6181908    PS.depth += 1
     
    6201910    cells_to_refine_by[0] = PS_split_point(PS, v)
    6211911    return refine_and_return_invariant(PS, S, cells_to_refine_by, 1)
    6221912
     1913cdef inline update_perm_stack(StabilizerChain *group, int level, int point,
     1914    int *perm_stack):
     1915    """
     1916    Ensure that perm_stack[level] is gamma_0^{-1}...gamma_{level-1}^{-1}, where
     1917    each gamma_i represents the coset representative at the ith level determined
     1918    by our current position in the search tree.
     1919   
     1920    For internal use within the automorphism group, canonical label and double
     1921    coset functions, to be called after refinement (level = depth after refinement).
     1922    """
     1923    cdef int n = group.degree
     1924    memcpy(perm_stack + n*level, perm_stack + n*(level-1), n*sizeof(int))
     1925    SC_compose_up_to_base(group, level-1, perm_stack[n*(level-1) + point], perm_stack + n*level)
     1926
     1927cdef int refine_by_orbits(PartitionStack *PS, StabilizerChain *SC, int *perm_stack, int *cells_to_refine_by, int *ctrb_len):
     1928    """
     1929    Given a stabilizer chain SC, refine the partition stack PS so that each cell
     1930    contains elements from at most one orbit, and sort the refined cells by
     1931    their minimal cell representatives in the orbit of the group.
     1932    """
     1933    cdef int start, end, level, gen_index, i, j, k, x, n = SC.degree
     1934    cdef int *gen, *min_cell_reps = SC.perm_scratch, *perm = perm_stack + n*PS.depth
     1935    cdef OrbitPartition *OP = SC.OP_scratch
     1936    cdef int invariant = 1
     1937    OP_clear(OP)
     1938    for level from PS.depth <= level < SC.base_size:
     1939        for gen_index from 0 <= gen_index < SC.num_gens[level]:
     1940            gen = SC.generators[level] + gen_index*n
     1941            for i from 0 <= i < n:
     1942                OP_join(OP, i, gen[i])
     1943    start = 0
     1944    while start < n:
     1945        i = 0
     1946        while 1:
     1947            x = PS.entries[start+i]
     1948            x = perm[x]
     1949            min_cell_reps[i] = OP.mcr[OP_find(OP, x)]
     1950            if PS.levels[start+i] <= PS.depth:
     1951                break
     1952            i += 1
     1953        invariant += sort_by_function(PS, start, min_cell_reps)
     1954        invariant += i
     1955        # update cells_to_refine_by
     1956        k = start
     1957        for j from start <= j < start+i:
     1958            if PS.levels[j] == PS.depth:
     1959                cells_to_refine_by[ctrb_len[0]] = k
     1960                ctrb_len[0] += 1
     1961                k = j+1
     1962        start += i+1
     1963   
     1964    return invariant
     1965
     1966cdef inline int split_point_and_refine_by_orbits(PartitionStack *PS, int v,
     1967    void *S, int (*refine_and_return_invariant)\
     1968         (PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len),
     1969    int *cells_to_refine_by, StabilizerChain *SC, int *perm_stack):
     1970    """ """
     1971    PS.depth += 1
     1972    PS_clear(PS)
     1973    cells_to_refine_by[0] = PS_split_point(PS, v)
     1974    update_perm_stack(SC, PS.depth, v, perm_stack)
     1975    return refine_also_by_orbits(PS, S, refine_and_return_invariant, cells_to_refine_by, 1, SC, perm_stack)
     1976
     1977cdef inline int refine_also_by_orbits(PartitionStack *PS, void *S,
     1978    int (*refine_and_return_invariant)\
     1979         (PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len),
     1980    int *cells_to_refine_by, int ctrb_len, StabilizerChain *SC, int *perm_stack):
     1981    """ """
     1982    cdef int inv, n = PS.degree
     1983    inv  = refine_by_orbits(PS, SC, perm_stack, cells_to_refine_by, &ctrb_len)
     1984    inv += refine_and_return_invariant(PS, S, cells_to_refine_by, ctrb_len)
     1985    return inv
     1986
     1987cdef int compute_relabeling(StabilizerChain *group, StabilizerChain *scratch_group,
     1988    int *permutation, int *relabeling):
     1989    """
     1990    Technically, compute the INVERSE of the relabeling
     1991    """
     1992    cdef int i, j, x, y, m, n = group.degree
     1993    cdef int *scratch = group.perm_scratch
     1994    if SC_new_base_nomalloc(scratch_group, group, permutation, n):
     1995        return 1
     1996    group = scratch_group
     1997    SC_identify(relabeling, n)
     1998    for i from 0 <= i < n:
     1999        m = n
     2000        for j from 0 <= j < group.orbit_sizes[i]:
     2001            x = relabeling[group.base_orbits[i][j]]
     2002            if x < m:
     2003                m = x
     2004                y = group.base_orbits[i][j]
     2005        SC_invert_perm(scratch, relabeling, n)
     2006        SC_compose_up_to_base(group, i, y, scratch) # doesn't use group.perm_scratch
     2007        SC_invert_perm(relabeling, scratch, n)
     2008    SC_invert_perm(scratch, relabeling, n)
     2009    memcpy(relabeling, scratch, n*sizeof(int))
     2010    return 0
     2011
     2012
     2013
  • sage/groups/perm_gps/partn_ref/double_coset.pxd

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/double_coset.pxd
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#
    55# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    66#                         http://www.gnu.org/licenses/
     
    1212
    1313from sage.rings.integer cimport Integer
    1414
    15 cdef int *double_coset( object, object, int **, int *, int,
    16     bint (*)(PartitionStack *, object),
    17     int (*)(PartitionStack *, object, int *, int),
    18     int (*)(int *, int *, object, object))
     15cdef int *double_coset( void *, void *, PartitionStack *, int *, int,
     16    bint (*)(PartitionStack *, void *),
     17    int (*)(PartitionStack *, void *, int *, int),
     18    int (*)(int *, int *, void *, void *),
     19    StabilizerChain *)
    1920
    2021
    2122
  • sage/groups/perm_gps/partn_ref/double_coset.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/double_coset.pyx
    a b  
    1919   
    2020    Signature:
    2121   
    22     \code{int refine_and_return_invariant(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len)}
     22    \code{int refine_and_return_invariant(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len)}
    2323   
    2424    This function should split up cells in the partition at the top of the
    2525    partition stack in such a way that any automorphism that respects the
     
    4141   
    4242    Signature:
    4343   
    44     \code{int compare_structures(int *gamma_1, int *gamma_2, object S1, object S2)}
     44    \code{int compare_structures(int *gamma_1, int *gamma_2, void *S1, void *S2)}
    4545   
    4646    This function must implement a total ordering on the set of objects of fixed
    4747    order. Return:
     
    6060   
    6161    Signature:
    6262   
    63     \code{bint all_children_are_equivalent(PartitionStack *PS, object S)}
     63    \code{bint all_children_are_equivalent(PartitionStack *PS, void *S)}
    6464   
    6565    This function must return False unless it is the case that each discrete
    6666    partition finer than the top of the partition stack is equivalent to the
     
    8383"""
    8484
    8585#*****************************************************************************
    86 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     86#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    8787#
    8888# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    8989#                         http://www.gnu.org/licenses/
     
    9696    elif a == b: return 0
    9797    else: return 1
    9898
    99 cdef inline bint stacks_are_equivalent(PartitionStack *PS1, PartitionStack *PS2):
    100     cdef int i, j, depth = min(PS1.depth, PS2.depth)
    101     for i from 0 <= i < PS1.degree:
    102         j = cmp(PS1.levels[i], PS2.levels[i])
    103         if j == 0: continue
    104         if ( (j < 0 and PS1.levels[i] <= depth and PS2.levels[i] > depth)
    105             or (j > 0 and PS2.levels[i] <= depth and PS1.levels[i] > depth) ):
    106             return 0
    107     return 1
    108 
    10999# Functions
    110100
    111 cdef int *double_coset(object S1, object S2, int **partition1, int *ordering2,
    112     int n, bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
     101cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, void *S):
     102    return 0
     103
     104cdef int refine_and_return_invariant_trivial(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len):
     105    return 0
     106
     107cdef int compare_perms(int *gamma_1, int *gamma_2, void *S1, void *S2):
     108    cdef list MS1 = <list> S1
     109    cdef list MS2 = <list> S2
     110    cdef int i, j
     111    for i from 0 <= i < len(MS1):
     112        j = cmp(MS1[gamma_1[i]], MS2[gamma_2[i]])
     113        if j != 0: return j
     114    return 0
     115
     116def coset_eq(list perm1=[0,1,2,3,4,5], list perm2=[1,2,3,4,5,0], list gens=[[1,2,3,4,5,0]]):
     117    """
     118    Given a group G generated by the given generators, tests whether the given
     119    permutations are in the same right coset of G. Tests nontrivial input group
     120    when using double_coset. If they are, return an element g so that
     121    g.perm1 = perm2 (composing left to right).
     122   
     123    TESTS::
     124
     125        sage: from sage.groups.perm_gps.partn_ref.double_coset import coset_eq
     126        sage: coset_eq()
     127        [5, 0, 1, 2, 3, 4]
     128        sage: gens = [[1,2,3,0]]
     129        sage: reps = [[0,1,2,3]]
     130        sage: for p in SymmetricGroup(4):
     131        ...     p = [a-1 for a in p.list()]
     132        ...     found = False
     133        ...     for r in reps:
     134        ...         if coset_eq(p, r, gens):
     135        ...             found = True
     136        ...             break
     137        ...     if not found:
     138        ...         reps.append(p)
     139        sage: len(reps)
     140        6
     141        sage: gens = [[1,0,2,3],[0,1,3,2]]
     142        sage: reps = [[0,1,2,3]]
     143        sage: for p in SymmetricGroup(4):
     144        ...     p = [a-1 for a in p.list()]
     145        ...     found = False
     146        ...     for r in reps:
     147        ...         if coset_eq(p, r, gens):
     148        ...             found = True
     149        ...             break
     150        ...     if not found:
     151        ...         reps.append(p)
     152        sage: len(reps)
     153        6
     154        sage: gens = [[1,2,0,3]]
     155        sage: reps = [[0,1,2,3]]
     156        sage: for p in SymmetricGroup(4):
     157        ...     p = [a-1 for a in p.list()]
     158        ...     found = False
     159        ...     for r in reps:
     160        ...         if coset_eq(p, r, gens):
     161        ...             found = True
     162        ...             break
     163        ...     if not found:
     164        ...         reps.append(p)
     165        sage: len(reps)
     166        8
     167
     168    """
     169    cdef int i, n = len(perm1)
     170    assert all(len(g) == n for g in gens+[perm2])
     171    cdef PartitionStack *part = PS_new(n, 1)
     172    cdef int *c_perm = <int *> sage_malloc(n * sizeof(int))
     173    cdef StabilizerChain *group = SC_new(n, 1)
     174    if part is NULL or c_perm is NULL or group is NULL:
     175        sage_free(c_perm)
     176        PS_dealloc(part)
     177        SC_dealloc(group)
     178        raise MemoryError
     179    for g in gens:
     180        for i from 0 <= i < n:
     181            c_perm[i] = g[i]
     182        SC_insert(group, 0, c_perm, 1)
     183    for i from 0 <= i < n:
     184        c_perm[i] = i
     185    cdef int *output = double_coset(<void *> perm1, <void *> perm2, part, c_perm, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_perms, group)
     186    sage_free(c_perm)
     187    PS_dealloc(part)
     188    SC_dealloc(group)
     189    if output is NULL:
     190        return False
     191    else:
     192        x = [output[i] for i from 0 <= i < n]
     193        sage_free(output)
     194        return x
     195
     196cdef int *double_coset(void *S1, void *S2, PartitionStack *partition1, int *ordering2,
     197    int n, bint (*all_children_are_equivalent)(PartitionStack *PS, void *S),
    113198    int (*refine_and_return_invariant)\
    114          (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
    115     int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2) ):
     199         (PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len),
     200    int (*compare_structures)(int *gamma_1, int *gamma_2, void *S1, void *S2),
     201    StabilizerChain *input_group):
    116202    """
    117203    Traverse the search space for double coset calculation.
    118204   
    119205    INPUT:
    120206    S1, S2 -- pointers to the structures
    121     partition1 -- array representing a partition of the points of S1
     207    partition1 -- PartitionStack of depth 0 and degree n,
     208        whose first partition is of the points of S1
    122209    ordering2 -- an ordering of the points of S2 representing a second partition
    123210    n -- the number of points (points are assumed to be 0,1,...,n-1)
    124211    all_children_are_equivalent -- pointer to a function
     
    143230        OUTPUT:
    144231        int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp),
    145232            such that the set of all structures is well-ordered
     233   
     234    NOTE:
     235    The partition ``partition1`` and the resulting partition from ``ordering2``
     236    *must* satisfy the property that in each cell, the smallest element occurs
     237    first!
     238   
    146239    OUTPUT:
    147240    If S1 and S2 are isomorphic, a pointer to an integer array representing an
    148241    isomorphism. Otherwise, a NULL pointer.
     
    168261    cdef bitset_t vertices_have_been_reduced
    169262    cdef int *permutation, *id_perm, *cells_to_refine_by
    170263    cdef int *vertices_determining_current_stack
     264    cdef int *perm_stack
     265    cdef StabilizerChain *group, *old_group, *tmp_gp
    171266   
    172267    cdef int *output
    173268   
     
    178273    if n == 0:
    179274        return NULL
    180275   
    181     indicators = <int *> sage_malloc(n * sizeof(int))
    182    
    183     fixed_points_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
    184     minimal_cell_reps_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
    185    
    186     vertices_to_split = <bitset_t *> sage_malloc( n * sizeof(bitset_t) )
    187     permutation = <int *> sage_malloc( n * sizeof(int) )
    188     id_perm = <int *> sage_malloc( n * sizeof(int) )
    189     cells_to_refine_by = <int *> sage_malloc( n * sizeof(int) )
    190     vertices_determining_current_stack = <int *> sage_malloc( n * sizeof(int) )
    191 
     276    cdef int *int_array = <int *> sage_malloc( 5*n * sizeof(int) )
    192277    output = <int *> sage_malloc( n * sizeof(int) )
    193    
     278    if input_group is not NULL:
     279        perm_stack = <int *> sage_malloc( n*n * sizeof(int) )
     280        group = SC_copy(input_group, n)
     281        old_group = SC_new(n)
     282        if perm_stack is NULL or group is NULL or old_group is NULL:
     283            mem_err = 1
     284        else:
     285            SC_identify(perm_stack, n)
     286    cdef bitset_t *bitset_array = <bitset_t *> sage_malloc( (n+2*len_of_fp_and_mcr) * sizeof(bitset_t) )
     287    left_ps = partition1
    194288    current_ps = PS_new(n, 0)
    195     left_ps = PS_new(n, 0)
    196289    orbits_of_subgroup = OP_new(n)
    197290   
    198     # Check for allocation failures:
    199     if indicators                         is NULL or \
    200        fixed_points_of_generators         is NULL or \
    201        minimal_cell_reps_of_generators    is NULL or \
    202        vertices_to_split                  is NULL or \
    203        permutation                        is NULL or \
    204        id_perm                            is NULL or \
    205        cells_to_refine_by                 is NULL or \
    206        vertices_determining_current_stack is NULL or \
    207        current_ps                         is NULL or \
    208        orbits_of_subgroup                 is NULL or \
    209        output                             is NULL:
    210         if indicators                         is not NULL:
    211             sage_free(indicators)
    212         if fixed_points_of_generators         is not NULL:
    213             sage_free(fixed_points_of_generators)
    214         if minimal_cell_reps_of_generators    is not NULL:
    215             sage_free(minimal_cell_reps_of_generators)
    216         if vertices_to_split                  is not NULL:
    217             sage_free(vertices_to_split)
    218         if permutation                        is not NULL:
    219             sage_free(permutation)
    220         if id_perm                            is not NULL:
    221             sage_free(id_perm)
    222         if cells_to_refine_by                 is not NULL:
    223             sage_free(cells_to_refine_by)
    224         if vertices_determining_current_stack is not NULL:
    225             sage_free(vertices_determining_current_stack)
    226         if current_ps is not NULL:
    227             PS_dealloc(current_ps)
    228         if orbits_of_subgroup is not NULL:
    229             OP_dealloc(orbits_of_subgroup)
    230         if output is not NULL:
    231             sage_free(output)
    232         raise MemoryError
     291    if int_array  is NULL or output is NULL or bitset_array is NULL or \
     292       current_ps is NULL or orbits_of_subgroup is NULL:
     293        mem_err = 1
    233294   
    234     # Initialize bitsets, checking for allocation failures:
    235     cdef bint succeeded = 1
    236     for i from 0 <= i < len_of_fp_and_mcr:
     295    if not mem_err:
     296        indicators                         = int_array
     297        permutation                        = int_array +   n
     298        id_perm                            = int_array + 2*n
     299        cells_to_refine_by                 = int_array + 3*n
     300        vertices_determining_current_stack = int_array + 4*n
     301       
     302        fixed_points_of_generators      = bitset_array
     303        minimal_cell_reps_of_generators = bitset_array + len_of_fp_and_mcr
     304        vertices_to_split               = bitset_array + 2*len_of_fp_and_mcr
    237305        try:
    238             bitset_init(fixed_points_of_generators[i], n)
    239         except MemoryError:
    240             succeeded = 0
    241             for j from 0 <= j < i:
    242                 bitset_free(fixed_points_of_generators[j])
    243                 bitset_free(minimal_cell_reps_of_generators[j])
    244             break
    245         try:
    246             bitset_init(minimal_cell_reps_of_generators[i], n)
    247         except MemoryError:
    248             succeeded = 0
    249             for j from 0 <= j < i:
    250                 bitset_free(fixed_points_of_generators[j])
    251                 bitset_free(minimal_cell_reps_of_generators[j])
    252             bitset_free(fixed_points_of_generators[i])
    253             break
    254     if succeeded:
    255         for i from 0 <= i < n:
    256             try:
    257                 bitset_init(vertices_to_split[i], n)
    258             except MemoryError:
    259                 succeeded = 0
    260                 for j from 0 <= j < i:
    261                     bitset_free(vertices_to_split[j])
    262                 for j from 0 <= j < len_of_fp_and_mcr:
    263                     bitset_free(fixed_points_of_generators[j])
    264                     bitset_free(minimal_cell_reps_of_generators[j])
    265                 break
    266     if succeeded:
    267         try:
     306            for i from 0 <= i < n+2*len_of_fp_and_mcr:
     307                bitset_init(bitset_array[i], n)
    268308            bitset_init(vertices_have_been_reduced, n)
    269309        except MemoryError:
    270             succeeded = 0
    271             for j from 0 <= j < n:
    272                 bitset_free(vertices_to_split[j])
    273             for j from 0 <= j < len_of_fp_and_mcr:
    274                 bitset_free(fixed_points_of_generators[j])
    275                 bitset_free(minimal_cell_reps_of_generators[j])
    276     if not succeeded:
    277         sage_free(indicators)
    278         sage_free(fixed_points_of_generators)
    279         sage_free(minimal_cell_reps_of_generators)
    280         sage_free(vertices_to_split)
    281         sage_free(permutation)
    282         sage_free(id_perm)
    283         sage_free(cells_to_refine_by)
    284         sage_free(vertices_determining_current_stack)
    285         PS_dealloc(current_ps)
    286         PS_dealloc(left_ps)
    287         OP_dealloc(orbits_of_subgroup)
    288         raise MemoryError
    289    
    290     bitset_zero(vertices_have_been_reduced)
    291    
    292     # set up the identity permutation
    293     for i from 0 <= i < n:
    294         id_perm[i] = i
    295     if ordering2 is NULL:
    296         ordering2 = id_perm
    297    
    298     # Copy data of partition to left_ps, and
    299     # ordering of that partition to current_ps.
    300     i = 0
    301     j = 0
    302     while partition1[i] is not NULL:
    303         k = 0
    304         while partition1[i][k] != -1:
    305             left_ps.entries[j+k] = partition1[i][k]
    306             left_ps.levels[j+k] = n
    307             current_ps.entries[j+k] = ordering2[j+k]
    308             current_ps.levels[j+k] = n
    309             k += 1
    310         left_ps.levels[j+k-1] = 0
    311         current_ps.levels[j+k-1] = 0
    312         PS_move_min_to_front(current_ps, j, j+k-1)
    313         j += k
    314         i += 1
    315     left_ps.levels[j-1] = -1
    316     current_ps.levels[j-1] = -1
    317     left_ps.depth = 0
    318     current_ps.depth = 0
    319     left_ps.degree = n
    320     current_ps.degree = n
    321    
    322     # default values of "infinity"
    323     for i from 0 <= i < n:
    324         indicators[i] = -1
     310            mem_err = 1
    325311   
    326312    cdef bint possible = 1
    327313    cdef bint unknown = 1
     314       
     315    if not mem_err:
     316        bitset_zero(vertices_have_been_reduced)
     317       
     318        # set up the identity permutation
     319        for i from 0 <= i < n:
     320            id_perm[i] = i
     321        if ordering2 is NULL:
     322            ordering2 = id_perm
     323       
     324        # Copy reordering of left_ps coming from ordering2 to current_ps.
     325        for i from 0 <= i < n:
     326            current_ps.entries[i] = ordering2[i] # memcpy faster?
     327            current_ps.levels[i] = left_ps.levels[i]
     328        # note that current_ps depth and degree already set
     329       
     330        # default values of "infinity"
     331        for i from 0 <= i < n:
     332            indicators[i] = -1
     333       
     334        # Our first refinement needs to check every cell of the partition,
     335        # so cells_to_refine_by needs to be a list of pointers to each cell.
     336        j = 1
     337        cells_to_refine_by[0] = 0
     338        for i from 0 < i < n:
     339            if left_ps.levels[i-1] == 0:
     340                cells_to_refine_by[j] = i
     341                j += 1
     342        if input_group is NULL:
     343            k = refine_and_return_invariant(left_ps, S1, cells_to_refine_by, j)
     344        else:
     345            k = refine_also_by_orbits(left_ps, S1, refine_and_return_invariant,
     346                cells_to_refine_by, j, group, perm_stack)
     347       
     348        j = 1
     349        cells_to_refine_by[0] = 0
     350        for i from 0 < i < n:
     351            if current_ps.levels[i-1] == 0:
     352                cells_to_refine_by[j] = i
     353                j += 1
     354        if input_group is NULL:
     355            j = refine_and_return_invariant(current_ps, S2, cells_to_refine_by, j)
     356        else:
     357            j = refine_also_by_orbits(current_ps, S2, refine_and_return_invariant,
     358                cells_to_refine_by, j, group, perm_stack)
     359       
     360        if k != j:
     361            possible = 0; unknown = 0
     362        elif not stacks_are_equivalent(left_ps, current_ps):
     363            possible = 0; unknown = 0
     364        else:
     365            PS_move_all_mins_to_front(current_ps)
     366
     367        first_ps = NULL
     368        # Refine down to a discrete partition
     369        while not PS_is_discrete(left_ps):
     370            i = left_ps.depth
     371            k = PS_first_smallest(left_ps, vertices_to_split[i]) # writes to vertices_to_split, but this is never used
     372            if input_group is not NULL:
     373                # split the point
     374                left_ps.depth += 1
     375                PS_clear(left_ps)
     376                cells_to_refine_by[0] = PS_split_point(left_ps, k)
     377                # update the base
     378                tmp_gp = group
     379                group = old_group
     380                old_group = tmp_gp
     381                if SC_insert_base_point_nomalloc(group, old_group, i, k):
     382                    mem_err = 1
     383                    break
     384                SC_identify(perm_stack + n*left_ps.depth, n)
     385                # do the refinements
     386                indicators[i] = refine_also_by_orbits(left_ps, S1, refine_and_return_invariant, cells_to_refine_by, 1, group, perm_stack)
     387            else:
     388                indicators[i] = split_point_and_refine(left_ps, k, S1, refine_and_return_invariant, cells_to_refine_by)
     389            bitset_unset(vertices_have_been_reduced, left_ps.depth)
     390    if not mem_err:
     391        while not PS_is_discrete(current_ps) and possible:
     392            i = current_ps.depth
     393            vertices_determining_current_stack[i] = PS_first_smallest(current_ps, vertices_to_split[i])
     394            if input_group is not NULL:
     395                if group.parents[i][perm_stack[n*i + vertices_determining_current_stack[i]]] == -1:
     396                    possible = 0
     397            while possible:
     398                i = current_ps.depth
     399                if input_group is not NULL:
     400                    # split the point
     401                    current_ps.depth += 1
     402                    PS_clear(current_ps)
     403                    cells_to_refine_by[0] = PS_split_point(current_ps, vertices_determining_current_stack[i])
     404                    # update the base
     405                    tmp_gp = group
     406                    group = old_group
     407                    old_group = tmp_gp
     408                    if SC_insert_base_point_nomalloc(group, old_group, i, vertices_determining_current_stack[i]):
     409                        mem_err = 1
     410                        break
     411                    # update perm_stack
     412                    update_perm_stack(group, current_ps.depth, vertices_determining_current_stack[i], perm_stack)
     413                    # do the refinements
     414                    j = refine_also_by_orbits(current_ps, S2, refine_and_return_invariant, cells_to_refine_by, 1, group, perm_stack)
     415                else:
     416                    j = split_point_and_refine(current_ps,
     417                        vertices_determining_current_stack[i], S2,
     418                        refine_and_return_invariant, cells_to_refine_by)
     419                if indicators[i] != j:
     420                    possible = 0
     421                elif not stacks_are_equivalent(left_ps, current_ps):
     422                    possible = 0
     423                else:
     424                    PS_move_all_mins_to_front(current_ps)
     425                    if not all_children_are_equivalent(current_ps, S2):
     426                        current_kids_are_same = current_ps.depth + 1
     427                    break
     428                current_ps.depth -= 1
     429                while current_ps.depth >= 0: # not possible, so look for another
     430                    i = current_ps.depth
     431                    j = vertices_determining_current_stack[i] + 1
     432                    j = bitset_next(vertices_to_split[i], j)
     433                    if j == -1:
     434                        current_ps.depth -= 1 # backtrack
     435                    else:
     436                        possible = 1
     437                        vertices_determining_current_stack[i] = j
     438                        break # found another
     439        if possible:
     440            if input_group is NULL:
     441                if compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
     442                    unknown = 0
     443            else:
     444                PS_get_perm_from(left_ps, current_ps, permutation)
     445                if SC_contains(group, 0, permutation, 0) and compare_structures(permutation, id_perm, S1, S2) == 0:
     446                    # TODO: might be slight optimization for containment using perm_stack
     447                    unknown = 0
     448            if unknown:
     449                first_meets_current = current_ps.depth
     450                first_kids_are_same = current_ps.depth
     451                first_ps = PS_copy(current_ps)
     452                if first_ps is NULL:
     453                    mem_err = 1
     454                current_ps.depth -= 1
    328455   
    329     # Our first refinement needs to check every cell of the partition,
    330     # so cells_to_refine_by needs to be a list of pointers to each cell.
    331     j = 1
    332     cells_to_refine_by[0] = 0
    333     for i from 0 < i < n:
    334         if left_ps.levels[i-1] == 0:
    335             cells_to_refine_by[j] = i
    336             j += 1
    337    
    338     k = refine_and_return_invariant(left_ps, S1, cells_to_refine_by, j)
    339    
    340     j = 1
    341     cells_to_refine_by[0] = 0
    342     for i from 0 < i < n:
    343         if current_ps.levels[i-1] == 0:
    344             cells_to_refine_by[j] = i
    345             j += 1
    346     j = refine_and_return_invariant(current_ps, S2, cells_to_refine_by, j)
    347     if k != j:
    348         possible = 0; unknown = 0
    349     elif not stacks_are_equivalent(left_ps, current_ps):
    350         possible = 0; unknown = 0
    351     else:
    352         PS_move_all_mins_to_front(current_ps)
     456    if mem_err:
     457        sage_free(int_array)
     458        if input_group is not NULL:
     459            sage_free(perm_stack)
     460            SC_dealloc(group)
     461            SC_dealloc(old_group)
     462        OP_dealloc(orbits_of_subgroup)
     463        sage_free(output)
     464        if bitset_array is not NULL:
     465            for i from 0 <= i < n+2*len_of_fp_and_mcr:
     466                bitset_free(bitset_array[i])
     467        bitset_free(vertices_have_been_reduced)
     468        sage_free(bitset_array)
     469        PS_dealloc(current_ps)
     470        raise MemoryError
    353471
    354     first_ps = NULL
    355     # Refine down to a discrete partition
    356     while not PS_is_discrete(left_ps):
    357         i = left_ps.depth
    358         k = PS_first_smallest(left_ps, vertices_to_split[i]) # writes to vertices_to_split, but this is never used
    359         indicators[i] = split_point_and_refine(left_ps, k, S1, refine_and_return_invariant, cells_to_refine_by)
    360         bitset_unset(vertices_have_been_reduced, left_ps.depth)
    361     while not PS_is_discrete(current_ps) and possible:
    362         i = current_ps.depth
    363         vertices_determining_current_stack[i] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
    364         while possible:
    365             i = current_ps.depth
    366             j = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by)
    367             if indicators[i] != j:
    368                 possible = 0
    369             elif not stacks_are_equivalent(left_ps, current_ps):
    370                 possible = 0
    371             else:
    372                 PS_move_all_mins_to_front(current_ps)
    373                 if not all_children_are_equivalent(current_ps, S2):
    374                     current_kids_are_same = current_ps.depth + 1
    375                 break
    376             current_ps.depth -= 1 # reset from split_point_and_refine
    377             while current_ps.depth >= 0: # not possible, so look for another
    378                 i = current_ps.depth
    379                 j = vertices_determining_current_stack[i] + 1
    380                 j = bitset_next(vertices_to_split[i], j)
    381                 if j == -1:
    382                     current_ps.depth -= 1 # backtrack
    383                 else:
    384                     possible = 1
    385                     vertices_determining_current_stack[i] = j
    386                     break # found another
    387    
    388     if possible:
    389         if compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
    390             unknown = 0
    391         else:
    392             first_meets_current = current_ps.depth
    393             first_kids_are_same = current_ps.depth
    394             first_ps = PS_copy(current_ps)
    395             current_ps.depth -= 1
    396    
    397472    # Main loop:
    398473    while possible and unknown and current_ps.depth != -1:
    399474       
     
    484559        while 1:
    485560            i = current_ps.depth
    486561            while 1:
    487                 k = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by)
     562                if input_group is not NULL:
     563                    k = split_point_and_refine_by_orbits(current_ps,
     564                        vertices_determining_current_stack[i], S2,
     565                        refine_and_return_invariant, cells_to_refine_by,
     566                        group, perm_stack)
     567                    update_perm_stack(group, current_ps.depth, vertices_determining_current_stack[i], perm_stack)
     568                else:
     569                    k = split_point_and_refine(current_ps,
     570                        vertices_determining_current_stack[i], S2,
     571                        refine_and_return_invariant, cells_to_refine_by)
    488572                PS_move_all_mins_to_front(current_ps)
    489573                if indicators[i] != k:
    490574                    possible = 0
    491575                elif not stacks_are_equivalent(left_ps, current_ps):
    492576                    possible = 0
     577                if PS_is_discrete(current_ps):
     578                    break
     579                vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
     580                if input_group is not NULL:
     581                    if group.parents[current_ps.depth][perm_stack[n*current_ps.depth + vertices_determining_current_stack[current_ps.depth]]] == -1:
     582                        possible = 0
    493583                if not possible:
    494584                    j = vertices_determining_current_stack[i] + 1
    495585                    j = bitset_next(vertices_to_split[i], j)
     
    504594                break
    505595            if PS_is_discrete(current_ps):
    506596                break
    507             vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
    508597            bitset_unset(vertices_have_been_reduced, current_ps.depth)
    509598            if not all_children_are_equivalent(current_ps, S2):
    510599                current_kids_are_same = current_ps.depth + 1
     
    512601        # III. Check for automorphisms and isomorphisms
    513602        automorphism = 0
    514603        if possible:
    515             PS_get_perm_from(current_ps, first_ps, permutation)
     604            PS_get_perm_from(first_ps, current_ps, permutation)
    516605            if compare_structures(permutation, id_perm, S2, S2) == 0:
    517                 automorphism = 1
     606                if input_group is NULL or SC_contains(group, 0, permutation, 0):
     607                    # TODO: might be slight optimization for containment using perm_stack
     608                    automorphism = 1
    518609        if not automorphism and possible:
    519610            # if we get here, discrete must be true
    520611            if current_ps.depth != left_ps.depth:
    521612                possible = 0
    522             elif compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
    523                 unknown = 0
    524                 break # main loop
     613            elif input_group is NULL:
     614                if compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
     615                    unknown = 0
     616                    break
     617                else:
     618                    possible = 0
    525619            else:
    526                 possible = 0
     620                PS_get_perm_from(left_ps, current_ps, permutation)
     621                if SC_contains(group, 0, permutation, 0) and compare_structures(permutation, id_perm, S1, S2) == 0:
     622                    # TODO: might be slight optimization for containment using perm_stack
     623                    unknown = 0
     624                    break
     625                else:
     626                    possible = 0
    527627        if automorphism:
    528628            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
    529629                index_in_fp_and_mcr += 1
     
    582682        output = NULL
    583683   
    584684    # Deallocate:
    585     for i from 0 <= i < len_of_fp_and_mcr:
    586         bitset_free(fixed_points_of_generators[i])
    587         bitset_free(minimal_cell_reps_of_generators[i])
    588     for i from 0 <= i < n:
    589         bitset_free(vertices_to_split[i])
     685    sage_free(int_array)
     686    OP_dealloc(orbits_of_subgroup)
     687    if bitset_array is not NULL:
     688        for i from 0 <= i < n+2*len_of_fp_and_mcr:
     689            bitset_free(bitset_array[i])
    590690    bitset_free(vertices_have_been_reduced)
    591     sage_free(indicators)
    592     sage_free(fixed_points_of_generators)
    593     sage_free(minimal_cell_reps_of_generators)
    594     sage_free(vertices_to_split)
    595     sage_free(permutation)
    596     sage_free(id_perm)
    597     sage_free(cells_to_refine_by)
    598     sage_free(vertices_determining_current_stack)
     691    sage_free(bitset_array)
     692    if input_group is not NULL:
     693        sage_free(perm_stack)
     694        SC_dealloc(group)
     695        SC_dealloc(old_group)
     696    PS_dealloc(first_ps)
    599697    PS_dealloc(current_ps)
    600     PS_dealloc(left_ps)
    601     OP_dealloc(orbits_of_subgroup)
    602     if first_ps is not NULL: PS_dealloc(first_ps)
    603    
    604698    return output
    605699
    606700
  • sage/groups/perm_gps/partn_ref/refinement_binary.pxd

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_binary.pxd
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#
    55# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    66#                         http://www.gnu.org/licenses/
     
    1111include 'data_structures_pxd.pxi' # includes bitsets
    1212
    1313from sage.rings.integer cimport Integer
    14 from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     14from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab
    1515from double_coset cimport double_coset
    1616
    1717cdef class BinaryCodeStruct:
     
    2222    cdef PartitionStack *word_ps
    2323    cdef int *alpha # length nwords + degree
    2424    cdef int *scratch # length 3*nwords + 3*degree + 2
    25     cdef aut_gp_and_can_lab_return *output
     25    cdef aut_gp_and_can_lab *output
    2626    cdef int (*ith_word)(BinaryCodeStruct self, int, bitset_s *)
    2727
    2828cdef class LinearBinaryCodeStruct(BinaryCodeStruct):
     
    3636    cdef bitset_s *scratch_bitsets # length 4*nwords + 1
    3737cdef int ith_word_nonlinear(BinaryCodeStruct, int, bitset_s *)
    3838
    39 cdef int refine_by_bip_degree(PartitionStack *, object, int *, int)
    40 cdef int compare_linear_codes(int *, int *, object, object)
    41 cdef int compare_nonlinear_codes(int *, int *, object, object)
    42 cdef bint all_children_are_equivalent(PartitionStack *, object)
     39cdef int refine_by_bip_degree(PartitionStack *, void *, int *, int)
     40cdef int compare_linear_codes(int *, int *, void *, void *)
     41cdef int compare_nonlinear_codes(int *, int *, void *, void *)
     42cdef bint all_children_are_equivalent(PartitionStack *, void *)
    4343cdef inline int word_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *)
    4444cdef inline int col_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *)
    45 cdef inline int sort_by_function(PartitionStack *, int, int *, int *, int *, int)
     45cdef inline int sort_by_function_codes(PartitionStack *, int, int *, int *, int *, int)
    4646
    4747
    4848
  • sage/groups/perm_gps/partn_ref/refinement_binary.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_binary.pyx
    a b  
    1616"""
    1717
    1818#*****************************************************************************
    19 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     19#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    2020#
    2121# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    2222#                         http://www.gnu.org/licenses/
     
    222222            17868913969917295853568000000
    223223 
    224224        """
    225         cdef int **part, i, j
     225        cdef int i, n = self.degree
     226        cdef PartitionStack *part
    226227        if partition is None:
    227             partition = [range(self.degree)]
    228         part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     228            part = PS_new(n, 1)
     229        else:
     230            part = PS_from_list(partition)
    229231        if part is NULL:
    230232            raise MemoryError
    231         for i from 0 <= i < len(partition):
    232             part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    233             if part[i] is NULL:
    234                 for j from 0 <= j < i:
    235                     sage_free(part[j])
    236                 sage_free(part)
    237                 raise MemoryError
    238             for j from 0 <= j < len(partition[i]):
    239                 part[i][j] = partition[i][j]
    240             part[i][len(partition[i])] = -1
    241         part[len(partition)] = NULL
     233
    242234        self.first_time = 1
    243235       
    244         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)
     236        self.output = get_aut_gp_and_can_lab(<void *>self, part, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes, 1, NULL)
    245237       
    246         for i from 0 <= i < len(partition):
    247             sage_free(part[i])
    248         sage_free(part)
     238        PS_dealloc(part)
    249239   
    250240    def automorphism_group(self):
    251241        """
     
    262252
    263253        """
    264254        cdef int i, j
    265         cdef object generators, base
     255        cdef list generators, base
    266256        cdef Integer order
    267257        if self.output is NULL:
    268258            self.run()
     
    270260        for i from 0 <= i < self.output.num_gens:
    271261            generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree])
    272262        order = Integer()
    273         mpz_set(order.value, self.output.order)
    274         base = [self.output.base[i] for i from 0 <= i < self.output.base_size]
     263        SC_order(self.output.group, 0, order.value)
     264        base = [self.output.group.base_orbits[i][0] for i from 0 <= i < self.output.group.base_size]
    275265        return generators, order, base
    276266   
    277267    def canonical_relabeling(self):
     
    311301            [0, 1, 2, 5, 3, 4]
    312302           
    313303        """
    314         cdef int **part, i, j
     304        cdef int i, n = self.degree
    315305        cdef int *output, *ordering
    316         partition = [range(self.degree)]
    317         part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     306        cdef PartitionStack *part
     307        part = PS_new(n, 1)
    318308        ordering = <int *> sage_malloc(self.degree * sizeof(int))
    319309        if part is NULL or ordering is NULL:
    320             if part is not NULL: sage_free(part)
     310            if part is not NULL: PS_dealloc(part)
    321311            if ordering is not NULL: sage_free(ordering)
    322312            raise MemoryError
    323         for i from 0 <= i < len(partition):
    324             part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    325             if part[i] is NULL:
    326                 for j from 0 <= j < i:
    327                     sage_free(part[j])
    328                 sage_free(part)
    329                 raise MemoryError
    330             for j from 0 <= j < len(partition[i]):
    331                 part[i][j] = partition[i][j]
    332             part[i][len(partition[i])] = -1
    333         part[len(partition)] = NULL
    334         for i from 0 <= i < self.degree:
     313        for i from 0 <= i < n:
    335314            ordering[i] = i
    336315        self.first_time = 1
    337316        other.first_time = 1
    338317       
    339         output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes)
     318        output = double_coset(<void *> self, <void *> other, part, ordering, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes, NULL)
    340319       
    341         for i from 0 <= i < len(partition):
    342             sage_free(part[i])
    343         sage_free(part)
     320        PS_dealloc(part)
    344321        sage_free(ordering)
    345322
    346323        if output is NULL:
    347324            return False
    348325        else:
    349             output_py = [output[i] for i from 0 <= i < self.degree]
     326            output_py = [output[i] for i from 0 <= i < n]
    350327            sage_free(output)
    351328            return output_py
    352329   
     
    361338        sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
    362339        sage_free(self.alpha); sage_free(self.scratch)
    363340        if self.output is not NULL:
    364             mpz_clear(self.output.order)
    365341            sage_free(self.output.generators)
    366             sage_free(self.output.base)
     342            SC_dealloc(self.output.group)
    367343            sage_free(self.output.relabeling)
    368344            sage_free(self.output)
    369345           
     
    468444        sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
    469445        sage_free(self.alpha); sage_free(self.scratch)
    470446        if self.output is not NULL:
    471             mpz_clear(self.output.order)
    472447            sage_free(self.output.generators)
    473             sage_free(self.output.base)
     448            SC_dealloc(self.output.group)
    474449            sage_free(self.output.relabeling)
    475450            sage_free(self.output)
    476451           
     
    513488            [2, 3, 4, 5, 0, 1]
    514489
    515490        """
    516         cdef int **part, i, j
     491        cdef int n = self.degree
     492        cdef PartitionStack *part
    517493        if partition is None:
    518             partition = [range(self.degree)]
    519         part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     494            part = PS_new(n, 1)
     495        else:
     496            part = PS_from_list(partition)
    520497        if part is NULL:
    521498            raise MemoryError
    522         for i from 0 <= i < len(partition):
    523             part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    524             if part[i] is NULL:
    525                 for j from 0 <= j < i:
    526                     sage_free(part[j])
    527                 sage_free(part)
    528                 raise MemoryError
    529             for j from 0 <= j < len(partition[i]):
    530                 part[i][j] = partition[i][j]
    531             part[i][len(partition[i])] = -1
    532         part[len(partition)] = NULL
    533499        self.first_time = 1
    534500       
    535         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)
     501        self.output = get_aut_gp_and_can_lab(<void *> self, part, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes, 1, NULL)
    536502       
    537         for i from 0 <= i < len(partition):
    538             sage_free(part[i])
    539         sage_free(part)
     503        PS_dealloc(part)
    540504   
    541505    def automorphism_group(self):
    542506        """
     
    559523
    560524        """
    561525        cdef int i, j
    562         cdef object generators, base
     526        cdef list generators, base
    563527        cdef Integer order
    564528        if self.output is NULL:
    565529            self.run()
     
    567531        for i from 0 <= i < self.output.num_gens:
    568532            generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree])
    569533        order = Integer()
    570         mpz_set(order.value, self.output.order)
    571         base = [self.output.base[i] for i from 0 <= i < self.output.base_size]
     534        SC_order(self.output.group, 0, order.value)
     535        base = [self.output.group.base_orbits[i][0] for i from 0 <= i < self.output.group.base_size]
    572536        return generators, order, base
    573537   
    574538    def canonical_relabeling(self):
     
    602566            [2, 3, 0, 1, 4, 5]
    603567           
    604568        """
    605         cdef int **part, i, j
     569        cdef int i, n = self.degree
    606570        cdef int *output, *ordering
    607         partition = [range(self.degree)]
    608         part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
    609         ordering = <int *> sage_malloc(self.degree * sizeof(int))
     571        cdef PartitionStack *part
     572        part = PS_new(n, 1)
     573        ordering = <int *> sage_malloc(n * sizeof(int))
    610574        if part is NULL or ordering is NULL:
    611             if part is not NULL: sage_free(part)
     575            if part is not NULL: PS_dealloc(part)
    612576            if ordering is not NULL: sage_free(ordering)
    613577            raise MemoryError
    614         for i from 0 <= i < len(partition):
    615             part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    616             if part[i] is NULL:
    617                 for j from 0 <= j < i:
    618                     sage_free(part[j])
    619                 sage_free(part)
    620                 raise MemoryError
    621             for j from 0 <= j < len(partition[i]):
    622                 part[i][j] = partition[i][j]
    623             part[i][len(partition[i])] = -1
    624         part[len(partition)] = NULL
    625         for i from 0 <= i < self.degree:
     578        for i from 0 <= i < n:
    626579            ordering[i] = i
    627580        self.first_time = 1
    628581        other.first_time = 1
    629582       
    630         output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes)
     583        output = double_coset(<void *> self, <void *> other, part, ordering, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes, NULL)
    631584       
    632         for i from 0 <= i < len(partition):
    633             sage_free(part[i])
    634         sage_free(part)
     585        PS_dealloc(part)
    635586        sage_free(ordering)
    636587
    637588        if output is NULL:
     
    646597    bitset_copy(word, &NBCS.words[i])
    647598    return 0
    648599
    649 cdef int refine_by_bip_degree(PartitionStack *col_ps, object S, int *cells_to_refine_by, int ctrb_len):
     600cdef int refine_by_bip_degree(PartitionStack *col_ps, void *S, int *cells_to_refine_by, int ctrb_len):
    650601    """
    651602    Refines the input partition by checking degrees of vertices to the given
    652603    cells in the associated bipartite graph (vertices split into columns and
     
    711662                # now, i points to the next cell (before refinement)
    712663                if necessary_to_split_cell:
    713664                    invariant += 8
    714                     first_largest_subcell = sort_by_function(col_ps, current_cell, col_degrees, col_counts, col_output, BCS.nwords+1)
     665                    first_largest_subcell = sort_by_function_codes(col_ps, current_cell, col_degrees, col_counts, col_output, BCS.nwords+1)
    715666                    invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps)
    716667                    invariant += first_largest_subcell
    717668                    against_index = current_cell_against
     
    746697                # now, i points to the next cell (before refinement)
    747698                if necessary_to_split_cell:
    748699                    invariant += 64
    749                     first_largest_subcell = sort_by_function(word_ps, current_cell, word_degrees, word_counts, word_output, BCS.degree+1)
     700                    first_largest_subcell = sort_by_function_codes(word_ps, current_cell, word_degrees, word_counts, word_output, BCS.degree+1)
    750701                    invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps)
    751702                    invariant += first_largest_subcell
    752703                    against_index = current_cell_against
     
    770721        current_cell_against += 1
    771722    return invariant
    772723
    773 cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S1, object S2):
     724cdef int compare_linear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2):
    774725    """
    775726    Compare gamma_1(S1) and gamma_2(S2).
    776727   
     
    843794                    return <int>bitset_check(&basis_2[i], gamma_2[cur_col]) - <int>bitset_check(&basis_1[i], gamma_1[cur_col])
    844795    return 0
    845796
    846 cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S1, object S2):
     797cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2):
    847798    """
    848799    Compare gamma_1(S1) and gamma_2(S2).
    849800   
     
    923874
    924875    return 0
    925876
    926 cdef bint all_children_are_equivalent(PartitionStack *col_ps, object S):
     877cdef bint all_children_are_equivalent(PartitionStack *col_ps, void *S):
    927878    """
    928879    Returns True if any refinement of the current partition results in the same
    929880    structure.
     
    1022973    bitset_free(word)
    1023974    return degree
    1024975
    1025 cdef inline int sort_by_function(PartitionStack *PS, int start, int *degrees, int *counts, int *output, int count_max):
     976cdef inline int sort_by_function_codes(PartitionStack *PS, int start, int *degrees, int *counts, int *output, int count_max):
    1026977    """
    1027978    A simple counting sort, given the degrees of vertices to a certain cell.
    1028979   
  • sage/groups/perm_gps/partn_ref/refinement_graphs.pxd

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_graphs.pxd
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#
    55# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    66#                         http://www.gnu.org/licenses/
     
    1414from sage.graphs.base.sparse_graph cimport SparseGraph
    1515from sage.graphs.base.dense_graph cimport DenseGraph
    1616from sage.rings.integer cimport Integer
    17 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     17from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab
    1818from double_coset cimport double_coset
    1919
    2020cdef class GraphStruct:
     
    2323    cdef bint use_indicator
    2424    cdef int *scratch # length 3n+1
    2525
    26 cdef int refine_by_degree(PartitionStack *, object, int *, int)
    27 cdef int compare_graphs(int *, int *, object, object)
    28 cdef bint all_children_are_equivalent(PartitionStack *, object)
     26cdef int refine_by_degree(PartitionStack *, void *, int *, int)
     27cdef int compare_graphs(int *, int *, void *, void *)
     28cdef bint all_children_are_equivalent(PartitionStack *, void *)
    2929cdef inline int degree(PartitionStack *, CGraph, int, int, bint)
    30 cdef inline int sort_by_function(PartitionStack *, int, int *)
    3130
    3231
    3332
  • sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_graphs.pyx
    a b  
    1212"""
    1313
    1414#*****************************************************************************
    15 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     15#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    1616#
    1717# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    1818#                         http://www.gnu.org/licenses/
     
    4949    {0: 1, 1: 2, 2: 0}
    5050
    5151    """
    52     cdef int **part
     52    cdef PartitionStack *part
    5353    cdef int *output, *ordering
    5454    cdef CGraph G
    5555    cdef GraphStruct GS1 = GraphStruct()
    5656    cdef GraphStruct GS2 = GraphStruct()
    5757    cdef GraphStruct GS
    58     cdef int i, j, n = -1
     58    cdef int i, j, k, n = -1, cell_len
     59    cdef list partition, cell
    5960   
    6061    from sage.graphs.all import Graph, DiGraph
    6162    from sage.graphs.generic_graph import GenericGraph
     
    8283                if first:
    8384                    partition = [[to[v] for v in cell] for cell in partn]
    8485            else:
     86                if first:
     87                    partition = partn
    8588                to = range(n)
    8689                frm = to
    8790            if sparse:
     
    104107            to = {}
    105108            for a in G.verts(): to[a]=a
    106109            frm = to
     110            if first:
     111                partition = partn
    107112        else:
    108113            raise TypeError("G must be a Sage graph.")
    109114        if first: frm1=frm;to1=to
     
    115120
    116121    if n == 0:
    117122        return {}
    118    
    119     part = <int **> sage_malloc((len(partn)+1) * sizeof(int *))
     123
     124    part = PS_from_list(partition)
    120125    ordering = <int *> sage_malloc(n * sizeof(int))
    121126    if part is NULL or ordering is NULL:
    122         if part is not NULL: sage_free(part)
     127        if part is not NULL: PS_dealloc(part)
    123128        if ordering is not NULL: sage_free(ordering)
    124129        raise MemoryError
    125     for i from 0 <= i < len(partn):
    126         part[i] = <int *> sage_malloc((len(partn[i])+1) * sizeof(int))
    127         if part[i] is NULL:
    128             for j from 0 <= j < i:
    129                 sage_free(part[j])
    130             sage_free(part)
    131             raise MemoryError
    132         for j from 0 <= j < len(partn[i]):
    133             part[i][j] = to1[partn[i][j]]
    134         part[i][len(partn[i])] = -1
    135     part[len(partn)] = NULL
    136130    for i from 0 <= i < n:
    137131        ordering[i] = to2[ordering2[i]]
    138132
     
    141135    if GS1.scratch is NULL or GS2.scratch is NULL:
    142136        if GS1.scratch is not NULL: sage_free(GS1.scratch)
    143137        if GS2.scratch is not NULL: sage_free(GS2.scratch)
    144         for j from 0 <= j < len(partition):
    145             sage_free(part[j])
    146         sage_free(part)
     138        PS_dealloc(part)
     139        sage_free(ordering)
    147140        raise MemoryError
    148141
    149     output = double_coset(GS1, GS2, part, ordering, n, &all_children_are_equivalent, &refine_by_degree, &compare_graphs)
     142    output = double_coset(<void *>GS1, <void *>GS2, part, ordering, n, &all_children_are_equivalent, &refine_by_degree, &compare_graphs, NULL)
    150143
    151     for i from 0 <= i < len(partn):
    152         sage_free(part[i])
    153     sage_free(part)
     144    PS_dealloc(part)
    154145    sage_free(ordering)
    155146    sage_free(GS1.scratch)
    156147    sage_free(GS2.scratch)
     
    350341    cdef CGraph G
    351342    cdef int i, j, n
    352343    cdef Integer I
    353     cdef aut_gp_and_can_lab_return *output
    354     cdef int **part
     344    cdef aut_gp_and_can_lab *output
     345    cdef PartitionStack *part
    355346    from sage.graphs.all import Graph, DiGraph
    356347    from sage.graphs.generic_graph import GenericGraph
    357348    from copy import copy
     
    387378    else:
    388379        raise TypeError("G must be a Sage graph.")
    389380   
    390     part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
    391     if part is NULL:
    392         raise MemoryError
    393     for i from 0 <= i < len(partition):
    394         part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    395         if part[i] is NULL:
    396             for j from 0 <= j < i:
    397                 sage_free(part[j])
    398             sage_free(part)
    399             raise MemoryError
    400         for j from 0 <= j < len(partition[i]):
    401             part[i][j] = partition[i][j]
    402         part[i][len(partition[i])] = -1
    403     part[len(partition)] = NULL
    404 
    405381    cdef GraphStruct GS = GraphStruct()
    406382    GS.G = G
    407383    GS.directed = 1 if dig else 0
    408384    GS.use_indicator = 1 if use_indicator_function else 0
     385
     386    if n == 0:
     387        return_tuple = [[]]
     388        if dict_rep:
     389            return_tuple.append({})
     390        if lab:
     391            if isinstance(G_in, GenericGraph):
     392                G_C = copy(G_in)
     393            else:
     394                if isinstance(G, SparseGraph):
     395                    G_C = SparseGraph(n)
     396                else:
     397                    G_C = DenseGraph(n)
     398            return_tuple.append(G_C)
     399        if certify:
     400            return_tuple.append({})
     401        if base:
     402            return_tuple.append([])
     403        if order:
     404            return_tuple.append(Integer(1))
     405        if len(return_tuple) == 1:
     406            return return_tuple[0]
     407        else:
     408            return tuple(return_tuple)
     409
    409410    GS.scratch = <int *> sage_malloc( (3*G.num_verts + 1) * sizeof(int) )
    410     if GS.scratch is NULL:
    411         for j from 0 <= j < len(partition):
    412             sage_free(part[j])
    413         sage_free(part)
     411    part = PS_from_list(partition)
     412    if GS.scratch is NULL or part is NULL:
     413        if part is not NULL: PS_dealloc(part)
     414        if GS.scratch is not NULL: sage_free(GS.scratch)
    414415        raise MemoryError
    415416
    416417    lab_new = lab or certify
    417     output = get_aut_gp_and_can_lab(GS, part, G.num_verts, &all_children_are_equivalent, &refine_by_degree, &compare_graphs, lab_new, base, order)
     418    output = get_aut_gp_and_can_lab(<void *>GS, part, G.num_verts, &all_children_are_equivalent, &refine_by_degree, &compare_graphs, lab, NULL)
    418419    sage_free( GS.scratch )
    419420    # prepare output
    420421    list_of_gens = []
     
    445446            dd[frm[i]] = output.relabeling[i]
    446447        return_tuple.append(dd)
    447448    if base:
    448         return_tuple.append([output.base[i] for i from 0 <= i < output.base_size])
     449        return_tuple.append([output.group.base_orbits[i][0] for i from 0 <= i < output.group.base_size])
    449450    if order:
    450451        I = Integer()
    451         mpz_set(I.value, output.order)
     452        SC_order(output.group, 0, I.value)
    452453        return_tuple.append(I)
    453     for i from 0 <= i < len(partition):
    454         sage_free(part[i])
    455     sage_free(part)
    456     mpz_clear(output.order)
     454    PS_dealloc(part)
    457455    sage_free(output.generators)
    458     if base:
    459         sage_free(output.base)
     456    SC_dealloc(output.group)
    460457    if lab_new:
    461458        sage_free(output.relabeling)
    462459    sage_free(output)
     
    465462    else:
    466463        return tuple(return_tuple)
    467464
    468 cdef int refine_by_degree(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     465cdef int refine_by_degree(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len):
    469466    r"""
    470467    Refines the input partition by checking degrees of vertices to the given
    471468    cells.
     
    584581    else:
    585582        return 0
    586583   
    587 cdef int compare_graphs(int *gamma_1, int *gamma_2, object S1, object S2):
     584cdef int compare_graphs(int *gamma_1, int *gamma_2, void *S1, void *S2):
    588585    r"""
    589586    Compare gamma_1(S1) and gamma_2(S2).
    590587
     
    611608                return -1
    612609    return 0
    613610
    614 cdef bint all_children_are_equivalent(PartitionStack *PS, object S):
     611cdef bint all_children_are_equivalent(PartitionStack *PS, void *S):
    615612    """
    616613    Return True if every refinement of the current partition results in the
    617614    same structure.
     
    681678                break
    682679    return num_arcs
    683680
    684 cdef inline int sort_by_function(PartitionStack *PS, int start, int *degrees):
    685     """
    686     A simple counting sort, given the degrees of vertices to a certain cell.
    687    
    688     INPUT:
    689     PS -- the partition stack to be checked
    690     start -- beginning index of the cell to be sorted
    691     degrees -- the values to be sorted by
    692 
    693     """
    694     cdef int n = PS.degree
    695     cdef int i, j, max, max_location
    696     cdef int *counts = degrees + n, *output = degrees + 2*n + 1
    697     for i from 0 <= i <= n:
    698         counts[i] = 0
    699     i = 0
    700     while PS.levels[i+start] > PS.depth:
    701         counts[degrees[i]] += 1
    702         i += 1
    703     counts[degrees[i]] += 1
    704     # i+start is the right endpoint of the cell now
    705     max = counts[0]
    706     max_location = 0
    707     for j from 0 < j <= n:
    708         if counts[j] > max:
    709             max = counts[j]
    710             max_location = j
    711         counts[j] += counts[j - 1]
    712     for j from i >= j >= 0:
    713         counts[degrees[j]] -= 1
    714         output[counts[degrees[j]]] = PS.entries[start+j]
    715     max_location = counts[max_location]+start
    716     for j from 0 <= j <= i:
    717         PS.entries[start+j] = output[j]
    718     j = 1
    719     while j <= n and counts[j] <= i:
    720         if counts[j] > 0:
    721             PS.levels[start + counts[j] - 1] = PS.depth
    722         PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1)
    723         j += 1
    724     return max_location
    725 
    726681def all_labeled_graphs(n):
    727682    """
    728683    Returns all labeled graphs on n vertices {0,1,...,n-1}. Used in
     
    773728    return Glist
    774729
    775730
    776 def random_tests(num=20, n_max=50, perms_per_graph=8):
     731def random_tests(num=10, n_max=60, perms_per_graph=5):
    777732    """
    778733    Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma
    779734    and random graphs G, and that isomorphic returns an isomorphism.
     
    796751    TESTS::
    797752
    798753        sage: import sage.groups.perm_gps.partn_ref.refinement_graphs
    799         sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests()  # long time (up to 25s on sage.math, 2011)
    800         All passed: 640 random tests on 40 graphs.
     754        sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests()  # long time
     755        All passed: 200 random tests on 20 graphs.
    801756    """
    802757    from sage.misc.prandom import random, randint
    803758    from sage.graphs.graph_generators import GraphGenerators
     
    982937    GS.G = G
    983938    GS.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
    984939    if GS.scratch is NULL:
    985         PS_clear(nu)
     940        PS_dealloc(nu)
    986941        raise MemoryError
    987942    GS.directed = directed
    988943    GS.use_indicator = 0
     
    991946    cdef int num_cells = len(partition)
    992947    cdef int *alpha = <int *>sage_malloc(n * sizeof(int))
    993948    if alpha is NULL:
    994         PS_clear(nu)
     949        PS_dealloc(nu)
    995950        sage_free(GS.scratch)
    996951        raise MemoryError
    997952    j = 0
     
    1000955        j += len(partition[i])
    1001956   
    1002957    # refine, and get the result
    1003     refine_by_degree(nu, GS, alpha, num_cells)
     958    refine_by_degree(nu, <void *>GS, alpha, num_cells)
    1004959   
    1005960    eq_part = []
    1006961    cell = []
     
    1010965            eq_part.append(cell)
    1011966            cell = []
    1012967   
    1013     PS_clear(nu)
     968    PS_dealloc(nu)
    1014969    sage_free(GS.scratch)
    1015970    sage_free(alpha)
    1016971
  • sage/groups/perm_gps/partn_ref/refinement_lists.pxd

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_lists.pxd
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#      Copyright (C) 2009 Nicolas Borie <nicolas.borie@math.u-psud.fr>
    55#
    66# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     
    1717from double_coset cimport double_coset
    1818
    1919# name of the three functions to customize
    20 cdef int refine_list(PartitionStack *, object, int *, int)
    21 cdef int compare_lists(int *, int *, object, object)
    22 cdef bint all_list_children_are_equivalent(PartitionStack *, object)
    23  No newline at end of file
     20cdef int refine_list(PartitionStack *, void *, int *, int)
     21cdef int compare_lists(int *, int *, void *, void *)
     22cdef bint all_list_children_are_equivalent(PartitionStack *, void *)
     23 No newline at end of file
  • sage/groups/perm_gps/partn_ref/refinement_lists.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_lists.pyx
    a b  
    77"""
    88
    99#*****************************************************************************
    10 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     10#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    1111#      Copyright (C) 2009 Nicolas Borie <nicolas.borie@math.u-psud.fr>
    1212#
    1313# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     
    2828        [1, 2, 0]
    2929
    3030    """
    31     cdef int **part, i, j
     31    cdef int i, n = len(self)
     32    cdef PartitionStack *part
    3233    cdef int *output, *ordering
    33     partition = [range(len(self))]
    34     part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     34    part = PS_new(n, 1)
    3535    ordering = <int *> sage_malloc((len(self)) * sizeof(int))
    3636    if part is NULL or ordering is NULL:
    37         if part is not NULL: sage_free(part)
     37        if part is not NULL: PS_dealloc(part)
    3838        if ordering is not NULL: sage_free(ordering)
    3939        raise MemoryError
    40     for i from 0 <= i < len(partition):
    41         part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    42         if part[i] is NULL:
    43             for j from 0 <= j < i:
    44                 sage_free(part[j])
    45             sage_free(part)
    46             raise MemoryError
    47         for j from 0 <= j < len(partition[i]):
    48             part[i][j] = partition[i][j]
    49         part[i][len(partition[i])] = -1
    50     part[len(partition)] = NULL
    5140    for i from 0 <= i < (len(self)):
    5241        ordering[i] = i
    5342       
    54     output = double_coset(self, other, part, ordering, (len(self)), &all_list_children_are_equivalent, &refine_list, &compare_lists)
     43    output = double_coset(<void *> self, <void *> other, part, ordering, (len(self)), &all_list_children_are_equivalent, &refine_list, &compare_lists, NULL)
    5544       
    56     for i from 0 <= i < len(partition):
    57         sage_free(part[i])
    58     sage_free(part)
     45    PS_dealloc(part)
    5946    sage_free(ordering)
    6047
    6148    if output is NULL:
     
    6552        sage_free(output)
    6653        return output_py
    6754
    68 cdef bint all_list_children_are_equivalent(PartitionStack *PS, object S):
     55cdef bint all_list_children_are_equivalent(PartitionStack *PS, void *S):
    6956    return 0
    7057
    71 cdef int refine_list(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     58cdef int refine_list(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len):
    7259    return 0
    7360
    74 cdef int compare_lists(int *gamma_1, int *gamma_2, object S1, object S2):
     61cdef int compare_lists(int *gamma_1, int *gamma_2, void *S1, void *S2):
    7562    r"""
    7663    Compare two lists according to the lexicographic order.
    7764    """
  • sage/groups/perm_gps/partn_ref/refinement_matrices.pxd

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_matrices.pxd
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#
    55# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    66#                         http://www.gnu.org/licenses/
     
    1111include 'data_structures_pxd.pxi' # includes bitsets
    1212
    1313from sage.rings.integer cimport Integer
    14 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     14from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab
    1515from refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree
    1616from refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent
    1717from double_coset cimport double_coset
     
    2424    cdef list symbols
    2525    cdef int nsymbols
    2626    cdef PartitionStack *temp_col_ps
    27     cdef aut_gp_and_can_lab_return *output
     27    cdef aut_gp_and_can_lab *output
    2828
    29 cdef int refine_matrix(PartitionStack *, object, int *, int)
    30 cdef int compare_matrices(int *, int *, object, object)
    31 cdef bint all_matrix_children_are_equivalent(PartitionStack *, object)
     29cdef int refine_matrix(PartitionStack *, void *, int *, int)
     30cdef int compare_matrices(int *, int *, void *, void *)
     31cdef bint all_matrix_children_are_equivalent(PartitionStack *, void *)
    3232
  • sage/groups/perm_gps/partn_ref/refinement_matrices.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_matrices.pyx
    a b  
    1616"""
    1717
    1818#*****************************************************************************
    19 #      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     19#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    2020#
    2121# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    2222#                         http://www.gnu.org/licenses/
     
    8181    def __dealloc__(self):
    8282        PS_dealloc(self.temp_col_ps)
    8383        if self.output is not NULL:
    84             mpz_clear(self.output.order)
    8584            sage_free(self.output.generators)
    86             sage_free(self.output.base)
     85            SC_dealloc(self.output.group)
    8786            sage_free(self.output.relabeling)
    8887            sage_free(self.output)
    8988
     
    148147            True
    149148
    150149        """
    151         cdef int **part, i, j
     150        cdef int i, n = self.degree
     151        cdef PartitionStack *part
    152152        cdef NonlinearBinaryCodeStruct S_temp
    153153        for i from 0 <= i < self.nsymbols:
    154154            S_temp = <NonlinearBinaryCodeStruct> self.symbol_structs[i]
    155155            S_temp.first_time = 1
    156156
    157157        if partition is None:
    158             partition = [range(self.degree)]
    159         part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     158            part = PS_new(n, 1)
     159        else:
     160            part = PS_from_list(partition)
    160161        if part is NULL:
    161162            raise MemoryError
    162         for i from 0 <= i < len(partition):
    163             part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    164             if part[i] is NULL:
    165                 for j from 0 <= j < i:
    166                     sage_free(part[j])
    167                 sage_free(part)
    168                 raise MemoryError
    169             for j from 0 <= j < len(partition[i]):
    170                 part[i][j] = partition[i][j]
    171             part[i][len(partition[i])] = -1
    172         part[len(partition)] = NULL
    173163       
    174         self.output = get_aut_gp_and_can_lab(self, part, self.degree, &all_matrix_children_are_equivalent, &refine_matrix, &compare_matrices, 1, 1, 1)
     164        self.output = get_aut_gp_and_can_lab(<void *> self, part, self.degree, &all_matrix_children_are_equivalent, &refine_matrix, &compare_matrices, 1, NULL)
    175165       
    176         for i from 0 <= i < len(partition):
    177             sage_free(part[i])
    178         sage_free(part)
     166        PS_dealloc(part)
    179167   
    180168
    181169    def automorphism_group(self):
     
    193181
    194182        """
    195183        cdef int i, j
    196         cdef object generators, base
     184        cdef list generators, base
    197185        cdef Integer order
    198186        if self.output is NULL:
    199187            self.run()
     
    201189        for i from 0 <= i < self.output.num_gens:
    202190            generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree])
    203191        order = Integer()
    204         mpz_set(order.value, self.output.order)
    205         base = [self.output.base[i] for i from 0 <= i < self.output.base_size]
     192        SC_order(self.output.group, 0, order.value)
     193        base = [self.output.group.base_orbits[i][0] for i from 0 <= i < self.output.group.base_size]
    206194        return generators, order, base
    207195   
    208196    def canonical_relabeling(self):
     
    234222            [0, 2, 4, 1, 3, 5]
    235223           
    236224        """
    237        
    238         cdef int **part, i, j
     225        cdef int i, j, n = self.degree
    239226        cdef int *output, *ordering
     227        cdef PartitionStack *part
    240228        cdef NonlinearBinaryCodeStruct S_temp
    241229        for i from 0 <= i < self.nsymbols:
    242230            S_temp = self.symbol_structs[i]
    243231            S_temp.first_time = 1
    244232            S_temp = other.symbol_structs[i]
    245233            S_temp.first_time = 1
    246         partition = [range(self.degree)]
    247         part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     234        part = PS_new(n, 1)
    248235        ordering = <int *> sage_malloc(self.degree * sizeof(int))
    249236        if part is NULL or ordering is NULL:
    250             if part is not NULL: sage_free(part)
     237            if part is not NULL: PS_dealloc(part)
    251238            if ordering is not NULL: sage_free(ordering)
    252239            raise MemoryError
    253         for i from 0 <= i < len(partition):
    254             part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    255             if part[i] is NULL:
    256                 for j from 0 <= j < i:
    257                     sage_free(part[j])
    258                 sage_free(part)
    259                 raise MemoryError
    260             for j from 0 <= j < len(partition[i]):
    261                 part[i][j] = partition[i][j]
    262             part[i][len(partition[i])] = -1
    263         part[len(partition)] = NULL
    264240        for i from 0 <= i < self.degree:
    265241            ordering[i] = i
    266242       
    267         output = double_coset(self, other, part, ordering, self.degree, &all_matrix_children_are_equivalent, &refine_matrix, &compare_matrices)
     243        output = double_coset(<void *> self, <void *> other, part, ordering, self.degree, &all_matrix_children_are_equivalent, &refine_matrix, &compare_matrices, NULL)
    268244       
    269         for i from 0 <= i < len(partition):
    270             sage_free(part[i])
    271         sage_free(part)
     245        PS_dealloc(part)
    272246        sage_free(ordering)
    273247
    274248        if output is NULL:
     
    278252            sage_free(output)
    279253            return output_py
    280254   
    281 cdef int refine_matrix(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     255cdef int refine_matrix(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len):
    282256    cdef MatrixStruct M = <MatrixStruct> S
    283257    cdef int i, temp_inv, invariant = 1
    284258    cdef bint changed = 1
    285259    while changed:
    286260        PS_copy_from_to(PS, M.temp_col_ps)
    287261        for BCS in M.symbol_structs:
    288             temp_inv = refine_by_bip_degree(PS, BCS, cells_to_refine_by, ctrb_len)
     262            temp_inv = refine_by_bip_degree(PS, <void *> BCS, cells_to_refine_by, ctrb_len)
    289263            invariant *= temp_inv + 1
    290         if (memcmp(PS.entries, M.temp_col_ps.entries, M.degree * sizeof(int)) == 0
    291          and memcmp(PS.levels, M.temp_col_ps.levels, M.degree * sizeof(int)) == 0):
     264        if memcmp(PS.entries, M.temp_col_ps.entries, 2*M.degree * sizeof(int)) == 0:
    292265            changed = 0
    293266    return invariant
    294267   
    295 cdef int compare_matrices(int *gamma_1, int *gamma_2, object S1, object S2):
     268cdef int compare_matrices(int *gamma_1, int *gamma_2, void *S1, void *S2):
    296269    cdef MatrixStruct MS1 = <MatrixStruct> S1
    297270    cdef MatrixStruct MS2 = <MatrixStruct> S2
    298271    M1 = MS1.matrix
     
    305278        MM2.set_column(i, M2.column(gamma_2[i]))
    306279    return cmp(sorted(MM1.rows()), sorted(MM2.rows()))
    307280
    308 cdef bint all_matrix_children_are_equivalent(PartitionStack *PS, object S):
     281cdef bint all_matrix_children_are_equivalent(PartitionStack *PS, void *S):
    309282    return 0
    310283
    311 def random_tests(n=15, nrows_max=50, ncols_max=50, nsymbols_max=20, perms_per_matrix=10, density_range=(.1,.9)):
     284def random_tests(n=10, nrows_max=50, ncols_max=50, nsymbols_max=10, perms_per_matrix=5, density_range=(.1,.9)):
    312285    """
    313286    Tests to make sure that C(gamma(M)) == C(M) for random permutations gamma
    314287    and random matrices M, and that M.is_isomorphic(gamma(M)) returns an
  • sage/groups/perm_gps/partn_ref/refinement_python.pxd

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_python.pxd
    a b  
    11
    22#*****************************************************************************
    3 #      Copyright (C) 2006 - 2009 Robert L. Miller <rlmillster@gmail.com>
     3#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    44#
    55# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    66#                         http://www.gnu.org/licenses/
     
    1313
    1414
    1515from sage.rings.integer cimport Integer
    16 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     16from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab
    1717from double_coset cimport double_coset
    1818
    1919
  • sage/groups/perm_gps/partn_ref/refinement_python.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/groups/perm_gps/partn_ref/refinement_python.pyx
    a b  
    2020"""
    2121
    2222#*****************************************************************************
    23 #      Copyright (C) 2006 - 2009 Robert L. Miller <rlmillster@gmail.com>
     23#      Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com>
    2424#
    2525# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    2626#                         http://www.gnu.org/licenses/
     
    376376        self.rari_fn = rari_fn
    377377        self.cs_fn = cs_fn
    378378
    379 cdef bint all_children_are_equivalent_python(PartitionStack *PS, object S):
     379cdef bint all_children_are_equivalent_python(PartitionStack *PS, void *S):
    380380    """
    381381    Python conversion of all_children_are_equivalent function.
    382382    """
    383383    cdef PythonPartitionStack Py_PS = PythonPartitionStack(PS.degree)
     384    cdef object S_obj = <object> S
    384385    PS_copy_from_to(PS, Py_PS.c_ps)
    385     return S.acae_fn(Py_PS, S.obj)
     386    return S_obj.acae_fn(Py_PS, S_obj.obj)
    386387
    387 cdef int refine_and_return_invariant_python(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     388cdef int refine_and_return_invariant_python(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len):
    388389    """
    389390    Python conversion of refine_and_return_invariant function.
    390391    """
    391392    cdef PythonPartitionStack Py_PS = PythonPartitionStack(PS.degree)
     393    cdef object S_obj = <object> S
    392394    PS_copy_from_to(PS, Py_PS.c_ps)
    393395    cdef int i
    394     cdef list ctrb_py = [cells_to_refine_by[i] for i from 0 <= i < PS.degree]
    395     return S.rari_fn(Py_PS, S.obj, ctrb_py)
     396    cdef list ctrb_py = [cells_to_refine_by[i] for i from 0 <= i < ctrb_len]
     397    return S_obj.rari_fn(Py_PS, S_obj.obj, ctrb_py)
    396398
    397 cdef int compare_structures_python(int *gamma_1, int *gamma_2, object S1, object S2):
     399cdef int compare_structures_python(int *gamma_1, int *gamma_2, void *S1, void *S2):
    398400    """
    399401    Python conversion of compare_structures function.
    400402    """
    401403    cdef int i
    402     cdef list gamma_1_py = [gamma_1[i] for i from 0 <= i < S1.degree]
    403     cdef list gamma_2_py = [gamma_2[i] for i from 0 <= i < S1.degree]
    404     return S1.cs_fn(gamma_1_py, gamma_2_py, S1.obj, S2.obj)
     404    cdef object S1_obj = <object> S1, S2_obj = <object> S2
     405    cdef list gamma_1_py = [gamma_1[i] for i from 0 <= i < S1_obj.degree]
     406    cdef list gamma_2_py = [gamma_2[i] for i from 0 <= i < S1_obj.degree]
     407    return S1_obj.cs_fn(gamma_1_py, gamma_2_py, S1_obj.obj, S2_obj.obj)
    405408
    406409def aut_gp_and_can_lab_python(S, partition, n,
    407410    all_children_are_equivalent,
     
    455458
    456459    """
    457460    obj_wrapper = PythonObjectWrapper(S, all_children_are_equivalent, refine_and_return_invariant, compare_structures, n)
    458     cdef aut_gp_and_can_lab_return *output
     461    cdef aut_gp_and_can_lab *output
    459462    cdef PythonPartitionStack Py_PS = PythonPartitionStack(n)
    460463    cdef int i, j
    461464    cdef Integer I
    462465
    463     cdef int **part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     466    cdef PartitionStack *part = PS_from_list(partition)
    464467    if part is NULL:
    465468        raise MemoryError
    466     for i from 0 <= i < len(partition):
    467         part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
    468         if part[i] is NULL:
    469             for j from 0 <= j < i:
    470                 sage_free(part[j])
    471             sage_free(part)
    472             raise MemoryError
    473         for j from 0 <= j < len(partition[i]):
    474             part[i][j] = partition[i][j]
    475         part[i][len(partition[i])] = -1
    476     part[len(partition)] = NULL
    477469
    478     output = get_aut_gp_and_can_lab(obj_wrapper, part, n,
     470    output = get_aut_gp_and_can_lab(<void *> obj_wrapper, part, n,
    479471        &all_children_are_equivalent_python,
    480472        &refine_and_return_invariant_python,
    481473        &compare_structures_python,
    482         canonical_label, base, order)
     474        canonical_label, NULL)
    483475
    484476    list_of_gens = []
    485477    for i from 0 <= i < output.num_gens:
     
    488480    if canonical_label:
    489481        return_tuple.append([output.relabeling[i] for i from 0 <= i < n])
    490482    if base:
    491         return_tuple.append([output.base[i] for i from 0 <= i < output.base_size])
     483        return_tuple.append([output.group.base_orbits[i][0] for i from 0 <= i < output.group.base_size])
    492484    if order:
    493485        I = Integer()
    494         mpz_set(I.value, output.order)
     486        SC_order(output.group, 0, I.value)
    495487        return_tuple.append(I)
    496     for i from 0 <= i < len(partition):
    497         sage_free(part[i])
    498     sage_free(part)
    499     mpz_clear(output.order)
     488    PS_dealloc(part)
    500489    sage_free(output.generators)
    501     if base:
    502         sage_free(output.base)
     490    SC_dealloc(output.group)
    503491    if canonical_label:
    504492        sage_free(output.relabeling)
    505493    sage_free(output)
     
    559547    obj_wrapper1 = PythonObjectWrapper(S1, all_children_are_equivalent, refine_and_return_invariant, compare_structures, n)
    560548    obj_wrapper2 = PythonObjectWrapper(S2, all_children_are_equivalent, refine_and_return_invariant, compare_structures, n)
    561549
    562     cdef int **part = <int **> sage_malloc((len(partition1)+1) * sizeof(int *))
     550    cdef PartitionStack *part = PS_from_list(partition1)
    563551    cdef int *ordering = <int *> sage_malloc(n * sizeof(int))
    564552    if part is NULL or ordering is NULL:
    565         if part is not NULL:
    566             sage_free(part)
    567         if ordering is not NULL:
    568             sage_free(ordering)
     553        if part is not NULL: PS_dealloc(part)
     554        if ordering is not NULL: sage_free(ordering)
    569555        raise MemoryError
    570     for i from 0 <= i < len(partition1):
    571         part[i] = <int *> sage_malloc((len(partition1[i])+1) * sizeof(int))
    572         if part[i] is NULL:
    573             for j from 0 <= j < i:
    574                 sage_free(part[j])
    575             sage_free(part)
    576             raise MemoryError
    577         for j from 0 <= j < len(partition1[i]):
    578             part[i][j] = partition1[i][j]
    579         part[i][len(partition1[i])] = -1
    580     part[len(partition1)] = NULL
    581556    for i from 0 <= i < n:
    582557        ordering[i] = ordering2[i]
    583558   
    584     cdef int *output = double_coset(obj_wrapper1, obj_wrapper2, part, ordering, n,
     559    cdef int *output = double_coset(<void *> obj_wrapper1, <void *> obj_wrapper2,
     560        part, ordering, n,
    585561        &all_children_are_equivalent_python,
    586562        &refine_and_return_invariant_python,
    587         &compare_structures_python)
     563        &compare_structures_python, NULL)
    588564
    589     for i from 0 <= i < len(partition1):
    590         sage_free(part[i])
    591     sage_free(part)
     565    PS_dealloc(part)
    592566    sage_free(ordering)
    593567   
    594568    if output is NULL:
  • sage/sets/disjoint_set.pyx

    diff -r 872cd84f620c -r 8ddd7727b639 sage/sets/disjoint_set.pyx
    a b  
    6060include '../groups/perm_gps/partn_ref/data_structures_pyx.pxi'
    6161
    6262import itertools
    63 from sage.rings.all import Integer
     63from sage.rings.integer import Integer
    6464from sage.structure.sage_object cimport SageObject
    6565
    6666def DisjointSet(arg):