CoCalc -- refinement

GitHub Repository: sagemath/sagesmc
Path: blob/master/src/sage/groups/perm_gps/partn_ref/refinement_sets.pyx
⁸⁸¹⁷ views
1
r"""
2
Partition backtrack functions for sets
3

4
DOCTEST:
5
    sage: import sage.groups.perm_gps.partn_ref.refinement_sets
6

7
REFERENCE:
8

9
    [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
10
        Vol. 30 (1981), pp. 45-87.
11

12
    [2] Leon, Jeffrey. Permutation Group Algorithms Based on Partitions, I:
13
        Theory and Algorithms. J. Symbolic Computation, Vol. 12 (1991), pp.
14
        533-583.
15

16
"""
17

18
#*****************************************************************************
19
#      Copyright (C) 2006 - 2011 Robert L. Miller <[email protected]>
20
#
21
# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
22
#                         http://www.gnu.org/licenses/
23
#*****************************************************************************
24

25
include 'data_structures_pyx.pxi' # includes bitsets
26

27
def set_stab_py(generators, sett, relab=False):
28
    r"""
29
    Compute the setwise stabilizer of a subset of [0..n-1] in a subgroup of S_n.
30

31
    EXAMPLES::
32

33
        sage: from sage.groups.perm_gps.partn_ref.refinement_sets import set_stab_py
34

35
    Degree 4 examples
36

37
    A four-cycle::
38

39
        sage: set_stab_py([[1,2,3,0]], [0])
40
        []
41
        sage: set_stab_py([[1,2,3,0]], [0,1])
42
        []
43
        sage: set_stab_py([[1,2,3,0]], [0,1,2])
44
        []
45
        sage: set_stab_py([[1,2,3,0]], [0,1,2,3])
46
        [[1, 2, 3, 0]]
47
        sage: set_stab_py([[1,2,3,0]], [0,2])
48
        [[2, 3, 0, 1]]
49

50
    Symmetric group::
51

52
        sage: set_stab_py([[1,0,2,3],[1,2,3,0]], [0])
53
        [[0, 1, 3, 2], [0, 2, 1, 3]]
54
        sage: set_stab_py([[1,0,2,3],[1,2,3,0]], [0,1])
55
        [[1, 0, 2, 3], [0, 1, 3, 2]]
56
        sage: set_stab_py([[1,0,2,3],[1,2,3,0]], [0,1,2,3])
57
        [[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]]
58
        sage: set_stab_py([[1,0,2,3],[1,2,3,0]], [0,3])
59
        [[3, 1, 2, 0], [0, 2, 1, 3]]
60

61
    Klein 4-group::
62

63
        sage: set_stab_py([[1,0,2,3],[0,1,3,2]], [0])
64
        [[0, 1, 3, 2]]
65
        sage: set_stab_py([[1,0,2,3],[0,1,3,2]], [0,1])
66
        [[0, 1, 3, 2], [1, 0, 2, 3]]
67
        sage: set_stab_py([[1,0,2,3],[0,1,3,2]], [0,2])
68
        []
69

70
    Dihedral group::
71

72
        sage: set_stab_py([[1,2,3,0],[0,3,2,1]], [0])
73
        [[0, 3, 2, 1]]
74
        sage: set_stab_py([[1,2,3,0],[0,3,2,1]], [0,1])
75
        [[1, 0, 3, 2]]
76
        sage: set_stab_py([[1,2,3,0],[0,3,2,1]], [0,2])
77
        [[2, 1, 0, 3], [0, 3, 2, 1]]
78
        sage: set_stab_py([[1,2,3,0],[0,3,2,1]], [1])
79
        [[2, 1, 0, 3]]
80
        sage: set_stab_py([[1,2,3,0],[0,3,2,1]], [1,2,3])
81
        [[0, 3, 2, 1]]
82

83
    Alternating group::
84

85
        sage: set_stab_py([[1,2,0,3],[0,2,3,1]], [0])
86
        [[0, 2, 3, 1]]
87
        sage: set_stab_py([[1,2,0,3],[0,2,3,1]], [1])
88
        [[2, 1, 3, 0]]
89
        sage: set_stab_py([[1,2,0,3],[0,2,3,1]], [2])
90
        [[1, 3, 2, 0]]
91
        sage: set_stab_py([[1,2,0,3],[0,2,3,1]], [3])
92
        [[1, 2, 0, 3]]
93
        sage: set_stab_py([[1,2,0,3],[0,2,3,1]], [0,1])
94
        [[1, 0, 3, 2]]
95
        sage: set_stab_py([[1,2,0,3],[0,2,3,1]], [0,2])
96
        [[2, 3, 0, 1]]
97
        sage: set_stab_py([[1,2,0,3],[0,2,3,1]], [0,3])
98
        [[3, 2, 1, 0]]
99

100
    Larger degree examples
101

102
    Dihedral group of degree 5::
103

104
        sage: set_stab_py([[1,2,3,4,0],[0,4,3,2,1]], [])
105
        [[0, 4, 3, 2, 1], [1, 0, 4, 3, 2]]
106
        sage: set_stab_py([[1,2,3,4,0],[0,4,3,2,1]], [0])
107
        [[0, 4, 3, 2, 1]]
108
        sage: set_stab_py([[1,2,3,4,0],[0,4,3,2,1]], [0,2])
109
        [[2, 1, 0, 4, 3]]
110

111
    Dihedral group of degree 6::
112

113
        sage: set_stab_py([[1,2,3,4,5,0],[0,5,4,3,2,1]], [])
114
        [[0, 5, 4, 3, 2, 1], [1, 0, 5, 4, 3, 2]]
115
        sage: set_stab_py([[1,2,3,4,5,0],[0,5,4,3,2,1]], [0])
116
        [[0, 5, 4, 3, 2, 1]]
117
        sage: set_stab_py([[1,2,3,4,5,0],[0,5,4,3,2,1]], [0,1])
118
        [[1, 0, 5, 4, 3, 2]]
119
        sage: set_stab_py([[1,2,3,4,5,0],[0,5,4,3,2,1]], [0,2])
120
        [[2, 1, 0, 5, 4, 3]]
121
        sage: set_stab_py([[1,2,3,4,5,0],[0,5,4,3,2,1]], [0,3])
122
        [[0, 5, 4, 3, 2, 1], [3, 2, 1, 0, 5, 4]]
123
        sage: set_stab_py([[1,2,3,4,5,0],[0,5,4,3,2,1]], [0,2,4])
124
        [[2, 1, 0, 5, 4, 3], [4, 3, 2, 1, 0, 5]]
125

126
    Canonical labels::
127

128
        sage: set_stab_py([[0,2,1,4,3,5,8,7,6],[8,7,6,3,5,4,2,1,0]], [0,1,3,5,6], True)
129
        ([], [7, 8, 6, 3, 4, 5, 2, 0, 1])
130
        sage: set_stab_py([[0,2,1,4,3,5,8,7,6],[8,7,6,3,5,4,2,1,0]], [0,3,5,6,8], True)
131
        ([], [2, 1, 0, 5, 4, 3, 7, 6, 8])
132

133
    """
