Coverage for /usr/lib/python3/dist-packages/sympy/combinatorics/perm_groups.py: 9%

2114 statements  

« prev     ^ index     » next       coverage.py v7.9.1, created at 2025-06-14 15:55 +0200

1from math import factorial as _factorial, log, prod 

2from itertools import chain, islice, product 

3 

4 

5from sympy.combinatorics import Permutation 

6from sympy.combinatorics.permutations import (_af_commutes_with, _af_invert, 

7 _af_rmul, _af_rmuln, _af_pow, Cycle) 

8from sympy.combinatorics.util import (_check_cycles_alt_sym, 

9 _distribute_gens_by_base, _orbits_transversals_from_bsgs, 

10 _handle_precomputed_bsgs, _base_ordering, _strong_gens_from_distr, 

11 _strip, _strip_af) 

12from sympy.core import Basic 

13from sympy.core.random import _randrange, randrange, choice 

14from sympy.core.symbol import Symbol 

15from sympy.core.sympify import _sympify 

16from sympy.functions.combinatorial.factorials import factorial 

17from sympy.ntheory import primefactors, sieve 

18from sympy.ntheory.factor_ import (factorint, multiplicity) 

19from sympy.ntheory.primetest import isprime 

20from sympy.utilities.iterables import has_variety, is_sequence, uniq 

21 

22rmul = Permutation.rmul_with_af 

23_af_new = Permutation._af_new 

24 

25 

26class PermutationGroup(Basic): 

27 r"""The class defining a Permutation group. 

28 

29 Explanation 

30 =========== 

31 

32 ``PermutationGroup([p1, p2, ..., pn])`` returns the permutation group 

33 generated by the list of permutations. This group can be supplied 

34 to Polyhedron if one desires to decorate the elements to which the 

35 indices of the permutation refer. 

36 

37 Examples 

38 ======== 

39 

40 >>> from sympy.combinatorics import Permutation, PermutationGroup 

41 >>> from sympy.combinatorics import Polyhedron 

42 

43 The permutations corresponding to motion of the front, right and 

44 bottom face of a $2 \times 2$ Rubik's cube are defined: 

45 

46 >>> F = Permutation(2, 19, 21, 8)(3, 17, 20, 10)(4, 6, 7, 5) 

47 >>> R = Permutation(1, 5, 21, 14)(3, 7, 23, 12)(8, 10, 11, 9) 

48 >>> D = Permutation(6, 18, 14, 10)(7, 19, 15, 11)(20, 22, 23, 21) 

49 

50 These are passed as permutations to PermutationGroup: 

51 

52 >>> G = PermutationGroup(F, R, D) 

53 >>> G.order() 

54 3674160 

55 

56 The group can be supplied to a Polyhedron in order to track the 

57 objects being moved. An example involving the $2 \times 2$ Rubik's cube is 

58 given there, but here is a simple demonstration: 

59 

60 >>> a = Permutation(2, 1) 

61 >>> b = Permutation(1, 0) 

62 >>> G = PermutationGroup(a, b) 

63 >>> P = Polyhedron(list('ABC'), pgroup=G) 

64 >>> P.corners 

65 (A, B, C) 

66 >>> P.rotate(0) # apply permutation 0 

67 >>> P.corners 

68 (A, C, B) 

69 >>> P.reset() 

70 >>> P.corners 

71 (A, B, C) 

72 

73 Or one can make a permutation as a product of selected permutations 

74 and apply them to an iterable directly: 

75 

76 >>> P10 = G.make_perm([0, 1]) 

77 >>> P10('ABC') 

78 ['C', 'A', 'B'] 

79 

80 See Also 

81 ======== 

82 

83 sympy.combinatorics.polyhedron.Polyhedron, 

84 sympy.combinatorics.permutations.Permutation 

85 

86 References 

87 ========== 

88 

89 .. [1] Holt, D., Eick, B., O'Brien, E. 

90 "Handbook of Computational Group Theory" 

91 

92 .. [2] Seress, A. 

93 "Permutation Group Algorithms" 

94 

95 .. [3] https://en.wikipedia.org/wiki/Schreier_vector 

96 

97 .. [4] https://en.wikipedia.org/wiki/Nielsen_transformation#Product_replacement_algorithm 

98 

99 .. [5] Frank Celler, Charles R.Leedham-Green, Scott H.Murray, 

100 Alice C.Niemeyer, and E.A.O'Brien. "Generating Random 

101 Elements of a Finite Group" 

102 

103 .. [6] https://en.wikipedia.org/wiki/Block_%28permutation_group_theory%29 

104 

105 .. [7] https://algorithmist.com/wiki/Union_find 

106 

107 .. [8] https://en.wikipedia.org/wiki/Multiply_transitive_group#Multiply_transitive_groups 

108 

109 .. [9] https://en.wikipedia.org/wiki/Center_%28group_theory%29 

110 

111 .. [10] https://en.wikipedia.org/wiki/Centralizer_and_normalizer 

112 

113 .. [11] https://groupprops.subwiki.org/wiki/Derived_subgroup 

114 

115 .. [12] https://en.wikipedia.org/wiki/Nilpotent_group 

116 

117 .. [13] https://www.math.colostate.edu/~hulpke/CGT/cgtnotes.pdf 

118 

119 .. [14] https://docs.gap-system.org/doc/ref/manual.pdf 

120 

121 """ 

122 is_group = True 

123 

124 def __new__(cls, *args, dups=True, **kwargs): 

125 """The default constructor. Accepts Cycle and Permutation forms. 

126 Removes duplicates unless ``dups`` keyword is ``False``. 

127 """ 

128 if not args: 

129 args = [Permutation()] 

130 else: 

131 args = list(args[0] if is_sequence(args[0]) else args) 

132 if not args: 

133 args = [Permutation()] 

134 if any(isinstance(a, Cycle) for a in args): 

135 args = [Permutation(a) for a in args] 

136 if has_variety(a.size for a in args): 

137 degree = kwargs.pop('degree', None) 

138 if degree is None: 

139 degree = max(a.size for a in args) 

140 for i in range(len(args)): 

141 if args[i].size != degree: 

142 args[i] = Permutation(args[i], size=degree) 

143 if dups: 

144 args = list(uniq([_af_new(list(a)) for a in args])) 

145 if len(args) > 1: 

146 args = [g for g in args if not g.is_identity] 

147 return Basic.__new__(cls, *args, **kwargs) 

148 

149 def __init__(self, *args, **kwargs): 

150 self._generators = list(self.args) 

151 self._order = None 

152 self._center = [] 

153 self._is_abelian = None 

154 self._is_transitive = None 

155 self._is_sym = None 

156 self._is_alt = None 

157 self._is_primitive = None 

158 self._is_nilpotent = None 

159 self._is_solvable = None 

160 self._is_trivial = None 

161 self._transitivity_degree = None 

162 self._max_div = None 

163 self._is_perfect = None 

164 self._is_cyclic = None 

165 self._is_dihedral = None 

166 self._r = len(self._generators) 

167 self._degree = self._generators[0].size 

168 

169 # these attributes are assigned after running schreier_sims 

170 self._base = [] 

171 self._strong_gens = [] 

172 self._strong_gens_slp = [] 

173 self._basic_orbits = [] 

174 self._transversals = [] 

175 self._transversal_slp = [] 

176 

177 # these attributes are assigned after running _random_pr_init 

178 self._random_gens = [] 

179 

180 # finite presentation of the group as an instance of `FpGroup` 

181 self._fp_presentation = None 

182 

183 def __getitem__(self, i): 

184 return self._generators[i] 

185 

186 def __contains__(self, i): 

187 """Return ``True`` if *i* is contained in PermutationGroup. 

188 

189 Examples 

190 ======== 

191 

192 >>> from sympy.combinatorics import Permutation, PermutationGroup 

193 >>> p = Permutation(1, 2, 3) 

194 >>> Permutation(3) in PermutationGroup(p) 

195 True 

196 

197 """ 

198 if not isinstance(i, Permutation): 

199 raise TypeError("A PermutationGroup contains only Permutations as " 

200 "elements, not elements of type %s" % type(i)) 

201 return self.contains(i) 

202 

203 def __len__(self): 

204 return len(self._generators) 

205 

206 def equals(self, other): 

207 """Return ``True`` if PermutationGroup generated by elements in the 

208 group are same i.e they represent the same PermutationGroup. 

209 

210 Examples 

211 ======== 

212 

213 >>> from sympy.combinatorics import Permutation, PermutationGroup 

214 >>> p = Permutation(0, 1, 2, 3, 4, 5) 

215 >>> G = PermutationGroup([p, p**2]) 

216 >>> H = PermutationGroup([p**2, p]) 

217 >>> G.generators == H.generators 

218 False 

219 >>> G.equals(H) 

220 True 

221 

222 """ 

223 if not isinstance(other, PermutationGroup): 

224 return False 

225 

226 set_self_gens = set(self.generators) 

227 set_other_gens = set(other.generators) 

228 

229 # before reaching the general case there are also certain 

230 # optimisation and obvious cases requiring less or no actual 

231 # computation. 

232 if set_self_gens == set_other_gens: 

233 return True 

234 

235 # in the most general case it will check that each generator of 

236 # one group belongs to the other PermutationGroup and vice-versa 

237 for gen1 in set_self_gens: 

238 if not other.contains(gen1): 

239 return False 

240 for gen2 in set_other_gens: 

241 if not self.contains(gen2): 

242 return False 

243 return True 

244 

245 def __mul__(self, other): 

246 """ 

247 Return the direct product of two permutation groups as a permutation 

248 group. 

249 

250 Explanation 

251 =========== 

252 

253 This implementation realizes the direct product by shifting the index 

254 set for the generators of the second group: so if we have ``G`` acting 

255 on ``n1`` points and ``H`` acting on ``n2`` points, ``G*H`` acts on 

256 ``n1 + n2`` points. 

257 

258 Examples 

259 ======== 

260 

261 >>> from sympy.combinatorics.named_groups import CyclicGroup 

262 >>> G = CyclicGroup(5) 

263 >>> H = G*G 

264 >>> H 

265 PermutationGroup([ 

266 (9)(0 1 2 3 4), 

267 (5 6 7 8 9)]) 

268 >>> H.order() 

269 25 

270 

271 """ 

272 if isinstance(other, Permutation): 

273 return Coset(other, self, dir='+') 

274 gens1 = [perm._array_form for perm in self.generators] 

275 gens2 = [perm._array_form for perm in other.generators] 

276 n1 = self._degree 

277 n2 = other._degree 

278 start = list(range(n1)) 

279 end = list(range(n1, n1 + n2)) 

280 for i in range(len(gens2)): 

281 gens2[i] = [x + n1 for x in gens2[i]] 

282 gens2 = [start + gen for gen in gens2] 

283 gens1 = [gen + end for gen in gens1] 

284 together = gens1 + gens2 

285 gens = [_af_new(x) for x in together] 

286 return PermutationGroup(gens) 

287 

288 def _random_pr_init(self, r, n, _random_prec_n=None): 

289 r"""Initialize random generators for the product replacement algorithm. 

290 

291 Explanation 

292 =========== 

293 

294 The implementation uses a modification of the original product 

295 replacement algorithm due to Leedham-Green, as described in [1], 

296 pp. 69-71; also, see [2], pp. 27-29 for a detailed theoretical 

297 analysis of the original product replacement algorithm, and [4]. 

298 

299 The product replacement algorithm is used for producing random, 

300 uniformly distributed elements of a group `G` with a set of generators 

301 `S`. For the initialization ``_random_pr_init``, a list ``R`` of 

302 `\max\{r, |S|\}` group generators is created as the attribute 

303 ``G._random_gens``, repeating elements of `S` if necessary, and the 

304 identity element of `G` is appended to ``R`` - we shall refer to this 

305 last element as the accumulator. Then the function ``random_pr()`` 

306 is called ``n`` times, randomizing the list ``R`` while preserving 

307 the generation of `G` by ``R``. The function ``random_pr()`` itself 

308 takes two random elements ``g, h`` among all elements of ``R`` but 

309 the accumulator and replaces ``g`` with a randomly chosen element 

310 from `\{gh, g(~h), hg, (~h)g\}`. Then the accumulator is multiplied 

311 by whatever ``g`` was replaced by. The new value of the accumulator is 

312 then returned by ``random_pr()``. 

313 

314 The elements returned will eventually (for ``n`` large enough) become 

315 uniformly distributed across `G` ([5]). For practical purposes however, 

316 the values ``n = 50, r = 11`` are suggested in [1]. 

317 

318 Notes 

319 ===== 

320 

321 THIS FUNCTION HAS SIDE EFFECTS: it changes the attribute 

322 self._random_gens 

323 

324 See Also 

325 ======== 

326 

327 random_pr 

328 

329 """ 

330 deg = self.degree 

331 random_gens = [x._array_form for x in self.generators] 

332 k = len(random_gens) 

333 if k < r: 

334 for i in range(k, r): 

335 random_gens.append(random_gens[i - k]) 

336 acc = list(range(deg)) 

337 random_gens.append(acc) 

338 self._random_gens = random_gens 

339 

340 # handle randomized input for testing purposes 

341 if _random_prec_n is None: 

342 for i in range(n): 

343 self.random_pr() 

344 else: 

345 for i in range(n): 

346 self.random_pr(_random_prec=_random_prec_n[i]) 

347 

348 def _union_find_merge(self, first, second, ranks, parents, not_rep): 

349 """Merges two classes in a union-find data structure. 

350 

351 Explanation 

352 =========== 

353 

354 Used in the implementation of Atkinson's algorithm as suggested in [1], 

355 pp. 83-87. The class merging process uses union by rank as an 

356 optimization. ([7]) 

357 

358 Notes 

359 ===== 

360 

361 THIS FUNCTION HAS SIDE EFFECTS: the list of class representatives, 

362 ``parents``, the list of class sizes, ``ranks``, and the list of 

363 elements that are not representatives, ``not_rep``, are changed due to 

364 class merging. 

365 

366 See Also 

367 ======== 

368 

369 minimal_block, _union_find_rep 

370 

371 References 

372 ========== 

373 

374 .. [1] Holt, D., Eick, B., O'Brien, E. 

375 "Handbook of computational group theory" 

376 

377 .. [7] https://algorithmist.com/wiki/Union_find 

378 

379 """ 

380 rep_first = self._union_find_rep(first, parents) 

381 rep_second = self._union_find_rep(second, parents) 

382 if rep_first != rep_second: 

383 # union by rank 

384 if ranks[rep_first] >= ranks[rep_second]: 

385 new_1, new_2 = rep_first, rep_second 

386 else: 

387 new_1, new_2 = rep_second, rep_first 

388 total_rank = ranks[new_1] + ranks[new_2] 

389 if total_rank > self.max_div: 

390 return -1 

391 parents[new_2] = new_1 

392 ranks[new_1] = total_rank 

393 not_rep.append(new_2) 

394 return 1 

395 return 0 

396 

397 def _union_find_rep(self, num, parents): 

398 """Find representative of a class in a union-find data structure. 

399 

400 Explanation 

401 =========== 

402 

403 Used in the implementation of Atkinson's algorithm as suggested in [1], 

404 pp. 83-87. After the representative of the class to which ``num`` 

405 belongs is found, path compression is performed as an optimization 

406 ([7]). 

407 

408 Notes 

409 ===== 

410 

411 THIS FUNCTION HAS SIDE EFFECTS: the list of class representatives, 

412 ``parents``, is altered due to path compression. 

413 

414 See Also 

415 ======== 

416 

417 minimal_block, _union_find_merge 

418 

419 References 

420 ========== 

421 

422 .. [1] Holt, D., Eick, B., O'Brien, E. 

423 "Handbook of computational group theory" 

424 

425 .. [7] https://algorithmist.com/wiki/Union_find 

426 

427 """ 

428 rep, parent = num, parents[num] 

429 while parent != rep: 

430 rep = parent 

431 parent = parents[rep] 

432 # path compression 

433 temp, parent = num, parents[num] 

434 while parent != rep: 

435 parents[temp] = rep 

436 temp = parent 

437 parent = parents[temp] 

438 return rep 

439 

440 @property 

441 def base(self): 

442 r"""Return a base from the Schreier-Sims algorithm. 

443 

444 Explanation 

445 =========== 

446 

447 For a permutation group `G`, a base is a sequence of points 

448 `B = (b_1, b_2, \dots, b_k)` such that no element of `G` apart 

449 from the identity fixes all the points in `B`. The concepts of 

450 a base and strong generating set and their applications are 

451 discussed in depth in [1], pp. 87-89 and [2], pp. 55-57. 

452 

453 An alternative way to think of `B` is that it gives the 

454 indices of the stabilizer cosets that contain more than the 

455 identity permutation. 

456 

457 Examples 

458 ======== 

459 

460 >>> from sympy.combinatorics import Permutation, PermutationGroup 

461 >>> G = PermutationGroup([Permutation(0, 1, 3)(2, 4)]) 

462 >>> G.base 

463 [0, 2] 

464 

465 See Also 

466 ======== 

467 

468 strong_gens, basic_transversals, basic_orbits, basic_stabilizers 

469 

470 """ 

471 if self._base == []: 

472 self.schreier_sims() 

473 return self._base 

474 

475 def baseswap(self, base, strong_gens, pos, randomized=False, 

476 transversals=None, basic_orbits=None, strong_gens_distr=None): 

477 r"""Swap two consecutive base points in base and strong generating set. 

478 

479 Explanation 

480 =========== 

481 

482 If a base for a group `G` is given by `(b_1, b_2, \dots, b_k)`, this 

483 function returns a base `(b_1, b_2, \dots, b_{i+1}, b_i, \dots, b_k)`, 

484 where `i` is given by ``pos``, and a strong generating set relative 

485 to that base. The original base and strong generating set are not 

486 modified. 

487 

488 The randomized version (default) is of Las Vegas type. 

489 

490 Parameters 

491 ========== 

492 

493 base, strong_gens 

494 The base and strong generating set. 

495 pos 

496 The position at which swapping is performed. 

497 randomized 

498 A switch between randomized and deterministic version. 

499 transversals 

500 The transversals for the basic orbits, if known. 

501 basic_orbits 

502 The basic orbits, if known. 

503 strong_gens_distr 

504 The strong generators distributed by basic stabilizers, if known. 

505 

506 Returns 

507 ======= 

508 

509 (base, strong_gens) 

510 ``base`` is the new base, and ``strong_gens`` is a generating set 

511 relative to it. 

512 

513 Examples 

514 ======== 

515 

516 >>> from sympy.combinatorics.named_groups import SymmetricGroup 

517 >>> from sympy.combinatorics.testutil import _verify_bsgs 

518 >>> from sympy.combinatorics.perm_groups import PermutationGroup 

519 >>> S = SymmetricGroup(4) 

520 >>> S.schreier_sims() 

521 >>> S.base 

522 [0, 1, 2] 

523 >>> base, gens = S.baseswap(S.base, S.strong_gens, 1, randomized=False) 

524 >>> base, gens 

525 ([0, 2, 1], 

526 [(0 1 2 3), (3)(0 1), (1 3 2), 

527 (2 3), (1 3)]) 

528 

529 check that base, gens is a BSGS 

530 

531 >>> S1 = PermutationGroup(gens) 

532 >>> _verify_bsgs(S1, base, gens) 

533 True 

534 

535 See Also 

536 ======== 

537 

538 schreier_sims 

539 

540 Notes 

541 ===== 

542 

543 The deterministic version of the algorithm is discussed in 

544 [1], pp. 102-103; the randomized version is discussed in [1], p.103, and 

545 [2], p.98. It is of Las Vegas type. 

546 Notice that [1] contains a mistake in the pseudocode and 

547 discussion of BASESWAP: on line 3 of the pseudocode, 

548 `|\beta_{i+1}^{\left\langle T\right\rangle}|` should be replaced by 

549 `|\beta_{i}^{\left\langle T\right\rangle}|`, and the same for the 

550 discussion of the algorithm. 

551 

552 """ 

553 # construct the basic orbits, generators for the stabilizer chain 

554 # and transversal elements from whatever was provided 

555 transversals, basic_orbits, strong_gens_distr = \ 

556 _handle_precomputed_bsgs(base, strong_gens, transversals, 

557 basic_orbits, strong_gens_distr) 

558 base_len = len(base) 

559 degree = self.degree 

560 # size of orbit of base[pos] under the stabilizer we seek to insert 

561 # in the stabilizer chain at position pos + 1 

562 size = len(basic_orbits[pos])*len(basic_orbits[pos + 1]) \ 

563 //len(_orbit(degree, strong_gens_distr[pos], base[pos + 1])) 

564 # initialize the wanted stabilizer by a subgroup 

565 if pos + 2 > base_len - 1: 

566 T = [] 

567 else: 

568 T = strong_gens_distr[pos + 2][:] 

569 # randomized version 

570 if randomized is True: 

571 stab_pos = PermutationGroup(strong_gens_distr[pos]) 

572 schreier_vector = stab_pos.schreier_vector(base[pos + 1]) 

573 # add random elements of the stabilizer until they generate it 

574 while len(_orbit(degree, T, base[pos])) != size: 

575 new = stab_pos.random_stab(base[pos + 1], 

576 schreier_vector=schreier_vector) 

577 T.append(new) 

578 # deterministic version 

579 else: 

580 Gamma = set(basic_orbits[pos]) 

581 Gamma.remove(base[pos]) 

582 if base[pos + 1] in Gamma: 

583 Gamma.remove(base[pos + 1]) 

584 # add elements of the stabilizer until they generate it by 

585 # ruling out member of the basic orbit of base[pos] along the way 

586 while len(_orbit(degree, T, base[pos])) != size: 

587 gamma = next(iter(Gamma)) 

588 x = transversals[pos][gamma] 

589 temp = x._array_form.index(base[pos + 1]) # (~x)(base[pos + 1]) 

590 if temp not in basic_orbits[pos + 1]: 

591 Gamma = Gamma - _orbit(degree, T, gamma) 

592 else: 

593 y = transversals[pos + 1][temp] 

594 el = rmul(x, y) 

595 if el(base[pos]) not in _orbit(degree, T, base[pos]): 

596 T.append(el) 

597 Gamma = Gamma - _orbit(degree, T, base[pos]) 

598 # build the new base and strong generating set 

599 strong_gens_new_distr = strong_gens_distr[:] 

600 strong_gens_new_distr[pos + 1] = T 

601 base_new = base[:] 

602 base_new[pos], base_new[pos + 1] = base_new[pos + 1], base_new[pos] 

603 strong_gens_new = _strong_gens_from_distr(strong_gens_new_distr) 

604 for gen in T: 

605 if gen not in strong_gens_new: 

606 strong_gens_new.append(gen) 

607 return base_new, strong_gens_new 

608 

609 @property 

610 def basic_orbits(self): 

611 r""" 

612 Return the basic orbits relative to a base and strong generating set. 

613 

614 Explanation 

615 =========== 

616 

617 If `(b_1, b_2, \dots, b_k)` is a base for a group `G`, and 

618 `G^{(i)} = G_{b_1, b_2, \dots, b_{i-1}}` is the ``i``-th basic stabilizer 

619 (so that `G^{(1)} = G`), the ``i``-th basic orbit relative to this base 

620 is the orbit of `b_i` under `G^{(i)}`. See [1], pp. 87-89 for more 

621 information. 

622 

623 Examples 

624 ======== 

625 

626 >>> from sympy.combinatorics.named_groups import SymmetricGroup 

627 >>> S = SymmetricGroup(4) 

628 >>> S.basic_orbits 

629 [[0, 1, 2, 3], [1, 2, 3], [2, 3]] 

630 

631 See Also 

632 ======== 

633 

634 base, strong_gens, basic_transversals, basic_stabilizers 

635 

636 """ 

637 if self._basic_orbits == []: 

638 self.schreier_sims() 

639 return self._basic_orbits 

640 

641 @property 

642 def basic_stabilizers(self): 

643 r""" 

644 Return a chain of stabilizers relative to a base and strong generating 

645 set. 

646 

647 Explanation 

648 =========== 

649 

650 The ``i``-th basic stabilizer `G^{(i)}` relative to a base 

651 `(b_1, b_2, \dots, b_k)` is `G_{b_1, b_2, \dots, b_{i-1}}`. For more 

652 information, see [1], pp. 87-89. 

653 

654 Examples 

655 ======== 

656 

657 >>> from sympy.combinatorics.named_groups import AlternatingGroup 

658 >>> A = AlternatingGroup(4) 

659 >>> A.schreier_sims() 

660 >>> A.base 

661 [0, 1] 

662 >>> for g in A.basic_stabilizers: 

663 ... print(g) 

664 ... 

665 PermutationGroup([ 

666 (3)(0 1 2), 

667 (1 2 3)]) 

668 PermutationGroup([ 

669 (1 2 3)]) 

670 

671 See Also 

672 ======== 

673 

674 base, strong_gens, basic_orbits, basic_transversals 

675 

676 """ 

677 

678 if self._transversals == []: 

679 self.schreier_sims() 

680 strong_gens = self._strong_gens 

681 base = self._base 

682 if not base: # e.g. if self is trivial 

683 return [] 

684 strong_gens_distr = _distribute_gens_by_base(base, strong_gens) 

