Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690

691

692

693

694

695

696

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722

723

724

725

726

727

728

729

730

731

732

733

734

735

736

737

738

739

740

741

742

743

744

745

746

747

748

749

750

751

752

753

754

755

756

757

758

759

760

761

762

763

764

765

766

767

768

769

770

771

772

773

774

775

776

777

778

779

780

781

782

783

784

785

786

787

788

789

790

791

792

793

794

795

796

797

798

799

800

801

802

803

804

805

806

807

808

809

810

811

812

813

814

815

816

817

818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

839

840

841

842

843

844

845

846

847

848

849

850

851

852

853

854

855

856

857

858

859

860

861

862

863

864

865

866

867

868

869

870

871

872

873

874

875

876

877

878

879

880

881

882

883

884

885

886

887

888

889

890

891

892

893

894

895

896

897

898

899

900

901

902

903

904

905

906

907

908

909

910

911

912

913

914

915

916

917

918

919

920

921

922

923

924

925

926

927

928

929

930

931

932

933

934

935

936

937

938

939

940

941

942

943

944

945

946

947

948

949

950

951

952

953

954

955

956

957

958

959

960

961

962

963

964

965

966

967

968

969

970

971

972

973

974

975

976

977

978

979

980

981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

999

1000

1001

1002

1003

1004

1005

1006

1007

1008

1009

1010

1011

1012

1013

1014

1015

1016

1017

1018

1019

1020

1021

1022

1023

1024

1025

1026

1027

1028

1029

1030

1031

1032

1033

1034

1035

1036

1037

1038

1039

1040

1041

1042

1043

1044

1045

1046

1047

1048

1049

1050

1051

1052

1053

1054

1055

1056

1057

1058

1059

1060

1061

1062

1063

1064

1065

1066

1067

1068

1069

1070

1071

1072

1073

1074

1075

1076

1077

1078

1079

1080

1081

1082

1083

1084

1085

1086

1087

1088

1089

1090

1091

1092

1093

1094

1095

1096

1097

1098

1099

1100

1101

1102

1103

1104

1105

1106

1107

1108

1109

1110

1111

1112

1113

1114

1115

1116

1117

1118

1119

1120

1121

1122

1123

1124

1125

1126

1127

1128

1129

1130

1131

1132

1133

1134

1135

1136

1137

1138

1139

1140

1141

1142

1143

1144

1145

1146

1147

1148

1149

1150

1151

1152

1153

1154

1155

1156

1157

1158

1159

1160

1161

1162

1163

1164

1165

1166

1167

1168

1169

1170

1171

1172

1173

1174

1175

1176

1177

1178

1179

1180

1181

1182

1183

1184

1185

1186

1187

1188

1189

1190

1191

1192

1193

1194

1195

1196

1197

1198

1199

1200

1201

1202

1203

1204

1205

1206

1207

1208

1209

1210

1211

1212

1213

1214

1215

1216

1217

1218

1219

1220

1221

1222

1223

1224

1225

1226

1227

1228

1229

1230

1231

1232

1233

1234

1235

1236

1237

1238

1239

1240

1241

1242

1243

1244

1245

1246

1247

1248

1249

1250

1251

1252

1253

1254

1255

1256

1257

1258

1259

1260

1261

1262

1263

1264

1265

1266

1267

1268

1269

1270

1271

1272

1273

1274

1275

1276

1277

1278

1279

1280

1281

1282

1283

1284

1285

1286

1287

1288

1289

1290

1291

1292

1293

1294

1295

1296

1297

1298

1299

1300

1301

1302

1303

1304

1305

1306

1307

1308

1309

1310

1311

1312

1313

1314

1315

1316

1317

1318

1319

1320

1321

1322

1323

1324

1325

1326

1327

1328

1329

1330

1331

1332

1333

1334

1335

1336

1337

1338

1339

1340

1341

1342

1343

1344

1345

1346

1347

1348

1349

1350

1351

1352

1353

1354

1355

1356

1357

1358

1359

1360

1361

1362

1363

1364

1365

1366

1367

1368

1369

1370

1371

1372

1373

1374

1375

1376

1377

1378

1379

1380

1381

1382

1383

1384

1385

1386

1387

1388

1389

1390

1391

1392

1393

1394

1395

1396

1397

1398

1399

1400

1401

1402

1403

1404

1405

1406

1407

1408

1409

1410

1411

1412

1413

1414

1415

1416

1417

1418

1419

1420

1421

1422

1423

1424

1425

1426

1427

1428

1429

1430

1431

1432

1433

1434

1435

1436

1437

1438

1439

1440

1441

1442

1443

1444

1445

1446

1447

1448

1449

1450

1451

1452

1453

1454

1455

1456

1457

1458

1459

1460

1461

1462

1463

1464

1465

1466

1467

1468

1469

1470

1471

1472

1473

1474

1475

1476

1477

1478

1479

1480

1481

1482

1483

1484

1485

1486

1487

1488

1489

1490

1491

1492

1493

1494

1495

1496

1497

1498

1499

1500

1501

1502

1503

1504

1505

1506

1507

1508

1509

1510

1511

1512

1513

1514

1515

1516

1517

1518

1519

1520

1521

1522

1523

1524

1525

1526

1527

1528

1529

1530

1531

1532

1533

1534

1535

1536

1537

1538

1539

1540

1541

1542

1543

1544

1545

1546

1547

1548

1549

1550

1551

1552

1553

1554

1555

1556

1557

1558

1559

1560

1561

1562

1563

1564

1565

1566

1567

1568

1569

1570

1571

1572

1573

1574

1575

1576

1577

1578

1579

1580

1581

1582

1583

1584

1585

1586

1587

1588

1589

1590

1591

1592

1593

1594

1595

1596

1597

1598

1599

1600

1601

1602

1603

1604

1605

1606

1607

1608

1609

1610

1611

1612

1613

1614

1615

1616

1617

1618

1619

1620

1621

1622

1623

1624

1625

1626

1627

1628

1629

1630

1631

1632

1633

1634

1635

1636

1637

1638

1639

1640

1641

1642

1643

1644

1645

1646

1647

1648

1649

1650

1651

1652

1653

1654

1655

1656

1657

1658

1659

1660

1661

1662

1663

1664

1665

1666

1667

1668

1669

1670

1671

1672

1673

1674

1675

1676

1677

1678

1679

1680

1681

1682

1683

1684

1685

1686

1687

1688

1689

1690

1691

1692

1693

1694

1695

1696

1697

1698

1699

1700

1701

1702

1703

1704

1705

1706

1707

1708

1709

1710

1711

1712

1713

1714

1715

r""" 

Data structures 

  

This module implements basic data structures essential to the rest of the 

partn_ref module. 

  

REFERENCES: 

  

[1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium, 

Vol. 30 (1981), pp. 45-87. 

[2] Fredman, M. and Saks, M. The cell probe complexity of dynamic data 

structures. Proceedings of the Twenty-First Annual ACM Symposium on 

Theory of Computing, pp. 345–354. May 1989. 

[3] Seress, Akos. Permutation Group Algorithms. Cambridge University Press, 

2003. 

  

""" 

  

#***************************************************************************** 

# Copyright (C) 2006 - 2011 Robert L. Miller <rlmillster@gmail.com> 

# 

# This program is free software: you can redistribute it and/or modify 

# it under the terms of the GNU General Public License as published by 

# the Free Software Foundation, either version 2 of the License, or 

# (at your option) any later version. 

# http://www.gnu.org/licenses/ 

#***************************************************************************** 

  

from libc.math cimport log, ceil 

from libc.string cimport memcpy, memset 

from libc.stdlib cimport rand 

from cysignals.memory cimport sig_malloc, sig_calloc, sig_realloc, sig_free 

  

include "sage/data_structures/bitset.pxi" 

from sage.rings.integer cimport Integer 

from sage.libs.flint.ulong_extras cimport n_is_prime 

  

  

# OrbitPartition (OP) 

  

cdef inline OrbitPartition *OP_new(int n): 

""" 

Allocate and return a pointer to a new OrbitPartition of degree n. Returns a 

null pointer in the case of an allocation failure. 

""" 

cdef int i 

cdef OrbitPartition *OP = <OrbitPartition *> \ 

sig_malloc(sizeof(OrbitPartition)) 

cdef int *int_array = <int *> sig_malloc( 4*n * sizeof(int) ) 

if OP is NULL or int_array is NULL: 

sig_free(OP) 

sig_free(int_array) 

return NULL 

OP.degree = n 

OP.num_cells = n 

OP.parent = int_array 

OP.rank = int_array + n 

OP.mcr = int_array + 2*n 

OP.size = int_array + 3*n 

OP_clear(OP) 

return OP 

  

cdef inline void OP_dealloc(OrbitPartition *OP): 

if OP is not NULL: 

sig_free(OP.parent) 

sig_free(OP) 

  

cdef OP_string(OrbitPartition *OP): 

""" 

Return a string representation of the OrbitPartition. 

""" 

cdef i,j 

s = "" 

for i from 0 <= i < OP.degree: 

s += " " 