134
    if len(generators) == 0:
135
        return []
136
    cdef int i, j, n = len(generators[0]), n_gens = len(generators)
137
    cdef StabilizerChain *supergroup = SC_new(n)
138
    cdef aut_gp_and_can_lab *stabilizer
139
    cdef int *gens = <int *> sage_malloc(n*n_gens * sizeof(int))
140
    cdef subset *subset_sett = <subset *> sage_malloc(sizeof(subset))
141
    if gens is NULL or supergroup is NULL or subset_sett is NULL:
142
        SC_dealloc(supergroup)
143
        sage_free(gens)
144
        sage_free(subset_sett)
145
        raise MemoryError
146
    bitset_init(&subset_sett.bits, n)
147
    subset_sett.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
148
    for i from 0 <= i < len(generators):
149
        for j from 0 <= j < n:
150
            gens[n*i + j] = generators[i][j]
151
    if SC_insert(supergroup, 0, gens, n_gens):
152
        SC_dealloc(supergroup)
153
        sage_free(gens)
154
        sage_free(subset_sett)
155
        raise MemoryError
156
    sage_free(gens)
157
    bitset_clear(&subset_sett.bits)
158
    for i in sett:
159
        bitset_add(&subset_sett.bits, i)
160
    stabilizer = set_stab(supergroup, subset_sett, relab)
161
    SC_dealloc(supergroup)
162
    bitset_free(&subset_sett.bits)
163
    sage_free(subset_sett.scratch)
164
    sage_free(subset_sett)
165
    if stabilizer is NULL:
166
        raise MemoryError
167
    stab_gens = []
168
    for i from 0 <= i < stabilizer.num_gens:
169
        stab_gens.append([stabilizer.generators[i*n+j] for j from 0 <= j < n])
170
    if relab:
171
        relabeling = [stabilizer.relabeling[j] for j from 0 <= j < n]
172
    deallocate_agcl_output(stabilizer)
173
    if relab:
174
        return stab_gens, relabeling
175
    return stab_gens
176

177
cdef aut_gp_and_can_lab *set_stab(StabilizerChain *supergroup, subset *sett, bint relab):
178
    r"""
179
    Computes the set stabilizer of ``sett`` within ``supergroup``. (Note that
180
    ``set`` is a reserved Python keyword.) If ``relab`` is specified then
181
    computes the canonical label of the set under the action of the group.
182
    """
183
    cdef aut_gp_and_can_lab *output
184
    cdef int n = supergroup.degree
185
    cdef PartitionStack *part = PS_new(n, 1)
186
    if part is NULL:
187
        return NULL
188
    output = get_aut_gp_and_can_lab(<void *> sett, part, n,
189
        &all_set_children_are_equivalent, &refine_set, &compare_sets, relab,
190
        supergroup, NULL, NULL)
191
    PS_dealloc(part)
192
    if output is NULL:
193
        return NULL
194
    return output