685 basic_stabilizers = [] 

686 for gens in strong_gens_distr: 

687 basic_stabilizers.append(PermutationGroup(gens)) 

688 return basic_stabilizers 

689 

690 @property 

691 def basic_transversals(self): 

692 """ 

693 Return basic transversals relative to a base and strong generating set. 

694 

695 Explanation 

696 =========== 

697 

698 The basic transversals are transversals of the basic orbits. They 

699 are provided as a list of dictionaries, each dictionary having 

700 keys - the elements of one of the basic orbits, and values - the 

701 corresponding transversal elements. See [1], pp. 87-89 for more 

702 information. 

703 

704 Examples 

705 ======== 

706 

707 >>> from sympy.combinatorics.named_groups import AlternatingGroup 

708 >>> A = AlternatingGroup(4) 

709 >>> A.basic_transversals 

710 [{0: (3), 1: (3)(0 1 2), 2: (3)(0 2 1), 3: (0 3 1)}, {1: (3), 2: (1 2 3), 3: (1 3 2)}] 

711 

712 See Also 

713 ======== 

714 

715 strong_gens, base, basic_orbits, basic_stabilizers 

716 

717 """ 

718 

719 if self._transversals == []: 

720 self.schreier_sims() 

721 return self._transversals 

722 

723 def composition_series(self): 

724 r""" 

725 Return the composition series for a group as a list 

726 of permutation groups. 

727 

728 Explanation 

729 =========== 

730 

731 The composition series for a group `G` is defined as a 

732 subnormal series `G = H_0 > H_1 > H_2 \ldots` A composition 

733 series is a subnormal series such that each factor group 

734 `H(i+1) / H(i)` is simple. 

735 A subnormal series is a composition series only if it is of 

736 maximum length. 

737 

738 The algorithm works as follows: 

739 Starting with the derived series the idea is to fill 

740 the gap between `G = der[i]` and `H = der[i+1]` for each 

741 `i` independently. Since, all subgroups of the abelian group 

742 `G/H` are normal so, first step is to take the generators 

743 `g` of `G` and add them to generators of `H` one by one. 

744 

745 The factor groups formed are not simple in general. Each 

746 group is obtained from the previous one by adding one 

747 generator `g`, if the previous group is denoted by `H` 

748 then the next group `K` is generated by `g` and `H`. 

749 The factor group `K/H` is cyclic and it's order is 

750 `K.order()//G.order()`. The series is then extended between 

751 `K` and `H` by groups generated by powers of `g` and `H`. 

752 The series formed is then prepended to the already existing 

753 series. 

754 

755 Examples 

756 ======== 

757 >>> from sympy.combinatorics.named_groups import SymmetricGroup 

758 >>> from sympy.combinatorics.named_groups import CyclicGroup 

759 >>> S = SymmetricGroup(12) 

760 >>> G = S.sylow_subgroup(2) 

761 >>> C = G.composition_series() 

762 >>> [H.order() for H in C] 

763 [1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1] 

764 >>> G = S.sylow_subgroup(3) 

765 >>> C = G.composition_series() 

766 >>> [H.order() for H in C] 

767 [243, 81, 27, 9, 3, 1] 

768 >>> G = CyclicGroup(12) 

769 >>> C = G.composition_series() 

770 >>> [H.order() for H in C] 

771 [12, 6, 3, 1] 

772 

773 """ 

774 der = self.derived_series() 

775 if not all(g.is_identity for g in der[-1].generators): 

776 raise NotImplementedError('Group should be solvable') 

777 series = [] 

778 

779 for i in range(len(der)-1): 

780 H = der[i+1] 

781 up_seg = [] 

782 for g in der[i].generators: 

783 K = PermutationGroup([g] + H.generators) 

784 order = K.order() // H.order() 

785 down_seg = [] 

786 for p, e in factorint(order).items(): 

787 for _ in range(e): 

788 down_seg.append(PermutationGroup([g] + H.generators)) 

789 g = g**p 

790 up_seg = down_seg + up_seg 

791 H = K 

792 up_seg[0] = der[i] 

793 series.extend(up_seg) 

794 series.append(der[-1]) 

795 return series 

796 

797 def coset_transversal(self, H): 

798 """Return a transversal of the right cosets of self by its subgroup H 

799 using the second method described in [1], Subsection 4.6.7 

800 

801 """ 

802 

803 if not H.is_subgroup(self): 

804 raise ValueError("The argument must be a subgroup") 

805 

806 if H.order() == 1: 

807 return self._elements 

808 

809 self._schreier_sims(base=H.base) # make G.base an extension of H.base 

810 

811 base = self.base 

812 base_ordering = _base_ordering(base, self.degree) 

813 identity = Permutation(self.degree - 1) 

814 

815 transversals = self.basic_transversals[:] 

816 # transversals is a list of dictionaries. Get rid of the keys 

817 # so that it is a list of lists and sort each list in 

818 # the increasing order of base[l]^x 

819 for l, t in enumerate(transversals): 

820 transversals[l] = sorted(t.values(), 

821 key = lambda x: base_ordering[base[l]^x]) 

822 

823 orbits = H.basic_orbits 

824 h_stabs = H.basic_stabilizers 

825 g_stabs = self.basic_stabilizers 

826 