j = OP_find(OP, i) 

s += "%d -> %d"%(i, j) 

return s 

  

def OP_represent(int n, merges, perm): 

""" 

Demonstration and testing. 

  

TESTS:: 

  

sage: from sage.groups.perm_gps.partn_ref.data_structures import OP_represent 

sage: OP_represent(9, [(0,1),(2,3),(3,4)], [1,2,0,4,3,6,7,5,8]) 

Allocating OrbitPartition... 

Allocation passed. 

Checking that each element reports itself as its root. 

Each element reports itself as its root. 

Merging: 

Merged 0 and 1. 

Merged 2 and 3. 

Merged 3 and 4. 

Done merging. 

Finding: 

0 -> 0, root: size=2, mcr=0, rank=1 

1 -> 0 

2 -> 2, root: size=3, mcr=2, rank=1 

3 -> 2 

4 -> 2 

5 -> 5, root: size=1, mcr=5, rank=0 

6 -> 6, root: size=1, mcr=6, rank=0 

7 -> 7, root: size=1, mcr=7, rank=0 

8 -> 8, root: size=1, mcr=8, rank=0 

Allocating array to test merge_perm. 

Allocation passed. 

Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8] 

Done merging. 

Finding: 

0 -> 0, root: size=5, mcr=0, rank=2 

1 -> 0 

2 -> 0 

3 -> 0 

4 -> 0 

5 -> 5, root: size=3, mcr=5, rank=1 

6 -> 5 

7 -> 5 

8 -> 8, root: size=1, mcr=8, rank=0 

Deallocating OrbitPartition. 

Done. 

  

""" 

cdef int i 

print("Allocating OrbitPartition...") 

cdef OrbitPartition *OP = OP_new(n) 

if OP is NULL: 

print("Allocation failed!") 

return 

print("Allocation passed.") 

print("Checking that each element reports itself as its root.") 

good = True 

for i from 0 <= i < n: 

if not OP_find(OP, i) == i: 

print("Failed at i = %d!" % i) 

good = False 

if good: print("Each element reports itself as its root.") 

print("Merging:") 

for i,j in merges: 

OP_join(OP, i, j) 

print("Merged %d and %d." % (i, j)) 

print("Done merging.") 

print("Finding:") 

for i from 0 <= i < n: 

j = OP_find(OP, i) 

s = "%d -> %d"%(i, j) 

if i == j: 

s += ", root: size=%d, mcr=%d, rank=%d"%\ 

(OP.size[i], OP.mcr[i], OP.rank[i]) 

print(s) 

print("Allocating array to test merge_perm.") 

cdef int *gamma = <int *> sig_malloc( n * sizeof(int) ) 

if gamma is NULL: 

print("Allocation failed!") 

OP_dealloc(OP) 

return 

print("Allocation passed.") 

for i from 0 <= i < n: 

gamma[i] = perm[i] 

print("Merging permutation: %s" % perm) 

OP_merge_list_perm(OP, gamma) 

print("Done merging.") 

print("Finding:") 

for i from 0 <= i < n: 

j = OP_find(OP, i) 

s = "%d -> %d"%(i, j) 

if i == j: 

s += ", root: size=%d, mcr=%d, rank=%d"%\ 

(OP.size[i], OP.mcr[i], OP.rank[i]) 

print(s) 

print("Deallocating OrbitPartition.") 

sig_free(gamma) 

OP_dealloc(OP) 

print("Done.") 

  

  

# PartitionStack (PS) 

  

cdef inline PartitionStack *PS_new(int n, bint unit_partition): 

""" 

Allocate and return a pointer to a new PartitionStack of degree n. Returns a 

null pointer in the case of an allocation failure. 

""" 

cdef int i 

cdef PartitionStack *PS = <PartitionStack *> \ 

sig_malloc(sizeof(PartitionStack)) 

cdef int *int_array = <int *> sig_malloc( 2*n * sizeof(int) ) 

if PS is NULL or int_array is NULL: 

sig_free(PS) 

sig_free(int_array) 

return NULL 

PS.entries = int_array 

PS.levels = int_array + n 

PS.depth = 0 

PS.degree = n 

if unit_partition: 

PS_unit_partition(PS) 

return PS 

  

cdef void PS_unit_partition(PartitionStack *PS): 

""" 

Set partition stack to a single partition with a single cell. 

""" 

cdef int i, n = PS.degree 

PS.depth = 0 

for i from 0 <= i < n-1: 

PS.entries[i] = i 

PS.levels[i] = n 

PS.entries[n-1] = n-1 

PS.levels[n-1] = -1 

  

cdef inline PartitionStack *PS_copy(PartitionStack *PS): 

""" 

Allocate and return a pointer to a copy of PartitionStack PS. Returns a null 

pointer in the case of an allocation failure. 

""" 

cdef int i, n = PS.degree 

  

cdef PartitionStack *PS2 = <PartitionStack *> \ 

sig_malloc(sizeof(PartitionStack)) 

cdef int *int_array = <int *> sig_malloc( 2*n * sizeof(int) ) 

if PS2 is NULL or int_array is NULL: 

sig_free(PS2) 

sig_free(int_array) 

return NULL 

PS2.entries = int_array 

PS2.levels = int_array + n 

PS_copy_from_to(PS, PS2) 

return PS2 

  

cdef inline void PS_dealloc(PartitionStack *PS): 

if PS is not NULL: 

sig_free(PS.entries) 

sig_free(PS) 

  

cdef PartitionStack *PS_from_list(list L): 

""" 

Allocate and return a pointer to a PartitionStack representing L. Returns a 

null pointer in the case of an allocation failure. 

""" 

cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0 

for cell from 0 <= cell < num_cells: 

n += len(L[cell]) 

cdef PartitionStack *PS = PS_new(n, 0) 

if PS is NULL: 

return NULL 

for cell from 0 <= cell < num_cells: 

cur_len = len(L[cell]) 

for i from 0 <= i < cur_len: 

PS.entries[cur_start + i] = L[cell][i] 

PS.levels[cur_start + i] = n 

PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1) 

cur_start += cur_len 

PS.levels[cur_start-1] = 0 

if n > 0: 

PS.levels[n-1] = -1 

PS.depth = 0 

PS.degree = n 

return PS 

  

cdef PS_print(PartitionStack *PS): 

""" 

Print a visual representation of PS. 

""" 

cdef int i 

for i from 0 <= i <= PS.depth: 

PS_print_partition(PS, i) 

  

cdef PS_print_partition(PartitionStack *PS, int k): 

""" 

Print the partition at depth k. 

""" 

s = '(' 

for i from 0 <= i < PS.degree: 

s += str(PS.entries[i]) 

if PS.levels[i] <= k: 

s += '|' 

else: 

s += ' ' 

s = s[:-1] + ')' 

print(s) 

  

cdef int PS_first_smallest(PartitionStack *PS, bitset_t b, int *second_pos=NULL, 

PartitionRefinement_generic partn_ref_alg=None): 

""" 

Find the first occurrence of the smallest cell of size greater than one, 

which is admissible (checked by the function ``test_allowance``). 

Its entries are stored to b and its minimum element is returned. 

""" 

cdef int i = 0, j = 0, location = 0, n = PS.degree 

bitset_zero(b) 

while 1: 

if PS.levels[i] <= PS.depth: 

if i != j and n > i - j + 1 and (partn_ref_alg is None or  

partn_ref_alg._minimization_allowed_on_col(PS.entries[j])): 

n = i - j + 1 

location = j 

j = i + 1 

if PS.levels[i] == -1: break 

i += 1 

# location now points to the beginning of the first, smallest, 

# nontrivial cell 

i = location 

while 1: 

bitset_flip(b, PS.entries[i]) 

if PS.levels[i] <= PS.depth: break 

i += 1 

  

if second_pos != NULL: 

if n==2: 

second_pos[0] = PS.entries[location+1] 

else: 

second_pos[0] = -1 

  

  

return PS.entries[location] 

  

  

cdef int PS_all_new_cells(PartitionStack *PS, bitset_t** nonsingletons_ptr): 

""" 

Suppose a cell ``C`` was split into ``a`` components at ``PS.level``. 

Set the rows of the matrix ``nonsingletons_ptr`` to the first 

``a-1`` components of ``C`` for all those ``C``. 

Return the number of rows of ``nonsingletons_ptr``. 

""" 

cdef int beg=0, end, n = PS.degree, count=0, i, n_1 = n-1 

cdef bint non_unit_partition = False 

cdef bitset_t scratch 

bitset_init(scratch, n) 

cdef bitset_t* nonsingletons = nonsingletons_ptr[0] 

  

while beg < n_1: 

end = beg 

while end!=n and PS.levels[end] > PS.depth: 

end+=1 

  

if end != n: 

if PS.levels[end] == PS.depth: 

bitset_zero(scratch) 

for i from beg <= i <= end: 

bitset_set(scratch, PS.entries[i]) 

count +=1 

nonsingletons = <bitset_t*> sig_realloc(nonsingletons, count * sizeof(bitset_t)) 

if nonsingletons is NULL: 