195

196
def sets_isom_py(generators, set1, set2):
197
    r"""
198
    Computes whether ``set1`` and ``set2`` are isomorphic under the action of
199
    the group generated by the generators given in list form.
200

201
    EXAMPLES::
202

203
        sage: from sage.groups.perm_gps.partn_ref.refinement_sets import sets_isom_py
204

205
    Degree 3 examples
206

207
    Trivial group::
208

209
        sage: sets_isom_py([], [0,1,2], [0,1])
210
        False
211
        sage: sets_isom_py([], [0,1,2], [0,2,1])
212
        [0, 1, 2]
213
        sage: sets_isom_py([[0,1,2]], [0,1,2], [0,2,1])
214
        [0, 1, 2]
215
        sage: sets_isom_py([[0,1,2]], [0], [0])
216
        [0, 1, 2]
217
        sage: sets_isom_py([[0,1,2]], [0], [1])
218
        False
219
        sage: sets_isom_py([[0,1,2]], [0], [2])
220
        False
221
        sage: sets_isom_py([[0,1,2]], [0,1], [1,0])
222
        [0, 1, 2]
223

224
    Three-cycle::
225

226
        sage: sets_isom_py([[1,2,0]], [0], [1])
227
        [1, 2, 0]
228
        sage: sets_isom_py([[1,2,0]], [0], [2])
229
        [2, 0, 1]
230
        sage: sets_isom_py([[1,2,0]], [0], [0])
231
        [0, 1, 2]
232
        sage: sets_isom_py([[1,2,0]], [0,1], [0])
233
        False
234
        sage: sets_isom_py([[1,2,0]], [0,1], [1])
235
        False
236
        sage: sets_isom_py([[1,2,0]], [0,1], [2])
237
        False
238
        sage: sets_isom_py([[1,2,0]], [0,1], [0,2])
239
        [2, 0, 1]
240
        sage: sets_isom_py([[1,2,0]], [0,1], [1,2])
241
        [1, 2, 0]
242
        sage: sets_isom_py([[1,2,0]], [0,1], [1,0])
243
        [0, 1, 2]
244
        sage: sets_isom_py([[1,2,0]], [0,2,1], [2,1,0])
245
        [0, 1, 2]
246

247
    Transposition::
248

249
        sage: sets_isom_py([[1,0,2]], [0], [])
250
        False
251
        sage: sets_isom_py([[1,0,2]], [0], [1,2])
252
        False
253
        sage: sets_isom_py([[1,0,2]], [0], [1])
254
        [1, 0, 2]
255
        sage: sets_isom_py([[1,0,2]], [0], [0])
256
        [0, 1, 2]
257
        sage: sets_isom_py([[1,0,2]], [0,1], [2])
258
        False
259
        sage: sets_isom_py([[1,0,2]], [0,2], [1,2])
260
        [1, 0, 2]
261
        sage: sets_isom_py([[1,0,2]], [0], [2])
262
        False
263
        sage: sets_isom_py([[1,0,2]], [0,1], [1,2])
264
        False
265

266
    Symmetric group S_3::
267

268
        sage: sets_isom_py([[1,0,2],[1,2,0]], [], [])
269
        [0, 1, 2]
270
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0], [])
271
        False
272
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0], [0])
273
        [0, 1, 2]
274
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0], [1])
275
        [1, 0, 2]
276
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0], [2])
277
        [2, 0, 1]
278
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0,2], [2])
279
        False
280
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0,2], [1,2])
281
        [1, 0, 2]
282
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0,2], [0,1])
283
        [0, 2, 1]
284
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0,2], [0,2])
285
        [0, 1, 2]
286
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0,2], [0,1,2])
287
        False
288
        sage: sets_isom_py([[1,0,2],[1,2,0]], [0,2,1], [0,1,2])
289
        [0, 1, 2]
290

291
    Degree 4 examples
292

293
    Trivial group::
294

295
        sage: sets_isom_py([[0,1,2,3]], [], [])
296
        [0, 1, 2, 3]
297
        sage: sets_isom_py([[0,1,2,3]], [0], [])
298
        False
299
        sage: sets_isom_py([[0,1,2,3]], [0], [1])
300
        False
301
        sage: sets_isom_py([[0,1,2,3]], [0], [2])
302
        False
303
        sage: sets_isom_py([[0,1,2,3]], [0], [3])
304
        False
305
        sage: sets_isom_py([[0,1,2,3]], [0], [0])
306
        [0, 1, 2, 3]
307
        sage: sets_isom_py([[0,1,2,3]], [0,1], [1,2])
308
        False
309
        sage: sets_isom_py([[0,1,2,3]], [0,1], [0,1])
310
        [0, 1, 2, 3]
311
        sage: sets_isom_py([[0,1,2,3]], [0,1,2,3], [0,1])
312
        False
313
        sage: sets_isom_py([[0,1,2,3]], [0,1,2,3], [0,1,2,3])
314
        [0, 1, 2, 3]
315

316
    Four-cycle::
317

318
        sage: sets_isom_py([[1,2,3,0]], [0,1,2,3], [0,1,2,3])
319
        [0, 1, 2, 3]
320
        sage: sets_isom_py([[1,2,3,0]], [], [])
321
        [0, 1, 2, 3]
322
        sage: sets_isom_py([[1,2,3,0]], [0], [0])
323
        [0, 1, 2, 3]
324
        sage: sets_isom_py([[1,2,3,0]], [0], [1])
325
        [1, 2, 3, 0]
326
        sage: sets_isom_py([[1,2,3,0]], [0], [2])
327
        [2, 3, 0, 1]
328
        sage: sets_isom_py([[1,2,3,0]], [0], [3])
329
        [3, 0, 1, 2]
330
        sage: sets_isom_py([[1,2,3,0]], [0,1], [3])
331
        False
332
        sage: sets_isom_py([[1,2,3,0]], [0,1], [])
333
        False
334
        sage: sets_isom_py([[1,2,3,0]], [0,1], [1,2,3])
335
        False
336
        sage: sets_isom_py([[1,2,3,0]], [0,1], [1,2])
337
        [1, 2, 3, 0]
338
        sage: sets_isom_py([[1,2,3,0]], [0,1], [2,3])
339
        [2, 3, 0, 1]
340
        sage: sets_isom_py([[1,2,3,0]], [0,1], [0,3])
341
        [3, 0, 1, 2]
342
        sage: sets_isom_py([[1,2,3,0]], [0,2], [0,2])
343
        [0, 1, 2, 3]
344
        sage: sets_isom_py([[1,2,3,0]], [0,2], [1,3])
345
        [3, 0, 1, 2]
346
        sage: sets_isom_py([[1,2,3,0]], [0,1,2], [1,2,3])
347
        [1, 2, 3, 0]
348
        sage: sets_isom_py([[1,2,3,0]], [0,1,2], [0,2,3])
349
        [2, 3, 0, 1]
350
        sage: sets_isom_py([[1,2,3,0]], [0,1,2], [0,1,3])
351
        [3, 0, 1, 2]
352
        sage: sets_isom_py([[1,2,3,0]], [0,1,2], [0,1,2])
353
        [0, 1, 2, 3]
354

355
    Two transpositions::
356

357
        sage: from sage.groups.perm_gps.partn_ref.refinement_sets import sets_isom_py
358
        sage: sets_isom_py([[2,3,0,1]], [0], [2])
359
        [2, 3, 0, 1]
360
        sage: sets_isom_py([[2,3,0,1]], [1], [3])
361
        [2, 3, 0, 1]
362
        sage: sets_isom_py([[2,3,0,1]], [1], [1])
363
        [0, 1, 2, 3]
364
        sage: sets_isom_py([[2,3,0,1]], [1], [2])
365
        False
366
        sage: sets_isom_py([[2,3,0,1]], [0,3], [1,2])
367
        [2, 3, 0, 1]
368
        sage: sets_isom_py([[2,3,0,1]], [0,3], [3,0])
369
        [0, 1, 2, 3]
370
        sage: sets_isom_py([[2,3,0,1]], [0,1,3], [0,2,3])
371
        False
372
        sage: sets_isom_py([[2,3,0,1]], [0,1,3], [1,2,3])
373
        [2, 3, 0, 1]
374

375

376
    """