827 indices = [x.order()//y.order() for x, y in zip(g_stabs, h_stabs)] 

828 

829 # T^(l) should be a right transversal of H^(l) in G^(l) for 

830 # 1<=l<=len(base). While H^(l) is the trivial group, T^(l) 

831 # contains all the elements of G^(l) so we might just as well 

832 # start with l = len(h_stabs)-1 

833 if len(g_stabs) > len(h_stabs): 

834 T = g_stabs[len(h_stabs)]._elements 

835 else: 

836 T = [identity] 

837 l = len(h_stabs)-1 

838 t_len = len(T) 

839 while l > -1: 

840 T_next = [] 

841 for u in transversals[l]: 

842 if u == identity: 

843 continue 

844 b = base_ordering[base[l]^u] 

845 for t in T: 

846 p = t*u 

847 if all(base_ordering[h^p] >= b for h in orbits[l]): 

848 T_next.append(p) 

849 if t_len + len(T_next) == indices[l]: 

850 break 

851 if t_len + len(T_next) == indices[l]: 

852 break 

853 T += T_next 

854 t_len += len(T_next) 

855 l -= 1 

856 T.remove(identity) 

857 T = [identity] + T 

858 return T 

859 

860 def _coset_representative(self, g, H): 

861 """Return the representative of Hg from the transversal that 

862 would be computed by ``self.coset_transversal(H)``. 

863 

864 """ 

865 if H.order() == 1: 

866 return g 

867 # The base of self must be an extension of H.base. 

868 if not(self.base[:len(H.base)] == H.base): 

869 self._schreier_sims(base=H.base) 

870 orbits = H.basic_orbits[:] 

871 h_transversals = [list(_.values()) for _ in H.basic_transversals] 

872 transversals = [list(_.values()) for _ in self.basic_transversals] 

873 base = self.base 

874 base_ordering = _base_ordering(base, self.degree) 

875 def step(l, x): 

876 gamma = sorted(orbits[l], key = lambda y: base_ordering[y^x])[0] 

877 i = [base[l]^h for h in h_transversals[l]].index(gamma) 

878 x = h_transversals[l][i]*x 

879 if l < len(orbits)-1: 

880 for u in transversals[l]: 

881 if base[l]^u == base[l]^x: 

882 break 

883 x = step(l+1, x*u**-1)*u 

884 return x 

885 return step(0, g) 

886 

887 def coset_table(self, H): 

888 """Return the standardised (right) coset table of self in H as 

889 a list of lists. 

890 """ 

891 # Maybe this should be made to return an instance of CosetTable 

892 # from fp_groups.py but the class would need to be changed first 

893 # to be compatible with PermutationGroups 

894 

895 if not H.is_subgroup(self): 

896 raise ValueError("The argument must be a subgroup") 

897 T = self.coset_transversal(H) 

898 n = len(T) 

899 

900 A = list(chain.from_iterable((gen, gen**-1) 

901 for gen in self.generators)) 

902 

903 table = [] 

904 for i in range(n): 

905 row = [self._coset_representative(T[i]*x, H) for x in A] 

906 row = [T.index(r) for r in row] 

907 table.append(row) 

908 

909 # standardize (this is the same as the algorithm used in coset_table) 

910 # If CosetTable is made compatible with PermutationGroups, this 

911 # should be replaced by table.standardize() 

912 A = range(len(A)) 

913 gamma = 1 

914 for alpha, a in product(range(n), A): 

915 beta = table[alpha][a] 

916 if beta >= gamma: 

917 if beta > gamma: 

918 for x in A: 

919 z = table[gamma][x] 

920 table[gamma][x] = table[beta][x] 

921 table[beta][x] = z 

922 for i in range(n): 

923 if table[i][x] == beta: 

924 table[i][x] = gamma 

925 elif table[i][x] == gamma: 

926 table[i][x] = beta 

927 gamma += 1 

928 if gamma >= n-1: 

929 return table 

930 

931 def center(self): 

932 r""" 

933 Return the center of a permutation group. 

934 

935 Explanation 

936 =========== 

937 

938 The center for a group `G` is defined as 

939 `Z(G) = \{z\in G | \forall g\in G, zg = gz \}`, 

940 the set of elements of `G` that commute with all elements of `G`. 

941 It is equal to the centralizer of `G` inside `G`, and is naturally a 

942 subgroup of `G` ([9]). 

943 

944 Examples 

945 ======== 

946 

947 >>> from sympy.combinatorics.named_groups import DihedralGroup 

948 >>> D = DihedralGroup(4) 

949 >>> G = D.center() 

950 >>> G.order() 

951 2 

952 

953 See Also 

954 ======== 

955 

956 centralizer 

957 

958 Notes 

959 ===== 

960 

961 This is a naive implementation that is a straightforward application 

962 of ``.centralizer()`` 

963 

964 """ 

965 return self.centralizer(self) 

966 

967 def centralizer(self, other): 

968 r""" 

969 Return the centralizer of a group/set/element. 

970 

971 Explanation 

972 =========== 

973 

974 The centralizer of a set of permutations ``S`` inside 

975 a group ``G`` is the set of elements of ``G`` that commute with all 

976 elements of ``S``:: 

977 

978 `C_G(S) = \{ g \in G | gs = sg \forall s \in S\}` ([10]) 

979 

980 Usually, ``S`` is a subset of ``G``, but if ``G`` is a proper subgroup of 

981 the full symmetric group, we allow for ``S`` to have elements outside 

982 ``G``. 

983 

984 It is naturally a subgroup of ``G``; the centralizer of a permutation 

985 group is equal to the centralizer of any set of generators for that 

986 group, since any element commuting with the generators commutes with 

987 any product of the generators. 

988 

989 Parameters 

990 ========== 

991 

992 other 

993 a permutation group/list of permutations/single permutation 

994 

995 Examples 

996 ======== 

997 

998 >>> from sympy.combinatorics.named_groups import (SymmetricGroup, 

999 ... CyclicGroup) 

1000 >>> S = SymmetricGroup(6) 

1001 >>> C = CyclicGroup(6) 

1002 >>> H = S.centralizer(C) 

1003 >>> H.is_subgroup(C) 

1004 True 

1005 

1006 See Also 

1007 ======== 

1008 

1009 subgroup_search 

1010 

1011 Notes 

1012 ===== 

1013 

1014 The implementation is an application of ``.subgroup_search()`` with 

1015 tests using a specific base for the group ``G``. 

1016 

1017 """ 

1018 if hasattr(other, 'generators'): 

1019 if other.is_trivial or self.is_trivial: 

1020 return self 

1021 degree = self.degree 

1022 identity = _af_new(list(range(degree))) 

1023 orbits = other.orbits() 

1024 num_orbits = len(orbits) 

1025 orbits.sort(key=lambda x: -len(x)) 

1026 long_base = [] 

1027 orbit_reps = [None]*num_orbits 

1028 orbit_reps_indices = [None]*num_orbits 

1029 orbit_descr = [None]*degree 

1030 for i in range(num_orbits): 

1031 orbit = list(orbits[i]) 

1032 orbit_reps[i] = orbit[0] 

1033 orbit_reps_indices[i] = len(long_base) 

1034 for point in orbit: 

1035 orbit_descr[point] = i 

1036 long_base = long_base + orbit 

1037 base, strong_gens = self.schreier_sims_incremental(base=long_base) 

1038 strong_gens_distr = _distribute_gens_by_base(base, strong_gens) 

1039 i = 0 

1040 for i in range(len(base)): 

1041 if strong_gens_distr[i] == [identity]: 

1042 break 

1043 base = base[:i] 

1044 base_len = i 

1045 for j in range(num_orbits): 

1046 if base[base_len - 1] in orbits[j]: 

1047 break 

1048 rel_orbits = orbits[: j + 1] 

1049 num_rel_orbits = len(rel_orbits) 

1050 transversals = [None]*num_rel_orbits 

1051 for j in range(num_rel_orbits): 

1052 rep = orbit_reps[j] 

1053 transversals[j] = dict( 

1054 other.orbit_transversal(rep, pairs=True)) 

1055 trivial_test = lambda x: True 

1056 tests = [None]*base_len 

1057 for l in range(base_len): 

1058 if base[l] in orbit_reps: 

1059 tests[l] = trivial_test 

1060 else: 

1061 def test(computed_words, l=l): 

1062 g = computed_words[l] 

1063 rep_orb_index = orbit_descr[base[l]] 

1064 rep = orbit_reps[rep_orb_index] 

1065 im = g._array_form[base[l]] 

1066 im_rep = g._array_form[rep] 

1067 tr_el = transversals[rep_orb_index][base[l]] 

1068 # using the definition of transversal, 

1069 # base[l]^g = rep^(tr_el*g); 

1070 # if g belongs to the centralizer, then 

1071 # base[l]^g = (rep^g)^tr_el 

1072 return im == tr_el._array_form[im_rep] 

1073 tests[l] = test 

1074 

1075 def prop(g): 

1076 return [rmul(g, gen) for gen in other.generators] == \ 

1077 [rmul(gen, g) for gen in other.generators] 

1078 return self.subgroup_search(prop, base=base, 

1079 strong_gens=strong_gens, tests=tests) 

1080 elif hasattr(other, '__getitem__'): 

1081 gens = list(other) 

1082 return self.centralizer(PermutationGroup(gens)) 

1083 elif hasattr(other, 'array_form'): 

1084 return self.centralizer(PermutationGroup([other])) 

1085 

1086 def commutator(self, G, H): 

1087 """ 

1088 Return the commutator of two subgroups. 

1089 

1090 Explanation 

1091 =========== 

1092 

1093 For a permutation group ``K`` and subgroups ``G``, ``H``, the 

1094 commutator of ``G`` and ``H`` is defined as the group generated 

1095 by all the commutators `[g, h] = hgh^{-1}g^{-1}` for ``g`` in ``G`` and 

1096 ``h`` in ``H``. It is naturally a subgroup of ``K`` ([1], p.27). 

1097 

1098 Examples 

1099 ======== 

1100 

1101 >>> from sympy.combinatorics.named_groups import (SymmetricGroup, 

1102 ... AlternatingGroup) 

1103 >>> S = SymmetricGroup(5) 

1104 >>> A = AlternatingGroup(5) 

1105 >>> G = S.commutator(S, A) 

1106 >>> G.is_subgroup(A) 

1107 True 

1108 

1109 See Also 

1110 ======== 

1111 

1112 derived_subgroup 

1113 

1114 Notes 

1115 ===== 

1116 

1117 The commutator of two subgroups `H, G` is equal to the normal closure 

1118 of the commutators of all the generators, i.e. `hgh^{-1}g^{-1}` for `h` 

1119 a generator of `H` and `g` a generator of `G` ([1], p.28) 

1120 

1121 """ 

1122 ggens = G.generators 

1123 hgens = H.generators 

1124 commutators = [] 

1125 for ggen in ggens: 

1126 for hgen in hgens: 

1127 commutator = rmul(hgen, ggen, ~hgen, ~ggen) 

1128 if commutator not in commutators: 

1129 commutators.append(commutator) 

1130 res = self.normal_closure(commutators) 

1131 return res 

1132 

1133 def coset_factor(self, g, factor_index=False): 

1134 """Return ``G``'s (self's) coset factorization of ``g`` 

1135 

1136 Explanation 

1137 =========== 

1138 

1139 If ``g`` is an element of ``G`` then it can be written as the product 

1140 of permutations drawn from the Schreier-Sims coset decomposition, 

1141 

1142 The permutations returned in ``f`` are those for which 

1143 the product gives ``g``: ``g = f[n]*...f[1]*f[0]`` where ``n = len(B)`` 

1144 and ``B = G.base``. f[i] is one of the permutations in 

1145 ``self._basic_orbits[i]``. 

1146 

1147 If factor_index==True, 

1148 returns a tuple ``[b[0],..,b[n]]``, where ``b[i]`` 

1149 belongs to ``self._basic_orbits[i]`` 

1150 

1151 Examples 

1152 ======== 

1153 

1154 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1155 >>> a = Permutation(0, 1, 3, 7, 6, 4)(2, 5) 

1156 >>> b = Permutation(0, 1, 3, 2)(4, 5, 7, 6) 

1157 >>> G = PermutationGroup([a, b]) 

1158 

1159 Define g: 

1160 

1161 >>> g = Permutation(7)(1, 2, 4)(3, 6, 5) 

1162 

1163 Confirm that it is an element of G: 

1164 

1165 >>> G.contains(g) 

1166 True 

1167 

1168 Thus, it can be written as a product of factors (up to 

1169 3) drawn from u. See below that a factor from u1 and u2 

1170 and the Identity permutation have been used: 

1171 

1172 >>> f = G.coset_factor(g) 

1173 >>> f[2]*f[1]*f[0] == g 

1174 True 

1175 >>> f1 = G.coset_factor(g, True); f1 

1176 [0, 4, 4] 

1177 >>> tr = G.basic_transversals 

1178 >>> f[0] == tr[0][f1[0]] 

1179 True 

1180 

1181 If g is not an element of G then [] is returned: 

1182 

1183 >>> c = Permutation(5, 6, 7) 

1184 >>> G.coset_factor(c) 

1185 [] 

1186 

1187 See Also 

1188 ======== 

1189 

1190 sympy.combinatorics.util._strip 

1191 

1192 """ 

1193 if isinstance(g, (Cycle, Permutation)): 

1194 g = g.list() 

1195 if len(g) != self._degree: 

1196 # this could either adjust the size or return [] immediately 

1197 # but we don't choose between the two and just signal a possible 

1198 # error 

1199 raise ValueError('g should be the same size as permutations of G') 

1200 I = list(range(self._degree)) 

1201 basic_orbits = self.basic_orbits 

1202 transversals = self._transversals 

1203 factors = [] 

1204 base = self.base 

1205 h = g 

1206 for i in range(len(base)): 

1207 beta = h[base[i]] 

1208 if beta == base[i]: 

1209 factors.append(beta) 

1210 continue 

1211 if beta not in basic_orbits[i]: 

1212 return [] 

1213 u = transversals[i][beta]._array_form 

1214 h = _af_rmul(_af_invert(u), h) 

1215 factors.append(beta) 

1216 if h != I: 

1217 return [] 

1218 if factor_index: 

1219 return factors 

1220 tr = self.basic_transversals 

1221 factors = [tr[i][factors[i]] for i in range(len(base))] 

1222 return factors 

1223 

1224 def generator_product(self, g, original=False): 

1225 r''' 

1226 Return a list of strong generators `[s1, \dots, sn]` 

1227 s.t `g = sn \times \dots \times s1`. If ``original=True``, make the 

1228 list contain only the original group generators 

1229 

1230 ''' 

1231 product = [] 

1232 if g.is_identity: 

1233 return [] 

1234 if g in self.strong_gens: 

1235 if not original or g in self.generators: 

1236 return [g] 

1237 else: 

1238 slp = self._strong_gens_slp[g] 

1239 for s in slp: 

1240 product.extend(self.generator_product(s, original=True)) 

1241 return product 

1242 elif g**-1 in self.strong_gens: 

1243 g = g**-1 

1244 if not original or g in self.generators: 

1245 return [g**-1] 

1246 else: 

1247 slp = self._strong_gens_slp[g] 

1248 for s in slp: 

1249 product.extend(self.generator_product(s, original=True)) 

1250 l = len(product) 

1251 product = [product[l-i-1]**-1 for i in range(l)] 

1252 return product 

1253 

1254 f = self.coset_factor(g, True) 

1255 for i, j in enumerate(f): 

1256 slp = self._transversal_slp[i][j] 

1257 for s in slp: 

1258 if not original: 

1259 product.append(self.strong_gens[s]) 

1260 else: 

1261 s = self.strong_gens[s] 

1262 product.extend(self.generator_product(s, original=True)) 

1263 return product 

1264 

1265 def coset_rank(self, g): 

1266 """rank using Schreier-Sims representation. 

1267 

1268 Explanation 

1269 =========== 

1270 

1271 The coset rank of ``g`` is the ordering number in which 

1272 it appears in the lexicographic listing according to the 

1273 coset decomposition 

1274 

1275 The ordering is the same as in G.generate(method='coset'). 

1276 If ``g`` does not belong to the group it returns None. 

1277 

1278 Examples 

1279 ======== 

1280 

1281 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1282 >>> a = Permutation(0, 1, 3, 7, 6, 4)(2, 5) 

1283 >>> b = Permutation(0, 1, 3, 2)(4, 5, 7, 6) 

1284 >>> G = PermutationGroup([a, b]) 

1285 >>> c = Permutation(7)(2, 4)(3, 5) 

1286 >>> G.coset_rank(c) 

1287 16 

1288 >>> G.coset_unrank(16) 

1289 (7)(2 4)(3 5) 

1290 

1291 See Also 

1292 ======== 

1293 

1294 coset_factor 

1295 

1296 """ 

1297 factors = self.coset_factor(g, True) 

1298 if not factors: 

1299 return None 

1300 rank = 0 

1301 b = 1 

1302 transversals = self._transversals 

1303 base = self._base 

1304 basic_orbits = self._basic_orbits 

1305 for i in range(len(base)): 

1306 k = factors[i] 

1307 j = basic_orbits[i].index(k) 

1308 rank += b*j 

1309 b = b*len(transversals[i]) 

1310 return rank 

1311 

1312 def coset_unrank(self, rank, af=False): 

1313 """unrank using Schreier-Sims representation 

1314 

1315 coset_unrank is the inverse operation of coset_rank 

1316 if 0 <= rank < order; otherwise it returns None. 

1317 

1318 """ 

1319 if rank < 0 or rank >= self.order(): 

1320 return None 

1321 base = self.base 

1322 transversals = self.basic_transversals 

1323 basic_orbits = self.basic_orbits 

1324 m = len(base) 

1325 v = [0]*m 

1326 for i in range(m): 

1327 rank, c = divmod(rank, len(transversals[i])) 

1328 v[i] = basic_orbits[i][c] 

1329 a = [transversals[i][v[i]]._array_form for i in range(m)] 

1330 h = _af_rmuln(*a) 

1331 if af: 

1332 return h 

1333 else: 

1334 return _af_new(h) 

1335 

1336 @property 

1337 def degree(self): 

1338 """Returns the size of the permutations in the group. 

1339 

1340 Explanation 

1341 =========== 

1342 

1343 The number of permutations comprising the group is given by 

1344 ``len(group)``; the number of permutations that can be generated 

1345 by the group is given by ``group.order()``. 

1346 

1347 Examples 

1348 ======== 

1349 

1350 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1351 >>> a = Permutation([1, 0, 2]) 

1352 >>> G = PermutationGroup([a]) 

1353 >>> G.degree 

1354 3 

1355 >>> len(G) 

1356 1 

1357 >>> G.order() 

1358 2 

1359 >>> list(G.generate()) 

1360 [(2), (2)(0 1)] 

1361 

1362 See Also 

1363 ======== 

1364 

1365 order 

1366 """ 

1367 return self._degree 

1368 

1369 @property 

1370 def identity(self): 

1371 ''' 

1372 Return the identity element of the permutation group. 

1373 

1374 ''' 

1375 return _af_new(list(range(self.degree))) 

1376 

1377 @property 

1378 def elements(self): 

1379 """Returns all the elements of the permutation group as a set 

1380 

1381 Examples 

1382 ======== 

1383 

1384 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1385 >>> p = PermutationGroup(Permutation(1, 3), Permutation(1, 2)) 

1386 >>> p.elements 

1387 {(1 2 3), (1 3 2), (1 3), (2 3), (3), (3)(1 2)} 

1388 

1389 """ 

1390 return set(self._elements) 

1391 

1392 @property 

1393 def _elements(self): 

1394 """Returns all the elements of the permutation group as a list 

1395 

1396 Examples 

1397 ======== 

1398 

1399 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1400 >>> p = PermutationGroup(Permutation(1, 3), Permutation(1, 2)) 

1401 >>> p._elements 

1402 [(3), (3)(1 2), (1 3), (2 3), (1 2 3), (1 3 2)] 

1403 

1404 """ 

1405 return list(islice(self.generate(), None)) 

1406 

1407 def derived_series(self): 

1408 r"""Return the derived series for the group. 

1409 

1410 Explanation 

1411 =========== 

1412 

1413 The derived series for a group `G` is defined as 

1414 `G = G_0 > G_1 > G_2 > \ldots` where `G_i = [G_{i-1}, G_{i-1}]`, 

1415 i.e. `G_i` is the derived subgroup of `G_{i-1}`, for 

1416 `i\in\mathbb{N}`. When we have `G_k = G_{k-1}` for some 

1417 `k\in\mathbb{N}`, the series terminates. 

1418 

1419 Returns 

1420 ======= 

1421 

1422 A list of permutation groups containing the members of the derived 

1423 series in the order `G = G_0, G_1, G_2, \ldots`. 

1424 

1425 Examples 

1426 ======== 

1427 

1428 >>> from sympy.combinatorics.named_groups import (SymmetricGroup, 

1429 ... AlternatingGroup, DihedralGroup) 

1430 >>> A = AlternatingGroup(5) 

1431 >>> len(A.derived_series()) 

1432 1 

1433 >>> S = SymmetricGroup(4) 

1434 >>> len(S.derived_series()) 

1435 4 

1436 >>> S.derived_series()[1].is_subgroup(AlternatingGroup(4)) 

1437 True 

1438 >>> S.derived_series()[2].is_subgroup(DihedralGroup(2)) 

1439 True 

1440 

1441 See Also 

1442 ======== 

1443 

1444 derived_subgroup 

1445 

1446 """ 

1447 res = [self] 

1448 current = self 

1449 nxt = self.derived_subgroup() 

1450 while not current.is_subgroup(nxt): 

1451 res.append(nxt) 

1452 current = nxt 

1453 nxt = nxt.derived_subgroup() 

1454 return res 

1455 

1456 def derived_subgroup(self): 

1457 r"""Compute the derived subgroup. 

1458 

1459 Explanation 

1460 =========== 

1461 

1462 The derived subgroup, or commutator subgroup is the subgroup generated 

1463 by all commutators `[g, h] = hgh^{-1}g^{-1}` for `g, h\in G` ; it is 

1464 equal to the normal closure of the set of commutators of the generators 

1465 ([1], p.28, [11]). 

1466 

1467 Examples 

1468 ======== 

1469 

1470 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1471 >>> a = Permutation([1, 0, 2, 4, 3]) 

1472 >>> b = Permutation([0, 1, 3, 2, 4]) 

1473 >>> G = PermutationGroup([a, b]) 

1474 >>> C = G.derived_subgroup() 

1475 >>> list(C.generate(af=True)) 

1476 [[0, 1, 2, 3, 4], [0, 1, 3, 4, 2], [0, 1, 4, 2, 3]] 

1477 

1478 See Also 

1479 ======== 

1480 

1481 derived_series 

1482 

1483 """ 

1484 r = self._r 

1485 gens = [p._array_form for p in self.generators] 

1486 set_commutators = set() 

1487 degree = self._degree 

1488 rng = list(range(degree)) 

1489 for i in range(r): 

1490 for j in range(r): 

1491 p1 = gens[i] 

1492 p2 = gens[j] 

1493 c = list(range(degree)) 

1494 for k in rng: 

1495 c[p2[p1[k]]] = p1[p2[k]] 

1496 ct = tuple(c) 

1497 if ct not in set_commutators: 

1498 set_commutators.add(ct) 

1499 cms = [_af_new(p) for p in set_commutators] 

1500 G2 = self.normal_closure(cms) 

1501 return G2 

1502 

1503 def generate(self, method="coset", af=False): 

1504 """Return iterator to generate the elements of the group. 

1505 

1506 Explanation 

1507 =========== 

1508 

1509 Iteration is done with one of these methods:: 

1510 

1511 method='coset' using the Schreier-Sims coset representation 

1512 method='dimino' using the Dimino method 

1513 

1514 If ``af = True`` it yields the array form of the permutations 

1515 

1516 Examples 

1517 ======== 

1518 

1519 >>> from sympy.combinatorics import PermutationGroup 

1520 >>> from sympy.combinatorics.polyhedron import tetrahedron 

1521 

1522 The permutation group given in the tetrahedron object is also 

1523 true groups: 

1524 

1525 >>> G = tetrahedron.pgroup 

1526 >>> G.is_group 

1527 True 

1528 

1529 Also the group generated by the permutations in the tetrahedron 

1530 pgroup -- even the first two -- is a proper group: 

1531 

1532 >>> H = PermutationGroup(G[0], G[1]) 

1533 >>> J = PermutationGroup(list(H.generate())); J 

1534 PermutationGroup([ 

1535 (0 1)(2 3), 

1536 (1 2 3), 

1537 (1 3 2), 

1538 (0 3 1), 

1539 (0 2 3), 

1540 (0 3)(1 2), 

1541 (0 1 3), 

1542 (3)(0 2 1), 

1543 (0 3 2), 

1544 (3)(0 1 2), 

1545 (0 2)(1 3)]) 

1546 >>> _.is_group 

1547 True 

1548 """ 

1549 if method == "coset": 

1550 return self.generate_schreier_sims(af) 

1551 elif method == "dimino": 

1552 return self.generate_dimino(af) 

1553 else: 

1554 raise NotImplementedError('No generation defined for %s' % method) 

1555 

1556 def generate_dimino(self, af=False): 

1557 """Yield group elements using Dimino's algorithm. 

1558 

1559 If ``af == True`` it yields the array form of the permutations. 

1560 

1561 Examples 

1562 ======== 

1563 

1564 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1565 >>> a = Permutation([0, 2, 1, 3]) 

1566 >>> b = Permutation([0, 2, 3, 1]) 

1567 >>> g = PermutationGroup([a, b]) 

1568 >>> list(g.generate_dimino(af=True)) 

1569 [[0, 1, 2, 3], [0, 2, 1, 3], [0, 2, 3, 1], 

1570 [0, 1, 3, 2], [0, 3, 2, 1], [0, 3, 1, 2]] 

1571 

1572 References 

1573 ========== 

1574 

1575 .. [1] The Implementation of Various Algorithms for Permutation Groups in 

1576 the Computer Algebra System: AXIOM, N.J. Doye, M.Sc. Thesis 

1577 

1578 """ 

1579 idn = list(range(self.degree)) 

1580 order = 0 

1581 element_list = [idn] 

1582 set_element_list = {tuple(idn)} 

1583 if af: 

1584 yield idn 

1585 else: 

1586 yield _af_new(idn) 

1587 gens = [p._array_form for p in self.generators] 

1588 

1589 for i in range(len(gens)): 

1590 # D elements of the subgroup G_i generated by gens[:i] 

1591 D = element_list[:] 

1592 N = [idn] 

1593 while N: 

1594 A = N 

1595 N = [] 

1596 for a in A: 

1597 for g in gens[:i + 1]: 

1598 ag = _af_rmul(a, g) 

1599 if tuple(ag) not in set_element_list: 

1600 # produce G_i*g 

1601 for d in D: 

1602 order += 1 

1603 ap = _af_rmul(d, ag) 

1604 if af: 

1605 yield ap 

1606 else: 

1607 p = _af_new(ap) 

1608 yield p 

1609 element_list.append(ap) 

1610 set_element_list.add(tuple(ap)) 

1611 N.append(ap) 

1612 self._order = len(element_list) 

1613 

1614 def generate_schreier_sims(self, af=False): 

1615 """Yield group elements using the Schreier-Sims representation 

1616 in coset_rank order 

1617 

1618 If ``af = True`` it yields the array form of the permutations 

1619 

1620 Examples 

1621 ======== 

1622 

1623 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1624 >>> a = Permutation([0, 2, 1, 3]) 

1625 >>> b = Permutation([0, 2, 3, 1]) 

1626 >>> g = PermutationGroup([a, b]) 

1627 >>> list(g.generate_schreier_sims(af=True)) 

1628 [[0, 1, 2, 3], [0, 2, 1, 3], [0, 3, 2, 1], 

1629 [0, 1, 3, 2], [0, 2, 3, 1], [0, 3, 1, 2]] 

1630 """ 

1631 

1632 n = self._degree 

1633 u = self.basic_transversals 

1634 basic_orbits = self._basic_orbits 

1635 if len(u) == 0: 

1636 for x in self.generators: 

1637 if af: 

1638 yield x._array_form 

1639 else: 

1640 yield x 

1641 return 

1642 if len(u) == 1: 

1643 for i in basic_orbits[0]: 

1644 if af: 

1645 yield u[0][i]._array_form 

1646 else: 

1647 yield u[0][i] 

1648 return 

1649 

1650 u = list(reversed(u)) 

1651 basic_orbits = basic_orbits[::-1] 

1652 # stg stack of group elements 

1653 stg = [list(range(n))] 

1654 posmax = [len(x) for x in u] 

1655 n1 = len(posmax) - 1 

1656 pos = [0]*n1 

1657 h = 0 

1658 while 1: 

1659 # backtrack when finished iterating over coset 

1660 if pos[h] >= posmax[h]: 

1661 if h == 0: 

1662 return 

1663 pos[h] = 0 

1664 h -= 1 

1665 stg.pop() 

1666 continue 

1667 p = _af_rmul(u[h][basic_orbits[h][pos[h]]]._array_form, stg[-1]) 

1668 pos[h] += 1 

1669 stg.append(p) 

1670 h += 1 

1671 if h == n1: 

1672 if af: 

1673 for i in basic_orbits[-1]: 

1674 p = _af_rmul(u[-1][i]._array_form, stg[-1]) 

1675 yield p 

1676 else: 

1677 for i in basic_orbits[-1]: 

1678 p = _af_rmul(u[-1][i]._array_form, stg[-1]) 

1679 p1 = _af_new(p) 

1680 yield p1 

1681 stg.pop() 

1682 h -= 1 

1683 

1684 @property 

1685 def generators(self): 

1686 """Returns the generators of the group. 

1687 

1688 Examples 

1689 ======== 

1690 

1691 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1692 >>> a = Permutation([0, 2, 1]) 

1693 >>> b = Permutation([1, 0, 2]) 

1694 >>> G = PermutationGroup([a, b]) 

1695 >>> G.generators 

1696 [(1 2), (2)(0 1)] 

1697 

1698 """ 

1699 return self._generators 

1700 

1701 def contains(self, g, strict=True): 

1702 """Test if permutation ``g`` belong to self, ``G``. 

1703 

1704 Explanation 

1705 =========== 

1706 

1707 If ``g`` is an element of ``G`` it can be written as a product 

1708 of factors drawn from the cosets of ``G``'s stabilizers. To see 

1709 if ``g`` is one of the actual generators defining the group use 

1710 ``G.has(g)``. 

1711 

1712 If ``strict`` is not ``True``, ``g`` will be resized, if necessary, 

1713 to match the size of permutations in ``self``. 

1714 

1715 Examples 

1716 ======== 

1717 

1718 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1719 

1720 >>> a = Permutation(1, 2) 

1721 >>> b = Permutation(2, 3, 1) 

1722 >>> G = PermutationGroup(a, b, degree=5) 

1723 >>> G.contains(G[0]) # trivial check 

1724 True 

1725 >>> elem = Permutation([[2, 3]], size=5) 

1726 >>> G.contains(elem) 

1727 True 

1728 >>> G.contains(Permutation(4)(0, 1, 2, 3)) 

1729 False 

1730 

1731 If strict is False, a permutation will be resized, if 

1732 necessary: 

1733 

1734 >>> H = PermutationGroup(Permutation(5)) 

1735 >>> H.contains(Permutation(3)) 

1736 False 

1737 >>> H.contains(Permutation(3), strict=False) 

1738 True 

1739 

1740 To test if a given permutation is present in the group: 

1741 

1742 >>> elem in G.generators 

1743 False 

1744 >>> G.has(elem) 

1745 False 

1746 

1747 See Also 

1748 ======== 

1749 

1750 coset_factor, sympy.core.basic.Basic.has, __contains__ 

1751 

1752 """ 

1753 if not isinstance(g, Permutation): 

1754 return False 

1755 if g.size != self.degree: 

1756 if strict: 

1757 return False 

1758 g = Permutation(g, size=self.degree) 

1759 if g in self.generators: 

1760 return True 

1761 return bool(self.coset_factor(g.array_form, True)) 

1762 

1763 @property 

1764 def is_perfect(self): 

1765 """Return ``True`` if the group is perfect. 

1766 A group is perfect if it equals to its derived subgroup. 

1767 

1768 Examples 

1769 ======== 

1770 

1771 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1772 >>> a = Permutation(1,2,3)(4,5) 

1773 >>> b = Permutation(1,2,3,4,5) 

1774 >>> G = PermutationGroup([a, b]) 

1775 >>> G.is_perfect 

1776 False 

1777 

1778 """ 

1779 if self._is_perfect is None: 

1780 self._is_perfect = self.equals(self.derived_subgroup()) 

1781 return self._is_perfect 

1782 

1783 @property 

1784 def is_abelian(self): 

1785 """Test if the group is Abelian. 

1786 

1787 Examples 

1788 ======== 

1789 

1790 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1791 >>> a = Permutation([0, 2, 1]) 

1792 >>> b = Permutation([1, 0, 2]) 

1793 >>> G = PermutationGroup([a, b]) 

1794 >>> G.is_abelian 

1795 False 

1796 >>> a = Permutation([0, 2, 1]) 

1797 >>> G = PermutationGroup([a]) 

1798 >>> G.is_abelian 

1799 True 

1800 

1801 """ 

1802 if self._is_abelian is not None: 

1803 return self._is_abelian 

1804 

1805 self._is_abelian = True 

1806 gens = [p._array_form for p in self.generators] 

1807 for x in gens: 

1808 for y in gens: 

1809 if y <= x: 

1810 continue 

1811 if not _af_commutes_with(x, y): 

1812 self._is_abelian = False 

1813 return False 

1814 return True 

1815 

1816 def abelian_invariants(self): 

1817 """ 

1818 Returns the abelian invariants for the given group. 

1819 Let ``G`` be a nontrivial finite abelian group. Then G is isomorphic to 

1820 the direct product of finitely many nontrivial cyclic groups of 

1821 prime-power order. 

1822 

1823 Explanation 

1824 =========== 

1825 

1826 The prime-powers that occur as the orders of the factors are uniquely 

1827 determined by G. More precisely, the primes that occur in the orders of the 

1828 factors in any such decomposition of ``G`` are exactly the primes that divide 

1829 ``|G|`` and for any such prime ``p``, if the orders of the factors that are 

1830 p-groups in one such decomposition of ``G`` are ``p^{t_1} >= p^{t_2} >= ... p^{t_r}``, 

1831 then the orders of the factors that are p-groups in any such decomposition of ``G`` 

1832 are ``p^{t_1} >= p^{t_2} >= ... p^{t_r}``. 

1833 

1834 The uniquely determined integers ``p^{t_1} >= p^{t_2} >= ... p^{t_r}``, taken 

1835 for all primes that divide ``|G|`` are called the invariants of the nontrivial 

1836 group ``G`` as suggested in ([14], p. 542). 

1837 

1838 Notes 

1839 ===== 

1840 

1841 We adopt the convention that the invariants of a trivial group are []. 

1842 

1843 Examples 

1844 ======== 

1845 

1846 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1847 >>> a = Permutation([0, 2, 1]) 

1848 >>> b = Permutation([1, 0, 2]) 

1849 >>> G = PermutationGroup([a, b]) 

1850 >>> G.abelian_invariants() 

1851 [2] 

1852 >>> from sympy.combinatorics import CyclicGroup 

1853 >>> G = CyclicGroup(7) 

1854 >>> G.abelian_invariants() 

1855 [7] 

1856 

1857 """ 

1858 if self.is_trivial: 

1859 return [] 

1860 gns = self.generators 

1861 inv = [] 

1862 G = self 

1863 H = G.derived_subgroup() 

1864 Hgens = H.generators 

1865 for p in primefactors(G.order()): 

1866 ranks = [] 

1867 while True: 

1868 pows = [] 

1869 for g in gns: 

1870 elm = g**p 

1871 if not H.contains(elm): 

1872 pows.append(elm) 

1873 K = PermutationGroup(Hgens + pows) if pows else H 

1874 r = G.order()//K.order() 

1875 G = K 

1876 gns = pows 

1877 if r == 1: 

1878 break 

1879 ranks.append(multiplicity(p, r)) 

1880 

1881 if ranks: 

1882 pows = [1]*ranks[0] 

1883 for i in ranks: 

1884 for j in range(i): 

1885 pows[j] = pows[j]*p 

1886 inv.extend(pows) 

1887 inv.sort() 

1888 return inv 

1889 

1890 def is_elementary(self, p): 

1891 """Return ``True`` if the group is elementary abelian. An elementary 

1892 abelian group is a finite abelian group, where every nontrivial 

1893 element has order `p`, where `p` is a prime. 

1894 

1895 Examples 

1896 ======== 

1897 

1898 >>> from sympy.combinatorics import Permutation, PermutationGroup 

1899 >>> a = Permutation([0, 2, 1]) 

1900 >>> G = PermutationGroup([a]) 

1901 >>> G.is_elementary(2) 

1902 True 

1903 >>> a = Permutation([0, 2, 1, 3]) 

1904 >>> b = Permutation([3, 1, 2, 0]) 

1905 >>> G = PermutationGroup([a, b]) 

1906 >>> G.is_elementary(2) 

1907 True 

1908 >>> G.is_elementary(3) 

1909 False 

1910 

1911 """ 

1912 return self.is_abelian and all(g.order() == p for g in self.generators) 

1913 

1914 def _eval_is_alt_sym_naive(self, only_sym=False, only_alt=False): 

1915 """A naive test using the group order.""" 

1916 if only_sym and only_alt: 

1917 raise ValueError( 

1918 "Both {} and {} cannot be set to True" 

1919 .format(only_sym, only_alt)) 

1920 

1921 n = self.degree 

1922 sym_order = _factorial(n) 

1923 order = self.order() 

1924 

1925 if order == sym_order: 

1926 self._is_sym = True 

1927 self._is_alt = False 

1928 if only_alt: 

1929 return False 

1930 return True 

1931 

1932 elif 2*order == sym_order: 

1933 self._is_sym = False 

1934 self._is_alt = True 

1935 if only_sym: 

1936 return False 

1937 return True 

1938 

1939 return False 

1940 

1941 def _eval_is_alt_sym_monte_carlo(self, eps=0.05, perms=None): 

1942 """A test using monte-carlo algorithm. 

1943 

1944 Parameters 

1945 ========== 

1946 

1947 eps : float, optional 

1948 The criterion for the incorrect ``False`` return. 

1949 

1950 perms : list[Permutation], optional 

1951 If explicitly given, it tests over the given candidates 

1952 for testing. 

1953 

1954 If ``None``, it randomly computes ``N_eps`` and chooses 

1955 ``N_eps`` sample of the permutation from the group. 

1956 

1957 See Also 

1958 ======== 

1959 

1960 _check_cycles_alt_sym 

1961 """ 

1962 if perms is None: 

1963 n = self.degree 

1964 if n < 17: 

1965 c_n = 0.34 

1966 else: 

1967 c_n = 0.57 

1968 d_n = (c_n*log(2))/log(n) 

1969 N_eps = int(-log(eps)/d_n) 

1970 

1971 perms = (self.random_pr() for i in range(N_eps)) 

1972 return self._eval_is_alt_sym_monte_carlo(perms=perms) 

1973 

1974 for perm in perms: 

1975 if _check_cycles_alt_sym(perm): 

1976 return True 

1977 return False 

1978 

1979 def is_alt_sym(self, eps=0.05, _random_prec=None): 

1980 r"""Monte Carlo test for the symmetric/alternating group for degrees 

1981 >= 8. 

1982 

1983 Explanation 

1984 =========== 

1985 

1986 More specifically, it is one-sided Monte Carlo with the 

1987 answer True (i.e., G is symmetric/alternating) guaranteed to be 

1988 correct, and the answer False being incorrect with probability eps. 

1989 

1990 For degree < 8, the order of the group is checked so the test 

1991 is deterministic. 

1992 

1993 Notes 

1994 ===== 

1995 

1996 The algorithm itself uses some nontrivial results from group theory and 

1997 number theory: 

1998 1) If a transitive group ``G`` of degree ``n`` contains an element 

1999 with a cycle of length ``n/2 < p < n-2`` for ``p`` a prime, ``G`` is the 

2000 symmetric or alternating group ([1], pp. 81-82) 

2001 2) The proportion of elements in the symmetric/alternating group having 

2002 the property described in 1) is approximately `\log(2)/\log(n)` 

2003 ([1], p.82; [2], pp. 226-227). 

2004 The helper function ``_check_cycles_alt_sym`` is used to 

2005 go over the cycles in a permutation and look for ones satisfying 1). 

2006 

2007 Examples 

2008 ======== 

2009 

2010 >>> from sympy.combinatorics.named_groups import DihedralGroup 

2011 >>> D = DihedralGroup(10) 

2012 >>> D.is_alt_sym() 

2013 False 

2014 

2015 See Also 

2016 ======== 

2017 

2018 _check_cycles_alt_sym 

2019 

2020 """ 

2021 if _random_prec is not None: 

2022 N_eps = _random_prec['N_eps'] 

2023 perms= (_random_prec[i] for i in range(N_eps)) 

2024 return self._eval_is_alt_sym_monte_carlo(perms=perms) 

2025 

2026 if self._is_sym or self._is_alt: 

2027 return True 

2028 if self._is_sym is False and self._is_alt is False: 

2029 return False 

2030 

2031 n = self.degree 

2032 if n < 8: 

2033 return self._eval_is_alt_sym_naive() 

2034 elif self.is_transitive(): 

2035 return self._eval_is_alt_sym_monte_carlo(eps=eps) 

2036 

2037 self._is_sym, self._is_alt = False, False 

2038 return False 

2039 

2040 @property 

2041 def is_nilpotent(self): 

2042 """Test if the group is nilpotent. 

2043 

2044 Explanation 

2045 =========== 

2046 

2047 A group `G` is nilpotent if it has a central series of finite length. 

2048 Alternatively, `G` is nilpotent if its lower central series terminates 

2049 with the trivial group. Every nilpotent group is also solvable 

2050 ([1], p.29, [12]). 

2051 

2052 Examples 

2053 ======== 

2054 

2055 >>> from sympy.combinatorics.named_groups import (SymmetricGroup, 

2056 ... CyclicGroup) 

2057 >>> C = CyclicGroup(6) 

2058 >>> C.is_nilpotent 

2059 True 

2060 >>> S = SymmetricGroup(5) 

2061 >>> S.is_nilpotent 

2062 False 

2063 

2064 See Also 

2065 ======== 

2066 

2067 lower_central_series, is_solvable 

2068 

2069 """ 

2070 if self._is_nilpotent is None: 

2071 lcs = self.lower_central_series() 

2072 terminator = lcs[len(lcs) - 1] 

2073 gens = terminator.generators 

2074 degree = self.degree 

2075 identity = _af_new(list(range(degree))) 

2076 if all(g == identity for g in gens): 

2077 self._is_solvable = True 

2078 self._is_nilpotent = True 

2079 return True 

2080 else: 

2081 self._is_nilpotent = False 

2082 return False 

2083 else: 

2084 return self._is_nilpotent 

2085 

2086 def is_normal(self, gr, strict=True): 

2087 """Test if ``G=self`` is a normal subgroup of ``gr``. 

2088 

2089 Explanation 

2090 =========== 

2091 

2092 G is normal in gr if 

2093 for each g2 in G, g1 in gr, ``g = g1*g2*g1**-1`` belongs to G 

2094 It is sufficient to check this for each g1 in gr.generators and 

2095 g2 in G.generators. 

2096 

2097 Examples 

2098 ======== 

2099 

2100 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2101 >>> a = Permutation([1, 2, 0]) 

2102 >>> b = Permutation([1, 0, 2]) 

2103 >>> G = PermutationGroup([a, b]) 

2104 >>> G1 = PermutationGroup([a, Permutation([2, 0, 1])]) 

2105 >>> G1.is_normal(G) 

2106 True 

2107 

2108 """ 

2109 if not self.is_subgroup(gr, strict=strict): 

2110 return False 

2111 d_self = self.degree 

2112 d_gr = gr.degree 

2113 if self.is_trivial and (d_self == d_gr or not strict): 

2114 return True 

2115 if self._is_abelian: 

2116 return True 

2117 new_self = self.copy() 

2118 if not strict and d_self != d_gr: 

2119 if d_self < d_gr: 

2120 new_self = PermGroup(new_self.generators + [Permutation(d_gr - 1)]) 

2121 else: 

2122 gr = PermGroup(gr.generators + [Permutation(d_self - 1)]) 

2123 gens2 = [p._array_form for p in new_self.generators] 

2124 gens1 = [p._array_form for p in gr.generators] 

2125 for g1 in gens1: 

2126 for g2 in gens2: 

2127 p = _af_rmuln(g1, g2, _af_invert(g1)) 

2128 if not new_self.coset_factor(p, True): 

2129 return False 

2130 return True 

2131 

2132 def is_primitive(self, randomized=True): 

2133 r"""Test if a group is primitive. 

2134 

2135 Explanation 

2136 =========== 

2137 

2138 A permutation group ``G`` acting on a set ``S`` is called primitive if 

2139 ``S`` contains no nontrivial block under the action of ``G`` 

2140 (a block is nontrivial if its cardinality is more than ``1``). 

2141 

2142 Notes 

2143 ===== 

2144 

2145 The algorithm is described in [1], p.83, and uses the function 

2146 minimal_block to search for blocks of the form `\{0, k\}` for ``k`` 

2147 ranging over representatives for the orbits of `G_0`, the stabilizer of 

2148 ``0``. This algorithm has complexity `O(n^2)` where ``n`` is the degree 

2149 of the group, and will perform badly if `G_0` is small. 

2150 

2151 There are two implementations offered: one finds `G_0` 

2152 deterministically using the function ``stabilizer``, and the other 

2153 (default) produces random elements of `G_0` using ``random_stab``, 

2154 hoping that they generate a subgroup of `G_0` with not too many more 

2155 orbits than `G_0` (this is suggested in [1], p.83). Behavior is changed 

2156 by the ``randomized`` flag. 

2157 

2158 Examples 

2159 ======== 

2160 

2161 >>> from sympy.combinatorics.named_groups import DihedralGroup 

2162 >>> D = DihedralGroup(10) 

2163 >>> D.is_primitive() 

2164 False 

2165 

2166 See Also 

2167 ======== 

2168 

2169 minimal_block, random_stab 

2170 

2171 """ 

2172 if self._is_primitive is not None: 

2173 return self._is_primitive 

2174 

2175 if self.is_transitive() is False: 

2176 return False 

2177 

2178 if randomized: 

2179 random_stab_gens = [] 

2180 v = self.schreier_vector(0) 

2181 for _ in range(len(self)): 

2182 random_stab_gens.append(self.random_stab(0, v)) 

2183 stab = PermutationGroup(random_stab_gens) 

2184 else: 

2185 stab = self.stabilizer(0) 

2186 orbits = stab.orbits() 

2187 for orb in orbits: 

2188 x = orb.pop() 

2189 if x != 0 and any(e != 0 for e in self.minimal_block([0, x])): 

2190 self._is_primitive = False 

2191 return False 

2192 self._is_primitive = True 

2193 return True 

2194 

2195 def minimal_blocks(self, randomized=True): 

2196 ''' 

2197 For a transitive group, return the list of all minimal 

2198 block systems. If a group is intransitive, return `False`. 

2199 

2200 Examples 

2201 ======== 

2202 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2203 >>> from sympy.combinatorics.named_groups import DihedralGroup 

2204 >>> DihedralGroup(6).minimal_blocks() 

2205 [[0, 1, 0, 1, 0, 1], [0, 1, 2, 0, 1, 2]] 

2206 >>> G = PermutationGroup(Permutation(1,2,5)) 

2207 >>> G.minimal_blocks() 

2208 False 

2209 

2210 See Also 

2211 ======== 

2212 

2213 minimal_block, is_transitive, is_primitive 

2214 

2215 ''' 

2216 def _number_blocks(blocks): 

2217 # number the blocks of a block system 

2218 # in order and return the number of 

2219 # blocks and the tuple with the 

2220 # reordering 

2221 n = len(blocks) 

2222 appeared = {} 

2223 m = 0 

2224 b = [None]*n 

2225 for i in range(n): 

2226 if blocks[i] not in appeared: 

2227 appeared[blocks[i]] = m 

2228 b[i] = m 

2229 m += 1 

2230 else: 

2231 b[i] = appeared[blocks[i]] 

2232 return tuple(b), m 

2233 

2234 if not self.is_transitive(): 

2235 return False 

2236 blocks = [] 

2237 num_blocks = [] 

2238 rep_blocks = [] 

2239 if randomized: 

2240 random_stab_gens = [] 

2241 v = self.schreier_vector(0) 

2242 for i in range(len(self)): 

2243 random_stab_gens.append(self.random_stab(0, v)) 

2244 stab = PermutationGroup(random_stab_gens) 

2245 else: 

2246 stab = self.stabilizer(0) 

2247 orbits = stab.orbits() 

2248 for orb in orbits: 

2249 x = orb.pop() 

2250 if x != 0: 

2251 block = self.minimal_block([0, x]) 

2252 num_block, _ = _number_blocks(block) 

2253 # a representative block (containing 0) 

2254 rep = {j for j in range(self.degree) if num_block[j] == 0} 

2255 # check if the system is minimal with 

2256 # respect to the already discovere ones 

2257 minimal = True 

2258 blocks_remove_mask = [False] * len(blocks) 

2259 for i, r in enumerate(rep_blocks): 

2260 if len(r) > len(rep) and rep.issubset(r): 

2261 # i-th block system is not minimal 

2262 blocks_remove_mask[i] = True 

2263 elif len(r) < len(rep) and r.issubset(rep): 

2264 # the system being checked is not minimal 

2265 minimal = False 

2266 break 

2267 # remove non-minimal representative blocks 

2268 blocks = [b for i, b in enumerate(blocks) if not blocks_remove_mask[i]] 

2269 num_blocks = [n for i, n in enumerate(num_blocks) if not blocks_remove_mask[i]] 

2270 rep_blocks = [r for i, r in enumerate(rep_blocks) if not blocks_remove_mask[i]] 

2271 

2272 if minimal and num_block not in num_blocks: 

2273 blocks.append(block) 

2274 num_blocks.append(num_block) 

2275 rep_blocks.append(rep) 

2276 return blocks 

2277 

2278 @property 

2279 def is_solvable(self): 

2280 """Test if the group is solvable. 

2281 

2282 ``G`` is solvable if its derived series terminates with the trivial 

2283 group ([1], p.29). 

2284 

2285 Examples 

2286 ======== 

2287 

2288 >>> from sympy.combinatorics.named_groups import SymmetricGroup 

2289 >>> S = SymmetricGroup(3) 

2290 >>> S.is_solvable 

2291 True 

2292 

2293 See Also 

2294 ======== 

2295 

2296 is_nilpotent, derived_series 

2297 

2298 """ 

2299 if self._is_solvable is None: 

2300 if self.order() % 2 != 0: 

2301 return True 

2302 ds = self.derived_series() 

2303 terminator = ds[len(ds) - 1] 

2304 gens = terminator.generators 

2305 degree = self.degree 

2306 identity = _af_new(list(range(degree))) 

2307 if all(g == identity for g in gens): 

2308 self._is_solvable = True 

2309 return True 

2310 else: 

2311 self._is_solvable = False 

2312 return False 

2313 else: 

2314 return self._is_solvable 

2315 

2316 def is_subgroup(self, G, strict=True): 

2317 """Return ``True`` if all elements of ``self`` belong to ``G``. 

2318 

2319 If ``strict`` is ``False`` then if ``self``'s degree is smaller 

2320 than ``G``'s, the elements will be resized to have the same degree. 

2321 

2322 Examples 

2323 ======== 

2324 

2325 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2326 >>> from sympy.combinatorics import SymmetricGroup, CyclicGroup 

2327 

2328 Testing is strict by default: the degree of each group must be the 

2329 same: 

2330 

2331 >>> p = Permutation(0, 1, 2, 3, 4, 5) 

2332 >>> G1 = PermutationGroup([Permutation(0, 1, 2), Permutation(0, 1)]) 

2333 >>> G2 = PermutationGroup([Permutation(0, 2), Permutation(0, 1, 2)]) 

2334 >>> G3 = PermutationGroup([p, p**2]) 

2335 >>> assert G1.order() == G2.order() == G3.order() == 6 

2336 >>> G1.is_subgroup(G2) 

2337 True 

2338 >>> G1.is_subgroup(G3) 

2339 False 

2340 >>> G3.is_subgroup(PermutationGroup(G3[1])) 

2341 False 

2342 >>> G3.is_subgroup(PermutationGroup(G3[0])) 

2343 True 

2344 

2345 To ignore the size, set ``strict`` to ``False``: 

2346 

2347 >>> S3 = SymmetricGroup(3) 

2348 >>> S5 = SymmetricGroup(5) 

2349 >>> S3.is_subgroup(S5, strict=False) 

2350 True 

2351 >>> C7 = CyclicGroup(7) 

2352 >>> G = S5*C7 

2353 >>> S5.is_subgroup(G, False) 

2354 True 

2355 >>> C7.is_subgroup(G, 0) 

2356 False 

2357 

2358 """ 

2359 if isinstance(G, SymmetricPermutationGroup): 

2360 if self.degree != G.degree: 

2361 return False 

2362 return True 

2363 if not isinstance(G, PermutationGroup): 

2364 return False 

2365 if self == G or self.generators[0]==Permutation(): 

2366 return True 

2367 if G.order() % self.order() != 0: 

2368 return False 

2369 if self.degree == G.degree or \ 

2370 (self.degree < G.degree and not strict): 

2371 gens = self.generators 

2372 else: 

2373 return False 

2374 return all(G.contains(g, strict=strict) for g in gens) 

2375 

2376 @property 

2377 def is_polycyclic(self): 

2378 """Return ``True`` if a group is polycyclic. A group is polycyclic if 

2379 it has a subnormal series with cyclic factors. For finite groups, 

2380 this is the same as if the group is solvable. 

2381 

2382 Examples 

2383 ======== 

2384 

2385 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2386 >>> a = Permutation([0, 2, 1, 3]) 

2387 >>> b = Permutation([2, 0, 1, 3]) 

2388 >>> G = PermutationGroup([a, b]) 

2389 >>> G.is_polycyclic 

2390 True 

2391 

2392 """ 

2393 return self.is_solvable 

2394 

2395 def is_transitive(self, strict=True): 

2396 """Test if the group is transitive. 

2397 

2398 Explanation 

2399 =========== 

2400 

2401 A group is transitive if it has a single orbit. 

2402 

2403 If ``strict`` is ``False`` the group is transitive if it has 

2404 a single orbit of length different from 1. 

2405 

2406 Examples 

2407 ======== 

2408 

2409 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2410 >>> a = Permutation([0, 2, 1, 3]) 

2411 >>> b = Permutation([2, 0, 1, 3]) 

2412 >>> G1 = PermutationGroup([a, b]) 

2413 >>> G1.is_transitive() 

2414 False 

2415 >>> G1.is_transitive(strict=False) 

2416 True 

2417 >>> c = Permutation([2, 3, 0, 1]) 

2418 >>> G2 = PermutationGroup([a, c]) 

2419 >>> G2.is_transitive() 

2420 True 

2421 >>> d = Permutation([1, 0, 2, 3]) 

2422 >>> e = Permutation([0, 1, 3, 2]) 

2423 >>> G3 = PermutationGroup([d, e]) 

2424 >>> G3.is_transitive() or G3.is_transitive(strict=False) 

2425 False 

2426 

2427 """ 

2428 if self._is_transitive: # strict or not, if True then True 

2429 return self._is_transitive 

2430 if strict: 

2431 if self._is_transitive is not None: # we only store strict=True 

2432 return self._is_transitive 

2433 

2434 ans = len(self.orbit(0)) == self.degree 

2435 self._is_transitive = ans 

2436 return ans 

2437 

2438 got_orb = False 

2439 for x in self.orbits(): 

2440 if len(x) > 1: 

2441 if got_orb: 

2442 return False 

2443 got_orb = True 

2444 return got_orb 

2445 

2446 @property 

2447 def is_trivial(self): 

2448 """Test if the group is the trivial group. 

2449 

2450 This is true if the group contains only the identity permutation. 

2451 

2452 Examples 

2453 ======== 

2454 

2455 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2456 >>> G = PermutationGroup([Permutation([0, 1, 2])]) 

2457 >>> G.is_trivial 

2458 True 

2459 

2460 """ 

2461 if self._is_trivial is None: 

2462 self._is_trivial = len(self) == 1 and self[0].is_Identity 

2463 return self._is_trivial 

2464 

2465 def lower_central_series(self): 

2466 r"""Return the lower central series for the group. 

2467 

2468 The lower central series for a group `G` is the series 

2469 `G = G_0 > G_1 > G_2 > \ldots` where 

2470 `G_k = [G, G_{k-1}]`, i.e. every term after the first is equal to the 

2471 commutator of `G` and the previous term in `G1` ([1], p.29). 

2472 

2473 Returns 

2474 ======= 

2475 

2476 A list of permutation groups in the order `G = G_0, G_1, G_2, \ldots` 

2477 

2478 Examples 

2479 ======== 

2480 

2481 >>> from sympy.combinatorics.named_groups import (AlternatingGroup, 

2482 ... DihedralGroup) 

2483 >>> A = AlternatingGroup(4) 

2484 >>> len(A.lower_central_series()) 

2485 2 

2486 >>> A.lower_central_series()[1].is_subgroup(DihedralGroup(2)) 

2487 True 

2488 

2489 See Also 

2490 ======== 

2491 

2492 commutator, derived_series 

2493 

2494 """ 

2495 res = [self] 

2496 current = self 

2497 nxt = self.commutator(self, current) 

2498 while not current.is_subgroup(nxt): 

2499 res.append(nxt) 

2500 current = nxt 

2501 nxt = self.commutator(self, current) 

2502 return res 

2503 

2504 @property 

2505 def max_div(self): 

2506 """Maximum proper divisor of the degree of a permutation group. 

2507 

2508 Explanation 

2509 =========== 

2510 

2511 Obviously, this is the degree divided by its minimal proper divisor 

2512 (larger than ``1``, if one exists). As it is guaranteed to be prime, 

2513 the ``sieve`` from ``sympy.ntheory`` is used. 

2514 This function is also used as an optimization tool for the functions 

2515 ``minimal_block`` and ``_union_find_merge``. 

2516 

2517 Examples 

2518 ======== 

2519 

2520 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2521 >>> G = PermutationGroup([Permutation([0, 2, 1, 3])]) 

2522 >>> G.max_div 

2523 2 

2524 

2525 See Also 

2526 ======== 

2527 

2528 minimal_block, _union_find_merge 

2529 

2530 """ 

2531 if self._max_div is not None: 

2532 return self._max_div 

2533 n = self.degree 

2534 if n == 1: 

2535 return 1 

2536 for x in sieve: 

2537 if n % x == 0: 

2538 d = n//x 

2539 self._max_div = d 

2540 return d 

2541 

2542 def minimal_block(self, points): 

2543 r"""For a transitive group, finds the block system generated by 

2544 ``points``. 

2545 

2546 Explanation 

2547 =========== 

2548 

2549 If a group ``G`` acts on a set ``S``, a nonempty subset ``B`` of ``S`` 

2550 is called a block under the action of ``G`` if for all ``g`` in ``G`` 

2551 we have ``gB = B`` (``g`` fixes ``B``) or ``gB`` and ``B`` have no 

2552 common points (``g`` moves ``B`` entirely). ([1], p.23; [6]). 

2553 

2554 The distinct translates ``gB`` of a block ``B`` for ``g`` in ``G`` 

2555 partition the set ``S`` and this set of translates is known as a block 

2556 system. Moreover, we obviously have that all blocks in the partition 

2557 have the same size, hence the block size divides ``|S|`` ([1], p.23). 

2558 A ``G``-congruence is an equivalence relation ``~`` on the set ``S`` 

2559 such that ``a ~ b`` implies ``g(a) ~ g(b)`` for all ``g`` in ``G``. 

2560 For a transitive group, the equivalence classes of a ``G``-congruence 

2561 and the blocks of a block system are the same thing ([1], p.23). 

2562 

2563 The algorithm below checks the group for transitivity, and then finds 

2564 the ``G``-congruence generated by the pairs ``(p_0, p_1), (p_0, p_2), 

2565 ..., (p_0,p_{k-1})`` which is the same as finding the maximal block 

2566 system (i.e., the one with minimum block size) such that 

2567 ``p_0, ..., p_{k-1}`` are in the same block ([1], p.83). 

2568 

2569 It is an implementation of Atkinson's algorithm, as suggested in [1], 

2570 and manipulates an equivalence relation on the set ``S`` using a 

2571 union-find data structure. The running time is just above 

2572 `O(|points||S|)`. ([1], pp. 83-87; [7]). 

2573 

2574 Examples 

2575 ======== 

2576 

2577 >>> from sympy.combinatorics.named_groups import DihedralGroup 

2578 >>> D = DihedralGroup(10) 

2579 >>> D.minimal_block([0, 5]) 

2580 [0, 1, 2, 3, 4, 0, 1, 2, 3, 4] 

2581 >>> D.minimal_block([0, 1]) 

2582 [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] 

2583 

2584 See Also 

2585 ======== 

2586 

2587 _union_find_rep, _union_find_merge, is_transitive, is_primitive 

2588 

2589 """ 

2590 if not self.is_transitive(): 

2591 return False 

2592 n = self.degree 

2593 gens = self.generators 

2594 # initialize the list of equivalence class representatives 

2595 parents = list(range(n)) 

2596 ranks = [1]*n 

2597 not_rep = [] 

2598 k = len(points) 

2599 # the block size must divide the degree of the group 

2600 if k > self.max_div: 

2601 return [0]*n 

2602 for i in range(k - 1): 

2603 parents[points[i + 1]] = points[0] 

2604 not_rep.append(points[i + 1]) 

2605 ranks[points[0]] = k 

2606 i = 0 

2607 len_not_rep = k - 1 

2608 while i < len_not_rep: 

2609 gamma = not_rep[i] 

2610 i += 1 

2611 for gen in gens: 

2612 # find has side effects: performs path compression on the list 

2613 # of representatives 

2614 delta = self._union_find_rep(gamma, parents) 

2615 # union has side effects: performs union by rank on the list 

2616 # of representatives 

2617 temp = self._union_find_merge(gen(gamma), gen(delta), ranks, 

2618 parents, not_rep) 

2619 if temp == -1: 

2620 return [0]*n 

2621 len_not_rep += temp 

2622 for i in range(n): 

2623 # force path compression to get the final state of the equivalence 

2624 # relation 

2625 self._union_find_rep(i, parents) 

2626 

2627 # rewrite result so that block representatives are minimal 

2628 new_reps = {} 

2629 return [new_reps.setdefault(r, i) for i, r in enumerate(parents)] 

2630 

2631 def conjugacy_class(self, x): 

2632 r"""Return the conjugacy class of an element in the group. 

2633 

2634 Explanation 

2635 =========== 

2636 

2637 The conjugacy class of an element ``g`` in a group ``G`` is the set of 

2638 elements ``x`` in ``G`` that are conjugate with ``g``, i.e. for which 

2639 

2640 ``g = xax^{-1}`` 

2641 

2642 for some ``a`` in ``G``. 

2643 

2644 Note that conjugacy is an equivalence relation, and therefore that 

2645 conjugacy classes are partitions of ``G``. For a list of all the 

2646 conjugacy classes of the group, use the conjugacy_classes() method. 

2647 

2648 In a permutation group, each conjugacy class corresponds to a particular 

2649 `cycle structure': for example, in ``S_3``, the conjugacy classes are: 

2650 

2651 * the identity class, ``{()}`` 

2652 * all transpositions, ``{(1 2), (1 3), (2 3)}`` 

2653 * all 3-cycles, ``{(1 2 3), (1 3 2)}`` 

2654 

2655 Examples 

2656 ======== 

2657 

2658 >>> from sympy.combinatorics import Permutation, SymmetricGroup 

2659 >>> S3 = SymmetricGroup(3) 

2660 >>> S3.conjugacy_class(Permutation(0, 1, 2)) 

2661 {(0 1 2), (0 2 1)} 

2662 

2663 Notes 

2664 ===== 

2665 

2666 This procedure computes the conjugacy class directly by finding the 

2667 orbit of the element under conjugation in G. This algorithm is only 

2668 feasible for permutation groups of relatively small order, but is like 

2669 the orbit() function itself in that respect. 

2670 """ 

2671 # Ref: "Computing the conjugacy classes of finite groups"; Butler, G. 

2672 # Groups '93 Galway/St Andrews; edited by Campbell, C. M. 

2673 new_class = {x} 

2674 last_iteration = new_class 

2675 

2676 while len(last_iteration) > 0: 

2677 this_iteration = set() 

2678 

2679 for y in last_iteration: 

2680 for s in self.generators: 

2681 conjugated = s * y * (~s) 

2682 if conjugated not in new_class: 

2683 this_iteration.add(conjugated) 

2684 

2685 new_class.update(last_iteration) 

2686 last_iteration = this_iteration 

2687 

2688 return new_class 

2689 

2690 

2691 def conjugacy_classes(self): 

2692 r"""Return the conjugacy classes of the group. 

2693 

2694 Explanation 

2695 =========== 

2696 

2697 As described in the documentation for the .conjugacy_class() function, 

2698 conjugacy is an equivalence relation on a group G which partitions the 

2699 set of elements. This method returns a list of all these conjugacy 

2700 classes of G. 

2701 

2702 Examples 

2703 ======== 

2704 

2705 >>> from sympy.combinatorics import SymmetricGroup 

2706 >>> SymmetricGroup(3).conjugacy_classes() 

2707 [{(2)}, {(0 1 2), (0 2 1)}, {(0 2), (1 2), (2)(0 1)}] 

2708 

2709 """ 

2710 identity = _af_new(list(range(self.degree))) 

2711 known_elements = {identity} 

2712 classes = [known_elements.copy()] 

2713 

2714 for x in self.generate(): 

2715 if x not in known_elements: 

2716 new_class = self.conjugacy_class(x) 

2717 classes.append(new_class) 

2718 known_elements.update(new_class) 

2719 

2720 return classes 

2721 

2722 def normal_closure(self, other, k=10): 

2723 r"""Return the normal closure of a subgroup/set of permutations. 

2724 

2725 Explanation 

2726 =========== 

2727 

2728 If ``S`` is a subset of a group ``G``, the normal closure of ``A`` in ``G`` 

2729 is defined as the intersection of all normal subgroups of ``G`` that 

2730 contain ``A`` ([1], p.14). Alternatively, it is the group generated by 

2731 the conjugates ``x^{-1}yx`` for ``x`` a generator of ``G`` and ``y`` a 

2732 generator of the subgroup ``\left\langle S\right\rangle`` generated by 

2733 ``S`` (for some chosen generating set for ``\left\langle S\right\rangle``) 

2734 ([1], p.73). 

2735 

2736 Parameters 

2737 ========== 

2738 

2739 other 

2740 a subgroup/list of permutations/single permutation 

2741 k 

2742 an implementation-specific parameter that determines the number 

2743 of conjugates that are adjoined to ``other`` at once 

2744 

2745 Examples 

2746 ======== 

2747 

2748 >>> from sympy.combinatorics.named_groups import (SymmetricGroup, 

2749 ... CyclicGroup, AlternatingGroup) 

2750 >>> S = SymmetricGroup(5) 

2751 >>> C = CyclicGroup(5) 

2752 >>> G = S.normal_closure(C) 

2753 >>> G.order() 

2754 60 

2755 >>> G.is_subgroup(AlternatingGroup(5)) 

2756 True 

2757 

2758 See Also 

2759 ======== 

2760 

2761 commutator, derived_subgroup, random_pr 

2762 

2763 Notes 

2764 ===== 

2765 

2766 The algorithm is described in [1], pp. 73-74; it makes use of the 

2767 generation of random elements for permutation groups by the product 

2768 replacement algorithm. 

2769 

2770 """ 

2771 if hasattr(other, 'generators'): 

2772 degree = self.degree 

2773 identity = _af_new(list(range(degree))) 

2774 

2775 if all(g == identity for g in other.generators): 

2776 return other 

2777 Z = PermutationGroup(other.generators[:]) 

2778 base, strong_gens = Z.schreier_sims_incremental() 

2779 strong_gens_distr = _distribute_gens_by_base(base, strong_gens) 

2780 basic_orbits, basic_transversals = \ 

2781 _orbits_transversals_from_bsgs(base, strong_gens_distr) 

2782 

2783 self._random_pr_init(r=10, n=20) 

2784 

2785 _loop = True 

2786 while _loop: 

2787 Z._random_pr_init(r=10, n=10) 

2788 for _ in range(k): 

2789 g = self.random_pr() 

2790 h = Z.random_pr() 

2791 conj = h^g 

2792 res = _strip(conj, base, basic_orbits, basic_transversals) 

2793 if res[0] != identity or res[1] != len(base) + 1: 

2794 gens = Z.generators 

2795 gens.append(conj) 

2796 Z = PermutationGroup(gens) 

2797 strong_gens.append(conj) 

2798 temp_base, temp_strong_gens = \ 

2799 Z.schreier_sims_incremental(base, strong_gens) 

2800 base, strong_gens = temp_base, temp_strong_gens 

2801 strong_gens_distr = \ 

2802 _distribute_gens_by_base(base, strong_gens) 

2803 basic_orbits, basic_transversals = \ 

2804 _orbits_transversals_from_bsgs(base, 

2805 strong_gens_distr) 

2806 _loop = False 

2807 for g in self.generators: 

2808 for h in Z.generators: 

2809 conj = h^g 

2810 res = _strip(conj, base, basic_orbits, 

2811 basic_transversals) 

2812 if res[0] != identity or res[1] != len(base) + 1: 

2813 _loop = True 

2814 break 

2815 if _loop: 

2816 break 

2817 return Z 

2818 elif hasattr(other, '__getitem__'): 

2819 return self.normal_closure(PermutationGroup(other)) 

2820 elif hasattr(other, 'array_form'): 

2821 return self.normal_closure(PermutationGroup([other])) 

2822 

2823 def orbit(self, alpha, action='tuples'): 

2824 r"""Compute the orbit of alpha `\{g(\alpha) | g \in G\}` as a set. 

2825 

2826 Explanation 

2827 =========== 

2828 

2829 The time complexity of the algorithm used here is `O(|Orb|*r)` where 

2830 `|Orb|` is the size of the orbit and ``r`` is the number of generators of 

2831 the group. For a more detailed analysis, see [1], p.78, [2], pp. 19-21. 

2832 Here alpha can be a single point, or a list of points. 

2833 

2834 If alpha is a single point, the ordinary orbit is computed. 

2835 if alpha is a list of points, there are three available options: 

2836 

2837 'union' - computes the union of the orbits of the points in the list 

2838 'tuples' - computes the orbit of the list interpreted as an ordered 

2839 tuple under the group action ( i.e., g((1,2,3)) = (g(1), g(2), g(3)) ) 

2840 'sets' - computes the orbit of the list interpreted as a sets 

2841 

2842 Examples 

2843 ======== 

2844 

2845 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2846 >>> a = Permutation([1, 2, 0, 4, 5, 6, 3]) 

2847 >>> G = PermutationGroup([a]) 

2848 >>> G.orbit(0) 

2849 {0, 1, 2} 

2850 >>> G.orbit([0, 4], 'union') 

2851 {0, 1, 2, 3, 4, 5, 6} 

2852 

2853 See Also 

2854 ======== 

2855 

2856 orbit_transversal 

2857 

2858 """ 

2859 return _orbit(self.degree, self.generators, alpha, action) 

2860 

2861 def orbit_rep(self, alpha, beta, schreier_vector=None): 

2862 """Return a group element which sends ``alpha`` to ``beta``. 

2863 

2864 Explanation 

2865 =========== 

2866 

2867 If ``beta`` is not in the orbit of ``alpha``, the function returns 

2868 ``False``. This implementation makes use of the schreier vector. 

2869 For a proof of correctness, see [1], p.80 

2870 

2871 Examples 

2872 ======== 

2873 

2874 >>> from sympy.combinatorics.named_groups import AlternatingGroup 

2875 >>> G = AlternatingGroup(5) 

2876 >>> G.orbit_rep(0, 4) 

2877 (0 4 1 2 3) 

2878 

2879 See Also 

2880 ======== 

2881 

2882 schreier_vector 

2883 

2884 """ 

2885 if schreier_vector is None: 

2886 schreier_vector = self.schreier_vector(alpha) 

2887 if schreier_vector[beta] is None: 

2888 return False 

2889 k = schreier_vector[beta] 

2890 gens = [x._array_form for x in self.generators] 

2891 a = [] 

2892 while k != -1: 

2893 a.append(gens[k]) 

2894 beta = gens[k].index(beta) # beta = (~gens[k])(beta) 

2895 k = schreier_vector[beta] 

2896 if a: 

2897 return _af_new(_af_rmuln(*a)) 

2898 else: 

2899 return _af_new(list(range(self._degree))) 

2900 

2901 def orbit_transversal(self, alpha, pairs=False): 

2902 r"""Computes a transversal for the orbit of ``alpha`` as a set. 

2903 

2904 Explanation 

2905 =========== 

2906 

2907 For a permutation group `G`, a transversal for the orbit 

2908 `Orb = \{g(\alpha) | g \in G\}` is a set 

2909 `\{g_\beta | g_\beta(\alpha) = \beta\}` for `\beta \in Orb`. 

2910 Note that there may be more than one possible transversal. 

2911 If ``pairs`` is set to ``True``, it returns the list of pairs 

2912 `(\beta, g_\beta)`. For a proof of correctness, see [1], p.79 

2913 

2914 Examples 

2915 ======== 

2916 

2917 >>> from sympy.combinatorics.named_groups import DihedralGroup 

2918 >>> G = DihedralGroup(6) 

2919 >>> G.orbit_transversal(0) 

2920 [(5), (0 1 2 3 4 5), (0 5)(1 4)(2 3), (0 2 4)(1 3 5), (5)(0 4)(1 3), (0 3)(1 4)(2 5)] 

2921 

2922 See Also 

2923 ======== 

2924 

2925 orbit 

2926 

2927 """ 

2928 return _orbit_transversal(self._degree, self.generators, alpha, pairs) 

2929 

2930 def orbits(self, rep=False): 

2931 """Return the orbits of ``self``, ordered according to lowest element 

2932 in each orbit. 

2933 

2934 Examples 

2935 ======== 

2936 

2937 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2938 >>> a = Permutation(1, 5)(2, 3)(4, 0, 6) 

2939 >>> b = Permutation(1, 5)(3, 4)(2, 6, 0) 

2940 >>> G = PermutationGroup([a, b]) 

2941 >>> G.orbits() 

2942 [{0, 2, 3, 4, 6}, {1, 5}] 

2943 """ 

2944 return _orbits(self._degree, self._generators) 

2945 

2946 def order(self): 

2947 """Return the order of the group: the number of permutations that 

2948 can be generated from elements of the group. 

2949 

2950 The number of permutations comprising the group is given by 

2951 ``len(group)``; the length of each permutation in the group is 

2952 given by ``group.size``. 

2953 

2954 Examples 

2955 ======== 

2956 

2957 >>> from sympy.combinatorics import Permutation, PermutationGroup 

2958 

2959 >>> a = Permutation([1, 0, 2]) 

2960 >>> G = PermutationGroup([a]) 

2961 >>> G.degree 

2962 3 

2963 >>> len(G) 

2964 1 

2965 >>> G.order() 

2966 2 

2967 >>> list(G.generate()) 

2968 [(2), (2)(0 1)] 

2969 

2970 >>> a = Permutation([0, 2, 1]) 

2971 >>> b = Permutation([1, 0, 2]) 

2972 >>> G = PermutationGroup([a, b]) 

2973 >>> G.order() 

2974 6 

2975 

2976 See Also 

2977 ======== 

2978 

2979 degree 

2980 

2981 """ 

2982 if self._order is not None: 

2983 return self._order 

2984 if self._is_sym: 

2985 n = self._degree 

2986 self._order = factorial(n) 

2987 return self._order 

2988 if self._is_alt: 

2989 n = self._degree 

2990 self._order = factorial(n)/2 

2991 return self._order 

2992 

2993 m = prod([len(x) for x in self.basic_transversals]) 

2994 self._order = m 

2995 return m 

2996 

2997 def index(self, H): 

2998 """ 

2999 Returns the index of a permutation group. 

3000 

3001 Examples 

3002 ======== 

3003 

3004 >>> from sympy.combinatorics import Permutation, PermutationGroup 

3005 >>> a = Permutation(1,2,3) 

3006 >>> b =Permutation(3) 

3007 >>> G = PermutationGroup([a]) 

3008 >>> H = PermutationGroup([b]) 

3009 >>> G.index(H) 

3010 3 

3011 

3012 """ 

3013 if H.is_subgroup(self): 

3014 return self.order()//H.order() 

3015 

3016 @property 

3017 def is_symmetric(self): 

3018 """Return ``True`` if the group is symmetric. 

3019 

3020 Examples 

3021 ======== 

3022 

3023 >>> from sympy.combinatorics import SymmetricGroup 

3024 >>> g = SymmetricGroup(5) 

3025 >>> g.is_symmetric 

3026 True 

3027 

3028 >>> from sympy.combinatorics import Permutation, PermutationGroup 

3029 >>> g = PermutationGroup( 

3030 ... Permutation(0, 1, 2, 3, 4), 

3031 ... Permutation(2, 3)) 

3032 >>> g.is_symmetric 

3033 True 

3034 

3035 Notes 

3036 ===== 

3037 

3038 This uses a naive test involving the computation of the full 

3039 group order. 

3040 If you need more quicker taxonomy for large groups, you can use 

3041 :meth:`PermutationGroup.is_alt_sym`. 

3042 However, :meth:`PermutationGroup.is_alt_sym` may not be accurate 

3043 and is not able to distinguish between an alternating group and 

3044 a symmetric group. 

3045 

3046 See Also 

3047 ======== 

3048 

3049 is_alt_sym 

3050 """ 

3051 _is_sym = self._is_sym 

3052 if _is_sym is not None: 

3053 return _is_sym 

3054 

3055 n = self.degree 

3056 if n >= 8: 

3057 if self.is_transitive(): 

3058 _is_alt_sym = self._eval_is_alt_sym_monte_carlo() 

3059 if _is_alt_sym: 

3060 if any(g.is_odd for g in self.generators): 

3061 self._is_sym, self._is_alt = True, False 

3062 return True 

3063 

3064 self._is_sym, self._is_alt = False, True 

3065 return False 

3066 

3067 return self._eval_is_alt_sym_naive(only_sym=True) 

3068 

3069 self._is_sym, self._is_alt = False, False 

3070 return False 

3071 

3072 return self._eval_is_alt_sym_naive(only_sym=True) 

3073 

3074 

3075 @property 

3076 def is_alternating(self): 

3077 """Return ``True`` if the group is alternating. 

3078 

3079 Examples 

3080 ======== 

3081 

3082 >>> from sympy.combinatorics import AlternatingGroup 

3083 >>> g = AlternatingGroup(5) 

3084 >>> g.is_alternating 

3085 True 

3086 

3087 >>> from sympy.combinatorics import Permutation, PermutationGroup 

3088 >>> g = PermutationGroup( 

3089 ... Permutation(0, 1, 2, 3, 4), 

3090 ... Permutation(2, 3, 4)) 

3091 >>> g.is_alternating 

3092 True 

3093 

3094 Notes 

3095 ===== 

3096 

3097 This uses a naive test involving the computation of the full 

3098 group order. 

3099 If you need more quicker taxonomy for large groups, you can use 

3100 :meth:`PermutationGroup.is_alt_sym`. 

3101 However, :meth:`PermutationGroup.is_alt_sym` may not be accurate 

3102 and is not able to distinguish between an alternating group and 

3103 a symmetric group. 

3104 

3105 See Also 

3106 ======== 

3107 

3108 is_alt_sym 

3109 """ 

3110 _is_alt = self._is_alt 

3111 if _is_alt is not None: 

3112 return _is_alt 

3113 

3114 n = self.degree 

3115 if n >= 8: 

3116 if self.is_transitive(): 

3117 _is_alt_sym = self._eval_is_alt_sym_monte_carlo() 

3118 if _is_alt_sym: 

3119 if all(g.is_even for g in self.generators): 

3120 self._is_sym, self._is_alt = False, True 

3121 return True 

3122 

3123 self._is_sym, self._is_alt = True, False 

3124 return False 

3125 

3126 return self._eval_is_alt_sym_naive(only_alt=True) 

3127 

3128 self._is_sym, self._is_alt = False, False 

3129 return False 

3130 

3131 return self._eval_is_alt_sym_naive(only_alt=True) 

3132 

3133 @classmethod 

3134 def _distinct_primes_lemma(cls, primes): 

3135 """Subroutine to test if there is only one cyclic group for the 

3136 order.""" 

3137 primes = sorted(primes) 

3138 l = len(primes) 

3139 for i in range(l): 

3140 for j in range(i+1, l): 

3141 if primes[j] % primes[i] == 1: 

3142 return None 

3143 return True 

3144 

3145 @property 

3146 def is_cyclic(self): 

3147 r""" 

3148 Return ``True`` if the group is Cyclic. 

3149 

3150 Examples 

3151 ======== 

3152 

3153 >>> from sympy.combinatorics.named_groups import AbelianGroup 

3154 >>> G = AbelianGroup(3, 4) 

3155 >>> G.is_cyclic 

3156 True 

3157 >>> G = AbelianGroup(4, 4) 

3158 >>> G.is_cyclic 

3159 False 

3160 

3161 Notes 

3162 ===== 

3163 

3164 If the order of a group $n$ can be factored into the distinct 

3165 primes $p_1, p_2, \dots , p_s$ and if 

3166 

3167 .. math:: 

3168 \forall i, j \in \{1, 2, \dots, s \}: 

3169 p_i \not \equiv 1 \pmod {p_j} 

3170 

3171 holds true, there is only one group of the order $n$ which 

3172 is a cyclic group [1]_. This is a generalization of the lemma 

3173 that the group of order $15, 35, \dots$ are cyclic. 

3174 

3175 And also, these additional lemmas can be used to test if a 

3176 group is cyclic if the order of the group is already found. 

3177 

3178 - If the group is abelian and the order of the group is 

3179 square-free, the group is cyclic. 

3180 - If the order of the group is less than $6$ and is not $4$, the 

3181 group is cyclic. 

3182 - If the order of the group is prime, the group is cyclic. 

3183 

3184 References 

3185 ========== 

3186 

3187 .. [1] 1978: John S. Rose: A Course on Group Theory, 

3188 Introduction to Finite Group Theory: 1.4 

3189 """ 

3190 if self._is_cyclic is not None: 

3191 return self._is_cyclic 

3192 

3193 if len(self.generators) == 1: 

3194 self._is_cyclic = True 

3195 self._is_abelian = True 

3196 return True 

3197 

3198 if self._is_abelian is False: 

3199 self._is_cyclic = False 

3200 return False 

3201 

3202 order = self.order() 

3203 

3204 if order < 6: 

3205 self._is_abelian = True 

3206 if order != 4: 

3207 self._is_cyclic = True 

3208 return True 

3209 

3210 factors = factorint(order) 

3211 if all(v == 1 for v in factors.values()): 

3212 if self._is_abelian: 

3213 self._is_cyclic = True 

3214 return True 

3215 

3216 primes = list(factors.keys()) 

3217 if PermutationGroup._distinct_primes_lemma(primes) is True: 

3218 self._is_cyclic = True 

3219 self._is_abelian = True 

3220 return True 

3221 

3222 if not self.is_abelian: 

3223 self._is_cyclic = False 

3224 return False 

3225 

3226 self._is_cyclic = all( 

3227 any(g**(order//p) != self.identity for g in self.generators) 

3228 for p, e in factors.items() if e > 1 

3229 ) 

3230 return self._is_cyclic 

3231 

3232 @property 

3233 def is_dihedral(self): 

3234 r""" 

3235 Return ``True`` if the group is dihedral. 

3236 

3237 Examples 

3238 ======== 

3239 

3240 >>> from sympy.combinatorics.perm_groups import PermutationGroup 

3241 >>> from sympy.combinatorics.permutations import Permutation 

3242 >>> from sympy.combinatorics.named_groups import SymmetricGroup, CyclicGroup 

3243 >>> G = PermutationGroup(Permutation(1, 6)(2, 5)(3, 4), Permutation(0, 1, 2, 3, 4, 5, 6)) 

3244 >>> G.is_dihedral 

3245 True 

3246 >>> G = SymmetricGroup(3) 

3247 >>> G.is_dihedral 

3248 True 

3249 >>> G = CyclicGroup(6) 

3250 >>> G.is_dihedral 

3251 False 

3252 

3253 References 

3254 ========== 

3255 

3256 .. [Di1] https://math.stackexchange.com/a/827273 

3257 .. [Di2] https://kconrad.math.uconn.edu/blurbs/grouptheory/dihedral.pdf 

3258 .. [Di3] https://kconrad.math.uconn.edu/blurbs/grouptheory/dihedral2.pdf 

3259 .. [Di4] https://en.wikipedia.org/wiki/Dihedral_group 

3260 """ 

3261 if self._is_dihedral is not None: 

3262 return self._is_dihedral 

3263 

3264 order = self.order() 

3265 

3266 if order % 2 == 1: 

3267 self._is_dihedral = False 

3268 return False 

3269 if order == 2: 

3270 self._is_dihedral = True 

3271 return True 

3272 if order == 4: 

3273 # The dihedral group of order 4 is the Klein 4-group. 

3274 self._is_dihedral = not self.is_cyclic 

3275 return self._is_dihedral 

3276 if self.is_abelian: 

3277 # The only abelian dihedral groups are the ones of orders 2 and 4. 

3278 self._is_dihedral = False 

3279 return False 

3280 

3281 # Now we know the group is of even order >= 6, and nonabelian. 

3282 n = order // 2 

3283 

3284 # Handle special cases where there are exactly two generators. 

3285 gens = self.generators 

3286 if len(gens) == 2: 

3287 x, y = gens 

3288 a, b = x.order(), y.order() 

3289 # Make a >= b 

3290 if a < b: 

3291 x, y, a, b = y, x, b, a 

3292 # Using Theorem 2.1 of [Di3]: 

3293 if a == 2 == b: 

3294 self._is_dihedral = True 

3295 return True 

3296 # Using Theorem 1.1 of [Di3]: 

3297 if a == n and b == 2 and y*x*y == ~x: 

3298 self._is_dihedral = True 

3299 return True 

3300 

3301 # Proceed with algorithm of [Di1] 

3302 # Find elements of orders 2 and n 

3303 order_2, order_n = [], [] 

3304 for p in self.elements: 

3305 k = p.order() 

3306 if k == 2: 

3307 order_2.append(p) 

3308 elif k == n: 

3309 order_n.append(p) 

3310 

3311 if len(order_2) != n + 1 - (n % 2): 

3312 self._is_dihedral = False 

3313 return False 

3314 

3315 if not order_n: 

3316 self._is_dihedral = False 

3317 return False 

3318 

3319 x = order_n[0] 

3320 # Want an element y of order 2 that is not a power of x 

3321 # (i.e. that is not the 180-deg rotation, when n is even). 

3322 y = order_2[0] 

3323 if n % 2 == 0 and y == x**(n//2): 

3324 y = order_2[1] 

3325 

3326 self._is_dihedral = (y*x*y == ~x) 

3327 return self._is_dihedral 

3328 

3329 def pointwise_stabilizer(self, points, incremental=True): 

3330 r"""Return the pointwise stabilizer for a set of points. 

3331 

3332 Explanation 

3333 =========== 

3334 

3335 For a permutation group `G` and a set of points 

3336 `\{p_1, p_2,\ldots, p_k\}`, the pointwise stabilizer of 

3337 `p_1, p_2, \ldots, p_k` is defined as 

3338 `G_{p_1,\ldots, p_k} = 

3339 \{g\in G | g(p_i) = p_i \forall i\in\{1, 2,\ldots,k\}\}` ([1],p20). 

3340 It is a subgroup of `G`. 

3341 

3342 Examples 

3343 ======== 

3344 

3345 >>> from sympy.combinatorics.named_groups import SymmetricGroup 

3346 >>> S = SymmetricGroup(7) 

3347 >>> Stab = S.pointwise_stabilizer([2, 3, 5]) 

3348 >>> Stab.is_subgroup(S.stabilizer(2).stabilizer(3).stabilizer(5)) 

3349 True 

3350 

3351 See Also 

3352 ======== 

3353 

3354 stabilizer, schreier_sims_incremental 

3355 

3356 Notes 

3357 ===== 

3358 

3359 When incremental == True, 

3360 rather than the obvious implementation using successive calls to 

3361 ``.stabilizer()``, this uses the incremental Schreier-Sims algorithm 

3362 to obtain a base with starting segment - the given points. 

3363 

3364 """ 

3365 if incremental: 

3366 base, strong_gens = self.schreier_sims_incremental(base=points) 

3367 stab_gens = [] 

3368 degree = self.degree 

3369 for gen in strong_gens: 

3370 if [gen(point) for point in points] == points: 

3371 stab_gens.append(gen) 

3372 if not stab_gens: 

3373 stab_gens = _af_new(list(range(degree))) 

3374 return PermutationGroup(stab_gens) 

3375 else: 

3376 gens = self._generators 

3377 degree = self.degree 

3378 for x in points: 

3379 gens = _stabilizer(degree, gens, x) 

3380 return PermutationGroup(gens) 

3381 

3382 def make_perm(self, n, seed=None): 

3383 """ 

3384 Multiply ``n`` randomly selected permutations from 

3385 pgroup together, starting with the identity 

3386 permutation. If ``n`` is a list of integers, those 

3387 integers will be used to select the permutations and they 

3388 will be applied in L to R order: make_perm((A, B, C)) will 

3389 give CBA(I) where I is the identity permutation. 

3390 

3391 ``seed`` is used to set the seed for the random selection 

3392 of permutations from pgroup. If this is a list of integers, 

3393 the corresponding permutations from pgroup will be selected 

3394 in the order give. This is mainly used for testing purposes. 

3395 

3396 Examples 

3397 ======== 

3398 

3399 >>> from sympy.combinatorics import Permutation, PermutationGroup 

3400 >>> a, b = [Permutation([1, 0, 3, 2]), Permutation([1, 3, 0, 2])] 

3401 >>> G = PermutationGroup([a, b]) 

3402 >>> G.make_perm(1, [0]) 

3403 (0 1)(2 3) 

3404 >>> G.make_perm(3, [0, 1, 0]) 

3405 (0 2 3 1) 

3406 >>> G.make_perm([0, 1, 0]) 

3407 (0 2 3 1) 

3408 

3409 See Also 

3410 ======== 

3411 

3412 random 

3413 """ 

3414 if is_sequence(n): 

3415 if seed is not None: 

3416 raise ValueError('If n is a sequence, seed should be None') 

3417 n, seed = len(n), n 

3418 else: 

3419 try: 

3420 n = int(n) 

3421 except TypeError: 

3422 raise ValueError('n must be an integer or a sequence.') 

3423 randomrange = _randrange(seed) 

3424 

3425 # start with the identity permutation 

3426 result = Permutation(list(range(self.degree))) 

3427 m = len(self) 

3428 for _ in range(n): 

3429 p = self[randomrange(m)] 

3430 result = rmul(result, p) 

3431 return result 

3432 

3433 def random(self, af=False): 

3434 """Return a random group element 

3435 """ 

3436 rank = randrange(self.order()) 

3437 return self.coset_unrank(rank, af) 

3438 

3439 def random_pr(self, gen_count=11, iterations=50, _random_prec=None): 

3440 """Return a random group element using product replacement. 

3441 

3442 Explanation 

3443 =========== 

3444 

3445 For the details of the product replacement algorithm, see 

3446 ``_random_pr_init`` In ``random_pr`` the actual 'product replacement' 

3447 is performed. Notice that if the attribute ``_random_gens`` 

3448 is empty, it needs to be initialized by ``_random_pr_init``. 

3449 

3450 See Also 

3451 ======== 

3452 

3453 _random_pr_init 

3454 

3455 """ 

3456 if self._random_gens == []: 

3457 self._random_pr_init(gen_count, iterations) 

3458 random_gens = self._random_gens 

3459 r = len(random_gens) - 1 

3460 

3461 # handle randomized input for testing purposes 

3462 if _random_prec is None: 

3463 s = randrange(r) 

3464 t = randrange(r - 1) 

3465 if t == s: 

3466 t = r - 1 

3467 x = choice([1, 2]) 

3468 e = choice([-1, 1]) 

3469 else: 

3470 s = _random_prec['s'] 

3471 t = _random_prec['t'] 

3472 if t == s: 

3473 t = r - 1 

3474 x = _random_prec['x'] 

3475 e = _random_prec['e'] 

3476 

3477 if x == 1: 

3478 random_gens[s] = _af_rmul(random_gens[s], _af_pow(random_gens[t], e)) 

3479 random_gens[r] = _af_rmul(random_gens[r], random_gens[s]) 

3480 else: 

3481 random_gens[s] = _af_rmul(_af_pow(random_gens[t], e), random_gens[s]) 

3482 random_gens[r] = _af_rmul(random_gens[s], random_gens[r]) 

3483 return _af_new(random_gens[r]) 

3484 

3485 def random_stab(self, alpha, schreier_vector=None, _random_prec=None): 

3486 """Random element from the stabilizer of ``alpha``. 

3487 

3488 The schreier vector for ``alpha`` is an optional argument used 

3489 for speeding up repeated calls. The algorithm is described in [1], p.81 

3490 

3491 See Also 

3492 ======== 

3493 

3494 random_pr, orbit_rep 

3495 

3496 """ 

3497 if schreier_vector is None: 

3498 schreier_vector = self.schreier_vector(alpha) 

3499 if _random_prec is None: 

3500 rand = self.random_pr() 

3501 else: 

3502 rand = _random_prec['rand'] 

3503 beta = rand(alpha) 

3504 h = self.orbit_rep(alpha, beta, schreier_vector) 

3505 return rmul(~h, rand) 

3506 

3507 def schreier_sims(self): 

3508 """Schreier-Sims algorithm. 

3509 

3510 Explanation 

3511 =========== 

3512 

3513 It computes the generators of the chain of stabilizers 

3514 `G > G_{b_1} > .. > G_{b1,..,b_r} > 1` 

3515 in which `G_{b_1,..,b_i}` stabilizes `b_1,..,b_i`, 

3516 and the corresponding ``s`` cosets. 

3517 An element of the group can be written as the product 

3518 `h_1*..*h_s`. 

3519 

3520 We use the incremental Schreier-Sims algorithm. 

3521 

3522 Examples 

3523 ======== 

3524 

3525 >>> from sympy.combinatorics import Permutation, PermutationGroup 

3526 >>> a = Permutation([0, 2, 1]) 

3527 >>> b = Permutation([1, 0, 2]) 

3528 >>> G = PermutationGroup([a, b]) 

3529 >>> G.schreier_sims() 

3530 >>> G.basic_transversals 

3531 [{0: (2)(0 1), 1: (2), 2: (1 2)}, 

3532 {0: (2), 2: (0 2)}] 

3533 """ 

3534 if self._transversals: 

3535 return 

3536 self._schreier_sims() 

3537 return 

3538 

3539 def _schreier_sims(self, base=None): 

3540 schreier = self.schreier_sims_incremental(base=base, slp_dict=True) 

3541 base, strong_gens = schreier[:2] 

3542 self._base = base 

3543 self._strong_gens = strong_gens 

3544 self._strong_gens_slp = schreier[2] 

3545 if not base: 

3546 self._transversals = [] 

3547 self._basic_orbits = [] 

3548 return 

3549 

3550 strong_gens_distr = _distribute_gens_by_base(base, strong_gens) 

3551 basic_orbits, transversals, slps = _orbits_transversals_from_bsgs(base,\ 

3552 strong_gens_distr, slp=True) 

3553 

3554 # rewrite the indices stored in slps in terms of strong_gens 

3555 for i, slp in enumerate(slps): 

3556 gens = strong_gens_distr[i] 

3557 for k in slp: 

3558 slp[k] = [strong_gens.index(gens[s]) for s in slp[k]] 

3559 

3560 self._transversals = transversals 

3561 self._basic_orbits = [sorted(x) for x in basic_orbits] 

3562 self._transversal_slp = slps 

3563 

3564 def schreier_sims_incremental(self, base=None, gens=None, slp_dict=False): 

3565 """Extend a sequence of points and generating set to a base and strong 

3566 generating set. 

3567 

3568 Parameters 

3569 ========== 

3570 

3571 base 

3572 The sequence of points to be extended to a base. Optional 

3573 parameter with default value ``[]``. 

3574 gens 

3575 The generating set to be extended to a strong generating set 

3576 relative to the base obtained. Optional parameter with default 

3577 value ``self.generators``. 

3578 

3579 slp_dict 

3580 If `True`, return a dictionary `{g: gens}` for each strong 

3581 generator `g` where `gens` is a list of strong generators 

3582 coming before `g` in `strong_gens`, such that the product 

3583 of the elements of `gens` is equal to `g`. 

3584 

3585 Returns 

3586 ======= 

3587 

3588 (base, strong_gens) 

3589 ``base`` is the base obtained, and ``strong_gens`` is the strong 

3590 generating set relative to it. The original parameters ``base``, 

3591 ``gens`` remain unchanged. 

3592 

3593 Examples 

3594 ======== 

3595 

3596 >>> from sympy.combinatorics.named_groups import AlternatingGroup 

3597 >>> from sympy.combinatorics.testutil import _verify_bsgs 

3598 >>> A = AlternatingGroup(7) 

3599 >>> base = [2, 3] 

3600 >>> seq = [2, 3] 

3601 >>> base, strong_gens = A.schreier_sims_incremental(base=seq) 

3602 >>> _verify_bsgs(A, base, strong_gens) 

3603 True 

3604 >>> base[:2] 

3605 [2, 3] 

3606 

3607 Notes 

3608 ===== 

3609 

3610 This version of the Schreier-Sims algorithm runs in polynomial time. 

3611 There are certain assumptions in the implementation - if the trivial 

3612 group is provided, ``base`` and ``gens`` are returned immediately, 

3613 as any sequence of points is a base for the trivial group. If the 

3614 identity is present in the generators ``gens``, it is removed as 

3615 it is a redundant generator. 

3616 The implementation is described in [1], pp. 90-93. 

3617 

3618 See Also 

3619 ======== 

3620 

3621 schreier_sims, schreier_sims_random 

3622 

3623 """ 

3624 if base is None: 

3625 base = [] 

3626 if gens is None: 

3627 gens = self.generators[:] 

3628 degree = self.degree 

3629 id_af = list(range(degree)) 

3630 # handle the trivial group 

3631 if len(gens) == 1 and gens[0].is_Identity: 

3632 if slp_dict: 

3633 return base, gens, {gens[0]: [gens[0]]} 

3634 return base, gens 

3635 # prevent side effects 

3636 _base, _gens = base[:], gens[:] 

3637 # remove the identity as a generator 

3638 _gens = [x for x in _gens if not x.is_Identity] 

3639 # make sure no generator fixes all base points 

3640 for gen in _gens: 

3641 if all(x == gen._array_form[x] for x in _base): 

3642 for new in id_af: 

3643 if gen._array_form[new] != new: 

3644 break 

3645 else: 

3646 assert None # can this ever happen? 

3647 _base.append(new) 

3648 # distribute generators according to basic stabilizers 

3649 strong_gens_distr = _distribute_gens_by_base(_base, _gens) 

3650 strong_gens_slp = [] 

3651 # initialize the basic stabilizers, basic orbits and basic transversals 

3652 orbs = {} 

3653 transversals = {} 

3654 slps = {} 

3655 base_len = len(_base) 

3656 for i in range(base_len): 

3657 transversals[i], slps[i] = _orbit_transversal(degree, strong_gens_distr[i], 

3658 _base[i], pairs=True, af=True, slp=True) 

3659 transversals[i] = dict(transversals[i]) 

3660 orbs[i] = list(transversals[i].keys()) 

3661 # main loop: amend the stabilizer chain until we have generators 

3662 # for all stabilizers 

3663 i = base_len - 1 

3664 while i >= 0: 

3665 # this flag is used to continue with the main loop from inside 

3666 # a nested loop 

3667 continue_i = False 

3668 # test the generators for being a strong generating set 

3669 db = {} 

3670 for beta, u_beta in list(transversals[i].items()): 

3671 for j, gen in enumerate(strong_gens_distr[i]): 

3672 gb = gen._array_form[beta] 

3673 u1 = transversals[i][gb] 

3674 g1 = _af_rmul(gen._array_form, u_beta) 

3675 slp = [(i, g) for g in slps[i][beta]] 

3676 slp = [(i, j)] + slp 

3677 if g1 != u1: 

3678 # test if the schreier generator is in the i+1-th 

3679 # would-be basic stabilizer 

3680 y = True 

3681 try: 

3682 u1_inv = db[gb] 

3683 except KeyError: 

3684 u1_inv = db[gb] = _af_invert(u1) 

3685 schreier_gen = _af_rmul(u1_inv, g1) 

3686 u1_inv_slp = slps[i][gb][:] 

3687 u1_inv_slp.reverse() 

3688 u1_inv_slp = [(i, (g,)) for g in u1_inv_slp] 

3689 slp = u1_inv_slp + slp 

3690 h, j, slp = _strip_af(schreier_gen, _base, orbs, transversals, i, slp=slp, slps=slps) 

3691 if j <= base_len: 

3692 # new strong generator h at level j 

3693 y = False 

3694 elif h: 

3695 # h fixes all base points 

3696 y = False 

3697 moved = 0 

3698 while h[moved] == moved: 

3699 moved += 1 

3700 _base.append(moved) 

3701 base_len += 1 

3702 strong_gens_distr.append([]) 

3703 if y is False: 

3704 # if a new strong generator is found, update the 

3705 # data structures and start over 

3706 h = _af_new(h) 

3707 strong_gens_slp.append((h, slp)) 

3708 for l in range(i + 1, j): 

3709 strong_gens_distr[l].append(h) 

3710 transversals[l], slps[l] =\ 

3711 _orbit_transversal(degree, strong_gens_distr[l], 

3712 _base[l], pairs=True, af=True, slp=True) 

3713 transversals[l] = dict(transversals[l]) 

3714 orbs[l] = list(transversals[l].keys()) 

3715 i = j - 1 

3716 # continue main loop using the flag 

3717 continue_i = True 

3718 if continue_i is True: 

3719 break 

3720 if continue_i is True: 

3721 break 

3722 if continue_i is True: 

3723 continue 

3724 i -= 1 

3725 

3726 strong_gens = _gens[:] 

3727 

3728 if slp_dict: 

3729 # create the list of the strong generators strong_gens and 

3730 # rewrite the indices of strong_gens_slp in terms of the 

3731 # elements of strong_gens 

3732 for k, slp in strong_gens_slp: 

3733 strong_gens.append(k) 

3734 for i in range(len(slp)): 

3735 s = slp[i] 

3736 if isinstance(s[1], tuple): 

3737 slp[i] = strong_gens_distr[s[0]][s[1][0]]**-1 

3738 else: 

3739 slp[i] = strong_gens_distr[s[0]][s[1]] 

3740 strong_gens_slp = dict(strong_gens_slp) 

3741 # add the original generators 

3742 for g in _gens: 

3743 strong_gens_slp[g] = [g] 

3744 return (_base, strong_gens, strong_gens_slp) 

3745 

3746 strong_gens.extend([k for k, _ in strong_gens_slp]) 

3747 return _base, strong_gens 

3748 

3749 def schreier_sims_random(self, base=None, gens=None, consec_succ=10, 

3750 _random_prec=None): 

3751 r"""Randomized Schreier-Sims algorithm. 

3752 

3753 Explanation 

3754 =========== 

3755 

3756 The randomized Schreier-Sims algorithm takes the sequence ``base`` 

3757 and the generating set ``gens``, and extends ``base`` to a base, and 

3758 ``gens`` to a strong generating set relative to that base with 

3759 probability of a wrong answer at most `2^{-consec\_succ}`, 

3760 provided the random generators are sufficiently random. 

3761 

3762 Parameters 

3763 ========== 

3764 

3765 base 

3766 The sequence to be extended to a base. 

3767 gens 

3768 The generating set to be extended to a strong generating set. 

3769 consec_succ 

3770 The parameter defining the probability of a wrong answer. 

3771 _random_prec 

3772 An internal parameter used for testing purposes. 

3773 

3774 Returns 

3775 ======= 

3776 

3777 (base, strong_gens) 

3778 ``base`` is the base and ``strong_gens`` is the strong generating 

3779 set relative to it. 

3780 

3781 Examples 

3782 ======== 

3783 

3784 >>> from sympy.combinatorics.testutil import _verify_bsgs 

3785 >>> from sympy.combinatorics.named_groups import SymmetricGroup 

3786 >>> S = SymmetricGroup(5) 

3787 >>> base, strong_gens = S.schreier_sims_random(consec_succ=5) 

3788 >>> _verify_bsgs(S, base, strong_gens) #doctest: +SKIP 

3789 True 

3790 

3791 Notes 

3792 ===== 

3793 

3794 The algorithm is described in detail in [1], pp. 97-98. It extends 

3795 the orbits ``orbs`` and the permutation groups ``stabs`` to 

3796 basic orbits and basic stabilizers for the base and strong generating 

3797 set produced in the end. 

3798 The idea of the extension process 

3799 is to "sift" random group elements through the stabilizer chain 

3800 and amend the stabilizers/orbits along the way when a sift 

3801 is not successful. 

3802 The helper function ``_strip`` is used to attempt 

3803 to decompose a random group element according to the current 

3804 state of the stabilizer chain and report whether the element was 

3805 fully decomposed (successful sift) or not (unsuccessful sift). In 

3806 the latter case, the level at which the sift failed is reported and 

3807 used to amend ``stabs``, ``base``, ``gens`` and ``orbs`` accordingly. 

3808 The halting condition is for ``consec_succ`` consecutive successful 

3809 sifts to pass. This makes sure that the current ``base`` and ``gens`` 

3810 form a BSGS with probability at least `1 - 1/\text{consec\_succ}`. 

3811 

3812 See Also 

3813 ======== 

3814 

3815 schreier_sims 

3816 

3817 """ 

3818 if base is None: 

3819 base = [] 

3820 if gens is None: 

3821 gens = self.generators 

3822 base_len = len(base) 

3823 n = self.degree 

3824 # make sure no generator fixes all base points 

3825 for gen in gens: 

3826 if all(gen(x) == x for x in base): 

3827 new = 0 

3828 while gen._array_form[new] == new: 

3829 new += 1 

3830 base.append(new) 

3831 base_len += 1 

3832 # distribute generators according to basic stabilizers 

3833 strong_gens_distr = _distribute_gens_by_base(base, gens) 

3834 # initialize the basic stabilizers, basic transversals and basic orbits 

3835 transversals = {} 

3836 orbs = {} 

3837 for i in range(base_len): 

3838 transversals[i] = dict(_orbit_transversal(n, strong_gens_distr[i], 

3839 base[i], pairs=True)) 

3840 orbs[i] = list(transversals[i].keys()) 

3841 # initialize the number of consecutive elements sifted 

3842 c = 0 

3843 # start sifting random elements while the number of consecutive sifts 

3844 # is less than consec_succ 

3845 while c < consec_succ: 

3846 if _random_prec is None: 

3847 g = self.random_pr() 

3848 else: 

3849 g = _random_prec['g'].pop() 

3850 h, j = _strip(g, base, orbs, transversals) 

3851 y = True 

3852 # determine whether a new base point is needed 

3853 if j <= base_len: 

3854 y = False 

3855 elif not h.is_Identity: 

3856 y = False 

3857 moved = 0 

3858 while h(moved) == moved: 

3859 moved += 1 

3860 base.append(moved) 

3861 base_len += 1 

3862 strong_gens_distr.append([]) 

3863 # if the element doesn't sift, amend the strong generators and 

3864 # associated stabilizers and orbits 

3865 if y is False: 

3866 for l in range(1, j): 

3867 strong_gens_distr[l].append(h) 

3868 transversals[l] = dict(_orbit_transversal(n, 

3869 strong_gens_distr[l], base[l], pairs=True)) 

3870 orbs[l] = list(transversals[l].keys()) 

3871 c = 0 

3872 else: 

3873 c += 1 

3874 # build the strong generating set 

3875 strong_gens = strong_gens_distr[0][:] 

3876 for gen in strong_gens_distr[1]: 

3877 if gen not in strong_gens: 

3878 strong_gens.append(gen) 

3879 return base, strong_gens 

3880 

3881 def schreier_vector(self, alpha): 

3882 """Computes the schreier vector for ``alpha``. 

3883 

3884 Explanation 

3885 =========== 

3886 

3887 The Schreier vector efficiently stores information 

3888 about the orbit of ``alpha``. It can later be used to quickly obtain 

3889 elements of the group that send ``alpha`` to a particular element 

3890 in the orbit. Notice that the Schreier vector depends on the order 

3891 in which the group generators are listed. For a definition, see [3]. 

3892 Since list indices start from zero, we adopt the convention to use 

3893 "None" instead of 0 to signify that an element does not belong 

3894 to the orbit. 

3895 For the algorithm and its correctness, see [2], pp.78-80. 

3896 

3897 Examples 

3898 ======== 

3899 

3900 >>> from sympy.combinatorics import Permutation, PermutationGroup 

3901 >>> a = Permutation([2, 4, 6, 3, 1, 5, 0]) 

3902 >>> b = Permutation([0, 1, 3, 5, 4, 6, 2]) 

3903 >>> G = PermutationGroup([a, b]) 

3904 >>> G.schreier_vector(0) 

3905 [-1, None, 0, 1, None, 1, 0] 

3906 

3907 See Also 

3908 ======== 

3909 

3910 orbit 

3911 

3912 """ 

3913 n = self.degree 

3914 v = [None]*n 

3915 v[alpha] = -1 

3916 orb = [alpha] 

3917 used = [False]*n 

3918 used[alpha] = True 

3919 gens = self.generators 

3920 r = len(gens) 

3921 for b in orb: 

3922 for i in range(r): 

3923 temp = gens[i]._array_form[b] 

3924 if used[temp] is False: 

3925 orb.append(temp) 

3926 used[temp] = True 

3927 v[temp] = i 

3928 return v 

3929 

3930 def stabilizer(self, alpha): 

3931 r"""Return the stabilizer subgroup of ``alpha``. 

3932 

3933 Explanation 

3934 =========== 

3935 

3936 The stabilizer of `\alpha` is the group `G_\alpha = 

3937 \{g \in G | g(\alpha) = \alpha\}`. 

3938 For a proof of correctness, see [1], p.79. 

3939 

3940 Examples 

3941 ======== 

3942 

3943 >>> from sympy.combinatorics.named_groups import DihedralGroup 

3944 >>> G = DihedralGroup(6) 

3945 >>> G.stabilizer(5) 

3946 PermutationGroup([ 

3947 (5)(0 4)(1 3)]) 

3948 

3949 See Also 

3950 ======== 

3951 

3952 orbit 

3953 

3954 """ 

3955 return PermGroup(_stabilizer(self._degree, self._generators, alpha)) 

3956 

3957 @property 

3958 def strong_gens(self): 

3959 r"""Return a strong generating set from the Schreier-Sims algorithm. 

3960 

3961 Explanation 

3962 =========== 

3963 

3964 A generating set `S = \{g_1, g_2, \dots, g_t\}` for a permutation group 

3965 `G` is a strong generating set relative to the sequence of points 

3966 (referred to as a "base") `(b_1, b_2, \dots, b_k)` if, for 

3967 `1 \leq i \leq k` we have that the intersection of the pointwise 

3968 stabilizer `G^{(i+1)} := G_{b_1, b_2, \dots, b_i}` with `S` generates 

3969 the pointwise stabilizer `G^{(i+1)}`. The concepts of a base and 

3970 strong generating set and their applications are discussed in depth 

3971 in [1], pp. 87-89 and [2], pp. 55-57. 

3972 

3973 Examples 

3974 ======== 

3975 

3976 >>> from sympy.combinatorics.named_groups import DihedralGroup 

3977 >>> D = DihedralGroup(4) 

3978 >>> D.strong_gens 

3979 [(0 1 2 3), (0 3)(1 2), (1 3)] 

3980 >>> D.base 

3981 [0, 1] 

3982 

3983 See Also 

3984 ======== 

3985 

3986 base, basic_transversals, basic_orbits, basic_stabilizers 

3987 

3988 """ 

3989 if self._strong_gens == []: 

3990 self.schreier_sims() 

3991 return self._strong_gens 

3992 

3993 def subgroup(self, gens): 

3994 """ 

3995 Return the subgroup generated by `gens` which is a list of 

3996 elements of the group 

3997 """ 

3998 

3999 if not all(g in self for g in gens): 

4000 raise ValueError("The group does not contain the supplied generators") 

4001 

4002 G = PermutationGroup(gens) 

4003 return G 

4004 

4005 def subgroup_search(self, prop, base=None, strong_gens=None, tests=None, 

4006 init_subgroup=None): 

4007 """Find the subgroup of all elements satisfying the property ``prop``. 

4008 

4009 Explanation 

4010 =========== 

4011 

4012 This is done by a depth-first search with respect to base images that 

4013 uses several tests to prune the search tree. 

4014 

4015 Parameters 

4016 ========== 

4017 

4018 prop 

4019 The property to be used. Has to be callable on group elements 

4020 and always return ``True`` or ``False``. It is assumed that 

4021 all group elements satisfying ``prop`` indeed form a subgroup. 

4022 base 

4023 A base for the supergroup. 

4024 strong_gens 

4025 A strong generating set for the supergroup. 

4026 tests 

4027 A list of callables of length equal to the length of ``base``. 

4028 These are used to rule out group elements by partial base images, 

4029 so that ``tests[l](g)`` returns False if the element ``g`` is known 

4030 not to satisfy prop base on where g sends the first ``l + 1`` base 

4031 points. 

4032 init_subgroup 

4033 if a subgroup of the sought group is 

4034 known in advance, it can be passed to the function as this 

4035 parameter. 

4036 

4037 Returns 

4038 ======= 

4039 

4040 res 

4041 The subgroup of all elements satisfying ``prop``. The generating 

4042 set for this group is guaranteed to be a strong generating set 

4043 relative to the base ``base``. 

4044 

4045 Examples 

4046 ======== 

4047 

4048 >>> from sympy.combinatorics.named_groups import (SymmetricGroup, 

4049 ... AlternatingGroup) 

4050 >>> from sympy.combinatorics.testutil import _verify_bsgs 

4051 >>> S = SymmetricGroup(7) 

4052 >>> prop_even = lambda x: x.is_even 

4053 >>> base, strong_gens = S.schreier_sims_incremental() 

4054 >>> G = S.subgroup_search(prop_even, base=base, strong_gens=strong_gens) 

4055 >>> G.is_subgroup(AlternatingGroup(7)) 

4056 True 

4057 >>> _verify_bsgs(G, base, G.generators) 

4058 True 

4059 

4060 Notes 

4061 ===== 

4062 

4063 This function is extremely lengthy and complicated and will require 

4064 some careful attention. The implementation is described in 

4065 [1], pp. 114-117, and the comments for the code here follow the lines 

4066 of the pseudocode in the book for clarity. 

4067 

4068 The complexity is exponential in general, since the search process by 

4069 itself visits all members of the supergroup. However, there are a lot 

4070 of tests which are used to prune the search tree, and users can define 

4071 their own tests via the ``tests`` parameter, so in practice, and for 

4072 some computations, it's not terrible. 

4073 

4074 A crucial part in the procedure is the frequent base change performed 

4075 (this is line 11 in the pseudocode) in order to obtain a new basic 

4076 stabilizer. The book mentiones that this can be done by using 

4077 ``.baseswap(...)``, however the current implementation uses a more 

4078 straightforward way to find the next basic stabilizer - calling the 

4079 function ``.stabilizer(...)`` on the previous basic stabilizer. 

4080 

4081 """ 

4082 # initialize BSGS and basic group properties 

4083 def get_reps(orbits): 

4084 # get the minimal element in the base ordering 

4085 return [min(orbit, key = lambda x: base_ordering[x]) \ 

4086 for orbit in orbits] 

4087 

4088 def update_nu(l): 

4089 temp_index = len(basic_orbits[l]) + 1 -\ 

4090 len(res_basic_orbits_init_base[l]) 

4091 # this corresponds to the element larger than all points 

4092 if temp_index >= len(sorted_orbits[l]): 

4093 nu[l] = base_ordering[degree] 

4094 else: 

4095 nu[l] = sorted_orbits[l][temp_index] 

4096 

4097 if base is None: 

4098 base, strong_gens = self.schreier_sims_incremental() 

4099 base_len = len(base) 

4100 degree = self.degree 

4101 identity = _af_new(list(range(degree))) 

4102 base_ordering = _base_ordering(base, degree) 

4103 # add an element larger than all points 

4104 base_ordering.append(degree) 

4105 # add an element smaller than all points 

4106 base_ordering.append(-1) 

4107 # compute BSGS-related structures 

4108 strong_gens_distr = _distribute_gens_by_base(base, strong_gens) 

4109 basic_orbits, transversals = _orbits_transversals_from_bsgs(base, 

4110 strong_gens_distr) 

4111 # handle subgroup initialization and tests 

4112 if init_subgroup is None: 

4113 init_subgroup = PermutationGroup([identity]) 

4114 if tests is None: 

4115 trivial_test = lambda x: True 

4116 tests = [] 

4117 for i in range(base_len): 

4118 tests.append(trivial_test) 

4119 # line 1: more initializations. 

4120 res = init_subgroup 

4121 f = base_len - 1 

4122 l = base_len - 1 

4123 # line 2: set the base for K to the base for G 

4124 res_base = base[:] 

4125 # line 3: compute BSGS and related structures for K 

4126 res_base, res_strong_gens = res.schreier_sims_incremental( 

4127 base=res_base) 

4128 res_strong_gens_distr = _distribute_gens_by_base(res_base, 

4129 res_strong_gens) 

4130 res_generators = res.generators 

4131 res_basic_orbits_init_base = \ 

4132 [_orbit(degree, res_strong_gens_distr[i], res_base[i])\ 

4133 for i in range(base_len)] 

4134 # initialize orbit representatives 

4135 orbit_reps = [None]*base_len 

4136 # line 4: orbit representatives for f-th basic stabilizer of K 

4137 orbits = _orbits(degree, res_strong_gens_distr[f]) 

4138 orbit_reps[f] = get_reps(orbits) 

4139 # line 5: remove the base point from the representatives to avoid 

4140 # getting the identity element as a generator for K 

4141 orbit_reps[f].remove(base[f]) 

4142 # line 6: more initializations 

4143 c = [0]*base_len 

4144 u = [identity]*base_len 

4145 sorted_orbits = [None]*base_len 

4146 for i in range(base_len): 

4147 sorted_orbits[i] = basic_orbits[i][:] 

4148 sorted_orbits[i].sort(key=lambda point: base_ordering[point]) 

4149 # line 7: initializations 

4150 mu = [None]*base_len 

4151 nu = [None]*base_len 

4152 # this corresponds to the element smaller than all points 

4153 mu[l] = degree + 1 

4154 update_nu(l) 

4155 # initialize computed words 

4156 computed_words = [identity]*base_len 

4157 # line 8: main loop 

4158 while True: 

4159 # apply all the tests 

4160 while l < base_len - 1 and \ 

4161 computed_words[l](base[l]) in orbit_reps[l] and \ 

4162 base_ordering[mu[l]] < \ 

4163 base_ordering[computed_words[l](base[l])] < \ 

4164 base_ordering[nu[l]] and \ 

4165 tests[l](computed_words): 

4166 # line 11: change the (partial) base of K 

4167 new_point = computed_words[l](base[l]) 

4168 res_base[l] = new_point 

4169 new_stab_gens = _stabilizer(degree, res_strong_gens_distr[l], 

4170 new_point) 

4171 res_strong_gens_distr[l + 1] = new_stab_gens 

4172 # line 12: calculate minimal orbit representatives for the 

4173 # l+1-th basic stabilizer 

4174 orbits = _orbits(degree, new_stab_gens) 

4175 orbit_reps[l + 1] = get_reps(orbits) 

4176 # line 13: amend sorted orbits 

4177 l += 1 

4178 temp_orbit = [computed_words[l - 1](point) for point 

4179 in basic_orbits[l]] 

4180 temp_orbit.sort(key=lambda point: base_ordering[point]) 

4181 sorted_orbits[l] = temp_orbit 

4182 # lines 14 and 15: update variables used minimality tests 

4183 new_mu = degree + 1 

4184 for i in range(l): 

4185 if base[l] in res_basic_orbits_init_base[i]: 

4186 candidate = computed_words[i](base[i]) 

4187 if base_ordering[candidate] > base_ordering[new_mu]: 

4188 new_mu = candidate 

4189 mu[l] = new_mu 

4190 update_nu(l) 

4191 # line 16: determine the new transversal element 

4192 c[l] = 0 

4193 temp_point = sorted_orbits[l][c[l]] 

4194 gamma = computed_words[l - 1]._array_form.index(temp_point) 

4195 u[l] = transversals[l][gamma] 

4196 # update computed words 

4197 computed_words[l] = rmul(computed_words[l - 1], u[l]) 

4198 # lines 17 & 18: apply the tests to the group element found 

4199 g = computed_words[l] 

4200 temp_point = g(base[l]) 

4201 if l == base_len - 1 and \ 

4202 base_ordering[mu[l]] < \ 

4203 base_ordering[temp_point] < base_ordering[nu[l]] and \ 

4204 temp_point in orbit_reps[l] and \ 

4205 tests[l](computed_words) and \ 

4206 prop(g): 

4207 # line 19: reset the base of K 

4208 res_generators.append(g) 

4209 res_base = base[:] 

4210 # line 20: recalculate basic orbits (and transversals) 

4211 res_strong_gens.append(g) 

4212 res_strong_gens_distr = _distribute_gens_by_base(res_base, 

4213 res_strong_gens) 

4214 res_basic_orbits_init_base = \ 

4215 [_orbit(degree, res_strong_gens_distr[i], res_base[i]) \ 

4216 for i in range(base_len)] 

4217 # line 21: recalculate orbit representatives 

4218 # line 22: reset the search depth 

4219 orbit_reps[f] = get_reps(orbits) 

4220 l = f 

4221 # line 23: go up the tree until in the first branch not fully 

4222 # searched 

4223 while l >= 0 and c[l] == len(basic_orbits[l]) - 1: 

4224 l = l - 1 

4225 # line 24: if the entire tree is traversed, return K 

4226 if l == -1: 

4227 return PermutationGroup(res_generators) 

4228 # lines 25-27: update orbit representatives 

4229 if l < f: 

4230 # line 26 

4231 f = l 

4232 c[l] = 0 

4233 # line 27 

4234 temp_orbits = _orbits(degree, res_strong_gens_distr[f]) 

4235 orbit_reps[f] = get_reps(temp_orbits) 

4236 # line 28: update variables used for minimality testing 

4237 mu[l] = degree + 1 

4238 temp_index = len(basic_orbits[l]) + 1 - \ 

4239 len(res_basic_orbits_init_base[l]) 

4240 if temp_index >= len(sorted_orbits[l]): 

4241 nu[l] = base_ordering[degree] 

4242 else: 

4243 nu[l] = sorted_orbits[l][temp_index] 

4244 # line 29: set the next element from the current branch and update 

4245 # accordingly 

4246 c[l] += 1 

4247 if l == 0: 

4248 gamma = sorted_orbits[l][c[l]] 

4249 else: 

4250 gamma = computed_words[l - 1]._array_form.index(sorted_orbits[l][c[l]]) 

4251 

4252 u[l] = transversals[l][gamma] 

4253 if l == 0: 

4254 computed_words[l] = u[l] 

4255 else: 

4256 computed_words[l] = rmul(computed_words[l - 1], u[l]) 

4257 

4258 @property 

4259 def transitivity_degree(self): 

4260 r"""Compute the degree of transitivity of the group. 

4261 

4262 Explanation 

4263 =========== 

4264 

4265 A permutation group `G` acting on `\Omega = \{0, 1, \dots, n-1\}` is 

4266 ``k``-fold transitive, if, for any `k` points 

4267 `(a_1, a_2, \dots, a_k) \in \Omega` and any `k` points 

4268 `(b_1, b_2, \dots, b_k) \in \Omega` there exists `g \in G` such that 

4269 `g(a_1) = b_1, g(a_2) = b_2, \dots, g(a_k) = b_k` 

4270 The degree of transitivity of `G` is the maximum ``k`` such that 

4271 `G` is ``k``-fold transitive. ([8]) 

4272 

4273 Examples 

4274 ======== 

4275 

4276 >>> from sympy.combinatorics import Permutation, PermutationGroup 

4277 >>> a = Permutation([1, 2, 0]) 

4278 >>> b = Permutation([1, 0, 2]) 

4279 >>> G = PermutationGroup([a, b]) 

4280 >>> G.transitivity_degree 

4281 3 

4282 

4283 See Also 

4284 ======== 

4285 

4286 is_transitive, orbit 

4287 

4288 """ 

4289 if self._transitivity_degree is None: 

4290 n = self.degree 

4291 G = self 

4292 # if G is k-transitive, a tuple (a_0,..,a_k) 

4293 # can be brought to (b_0,...,b_(k-1), b_k) 

4294 # where b_0,...,b_(k-1) are fixed points; 

4295 # consider the group G_k which stabilizes b_0,...,b_(k-1) 

4296 # if G_k is transitive on the subset excluding b_0,...,b_(k-1) 

4297 # then G is (k+1)-transitive 

4298 for i in range(n): 

4299 orb = G.orbit(i) 

4300 if len(orb) != n - i: 

4301 self._transitivity_degree = i 

4302 return i 

4303 G = G.stabilizer(i) 

4304 self._transitivity_degree = n 

4305 return n 

4306 else: 

4307 return self._transitivity_degree 

4308 

4309 def _p_elements_group(self, p): 

4310 ''' 

4311 For an abelian p-group, return the subgroup consisting of 

4312 all elements of order p (and the identity) 

4313 

4314 ''' 

4315 gens = self.generators[:] 

4316 gens = sorted(gens, key=lambda x: x.order(), reverse=True) 

4317 gens_p = [g**(g.order()/p) for g in gens] 

4318 gens_r = [] 

4319 for i in range(len(gens)): 

4320 x = gens[i] 

4321 x_order = x.order() 

4322 # x_p has order p 

4323 x_p = x**(x_order/p) 

4324 if i > 0: 

4325 P = PermutationGroup(gens_p[:i]) 

4326 else: 

4327 P = PermutationGroup(self.identity) 

4328 if x**(x_order/p) not in P: 

4329 gens_r.append(x**(x_order/p)) 

4330 else: 

4331 # replace x by an element of order (x.order()/p) 

4332 # so that gens still generates G 

4333 g = P.generator_product(x_p, original=True) 

4334 for s in g: 

4335 x = x*s**-1 

4336 x_order = x_order/p 

4337 # insert x to gens so that the sorting is preserved 

4338 del gens[i] 

4339 del gens_p[i] 

4340 j = i - 1 

4341 while j < len(gens) and gens[j].order() >= x_order: 

4342 j += 1 

4343 gens = gens[:j] + [x] + gens[j:] 

4344 gens_p = gens_p[:j] + [x] + gens_p[j:] 

4345 return PermutationGroup(gens_r) 

4346 

4347 def _sylow_alt_sym(self, p): 

4348 ''' 

4349 Return a p-Sylow subgroup of a symmetric or an 

4350 alternating group. 

4351 

4352 Explanation 

4353 =========== 

4354 

4355 The algorithm for this is hinted at in [1], Chapter 4, 

4356 Exercise 4. 

4357 

4358 For Sym(n) with n = p^i, the idea is as follows. Partition 

4359 the interval [0..n-1] into p equal parts, each of length p^(i-1): 

4360 [0..p^(i-1)-1], [p^(i-1)..2*p^(i-1)-1]...[(p-1)*p^(i-1)..p^i-1]. 

4361 Find a p-Sylow subgroup of Sym(p^(i-1)) (treated as a subgroup 

4362 of ``self``) acting on each of the parts. Call the subgroups 

4363 P_1, P_2...P_p. The generators for the subgroups P_2...P_p 

4364 can be obtained from those of P_1 by applying a "shifting" 

4365 permutation to them, that is, a permutation mapping [0..p^(i-1)-1] 

4366 to the second part (the other parts are obtained by using the shift 

4367 multiple times). The union of this permutation and the generators 

4368 of P_1 is a p-Sylow subgroup of ``self``. 

4369 

4370 For n not equal to a power of p, partition 

4371 [0..n-1] in accordance with how n would be written in base p. 

4372 E.g. for p=2 and n=11, 11 = 2^3 + 2^2 + 1 so the partition 

4373 is [[0..7], [8..9], {10}]. To generate a p-Sylow subgroup, 

4374 take the union of the generators for each of the parts. 

4375 For the above example, {(0 1), (0 2)(1 3), (0 4), (1 5)(2 7)} 

4376 from the first part, {(8 9)} from the second part and 

4377 nothing from the third. This gives 4 generators in total, and 

4378 the subgroup they generate is p-Sylow. 

4379 

4380 Alternating groups are treated the same except when p=2. In this 

4381 case, (0 1)(s s+1) should be added for an appropriate s (the start 

4382 of a part) for each part in the partitions. 

4383 

4384 See Also 

4385 ======== 

4386 

4387 sylow_subgroup, is_alt_sym 

4388 

4389 ''' 

4390 n = self.degree 

4391 gens = [] 

4392 identity = Permutation(n-1) 

4393 # the case of 2-sylow subgroups of alternating groups 

4394 # needs special treatment 

4395 alt = p == 2 and all(g.is_even for g in self.generators) 

4396 

4397 # find the presentation of n in base p 

4398 coeffs = [] 

4399 m = n 

4400 while m > 0: 

4401 coeffs.append(m % p) 

4402 m = m // p 

4403 

4404 power = len(coeffs)-1 

4405 # for a symmetric group, gens[:i] is the generating 

4406 # set for a p-Sylow subgroup on [0..p**(i-1)-1]. For 

4407 # alternating groups, the same is given by gens[:2*(i-1)] 

4408 for i in range(1, power+1): 

4409 if i == 1 and alt: 

4410 # (0 1) shouldn't be added for alternating groups 

4411 continue 

4412 gen = Permutation([(j + p**(i-1)) % p**i for j in range(p**i)]) 

4413 gens.append(identity*gen) 

4414 if alt: 

4415 gen = Permutation(0, 1)*gen*Permutation(0, 1)*gen 

4416 gens.append(gen) 

4417 

4418 # the first point in the current part (see the algorithm 

4419 # description in the docstring) 

4420 start = 0 

4421 

4422 while power > 0: 

4423 a = coeffs[power] 

4424 

4425 # make the permutation shifting the start of the first 

4426 # part ([0..p^i-1] for some i) to the current one 

4427 for _ in range(a): 

4428 shift = Permutation() 

4429 if start > 0: 

4430 for i in range(p**power): 

4431 shift = shift(i, start + i) 

4432 

4433 if alt: 

4434 gen = Permutation(0, 1)*shift*Permutation(0, 1)*shift 

4435 gens.append(gen) 

4436 j = 2*(power - 1) 

4437 else: 

4438 j = power 

4439 

4440 for i, gen in enumerate(gens[:j]): 

4441 if alt and i % 2 == 1: 

4442 continue 

4443 # shift the generator to the start of the 

4444 # partition part 

4445 gen = shift*gen*shift 

4446 gens.append(gen) 

4447 

4448 start += p**power 

4449 power = power-1 

4450 

4451 return gens 

4452 

4453 def sylow_subgroup(self, p): 

4454 ''' 

4455 Return a p-Sylow subgroup of the group. 

4456 

4457 The algorithm is described in [1], Chapter 4, Section 7 

4458 

4459 Examples 

4460 ======== 

4461 >>> from sympy.combinatorics.named_groups import DihedralGroup 

4462 >>> from sympy.combinatorics.named_groups import SymmetricGroup 

4463 >>> from sympy.combinatorics.named_groups import AlternatingGroup 

4464 

4465 >>> D = DihedralGroup(6) 

4466 >>> S = D.sylow_subgroup(2) 

4467 >>> S.order() 

4468 4 

4469 >>> G = SymmetricGroup(6) 

4470 >>> S = G.sylow_subgroup(5) 

4471 >>> S.order() 

4472 5 

4473 

4474 >>> G1 = AlternatingGroup(3) 

4475 >>> G2 = AlternatingGroup(5) 

4476 >>> G3 = AlternatingGroup(9) 

4477 

4478 >>> S1 = G1.sylow_subgroup(3) 

4479 >>> S2 = G2.sylow_subgroup(3) 

4480 >>> S3 = G3.sylow_subgroup(3) 

4481 

4482 >>> len1 = len(S1.lower_central_series()) 

4483 >>> len2 = len(S2.lower_central_series()) 

4484 >>> len3 = len(S3.lower_central_series()) 

4485 

4486 >>> len1 == len2 

4487 True 

4488 >>> len1 < len3 

4489 True 

4490 

4491 ''' 

4492 from sympy.combinatorics.homomorphisms import ( 

4493 orbit_homomorphism, block_homomorphism) 

4494 

4495 if not isprime(p): 

4496 raise ValueError("p must be a prime") 

4497 

4498 def is_p_group(G): 

4499 # check if the order of G is a power of p 

4500 # and return the power 

4501 m = G.order() 

4502 n = 0 

4503 while m % p == 0: 

4504 m = m/p 

4505 n += 1 

4506 if m == 1: 

4507 return True, n 

4508 return False, n 

4509 

4510 def _sylow_reduce(mu, nu): 

4511 # reduction based on two homomorphisms 

4512 # mu and nu with trivially intersecting 

4513 # kernels 

4514 Q = mu.image().sylow_subgroup(p) 

4515 Q = mu.invert_subgroup(Q) 

4516 nu = nu.restrict_to(Q) 

4517 R = nu.image().sylow_subgroup(p) 

4518 return nu.invert_subgroup(R) 

4519 

4520 order = self.order() 

4521 if order % p != 0: 

4522 return PermutationGroup([self.identity]) 

4523 p_group, n = is_p_group(self) 

4524 if p_group: 

4525 return self 

4526 

4527 if self.is_alt_sym(): 

4528 return PermutationGroup(self._sylow_alt_sym(p)) 

4529 

4530 # if there is a non-trivial orbit with size not divisible 

4531 # by p, the sylow subgroup is contained in its stabilizer 

4532 # (by orbit-stabilizer theorem) 

4533 orbits = self.orbits() 

4534 non_p_orbits = [o for o in orbits if len(o) % p != 0 and len(o) != 1] 

4535 if non_p_orbits: 

4536 G = self.stabilizer(list(non_p_orbits[0]).pop()) 

4537 return G.sylow_subgroup(p) 

4538 

4539 if not self.is_transitive(): 

4540 # apply _sylow_reduce to orbit actions 

4541 orbits = sorted(orbits, key=len) 

4542 omega1 = orbits.pop() 

4543 omega2 = orbits[0].union(*orbits) 

4544 mu = orbit_homomorphism(self, omega1) 

4545 nu = orbit_homomorphism(self, omega2) 

4546 return _sylow_reduce(mu, nu) 

4547 

4548 blocks = self.minimal_blocks() 

4549 if len(blocks) > 1: 

4550 # apply _sylow_reduce to block system actions 

4551 mu = block_homomorphism(self, blocks[0]) 

4552 nu = block_homomorphism(self, blocks[1]) 

4553 return _sylow_reduce(mu, nu) 

4554 elif len(blocks) == 1: 

4555 block = list(blocks)[0] 

4556 if any(e != 0 for e in block): 

4557 # self is imprimitive 

4558 mu = block_homomorphism(self, block) 

4559 if not is_p_group(mu.image())[0]: 

4560 S = mu.image().sylow_subgroup(p) 

4561 return mu.invert_subgroup(S).sylow_subgroup(p) 

4562 

4563 # find an element of order p 

4564 g = self.random() 

4565 g_order = g.order() 

4566 while g_order % p != 0 or g_order == 0: 

4567 g = self.random() 

4568 g_order = g.order() 

4569 g = g**(g_order // p) 

4570 if order % p**2 != 0: 

4571 return PermutationGroup(g) 

4572 

4573 C = self.centralizer(g) 

4574 while C.order() % p**n != 0: 

4575 S = C.sylow_subgroup(p) 

4576 s_order = S.order() 

4577 Z = S.center() 

4578 P = Z._p_elements_group(p) 

4579 h = P.random() 

4580 C_h = self.centralizer(h) 

4581 while C_h.order() % p*s_order != 0: 

4582 h = P.random() 

4583 C_h = self.centralizer(h) 

4584 C = C_h 

4585 

4586 return C.sylow_subgroup(p) 

4587 

4588 def _block_verify(self, L, alpha): 

4589 delta = sorted(self.orbit(alpha)) 

4590 # p[i] will be the number of the block 

4591 # delta[i] belongs to 

4592 p = [-1]*len(delta) 

4593 blocks = [-1]*len(delta) 

4594 

4595 B = [[]] # future list of blocks 

4596 u = [0]*len(delta) # u[i] in L s.t. alpha^u[i] = B[0][i] 

4597 

4598 t = L.orbit_transversal(alpha, pairs=True) 

4599 for a, beta in t: 

4600 B[0].append(a) 

4601 i_a = delta.index(a) 

4602 p[i_a] = 0 

4603 blocks[i_a] = alpha 

4604 u[i_a] = beta 

4605 

4606 rho = 0 

4607 m = 0 # number of blocks - 1 

4608 

4609 while rho <= m: 

4610 beta = B[rho][0] 

4611 for g in self.generators: 

4612 d = beta^g 

4613 i_d = delta.index(d) 

4614 sigma = p[i_d] 

4615 if sigma < 0: 

4616 # define a new block 

4617 m += 1 

4618 sigma = m 

4619 u[i_d] = u[delta.index(beta)]*g 

4620 p[i_d] = sigma 

4621 rep = d 

4622 blocks[i_d] = rep 

4623 newb = [rep] 

4624 for gamma in B[rho][1:]: 

4625 i_gamma = delta.index(gamma) 

4626 d = gamma^g 

4627 i_d = delta.index(d) 

4628 if p[i_d] < 0: 

4629 u[i_d] = u[i_gamma]*g 

4630 p[i_d] = sigma 

4631 blocks[i_d] = rep 

4632 newb.append(d) 

4633 else: 

4634 # B[rho] is not a block 

4635 s = u[i_gamma]*g*u[i_d]**(-1) 

4636 return False, s 

4637 

4638 B.append(newb) 

4639 else: 

4640 for h in B[rho][1:]: 

4641 if h^g not in B[sigma]: 

4642 # B[rho] is not a block 

4643 s = u[delta.index(beta)]*g*u[i_d]**(-1) 

4644 return False, s 

4645 rho += 1 

4646 

4647 return True, blocks 

4648 

4649 def _verify(H, K, phi, z, alpha): 

4650 ''' 

4651 Return a list of relators ``rels`` in generators ``gens`_h` that 

4652 are mapped to ``H.generators`` by ``phi`` so that given a finite 

4653 presentation <gens_k | rels_k> of ``K`` on a subset of ``gens_h`` 

4654 <gens_h | rels_k + rels> is a finite presentation of ``H``. 

4655 

4656 Explanation 

4657 =========== 

4658 

4659 ``H`` should be generated by the union of ``K.generators`` and ``z`` 

4660 (a single generator), and ``H.stabilizer(alpha) == K``; ``phi`` is a 

4661 canonical injection from a free group into a permutation group 

4662 containing ``H``. 

4663 

4664 The algorithm is described in [1], Chapter 6. 

4665 

4666 Examples 

4667 ======== 

4668 

4669 >>> from sympy.combinatorics import free_group, Permutation, PermutationGroup 

4670 >>> from sympy.combinatorics.homomorphisms import homomorphism 

4671 >>> from sympy.combinatorics.fp_groups import FpGroup 

4672 

4673 >>> H = PermutationGroup(Permutation(0, 2), Permutation (1, 5)) 

4674 >>> K = PermutationGroup(Permutation(5)(0, 2)) 

4675 >>> F = free_group("x_0 x_1")[0] 

4676 >>> gens = F.generators 

4677 >>> phi = homomorphism(F, H, F.generators, H.generators) 

4678 >>> rels_k = [gens[0]**2] # relators for presentation of K 

4679 >>> z= Permutation(1, 5) 

4680 >>> check, rels_h = H._verify(K, phi, z, 1) 

4681 >>> check 

4682 True 

4683 >>> rels = rels_k + rels_h 

4684 >>> G = FpGroup(F, rels) # presentation of H 

4685 >>> G.order() == H.order() 

4686 True 

4687 

4688 See also 

4689 ======== 

4690 

4691 strong_presentation, presentation, stabilizer 

4692 

4693 ''' 

4694 

4695 orbit = H.orbit(alpha) 

4696 beta = alpha^(z**-1) 

4697 

4698 K_beta = K.stabilizer(beta) 

4699 

4700 # orbit representatives of K_beta 

4701 gammas = [alpha, beta] 

4702 orbits = list({tuple(K_beta.orbit(o)) for o in orbit}) 

4703 orbit_reps = [orb[0] for orb in orbits] 

4704 for rep in orbit_reps: 

4705 if rep not in gammas: 

4706 gammas.append(rep) 

4707 

4708 # orbit transversal of K 

4709 betas = [alpha, beta] 

4710 transversal = {alpha: phi.invert(H.identity), beta: phi.invert(z**-1)} 

4711 

4712 for s, g in K.orbit_transversal(beta, pairs=True): 

4713 if s not in transversal: 

4714 transversal[s] = transversal[beta]*phi.invert(g) 

4715 

4716 

4717 union = K.orbit(alpha).union(K.orbit(beta)) 

4718 while (len(union) < len(orbit)): 

4719 for gamma in gammas: 

4720 if gamma in union: 

4721 r = gamma^z 

4722 if r not in union: 

4723 betas.append(r) 

4724 transversal[r] = transversal[gamma]*phi.invert(z) 

4725 for s, g in K.orbit_transversal(r, pairs=True): 

4726 if s not in transversal: 

4727 transversal[s] = transversal[r]*phi.invert(g) 

4728 union = union.union(K.orbit(r)) 

4729 break 

4730 

4731 # compute relators 

4732 rels = [] 

4733 

4734 for b in betas: 

4735 k_gens = K.stabilizer(b).generators 

4736 for y in k_gens: 

4737 new_rel = transversal[b] 

4738 gens = K.generator_product(y, original=True) 

4739 for g in gens[::-1]: 

4740 new_rel = new_rel*phi.invert(g) 

4741 new_rel = new_rel*transversal[b]**-1 

4742 

4743 perm = phi(new_rel) 

4744 try: 

4745 gens = K.generator_product(perm, original=True) 

4746 except ValueError: 

4747 return False, perm 

4748 for g in gens: 

4749 new_rel = new_rel*phi.invert(g)**-1 

4750 if new_rel not in rels: 

4751 rels.append(new_rel) 

4752 

4753 for gamma in gammas: 

4754 new_rel = transversal[gamma]*phi.invert(z)*transversal[gamma^z]**-1 

4755 perm = phi(new_rel) 

4756 try: 

4757 gens = K.generator_product(perm, original=True) 

4758 except ValueError: 

4759 return False, perm 

4760 for g in gens: 

4761 new_rel = new_rel*phi.invert(g)**-1 

4762 if new_rel not in rels: 

4763 rels.append(new_rel) 

4764 

4765 return True, rels 

4766 

4767 def strong_presentation(self): 

4768 ''' 

4769 Return a strong finite presentation of group. The generators 

4770 of the returned group are in the same order as the strong 

4771 generators of group. 

4772 

4773 The algorithm is based on Sims' Verify algorithm described 

4774 in [1], Chapter 6. 

4775 

4776 Examples 

4777 ======== 

4778 

4779 >>> from sympy.combinatorics.named_groups import DihedralGroup 

4780 >>> P = DihedralGroup(4) 

4781 >>> G = P.strong_presentation() 

4782 >>> P.order() == G.order() 

4783 True 

4784 

4785 See Also 

4786 ======== 

4787 

4788 presentation, _verify 

4789 

4790 ''' 

4791 from sympy.combinatorics.fp_groups import (FpGroup, 

4792 simplify_presentation) 

4793 from sympy.combinatorics.free_groups import free_group 

4794 from sympy.combinatorics.homomorphisms import (block_homomorphism, 

4795 homomorphism, GroupHomomorphism) 

4796 

4797 strong_gens = self.strong_gens[:] 

4798 stabs = self.basic_stabilizers[:] 

4799 base = self.base[:] 

4800 

4801 # injection from a free group on len(strong_gens) 

4802 # generators into G 

4803 gen_syms = [('x_%d'%i) for i in range(len(strong_gens))] 

4804 F = free_group(', '.join(gen_syms))[0] 

4805 phi = homomorphism(F, self, F.generators, strong_gens) 

4806 

4807 H = PermutationGroup(self.identity) 

4808 while stabs: 

4809 alpha = base.pop() 

4810 K = H 

4811 H = stabs.pop() 

4812 new_gens = [g for g in H.generators if g not in K] 

4813 

4814 if K.order() == 1: 

4815 z = new_gens.pop() 

4816 rels = [F.generators[-1]**z.order()] 

4817 intermediate_gens = [z] 

4818 K = PermutationGroup(intermediate_gens) 

4819 

4820 # add generators one at a time building up from K to H 

4821 while new_gens: 

4822 z = new_gens.pop() 

4823 intermediate_gens = [z] + intermediate_gens 

4824 K_s = PermutationGroup(intermediate_gens) 

4825 orbit = K_s.orbit(alpha) 

4826 orbit_k = K.orbit(alpha) 

4827 

4828 # split into cases based on the orbit of K_s 

4829 if orbit_k == orbit: 

4830 if z in K: 

4831 rel = phi.invert(z) 

4832 perm = z 

4833 else: 

4834 t = K.orbit_rep(alpha, alpha^z) 

4835 rel = phi.invert(z)*phi.invert(t)**-1 

4836 perm = z*t**-1 

4837 for g in K.generator_product(perm, original=True): 

4838 rel = rel*phi.invert(g)**-1 

4839 new_rels = [rel] 

4840 elif len(orbit_k) == 1: 

4841 # `success` is always true because `strong_gens` 

4842 # and `base` are already a verified BSGS. Later 

4843 # this could be changed to start with a randomly 

4844 # generated (potential) BSGS, and then new elements 

4845 # would have to be appended to it when `success` 

4846 # is false. 

4847 success, new_rels = K_s._verify(K, phi, z, alpha) 

4848 else: 

4849 # K.orbit(alpha) should be a block 

4850 # under the action of K_s on K_s.orbit(alpha) 

4851 check, block = K_s._block_verify(K, alpha) 

4852 if check: 

4853 # apply _verify to the action of K_s 

4854 # on the block system; for convenience, 

4855 # add the blocks as additional points 

4856 # that K_s should act on 

4857 t = block_homomorphism(K_s, block) 

4858 m = t.codomain.degree # number of blocks 

4859 d = K_s.degree 

4860 

4861 # conjugating with p will shift 

4862 # permutations in t.image() to 

4863 # higher numbers, e.g. 

4864 # p*(0 1)*p = (m m+1) 

4865 p = Permutation() 

4866 for i in range(m): 

4867 p *= Permutation(i, i+d) 

4868 

4869 t_img = t.images 

4870 # combine generators of K_s with their 

4871 # action on the block system 

4872 images = {g: g*p*t_img[g]*p for g in t_img} 

4873 for g in self.strong_gens[:-len(K_s.generators)]: 

4874 images[g] = g 

4875 K_s_act = PermutationGroup(list(images.values())) 

4876 f = GroupHomomorphism(self, K_s_act, images) 

4877 

4878 K_act = PermutationGroup([f(g) for g in K.generators]) 

4879 success, new_rels = K_s_act._verify(K_act, f.compose(phi), f(z), d) 

4880 

4881 for n in new_rels: 

4882 if n not in rels: 

4883 rels.append(n) 

4884 K = K_s 

4885 

4886 group = FpGroup(F, rels) 

4887 return simplify_presentation(group) 

4888 

4889 def presentation(self, eliminate_gens=True): 

4890 ''' 

4891 Return an `FpGroup` presentation of the group. 

4892 

4893 The algorithm is described in [1], Chapter 6.1. 

4894 

4895 ''' 

4896 from sympy.combinatorics.fp_groups import (FpGroup, 

4897 simplify_presentation) 

4898 from sympy.combinatorics.coset_table import CosetTable 

4899 from sympy.combinatorics.free_groups import free_group 

4900 from sympy.combinatorics.homomorphisms import homomorphism 

4901 

4902 if self._fp_presentation: 

4903 return self._fp_presentation 

4904 

4905 def _factor_group_by_rels(G, rels): 

4906 if isinstance(G, FpGroup): 

4907 rels.extend(G.relators) 

4908 return FpGroup(G.free_group, list(set(rels))) 

4909 return FpGroup(G, rels) 

4910 

4911 gens = self.generators 

4912 len_g = len(gens) 

4913 

4914 if len_g == 1: 

4915 order = gens[0].order() 

4916 # handle the trivial group 

4917 if order == 1: 

4918 return free_group([])[0] 

4919 F, x = free_group('x') 

4920 return FpGroup(F, [x**order]) 

4921 

4922 if self.order() > 20: 

4923 half_gens = self.generators[0:(len_g+1)//2] 

4924 else: 

4925 half_gens = [] 

4926 H = PermutationGroup(half_gens) 

4927 H_p = H.presentation() 

4928 

4929 len_h = len(H_p.generators) 

4930 

4931 C = self.coset_table(H) 

4932 n = len(C) # subgroup index 

4933 

4934 gen_syms = [('x_%d'%i) for i in range(len(gens))] 

4935 F = free_group(', '.join(gen_syms))[0] 

4936 

4937 # mapping generators of H_p to those of F 

4938 images = [F.generators[i] for i in range(len_h)] 

4939 R = homomorphism(H_p, F, H_p.generators, images, check=False) 

4940 

4941 # rewrite relators 

4942 rels = R(H_p.relators) 

4943 G_p = FpGroup(F, rels) 

4944 

4945 # injective homomorphism from G_p into self 

4946 T = homomorphism(G_p, self, G_p.generators, gens) 

4947 

4948 C_p = CosetTable(G_p, []) 

4949 

4950 C_p.table = [[None]*(2*len_g) for i in range(n)] 

4951 

4952 # initiate the coset transversal 

4953 transversal = [None]*n 

4954 transversal[0] = G_p.identity 

4955 

4956 # fill in the coset table as much as possible 

4957 for i in range(2*len_h): 

4958 C_p.table[0][i] = 0 

4959 

4960 gamma = 1 

4961 for alpha, x in product(range(n), range(2*len_g)): 

4962 beta = C[alpha][x] 

4963 if beta == gamma: 

4964 gen = G_p.generators[x//2]**((-1)**(x % 2)) 

4965 transversal[beta] = transversal[alpha]*gen 

4966 C_p.table[alpha][x] = beta 

4967 C_p.table[beta][x + (-1)**(x % 2)] = alpha 

4968 gamma += 1 

4969 if gamma == n: 

4970 break 

4971 

4972 C_p.p = list(range(n)) 

4973 beta = x = 0 

4974 

4975 while not C_p.is_complete(): 

4976 # find the first undefined entry 

4977 while C_p.table[beta][x] == C[beta][x]: 

4978 x = (x + 1) % (2*len_g) 

4979 if x == 0: 

4980 beta = (beta + 1) % n 

4981 

4982 # define a new relator 

4983 gen = G_p.generators[x//2]**((-1)**(x % 2)) 

4984 new_rel = transversal[beta]*gen*transversal[C[beta][x]]**-1 

4985 perm = T(new_rel) 

4986 nxt = G_p.identity 

4987 for s in H.generator_product(perm, original=True): 

4988 nxt = nxt*T.invert(s)**-1 

4989 new_rel = new_rel*nxt 

4990 

4991 # continue coset enumeration 

4992 G_p = _factor_group_by_rels(G_p, [new_rel]) 

4993 C_p.scan_and_fill(0, new_rel) 

4994 C_p = G_p.coset_enumeration([], strategy="coset_table", 

4995 draft=C_p, max_cosets=n, incomplete=True) 

4996 

4997 self._fp_presentation = simplify_presentation(G_p) 

4998 return self._fp_presentation 

4999 

5000 def polycyclic_group(self): 

5001 """ 

5002 Return the PolycyclicGroup instance with below parameters: 

5003 

5004 Explanation 

5005 =========== 

5006 

5007 * pc_sequence : Polycyclic sequence is formed by collecting all 

5008 the missing generators between the adjacent groups in the 

5009 derived series of given permutation group. 

5010 

5011 * pc_series : Polycyclic series is formed by adding all the missing 

5012 generators of ``der[i+1]`` in ``der[i]``, where ``der`` represents 

5013 the derived series. 

5014 

5015 * relative_order : A list, computed by the ratio of adjacent groups in 

5016 pc_series. 

5017 

5018 """ 

5019 from sympy.combinatorics.pc_groups import PolycyclicGroup 

5020 if not self.is_polycyclic: 

5021 raise ValueError("The group must be solvable") 

5022 

5023 der = self.derived_series() 

5024 pc_series = [] 

5025 pc_sequence = [] 

5026 relative_order = [] 

5027 pc_series.append(der[-1]) 

5028 der.reverse() 

5029 

5030 for i in range(len(der)-1): 

5031 H = der[i] 

5032 for g in der[i+1].generators: 

5033 if g not in H: 

5034 H = PermutationGroup([g] + H.generators) 

5035 pc_series.insert(0, H) 

5036 pc_sequence.insert(0, g) 

5037 

5038 G1 = pc_series[0].order() 

5039 G2 = pc_series[1].order() 

5040 relative_order.insert(0, G1 // G2) 

5041 

5042 return PolycyclicGroup(pc_sequence, pc_series, relative_order, collector=None) 

5043 

5044 

5045def _orbit(degree, generators, alpha, action='tuples'): 

5046 r"""Compute the orbit of alpha `\{g(\alpha) | g \in G\}` as a set. 

5047 

5048 Explanation 

5049 =========== 

5050 

5051 The time complexity of the algorithm used here is `O(|Orb|*r)` where 

5052 `|Orb|` is the size of the orbit and ``r`` is the number of generators of 

5053 the group. For a more detailed analysis, see [1], p.78, [2], pp. 19-21. 

5054 Here alpha can be a single point, or a list of points. 

5055 

5056 If alpha is a single point, the ordinary orbit is computed. 

5057 if alpha is a list of points, there are three available options: 

5058 

5059 'union' - computes the union of the orbits of the points in the list 

5060 'tuples' - computes the orbit of the list interpreted as an ordered 

5061 tuple under the group action ( i.e., g((1, 2, 3)) = (g(1), g(2), g(3)) ) 

5062 'sets' - computes the orbit of the list interpreted as a sets 

5063 

5064 Examples 

5065 ======== 

5066 

5067 >>> from sympy.combinatorics import Permutation, PermutationGroup 

5068 >>> from sympy.combinatorics.perm_groups import _orbit 

5069 >>> a = Permutation([1, 2, 0, 4, 5, 6, 3]) 

5070 >>> G = PermutationGroup([a]) 

5071 >>> _orbit(G.degree, G.generators, 0) 

5072 {0, 1, 2} 

5073 >>> _orbit(G.degree, G.generators, [0, 4], 'union') 

5074 {0, 1, 2, 3, 4, 5, 6} 

5075 

5076 See Also 

5077 ======== 

5078 

5079 orbit, orbit_transversal 

5080 

5081 """ 

5082 if not hasattr(alpha, '__getitem__'): 

5083 alpha = [alpha] 

5084 

5085 gens = [x._array_form for x in generators] 

5086 if len(alpha) == 1 or action == 'union': 

5087 orb = alpha 

5088 used = [False]*degree 

5089 for el in alpha: 

5090 used[el] = True 

5091 for b in orb: 

5092 for gen in gens: 

5093 temp = gen[b] 

5094 if used[temp] == False: 

5095 orb.append(temp) 

5096 used[temp] = True 

5097 return set(orb) 

5098 elif action == 'tuples': 

5099 alpha = tuple(alpha) 

5100 orb = [alpha] 

5101 used = {alpha} 

5102 for b in orb: 

5103 for gen in gens: 

5104 temp = tuple([gen[x] for x in b]) 

5105 if temp not in used: 

5106 orb.append(temp) 

5107 used.add(temp) 

5108 return set(orb) 

5109 elif action == 'sets': 

5110 alpha = frozenset(alpha) 

5111 orb = [alpha] 

5112 used = {alpha} 

5113 for b in orb: 

5114 for gen in gens: 

5115 temp = frozenset([gen[x] for x in b]) 

5116 if temp not in used: 

5117 orb.append(temp) 

5118 used.add(temp) 

5119 return {tuple(x) for x in orb} 

5120 

5121 

5122def _orbits(degree, generators): 

5123 """Compute the orbits of G. 

5124 

5125 If ``rep=False`` it returns a list of sets else it returns a list of 

5126 representatives of the orbits 

5127 

5128 Examples 

5129 ======== 

5130 

5131 >>> from sympy.combinatorics import Permutation 

5132 >>> from sympy.combinatorics.perm_groups import _orbits 

5133 >>> a = Permutation([0, 2, 1]) 

5134 >>> b = Permutation([1, 0, 2]) 

5135 >>> _orbits(a.size, [a, b]) 

5136 [{0, 1, 2}] 

5137 """ 

5138 

5139 orbs = [] 

5140 sorted_I = list(range(degree)) 

5141 I = set(sorted_I) 

5142 while I: 

5143 i = sorted_I[0] 

5144 orb = _orbit(degree, generators, i) 

5145 orbs.append(orb) 

5146 # remove all indices that are in this orbit 

5147 I -= orb 

5148 sorted_I = [i for i in sorted_I if i not in orb] 

5149 return orbs 

5150 

5151 

5152def _orbit_transversal(degree, generators, alpha, pairs, af=False, slp=False): 

5153 r"""Computes a transversal for the orbit of ``alpha`` as a set. 

5154 

5155 Explanation 

5156 =========== 

5157 

5158 generators generators of the group ``G`` 

5159 

5160 For a permutation group ``G``, a transversal for the orbit 

5161 `Orb = \{g(\alpha) | g \in G\}` is a set 

5162 `\{g_\beta | g_\beta(\alpha) = \beta\}` for `\beta \in Orb`. 

5163 Note that there may be more than one possible transversal. 

5164 If ``pairs`` is set to ``True``, it returns the list of pairs 

5165 `(\beta, g_\beta)`. For a proof of correctness, see [1], p.79 

5166 

5167 if ``af`` is ``True``, the transversal elements are given in 

5168 array form. 

5169 

5170 If `slp` is `True`, a dictionary `{beta: slp_beta}` is returned 

5171 for `\beta \in Orb` where `slp_beta` is a list of indices of the 

5172 generators in `generators` s.t. if `slp_beta = [i_1 \dots i_n]` 

5173 `g_\beta = generators[i_n] \times \dots \times generators[i_1]`. 

5174 

5175 Examples 

5176 ======== 

5177 

5178 >>> from sympy.combinatorics.named_groups import DihedralGroup 

5179 >>> from sympy.combinatorics.perm_groups import _orbit_transversal 

5180 >>> G = DihedralGroup(6) 

5181 >>> _orbit_transversal(G.degree, G.generators, 0, False) 

5182 [(5), (0 1 2 3 4 5), (0 5)(1 4)(2 3), (0 2 4)(1 3 5), (5)(0 4)(1 3), (0 3)(1 4)(2 5)] 

5183 """ 

5184 

5185 tr = [(alpha, list(range(degree)))] 

5186 slp_dict = {alpha: []} 

5187 used = [False]*degree 

5188 used[alpha] = True 

5189 gens = [x._array_form for x in generators] 

5190 for x, px in tr: 

5191 px_slp = slp_dict[x] 

5192 for gen in gens: 

5193 temp = gen[x] 

5194 if used[temp] == False: 

5195 slp_dict[temp] = [gens.index(gen)] + px_slp 

5196 tr.append((temp, _af_rmul(gen, px))) 

5197 used[temp] = True 

5198 if pairs: 

5199 if not af: 

5200 tr = [(x, _af_new(y)) for x, y in tr] 

5201 if not slp: 

5202 return tr 

5203 return tr, slp_dict 

5204 

5205 if af: 

5206 tr = [y for _, y in tr] 

5207 if not slp: 

5208 return tr 

5209 return tr, slp_dict 

5210 

5211 tr = [_af_new(y) for _, y in tr] 

5212 if not slp: 

5213 return tr 

5214 return tr, slp_dict 

5215 

5216 

5217def _stabilizer(degree, generators, alpha): 

5218 r"""Return the stabilizer subgroup of ``alpha``. 

5219 

5220 Explanation 

5221 =========== 

5222 

5223 The stabilizer of `\alpha` is the group `G_\alpha = 

5224 \{g \in G | g(\alpha) = \alpha\}`. 

5225 For a proof of correctness, see [1], p.79. 

5226 

5227 degree : degree of G 

5228 generators : generators of G 

5229 

5230 Examples 

5231 ======== 

5232 

5233 >>> from sympy.combinatorics.perm_groups import _stabilizer 

5234 >>> from sympy.combinatorics.named_groups import DihedralGroup 

5235 >>> G = DihedralGroup(6) 

5236 >>> _stabilizer(G.degree, G.generators, 5) 

5237 [(5)(0 4)(1 3), (5)] 

5238 

5239 See Also 

5240 ======== 

5241 

5242 orbit 

5243 

5244 """ 

5245 orb = [alpha] 

5246 table = {alpha: list(range(degree))} 

5247 table_inv = {alpha: list(range(degree))} 

5248 used = [False]*degree 

5249 used[alpha] = True 

5250 gens = [x._array_form for x in generators] 

5251 stab_gens = [] 

5252 for b in orb: 

5253 for gen in gens: 

5254 temp = gen[b] 

5255 if used[temp] is False: 

5256 gen_temp = _af_rmul(gen, table[b]) 

5257 orb.append(temp) 

5258 table[temp] = gen_temp 

5259 table_inv[temp] = _af_invert(gen_temp) 

5260 used[temp] = True 

5261 else: 

5262 schreier_gen = _af_rmuln(table_inv[temp], gen, table[b]) 

5263 if schreier_gen not in stab_gens: 

5264 stab_gens.append(schreier_gen) 

5265 return [_af_new(x) for x in stab_gens] 

5266 

5267 

5268PermGroup = PermutationGroup 

5269 

5270 

5271class SymmetricPermutationGroup(Basic): 

5272 """ 

5273 The class defining the lazy form of SymmetricGroup. 

5274 

5275 deg : int 

5276 

5277 """ 

5278 def __new__(cls, deg): 

5279 deg = _sympify(deg) 

5280 obj = Basic.__new__(cls, deg) 

5281 return obj 

5282 

5283 def __init__(self, *args, **kwargs): 

5284 self._deg = self.args[0] 

5285 self._order = None 

5286 

5287 def __contains__(self, i): 

5288 """Return ``True`` if *i* is contained in SymmetricPermutationGroup. 

5289 

5290 Examples 

5291 ======== 

5292 

5293 >>> from sympy.combinatorics import Permutation, SymmetricPermutationGroup 

5294 >>> G = SymmetricPermutationGroup(4) 

5295 >>> Permutation(1, 2, 3) in G 

5296 True 

5297 

5298 """ 

5299 if not isinstance(i, Permutation): 

5300 raise TypeError("A SymmetricPermutationGroup contains only Permutations as " 

5301 "elements, not elements of type %s" % type(i)) 

5302 return i.size == self.degree 

5303 

5304 def order(self): 

5305 """ 

5306 Return the order of the SymmetricPermutationGroup. 

5307 

5308 Examples 

5309 ======== 

5310 

5311 >>> from sympy.combinatorics import SymmetricPermutationGroup 

5312 >>> G = SymmetricPermutationGroup(4) 

5313 >>> G.order() 

5314 24 

5315 """ 

5316 if self._order is not None: 

5317 return self._order 

5318 n = self._deg 

5319 self._order = factorial(n) 

5320 return self._order 

5321 

5322 @property 

5323 def degree(self): 

5324 """ 

5325 Return the degree of the SymmetricPermutationGroup. 

5326 

5327 Examples 

5328 ======== 

5329 

5330 >>> from sympy.combinatorics import SymmetricPermutationGroup 

5331 >>> G = SymmetricPermutationGroup(4) 

5332 >>> G.degree 

5333 4 

5334 

5335 """ 

5336 return self._deg 

5337 

5338 @property 

5339 def identity(self): 

5340 ''' 

5341 Return the identity element of the SymmetricPermutationGroup. 

5342 

5343 Examples 

5344 ======== 

5345 

5346 >>> from sympy.combinatorics import SymmetricPermutationGroup 

5347 >>> G = SymmetricPermutationGroup(4) 

5348 >>> G.identity() 

5349 (3) 

5350 

5351 ''' 

5352 return _af_new(list(range(self._deg))) 

5353 

5354 

5355class Coset(Basic): 

5356 """A left coset of a permutation group with respect to an element. 

5357 

5358 Parameters 

5359 ========== 

5360 

5361 g : Permutation 

5362 

5363 H : PermutationGroup 

5364 

5365 dir : "+" or "-", If not specified by default it will be "+" 

5366 here ``dir`` specified the type of coset "+" represent the 

5367 right coset and "-" represent the left coset. 

5368 

5369 G : PermutationGroup, optional 

5370 The group which contains *H* as its subgroup and *g* as its 

5371 element. 

5372 

5373 If not specified, it would automatically become a symmetric 

5374 group ``SymmetricPermutationGroup(g.size)`` and 

5375 ``SymmetricPermutationGroup(H.degree)`` if ``g.size`` and ``H.degree`` 

5376 are matching.``SymmetricPermutationGroup`` is a lazy form of SymmetricGroup 

5377 used for representation purpose. 

5378 

5379 """ 

5380 

5381 def __new__(cls, g, H, G=None, dir="+"): 

5382 g = _sympify(g) 

5383 if not isinstance(g, Permutation): 

5384 raise NotImplementedError 

5385 

5386 H = _sympify(H) 

5387 if not isinstance(H, PermutationGroup): 

5388 raise NotImplementedError 

5389 

5390 if G is not None: 

5391 G = _sympify(G) 

5392 if not isinstance(G, (PermutationGroup, SymmetricPermutationGroup)): 

5393 raise NotImplementedError 

5394 if not H.is_subgroup(G): 

5395 raise ValueError("{} must be a subgroup of {}.".format(H, G)) 

5396 if g not in G: 

5397 raise ValueError("{} must be an element of {}.".format(g, G)) 

5398 else: 

5399 g_size = g.size 

5400 h_degree = H.degree 

5401 if g_size != h_degree: 

5402 raise ValueError( 

5403 "The size of the permutation {} and the degree of " 

5404 "the permutation group {} should be matching " 

5405 .format(g, H)) 

5406 G = SymmetricPermutationGroup(g.size) 

5407 

5408 if isinstance(dir, str): 

5409 dir = Symbol(dir) 

5410 elif not isinstance(dir, Symbol): 

5411 raise TypeError("dir must be of type basestring or " 

5412 "Symbol, not %s" % type(dir)) 

5413 if str(dir) not in ('+', '-'): 

5414 raise ValueError("dir must be one of '+' or '-' not %s" % dir) 

5415 obj = Basic.__new__(cls, g, H, G, dir) 

5416 return obj 

5417 

5418 def __init__(self, *args, **kwargs): 

5419 self._dir = self.args[3] 

5420 

5421 @property 

5422 def is_left_coset(self): 

5423 """ 

5424 Check if the coset is left coset that is ``gH``. 

5425 

5426 Examples 

5427 ======== 

5428 

5429 >>> from sympy.combinatorics import Permutation, PermutationGroup, Coset 

5430 >>> a = Permutation(1, 2) 

5431 >>> b = Permutation(0, 1) 

5432 >>> G = PermutationGroup([a, b]) 

5433 >>> cst = Coset(a, G, dir="-") 

5434 >>> cst.is_left_coset 

5435 True 

5436 

5437 """ 

5438 return str(self._dir) == '-' 

5439 

5440 @property 

5441 def is_right_coset(self): 

5442 """ 

5443 Check if the coset is right coset that is ``Hg``. 

5444 

5445 Examples 

5446 ======== 

5447 

5448 >>> from sympy.combinatorics import Permutation, PermutationGroup, Coset 

5449 >>> a = Permutation(1, 2) 

5450 >>> b = Permutation(0, 1) 

5451 >>> G = PermutationGroup([a, b]) 

5452 >>> cst = Coset(a, G, dir="+") 

5453 >>> cst.is_right_coset 

5454 True 

5455 

5456 """ 

5457 return str(self._dir) == '+' 

5458 

5459 def as_list(self): 

5460 """ 

5461 Return all the elements of coset in the form of list. 

5462 """ 

5463 g = self.args[0] 

5464 H = self.args[1] 

5465 cst = [] 

5466 if str(self._dir) == '+': 

5467 for h in H.elements: 

5468 cst.append(h*g) 

5469 else: 

5470 for h in H.elements: 

5471 cst.append(g*h) 

5472 return cst