raise MemoryError("Memory error in PS_all_new_cells") 

bitset_init(nonsingletons[count-1], n) 

bitset_copy(nonsingletons[count-1], scratch) 

else: 

if beg==0: 

nonsingletons = <bitset_t*> sig_realloc(nonsingletons, sizeof(bitset_t)) 

if nonsingletons is NULL: 

raise MemoryError("Memory error in PS_all_new_cells") 

bitset_init(nonsingletons[0], n) 

bitset_zero(scratch) 

bitset_complement(nonsingletons[0], scratch) 

count = 1 

beg = end+1 

nonsingletons_ptr[0] = nonsingletons 

return count 

  

cdef int PS_find_element(PartitionStack *PS, bitset_t b, int x) except -1: 

""" 

Find the cell containing x, store its entries to b and return the location 

of the beginning of the cell. 

""" 

cdef int i, location, n = PS.degree 

bitset_zero(b) 

for i from 0 <= i < n: 

if PS.entries[i] == x: 

location = i 

break 

else: 

raise ValueError("element not found") 

while location > 0 and PS.levels[location-1] > PS.depth: 

location -= 1 

i = 0 

while 1: 

bitset_set(b, PS.entries[location+i]) 

if PS.levels[location+i] <= PS.depth: 

break 

i += 1 

return location 

  

cdef list PS_singletons(PartitionStack * part): 

""" 

Return the list of all singletons in the PartitionStack. 

""" 

cdef list l = [] 

cdef int i 

  

if part.levels[0] <= part.depth: 

l.append(0) 

  

for i in range(1, part.degree): 

if part.levels[i] <= part.depth and part.levels[i - 1] <= part.depth: 

l.append(i) 

  

return l 

  

def PS_represent(partition, splits): 

""" 

Demonstration and testing. 

  

TESTS:: 

  

sage: from sage.groups.perm_gps.partn_ref.data_structures import PS_represent 

sage: PS_represent([[6],[3,4,8,7],[1,9,5],[0,2]], [6,1,8,2]) 

Allocating PartitionStack... 

Allocation passed: 

(0 1 2 3 4 5 6 7 8 9) 

Checking that entries are in order and correct level. 

Everything seems in order, deallocating. 

Deallocated. 

Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]]. 

PartitionStack's data: 

entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2] 

levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1] 

depth = 0, degree = 10 

(6|3 4 8 7|1 9 5|0 2) 

Checking PS_is_discrete: 

False 

Checking PS_num_cells: 

4 

Checking PS_is_mcr, min cell reps are: 

[6, 3, 1, 0] 

Checking PS_is_fixed, fixed elements are: 

[6] 

Copying PartitionStack: 

(6|3 4 8 7|1 9 5|0 2) 

Checking for consistency. 

Everything is consistent. 

Clearing copy: 

(0 3 4 8 7 1 9 5 6 2) 

Splitting point 6 from original: 

0 

(6|3 4 8 7|1 9 5|0 2) 

Splitting point 1 from original: 

5 

(6|3 4 8 7|1|5 9|0 2) 

Splitting point 8 from original: 

1 

(6|8|3 4 7|1|5 9|0 2) 

Splitting point 2 from original: 

8 

(6|8|3 4 7|1|5 9|2|0) 

Getting permutation from PS2->PS: 

[6, 1, 0, 8, 3, 9, 2, 7, 4, 5] 

Finding first smallest: 

Minimal element is 5, bitset is: 

0000010001 

Finding element 1: 

Location is: 5 

Bitset is: 

0100000000 

Deallocating PartitionStacks. 

Done. 

  

""" 

cdef int i, n = sum([len(cell) for cell in partition]) 

cdef int *gamma 

cdef bitset_t b 

print("Allocating PartitionStack...") 

cdef PartitionStack *PS = PS_new(n, 1) 

cdef PartitionStack *PS2 

if PS is NULL: 

print("Allocation failed!") 

return 

print("Allocation passed:") 

PS_print(PS) 

print("Checking that entries are in order and correct level.") 

good = True 

for i from 0 <= i < n-1: 

if not (PS.entries[i] == i and PS.levels[i] == n): 

print("Failed at i = %d!" % i) 

print(PS.entries[i], PS.levels[i], i, n) 

good = False 

if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1): 

print("Failed at i = %d!"%(n-1)) 

good = False 

if not PS.degree == n or not PS.depth == 0: 

print("Incorrect degree or depth!") 

good = False 

if good: print("Everything seems in order, deallocating.") 

PS_dealloc(PS) 

print("Deallocated.") 

print("Creating PartitionStack from partition %s."%partition) 

PS = PS_from_list(partition) 

print("PartitionStack's data:") 

print("entries -> %s"%[PS.entries[i] for i from 0 <= i < n]) 

print("levels -> %s"%[PS.levels[i] for i from 0 <= i < n]) 

print("depth = %d, degree = %d" % (PS.depth,PS.degree)) 

PS_print(PS) 

print("Checking PS_is_discrete:") 

print("True" if PS_is_discrete(PS) else "False") 

print("Checking PS_num_cells:") 

print(PS_num_cells(PS)) 

print("Checking PS_is_mcr, min cell reps are:") 

L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)] 

print(L) 

print("Checking PS_is_fixed, fixed elements are:") 

print([PS.entries[l] for l in L if PS_is_fixed(PS, l)]) 

print("Copying PartitionStack:") 

PS2 = PS_copy(PS) 

PS_print(PS2) 

print("Checking for consistency.") 

good = True 

for i from 0 <= i < n: 

if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]: 

print("Failed at i = %d!"%i) 

good = False 

if PS.degree != PS2.degree or PS.depth != PS2.depth: 

print("Failure with degree or depth!") 

good = False 

if good: 

print("Everything is consistent.") 

print("Clearing copy:") 

PS_clear(PS2) 

PS_print(PS2) 

for s in splits: 

print("Splitting point %d from original:"%s) 

print(PS_split_point(PS, s)) 

PS_print(PS) 

print("Getting permutation from PS2->PS:") 

gamma = <int *> sig_malloc(n * sizeof(int)) 

PS_get_perm_from(PS, PS2, gamma) 

print([gamma[i] for i from 0 <= i < n]) 

sig_free(gamma) 

print("Finding first smallest:") 

bitset_init(b, n) 

i = PS_first_smallest(PS, b) 

print("Minimal element is %d, bitset is:"%i) 

print(bitset_string(b)) 

bitset_free(b) 

print("Finding element 1:") 

bitset_init(b, n) 

print("Location is: {}".format(PS_find_element(PS, b, 1))) 

print("Bitset is:") 

print(bitset_string(b)) 

bitset_free(b) 

print("Deallocating PartitionStacks.") 

PS_dealloc(PS) 

PS_dealloc(PS2) 

print("Done.") 

  

  

# StabilizerChain (SC) 

  

cdef enum: 

default_num_gens = 8 

default_num_bits = 64 

  

cdef StabilizerChain *SC_new(int n, bint init_gens=True): 

""" 

Allocate and return a pointer to a new StabilizerChain of degree n. Returns 

a null pointer in the case of an allocation failure. 

""" 

cdef int i 

cdef StabilizerChain *SC = <StabilizerChain *> \ 

sig_calloc(1, sizeof(StabilizerChain)) 

cdef int *array1 

cdef int *array2 

cdef int *array3 

cdef bint mem_err = 0 

if SC is NULL: 

return NULL 

SC.degree = n 

SC.base_size = 0 

if n == 0: 

# All internal pointers have been initialized to NULL by sig_calloc 

return SC 

  

# first level allocations 

cdef int *int_array = <int *> sig_malloc( (3*n*n + 6*n + 1) * sizeof(int) ) 

cdef int **int_ptrs = <int **> sig_calloc( 5*n, sizeof(int *) ) 

SC.OP_scratch = OP_new(n) 

# bitset_init without the MemoryError: 

cdef long limbs = (default_num_bits - 1)/(8*sizeof(unsigned long)) + 1 

SC.gen_used.size = default_num_bits 

SC.gen_is_id.size = default_num_bits 

SC.gen_used.limbs = limbs 

SC.gen_is_id.limbs = limbs 

SC.gen_used.bits = <mp_limb_t*>sig_malloc(limbs * sizeof(mp_limb_t)) 

SC.gen_is_id.bits = <mp_limb_t*>sig_malloc(limbs * sizeof(mp_limb_t)) 

  

# check for allocation failures 

if int_array is NULL or int_ptrs is NULL or \ 

SC.gen_used.bits is NULL or SC.gen_is_id.bits is NULL or \ 

SC.OP_scratch is NULL: 

sig_free(int_array) 

sig_free(int_ptrs) 

SC_dealloc(SC) 

return NULL 

  

SC.gen_used.bits[limbs-1] = 0 

SC.gen_is_id.bits[limbs-1] = 0 

  

SC.orbit_sizes = int_array 

SC.num_gens = int_array + n 

SC.array_size = int_array + 2*n 

SC.perm_scratch = int_array + 3*n # perm_scratch is length 3*n+1 for sorting 

int_array += 6*n + 1 

  

SC.generators = int_ptrs 