377
    from sage.misc.misc import uniq
378
    set1 = uniq(set1)
379
    set2 = uniq(set2)
380
    if len(generators) == 0:
381
        if set1 == set2:
382
            return range(max(set1)+1)
383
        else:
384
            return False
385
    cdef int i, j, n = len(generators[0]), n_gens = len(generators)
386
    cdef StabilizerChain *supergroup = SC_new(n)
387
    cdef int *gens = <int *> sage_malloc(n*n_gens * sizeof(int))
388
    cdef int *isom = <int *> sage_malloc(n * sizeof(int))
389
    cdef subset *subset_sett1 = <subset *> sage_malloc(sizeof(subset))
390
    cdef subset *subset_sett2 = <subset *> sage_malloc(sizeof(subset))
391
    bitset_init(&subset_sett1.bits, n)
392
    bitset_init(&subset_sett2.bits, n)
393
    subset_sett1.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
394
    subset_sett2.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
395
    for i from 0 <= i < len(generators):
396
        for j from 0 <= j < n:
397
            gens[n*i + j] = generators[i][j]
398
    if SC_insert(supergroup, 0, gens, n_gens):
399
        raise MemoryError
400
    sage_free(gens)
401
    bitset_clear(&subset_sett1.bits)
402
    bitset_clear(&subset_sett2.bits)
403
    for i in set1:
404
        bitset_add(&subset_sett1.bits, i)
405
    for i in set2:
406
        bitset_add(&subset_sett2.bits, i)
407
    cdef bint isomorphic = sets_isom(supergroup, subset_sett1, subset_sett2, isom)
408
    SC_dealloc(supergroup)
409
    bitset_free(&subset_sett1.bits)
410
    bitset_free(&subset_sett2.bits)
