Ticket #10804: trac_10804.patch
| File trac_10804.patch, 191.3 KB (added by , 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 435 435 436 436 Extension('sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label', 437 437 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"]), 439 442 440 443 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"]), 442 449 443 450 Extension('sage.groups.perm_gps.partn_ref.refinement_binary', 444 451 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"]), 446 456 447 457 Extension('sage.groups.perm_gps.partn_ref.refinement_graphs', 448 458 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"]), 450 463 451 464 Extension('sage.groups.perm_gps.partn_ref.refinement_lists', 452 465 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"]), 454 470 455 471 Extension('sage.groups.perm_gps.partn_ref.refinement_matrices', 456 472 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"]), 458 477 459 478 Extension('sage.groups.perm_gps.partn_ref.refinement_python', 460 479 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"]), 462 484 463 485 ################################ 464 486 ## … … 1570 1592 ################################ 1571 1593 1572 1594 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"]), 1574 1600 1575 1601 ################################ 1576 1602 ## -
sage/graphs/matchpoly.pyx
diff -r 872cd84f620c -r 8ddd7727b639 sage/graphs/matchpoly.pyx
a b 266 266 267 267 # Computing the signless matching polynomial 268 268 269 _sig_on269 sig_on() 270 270 delete_and_add(edges, nverts, nedges, nverts, 0, pol) 271 _sig_off271 sig_off() 272 272 273 273 # Building the actual matching polynomial 274 274 -
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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # 5 5 # Distributed under the terms of the GNU General Public License (GPL) 6 6 # http://www.gnu.org/licenses/ … … 12 12 13 13 from sage.rings.integer cimport Integer 14 14 15 cdef struct aut_gp_and_can_lab _return:15 cdef struct aut_gp_and_can_lab: 16 16 int *generators 17 17 int num_gens 18 StabilizerChain *group 18 19 int *relabeling 19 int *base20 int base_size21 mpz_t order22 20 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) 21 cdef 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 27 26 28 27 29 28 -
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 36 36 37 37 Signature: 38 38 39 \code{int refine_and_return_invariant(PartitionStack *PS, objectS, 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)} 40 40 41 41 This function should split up cells in the partition at the top of the 42 42 partition stack in such a way that any automorphism that respects the … … 58 58 59 59 Signature: 60 60 61 \code{int compare_structures(int *gamma_1, int *gamma_2, object S1, objectS2)}61 \code{int compare_structures(int *gamma_1, int *gamma_2, void *S1, void *S2)} 62 62 63 63 This function must implement a total ordering on the set of objects of fixed 64 64 order. Return: … … 77 77 78 78 Signature: 79 79 80 \code{bint all_children_are_equivalent(PartitionStack *PS, objectS)}80 \code{bint all_children_are_equivalent(PartitionStack *PS, void *S)} 81 81 82 82 This function must return False unless it is the case that each discrete 83 83 partition finer than the top of the partition stack is equivalent to the … … 96 96 """ 97 97 98 98 #***************************************************************************** 99 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>99 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 100 100 # 101 101 # Distributed under the terms of the GNU General Public License (GPL) 102 102 # http://www.gnu.org/licenses/ 103 103 #***************************************************************************** 104 104 105 105 include 'data_structures_pyx.pxi' # includes bitsets 106 include 'sage/ext/interrupt.pxi' 106 107 107 108 cdef inline int cmp(int a, int b): 108 109 if a < b: return -1 … … 111 112 112 113 # Functions 113 114 114 cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, objectS):115 cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, void *S): 115 116 return 0 116 117 117 cdef int refine_and_return_invariant_trivial(PartitionStack *PS, objectS, int *cells_to_refine_by, int ctrb_len):118 cdef int refine_and_return_invariant_trivial(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len): 118 119 return 0 119 120 120 cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S1, objectS2):121 cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, void *S1, void *S2): 121 122 return 0 122 123 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):124 def 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): 125 126 """ 126 127 sage: tttt = sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label.test_get_aut_gp_and_can_lab_trivially 127 128 sage: tttt() … … 142 143 1 143 144 144 145 """ 145 cdef int i, j 146 cdef aut_gp_and_can_lab_return *output 146 cdef aut_gp_and_can_lab *output 147 147 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) 157 155 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) 162 158 sage_free(output.generators) 163 if base:164 sage_free(output.base)165 159 if canonical_label: 166 160 sage_free(output.relabeling) 167 161 sage_free(output) 168 162 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), 163 def 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 206 cdef 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 215 def 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 289 cdef 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), 171 292 int (*refine_and_return_invariant)\ 172 (PartitionStack *PS, objectS, int *cells_to_refine_by, int ctrb_len),173 int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, objectS2),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: 175 296 """ 176 297 Traverse the search space for subgroup/canonical label calculation. 177 298 178 299 INPUT: 179 300 S -- pointer to the structure 180 partition -- arrayrepresenting a partition of the points301 partition -- PartitionStack representing a partition of the points 181 302 len_partition -- length of the partition 182 303 n -- the number of points (points are assumed to be 0,1,...,n-1) 183 304 canonical_label -- whether to search for canonical label - if True, return 184 305 the permutation taking S to its canonical label 185 base -- if True, return the base for which the given generating set is a186 strong generating set187 order -- if True, return the order of the automorphism group188 306 all_children_are_equivalent -- pointer to a function 189 307 INPUT: 190 308 PS -- pointer to a partition stack … … 207 325 OUTPUT: 208 326 int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp), 209 327 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 210 333 OUTPUT: 211 pointer to a aut_gp_and_can_lab _returnstruct334 pointer to a aut_gp_and_can_lab struct 212 335 213 336 """ 214 337 cdef PartitionStack *current_ps, *first_ps, *label_ps … … 222 345 cdef int label_and_current_indicator_same = -1 223 346 cdef int compared_current_and_label_indicators 224 347 225 cdef OrbitPartition *orbits_of_subgroup 226 cdef mpz_t subgroup_size 348 cdef OrbitPartition *orbits_of_subgroup, *orbits_of_permutation 227 349 cdef int subgroup_primary_orbit_size = 0 228 350 cdef int minimal_in_primary_orbit 229 351 cdef int size_of_generator_array = 2*n*n … … 235 357 236 358 cdef bitset_t *vertices_to_split 237 359 cdef bitset_t vertices_have_been_reduced 238 cdef int *permutation, * id_perm, *cells_to_refine_by360 cdef int *permutation, *label_perm, *id_perm, *cells_to_refine_by 239 361 cdef int *vertices_determining_current_stack 362 cdef int *perm_stack 363 cdef StabilizerChain *group = NULL, *old_group, *tmp_gp 240 364 241 365 cdef int i, j, k 242 366 cdef bint discrete, automorphism, update_label 243 367 cdef bint backtrack, new_vertex, narrow, mem_err = 0 244 368 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)) 246 370 if output is NULL: 247 371 raise MemoryError 248 372 output.group = SC_new(n) 373 if output.group is NULL: 374 sage_free(output) 375 raise MemoryError 249 376 if n == 0: 250 377 output.generators = NULL 251 378 output.num_gens = 0 252 379 output.relabeling = NULL 253 output.base = NULL254 output.base_size = 0255 mpz_init_set_si(output.order, 1)256 380 return output 257 381 258 382 # Allocate: 259 mpz_init_set_si(subgroup_size, 1)260 261 383 output.generators = <int *> sage_malloc( size_of_generator_array * sizeof(int) ) 262 384 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) 266 402 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 269 409 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 345 422 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) 376 425 bitset_init(vertices_have_been_reduced, n) 377 426 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 390 437 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 428 442 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) 456 460 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 458 500 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 463 511 464 first_meets_current = current_ps.depth465 first_kids_are_same = current_ps.depth466 first_and_current_indicator_same = current_ps.depth467 first_ps = PS_copy(current_ps)468 if canonical_label:469 label_meets_current = current_ps.depth470 label_and_current_indicator_same = current_ps.depth471 compared_current_and_label_indicators = 0472 label_ps = PS_copy(current_ps)473 current_ps.depth -= 1474 475 512 # Main loop: 476 513 while current_ps.depth != -1: 477 514 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 478 538 # If necessary, update label_meets_current 479 539 if canonical_label and label_meets_current > current_ps.depth: 480 540 label_meets_current = current_ps.depth … … 517 577 # are looking at a new primary orbit (corresponding to a 518 578 # larger subgroup in the stabilizer chain). 519 579 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] 524 581 while 1: 525 582 i = vertices_determining_current_stack[current_ps.depth] 526 583 # 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: 530 586 subgroup_primary_orbit_size += 1 531 587 # Look for a new point to split. 532 588 i += 1 … … 550 606 j += 1 551 607 if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1: 552 608 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)555 609 subgroup_primary_orbit_size = 0 556 610 current_ps.depth -= 1 557 611 break … … 569 623 # II. Refine down to a discrete partition, or until 570 624 # we leave the part of the tree we are interested in 571 625 discrete = 0 626 backtrack = 0 572 627 while 1: 573 628 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) 575 638 PS_move_all_mins_to_front(current_ps) 576 639 if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]: 577 640 first_and_current_indicator_same += 1 … … 586 649 if compared_current_and_label_indicators > 0: 587 650 label_indicators[i] = current_indicators[i] 588 651 if first_and_current_indicator_same < current_ps.depth and (not canonical_label or compared_current_and_label_indicators < 0): 652 backtrack = 1 589 653 break 590 654 if PS_is_discrete(current_ps): 591 655 discrete = 1 … … 597 661 598 662 # III. Check for automorphisms and labels 599 663 automorphism = 0 600 backtrack = 1 - discrete601 664 if discrete: 602 665 if current_ps.depth == first_and_current_indicator_same: 603 666 PS_get_perm_from(current_ps, first_ps, permutation) 604 667 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 606 671 if not automorphism: 607 672 if (not canonical_label) or (compared_current_and_label_indicators < 0): 608 673 backtrack = 1 … … 612 677 if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth): 613 678 update_label = 1 614 679 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) 616 687 if i > 0: 617 688 update_label = 1 618 689 elif i < 0: 619 690 backtrack = 1 620 691 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) 622 697 automorphism = 1 623 698 if update_label: 624 699 PS_copy_from_to(current_ps, label_ps) 700 if input_group is not NULL: 701 memcpy(label_perm, permutation, n*sizeof(int)) 625 702 compared_current_and_label_indicators = 0 626 703 label_meets_current = current_ps.depth 627 704 label_and_current_indicator_same = current_ps.depth … … 632 709 index_in_fp_and_mcr += 1 633 710 bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr]) 634 711 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]) 635 716 for i from 0 <= i < n: 636 717 if permutation[i] == i: 637 718 bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i) 638 719 bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i) 639 720 else: 640 721 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]) 650 724 if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser 725 # add permutation as a generator of the automorphism group 651 726 if n*output.num_gens == size_of_generator_array: 652 727 # must double its size 653 728 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) ) 655 730 if output.generators is NULL: 656 731 mem_err = True 657 break# main loop732 continue # main loop 658 733 j = n*output.num_gens 659 734 for i from 0 <= i < n: 660 735 output.generators[j + i] = permutation[i] 661 736 output.num_gens += 1 737 if SC_update_tree(output.group, first_meets_current, permutation, 1): 738 mem_err = True 739 continue # main loop 662 740 if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit: 663 741 current_ps.depth = first_meets_current 664 742 continue # main loop … … 691 769 692 770 # End of main loop. 693 771 694 mpz_init_set(output.order, subgroup_size)695 772 if canonical_label: 696 773 for i from 0 <= i < n: 697 774 output.relabeling[label_ps.entries[i]] = i 698 775 699 776 # 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]) 705 789 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) 719 791 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 734 794 735 795 736 796 -
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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # 5 5 # Distributed under the terms of the GNU General Public License (GPL) 6 6 # http://www.gnu.org/licenses/ … … 10 10 include '../../../ext/stdsage.pxi' 11 11 include '../../../misc/bitset_pxd.pxi' 12 12 13 cdef extern from "stdlib.h": 14 int rand() 15 16 cdef extern from "FLINT/long_extras.h": 17 int z_isprime(unsigned long n) 18 13 19 cdef struct OrbitPartition: 14 20 # Disjoint set data structure for representing the orbits of the generators 15 21 # so far found. Also keeps track of the minimum elements of the cells and … … 31 37 int depth 32 38 int degree 33 39 40 cdef struct StabilizerChain: 41 # A representation of a permutation group acting on 0, 1, ..., degree-1. 42 int degree 43 int base_size 34 44 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 1 1 r""" 2 2 Data structures 3 3 4 This module implements TODO... 4 This module implements basic data structures essential to the rest of the 5 partn_ref module. 5 6 6 REFERENCE :7 REFERENCES: 7 8 8 9 [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium, 9 10 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. 10 16 11 17 """ 12 18 13 19 #***************************************************************************** 14 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>20 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 15 21 # 16 22 # Distributed under the terms of the GNU General Public License (GPL) 17 23 # http://www.gnu.org/licenses/ … … 19 25 20 26 include '../../../misc/bitset.pxi' 21 27 22 cdef inline int min(int a, int b):23 if a < b: return a24 else: return b28 cdef extern from "math.h": 29 float log(float x) 30 float ceil(float x) 25 31 26 cdef inline int max(int a, int b): 27 if a > b: return a 28 else: return b 32 from sage.groups.perm_gps.permgroup import PermutationGroup 33 from sage.rings.integer cimport Integer 29 34 30 35 # OrbitPartitions 31 36 … … 37 42 cdef int i 38 43 cdef OrbitPartition *OP = <OrbitPartition *> \ 39 44 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 42 50 OP.degree = n 43 51 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 59 cdef inline OP_clear(OrbitPartition *OP): 60 cdef int i, n = OP.degree 55 61 for i from 0 <= i < n: 56 62 OP.parent[i] = i 57 63 OP.rank[i] = 0 58 64 OP.mcr[i] = i 59 65 OP.size[i] = 1 60 return OP61 66 62 67 cdef 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) 67 70 sage_free(OP) 68 return 069 71 70 72 cdef inline int OP_find(OrbitPartition *OP, int n): 71 73 """ … … 227 229 cdef int i 228 230 cdef PartitionStack *PS = <PartitionStack *> \ 229 231 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) 237 236 return NULL 238 PS.depth = 0 237 PS.entries = int_array 238 PS.levels = int_array + n 239 PS.depth = 0 239 240 PS.degree = n 240 241 if unit_partition: 241 242 for i from 0 <= i < n-1: … … 250 251 Allocate and return a pointer to a copy of PartitionStack PS. Returns a null 251 252 pointer in the case of an allocation failure. 252 253 """ 253 cdef int i 254 cdef int i, n = PS.degree 255 254 256 cdef PartitionStack *PS2 = <PartitionStack *> \ 255 257 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) 263 262 return NULL 263 PS2.entries = int_array 264 PS2.levels = int_array + n 264 265 PS_copy_from_to(PS, PS2) 265 266 return PS2 266 267 … … 270 271 """ 271 272 PS2.depth = PS.depth 272 273 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) ) 276 275 277 cdef inline PartitionStack *PS_from_list(object L):276 cdef PartitionStack *PS_from_list(list L): 278 277 """ 279 278 Allocate and return a pointer to a PartitionStack representing L. Returns a 280 279 null pointer in the case of an allocation failure. … … 283 282 for cell from 0 <= cell < num_cells: 284 283 n += len(L[cell]) 285 284 cdef PartitionStack *PS = PS_new(n, 0) 286 cell = 0287 285 if PS is NULL: 288 return PS289 while 1:286 return NULL 287 for cell from 0 <= cell < num_cells: 290 288 cur_len = len(L[cell]) 291 289 for i from 0 <= i < cur_len: 292 290 PS.entries[cur_start + i] = L[cell][i] 293 291 PS.levels[cur_start + i] = n 294 292 PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1) 295 293 cur_start += cur_len 296 cell += 1297 if cell == num_cells:298 PS.levels[cur_start-1] = -1299 break300 294 PS.levels[cur_start-1] = 0 295 if n > 0: 296 PS.levels[n-1] = -1 297 PS.depth = 0 298 PS.degree = n 301 299 return PS 302 300 303 cdef inline intPS_dealloc(PartitionStack *PS):304 sage_free(PS.entries)305 sage_free(PS.levels)301 cdef inline PS_dealloc(PartitionStack *PS): 302 if PS is not NULL: 303 sage_free(PS.entries) 306 304 sage_free(PS) 307 return 0308 305 309 cdef inline intPS_print(PartitionStack *PS):306 cdef PS_print(PartitionStack *PS): 310 307 """ 311 308 Print a visual representation of PS. 312 309 """ … … 315 312 PS_print_partition(PS, i) 316 313 print 317 314 318 cdef inline intPS_print_partition(PartitionStack *PS, int k):315 cdef PS_print_partition(PartitionStack *PS, int k): 319 316 """ 320 317 Print the partition at depth k. 321 318 """ … … 349 346 ncells += 1 350 347 return ncells 351 348 352 cdef inline intPS_move_min_to_front(PartitionStack *PS, int start, int end):349 cdef inline PS_move_min_to_front(PartitionStack *PS, int start, int end): 353 350 """ 354 351 Makes sure that the first element of the segment of entries i with 355 352 start <= i <= end is minimal. … … 383 380 """ 384 381 cdef int i, cur_start = 0 385 382 for i from 0 <= i < PS.degree: 386 if PS.levels[i] >= PS.depth:383 if PS.levels[i] == PS.depth: 387 384 PS.levels[i] += 1 388 385 if PS.levels[i] < PS.depth: 389 386 PS_move_min_to_front(PS, cur_start, i) … … 399 396 PS_move_min_to_front(PS, cur_start, i) 400 397 cur_start = i+1 401 398 402 cdef in line int PS_split_point(PartitionStack *PS, int v):399 cdef int PS_split_point(PartitionStack *PS, int v): 403 400 """ 404 401 Detaches the point v from the cell it is in, putting the singleton cell 405 402 of just v in front. Returns the position where v is now located. … … 430 427 PS.levels[i] = PS.depth 431 428 return i 432 429 433 cdef in line int PS_first_smallest(PartitionStack *PS, bitset_t b):430 cdef int PS_first_smallest(PartitionStack *PS, bitset_t b): 434 431 """ 435 432 Find the first occurrence of the smallest cell of size greater than one, 436 433 store its entries to b, and return its minimum element. … … 454 451 i += 1 455 452 return PS.entries[location] 456 453 454 cdef 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 457 475 cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma): 458 476 """ 459 477 Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i] … … 463 481 for i from 0 <= i < PS1.degree: 464 482 gamma[PS2.entries[i]] = PS1.entries[i] 465 483 484 cdef 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 466 494 def PS_represent(partition, splits): 467 495 """ 468 496 Demonstration and testing. … … 514 542 Finding first smallest: 515 543 Minimal element is 5, bitset is: 516 544 0000010001 545 Finding element 1: 546 Location is: 5 547 Bitset is: 548 0100000000 517 549 Deallocating PartitionStacks. 518 550 Done. 519 551 … … 591 623 print "Minimal element is %d, bitset is:"%i 592 624 print bitset_string(b) 593 625 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) 594 632 print "Deallocating PartitionStacks." 595 633 PS_dealloc(PS) 596 634 PS_dealloc(PS2) 597 635 print "Done." 598 636 637 # StabilizerChains 638 639 cdef enum: 640 default_num_gens = 8 641 default_num_bits = 64 642 643 cdef 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 723 cdef 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 764 cdef 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 808 cdef 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 827 cdef 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 861 cdef 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 875 cdef 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 903 cdef 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 917 cdef 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 925 cdef 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 933 cdef 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 941 cdef inline SC_identify(int *perm, int degree): 942 cdef int i 943 for i from 0 <= i < degree: 944 perm[i] = i 945 946 cdef 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 961 cdef 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 975 cdef 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 993 cdef 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 1016 cdef 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 1028 cdef 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 1046 cdef 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 1056 cdef 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 1082 cdef 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 1100 cdef 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 1117 cdef 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 1157 cdef 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 1204 cdef 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 1307 cdef 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 1324 cdef 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 1327 cdef 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 1339 cdef 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 1360 cdef 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 1371 cdef 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 1444 def 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 599 1844 # Functions 600 1845 601 cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S, 1846 cdef 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 1889 cdef inline int split_point_and_refine(PartitionStack *PS, int v, void *S, 602 1890 int (*refine_and_return_invariant)\ 603 (PartitionStack *PS, objectS, int *cells_to_refine_by, int ctrb_len),1891 (PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len), 604 1892 int *cells_to_refine_by): 605 1893 """ 606 1894 Make the partition stack one longer by copying the last partition in the … … 613 1901 S -- the structure 614 1902 refine_and_return_invariant -- the refinement function provided 615 1903 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 616 1906 617 1907 """ 618 1908 PS.depth += 1 … … 620 1910 cells_to_refine_by[0] = PS_split_point(PS, v) 621 1911 return refine_and_return_invariant(PS, S, cells_to_refine_by, 1) 622 1912 1913 cdef 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 1927 cdef 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 1966 cdef 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 1977 cdef 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 1987 cdef 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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # 5 5 # Distributed under the terms of the GNU General Public License (GPL) 6 6 # http://www.gnu.org/licenses/ … … 12 12 13 13 from sage.rings.integer cimport Integer 14 14 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)) 15 cdef 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 *) 19 20 20 21 21 22 -
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 19 19 20 20 Signature: 21 21 22 \code{int refine_and_return_invariant(PartitionStack *PS, objectS, 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)} 23 23 24 24 This function should split up cells in the partition at the top of the 25 25 partition stack in such a way that any automorphism that respects the … … 41 41 42 42 Signature: 43 43 44 \code{int compare_structures(int *gamma_1, int *gamma_2, object S1, objectS2)}44 \code{int compare_structures(int *gamma_1, int *gamma_2, void *S1, void *S2)} 45 45 46 46 This function must implement a total ordering on the set of objects of fixed 47 47 order. Return: … … 60 60 61 61 Signature: 62 62 63 \code{bint all_children_are_equivalent(PartitionStack *PS, objectS)}63 \code{bint all_children_are_equivalent(PartitionStack *PS, void *S)} 64 64 65 65 This function must return False unless it is the case that each discrete 66 66 partition finer than the top of the partition stack is equivalent to the … … 83 83 """ 84 84 85 85 #***************************************************************************** 86 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>86 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 87 87 # 88 88 # Distributed under the terms of the GNU General Public License (GPL) 89 89 # http://www.gnu.org/licenses/ … … 96 96 elif a == b: return 0 97 97 else: return 1 98 98 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: continue104 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 0107 return 1108 109 99 # Functions 110 100 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), 101 cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, void *S): 102 return 0 103 104 cdef int refine_and_return_invariant_trivial(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len): 105 return 0 106 107 cdef 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 116 def 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 196 cdef int *double_coset(void *S1, void *S2, PartitionStack *partition1, int *ordering2, 197 int n, bint (*all_children_are_equivalent)(PartitionStack *PS, void *S), 113 198 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): 116 202 """ 117 203 Traverse the search space for double coset calculation. 118 204 119 205 INPUT: 120 206 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 122 209 ordering2 -- an ordering of the points of S2 representing a second partition 123 210 n -- the number of points (points are assumed to be 0,1,...,n-1) 124 211 all_children_are_equivalent -- pointer to a function … … 143 230 OUTPUT: 144 231 int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp), 145 232 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 146 239 OUTPUT: 147 240 If S1 and S2 are isomorphic, a pointer to an integer array representing an 148 241 isomorphism. Otherwise, a NULL pointer. … … 168 261 cdef bitset_t vertices_have_been_reduced 169 262 cdef int *permutation, *id_perm, *cells_to_refine_by 170 263 cdef int *vertices_determining_current_stack 264 cdef int *perm_stack 265 cdef StabilizerChain *group, *old_group, *tmp_gp 171 266 172 267 cdef int *output 173 268 … … 178 273 if n == 0: 179 274 return NULL 180 275 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) ) 192 277 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 194 288 current_ps = PS_new(n, 0) 195 left_ps = PS_new(n, 0)196 289 orbits_of_subgroup = OP_new(n) 197 290 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 233 294 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 237 305 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) 268 308 bitset_init(vertices_have_been_reduced, n) 269 309 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 325 311 326 312 cdef bint possible = 1 327 313 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 328 455 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 353 471 354 first_ps = NULL355 # Refine down to a discrete partition356 while not PS_is_discrete(left_ps):357 i = left_ps.depth358 k = PS_first_smallest(left_ps, vertices_to_split[i]) # writes to vertices_to_split, but this is never used359 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.depth363 vertices_determining_current_stack[i] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])364 while possible:365 i = current_ps.depth366 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 = 0369 elif not stacks_are_equivalent(left_ps, current_ps):370 possible = 0371 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 + 1375 break376 current_ps.depth -= 1 # reset from split_point_and_refine377 while current_ps.depth >= 0: # not possible, so look for another378 i = current_ps.depth379 j = vertices_determining_current_stack[i] + 1380 j = bitset_next(vertices_to_split[i], j)381 if j == -1:382 current_ps.depth -= 1 # backtrack383 else:384 possible = 1385 vertices_determining_current_stack[i] = j386 break # found another387 388 if possible:389 if compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:390 unknown = 0391 else:392 first_meets_current = current_ps.depth393 first_kids_are_same = current_ps.depth394 first_ps = PS_copy(current_ps)395 current_ps.depth -= 1396 397 472 # Main loop: 398 473 while possible and unknown and current_ps.depth != -1: 399 474 … … 484 559 while 1: 485 560 i = current_ps.depth 486 561 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) 488 572 PS_move_all_mins_to_front(current_ps) 489 573 if indicators[i] != k: 490 574 possible = 0 491 575 elif not stacks_are_equivalent(left_ps, current_ps): 492 576 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 493 583 if not possible: 494 584 j = vertices_determining_current_stack[i] + 1 495 585 j = bitset_next(vertices_to_split[i], j) … … 504 594 break 505 595 if PS_is_discrete(current_ps): 506 596 break 507 vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])508 597 bitset_unset(vertices_have_been_reduced, current_ps.depth) 509 598 if not all_children_are_equivalent(current_ps, S2): 510 599 current_kids_are_same = current_ps.depth + 1 … … 512 601 # III. Check for automorphisms and isomorphisms 513 602 automorphism = 0 514 603 if possible: 515 PS_get_perm_from( current_ps, first_ps, permutation)604 PS_get_perm_from(first_ps, current_ps, permutation) 516 605 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 518 609 if not automorphism and possible: 519 610 # if we get here, discrete must be true 520 611 if current_ps.depth != left_ps.depth: 521 612 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 525 619 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 527 627 if automorphism: 528 628 if index_in_fp_and_mcr < len_of_fp_and_mcr - 1: 529 629 index_in_fp_and_mcr += 1 … … 582 682 output = NULL 583 683 584 684 # 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]) 590 690 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) 599 697 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 604 698 return output 605 699 606 700 -
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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # 5 5 # Distributed under the terms of the GNU General Public License (GPL) 6 6 # http://www.gnu.org/licenses/ … … 11 11 include 'data_structures_pxd.pxi' # includes bitsets 12 12 13 13 from 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 _return14 from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab 15 15 from double_coset cimport double_coset 16 16 17 17 cdef class BinaryCodeStruct: … … 22 22 cdef PartitionStack *word_ps 23 23 cdef int *alpha # length nwords + degree 24 24 cdef int *scratch # length 3*nwords + 3*degree + 2 25 cdef aut_gp_and_can_lab _return*output25 cdef aut_gp_and_can_lab *output 26 26 cdef int (*ith_word)(BinaryCodeStruct self, int, bitset_s *) 27 27 28 28 cdef class LinearBinaryCodeStruct(BinaryCodeStruct): … … 36 36 cdef bitset_s *scratch_bitsets # length 4*nwords + 1 37 37 cdef int ith_word_nonlinear(BinaryCodeStruct, int, bitset_s *) 38 38 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)39 cdef int refine_by_bip_degree(PartitionStack *, void *, int *, int) 40 cdef int compare_linear_codes(int *, int *, void *, void *) 41 cdef int compare_nonlinear_codes(int *, int *, void *, void *) 42 cdef bint all_children_are_equivalent(PartitionStack *, void *) 43 43 cdef inline int word_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *) 44 44 cdef inline int col_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *) 45 cdef inline int sort_by_function (PartitionStack *, int, int *, int *, int *, int)45 cdef inline int sort_by_function_codes(PartitionStack *, int, int *, int *, int *, int) 46 46 47 47 48 48 -
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 16 16 """ 17 17 18 18 #***************************************************************************** 19 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>19 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 20 20 # 21 21 # Distributed under the terms of the GNU General Public License (GPL) 22 22 # http://www.gnu.org/licenses/ … … 222 222 17868913969917295853568000000 223 223 224 224 """ 225 cdef int **part, i, j 225 cdef int i, n = self.degree 226 cdef PartitionStack *part 226 227 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) 229 231 if part is NULL: 230 232 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 242 234 self.first_time = 1 243 235 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) 245 237 246 for i from 0 <= i < len(partition): 247 sage_free(part[i]) 248 sage_free(part) 238 PS_dealloc(part) 249 239 250 240 def automorphism_group(self): 251 241 """ … … 262 252 263 253 """ 264 254 cdef int i, j 265 cdef object generators, base255 cdef list generators, base 266 256 cdef Integer order 267 257 if self.output is NULL: 268 258 self.run() … … 270 260 for i from 0 <= i < self.output.num_gens: 271 261 generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree]) 272 262 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] 275 265 return generators, order, base 276 266 277 267 def canonical_relabeling(self): … … 311 301 [0, 1, 2, 5, 3, 4] 312 302 313 303 """ 314 cdef int **part, i, j304 cdef int i, n = self.degree 315 305 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) 318 308 ordering = <int *> sage_malloc(self.degree * sizeof(int)) 319 309 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) 321 311 if ordering is not NULL: sage_free(ordering) 322 312 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: 335 314 ordering[i] = i 336 315 self.first_time = 1 337 316 other.first_time = 1 338 317 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) 340 319 341 for i from 0 <= i < len(partition): 342 sage_free(part[i]) 343 sage_free(part) 320 PS_dealloc(part) 344 321 sage_free(ordering) 345 322 346 323 if output is NULL: 347 324 return False 348 325 else: 349 output_py = [output[i] for i from 0 <= i < self.degree]326 output_py = [output[i] for i from 0 <= i < n] 350 327 sage_free(output) 351 328 return output_py 352 329 … … 361 338 sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps) 362 339 sage_free(self.alpha); sage_free(self.scratch) 363 340 if self.output is not NULL: 364 mpz_clear(self.output.order)365 341 sage_free(self.output.generators) 366 sage_free(self.output.base)342 SC_dealloc(self.output.group) 367 343 sage_free(self.output.relabeling) 368 344 sage_free(self.output) 369 345 … … 468 444 sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps) 469 445 sage_free(self.alpha); sage_free(self.scratch) 470 446 if self.output is not NULL: 471 mpz_clear(self.output.order)472 447 sage_free(self.output.generators) 473 sage_free(self.output.base)448 SC_dealloc(self.output.group) 474 449 sage_free(self.output.relabeling) 475 450 sage_free(self.output) 476 451 … … 513 488 [2, 3, 4, 5, 0, 1] 514 489 515 490 """ 516 cdef int **part, i, j 491 cdef int n = self.degree 492 cdef PartitionStack *part 517 493 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) 520 497 if part is NULL: 521 498 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 MemoryError529 for j from 0 <= j < len(partition[i]):530 part[i][j] = partition[i][j]531 part[i][len(partition[i])] = -1532 part[len(partition)] = NULL533 499 self.first_time = 1 534 500 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) 536 502 537 for i from 0 <= i < len(partition): 538 sage_free(part[i]) 539 sage_free(part) 503 PS_dealloc(part) 540 504 541 505 def automorphism_group(self): 542 506 """ … … 559 523 560 524 """ 561 525 cdef int i, j 562 cdef object generators, base526 cdef list generators, base 563 527 cdef Integer order 564 528 if self.output is NULL: 565 529 self.run() … … 567 531 for i from 0 <= i < self.output.num_gens: 568 532 generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree]) 569 533 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] 572 536 return generators, order, base 573 537 574 538 def canonical_relabeling(self): … … 602 566 [2, 3, 0, 1, 4, 5] 603 567 604 568 """ 605 cdef int **part, i, j569 cdef int i, n = self.degree 606 570 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)) 610 574 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) 612 576 if ordering is not NULL: sage_free(ordering) 613 577 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: 626 579 ordering[i] = i 627 580 self.first_time = 1 628 581 other.first_time = 1 629 582 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) 631 584 632 for i from 0 <= i < len(partition): 633 sage_free(part[i]) 634 sage_free(part) 585 PS_dealloc(part) 635 586 sage_free(ordering) 636 587 637 588 if output is NULL: … … 646 597 bitset_copy(word, &NBCS.words[i]) 647 598 return 0 648 599 649 cdef int refine_by_bip_degree(PartitionStack *col_ps, objectS, int *cells_to_refine_by, int ctrb_len):600 cdef int refine_by_bip_degree(PartitionStack *col_ps, void *S, int *cells_to_refine_by, int ctrb_len): 650 601 """ 651 602 Refines the input partition by checking degrees of vertices to the given 652 603 cells in the associated bipartite graph (vertices split into columns and … … 711 662 # now, i points to the next cell (before refinement) 712 663 if necessary_to_split_cell: 713 664 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) 715 666 invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps) 716 667 invariant += first_largest_subcell 717 668 against_index = current_cell_against … … 746 697 # now, i points to the next cell (before refinement) 747 698 if necessary_to_split_cell: 748 699 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) 750 701 invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps) 751 702 invariant += first_largest_subcell 752 703 against_index = current_cell_against … … 770 721 current_cell_against += 1 771 722 return invariant 772 723 773 cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S1, objectS2):724 cdef int compare_linear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2): 774 725 """ 775 726 Compare gamma_1(S1) and gamma_2(S2). 776 727 … … 843 794 return <int>bitset_check(&basis_2[i], gamma_2[cur_col]) - <int>bitset_check(&basis_1[i], gamma_1[cur_col]) 844 795 return 0 845 796 846 cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S1, objectS2):797 cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2): 847 798 """ 848 799 Compare gamma_1(S1) and gamma_2(S2). 849 800 … … 923 874 924 875 return 0 925 876 926 cdef bint all_children_are_equivalent(PartitionStack *col_ps, objectS):877 cdef bint all_children_are_equivalent(PartitionStack *col_ps, void *S): 927 878 """ 928 879 Returns True if any refinement of the current partition results in the same 929 880 structure. … … 1022 973 bitset_free(word) 1023 974 return degree 1024 975 1025 cdef inline int sort_by_function (PartitionStack *PS, int start, int *degrees, int *counts, int *output, int count_max):976 cdef inline int sort_by_function_codes(PartitionStack *PS, int start, int *degrees, int *counts, int *output, int count_max): 1026 977 """ 1027 978 A simple counting sort, given the degrees of vertices to a certain cell. 1028 979 -
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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # 5 5 # Distributed under the terms of the GNU General Public License (GPL) 6 6 # http://www.gnu.org/licenses/ … … 14 14 from sage.graphs.base.sparse_graph cimport SparseGraph 15 15 from sage.graphs.base.dense_graph cimport DenseGraph 16 16 from sage.rings.integer cimport Integer 17 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab _return17 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab 18 18 from double_coset cimport double_coset 19 19 20 20 cdef class GraphStruct: … … 23 23 cdef bint use_indicator 24 24 cdef int *scratch # length 3n+1 25 25 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)26 cdef int refine_by_degree(PartitionStack *, void *, int *, int) 27 cdef int compare_graphs(int *, int *, void *, void *) 28 cdef bint all_children_are_equivalent(PartitionStack *, void *) 29 29 cdef inline int degree(PartitionStack *, CGraph, int, int, bint) 30 cdef inline int sort_by_function(PartitionStack *, int, int *)31 30 32 31 33 32 -
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 12 12 """ 13 13 14 14 #***************************************************************************** 15 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>15 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 16 16 # 17 17 # Distributed under the terms of the GNU General Public License (GPL) 18 18 # http://www.gnu.org/licenses/ … … 49 49 {0: 1, 1: 2, 2: 0} 50 50 51 51 """ 52 cdef int **part52 cdef PartitionStack *part 53 53 cdef int *output, *ordering 54 54 cdef CGraph G 55 55 cdef GraphStruct GS1 = GraphStruct() 56 56 cdef GraphStruct GS2 = GraphStruct() 57 57 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 59 60 60 61 from sage.graphs.all import Graph, DiGraph 61 62 from sage.graphs.generic_graph import GenericGraph … … 82 83 if first: 83 84 partition = [[to[v] for v in cell] for cell in partn] 84 85 else: 86 if first: 87 partition = partn 85 88 to = range(n) 86 89 frm = to 87 90 if sparse: … … 104 107 to = {} 105 108 for a in G.verts(): to[a]=a 106 109 frm = to 110 if first: 111 partition = partn 107 112 else: 108 113 raise TypeError("G must be a Sage graph.") 109 114 if first: frm1=frm;to1=to … … 115 120 116 121 if n == 0: 117 122 return {} 118 119 part = <int **> sage_malloc((len(partn)+1) * sizeof(int *))123 124 part = PS_from_list(partition) 120 125 ordering = <int *> sage_malloc(n * sizeof(int)) 121 126 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) 123 128 if ordering is not NULL: sage_free(ordering) 124 129 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 MemoryError132 for j from 0 <= j < len(partn[i]):133 part[i][j] = to1[partn[i][j]]134 part[i][len(partn[i])] = -1135 part[len(partn)] = NULL136 130 for i from 0 <= i < n: 137 131 ordering[i] = to2[ordering2[i]] 138 132 … … 141 135 if GS1.scratch is NULL or GS2.scratch is NULL: 142 136 if GS1.scratch is not NULL: sage_free(GS1.scratch) 143 137 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) 147 140 raise MemoryError 148 141 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) 150 143 151 for i from 0 <= i < len(partn): 152 sage_free(part[i]) 153 sage_free(part) 144 PS_dealloc(part) 154 145 sage_free(ordering) 155 146 sage_free(GS1.scratch) 156 147 sage_free(GS2.scratch) … … 350 341 cdef CGraph G 351 342 cdef int i, j, n 352 343 cdef Integer I 353 cdef aut_gp_and_can_lab _return*output354 cdef int **part344 cdef aut_gp_and_can_lab *output 345 cdef PartitionStack *part 355 346 from sage.graphs.all import Graph, DiGraph 356 347 from sage.graphs.generic_graph import GenericGraph 357 348 from copy import copy … … 387 378 else: 388 379 raise TypeError("G must be a Sage graph.") 389 380 390 part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))391 if part is NULL:392 raise MemoryError393 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 MemoryError400 for j from 0 <= j < len(partition[i]):401 part[i][j] = partition[i][j]402 part[i][len(partition[i])] = -1403 part[len(partition)] = NULL404 405 381 cdef GraphStruct GS = GraphStruct() 406 382 GS.G = G 407 383 GS.directed = 1 if dig else 0 408 384 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 409 410 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) 414 415 raise MemoryError 415 416 416 417 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) 418 419 sage_free( GS.scratch ) 419 420 # prepare output 420 421 list_of_gens = [] … … 445 446 dd[frm[i]] = output.relabeling[i] 446 447 return_tuple.append(dd) 447 448 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]) 449 450 if order: 450 451 I = Integer() 451 mpz_set(I.value, output.order)452 SC_order(output.group, 0, I.value) 452 453 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) 457 455 sage_free(output.generators) 458 if base: 459 sage_free(output.base) 456 SC_dealloc(output.group) 460 457 if lab_new: 461 458 sage_free(output.relabeling) 462 459 sage_free(output) … … 465 462 else: 466 463 return tuple(return_tuple) 467 464 468 cdef int refine_by_degree(PartitionStack *PS, objectS, int *cells_to_refine_by, int ctrb_len):465 cdef int refine_by_degree(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len): 469 466 r""" 470 467 Refines the input partition by checking degrees of vertices to the given 471 468 cells. … … 584 581 else: 585 582 return 0 586 583 587 cdef int compare_graphs(int *gamma_1, int *gamma_2, object S1, objectS2):584 cdef int compare_graphs(int *gamma_1, int *gamma_2, void *S1, void *S2): 588 585 r""" 589 586 Compare gamma_1(S1) and gamma_2(S2). 590 587 … … 611 608 return -1 612 609 return 0 613 610 614 cdef bint all_children_are_equivalent(PartitionStack *PS, objectS):611 cdef bint all_children_are_equivalent(PartitionStack *PS, void *S): 615 612 """ 616 613 Return True if every refinement of the current partition results in the 617 614 same structure. … … 681 678 break 682 679 return num_arcs 683 680 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 checked690 start -- beginning index of the cell to be sorted691 degrees -- the values to be sorted by692 693 """694 cdef int n = PS.degree695 cdef int i, j, max, max_location696 cdef int *counts = degrees + n, *output = degrees + 2*n + 1697 for i from 0 <= i <= n:698 counts[i] = 0699 i = 0700 while PS.levels[i+start] > PS.depth:701 counts[degrees[i]] += 1702 i += 1703 counts[degrees[i]] += 1704 # i+start is the right endpoint of the cell now705 max = counts[0]706 max_location = 0707 for j from 0 < j <= n:708 if counts[j] > max:709 max = counts[j]710 max_location = j711 counts[j] += counts[j - 1]712 for j from i >= j >= 0:713 counts[degrees[j]] -= 1714 output[counts[degrees[j]]] = PS.entries[start+j]715 max_location = counts[max_location]+start716 for j from 0 <= j <= i:717 PS.entries[start+j] = output[j]718 j = 1719 while j <= n and counts[j] <= i:720 if counts[j] > 0:721 PS.levels[start + counts[j] - 1] = PS.depth722 PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1)723 j += 1724 return max_location725 726 681 def all_labeled_graphs(n): 727 682 """ 728 683 Returns all labeled graphs on n vertices {0,1,...,n-1}. Used in … … 773 728 return Glist 774 729 775 730 776 def random_tests(num= 20, n_max=50, perms_per_graph=8):731 def random_tests(num=10, n_max=60, perms_per_graph=5): 777 732 """ 778 733 Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma 779 734 and random graphs G, and that isomorphic returns an isomorphism. … … 796 751 TESTS:: 797 752 798 753 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. 801 756 """ 802 757 from sage.misc.prandom import random, randint 803 758 from sage.graphs.graph_generators import GraphGenerators … … 982 937 GS.G = G 983 938 GS.scratch = <int *> sage_malloc((3*n+1) * sizeof(int)) 984 939 if GS.scratch is NULL: 985 PS_ clear(nu)940 PS_dealloc(nu) 986 941 raise MemoryError 987 942 GS.directed = directed 988 943 GS.use_indicator = 0 … … 991 946 cdef int num_cells = len(partition) 992 947 cdef int *alpha = <int *>sage_malloc(n * sizeof(int)) 993 948 if alpha is NULL: 994 PS_ clear(nu)949 PS_dealloc(nu) 995 950 sage_free(GS.scratch) 996 951 raise MemoryError 997 952 j = 0 … … 1000 955 j += len(partition[i]) 1001 956 1002 957 # refine, and get the result 1003 refine_by_degree(nu, GS, alpha, num_cells)958 refine_by_degree(nu, <void *>GS, alpha, num_cells) 1004 959 1005 960 eq_part = [] 1006 961 cell = [] … … 1010 965 eq_part.append(cell) 1011 966 cell = [] 1012 967 1013 PS_ clear(nu)968 PS_dealloc(nu) 1014 969 sage_free(GS.scratch) 1015 970 sage_free(alpha) 1016 971 -
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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # Copyright (C) 2009 Nicolas Borie <nicolas.borie@math.u-psud.fr> 5 5 # 6 6 # Distributed under the terms of the GNU General Public License (GPL) … … 17 17 from double_coset cimport double_coset 18 18 19 19 # 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 20 cdef int refine_list(PartitionStack *, void *, int *, int) 21 cdef int compare_lists(int *, int *, void *, void *) 22 cdef 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 7 7 """ 8 8 9 9 #***************************************************************************** 10 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>10 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 11 11 # Copyright (C) 2009 Nicolas Borie <nicolas.borie@math.u-psud.fr> 12 12 # 13 13 # Distributed under the terms of the GNU General Public License (GPL) … … 28 28 [1, 2, 0] 29 29 30 30 """ 31 cdef int **part, i, j 31 cdef int i, n = len(self) 32 cdef PartitionStack *part 32 33 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) 35 35 ordering = <int *> sage_malloc((len(self)) * sizeof(int)) 36 36 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) 38 38 if ordering is not NULL: sage_free(ordering) 39 39 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 MemoryError47 for j from 0 <= j < len(partition[i]):48 part[i][j] = partition[i][j]49 part[i][len(partition[i])] = -150 part[len(partition)] = NULL51 40 for i from 0 <= i < (len(self)): 52 41 ordering[i] = i 53 42 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) 55 44 56 for i from 0 <= i < len(partition): 57 sage_free(part[i]) 58 sage_free(part) 45 PS_dealloc(part) 59 46 sage_free(ordering) 60 47 61 48 if output is NULL: … … 65 52 sage_free(output) 66 53 return output_py 67 54 68 cdef bint all_list_children_are_equivalent(PartitionStack *PS, objectS):55 cdef bint all_list_children_are_equivalent(PartitionStack *PS, void *S): 69 56 return 0 70 57 71 cdef int refine_list(PartitionStack *PS, objectS, int *cells_to_refine_by, int ctrb_len):58 cdef int refine_list(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len): 72 59 return 0 73 60 74 cdef int compare_lists(int *gamma_1, int *gamma_2, object S1, objectS2):61 cdef int compare_lists(int *gamma_1, int *gamma_2, void *S1, void *S2): 75 62 r""" 76 63 Compare two lists according to the lexicographic order. 77 64 """ -
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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # 5 5 # Distributed under the terms of the GNU General Public License (GPL) 6 6 # http://www.gnu.org/licenses/ … … 11 11 include 'data_structures_pxd.pxi' # includes bitsets 12 12 13 13 from sage.rings.integer cimport Integer 14 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab _return14 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab 15 15 from refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree 16 16 from refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent 17 17 from double_coset cimport double_coset … … 24 24 cdef list symbols 25 25 cdef int nsymbols 26 26 cdef PartitionStack *temp_col_ps 27 cdef aut_gp_and_can_lab _return*output27 cdef aut_gp_and_can_lab *output 28 28 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)29 cdef int refine_matrix(PartitionStack *, void *, int *, int) 30 cdef int compare_matrices(int *, int *, void *, void *) 31 cdef bint all_matrix_children_are_equivalent(PartitionStack *, void *) 32 32 -
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 16 16 """ 17 17 18 18 #***************************************************************************** 19 # Copyright (C) 2006 - 20 08Robert L. Miller <rlmillster@gmail.com>19 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 20 20 # 21 21 # Distributed under the terms of the GNU General Public License (GPL) 22 22 # http://www.gnu.org/licenses/ … … 81 81 def __dealloc__(self): 82 82 PS_dealloc(self.temp_col_ps) 83 83 if self.output is not NULL: 84 mpz_clear(self.output.order)85 84 sage_free(self.output.generators) 86 sage_free(self.output.base)85 SC_dealloc(self.output.group) 87 86 sage_free(self.output.relabeling) 88 87 sage_free(self.output) 89 88 … … 148 147 True 149 148 150 149 """ 151 cdef int **part, i, j 150 cdef int i, n = self.degree 151 cdef PartitionStack *part 152 152 cdef NonlinearBinaryCodeStruct S_temp 153 153 for i from 0 <= i < self.nsymbols: 154 154 S_temp = <NonlinearBinaryCodeStruct> self.symbol_structs[i] 155 155 S_temp.first_time = 1 156 156 157 157 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) 160 161 if part is NULL: 161 162 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 MemoryError169 for j from 0 <= j < len(partition[i]):170 part[i][j] = partition[i][j]171 part[i][len(partition[i])] = -1172 part[len(partition)] = NULL173 163 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) 175 165 176 for i from 0 <= i < len(partition): 177 sage_free(part[i]) 178 sage_free(part) 166 PS_dealloc(part) 179 167 180 168 181 169 def automorphism_group(self): … … 193 181 194 182 """ 195 183 cdef int i, j 196 cdef object generators, base184 cdef list generators, base 197 185 cdef Integer order 198 186 if self.output is NULL: 199 187 self.run() … … 201 189 for i from 0 <= i < self.output.num_gens: 202 190 generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree]) 203 191 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] 206 194 return generators, order, base 207 195 208 196 def canonical_relabeling(self): … … 234 222 [0, 2, 4, 1, 3, 5] 235 223 236 224 """ 237 238 cdef int **part, i, j 225 cdef int i, j, n = self.degree 239 226 cdef int *output, *ordering 227 cdef PartitionStack *part 240 228 cdef NonlinearBinaryCodeStruct S_temp 241 229 for i from 0 <= i < self.nsymbols: 242 230 S_temp = self.symbol_structs[i] 243 231 S_temp.first_time = 1 244 232 S_temp = other.symbol_structs[i] 245 233 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) 248 235 ordering = <int *> sage_malloc(self.degree * sizeof(int)) 249 236 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) 251 238 if ordering is not NULL: sage_free(ordering) 252 239 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 MemoryError260 for j from 0 <= j < len(partition[i]):261 part[i][j] = partition[i][j]262 part[i][len(partition[i])] = -1263 part[len(partition)] = NULL264 240 for i from 0 <= i < self.degree: 265 241 ordering[i] = i 266 242 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) 268 244 269 for i from 0 <= i < len(partition): 270 sage_free(part[i]) 271 sage_free(part) 245 PS_dealloc(part) 272 246 sage_free(ordering) 273 247 274 248 if output is NULL: … … 278 252 sage_free(output) 279 253 return output_py 280 254 281 cdef int refine_matrix(PartitionStack *PS, objectS, int *cells_to_refine_by, int ctrb_len):255 cdef int refine_matrix(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len): 282 256 cdef MatrixStruct M = <MatrixStruct> S 283 257 cdef int i, temp_inv, invariant = 1 284 258 cdef bint changed = 1 285 259 while changed: 286 260 PS_copy_from_to(PS, M.temp_col_ps) 287 261 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) 289 263 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: 292 265 changed = 0 293 266 return invariant 294 267 295 cdef int compare_matrices(int *gamma_1, int *gamma_2, object S1, objectS2):268 cdef int compare_matrices(int *gamma_1, int *gamma_2, void *S1, void *S2): 296 269 cdef MatrixStruct MS1 = <MatrixStruct> S1 297 270 cdef MatrixStruct MS2 = <MatrixStruct> S2 298 271 M1 = MS1.matrix … … 305 278 MM2.set_column(i, M2.column(gamma_2[i])) 306 279 return cmp(sorted(MM1.rows()), sorted(MM2.rows())) 307 280 308 cdef bint all_matrix_children_are_equivalent(PartitionStack *PS, objectS):281 cdef bint all_matrix_children_are_equivalent(PartitionStack *PS, void *S): 309 282 return 0 310 283 311 def random_tests(n=1 5, nrows_max=50, ncols_max=50, nsymbols_max=20, perms_per_matrix=10, density_range=(.1,.9)):284 def random_tests(n=10, nrows_max=50, ncols_max=50, nsymbols_max=10, perms_per_matrix=5, density_range=(.1,.9)): 312 285 """ 313 286 Tests to make sure that C(gamma(M)) == C(M) for random permutations gamma 314 287 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 1 1 2 2 #***************************************************************************** 3 # Copyright (C) 2006 - 20 09Robert L. Miller <rlmillster@gmail.com>3 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 4 4 # 5 5 # Distributed under the terms of the GNU General Public License (GPL) 6 6 # http://www.gnu.org/licenses/ … … 13 13 14 14 15 15 from sage.rings.integer cimport Integer 16 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab _return16 from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab 17 17 from double_coset cimport double_coset 18 18 19 19 -
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 20 20 """ 21 21 22 22 #***************************************************************************** 23 # Copyright (C) 2006 - 20 09Robert L. Miller <rlmillster@gmail.com>23 # Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 24 24 # 25 25 # Distributed under the terms of the GNU General Public License (GPL) 26 26 # http://www.gnu.org/licenses/ … … 376 376 self.rari_fn = rari_fn 377 377 self.cs_fn = cs_fn 378 378 379 cdef bint all_children_are_equivalent_python(PartitionStack *PS, objectS):379 cdef bint all_children_are_equivalent_python(PartitionStack *PS, void *S): 380 380 """ 381 381 Python conversion of all_children_are_equivalent function. 382 382 """ 383 383 cdef PythonPartitionStack Py_PS = PythonPartitionStack(PS.degree) 384 cdef object S_obj = <object> S 384 385 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) 386 387 387 cdef int refine_and_return_invariant_python(PartitionStack *PS, objectS, int *cells_to_refine_by, int ctrb_len):388 cdef int refine_and_return_invariant_python(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len): 388 389 """ 389 390 Python conversion of refine_and_return_invariant function. 390 391 """ 391 392 cdef PythonPartitionStack Py_PS = PythonPartitionStack(PS.degree) 393 cdef object S_obj = <object> S 392 394 PS_copy_from_to(PS, Py_PS.c_ps) 393 395 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) 396 398 397 cdef int compare_structures_python(int *gamma_1, int *gamma_2, object S1, objectS2):399 cdef int compare_structures_python(int *gamma_1, int *gamma_2, void *S1, void *S2): 398 400 """ 399 401 Python conversion of compare_structures function. 400 402 """ 401 403 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) 405 408 406 409 def aut_gp_and_can_lab_python(S, partition, n, 407 410 all_children_are_equivalent, … … 455 458 456 459 """ 457 460 obj_wrapper = PythonObjectWrapper(S, all_children_are_equivalent, refine_and_return_invariant, compare_structures, n) 458 cdef aut_gp_and_can_lab _return*output461 cdef aut_gp_and_can_lab *output 459 462 cdef PythonPartitionStack Py_PS = PythonPartitionStack(n) 460 463 cdef int i, j 461 464 cdef Integer I 462 465 463 cdef int **part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))466 cdef PartitionStack *part = PS_from_list(partition) 464 467 if part is NULL: 465 468 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 MemoryError473 for j from 0 <= j < len(partition[i]):474 part[i][j] = partition[i][j]475 part[i][len(partition[i])] = -1476 part[len(partition)] = NULL477 469 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, 479 471 &all_children_are_equivalent_python, 480 472 &refine_and_return_invariant_python, 481 473 &compare_structures_python, 482 canonical_label, base, order)474 canonical_label, NULL) 483 475 484 476 list_of_gens = [] 485 477 for i from 0 <= i < output.num_gens: … … 488 480 if canonical_label: 489 481 return_tuple.append([output.relabeling[i] for i from 0 <= i < n]) 490 482 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]) 492 484 if order: 493 485 I = Integer() 494 mpz_set(I.value, output.order)486 SC_order(output.group, 0, I.value) 495 487 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) 500 489 sage_free(output.generators) 501 if base: 502 sage_free(output.base) 490 SC_dealloc(output.group) 503 491 if canonical_label: 504 492 sage_free(output.relabeling) 505 493 sage_free(output) … … 559 547 obj_wrapper1 = PythonObjectWrapper(S1, all_children_are_equivalent, refine_and_return_invariant, compare_structures, n) 560 548 obj_wrapper2 = PythonObjectWrapper(S2, all_children_are_equivalent, refine_and_return_invariant, compare_structures, n) 561 549 562 cdef int **part = <int **> sage_malloc((len(partition1)+1) * sizeof(int *))550 cdef PartitionStack *part = PS_from_list(partition1) 563 551 cdef int *ordering = <int *> sage_malloc(n * sizeof(int)) 564 552 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) 569 555 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 MemoryError577 for j from 0 <= j < len(partition1[i]):578 part[i][j] = partition1[i][j]579 part[i][len(partition1[i])] = -1580 part[len(partition1)] = NULL581 556 for i from 0 <= i < n: 582 557 ordering[i] = ordering2[i] 583 558 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, 585 561 &all_children_are_equivalent_python, 586 562 &refine_and_return_invariant_python, 587 &compare_structures_python )563 &compare_structures_python, NULL) 588 564 589 for i from 0 <= i < len(partition1): 590 sage_free(part[i]) 591 sage_free(part) 565 PS_dealloc(part) 592 566 sage_free(ordering) 593 567 594 568 if output is NULL: -
sage/sets/disjoint_set.pyx
diff -r 872cd84f620c -r 8ddd7727b639 sage/sets/disjoint_set.pyx
a b 60 60 include '../groups/perm_gps/partn_ref/data_structures_pyx.pxi' 61 61 62 62 import itertools 63 from sage.rings. allimport Integer63 from sage.rings.integer import Integer 64 64 from sage.structure.sage_object cimport SageObject 65 65 66 66 def DisjointSet(arg):