SC.gen_inverses = int_ptrs + n 

SC.base_orbits = int_ptrs + 2*n 

SC.parents = int_ptrs + 3*n 

SC.labels = int_ptrs + 4*n 

for i from 0 <= i < n: 

SC.base_orbits[i] = int_array 

SC.parents[i] = int_array + n 

SC.labels[i] = int_array + 2*n 

int_array += 3*n 

  

# second level allocations 

if init_gens: 

for i from 0 <= i < n: 

SC.array_size[i] = default_num_gens 

SC.generators[i] = <int *> sig_malloc( default_num_gens*n * sizeof(int) ) 

SC.gen_inverses[i] = <int *> sig_malloc( default_num_gens*n * sizeof(int) ) 

if SC.generators[i] is NULL or SC.gen_inverses[i] is NULL: 

SC_dealloc(SC) 

return NULL 

  

return SC 

  

cdef inline int SC_realloc_gens(StabilizerChain *SC, int level, int size): 

""" 

Reallocate generator array at level `level` to size `size`. 

  

Returns 1 in case of an allocation failure. 

""" 

cdef int *temp 

cdef int n = SC.degree 

  

temp = <int *> sig_realloc( SC.generators[level], n * size * sizeof(int) ) 

if temp is NULL: return 1 

SC.generators[level] = temp 

  

temp = <int *> sig_realloc( SC.gen_inverses[level], n * size * sizeof(int) ) 

if temp is NULL: return 1 

SC.gen_inverses[level] = temp 

  

SC.array_size[level] = size 

return 0 

  

cdef inline void SC_dealloc(StabilizerChain *SC): 

cdef int i, n 

if SC is not NULL: 

n = SC.degree 

if SC.generators is not NULL: 

for i from 0 <= i < n: 

sig_free(SC.generators[i]) 

sig_free(SC.gen_inverses[i]) 

sig_free(SC.generators) # frees int_ptrs 

sig_free(SC.orbit_sizes) # frees int_array 

sig_free(SC.gen_used.bits) 

sig_free(SC.gen_is_id.bits) 

OP_dealloc(SC.OP_scratch) 

sig_free(SC) 

  

cdef StabilizerChain *SC_symmetric_group(int n): 

""" 

Returns a stabilizer chain for the symmetric group on {0, 1, ..., n-1}. 

  

Returns NULL in the case of an allocation failure. 

""" 

cdef int i, j, b 

cdef StabilizerChain *SC = SC_new(n, False) 

if SC is NULL: 

return NULL 

SC.base_size = n-1 

for i from 0 <= i < n-1: 

SC.array_size[i] = n-i-1 

SC.array_size[n-1] = default_num_gens 

for i from 0 <= i < n: 

SC.generators[i] = <int *> sig_malloc( SC.array_size[i]*n * sizeof(int) ) 

SC.gen_inverses[i] = <int *> sig_malloc( SC.array_size[i]*n * sizeof(int) ) 

if SC.generators[i] is NULL or SC.gen_inverses[i] is NULL: 

SC_dealloc(SC) 

return NULL 

cdef int *id_perm = SC.perm_scratch 

for i from 0 <= i < n: 

id_perm[i] = i 

for i from 0 <= i < n-1: 

b = i 

SC.orbit_sizes[i] = n-i 

SC.num_gens[i] = n-i-1 

for j from 0 <= j < i: 

SC.parents[i][j] = -1 

for j from 0 <= j < n-i: 

SC.base_orbits[i][j] = i+j 

SC.parents[i][i+j] = b 

SC.labels[i][i+j] = j 

for j from 0 <= j < n-i-1: 

#j-th generator sends i+j+1 to b 

memcpy(SC.generators[i] + n*j, id_perm, n * sizeof(int) ) 

SC.generators[i][n*j + i+j+1] = b 

SC.generators[i][n*j + b] = i+j+1 

memcpy(SC.gen_inverses[i] + n*j, SC.generators[i] + n*j, n * sizeof(int) ) 

return SC 

  

cdef StabilizerChain *SC_alternating_group(int n): 

""" 

Returns a stabilizer chain for the alternating group on {0, 1, ..., n-1}. 

  

Returns NULL in the case of an allocation failure. 

""" 

cdef int i, j, b 

cdef StabilizerChain *SC = SC_new(n, False) 

if SC is NULL: 

return NULL 

SC.base_size = n-2 

for i from 0 <= i < n-2: 

SC.array_size[i] = n-i-1 

SC.array_size[n-2] = default_num_gens 

SC.array_size[n-1] = default_num_gens 

for i from 0 <= i < n: 

SC.generators[i] = <int *> sig_malloc( SC.array_size[i]*n * sizeof(int) ) 

SC.gen_inverses[i] = <int *> sig_malloc( SC.array_size[i]*n * sizeof(int) ) 

if SC.generators[i] is NULL or SC.gen_inverses[i] is NULL: 

SC_dealloc(SC) 

return NULL 

cdef int *id_perm = SC.perm_scratch 

for i from 0 <= i < n: 

id_perm[i] = i 

for i from 0 <= i < n-2: 

b = i 

SC.orbit_sizes[i] = n-i 

SC.num_gens[i] = n-i-2 

for j from 0 <= j < i: 

SC.parents[i][j] = -1 

for j from 0 <= j < n-i: 

SC.base_orbits[i][j] = i+j 

SC.parents[i][i+j] = b 

SC.labels[i][i+j] = j 

SC.labels[i][n-1] = -(n-i-2) 

for j from 0 <= j < n-i-2: 

#j-th generator sends i+j+1 to b, i+j+2 to i+j+1, and b to i+j+2 

memcpy(SC.generators[i] + n*j, id_perm, n * sizeof(int) ) 

SC.generators[i][n*j + i+j+1] = b 

SC.generators[i][n*j + b ] = i+j+2 

SC.generators[i][n*j + i+j+2] = i+j+1 

SC_invert_perm(SC.gen_inverses[i] + n*j, SC.generators[i] + n*j, n) 

return SC 

  

cdef int SC_realloc_bitsets(StabilizerChain *SC, long size): 

""" 

If size is larger than current allocation, double the size of the bitsets 

until it is not. 

  

Returns 1 in case of an allocation failure. 

""" 

cdef unsigned long size_old = SC.gen_used.size 

if size <= size_old: return 0 

cdef unsigned long new_size = size_old 

while new_size < size: 

new_size *= 2 

cdef unsigned long limbs_old = SC.gen_used.limbs 

cdef long limbs = (new_size - 1)/(8*sizeof(unsigned long)) + 1 

cdef mp_limb_t *tmp = <mp_limb_t*> sig_realloc(SC.gen_used.bits, limbs * sizeof(mp_limb_t)) 

if tmp is not NULL: 

SC.gen_used.bits = tmp 

else: 

return 1 

tmp = <mp_limb_t*> sig_realloc(SC.gen_is_id.bits, limbs * sizeof(mp_limb_t)) 

if tmp is not NULL: 

SC.gen_is_id.bits = tmp 

else: 

return 1 

SC.gen_used.limbs = limbs 

SC.gen_is_id.limbs = limbs 

SC.gen_used.size = new_size 

SC.gen_is_id.size = new_size 

SC.gen_used.bits[size_old >> index_shift] &= limb_lower_bits_down(size_old) 

memset(SC.gen_used.bits + (size_old >> index_shift) + 1, 0, (limbs - (size_old >> index_shift) - 1) * sizeof(unsigned long)) 

SC.gen_is_id.bits[size_old >> index_shift] &= limb_lower_bits_down(size_old) 

memset(SC.gen_is_id.bits + (size_old >> index_shift) + 1, 0, (limbs - (size_old >> index_shift) - 1) * sizeof(unsigned long)) 

return 0 

  

cdef StabilizerChain *SC_copy(StabilizerChain *SC, int level): 

""" 

Creates a copy of the first `level` levels of SC. Must have 0 < level. 

  

Returns a null pointer in case of allocation failure. 

""" 

cdef int i, n = SC.degree 

cdef StabilizerChain *SCC = SC_new(n, False) 

if SCC is NULL: 

return NULL 

level = min(level, SC.base_size) 

for i from 0 <= i < level: 

SCC.generators[i] = <int *> sig_malloc( SC.array_size[i]*n * sizeof(int) ) 

SCC.gen_inverses[i] = <int *> sig_malloc( SC.array_size[i]*n * sizeof(int) ) 

if SCC.generators[i] is NULL or SCC.gen_inverses[i] is NULL: 

SC_dealloc(SCC) 

return NULL 

SCC.array_size[i] = SC.array_size[i] 

for i from level <= i < n: 

SCC.generators[i] = <int *> sig_malloc( default_num_gens*n * sizeof(int) ) 

SCC.gen_inverses[i] = <int *> sig_malloc( default_num_gens*n * sizeof(int) ) 

if SCC.generators[i] is NULL or SCC.gen_inverses[i] is NULL: 

SC_dealloc(SCC) 

return NULL 

SCC.array_size[i] = default_num_gens 

SC_copy_nomalloc(SCC, SC, level) # no chance for memory error here... 