411
    sage_free(subset_sett1.scratch)
412
    sage_free(subset_sett2.scratch)
413
    sage_free(subset_sett1)
414
    sage_free(subset_sett2)
415
    if isomorphic:
416
        output_py = [isom[i] for i from 0 <= i < n]
417
    else:
418
        output_py = False
419
    sage_free(isom)
420
    return output_py
421

422
cdef int sets_isom(StabilizerChain *supergroup, subset *set1, subset *set2, int *isom) except -1:
423
    r"""
424
    Underlying C function for testing two sets for isomorphism.
425
    """
426
    cdef int n = supergroup.degree
427
    cdef bint x
428
    cdef PartitionStack *part = PS_new(n, 1)
429
    if part is NULL:
430
        raise MemoryError
431
    x = double_coset(set1, set2, part, NULL, n,
432
        &all_set_children_are_equivalent, &refine_set, &compare_sets,
433
        supergroup, NULL, isom)
434
    PS_dealloc(part)
435
    return x
436

437
cdef bint all_set_children_are_equivalent(PartitionStack *PS, void *S):
438
    return 0
439

440
cdef int refine_set(PartitionStack *PS, void *S, int *cells_to_refine_by, int ctrb_len):
441
    """
442
    Given a set S, refine the partition stack PS so that each cell contains
443
    elements which are all either in the set or not in the set. If the depth is
444
    positive we do nothing since the cells will have already been split at an
445
    earlier level.
446
    """
447
    if PS.depth > 0:
448
        return 0
449
    cdef subset *subset1 = <subset *> S
450
    cdef int *scratch = subset1.scratch
451
    cdef int start, i, n = PS.degree, x
452
    start = 0
453
    while start < n:
454
        i = 0
455
        while 1:
456
            scratch[i] = bitset_in(&subset1.bits, PS.entries[start+i])
457
            if PS.levels[start+i] <= PS.depth:
458
                break
459
            i += 1
460
        sort_by_function(PS, start, scratch)
461
        start += i+1
462
    return 0
463

464
cdef inline int _bint_cmp(bint a, bint b):
465
    return (<int> b) - (<int> a)
466

467
cdef int compare_sets(int *gamma_1, int *gamma_2, void *S1, void *S2, int degree):
468
    r"""
469
    Compare two sets according to the lexicographic order.
470
    """
471
    cdef subset *subset1 = <subset *> S1
472
    cdef subset *subset2 = <subset *> S2
473
    cdef bitset_s set1 = subset1.bits
474
    cdef bitset_s set2 = subset2.bits
475
    cdef int i, j
476
    for i from 0 <= i < degree:
477
        j = _bint_cmp(bitset_in(&set1, gamma_1[i]), bitset_in(&set2, gamma_2[i]))
478
        if j != 0: return j
479
    return 0
480

481
cdef void *allocate_subset(int n):
482
    r"""
483
    Allocates a subset struct of degree n.
484
    """
485
    cdef subset *set1 = <subset *> sage_malloc(sizeof(subset))
486
    cdef int *scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
487
    if set1 is NULL or scratch is NULL:
488
        sage_free(set1)
489
        sage_free(scratch)
490
        return NULL
491
    try:
492
        bitset_init(&set1.bits, n)
493
    except MemoryError:
494
        sage_free(set1)
495
        sage_free(scratch)
496
        return NULL
497
    set1.scratch = scratch
498
    return <void *> set1
499

500
cdef void free_subset(void *child):
501
    r"""
502
    Deallocates a subset struct.
503
    """
504
    cdef subset *set1 = <subset *> child
505
    if set1 is not NULL:
506
        sage_free(set1.scratch)
507
        bitset_free(&set1.bits)
508
    sage_free(set1)
509

510
cdef void *allocate_sgd(int degree):
511
    r"""
512
    Allocates the data part of an iterator which generates augmentations, i.e.,
513
    elements to add to the set.
514
    """
515
    cdef subset_generator_data *sgd = <subset_generator_data *> sage_malloc(sizeof(subset_generator_data))
516
    sgd.orbits = OP_new(degree)
517
    if sgd is NULL or sgd.orbits is NULL:
518
        deallocate_sgd(sgd)
519
        return NULL
520
    return <void *> sgd
521

522
cdef void deallocate_sgd(void *data):
523
    r"""
524
    Deallocates the data part of the augmentation iterator.
525
    """
526
    cdef subset_generator_data *sgd = <subset_generator_data *> data
527
    if sgd is not NULL:
528
        OP_dealloc(sgd.orbits)
529
    sage_free(sgd)
530

531
cdef void *subset_generator_next(void *data, int *degree, bint *mem_err):
532
    r"""
533
    Returns the next element to consider adding to the set.
534
    """
535
    cdef subset_generator_data *sgd = <subset_generator_data *> data
536
    while 1:
537
        sgd.cur_point += 1
538
        if sgd.cur_point == sgd.orbits.degree:
539
            break
540
        if OP_find(sgd.orbits, sgd.cur_point) == sgd.cur_point and \
541
          not bitset_in(&sgd.bits, sgd.cur_point):
542
            break
543
    if sgd.cur_point == sgd.orbits.degree or mem_err[0]:
544
        return NULL
545
    return <void *> &sgd.cur_point
546

547
cdef int generate_child_subsets(void *S, aut_gp_and_can_lab *group, iterator *child_iterator):
548
    r"""
549
    Sets up an iterator of augmentations, i.e., elements to add to the given set.
550
    """
551
    cdef subset *subset1 = <subset *> S
552
    cdef bitset_s set1 = subset1.bits
553
    cdef int i, j, n = group.group.degree
554
    cdef subset_generator_data *sgd = <subset_generator_data *> child_iterator.data
555
    OP_clear(sgd.orbits)
556
    for i from 0 <= i < group.num_gens:
557
        for j from 0 <= j < n:
558
            OP_join(sgd.orbits, j, group.generators[n*i + j])
559
    i = bitset_first(&subset1.bits)
560
    j = bitset_next(&subset1.bits, i+1)
561
    while j != -1:
562
        OP_join(sgd.orbits, i, j)
563
        j = bitset_next(&subset1.bits, j+1)
564
    sgd.cur_point = -1
565
    sgd.bits = subset1.bits
566
    return 0
567

568
cdef void *apply_subset_aug(void *parent, void *aug, void *child, int *degree, bint *mem_err):
569
    r"""
570
    Adds the element represented by ``aug`` to ``parent``, storing the result to
571
    ``child``.
572
    """
573
    cdef subset *set1 = <subset *> child, *par_set = <subset *> parent
574
    cdef bitset_s parbits = par_set.bits
575
    cdef int add_pt = (<int *> aug)[0], n = parbits.size
576
    bitset_copy(&set1.bits, &parbits)
577
    bitset_add(&set1.bits, add_pt)
578
    degree[0] = n
579
    return <void *> set1
580

581
cdef void free_subset_aug(void *aug):
582
    return
583

584
cdef void *canonical_set_parent(void *child, void *parent, int *permutation, int *degree, bint *mem_err):
585
    r"""
586
    Determines the canonical parent of the set ``child`` by applying
587
    ``permutation``, deleting the largest element in lexicographic order, and
588
    storing the result to ``parent``.
589
    """
590
    cdef subset *set1 = <subset *> child
591
    cdef bitset_t can_par
592
    cdef int i, max_in_can_lab, max_loc, n = set1.bits.size
593
    cdef subset *par
594
    if parent is NULL:
595
        par = <subset *> allocate_subset(n)
596
        if par is NULL:
597
            mem_err[0] = 1
598
    else:
599
        par = <subset *> parent
600
    if mem_err[0]:
601
        return NULL
602
    i = bitset_first(&set1.bits)
603
    max_in_can_lab = permutation[i]
604
    max_loc = i
605
    while i != -1:
606
        if max_in_can_lab < permutation[i]:
607
            max_in_can_lab = permutation[i]
608
            max_loc = i
609
        i = bitset_next(&set1.bits, i+1)
610
    bitset_copy(&par.bits, &set1.bits)
611
    bitset_discard(&par.bits, max_loc)
612
    degree[0] = n
613
    return <void *> par
614

615
cdef iterator *allocate_subset_gen(int degree, int max_size):
616
    r"""
617
    Allocates the generator of subsets.
618
    """
619
    cdef iterator *subset_gen = <iterator *> sage_malloc(sizeof(iterator))
620
    if subset_gen is not NULL:
621
        if allocate_subset_gen_2(degree, max_size, subset_gen):
622
            sage_free(subset_gen)
623
            subset_gen = NULL
624
    return subset_gen
625

626
cdef int allocate_subset_gen_2(int degree, int max_size, iterator *it):
627
    r"""
628
    Given an already allocated iterator, allocates the generator of subsets.
629
    """
630
    cdef canonical_generator_data *cgd = allocate_cgd(max_size + 1, degree)
631
    if cgd is NULL:
632
        return 1
633
    cdef int i, j
634
    for i from 0 <= i < max_size + 1:
635
        cgd.object_stack[i] = allocate_subset(degree)
636
        cgd.parent_stack[i] = allocate_subset(degree)
637
        cgd.iterator_stack[i].data = allocate_sgd(degree)
638
        cgd.iterator_stack[i].next = &subset_generator_next
639
        if cgd.iterator_stack[i].data is NULL or \
640
           cgd.object_stack[i]        is NULL or \
641
           cgd.parent_stack[i]        is NULL:
642
            for j from 0 <= j <= i:
643
                deallocate_sgd(cgd.iterator_stack[i].data)
644
                free_subset(cgd.object_stack[i])
645
                free_subset(cgd.parent_stack[i])
646
            deallocate_cgd(cgd)
647
            return 1
648
    it.data = <void *> cgd
649
    it.next = &canonical_generator_next
650
    return 0
651

652
cdef void free_subset_gen(iterator *subset_gen):
653
    r"""
654
    Frees the iterator of subsets.
655
    """