return SCC 

  

cdef int SC_copy_nomalloc(StabilizerChain *SC_dest, StabilizerChain *SC, int level): 

cdef int i, n = SC.degree 

level = min(level, SC.base_size) 

SC_dest.base_size = level 

memcpy(SC_dest.orbit_sizes, SC.orbit_sizes, 2*n * sizeof(int) ) # copies orbit_sizes, num_gens 

memcpy(SC_dest.base_orbits[0], SC.base_orbits[0], 3*n*n * sizeof(int) ) # copies base_orbits, parents, labels 

for i from 0 <= i < level: 

if SC.num_gens[i] > SC_dest.array_size[i]: 

if SC_realloc_gens(SC_dest, i, max(SC.num_gens[i], 2*SC_dest.array_size[i])): 

return 1 

memcpy(SC_dest.generators[i], SC.generators[i], SC.num_gens[i]*n * sizeof(int) ) 

memcpy(SC_dest.gen_inverses[i], SC.gen_inverses[i], SC.num_gens[i]*n * sizeof(int) ) 

return 0 

  

cdef SC_print_level(StabilizerChain *SC, int level): 

cdef int i, j, n = SC.degree 

if level < SC.base_size: 

print('/ level {}'.format(level)) 

print('| orbit {}'.format([SC.base_orbits[level][i] for i from 0 <= i < SC.orbit_sizes[level]])) 

print('| parents {}'.format([SC.parents [level][i] for i from 0 <= i < n])) 

print('| labels {}'.format([SC.labels [level][i] for i from 0 <= i < n])) 

print('|') 

print('| generators {}'.format([[SC.generators [level][n*i + j] for j from 0 <= j < n] for i from 0 <= i < SC.num_gens[level]])) 

print('\ inverses {}'.format([[SC.gen_inverses[level][n*i + j] for j from 0 <= j < n] for i from 0 <= i < SC.num_gens[level]])) 

else: 

print('/ level {}'.format(level)) 

print('|') 

print('\ base_size {}'.format(SC.base_size)) 

  

cdef StabilizerChain *SC_new_base(StabilizerChain *SC, int *base, int base_len): 

""" 

Create a new stabilizer chain whose base starts with the given base, and 

which represents the same permutation group. Original StabilizerChain is 

unmodified. 

  

Use SC_cleanup to remove redundant base points. 

  

Returns a null pointer in case of an allocation failure. 

""" 

cdef StabilizerChain *NEW = SC_new(SC.degree) 

if NEW is NULL: 

return NULL 

if SC_new_base_nomalloc(NEW, SC, base, base_len): 

SC_dealloc(NEW) 

return NULL 

return NEW 

  

cdef int SC_new_base_nomalloc(StabilizerChain *SC_dest, StabilizerChain *SC, int *base, int base_len): 

cdef int i, n = SC.degree 

SC_dest.base_size = 0 

for i from 0 <= i < base_len: 

SC_add_base_point(SC_dest, base[i]) 

if SC_update(SC_dest, SC, 0): 

SC_dealloc(SC_dest) 

return 1 

return 0 

  

cdef int SC_update(StabilizerChain *dest, StabilizerChain *source, int level): 

cdef mpz_t src_order, dst_order 

cdef int *perm = dest.perm_scratch 

mpz_init(src_order) 

mpz_init(dst_order) 

SC_order(source, level, src_order) 

SC_order(dest, level, dst_order) 

cdef int i, first_moved, b 

while mpz_cmp(dst_order, src_order): 

SC_random_element(source, level, perm) 

i = level 

while i < dest.base_size: 

b = dest.base_orbits[i][0] 

if perm[b] != b: 

break 

i += 1 

else: 

for b from 0 <= b < dest.degree: 

if perm[b] != b: 

break 

else: 

continue 

SC_add_base_point(dest, b) 

first_moved = i 

for i from level <= i <= first_moved: 

if SC_insert_and_sift(dest, i, perm, 1, 0): # don't sift! 

mpz_clear(dst_order) 

mpz_clear(src_order) 

return 1 

SC_order(dest, level, dst_order) 

mpz_clear(src_order) 

mpz_clear(dst_order) 

return 0 

  

cdef StabilizerChain *SC_insert_base_point(StabilizerChain *SC, int level, int p): 

""" 

Insert the point ``p`` as a base point on level ``level``. Return a new 

StabilizerChain with this new base. Original StabilizerChain is unmodified. 

  

Use SC_cleanup to remove redundant base points. 

  

Returns a null pointer in case of an allocation failure. 

""" 

cdef int i, b, n = SC.degree 

cdef StabilizerChain *NEW 

if level > 0: 

NEW = SC_copy(SC, level) 

else: 

NEW = SC_new(n) 

if NEW is NULL: 

return NULL 

SC_add_base_point(NEW, p) 

for i from level <= i < SC.base_size: 

b = SC.base_orbits[i][0] 

if b != p: 

SC_add_base_point(NEW, b) 

if SC_update(NEW, SC, level): 

SC_dealloc(NEW) 

return NULL 

return NEW 

  

cdef int SC_insert_base_point_nomalloc(StabilizerChain *SC_dest, StabilizerChain *SC, int level, int p): 

cdef int i, b, n = SC.degree 

SC_copy_nomalloc(SC_dest, SC, level) 

SC_add_base_point(SC_dest, p) 

for i from level <= i < SC.base_size: 

b = SC.base_orbits[i][0] 

if b != p: 

SC_add_base_point(SC_dest, b) 

if SC_update(SC_dest, SC, level): 

return 1 

return 0 

  

cdef int SC_re_tree(StabilizerChain *SC, int level, int *perm, int x): 

""" 

Return values: 

0 - No errors. 

1 - Allocation failure. 

""" 

cdef int *gen 

cdef int *gen_inv 

cdef int i, b, gen_index, error, n = SC.degree 

  

# make sure we have room for the new generator: 

if SC.array_size[level] == SC.num_gens[level]: 

if SC_realloc_gens(SC, level, 2*SC.array_size[level]): 

return 1 

cdef int *new_gen = SC.generators [level] + n*SC.num_gens[level] 

cdef int *new_gen_inv = SC.gen_inverses[level] + n*SC.num_gens[level] 

  

# new generator is perm^(-1) * (path from x to base) (left to right composition) 

SC_invert_perm(new_gen, perm, n) 

SC_compose_up_to_base(SC, level, x, new_gen) 

SC_invert_perm(new_gen_inv, new_gen, n) 

SC.num_gens[level] += 1 

  

# now that we have our generators, regenerate the tree, breadth-first 

b = SC.base_orbits[level][0] 

for i from 0 <= i < n: 

SC.parents[level][i] = -1 

SC.parents[level][b] = b 

i = 0 

SC.orbit_sizes[level] = 1 

while i < SC.orbit_sizes[level]: 

x = SC.base_orbits[level][i] 

for gen_index from SC.num_gens[level] > gen_index >= 0: 

gen_inv = SC.gen_inverses[level] + n*gen_index 

SC_scan(SC, level, x, gen_index, gen_inv, 1) 

for gen_index from 0 <= gen_index < SC.num_gens[level]: 

gen = SC.generators [level] + n*gen_index 

SC_scan(SC, level, x, gen_index, gen, -1) 

i += 1 

return 0 

  

cdef int SC_sift(StabilizerChain *SC, int level, int x, int *gens, int num_gens, int *new_gens): 

""" 

Apply Schreier's subgroup lemma[1] as follows. Given a level, a point x, and 

a generator s, find the coset traversal element r coming from x. 

Try inserting rsR(rs)^(-1) (composing left to right) on level+1, where 

R(g) is the coset representative of g. 

  

level - which level we are currently working on 

x - the object representing r 

gens - points to the generators to sift by 

num_gens - how many of these there are 

new_gens - space of size at least num_gens*n for the sifted perms to go 

  

Returns 1 in case of an allocation failure. 

""" 

cdef int n = SC.degree 

if num_gens == 0: 

return 0 

  

# copy a representative taking base to the point x to each of these 

cdef int i 

cdef int *temp = SC.gen_inverses[level] + n*SC.num_gens[level] # one more scratch space 

# (available since num_gens > 0) 

cdef int *rep_inv = temp 

SC_identify(rep_inv, n) 

SC_compose_up_to_base(SC, level, x, rep_inv) 

SC_invert_perm(new_gens, rep_inv, n) 

for i from 0 < i < num_gens: 

memcpy(new_gens + n*i, new_gens, n*sizeof(int)) 

  

# post-compose each one with a generator 

for i from 0 <= i < num_gens: 

SC_mult_perms(new_gens + n*i, new_gens + n*i, gens + n*i, n) 

  

# for each one, look up the representative it now gives, and post-compose with that inverse 

cdef int y, b = SC.base_orbits[level][0] 

cdef int *perm 

cdef int *perm_rep_inv = temp 

cdef int j 

for i from 0 <= i < num_gens: 

perm = new_gens + n*i # this is now rs 

y = perm[b] 

SC_identify(perm_rep_inv, n) 