656
    if subset_gen is NULL: return
657
    cdef canonical_generator_data *cgd = <canonical_generator_data *> subset_gen.data
658
    deallocate_cgd(cgd)
659
    sage_free(subset_gen)
660

661
cdef iterator *setup_set_gen(iterator *subset_gen, int degree, int max_size):
662
    r"""
663
    Initiates the iterator of subsets.
664
    """
665
    cdef subset *empty_set
666
    cdef iterator *subset_iterator = setup_canonical_generator(degree,
667
        &all_set_children_are_equivalent,
668
        &refine_set,
669
        &compare_sets,
670
        &generate_child_subsets,
671
        &apply_subset_aug,
672
        &free_subset,
673
        &deallocate_sgd,
674
        &free_subset_aug,
675
        &canonical_set_parent,
676
        max_size+1, 0, subset_gen)
677
    if subset_iterator is not NULL:
678
        empty_set = <subset *> (<canonical_generator_data *> subset_gen.data).object_stack[0]
679
        bitset_clear(&empty_set.bits)
680
    return subset_iterator
681

682
def sets_modulo_perm_group(list generators, int max_size, bint indicate_mem_err = 1):
683
    r"""
684
    Given generators of a permutation group, list subsets up to permutations in
685
    the group.
686

687
    INPUT:
688

689
        - ``generators`` - (list of lists) list of generators in list form
690
        - ``max_size`` - (int) maximum size of subsets to be generated
691
        - ``indicate_mem_err`` - (bool) whether to raise an error
692
            if we run out of memory, or simply append a MemoryError
693
            instance to the end of the output
694

695
    EXAMPLES::
696

697
        sage: from sage.groups.perm_gps.partn_ref.refinement_sets import sets_modulo_perm_group
698
        sage: sets_modulo_perm_group([], 0)
699
        [[]]
700
        sage: sets_modulo_perm_group([], 1)
701
        [[0], []]
702
        sage: sets_modulo_perm_group([], 2)
703
        [[0, 1], [0], []]
704
        sage: sets_modulo_perm_group([], 3)
705
        [[0, 1, 2], [0, 1], [0], []]
706
        sage: sets_modulo_perm_group([], 4)
707
        [[0, 1, 2, 3], [0, 1, 2], [0, 1], [0], []]
708
        sage: len(sets_modulo_perm_group([], 99))
709
        100
710

711
    ::
712

713
        sage: sets_modulo_perm_group([[1,2,0]], 4)
714
        [[0, 1, 2], [0, 1], [0], []]
715
        sage: sets_modulo_perm_group([[1,2,0]], 3)
716
        [[0, 1, 2], [0, 1], [0], []]
717
        sage: sets_modulo_perm_group([[1,2,0]], 2)
718
        [[0, 1], [0], []]
719
        sage: sets_modulo_perm_group([[1,2,0]], 1)
720
        [[0], []]
721
        sage: sets_modulo_perm_group([[1,2,0]], 0)
722
        [[]]
723
        sage: sets_modulo_perm_group([[0,1,2]], 3)
724
        [[0, 1, 2], [0, 1], [0, 2], [0], [1, 2], [1], [2], []]
725
        sage: sets_modulo_perm_group([[1,0,2]], 3)
726
        [[0, 1, 2], [0, 1], [0, 2], [0], [2], []]
727
        sage: sets_modulo_perm_group([[1,0,2],[1,2,0]], 3)
728
        [[0, 1, 2], [0, 1], [1], []]
729

730
    ::
731

732
        sage: sets_modulo_perm_group([[1,2,3,0]], 4)
733
        [[0, 1, 2, 3], [0, 1, 2], [0, 1], [0, 2], [0], []]
734
        sage: sets_modulo_perm_group([[1,2,3,0]], 5)
735
        [[0, 1, 2, 3], [0, 1, 2], [0, 1], [0, 2], [0], []]
736
        sage: sets_modulo_perm_group([[1,2,3,0]], 3)
737
        [[0, 1, 2], [0, 1], [0, 2], [0], []]
738
        sage: sets_modulo_perm_group([[1,2,3,0]], 2)
739
        [[0, 1], [0, 2], [0], []]
740
        sage: sets_modulo_perm_group([[1,2,3,0]], 1)
741
        [[0], []]
742
        sage: sets_modulo_perm_group([[1,2,3,0]], 0)
743
        [[]]
744
        sage: sets_modulo_perm_group([[0,1,3,2],[1,0,2,3]], 4)
745
        [[0, 1, 2, 3], [0, 1, 2], [0, 1], [0, 2, 3], [0, 2], [0], [2, 3], [2], []]
746
        sage: sets_modulo_perm_group([[1,0,2,3],[1,2,0,3]], 4)
747
        [[0, 1, 2, 3], [0, 1, 2], [0, 1, 3], [0, 1], [1, 3], [1], [3], []]
748
        sage: sets_modulo_perm_group([[1,2,0,3],[0,2,3,1]], 4)
749
        [[0, 1, 2, 3], [0, 1, 2], [0, 1], [1], []]
750
        sage: sets_modulo_perm_group([[1,0,2,3],[1,2,3,0]], 4)
751
        [[0, 1, 2, 3], [0, 1, 2], [1, 2], [2], []]
752
        sage: L = list(powerset(range(4)))
753
        sage: L.sort()
754
        sage: L == sorted(sets_modulo_perm_group([[0,1,2,3]], 4))
755
        True
756

757
    ::
758

759
        sage: sets_modulo_perm_group([[1,2,3,4,0]], 5)
760
        [[0, 1, 2, 3, 4], [0, 1, 2, 3], [0, 1, 2], [0, 1], [0, 2, 3], [0, 2], [0], []]
761
        sage: sets_modulo_perm_group([[1,0,2,3,4],[0,1,3,4,2]], 5)
762
        [[0, 1, 2, 3, 4], [0, 1, 2, 3], [0, 1, 2], [0, 1], [0, 2, 3, 4], [0, 2, 3], [0, 2], [0], [2, 3, 4], [2, 3], [2], []]
763
        sage: L = list(powerset(range(5)))
764
        sage: L.sort()
765
        sage: L == sorted(sets_modulo_perm_group([[0,1,2,3,4]], 5))
766
        True
767
        sage: sets_modulo_perm_group([[1,0,2,3,4],[1,2,3,4,0]], 5)
768
        [[0, 1, 2, 3, 4], [0, 1, 2, 3], [1, 2, 3], [2, 3], [3], []]
769
        sage: sets_modulo_perm_group([[1,2,0,3,4],[0,2,3,1,4],[0,1,3,4,2]], 5)
770
        [[0, 1, 2, 3, 4], [0, 1, 2, 3], [1, 2, 3], [1, 2], [2], []]
771

772
    ::
773

774
        sage: X = sets_modulo_perm_group([[1,2,3,4,5,0]], 6)
775
        sage: [a for a in X if len(a) == 0]
776
        [[]]
777
        sage: [a for a in X if len(a) == 1]
778
        [[0]]
779
        sage: [a for a in X if len(a) == 2]
780
        [[0, 1], [0, 2], [0, 3]]
781
        sage: [a for a in X if len(a) == 3]
782
        [[0, 1, 2], [0, 1, 3], [0, 2, 3], [0, 2, 4]]
783
        sage: [a for a in X if len(a) == 4]
784
        [[0, 1, 2, 3], [0, 1, 3, 4], [0, 2, 3, 4]]
785
        sage: [a for a in X if len(a) == 5]
786
        [[0, 1, 2, 3, 4]]
787
        sage: [a for a in X if len(a) == 6]
788
        [[0, 1, 2, 3, 4, 5]]
789

790
    ::
791

792
        sage: X = sets_modulo_perm_group([[0,2,1,4,3,5,8,7,6],[8,7,6,3,5,4,2,1,0]], 9)
793
        sage: len(X)
794
        74
795

796
    """
797
    cdef list out_list = []
798
    cdef int i
799
    if max_size == 0:
800
        return [[]]
801
    if len(generators) == 0:
802
        ll = []
803
        for i in range(max_size,-1,-1):
804
            ll.append(range(i))
805
        return ll
806
    cdef int n = len(generators[0]), n_gens = len(generators)
807
    cdef iterator *subset_iterator
808
    cdef subset *thing
809

810
    cdef StabilizerChain *group = SC_new(n)
811
    cdef int *gens = <int *> sage_malloc(n*n_gens * sizeof(int))
812
    if group is NULL or gens is NULL:
813
        SC_dealloc(group)
814
        sage_free(gens)
815
        raise MemoryError
816
    for i from 0 <= i < len(generators):
817
        for j from 0 <= j < n:
818
            gens[n*i + j] = generators[i][j]
819
    if SC_insert(group, 0, gens, n_gens):
820
        SC_dealloc(group)
821
        sage_free(gens)
822
        raise MemoryError
823
    sage_free(gens)
824

825
    cdef iterator *subset_gen = allocate_subset_gen(n, max_size)
826
    if subset_gen is NULL:
827
        SC_dealloc(group)
828
        raise MemoryError
829
    subset_iterator = setup_set_gen(subset_gen, n, max_size)
830
    cdef bint mem_err = 0
831
    if subset_iterator is NULL:
832
        SC_dealloc(group)
833
        free_subset_gen(subset_gen)
834
        mem_err = 1
835
    else:
836
        start_canonical_generator(group, NULL, n, subset_gen)
837
    while not mem_err:
838
        thing = <subset *> subset_iterator.next(subset_iterator.data, NULL, &mem_err)
839
        if thing is NULL: break
840
        out_list.append( bitset_list(&thing.bits) )
841
    free_subset_gen(subset_gen)
842
    SC_dealloc(group)
843
    if mem_err:
844
        if indicate_mem_err:
845
            raise MemoryError
846
        else:
847
            out_list.append(MemoryError())
848
    return out_list
849

850

851

852

853

854

855

856
Product

Resources

Company