SC_compose_up_to_base(SC, level, y, perm_rep_inv) 

SC_mult_perms(perm, perm, perm_rep_inv, n) 

return SC_insert(SC, level+1, new_gens, num_gens) 

  

cdef int SC_insert_and_sift(StabilizerChain *SC, int level, int *pi, int num_perms, bint sift): 

cdef int i, j, b, n = SC.degree 

cdef int perm_gen_index 

cdef int max_orbit_size, max_orbit_place 

if sift: 

if SC_realloc_bitsets(SC, num_perms): 

return 1 

bitset_clear(&SC.gen_used) 

bitset_clear(&SC.gen_is_id) 

b = -1 

for perm_gen_index from 0 <= perm_gen_index < num_perms: 

for i from 0 <= i < n: 

if pi[n*perm_gen_index + i] != i: 

b = i 

break 

else: 

bitset_set(&SC.gen_is_id, perm_gen_index) 

if b != -1: break 

if b == -1: return 0 

if sift and level == SC.base_size: 

SC_add_base_point(SC, b) 

else: 

b = SC.base_orbits[level][0] 

# Now b is the base element at level `level` 

  

# Record the old orbit elements and the old generators (see sifting phase) 

cdef int old_num_gens = SC.num_gens[level] 

cdef int old_num_points = SC.orbit_sizes[level] 

  

# Add new points to the tree: 

cdef int x 

cdef int *perm 

cdef int start_over = 1 

cdef int error 

cdef int re_treed = 0 

while start_over: 

start_over = 0 

for i from 0 <= i < SC.orbit_sizes[level]: 

x = SC.base_orbits[level][i] 

for perm_gen_index from 0 <= perm_gen_index < num_perms: 

if sift and bitset_check(&SC.gen_is_id, perm_gen_index): continue 

perm = pi + n*perm_gen_index 

if SC.parents[level][perm[x]] == -1: 

# now we have an x which maps to a new point under perm, 

re_treed = 1 

if sift: bitset_set(&SC.gen_used, perm_gen_index) 

if SC_re_tree(SC, level, perm, x): 

return 1 

start_over = 1 # we must look anew 

break 

if start_over: break 

if not re_treed: continue 

for perm_gen_index from 0 <= perm_gen_index < old_num_gens: 

perm = SC.generators[level] + n*perm_gen_index 

if SC.parents[level][perm[x]] == -1: 

# now we have an x which maps to a new point under perm, 

if SC_re_tree(SC, level, perm, x): 

return 1 

start_over = 1 # we must look anew 

break 

if start_over: break 

for j from level < j < SC.base_size: 

for perm_gen_index from 0 <= perm_gen_index < SC.num_gens[j]: 

perm = SC.generators[j] + n*perm_gen_index 

if SC.parents[level][perm[x]] == -1: 

# now we have an x which maps to a new point under perm, 

if SC_re_tree(SC, level, perm, x): 

return 1 

start_over = 1 # we must look anew 

break 

if not sift: 

return 0 

  

# store the rest of the unused generators in the generator array, after the added ones 

cdef int new_size 

cdef int unused_gens = 0 

for perm_gen_index from 0 <= perm_gen_index < num_perms: 

if not bitset_check(&SC.gen_used, perm_gen_index) and not bitset_check(&SC.gen_is_id, perm_gen_index): 

unused_gens += 1 

if 2*(SC.num_gens[level] + unused_gens) > SC.array_size[level]: 

new_size = max(2*(SC.num_gens[level] + unused_gens), 2*SC.array_size[level]) 

if SC_realloc_gens(SC, level, new_size): 

return 1 

i = 0 

for perm_gen_index from 0 <= perm_gen_index < num_perms: 

if not bitset_check(&SC.gen_used, perm_gen_index) and not bitset_check(&SC.gen_is_id, perm_gen_index): 

memcpy(SC.generators[level] + n*(SC.num_gens[level]+i), pi + n*perm_gen_index, n*sizeof(int)) 

i += 1 

  

# Sift: 

cdef int *gens = SC.generators[level] 

cdef int total_gens = SC.num_gens[level] + unused_gens 

cdef int section, gens_in_section 

for i from 0 <= i < SC.orbit_sizes[level]: 

x = SC.base_orbits[level][i] 

section = 0 

while section*n < total_gens: 

gens_in_section = min(n, total_gens - section*n) 

if SC_sift(SC, level, x, gens + n*n*section, gens_in_section, gens + n*total_gens): 

return 1 

section += 1 

return 0 

  

cdef bint SC_is_giant(int n, int num_perms, int *perms, float p, bitset_t support): 

""" 

Test whether the group generated by the input permutations is a giant, i.e., 

the alternating or symmetric group. 

  

If the group is not a giant, this routine will return False. This could also 

indicate an allocation failure. 

  

If the group is a giant, this routine will return True with approximate 

probability p. It will set `support' to the support of the group in this 

case. Use bitset_len to get the size of support. 

  

The bitset `support' must be initialized. Must have 0 <= p < 1. 

  

Running time is roughly O(-ln(1-p)*n*ln(m)) where m <= n is the size of the 

support of the group. 

""" 

cdef int i, j, num_steps, m = 1, support_root 

cdef unsigned long q 

cdef int *gen 

cdef int *perm = <int *> sig_malloc(n*sizeof(int)) 

cdef OrbitPartition *OP = OP_new(n) 

if OP is NULL or perm is NULL: 

OP_dealloc(OP) 

sig_free(perm) 

return False 

  

# giants are transitive 

for j from 0 <= j < num_perms: 

gen = perms + n*j 

for i from 0 <= i < n: 

OP_join(OP, i, gen[i]) 

for i from 0 <= i < n: 

if OP.parent[i] == i: 

if OP.size[i] != 1: 

if m != 1: 

m = 1 

break 

else: 

m = OP.size[i] 

support_root = i 

if m == 1: 

OP_dealloc(OP) 

sig_free(perm) 

return False 

bitset_zero(support) 

for i from 0 <= i < n: 

if OP_find(OP, i) == support_root: 

bitset_set(support, i) 

  

# get a bit lost in the group, so our random elements are more random: 

SC_identify(perm, n) 

for i from 0 <= i < 10: 

SC_mult_perms(perm, perm, perms + n*(rand()%num_perms), n) 

  

# look for elements with cycles of prime length q, m/2 < q < m-2 

num_steps = <int> ceil(-log(1-p)*log(m)/log(2)) 

for j from 0 <= j < num_steps: 

OP_clear(OP) 

for i from 0 <= i < n: 

OP_join(OP, i, perm[i]) 

for i from 0 <= i < n: 

if OP.parent[i] == i: 

q = OP.size[i] 

if m < 2*q and q < m-2: 

if n_is_prime(q): 

sig_free(perm) 

OP_dealloc(OP) 

return True 

SC_mult_perms(perm, perm, perms + n*(rand()%num_perms), n) 

OP_dealloc(OP) 

sig_free(perm) 

return False 

  

def SC_test_list_perms(list L, int n, int limit, bint gap, bint limit_complain, bint test_contains): 

""" 

Test that the permutation group generated by list perms in L of degree n 

is of the correct order, by comparing with GAP. Don't test if the group is 

of size greater than limit. 

  

TESTS:: 

  

sage: from sage.groups.perm_gps.partn_ref.data_structures import SC_test_list_perms 

sage: limit = 10^7 

sage: def test_Sn_on_m_points(n, m, gap, contains): 

....: perm1 = [1,0] + list(range(2, m)) 

....: perm2 = [(i+1)%n for i in range( n )] + list(range(n,m)) 

....: SC_test_list_perms([perm1, perm2], m, limit, gap, 0, contains) 

sage: for i in range(2,9): 

....: test_Sn_on_m_points(i,i,1,0) 

sage: for i in range(2,9): 

....: test_Sn_on_m_points(i,i,0,1) 

sage: for i in range(2,9): # long time 

....: test_Sn_on_m_points(i,i,1,1) # long time 

sage: test_Sn_on_m_points(8,8,1,1) 

sage: def test_stab_chain_fns_1(n, gap, contains): 

....: perm1 = sum([[2*i+1,2*i] for i in range(n)], []) 

....: perm2 = [(i+1)%(2*n) for i in range( 2*n)] 

....: SC_test_list_perms([perm1, perm2], 2*n, limit, gap, 0, contains) 

sage: for n in range(1,11): 

....: test_stab_chain_fns_1(n, 1, 0) 

sage: for n in range(1,11): 

....: test_stab_chain_fns_1(n, 0, 1) 

sage: for n in range(1,9): # long time 

....: test_stab_chain_fns_1(n, 1, 1) # long time 

sage: test_stab_chain_fns_1(11, 1, 1) 

sage: def test_stab_chain_fns_2(n, gap, contains): 

....: perms = [] 

....: for p,e in factor(n): 

....: perm1 = [(p*(i//p)) + ((i+1)%p) for i in range(n)] 

....: perms.append(perm1) 

....: SC_test_list_perms(perms, n, limit, gap, 0, contains) 

sage: for n in range(2,11): 

....: test_stab_chain_fns_2(n, 1, 0) 

sage: for n in range(2,11): 

....: test_stab_chain_fns_2(n, 0, 1) 

sage: for n in range(2,11): # long time 

....: test_stab_chain_fns_2(n, 1, 1) # long time 

sage: test_stab_chain_fns_2(11, 1, 1) 

sage: def test_stab_chain_fns_3(n, gap, contains): 

....: perm1 = [(-i)%n for i in range( n )] 

....: perm2 = [(i+1)%n for i in range( n )] 

....: SC_test_list_perms([perm1, perm2], n, limit, gap, 0, contains) 

sage: for n in range(2,20): 

....: test_stab_chain_fns_3(n, 1, 0) 

sage: for n in range(2,20): 

....: test_stab_chain_fns_3(n, 0, 1) 

sage: for n in range(2,14): # long time 

....: test_stab_chain_fns_3(n, 1, 1) # long time 

sage: test_stab_chain_fns_3(20, 1, 1) 

sage: def test_stab_chain_fns_4(n, g, gap, contains): 

....: perms = [] 

....: for _ in range(g): 

....: perm = list(range(n)) 

....: shuffle(perm) 

....: perms.append(perm) 

....: SC_test_list_perms(perms, n, limit, gap, 0, contains) 

sage: for n in range(4,9): # long time 

....: test_stab_chain_fns_4(n, 1, 1, 0) # long time 

....: test_stab_chain_fns_4(n, 2, 1, 0) # long time 

....: test_stab_chain_fns_4(n, 2, 1, 0) # long time 

....: test_stab_chain_fns_4(n, 2, 1, 0) # long time 

....: test_stab_chain_fns_4(n, 2, 1, 0) # long time 

....: test_stab_chain_fns_4(n, 3, 1, 0) # long time 

sage: for n in range(4,9): 

....: test_stab_chain_fns_4(n, 1, 0, 1) 

....: for j in range(6): 

....: test_stab_chain_fns_4(n, 2, 0, 1) 

....: test_stab_chain_fns_4(n, 3, 0, 1) 

sage: for n in range(4,8): # long time 

....: test_stab_chain_fns_4(n, 1, 1, 1) # long time 

....: test_stab_chain_fns_4(n, 2, 1, 1) # long time 

....: test_stab_chain_fns_4(n, 2, 1, 1) # long time 

....: test_stab_chain_fns_4(n, 3, 1, 1) # long time 

sage: test_stab_chain_fns_4(8, 2, 1, 1) 

sage: def test_stab_chain_fns_5(n, gap, contains): 

....: perms = [] 

....: m = n//3 

....: perm1 = list(range(2*m)) 

....: shuffle(perm1) 

....: perm1 += list(range(2*m,n)) 

....: perm2 = list(range(m,n)) 

....: shuffle(perm2) 

....: perm2 = list(range(m)) + perm2 

....: SC_test_list_perms([perm1, perm2], n, limit, gap, 0, contains) 

sage: for n in [4..9]: # long time 

....: for _ in range(2): # long time 

....: test_stab_chain_fns_5(n, 1, 0) # long time 

sage: for n in [4..8]: # long time 

....: test_stab_chain_fns_5(n, 0, 1) # long time 

sage: for n in [4..9]: # long time 

....: test_stab_chain_fns_5(n, 1, 1) # long time 

sage: def random_perm(x): 

....: shuffle(x) 

....: return x 

sage: def test_stab_chain_fns_6(m,n,k, gap, contains): 

....: perms = [] 

....: for i in range(k): 

....: perm = sum([random_perm(list(range(i*(n//m),min(n,(i+1)*(n//m))))) for i in range(m)], []) 

....: perms.append(perm) 

....: SC_test_list_perms(perms, m*(n//m), limit, gap, 0, contains) 

sage: for m in range(2,9): # long time 

....: for n in range(m,3*m): # long time 

....: for k in range(1,3): # long time 

....: test_stab_chain_fns_6(m,n,k, 1, 0) # long time 

sage: for m in range(2,10): 

....: for n in range(m,4*m): 

....: for k in range(1,3): 

....: test_stab_chain_fns_6(m,n,k, 0, 1) 

sage: test_stab_chain_fns_6(10,20,2, 1, 1) 

sage: test_stab_chain_fns_6(8,16,2, 1, 1) 

sage: test_stab_chain_fns_6(6,36,2, 1, 1) 

sage: test_stab_chain_fns_6(4,40,3, 1, 1) 

sage: test_stab_chain_fns_6(4,40,2, 1, 1) 

sage: def test_stab_chain_fns_7(n, cop, gap, contains): 

....: perms = [] 

....: for i in range(0,n//2,2): 

....: p = list(range(n)) 

....: p[i] = i+1 

....: p[i+1] = i 

....: if cop: 

....: perms.append([c for c in p]) 

....: else: 

....: perms.append(p) 

....: SC_test_list_perms(perms, n, limit, gap, 0, contains) 

sage: for n in [6..14]: 

....: test_stab_chain_fns_7(n, 1, 1, 0) 

....: test_stab_chain_fns_7(n, 0, 1, 0) 

sage: for n in [6..30]: 

....: test_stab_chain_fns_7(n, 1, 0, 1) 

....: test_stab_chain_fns_7(n, 0, 0, 1) 

sage: for n in [6..14]: # long time 

....: test_stab_chain_fns_7(n, 1, 1, 1) # long time 

....: test_stab_chain_fns_7(n, 0, 1, 1) # long time 

sage: test_stab_chain_fns_7(20, 1, 1, 1) 

sage: test_stab_chain_fns_7(20, 0, 1, 1) 

  

""" 

if gap: 

from sage.groups.perm_gps.permgroup import PermutationGroup 

from sage.groups.perm_gps.permgroup_element import PermutationGroupElement 

from sage.misc.prandom import shuffle 

  

cdef StabilizerChain *SC 

cdef StabilizerChain *SCC 

cdef StabilizerChain *SCCC 

cdef StabilizerChain *SC_nb 

cdef Integer order, order2 

cdef int i, j, m, SC_says 

cdef bitset_t giant_support 

if gap: 

G = PermutationGroup([[i+1 for i in p] for p in L]) 

if G.order() > limit: 

if limit_complain: print('TOO BIG') 

return 

SC = SC_new(n) 

cdef int *perm = <int *>sig_malloc(n * (len(L)+3) * sizeof(int)) 

try: 

bitset_init(giant_support, n) 

except MemoryError: 

sig_free(perm) 

SC_dealloc(SC) 

raise MemoryError 

if perm is NULL or SC is NULL: 

bitset_free(giant_support) 

sig_free(perm) 

SC_dealloc(SC) 

raise MemoryError 

cdef int *perm2 = perm + n 

cdef int *perm3 = perm + 2*n 

for Lperm in L: 

for i from 0 <= i < n: 

perm[i] = Lperm[i] 

if SC_insert(SC, 0, perm, 1): 

bitset_free(giant_support) 

sig_free(perm) 

SC_dealloc(SC) 

raise MemoryError 

SCC = SC_copy(SC, n) 

SCCC = SC_insert_base_point(SC, 0, n-1) 

for i from 0 <= i < n: 

perm[i] = n-i-1 

SC_nb = SC_new_base(SC, perm, n) 

if SCC is NULL or SCCC is NULL or SC_nb is NULL: 

bitset_free(giant_support) 

sig_free(perm) 

SC_dealloc(SC) 

SC_dealloc(SCC) 

SC_dealloc(SCCC) 

SC_dealloc(SC_nb) 

raise MemoryError 

giant = False 

try: 

order = Integer(0) 

SC_order(SC,0,order.value) 

j = 0 

for Lperm in L: 

for i from 0 <= i < n: 

perm[n*j + i] = Lperm[i] 

j += 1 

if SC_is_giant(n, len(L), perm, 0.9, giant_support): 

giant = True 

m = bitset_len(giant_support) 

from sage.arith.all import factorial 

if not (order == factorial(m) or order == factorial(m)/2): 

print("SC_is_giant failed: %s %s"%(str(L), order)) 

raise AssertionError 

if order == factorial(n): 

SC_dealloc(SC) 

SC = SC_symmetric_group(n) 

SC_order(SC,0,order.value) 

if not order == factorial(n): 

print("SC_symmetric_group failed: %s %s"%(str(L), order)) 

raise AssertionError 

elif order == factorial(n)/2: 

SC_dealloc(SC) 

SC = SC_alternating_group(n) 

SC_order(SC,0,order.value) 

if not order == factorial(n)/2: 

print("SC_alternating_group failed: %s %s"%(str(L), order)) 

raise AssertionError 

order2 = Integer(0) 

SC_order(SCC,0,order2.value) 

if order != order2: 

print("FAIL {}".format(L)) 

print('SC_copy(n) does not agree with order of original {} {}'.format(order, order2)) 

raise AssertionError 

SC_order(SCCC,0,order2.value) 

if order != order2: 

print("FAIL {}".format(L)) 

print('does not agree with order of inserted base point {} {}'.format(order, order2)) 

raise AssertionError 

SC_order(SC_nb,0,order2.value) 

if order != order2: 

print("FAIL {}".format(L)) 

print('does not agree with order of new base {} {}'.format(order, order2)) 

raise AssertionError 

if test_contains: 

for i from 0 <= i < 3: 

SC_random_element(SC, 0, perm) 

if not SC_contains(SC, 0, perm, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=0) does not') 

raise AssertionError 

if not SC_contains(SC, 0, perm, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=1) does not') 

raise AssertionError 

if not SC_contains(SCC, 0, perm, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=0) does not on copy') 

raise AssertionError 

if not SC_contains(SCC, 0, perm, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=1) does not on copy') 

raise AssertionError 

if not SC_contains(SCCC, 0, perm, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=0) does not on inserted base point') 

raise AssertionError 

if not SC_contains(SCCC, 0, perm, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=1) does not on inserted base point') 

raise AssertionError 

if not SC_contains(SC_nb, 0, perm, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=0) does not on new base') 

raise AssertionError 

if not SC_contains(SC_nb, 0, perm, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element, SC_contains(modify=1) does not on new base') 

raise AssertionError 

SC_random_element(SCC, 0, perm2) 

if not SC_contains(SC, 0, perm2, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element of copy, SC_contains(modify=0) does not') 

raise AssertionError 

if not SC_contains(SC, 0, perm2, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element of copy, SC_contains(modify=1) does not') 

raise AssertionError 

SC_random_element(SCCC, 0, perm3) 

if not SC_contains(SC, 0, perm3, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element of inserted base point, SC_contains(modify=0) does not') 

raise AssertionError 

if not SC_contains(SC, 0, perm3, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element of inserted base point, SC_contains(modify=1) does not') 

raise AssertionError 

SC_random_element(SC_nb, 0, perm3) 

if not SC_contains(SC, 0, perm3, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element of new base, SC_contains(modify=0) does not') 

raise AssertionError 

if not SC_contains(SC, 0, perm3, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format([perm[ii] for ii in range(n)])) 

print('SC_random_element says it is an element of new base, SC_contains(modify=1) does not') 

raise AssertionError 

if gap: 

order = Integer(0) 

SC_order(SC,0,order.value) 

j = 0 

for Lperm in L: 

for i from 0 <= i < n: 

perm[n*j + i] = Lperm[i] 

j += 1 

if SC_is_giant(n, len(L), perm, 0.9, giant_support): 

from sage.arith.all import factorial 

m = bitset_len(giant_support) 

if order != factorial(m) and order != factorial(m)/2: 

print("SC_is_giant failed: %s %s"%(str(L), order)) 

raise AssertionError 

if order != G.order(): 

print("FAIL {}".format(L)) 

print('order was computed to be {}'.format(order)) 

print('GAP says it is {}'.format(G.order())) 

raise AssertionError 

if test_contains: 

for i from 0 <= i < 3: 

permy = G.random_element() 

for j from 0 <= j < n: 

perm[j] = permy(j+1)-1 

if not SC_contains(SC, 0, perm, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format(permy)) 

print('GAP says it is an element, SC_contains(modify=0) does not') 

raise AssertionError 

if not SC_contains(SC, 0, perm, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format(permy)) 

print('GAP says it is an element, SC_contains(modify=1) does not') 

raise AssertionError 

permy = list(xrange(1, n + 1)) 

shuffle(permy) 

gap_says = (PermutationGroupElement(permy) in G) 

for j from 0 <= j < n: 

perm[j] = permy[j]-1 

SC_says = SC_contains(SC, 0, perm, 0) 

if bool(SC_says) != bool(gap_says): 

print("FAIL {}".format(L)) 

print('element {}'.format(permy)) 

print('GAP says %d, SC_contains(modify=0) says %d'%(gap_says, SC_says)) 

raise AssertionError 

SC_says = SC_contains(SC, 0, perm, 1) 

if bool(SC_says) != bool(gap_says): 

print("FAIL {}".format(L)) 

print('element {}'.format(permy)) 

print('GAP says %d, SC_contains(modify=0) says %d'%(gap_says, SC_says)) 

raise AssertionError 

SC_random_element(SC, 0, perm) 

for j from 0 <= j < n: 

permy[j] = perm[j]+1 

gap_says = (PermutationGroupElement(permy) in G) 

if not SC_contains(SC, 0, perm, 0): 

print("FAIL {}".format(L)) 

print('element {}'.format(permy)) 

print('random element not contained in group, modify=false') 

raise AssertionError 

if not SC_contains(SC, 0, perm, 1): 

print("FAIL {}".format(L)) 

print('element {}'.format(permy)) 

print('random element not contained in group, modify=true') 

raise AssertionError 

if not gap_says: 

print("FAIL {}".format(L)) 

print('element {}'.format(permy)) 

print('random element not contained in group, according to GAP') 

raise AssertionError 

except Exception: 

if giant: 

print("detected giant!") 

raise 

finally: 

bitset_free(giant_support) 

sig_free(perm) 

SC_dealloc(SC) 

SC_dealloc(SCC) 

SC_dealloc(SCCC) 

SC_dealloc(SC_nb) 

  

# Functions 

  

cdef int sort_by_function(PartitionStack *PS, int start, int *degrees): 

""" 

A simple counting sort, given the degrees of vertices to a certain cell. 

  

INPUT: 

PS -- the partition stack to be checked 

start -- beginning index of the cell to be sorted 

degrees -- the values to be sorted by, must have extra scratch space for a 

total of 3*n+1 

  

""" 

cdef int n = PS.degree 

cdef int i, j, max, max_location 

cdef int *counts = degrees + n 

cdef int *output = degrees + 2*n + 1 

for i from 0 <= i <= n: 

counts[i] = 0 

i = 0 

while PS.levels[i+start] > PS.depth: 

counts[degrees[i]] += 1 

i += 1 

counts[degrees[i]] += 1 

# i+start is the right endpoint of the cell now 

max = counts[0] 

max_location = 0 

for j from 0 < j <= n: 

if counts[j] > max: 

max = counts[j] 

max_location = j 

counts[j] += counts[j - 1] 

for j from i >= j >= 0: 

counts[degrees[j]] -= 1 

output[counts[degrees[j]]] = PS.entries[start+j] 

max_location = counts[max_location]+start 

for j from 0 <= j <= i: 

PS.entries[start+j] = output[j] 

j = 1 

while j <= n and counts[j] <= i: 

if counts[j] > 0: 

PS.levels[start + counts[j] - 1] = PS.depth 

PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1) 

j += 1 

return max_location 

  

cdef int refine_by_orbits(PartitionStack *PS, StabilizerChain *SC, int *perm_stack, int *cells_to_refine_by, int *ctrb_len): 

""" 

Given a stabilizer chain SC, refine the partition stack PS so that each cell 

contains elements from at most one orbit, and sort the refined cells by 

their minimal cell representatives in the orbit of the group. 

""" 

cdef int start, level, gen_index, i, j, k, x, n = SC.degree 

cdef int *gen 

cdef int *min_cell_reps = SC.perm_scratch 

cdef int *perm = perm_stack + n*PS.depth 

cdef OrbitPartition *OP = SC.OP_scratch 

cdef int invariant = 1 

OP_clear(OP) 

for level from PS.depth <= level < SC.base_size: 

for gen_index from 0 <= gen_index < SC.num_gens[level]: 

gen = SC.generators[level] + gen_index*n 

for i from 0 <= i < n: 

OP_join(OP, i, gen[i]) 

start = 0 

while start < n: 

i = 0 

while 1: 

x = PS.entries[start+i] 

x = perm[x] 

min_cell_reps[i] = OP.mcr[OP_find(OP, x)] 

i += 1 

if PS.levels[start+i-1] <= PS.depth: 

break 

invariant += sort_by_function(PS, start, min_cell_reps) 

invariant += i 

# update cells_to_refine_by 

k = start 

for j from start <= j < start+i: 

if PS.levels[j] == PS.depth: 

k = j+1 

cells_to_refine_by[ctrb_len[0]] = k 

ctrb_len[0] += 1 

start += i 

return invariant 

  

cdef int compute_relabeling(StabilizerChain *group, StabilizerChain *scratch_group, 

int *permutation, int *relabeling): 

""" 

Technically, compute the INVERSE of the relabeling 

""" 

cdef int i, j, x, y, m, n = group.degree, orbit_element 

cdef int *scratch = group.perm_scratch 

if SC_new_base_nomalloc(scratch_group, group, permutation, n): 

return 1 

group = scratch_group 

SC_identify(relabeling, n) 

for i from 0 <= i < n: 

m = n 

for j from 0 <= j < group.orbit_sizes[i]: 

orbit_element = group.base_orbits[i][j] 

x = relabeling[orbit_element] 

if x < m: 

m = x 

y = orbit_element 

SC_invert_perm(scratch, relabeling, n) 

SC_compose_up_to_base(group, i, y, scratch) 

SC_invert_perm(relabeling, scratch, n) 

SC_invert_perm(scratch, relabeling, n) 

memcpy(relabeling, scratch, n*sizeof(int)) 

return 0