From ebf68d9e4ed11c3c7b0e7c561ae0853cd9e614d3 Mon Sep 17 00:00:00 2001
From: Frederik Wunden <frederik@familie-wunden.de>
Date: Thu, 5 Nov 2026 09:03:52 +0000
Subject: [PATCH 1/2] Add simplification, building and restriction functions

---
 lib/Wunden/BUILDING_GENERAL.md | 278 +++++++++++++++++++++++++
 lib/Wunden/README.md           |  78 +++++++
 lib/Wunden/RESTRICT_NONFREE.md | 270 ++++++++++++++++++++++++
 lib/Wunden/balance.gi          |  51 +++++
 lib/Wunden/building_general.gi | 152 ++++++++++++++
 lib/Wunden/decsWunden.gd       |  22 +-
 lib/Wunden/forget_action.gi    |  53 +++++
 lib/Wunden/initWunden.gi       |  12 +-
 lib/Wunden/lattice.gi          | 238 +++++++++++++++++++++
 lib/Wunden/reduce.gi           |  58 ++++++
 lib/Wunden/replace.gi          | 365 +++++++++++++++++++++++++++++++++
 lib/Wunden/restrict_nonfree.gi | 184 +++++++++++++++++
 lib/Wunden/simplify.gi         | 123 +++++++++++
 lib/Wunden/star.gi             |  73 +++++++
 14 files changed, 1955 insertions(+), 2 deletions(-)
 create mode 100644 lib/Wunden/BUILDING_GENERAL.md
 create mode 100644 lib/Wunden/README.md
 create mode 100644 lib/Wunden/RESTRICT_NONFREE.md
 create mode 100644 lib/Wunden/balance.gi
 create mode 100644 lib/Wunden/building_general.gi
 create mode 100644 lib/Wunden/forget_action.gi
 create mode 100644 lib/Wunden/lattice.gi
 create mode 100644 lib/Wunden/reduce.gi
 create mode 100644 lib/Wunden/replace.gi
 create mode 100644 lib/Wunden/restrict_nonfree.gi
 create mode 100644 lib/Wunden/simplify.gi
 create mode 100644 lib/Wunden/star.gi

diff --git a/lib/Wunden/BUILDING_GENERAL.md b/lib/Wunden/BUILDING_GENERAL.md
new file mode 100644
index 00000000..5b33acf8
--- /dev/null
+++ b/lib/Wunden/BUILDING_GENERAL.md
@@ -0,0 +1,278 @@
+> Package adaptation: implementations load automatically with HAP. References to the original standalone files below describe their provenance.
+
+# The building for SL(n, Z[1/p])
+
+[building_general.gi](building_general.gi) generalizes the construction to any
+integer $n\geq2$ and prime $p$. It returns a `HapNonFreeResolution`, with the
+same row-lattice convention and HAP interface as [original rank-three test fixture](../../tst/Wunden/building_sl3_reference.g).
+Neither `buildin.g` nor `CONSTRUCTION.md` is changed by this implementation.
+
+```gap
+LoadPackage("hap");
+B := BruhatTitsBuildingSL(6,2);;
+IsHapNonFreeResolution(B);          # true
+List([0..5],B!.dimension);          # [6,15,20,15,6,1]
+Length(B);                         # 5
+H := B!.stabilizer(2,1);;           # stabilizer of the first triangle
+GeneratorsOfGroup(H);
+IdentityMat(6) in H;               # true
+B!.boundary(2,1);
+B!.cellRepresentatives[3][1];      # its three vertex lattice bases
+```
+
+The files can be loaded together. They use different constructor and attribute
+names. For the old case, call `BruhatTitsBuildingSL(3,2)`; its cell numbering,
+orientations, boundaries, and stabilizer subgroups agree with the original.
+The generating lists need not be identical.
+
+## 1. The fundamental chamber
+
+The vertices are homothety classes of full-rank $\mathbb{Z}_p$-lattices in
+$\mathbb{Q}_p^n$. A row basis $a$ represents $L(a)=\mathbb{Z}_p^n a$.
+Choose
+
+$$
+a_i=\operatorname{diag}(\underbrace{1,\ldots,1}_{n-i},
+                        \underbrace{p,\ldots,p}_{i}),
+\qquad 0\leq i<n.
+$$
+
+Their lattices satisfy
+
+$$
+pL_0\subsetneq L_{n-1}\subsetneq\cdots\subsetneq L_1\subsetneq L_0.
+$$
+
+The vertices $v_i=[L_i]$ form an $(n-1)$-simplex, a fundamental chamber for
+$G=\mathrm{SL}_n(\mathbb{Z}[1/p])$. Every nonempty face supplies one cell-orbit
+representative. There are $\binom{n}{d+1}$ representatives of dimension $d$,
+and $2^n-1$ in total. In particular, $n=6$ requires only 63 representatives.
+As in the original file, these represent an infinite building, not a finite ball.
+
+The code generates the face lists with
+
+```gap
+cells := List([1..n],k -> Combinations([1..n],k));
+```
+
+GAP index `i+1` represents mathematical vertex $v_i$. Combinations are ordered
+lexicographically, and every face is oriented by increasing vertex index.
+
+The type $v_p(\det a)\bmod n$ is preserved by $G$. Every chamber contains one
+vertex of each type, so a cell stabilizer fixes all its vertices and preserves
+orientation. Thus the HAP orientation function `action(d,k,g)` is always 1.
+
+## 2. Boundaries and the HAP object
+
+For an oriented face, the boundary is the usual alternating sum:
+
+$$
+\partial[v_{i_0},\ldots,v_{i_d}]
+=\sum_{j=0}^d(-1)^j[v_{i_0},\ldots,\widehat{v_{i_j}},\ldots,v_{i_d}].
+$$
+
+Each face is already a chosen representative, so every attaching group
+element is the identity, at `B!.elts[1]`. Boundary terms have HAP format
+`[signedCellIndex,elementIndex]`; terms are returned in ascending face-index
+order to agree with `buildin.g`. Negative input cell indices negate the boundary.
+
+The fields `group`, `dimension`, `boundary`, `stabilizer`, `action`, `elts`,
+`homotopy`, `cellRepresentatives`, and `properties` have the same meaning as in
+the original construction. Stabilizers are now computed on their first use
+and cached. Their generators are actual rational matrices.
+
+HAP's left action corresponds to $g\cdot[L]=[Lg^{-1}]$ for row lattices.
+The group is generated by all $I+E_{ij}$, $i\ne j$, together with
+$\operatorname{diag}(p,p^{-1},1,\ldots,1)$. Conjugating elementary matrices
+and taking powers gives arbitrary coefficients in $\mathbb{Z}[1/p]$.
+For $n=2$, conjugation scales the off-diagonal coefficient by $p^2$;
+taking the $p$th power also supplies a coefficient $1/p$.
+
+## 3. Reducing a cell stabilizer to an integral flag stabilizer
+
+The stabilizer of $v_i$ is $a_i^{-1}\mathrm{SL}_n(\mathbb{Z})a_i$.
+A face stabilizer is the intersection of its vertex stabilizers.
+Instead of asking GAP to intersect infinite matrix groups, the code rotates
+the face to one containing $v_0$.
+
+Let $r$ have entry $p$ in position $(1,n)$ and entries 1 immediately below
+the diagonal; all its other entries are zero. It satisfies
+
+$$
+L_i r=L_{i+1}\quad(i<n-1),\qquad L_{n-1}r=pL_0.
+$$
+
+If the smallest vertex label in a face is $s$, right multiplication by
+$r^{-s}$ subtracts $s$ from all its labels without wraparound. The normalized
+face contains $v_0$, so its stabilizer is a subgroup of
+$\mathrm{SL}_n(\mathbb{Z})$. If that stabilizer is $H$, the original one is
+$r^{-s}Hr^s$. Although $r$ is not in $G$, it is in
+$\mathrm{GL}_n(\mathbb{Z}[1/p])$ and normalizes $G$ by conjugation.
+
+For a normalized nonzero label $j$, the lattice quotient
+$L_j/pL_0$ is the prefix subspace
+
+$$
+W_{n-j}=\langle e_1,\ldots,e_{n-j}\rangle\subset\mathbb{F}_p^n.
+$$
+
+The list `cuts` contains these dimensions. A row matrix preserves $W_c$
+exactly when its entries in positions $i\leq c<j$ vanish modulo $p$.
+For several cuts this gives a block lower-triangular finite flag stabilizer
+$P\leq\mathrm{SL}_n(\mathbb{F}_p)$. The integral stabilizer is its full preimage:
+
+$$
+H=\rho_p^{-1}(P),\qquad
+\rho_p:\mathrm{SL}_n(\mathbb{Z})\longrightarrow\mathrm{SL}_n(\mathbb{F}_p).
+$$
+
+For example, for $n=3$ the first edge gives `cuts=[2]`, the second gives
+`cuts=[1]`, and the chamber gives `cuts=[1,2]`. These are precisely the plane,
+line, and full flag used in the original code.
+
+## 4. Why the higher-rank stabilizer algorithm changes
+
+The original file uses Schreier generators by enumerating a finite flag orbit.
+That works well for $n=3,p=2$, where the largest orbit has 21 points. For the
+full flag in dimension 6 the orbit size is
+
+$$
+\prod_{i=1}^{6}\frac{p^i-1}{p-1}.
+$$
+
+Already at $p=2$ this is 615,195. Applying every elementary generator to all
+these flags would overwhelm the otherwise small chamber construction.
+
+For $n\geq3$ the new code instead builds
+
+$$
+H=\langle\text{generators of }\Gamma_n(p),
+        \ \text{integer lifts of generators of }P\rangle,
+\qquad \Gamma_n(p)=\ker\rho_p.
+$$
+
+This equality follows directly: any element of $H$ has the same finite image
+as a word in the chosen lifts, and dividing by that word leaves an element
+of the kernel. It is essential to include the kernel; the finite flag group
+alone does not describe the infinite stabilizer.
+
+### 4.1 Principal congruence generators
+
+For $n\geq3$, the code uses the generating set in Proposition 1.10 of
+[Detinko–Flannery–Hulpke, *Algorithms for arithmetic groups with the congruence
+subgroup property*](https://arxiv.org/pdf/1906.10423), attributed there to
+Sury–Venkataramana. Put $t_{ij}(p)=I+pE_{ij}$ and
+
+$$
+\Sigma=\{I\}\ \cup\ \{\text{coordinate transpositions}\}
+\ \cup\ \{I-2E_{kk}-2E_{k+1,k+1}+E_{k+1,k}:1\leq k<n\}.
+$$
+
+The matrices $c^{-1}t_{ij}(p)c$, for $i<j$ and $c\in\Sigma$, generate
+$\Gamma_n(p)$. `principalGenerators` implements this formula and removes
+duplicates. The formula applies to $p=2$ as well as odd primes, but is not
+used for $n=2$. Its number of candidate generators depends only on $n$.
+Coordinate transpositions have determinant $-1$, which is harmless:
+conjugation still preserves determinant 1 and the congruence condition.
+
+### 4.2 Lifting the finite flag group
+
+The finite block lower-triangular group is generated by:
+
+- Elementary matrices $I+E_{ij}$ in the permitted positions. Their powers
+  provide all coefficients because the field is the prime field.
+- The determinant-one diagonal torus.
+
+The permitted positions are exactly those for which no cut satisfies
+$i\leq c<j$. This includes both directions within a diagonal block and
+the lower off-block entries. Elementary matrices inside each block generate
+its special linear part. Adding the diagonal torus gives all possible block
+determinants with total product 1.
+
+For odd $p$, choose a primitive residue $u$ modulo $p$, and choose an integer
+$v$ with $uv\equiv1\pmod{p^2}$. On coordinates $(i,n)$ the code uses the block
+
+$$
+\begin{pmatrix}u&p\\(uv-1)/p&v\end{pmatrix}.
+$$
+
+Its determinant is exactly 1, and its reduction is
+$\operatorname{diag}(u,u^{-1})$. The lower off-diagonal entry is divisible by
+$p$ because the inverse was computed modulo $p^2$, not just $p$.
+Using these blocks for $i=1,\ldots,n-1$ generates the entire finite diagonal
+torus. When $p=2$ that torus is trivial and needs no generators.
+
+## 5. Rank two: Schreier's lemma remains explicit
+
+For $n=2$ the chamber is an edge. Its two vertex stabilizers are conjugates
+of $\mathrm{SL}_2(\mathbb{Z})$. The edge stabilizer is
+
+$$
+\{g\in\mathrm{SL}_2(\mathbb{Z}):g_{12}\equiv0\pmod p\}.
+$$
+
+It fixes the row line $\langle(1,0)\rangle$ modulo $p$. Representatives for
+its $p+1$ right cosets are
+
+$$
+r_t=\begin{pmatrix}1&t\\0&1\end{pmatrix}\quad(0\leq t<p),
+\qquad r_\infty=\begin{pmatrix}0&1\\-1&0\end{pmatrix}.
+$$
+
+The first rows represent the projective points $[1:t]$ and $[0:1]$.
+For each representative $r$ and each elementary generator $s$, the code
+identifies the projective point represented by the first row of $rs$.
+If its chosen representative is $r'$, then $rs(r')^{-1}$ fixes the initial
+line. Schreier's lemma says these matrices generate the whole stabilizer.
+This is the same idea as in `buildin.g`, with an explicit transversal rather
+than a search for new orbit points. The computation has $2(p+1)$ transitions.
+
+## 6. Membership, caching, and size limits
+
+`DeclareAttribute` and `InstallMethod` are used as before, now with the distinct
+attribute name `BTBuildingLatticeBases`. For a rational determinant-one matrix,
+membership in a cell stabilizer is tested by checking that $aga^{-1}$ is
+integral for every vertex basis $a$ of the cell. Since those bases are
+$p$-integral diagonal matrices, the test also implies that $g$ has entries in
+$\mathbb{Z}[1/p]$.
+
+The principal congruence generating set is computed once per building object.
+Integral flag generating lists are cached by their cuts. Each requested cell
+stabilizer is then conjugated, given its membership attribute, and cached.
+The constructor does not enumerate finite flag orbits for $n\geq3$.
+
+There is no hard upper bound of 6 in the code, but the number of representative
+faces is exponential in $n$. The test suite exercises $2\leq n\leq6$, with
+multiple primes including odd primes. Large $p$ increases integer arithmetic
+costs and, in rank two, the explicit Schreier transversal size. Operations
+you subsequently ask GAP to perform on infinite stabilizers may have their own
+costs; an efficient constructor does not make every such operation efficient.
+
+## 7. Validation and limits
+
+Run:
+
+```sh
+gap -q -T -c 'LoadPackage("hap"); ReadPackage("hap", "tst/Wunden/test_building_general.g"); QUIT;'
+```
+
+For the package installation on this machine:
+
+```sh
+gap -q -T -c 'LoadPackage("hap"); ReadPackage("hap", "tst/Wunden/test_building_general.g"); QUIT;'
+```
+
+Checks cover all cell stabilizers for each tested pair, determinant and lattice
+preservation, face-stabilizer inclusions, signed boundaries, boundary squared
+zero, augmentation, cache reuse, and negative membership examples. They also
+check finite parabolic image orders in smaller dimensions and for every tested
+full chamber, and the full trace-zero image of the congruence generators modulo
+$p^2$. The case $(3,2)$ is compared directly with the original file.
+
+These finite checks supplement the generation theorem; agreement of finite
+images alone would not prove equality of infinite subgroups.
+
+As before, `homotopy=fail`: there is no explicit contracting homotopy.
+Stabilizers are infinite, so HAP routines requiring finite stabilizers do not
+apply directly. The construction models Moss's building $B_n(p)$, not his
+larger space $Y$ involving the real symmetric space.
diff --git a/lib/Wunden/README.md b/lib/Wunden/README.md
new file mode 100644
index 00000000..cc70c3e4
--- /dev/null
+++ b/lib/Wunden/README.md
@@ -0,0 +1,78 @@
+# Equivariant G-complex simplification
+
+Load HAP in a fresh GAP session and use `ReduceGComplex(L)` to return a
+reduced complex. The input uses the HAP G-complex interface: `dimension`,
+`boundary`, `stabilizer`, `action`, `elts`, and `group`.
+
+The algorithm repeatedly replaces balanced stars that pass the implemented
+disk-collapse criterion. `HomeomorphicToDisk` implements that criterion;
+it is not a general algorithm for recognizing arbitrary topological disks.
+The replacement does not construct a contracting homotopy (`homotopy = fail`).
+
+Progress messages are off by default. Enable them with:
+
+```gap
+SetInfoLevel(InfoWunden, 1);
+```
+
+Declarations are in `decsWunden.gd`; `initWunden.gi` loads the seven modules
+in dependency order. 
+
+Run the small regression suites in a fresh GAP session:
+
+```gap
+LoadPackage("hap");
+ReadPackage("hap", "tst/Wunden/core.g");
+ReadPackage("hap", "tst/Wunden/integration.g");
+```
+
+Optional, slower homology regressions:
+
+```gap
+ReadPackage("hap", "tst/Wunden/sl3zs.g");
+ReadPackage("hap", "tst/Wunden/bianchi.g");
+```
+
+These test scripts assign global example variables; use a disposable session.
+The examples are regression checks, not a proof for all input complexes.
+
+## Bruhat--Tits buildings and restriction
+
+```gap
+LoadPackage("hap");
+B := BruhatTitsBuildingSL(3, 2);;
+# K must be a finite-index subgroup of B!.group.
+S := RestrictNonFreeResolution(B, K);;
+```
+
+`BruhatTitsBuildingSL(n,p)` requires an integer `n >= 2` and prime `p`.
+It supplies the chamber cell orbits, signed boundaries, and infinite cell
+stabilizers for SL(n,Z[1/p]). Stabilizer membership is tested by preservation
+of the stored lattice bases (`BTBuildingLatticeBases`). No contracting
+homotopy is supplied. Routines requiring finite stabilizers do not apply.
+
+`RestrictNonFreeResolution(R,K)` requires finite index, terminating subgroup
+membership tests, and the same group-element representation for K and R.
+The input must have a finite nonnegative recorded length. Infinite-index
+inputs may cause enumeration not to terminate; subgroup inclusion is only
+checked automatically when the ambient group is known finite.
+
+The output includes `toOriginal`, `fromOriginal`, `restrictionCells`, and
+`cosetRepresentatives`. It transports an available contracting homotopy and
+preserves `fail` when none is available. Restriction may extend R!.elts and
+populate caches. `RestrictionCellStabilizerData` supports exact membership
+in the restricted cell stabilizers.
+
+Run the regression suites in a disposable session:
+
+```gap
+ReadPackage("hap", "tst/Wunden/test_building_general.g");
+ReadPackage("hap", "tst/Wunden/test_restrict_nonfree.g");
+```
+
+The building suite reads `building_sl3_reference.g`, a test-only copy of the
+original rank-three constructor, for independent comparison. That constructor
+is not loaded as part of the package implementation. The tests assign global
+names and are currently run manually rather than through HAP's test runner.
+Do not also read the original standalone implementation scripts in the same
+session: their global assignments conflict with the package declarations.
diff --git a/lib/Wunden/RESTRICT_NONFREE.md b/lib/Wunden/RESTRICT_NONFREE.md
new file mode 100644
index 00000000..e3add19b
--- /dev/null
+++ b/lib/Wunden/RESTRICT_NONFREE.md
@@ -0,0 +1,270 @@
+> Package adaptation: implementations load automatically with HAP. References to the original standalone files below describe their provenance.
+
+# Restricting a non-free resolution using left cosets
+
+[restrict_nonfree.gi](restrict_nonfree.gi) supplies:
+
+```gap
+LoadPackage("hap");
+S := RestrictNonFreeResolution(R,K);;
+```
+
+Here `R` is a `HapNonFreeResolution` for `G := R!.group`, and `K` is a
+**finite-index subgroup of G**, represented by elements of the same ambient
+group. The output is a `HapNonFreeResolution` for `K` on the same underlying
+chain complex. HAP and its dependencies must be available.
+
+For a building constructed earlier, the call is simply
+
+```gap
+# K is your finite-index subgroup of B63!.group.
+S := RestrictNonFreeResolution(B63,K);;
+List([0..Length(S)],S!.dimension);
+S!.boundary(1,1);
+S!.stabilizer(1,1);
+```
+
+The implementation is separate from the building constructors and does not
+change `buildin.g`, `building_general.g`, or `CONSTRUCTION.md`.
+
+## 1. Why changing `R!.group` is not enough
+
+Suppose an original representative cell is $\sigma$, with stabilizer
+$H=\operatorname{Stab}_G(\sigma)$. Its translates form one $G$-orbit, but
+may form several $K$-orbits. Each new orbit needs a representative, a new
+stabilizer, and boundaries expressed in the new representatives.
+
+For example, let $S_3$ act on a filled triangle. There is one orbit in each
+dimension. Restricting to the subgroup generated by a transposition gives
+two vertex orbits, two edge orbits, and one triangle orbit. Restricting
+instead to the trivial subgroup gives three vertices, three edges, and
+one triangle. The underlying triangle is the same in each case.
+
+In general the new cell counts are neither the old counts nor necessarily
+the old counts multiplied by $[G:K]$.
+
+## 2. Left cosets throughout
+
+The procedure enumerates the finite set of **left cosets**
+
+$$
+G/K=\{xK:x\in G\}.
+$$
+
+The list `S!.cosetRepresentatives` stores matrices or group elements $x$,
+not coset objects. To recognize a coset, the code uses
+
+$$
+yK=xK\quad\Longleftrightarrow\quad x^{-1}y\in K.
+$$
+
+Enumeration starts with $K$ and repeatedly applies the generators of $G$
+by **left multiplication**: $xK\mapsto gxK$. Only subgroup membership
+and group arithmetic are needed. Positive generators suffice because
+they act as permutations of a finite coset set, so their inverse actions
+are positive powers.
+
+For an original stabilizer $H$, the code then computes its orbits on
+$G/K$, again by left multiplication. An orbit represented by $xK$
+corresponds to the new cell representative
+
+$$
+\tau=x^{-1}\sigma.
+$$
+
+The inverse here is essential. If $yK=hxK$ for $h\in H$, then
+$y=hxk$ for some $k\in K$, and
+$y^{-1}\sigma=k^{-1}x^{-1}\sigma$ up to the orientation sign of $h$.
+Thus the two cosets give cells in the same $K$-orbit. Conversely, cells
+in the same $K$-orbit give cosets in the same $H$-orbit.
+
+Each original cell can split into a different number of new orbits.
+The code records the new indices explicitly rather than assuming a common
+number of translates for all original cells.
+
+## 3. Stabilizers via Schreier's lemma
+
+The new cell stabilizer is
+
+$$
+J=\operatorname{Stab}_K(x^{-1}\sigma)
+ =K\cap x^{-1}Hx.
+$$
+
+It is unnecessary to ask GAP to intersect two infinite groups. Instead,
+Schreier's lemma computes generators from the finite $H$-orbit of $xK$.
+
+For every visited coset $T_qK$, the code stores a witness $w_q\in H$ with
+
+$$
+w_qxK=T_qK.
+$$
+
+For a generator $a\in H$, suppose $aT_qK=T_rK$. Then
+
+$$
+z=w_r^{-1}aw_q
+$$
+
+fixes the coset $xK$. These elements generate its stabilizer
+$H\cap xKx^{-1}$. Conjugating them by $x^{-1}$ gives generators
+$x^{-1}zx$ of the required subgroup $J$.
+
+This works whether $H$ is finite or infinite, provided its generators and
+the required membership tests are available. The new stabilizers are actual
+GAP subgroups of `K`, with explicit generators.
+
+An attribute `RestrictionCellStabilizerData` stores $K$, $H$, and the
+translate $t=x^{-1}$. The installed membership method then uses the exact test
+
+$$
+g\in J\quad\Longleftrightarrow\quad
+g\in K\ \text{and}\ t^{-1}gt\in H.
+$$
+
+As in the building code, `DeclareAttribute` supplies the getter, setter and
+tester, while `InstallMethod` provides membership for groups carrying that
+attribute. This avoids having to rediscover a membership algorithm for each
+new infinite matrix subgroup.
+
+## 4. Rewriting boundaries, including orientations
+
+HAP represents a chain term by `[signedCellIndex,elementIndex]`.
+An original term $y\sigma$ must be expressed as a subgroup translate of
+one of the new representatives.
+
+Find the left coset $y^{-1}K$. Suppose its orbit has base $xK$ and stored
+witness $w$, so that
+
+$$
+y^{-1}K=wxK.
+$$
+
+Consequently $k=ywx\in K$, and
+
+$$
+y=kx^{-1}w^{-1}.
+$$
+
+If $\chi_\sigma:H\to\{1,-1\}$ is the original stabilizer's orientation
+character, the chain term becomes
+
+$$
+y\sigma=\chi_\sigma(w^{-1})\,k\,(x^{-1}\sigma).
+$$
+
+The code obtains this sign from `R!.action` evaluated on the stabilizer
+element $w^{-1}$. Ignoring this sign would give wrong boundaries for cells
+whose stabilizers reverse orientation. The building constructors have
+trivial orientation characters, but this restriction procedure does not
+assume that.
+
+For the representative $\tau=t\sigma$, compute
+$\partial\tau=t\partial\sigma$ and rewrite every resulting term by the
+formula above. All group elements in the returned boundary belong to `K`.
+The orientation character of $J$ is transported by conjugation:
+
+$$
+\chi_\tau(k)=\chi_\sigma(t^{-1}kt),\qquad k\in J.
+$$
+
+For elements outside a cell stabilizer, `S!.action` uses the corresponding
+conjugated call to the input's orientation function. The character on the
+stabilizer is the part needed for the chain-module relations.
+
+## 5. Chain conversion and contracting homotopies
+
+The output includes two conversion functions:
+
+```gap
+S!.fromOriginal(n,word);  # R-chain -> S-chain in degree n
+S!.toOriginal(n,word);    # S-chain -> R-chain in degree n
+```
+
+Both take a list of HAP terms. For example,
+
+```gap
+w := S!.toOriginal(1,[[1,1]]);;
+S!.fromOriginal(1,w);
+```
+
+The resulting words represent the same chain; they need not be identical
+lists, since stabilizers give relations between different translated terms.
+The routines preserve terms rather than combining all cancellations.
+
+If `R!.homotopy = fail`, the output also has `homotopy = fail`.
+If an explicit contracting homotopy is present, it is transported through
+these conversions. It remains available in the degrees where the original
+homotopy is defined, including its augmentation component if supplied.
+No new contracting homotopy is invented.
+
+## 6. Additional output fields
+
+| Field | Meaning |
+| --- | --- |
+| `S!.originalResolution` | The input R |
+| `S!.cosetRepresentatives` | Representatives x of the left cosets xK |
+| `S!.restrictionCells[n+1][j]` | Data for the jth new representative in degree n |
+| `.originalCell` in that record | Its original cell index |
+| `.translate` | The group element t such that the new cell is t times the original |
+| `.stabilizer` | Its subgroup J of K |
+
+Cell data are constructed by degree on demand. Call `S!.dimension(n)` before
+inspecting `S!.restrictionCells[n+1]` directly. The usual HAP functions
+`boundary`, `stabilizer`, and `action` also prepare the needed degree.
+
+The input's recorded length and characteristic are preserved.
+`reduced` is conservatively set to false. Geometric fields specific to a
+particular constructor are not copied: the new representatives are described
+by their original cell and translate instead.
+
+## 7. Preconditions and computation costs
+
+- `K` must actually be a finite-index subgroup of `R!.group`, in the same
+  element representation. For finite ambient groups this is checked with
+  `IsSubgroup`. For infinite groups it is the caller's precondition.
+  Different identity elements are rejected before coset enumeration. In the
+  locally installed HAP data, `ContractibleGcomplex("SL3Z")` uses 6-by-6
+  matrices, whereas `BruhatTitsBuildingSL(3,p)` uses 3-by-3 matrices. Those
+  two objects cannot be combined directly: an explicit representation map
+  is needed to transport the subgroup or resolution first. Naming the same
+  abstract group does not make the matrix representations interchangeable.
+- Membership in `K`, and in the original cell stabilizers, must be computable
+  by GAP. The procedure can use specialized membership methods such as those
+  already installed for the building stabilizers; it cannot automatically
+  solve membership for every arbitrary infinite matrix group.
+- The ambient group and original stabilizers must supply finite generating
+  lists. The resolution must have a finite recorded length and finitely many
+  generators in each recorded degree.
+- When a needed element is absent from `R!.elts`, it is appended so the
+  original orientation/homotopy functions can refer to it. Existing indices
+  are preserved. Thus `R!.elts` may grow; the input's group, boundaries and
+  cell orbit data are not replaced.
+
+The initial enumeration visits $m=[G:K]$ left cosets. Coset lookup uses a
+linear scan with membership tests, so this implementation is intended for
+manageable finite indices. For each original cell, its stabilizer acts on
+at most these same $m$ cosets. Finitely many cosets guarantee termination of
+the orbit searches, assuming the underlying group operations terminate.
+Infinite-index inputs are outside the contract and may enumerate forever.
+
+## 8. Tests
+
+Run [test_restrict_nonfree.g](../../tst/Wunden/test_restrict_nonfree.g):
+
+```sh
+gap -q -T -c 'LoadPackage("hap"); ReadPackage("hap", "tst/Wunden/test_restrict_nonfree.g"); QUIT;'
+```
+
+The finite tests restrict the triangle's $S_3$ action to the whole group,
+the trivial subgroup, the alternating subgroup, and a nonnormal order-two
+subgroup. They check orbit counts, actual generated stabilizers, left-coset
+uniqueness, orientation signs, boundary transport, boundary squared zero,
+and the transported contraction. Repeating restriction checks the case where
+different old cells split into different numbers of orbits.
+
+The infinite-group tests restrict the buildings for
+$\mathrm{SL}_2(\mathbb{Z}[1/2])$ and
+$\mathrm{SL}_3(\mathbb{Z}[1/3])$ to finite-index congruence subgroups. They
+check subgroup membership, lattice preservation, face inclusions, finite
+congruence images, and conversion of translated cells.
diff --git a/lib/Wunden/balance.gi b/lib/Wunden/balance.gi
new file mode 100644
index 00000000..f4bf6138
--- /dev/null
+++ b/lib/Wunden/balance.gi
@@ -0,0 +1,51 @@
+# A star has a unique center, so its setwise stabilizer is the center's
+# stabilizer. Check balance using subgroup generators and cell cosets;
+# L!.elts is only a table of elements used in boundaries, not the whole group.
+InstallGlobalFunction(PreIsGBalanced, function(L,e,n,k)
+local U, StabU, d, orbit, cells, first, stab, h, cell, g;
+
+U:=kStar(L,e,n,k);
+StabU:=L!.stabilizer(n,e);
+for d in [2..k+1] do
+    for orbit in Set(List(U[d],cell->cell[1])) do
+        cells:=Filtered(U[d],cell->cell[1]=orbit);
+        first:=cells[1][2];
+        stab:=cells[1][3];
+
+        # An element carrying h*f to l*f lies in l*Stab(f)*h^-1.
+        # Every such element must fix the center, otherwise the translated
+        # star overlaps this star without being equal to it.
+        for g in GeneratorsOfGroup(stab) do
+            if not first*g*Inverse(first) in StabU then
+                return [false,StabU];
+            fi;
+        od;
+        for cell in cells do
+            if not cell[2]*Inverse(first) in StabU then
+                return [false,StabU];
+            fi;
+        od;
+
+        # Each center-stabilizer generator must permute these finitely
+        # many translated cells. With left translates, equality means
+        # l^-1*g*h lies in the stabilizer of the orbit representative.
+        for g in GeneratorsOfGroup(StabU) do
+            for cell in cells do
+                h:=cell[2];
+                if not ForAny(cells,other->Inverse(other[2])*g*h in stab) then
+                    return [false,StabU];
+                fi;
+            od;
+        od;
+    od;
+od;
+return [true,StabU];
+end);
+
+InstallGlobalFunction(IsGBalanced, function(L,e,n,k)
+return PreIsGBalanced(L,e,n,k)[1];
+end);
+
+InstallGlobalFunction(Stabstar, function(L,e,n,k)
+return PreIsGBalanced(L,e,n,k)[2];
+end);
diff --git a/lib/Wunden/building_general.gi b/lib/Wunden/building_general.gi
new file mode 100644
index 00000000..116bc5f9
--- /dev/null
+++ b/lib/Wunden/building_general.gi
@@ -0,0 +1,152 @@
+# B_n(p) for SL(n,Z[1/p]), using the conventions of buildin.g.
+# After LoadPackage("hap"); B := BruhatTitsBuildingSL(6,2);
+
+# Exact membership in a cell stabilizer; no enumeration of infinite groups.
+InstallMethod(\in, [IsMatrix, IsGroup and HasBTBuildingLatticeBases],
+    function(g,h)
+        local bases, n;
+        bases := BTBuildingLatticeBases(h); n := Length(bases[1]);
+        return Length(g)=n and ForAll(g,r -> Length(r)=n)
+            and ForAll(Concatenation(g),IsRat) and DeterminantMat(g)=1
+            and ForAll(bases,a -> ForAll(Concatenation(a*g*a^-1),IsInt));
+    end);
+
+InstallGlobalFunction(BruhatTitsBuildingSL, function(n,p)
+    local one, elementary, pairs, gens, a, cells, rotation, powers,
+          kernel, torus, flagCache, stabilizers, principalGenerators,
+          sl2EdgeGenerators, flagGenerators, dimension, checkCell,
+          boundary, stabilizer, u, v, i, g, diagonal;
+    if not IsInt(n) or n < 2 then Error("n must be an integer >= 2"); fi;
+    if not IsInt(p) or not IsPrimeInt(p) then Error("p must be prime"); fi;
+    one := IdentityMat(n);
+    elementary := function(i,j,x)
+        local e;
+        e := StructuralCopy(one); e[i][j] := x; return e;
+    end;
+    pairs := Filtered(Tuples([1..n],2),ij -> ij[1]<>ij[2]);
+    gens := List(pairs,ij -> elementary(ij[1],ij[2],1));
+
+    # Row bases a_i=diag(1,...,1,p,...,p), with i trailing p's.
+    # Degree d representatives are the (d+1)-subsets, oriented increasingly.
+    a := List([0..n-1],i -> DiagonalMat(
+        Concatenation(List([1..n-i],j -> 1),List([1..i],j -> p))));
+    cells := List([1..n],k -> Combinations([1..n],k));
+    rotation := NullMat(n,n); rotation[1][n] := p;
+    for i in [2..n] do rotation[i][i-1] := 1; od;
+    powers := List([0..n-1],i -> rotation^i);
+    kernel := fail; flagCache := rec();
+    stabilizers := List([1..n],i -> []);
+
+    # Generators of Gamma_n(p), n>=3. Sury--Venkataramana's formula,
+    # as stated in Detinko--Flannery--Hulpke, Proposition 1.10:
+    # https://arxiv.org/pdf/1906.10423 (see BUILDING_GENERAL.md).
+    principalGenerators := function()
+        local sigma, upper, ij, c, k;
+        sigma := [one];
+        for ij in Combinations([1..n],2) do
+            c := StructuralCopy(one);
+            c[ij[1]] := one[ij[2]]; c[ij[2]] := one[ij[1]];
+            Add(sigma,c);                    # coordinate transposition
+        od;
+        for k in [1..n-1] do
+            c := StructuralCopy(one);
+            c[k][k] := -1; c[k+1][k+1] := -1; c[k+1][k] := 1;
+            Add(sigma,c);
+        od;
+        upper := List(Combinations([1..n],2),ij -> elementary(ij[1],ij[2],p));
+        return Set(Concatenation(List(sigma,c -> List(upper,t -> c^-1*t*c))));
+    end;
+
+    # Lifts of generators of the determinant-one diagonal torus mod p.
+    # The 2x2 block has determinant 1 over Z, not merely modulo p.
+    torus := [];
+    if p > 2 and n > 2 then
+        u := PrimitiveRootMod(p); v := PowerModInt(u,-1,p^2);
+        for i in [1..n-1] do
+            g := StructuralCopy(one);
+            g[i][i] := u; g[i][n] := p;
+            g[n][i] := (u*v-1)/p; g[n][n] := v;
+            Add(torus,g);
+        od;
+    fi;
+
+    # Rank 2: Schreier generators for the row-line stabilizer in SL(2,Z).
+    # Transversal rows represent [1,t] (0<=t<p) and [0,1].
+    sl2EdgeGenerators := function()
+        local reps, result, r, s, x, pos;
+        reps := List([0..p-1],t -> [[1,t],[0,1]]);
+        Add(reps,[[0,1],[-1,0]]); result := [];
+        for r in reps do
+            for s in gens do
+                x := r*s;
+                if x[1][1] mod p=0 then pos := p+1;
+                else pos := (x[1][2]*PowerModInt(x[1][1],-1,p)) mod p+1; fi;
+                x := x*reps[pos]^-1;
+                if x<>one then AddSet(result,x); fi;
+            od;
+        od;
+        return result;
+    end;
+
+    # Full preimage of a standard row-flag stabilizer in SL(n,F_p).
+    # cuts lists the dimensions of the prefix subspaces in that flag.
+    flagGenerators := function(cuts)
+        local key, allowed;
+        if IsEmpty(cuts) then return gens; fi;
+        key := String(cuts);
+        if not IsBound(flagCache.(key)) then
+            if n=2 then flagCache.(key) := sl2EdgeGenerators();
+            else
+                if kernel=fail then kernel := principalGenerators(); fi;
+                allowed := Filtered(pairs,ij ->
+                    ForAll(cuts,d -> not (ij[1]<=d and d<ij[2])));
+                flagCache.(key) := Union(kernel,torus,
+                    List(allowed,ij -> elementary(ij[1],ij[2],1)));
+            fi;
+        fi;
+        return flagCache.(key);
+    end;
+
+    dimension := function(d)
+        if IsInt(d) and d in [0..n-1] then return Length(cells[d+1]); fi;
+        return 0;
+    end;
+    checkCell := function(d,k)
+        if not IsInt(k) or k=0 or AbsInt(k)>dimension(d) then
+            Error("Cell index out of range");
+        fi;
+    end;
+    boundary := function(d,k)
+        local c;
+        checkCell(d,k);
+        if d=0 then return []; fi;
+        c := cells[d+1][AbsInt(k)];
+        return Reversed(List([1..Length(c)],j ->
+            [SignInt(k)*(-1)^(j-1)*Position(cells[d],Difference(c,[c[j]])),1]));
+    end;
+    stabilizer := function(d,k)
+        local c, shift, cuts, r, h;
+        checkCell(d,k); k := AbsInt(k);
+        if not IsBound(stabilizers[d+1][k]) then
+            c := cells[d+1][k]; shift := c[1]-1;
+            # Rotate the first vertex to v0, compute its mod-p flag,
+            # then conjugate the integral stabilizer back to the given cell.
+            cuts := Set(List(c{[2..Length(c)]},j -> n-(j-c[1])));
+            r := powers[shift+1];
+            h := Group(List(flagGenerators(cuts),g -> r^-1*g*r));
+            SetBTBuildingLatticeBases(h,a{c}); SetIsFinite(h,false);
+            stabilizers[d+1][k] := h;
+        fi;
+        return stabilizers[d+1][k];
+    end;
+    diagonal := List([1..n],i -> 1); diagonal[1] := p; diagonal[2] := 1/p;
+    return Objectify(HapNonFreeResolution,rec(
+        group := Group(Concatenation(gens,[DiagonalMat(diagonal)])),
+        elts := [one], dimension := dimension, boundary := boundary,
+        stabilizer := stabilizer,
+        action := function(d,k,g) return 1; end,
+        homotopy := fail,
+        cellRepresentatives := List(cells,level -> List(level,c -> a{c})),
+        properties := [["length",n-1],["characteristic",0],["type","resolution"]]
+    ));
+end);
diff --git a/lib/Wunden/decsWunden.gd b/lib/Wunden/decsWunden.gd
index 36bb8918..81d83a81 100644
--- a/lib/Wunden/decsWunden.gd
+++ b/lib/Wunden/decsWunden.gd
@@ -1 +1,21 @@
-#This file contains all function and operation declarations
+# Declarations for equivariant G-complex simplification.
+DeclareInfoClass("InfoWunden");
+DeclareGlobalFunction("kStar");
+DeclareGlobalFunction("ForgetG");
+DeclareGlobalFunction("PreIsGBalanced");
+DeclareGlobalFunction("IsGBalanced");
+DeclareGlobalFunction("Stabstar");
+DeclareGlobalFunction("EntireBoundary");
+DeclareGlobalFunction("Cobound");
+DeclareGlobalFunction("ShouldRemove");
+DeclareGlobalFunction("JoinCells");
+DeclareGlobalFunction("SimplifyLattice");
+DeclareGlobalFunction("HomeomorphicToDisk");
+DeclareGlobalFunction("Replace2");
+DeclareGlobalFunction("ReduceGComplex");
+
+# Bruhat--Tits buildings and finite-index restriction.
+DeclareAttribute("BTBuildingLatticeBases", IsGroup);
+DeclareAttribute("RestrictionCellStabilizerData", IsGroup);
+DeclareGlobalFunction("BruhatTitsBuildingSL");
+DeclareGlobalFunction("RestrictNonFreeResolution");
diff --git a/lib/Wunden/forget_action.gi b/lib/Wunden/forget_action.gi
new file mode 100644
index 00000000..1ec3bf44
--- /dev/null
+++ b/lib/Wunden/forget_action.gi
@@ -0,0 +1,53 @@
+InstallGlobalFunction(ForgetG, function(L,U,n)
+local FL, Y, k, d, cell, boundarycell, dcell, orbits, grouped, LO, i, g, g1, h, gh, t, r, z, Le, Stab, glist, j;
+
+k:=Length(U)-1;
+FL:=[];
+FL[1]:=[[1,0]];
+
+for d in [2..k+1] do
+    Le:=Length(U[d]);
+    FL[d]:=List([1..Le], i -> []);
+    for r in [1..Le] do
+        cell:=U[d][r];
+        t:=0; #t counts how many cells are in the boundary
+        FL[d][r][1]:=1; #Platzhalter
+        boundarycell:=L!.boundary(d+n-1,cell[1]);
+        orbits:=Set(List(U[d-1],x->x[1]));
+        grouped:=List(orbits, y -> Filtered(U[d-1], x -> x[1] = y));
+        for dcell in boundarycell do
+            Stab:=L!.stabilizer(d+n-2,AbsoluteValue(dcell[1]));
+            h:=L!.elts[dcell[2]];
+            g:=cell[2];
+            gh:=g*h;
+            LO:=Length(orbits);
+            z:=0;
+            for i in [1..LO] do
+                if i>1 then
+                    z:=z+Length(grouped[i-1]); #counts number of cell
+                fi;
+                if AbsoluteValue(dcell[1])=orbits[i] then
+                    glist:=(List(grouped[i],y->y[2]));
+                    for j in [1..Length(glist)] do
+                        g1:=glist[j];
+                        if Inverse(g1)*gh in Stab then #getauscht von g1*Inverse(gh) zu Inverse(g1)*gh
+                            Add(FL[d][r],z+j);
+                            t:=t+1;
+                        fi;
+                    od;
+                fi;
+            od;
+        od;
+        FL[d][r][1]:=t;
+    od;
+od;
+
+Add(FL,[]);
+#Y:=RegularCWComplex(FL);
+#cobnd:=Y!.coboundaries;
+
+#return [FL,cobnd];
+return FL;
+end);
+###################################
+
diff --git a/lib/Wunden/initWunden.gi b/lib/Wunden/initWunden.gi
index 1dafc032..f978021c 100644
--- a/lib/Wunden/initWunden.gi
+++ b/lib/Wunden/initWunden.gi
@@ -1 +1,11 @@
-#This file will read all of the .gi files
+# Load implementations in dependency order.
+ReadPackageHap("lib/Wunden/star.gi");
+ReadPackageHap("lib/Wunden/forget_action.gi");
+ReadPackageHap("lib/Wunden/balance.gi");
+ReadPackageHap("lib/Wunden/lattice.gi");
+ReadPackageHap("lib/Wunden/simplify.gi");
+ReadPackageHap("lib/Wunden/replace.gi");
+ReadPackageHap("lib/Wunden/reduce.gi");
+
+ReadPackageHap("lib/Wunden/building_general.gi");
+ReadPackageHap("lib/Wunden/restrict_nonfree.gi");
diff --git a/lib/Wunden/lattice.gi b/lib/Wunden/lattice.gi
new file mode 100644
index 00000000..d9d06de6
--- /dev/null
+++ b/lib/Wunden/lattice.gi
@@ -0,0 +1,238 @@
+InstallGlobalFunction(EntireBoundary, function(bnd,d,n)
+local V, U, D, tmp, v , cells;
+
+if d=0 then return []; fi;
+
+D:=d+1;
+V:=StructuralCopy(bnd[D][n]);
+V:=V{[2..Length(V)]};
+V:=SSortedList(V);
+cells:=List(V,i->[D-2,i]);
+D:=D-1;
+
+while D>1 do
+tmp:=[];
+for v in V do
+U:=StructuralCopy(bnd[D][v]);
+U:=U{[2..Length(U)]};
+Append(tmp,U);
+Append(cells,List(U,i->[D-2,i]));   #Changed D-1 to D-2  (May 2022)
+od;
+tmp:=SSortedList(tmp);
+cells:=SSortedList(cells);
+V:=tmp;
+D:=D-1;
+od;
+
+return cells;
+end);
+
+InstallGlobalFunction(Cobound, function(bnd)
+local dim, cobnd, b, k, i, j, n;
+
+dim:=PositionProperty(bnd,IsEmpty)-2;
+cobnd:=[];; #Cobnd[n+1] contains the info on n-cells.
+for n in [0..dim] do
+  #k:=2*(n+1)+1;#k:=1+2^(n+1);
+  cobnd[n+1]:=List(bnd[n+1],i->[0]);
+  for j in [1..Length(bnd[n+2])] do
+    b:=bnd[n+2][j];
+    k:=Length(b);
+    for i in b{[2..k]} do
+      cobnd[n+1][i][1]:=cobnd[n+1][i][1]+1;
+      Add(cobnd[n+1][i],j);
+    od;
+  od;
+od;
+cobnd[dim+1]:=List(bnd[dim+1],a->[0]);
+return cobnd;
+end);
+
+InstallGlobalFunction(ShouldRemove, function(FL,d1,n)
+local d2, cob, V1, V2, V3, bnd, cobnd, Y;
+
+bnd := FL;
+
+#Y := RegularCWComplex(FL);
+cobnd := Cobound(bnd);
+
+if n > Length(bnd[d1]) then
+    return false;
+fi;
+
+if not (cobnd[d1][n][1] = 2 and bnd[d1][n][1] > 0) then
+    return false;
+fi;
+
+d2 := d1+1;
+cob := 1*cobnd[d1][n];
+
+if Length(cob) < 3 then
+    return false;
+fi;
+
+if d2 > Length(cobnd) then
+    return false;
+fi;
+
+if cob[2] > Length(cobnd[d2]) or cob[3] > Length(cobnd[d2]) then
+    return false;
+fi;
+
+if not SortedList(cobnd[d2][cob[2]]) = SortedList(cobnd[d2][cob[3]]) then
+    return false;
+fi;
+
+
+# Print("d1=", d1, " n=", n, "\n");
+# Print("d2=", d2, "\n");
+# Print("cob=", cob, "\n");
+# Print("Length(bnd)=", Length(bnd), "\n");
+# Print("Length(bnd[d1])=", Length(bnd[d1]), "\n");
+# Print("Length(bnd[d2])=", Length(bnd[d2]), "\n");
+
+V1:=EntireBoundary(bnd,d1-1,n);
+V2:=EntireBoundary(bnd,d1,cobnd[d1][n][2]);
+V3:=EntireBoundary(bnd,d1,cobnd[d1][n][3]);
+
+if not 1+Size(V1) = Size(Intersection(V2,V3)) then
+    return false;
+fi;
+
+return true;
+end);
+
+
+InstallGlobalFunction(JoinCells, function(FL,d1,n)
+local bnd, cobnd, cob, d, d2, d3, m, t, poss, s, pos, u, perm, cnt, F, x, i;
+
+d2:=d1+1;
+d3:=d1+2;
+bnd:=FL;
+cobnd:=Cobound(bnd);
+#cobnd:=Y!.coboundaries;
+cob:=1*cobnd[d1][n];
+##
+##REDUCE COBOUNDARIES OF BOUNDARIES OF nTH CELL
+if d1>1 then
+d:=d1-1;;
+for m in bnd[d1][n]{[2..Length(bnd[d1][n])]} do
+    t:=1*cobnd[d][m]{[2..Length(cobnd[d][m])]};;
+    poss:=Position(t,n);;
+    Remove(t,poss);;
+    cobnd[d][m]:=Concatenation([Length(t)],t);;
+od;;
+fi;;
+#Print("cbnd ",cobnd,"\n");
+##
+##COBOUNDARIES OF BOUNDARIES REDUCED
+
+##
+##REMOVE nTH CELL, ITS COBOUNDARY, AND ADJUST ITS PRESENCE IN BOUNDARIES
+##OF ITS COBOUNDARIES
+bnd[d1][n]:=[0];;
+cobnd[d1][n]:=[0];;
+#Print("cobnd",cobnd,"\n");
+#Print("bnd",bnd,"\n");
+
+s:=bnd[d2][cob[2]];;
+s:=s{[2..Length(s)]};;
+pos:=Position(s,n);;
+Remove(s,pos);;
+
+
+t:=bnd[d2][cob[3]];;
+t:=t{[2..Length(t)]};;
+pos:=Position(t,n);;
+Remove(t,pos);;
+
+
+u:=SortedList(Concatenation(s,t));;
+
+if IsEmpty(u) then
+    bnd[d2][cob[2]]:=[1,0];;
+else
+    bnd[d2][cob[2]]:=Concatenation([Length(u)],u);;
+fi;;
+
+##
+##nTH CELL AND ITS PRESENCE REMOVED
+
+
+##
+##FOR SECOND CELL OF DIMENSION n+1  REDUCE THE BOUNDARIES OF ITS 
+##COBOUNDARIES
+
+for m in cobnd[d2][cob[3]]{[2..Length(cobnd[d2][cob[3]])]} do
+t:=bnd[d3][m]{[2..Length(bnd[d3][m])]};;
+
+pos:=Position(t,cob[3]);;
+Remove(t,pos);;
+bnd[d3][m]:=Concatenation([Length(t)],t);;
+od;;
+
+
+##
+##SECOND CELL  BOUNDARIES OF ITS COBOUNDARIES REDUCED
+
+##
+##REMOVE PRESENCE OF SECOND CELL IN COBOUNDARIES OF ITS BOUNDARIES
+for m in bnd[d2][cob[3]]{[2..Length(bnd[d2][cob[3]])]} do
+if cobnd[d1][m][1]>0 then
+
+t:=cobnd[d1][m]{[2..Length(cobnd[d1][m])]};;
+
+pos:=Position(t,cob[3]);;
+t[pos]:=cob[2];;
+t:=SSortedList(t);;
+cobnd[d1][m]:=Concatenation([Length(t)],t);;
+
+fi;;
+od;;
+
+
+bnd[d2][cob[3]]:=[0];;
+
+cobnd[d2][cob[3]]:=[0];;
+
+perm:=[];  ##perm[d][i] will be the new position of
+           ##the old i-th cell of dimension d;
+
+for d in [1..Length(bnd)] do
+perm[d]:=[];
+cnt:=0;
+for n in [1..Length(bnd[d])] do
+if bnd[d][n][1]=0 then cnt:=cnt+1;
+else
+perm[d][n]:=n-cnt;
+fi;
+od;
+od;
+
+
+
+for d in [1..Length(bnd)] do
+F:=1*Filtered([1..Length(bnd[d])],i->not bnd[d][i][1]=0);
+bnd[d]:=1*((bnd[d]){F});
+
+if d>1 then
+d1:=d-1;
+for x in bnd[d] do
+    if x[2]>0 then
+        for i in [2..Length(x)] do
+            x[i]:=1*perm[d1][x[i]];
+        od;
+    fi;
+od;
+fi;
+od;
+
+#if all lower-dimensional cells have been deleted, move all the remaining cells downwards
+while Length(bnd) > 0 and IsEmpty(bnd[1]) do
+    Remove(bnd, 1);
+od;
+
+return bnd;
+end);
+#################################
+
diff --git a/lib/Wunden/reduce.gi b/lib/Wunden/reduce.gi
new file mode 100644
index 00000000..4973f5ae
--- /dev/null
+++ b/lib/Wunden/reduce.gi
@@ -0,0 +1,58 @@
+#Iteratively simplifies a G-complex L by repeatedly finding a cell e in some
+#dimension n whose k-star is
+#   (i)  G-balanced      -- IsGBalanced(L,e,n,k) = true, and
+#   (ii) a disk           -- HomeomorphicToDisk(L,e,n,k) = true,
+#and replacing that star by a single new cell via Replace2(L,e,n,k).
+#
+#The search restarts from dimension 0 after every successful replacement,
+#since Replace2 renumbers cells across several dimensions at once and it is
+#not safe to keep old indices around. This continues until a full pass over
+#the whole complex finds nothing left to replace.
+#
+
+InstallGlobalFunction(ReduceGComplex, function(L)
+local CurrentTopDim, Progress, n, e, k, dim, kmax;
+ 
+# L:=arg[1];
+# if Length(arg)>1 and arg[2]>0 then
+#     MaxK:=arg[2];
+# else
+#     MaxK:=infinity;
+# fi;
+ 
+#Highest dimension in which L currently has any cells.
+CurrentTopDim:=function()
+local d;
+    d:=0;
+    while L!.dimension(d+1) > 0 do
+        d:=d+1;
+    od;
+    return d;
+end;
+ 
+repeat
+    Progress:=false;
+    n:=0;
+    while n <= CurrentTopDim() and not Progress do
+        dim:=L!.dimension(n);
+        e:=1;
+        while e <= dim and not Progress do
+            kmax:=CurrentTopDim()-n;
+            k:=1;
+            while k <= kmax and not Progress do
+                if IsGBalanced(L,e,n,k) and HomeomorphicToDisk(L,e,n,k) then
+                    Info(InfoWunden, 1, "ReduceGComplex: replacing the ",k,
+                          "-star of cell ",e," in dimension ",n);
+                    L:=Replace2(L,e,n,k);
+                    Progress:=true;
+                fi;
+                k:=k+1;
+            od;
+            e:=e+1;
+        od;
+        n:=n+1;
+    od;
+until not Progress;
+ 
+return L;
+end);
diff --git a/lib/Wunden/replace.gi b/lib/Wunden/replace.gi
new file mode 100644
index 00000000..1d9a2345
--- /dev/null
+++ b/lib/Wunden/replace.gi
@@ -0,0 +1,365 @@
+#Replace the k-star of the e-th n-dimensional cell in L by a single cell
+InstallGlobalFunction(Replace2, function(L,e,n,k)
+local perm, i, j, d, l, topdim, matched, bndnk1, cnt, orbits, bool, StabU, sign, Elts, g, g1, lg, poslg, c, cell, pos,
+newbnddup, newbnd, bnd, bndcell, h, instar, candidates, cellstar, incell, Boundary, 
+Stab, Dim, action, oldpos, oldindex, RemapBoundary, oldi, signs, found,
+kstar, ReduceTermOld, NormalizeBoundaryOld, matches,
+coefficients, original, expanded, coefficient, rep, twist;
+
+Elts := ShallowCopy(L!.elts);
+kstar:=kStar(L,e,n,k);
+
+pos := function(g)
+local p;
+    p := Position(Elts, g);
+    if p = fail then
+        Add(Elts, g);
+        return Length(Elts);
+    fi;
+    return p;
+end;
+
+oldpos := function(g)
+local p;
+    p := Position(L!.elts, g);
+    if p = fail then
+        Add(L!.elts, g);
+        return Length(L!.elts);
+    fi;
+    return p;
+end;
+
+#ki 09/07
+## ReduceTermOld(degree,j,h): j is an OLD (pre-perm) cell index,
+## h is an index into Elts. Returns [rep, eps] where rep is the smallest
+## index such that Elts[h] and Elts[rep] lie in the same left coset of
+## Stab(degree,j).
+ReduceTermOld := function(degree, j, h)
+local stab, p, t;
+    stab := L!.stabilizer(degree, j);
+    for p in [1..Length(Elts)] do
+        t := Inverse(Elts[h]) * Elts[p];
+        if t in stab then
+            return [ p, L!.action(degree, j, oldpos(t)) ];
+        fi;
+    od;
+    Error("ReduceTermOld: no representative found (should be unreachable)");
+end;
+
+## NormalizeBoundaryOld(bnd,degree): takes a raw list of [±j,h] terms (old
+## indices), groups them by (j, canonical left-coset rep), sums signed
+## multiplicities (folding in the eps twist), drops classes that cancel
+## to zero, and expands surviving classes back into repeated [±j,rep]
+## entries. This is what correctly cancels a coset that shows up twice
+## with opposite signs -- which the old same-signed-value dedup could
+## not do.
+NormalizeBoundaryOld := function(bnd, degree)
+local out, term, j, sgn, h, red, rep, eps, entry, matched, result, cnt;
+    out := [];
+    for term in bnd do
+        j   := AbsoluteValue(term[1]);
+        sgn := SignInt(term[1]);
+        h   := term[2];
+        red := ReduceTermOld(degree, j, h);
+        rep := red[1];
+        eps := red[2];
+        matched := false;
+        for entry in out do
+            if entry[1] = j and entry[2] = rep then
+                entry[3] := entry[3] + sgn*eps;
+                matched := true;
+                break;
+            fi;
+        od;
+        if not matched then
+            Add(out, [ j, rep, sgn*eps ]);
+        fi;
+    od;
+    out := Filtered(out, x -> x[3] <> 0);
+    result := [];
+    for entry in out do
+        for cnt in [1..AbsoluteValue(entry[3])] do
+            Add(result, [ SignInt(entry[3]) * entry[1], entry[2] ]);
+        od;
+    od;
+    return result;
+end;
+#09/07 ging bis hier
+
+perm:=[];  ##perm[i][j] will be the new position of
+           ##the old j-th cell of dimension n+i-1;
+for i in [1..k+1] do
+    if kstar[i]=[] then break;
+    fi;
+    perm[i]:=[];
+    cnt:=0;
+    orbits:=List(kstar[i],x->x[1]);
+    for j in [1..L!.dimension(n+i-1)] do
+        if j in orbits then cnt:=cnt+1;
+        else perm[i][j]:=j-cnt;
+        fi;
+    od;
+od;
+#nrremovedtopcells=cnt
+topdim:=L!.dimension(n+k)-cnt+1;
+
+#ki #???Fehler: d hat keine Funktion
+oldindex := function(p, i, d)
+local q;
+    q := Position(p, i);
+    if q = fail then
+        Error("No old index corresponding to new index ",
+              i, " in dimension ", d);
+    fi;
+    return q;
+end;
+
+signs:=SimplifyLattice(L,kstar,n)[2];
+
+newbnddup:=[];
+#signs:=SimplifyLattice(L,kstar,n)[2];
+for c in [1..Length(kstar[k+1])] do
+    cell:=kstar[k+1][c];
+    bndcell:=L!.boundary(n+k,cell[1]);
+    for bnd in bndcell do
+        h:=cell[2]*Elts[bnd[2]];
+        instar:=false;
+        candidates:=Filtered(kstar[k],y->y[1]=AbsoluteValue(bnd[1]));
+        for cellstar in candidates do
+            if Inverse(h)*cellstar[2] in cellstar[3] then #changed from h*Inverse(cellstar[2]) to the way it is now
+                instar:=true;
+                break;
+            fi;
+        od;
+        if not instar then
+            #sign:=L!.action(n+k,cell[1],pos(cell[2]));
+            #sign:=SignInt(SimplifyLattice(L,kstar,n)[2][c]);
+            Add(newbnddup,[signs[c]*bnd[1],pos(h)]);
+        fi;
+    od;
+od;
+#Print(newbnddup, "\n");
+
+# newbnd:=[];
+# for cell in newbnddup do
+#     bool:=true;
+#     for incell in Filtered(newbnd,y->y[1]=cell[1]) do
+#         if Inverse(Elts[incell[2]])*Elts[cell[2]] in L!.stabilizer(n+k-1,AbsoluteValue(cell[1])) then
+#         # changed from if Elts[incell[2]]*Inverse(Elts[cell[2]]) in L!.stabilizer(n+k-1,AbsoluteValue(cell[1])) then
+#             bool:=false;
+#             break;
+#         fi;
+#     od;
+#     if bool then
+#         Add(newbnd,cell);
+#     fi;
+# od;
+newbnd := NormalizeBoundaryOld(newbnddup,n+k-1); #ki
+
+
+newbnd:=List(newbnd, 
+    cell -> [SignInt(cell[1]) * perm[k][AbsoluteValue(cell[1])],cell[2]]);
+
+#ki
+RemapBoundary := function(b, p, d, i)
+local term, out, newidx;
+    out := [];
+    for term in b do
+        newidx := p[AbsoluteValue(term[1])];
+        if newidx = fail then
+            Error("Boundary remap failed at d=", d, " i=", i,
+                  " old face=", AbsoluteValue(term[1]));
+        fi;
+        Add(out, [ SignInt(term[1]) * newidx, term[2] ]);
+    od;
+    return out;
+end;
+
+
+Boundary := function(d,i)
+local oldi, b;
+
+    if d<n or d>n+k+1 then
+        return L!.boundary(d,i);
+
+    elif d=n then
+        # boundaries in degree n go to degree n-1, which was not changed
+        oldi := oldindex(perm[d-n+1], i, d);
+        return L!.boundary(d, oldi);
+
+    elif d>n and d<n+k then
+        oldi := oldindex(perm[d-n+1], i, d);
+        b := L!.boundary(d, oldi);
+        return RemapBoundary(b, perm[d-n], d, i);
+
+    elif d=n+k then
+        if AbsoluteValue(i)<topdim then
+            oldi := oldindex(perm[d-n+1], AbsoluteValue(i), d);
+            b := L!.boundary(d, oldi);
+            b := RemapBoundary(b, perm[d-n], d, AbsoluteValue(i));
+            if i<0 then
+                b := List(b, t -> [-t[1],t[2]]);
+            fi;
+            return b;
+        elif AbsoluteValue(i)=topdim then
+            if i<0 then
+                return List(newbnd, t -> [-t[1],t[2]]);
+            fi;
+            return newbnd;
+        fi;
+
+    elif d=n+k+1 then
+        if i<0 then
+            return List(bndnk1[-i], t -> [-t[1],t[2]]);
+        fi;
+        return bndnk1[i];
+    fi;
+end;
+
+StabU:=Stabstar(L,e,n,k);
+Stab:=function(d,i)
+    if d<n or d>n+k+1 then
+        return L!.stabilizer(d,i);
+    elif d>=n and d<n+k then
+        return L!.stabilizer(d,Position(perm[d-n+1],AbsoluteValue(i)));
+    elif d=n+k then
+        if AbsoluteValue(i)<topdim then
+            return L!.stabilizer(d,Position(perm[d-n+1],AbsoluteValue(i)));
+        elif AbsoluteValue(i)=topdim then
+            return StabU;
+        fi;
+    elif d=n+k+1 then
+        return L!.stabilizer(d,i);
+    fi;
+end; 
+
+
+
+Dim:=function(d)
+local p;
+    if d<n or d>=n+k+1 then
+        return L!.dimension(d);
+    elif d>=n and d<n+k then
+        p := Compacted(perm[d-n+1]);
+        if Length(p)=0 then
+            return 0;
+        fi;
+        return Maximum(p);
+    elif d=n+k then
+        p := Compacted(perm[d-n+1]);
+        if Length(p)=0 then
+            return 1;
+        fi;
+        return Maximum(p)+1;
+    fi;
+end;
+
+#changes start here
+
+
+
+action := function(d,i,j)
+local g, h, representative, cell, cell2;
+    g := Elts[j];
+
+    if d = n+k and AbsoluteValue(i) = Dim(d) then
+        if not (g in StabU) then
+            # FreeGResolution also asks for signs outside the stabilizer.
+            # Use the same coset decomposition as HAP's ExtendScalars:
+            # g = representative * h, and evaluate the character on h.
+            representative:=CanonicalRightCountableCosetElement(StabU,Inverse(g))^-1;
+            g:=Inverse(representative)*g;
+        fi;
+        # Read the orientation from a constituent top cell. This also works
+        # when the new cell's cellular boundary is zero.
+        cell := kstar[k+1][1];
+        for cell2 in kstar[k+1] do
+            if cell2[1] = cell[1] then
+                h := Inverse(cell2[2])*g*cell[2];
+                if h in cell[3] then
+                    return signs[1]*signs[Position(kstar[k+1],cell2)]*
+                           L!.action(n+k,cell[1],oldpos(h));
+                fi;
+            fi;
+        od;
+        Error("Replace2: star stabilizer does not preserve its top cells");
+
+    elif d >= n and d <= n+k then
+        return L!.action(d, oldindex(perm[d-n+1], AbsoluteValue(i), d), oldpos(g));
+
+    else
+        return L!.action(d, AbsoluteValue(i), oldpos(g));
+    fi;
+end;
+
+
+# An attaching boundary contains the entire oriented chain of each star
+# translate, not one copy of the new cell per constituent face. Collect
+# coefficients modulo the star stabilizer (including its orientation),
+# then divide by the number of constituent top cells. Preserve genuine
+# multiplicities, and verify by expanding back into the old cell module.
+# This must follow the definition of action, used for orientation twists.
+orbits:=List(kstar[k+1],x->x[1]);
+bndnk1:=[];
+for j in [1..L!.dimension(n+k+1)] do
+    bndnk1[j]:=[];
+    coefficients:=[];
+    original:=[];
+    for cell in L!.boundary(n+k+1,j) do
+        if AbsoluteValue(cell[1]) in orbits then
+            Add(original,cell);
+            c:=Position(orbits,AbsoluteValue(cell[1]));
+            lg:=L!.elts[cell[2]]*Inverse(kstar[k+1][c][2]);
+            poslg:=pos(lg);
+            rep:=First([1..poslg],p->Inverse(lg)*Elts[p] in StabU);
+            twist:=action(n+k,topdim,pos(Inverse(lg)*Elts[rep]));
+            if not IsBound(coefficients[rep]) then coefficients[rep]:=0; fi;
+            coefficients[rep]:=coefficients[rep]+SignInt(cell[1])*signs[c]*twist;
+        else
+            Add(bndnk1[j],[SignInt(cell[1])*perm[k+1][AbsoluteValue(cell[1])],cell[2]]);
+        fi;
+    od;
+    expanded:=[];
+    for rep in [1..Length(coefficients)] do
+        if IsBound(coefficients[rep]) and coefficients[rep]<>0 then
+            coefficient:=coefficients[rep]/Length(kstar[k+1]);
+            if not IsInt(coefficient) then
+                Error("Replace2: attaching boundary is not a whole star chain in dimension ",
+                      n+k+1," cell ",j);
+            fi;
+            for cnt in [1..AbsoluteValue(coefficient)] do
+                Add(bndnk1[j],[SignInt(coefficient)*topdim,rep]);
+                for c in [1..Length(kstar[k+1])] do
+                    Add(expanded,[SignInt(coefficient)*signs[c]*kstar[k+1][c][1],
+                                  pos(Elts[rep]*kstar[k+1][c][2])]);
+                od;
+            od;
+        fi;
+    od;
+    if SortedList(NormalizeBoundaryOld(original,n+k)) <>
+       SortedList(NormalizeBoundaryOld(expanded,n+k)) then
+        Error("Replace2: attaching boundary does not expand to the original star chain in dimension ",
+              n+k+1," cell ",j);
+    fi;
+od;
+
+
+return Objectify(HapNonFreeResolution,
+        rec(
+            dimension := Dim,
+            boundary := Boundary,
+            stabilizer := Stab,
+            action := action,
+            elts := Elts,
+            group := L!.group,
+            homotopy := fail,
+            properties := [
+                ["characteristic",0],
+                ["type","resolution"],
+                ["length",1000]
+            ]
+        )
+    );
+end);
+###############################
+
diff --git a/lib/Wunden/restrict_nonfree.gi b/lib/Wunden/restrict_nonfree.gi
new file mode 100644
index 00000000..75c0439e
--- /dev/null
+++ b/lib/Wunden/restrict_nonfree.gi
@@ -0,0 +1,184 @@
+# Restrict a HAP non-free resolution to a finite-index subgroup.
+# After LoadPackage("hap"); S := RestrictNonFreeResolution(R,K);
+
+# Exact membership in K intersect t*H*t^-1, even when H is infinite.
+InstallMethod(\in,
+    [IsMultiplicativeElementWithInverse,
+     IsGroup and HasRestrictionCellStabilizerData],
+    function(g,h)
+        local data;
+        data := RestrictionCellStabilizerData(h);
+        return g in data.subgroup and
+            data.translate^-1*g*data.translate in data.original;
+    end);
+
+InstallGlobalFunction(RestrictNonFreeResolution, function(R,K)
+    local G, one, length, top, T, cosetPosition, q, g, x,
+          data, records, elts, originalIndex, subgroupIndex, prepare,
+          dimension, checkCell, fromOriginal, toOriginal, boundary,
+          stabilizer, action, homotopy, boundaryCache;
+    if not IsHapNonFreeResolution(R) or not IsGroup(K) then
+        Error("Expected a HapNonFreeResolution and a subgroup");
+    fi;
+    G := R!.group; one := One(G); length := Length(R);
+    if one <> One(K) then
+        Error("R!.group and K must use the same element representation; ",
+            "their identity elements differ. Transport K or R first.");
+    fi;
+    if not IsInt(length) or length < 0 then
+        Error("The input must have a finite nonnegative recorded length");
+    fi;
+    # Enumerate LEFT cosets x*K using membership in K. Finite index and
+    # terminating membership tests in K are preconditions.
+    if HasIsFinite(G) and IsFinite(G) then
+        if not IsSubgroup(G,K) then Error("K must be a subgroup of R!.group"); fi;
+    fi;
+    T := [one];
+    cosetPosition := g -> PositionProperty(T,t -> t^-1*g in K);
+    q := 1;
+    while q <= Length(T) do
+        for g in GeneratorsOfGroup(G) do
+            x := g*T[q];
+            if cosetPosition(x)=fail then Add(T,x); fi;
+        od;
+        q := q+1;
+    od;
+    top := length;
+    while top > 0 and R!.dimension(top)=0 do top := top-1; od;
+    data := List([0..top],n -> []);
+    records := List([0..top],n -> []);
+    boundaryCache := List([0..top],n -> []);
+    elts := [one];
+
+    originalIndex := function(g)
+        local pos;
+        pos := Position(R!.elts,g);
+        if pos=fail then Add(R!.elts,g); pos := Length(R!.elts); fi;
+        return pos;
+    end;
+    subgroupIndex := function(g)
+        local pos;
+        pos := Position(elts,g);
+        if pos=fail then Add(elts,g); pos := Length(elts); fi;
+        return pos;
+    end;
+
+    # Orbits of H on the LEFT cosets G/K. For a base coset x*K,
+    # witness[q]*x*K = T[q]*K; the new cell is x^-1*e_i.
+    prepare := function(n)
+        local i, H, generators, owner, witness, base, orbit, cursor, q,
+              next, s, t, z, h, loops, number, entry;
+        if n<0 or n>top or not IsEmpty(data[n+1]) then return; fi;
+        for i in [1..R!.dimension(n)] do
+            H := R!.stabilizer(n,i); generators := GeneratorsOfGroup(H);
+            owner := []; witness := [];
+            for base in [1..Length(T)] do
+                if not IsBound(owner[base]) then
+                    number := Length(records[n+1])+1; t := T[base]^-1;
+                    owner[base] := number; witness[base] := one;
+                    orbit := [base]; cursor := 1; loops := [];
+                    while cursor <= Length(orbit) do
+                        q := orbit[cursor];
+                        for s in generators do
+                            next := cosetPosition(s*T[q]);
+                            if not IsBound(owner[next]) then
+                                owner[next] := number;
+                                witness[next] := s*witness[q];
+                                Add(orbit,next);
+                            fi;
+                            z := witness[next]^-1*s*witness[q];
+                            if z<>one then AddSet(loops,t*z*t^-1); fi;
+                        od;
+                        cursor := cursor+1;
+                    od;
+                    # Schreier's lemma: these generate K intersect t*H*t^-1.
+                    h := SubgroupNC(K,loops);
+                    SetRestrictionCellStabilizerData(h,
+                        rec(subgroup := K, original := H, translate := t));
+                    if HasIsFinite(H) and not IsFinite(H) then
+                        SetIsFinite(h,false);
+                    fi;
+                    entry := rec(originalCell := i, translate := t, stabilizer := h);
+                    Add(records[n+1],entry);
+                fi;
+            od;
+            data[n+1][i] := rec(owner := owner, witness := witness);
+        od;
+    end;
+    dimension := function(n)
+        if not IsInt(n) or n<0 or n>top then return 0; fi;
+        prepare(n); return Length(records[n+1]);
+    end;
+    checkCell := function(n,k)
+        if not IsInt(k) or k=0 or AbsInt(k)>dimension(n) then
+            Error("Cell index out of range");
+        fi;
+    end;
+
+    # Rewrite y*e_i as sign*h*(t*e_i), where y=h*t*s and s fixes e_i.
+    # Only the stabilizer's orientation character is needed here.
+    fromOriginal := function(n,word)
+        local result, term, i, y, q, info, k, entry, s, h, sign;
+        if IsEmpty(word) then return []; fi;
+        if n<0 or n>top then Error("Degree out of range"); fi;
+        prepare(n); result := [];
+        for term in word do
+            i := AbsInt(term[1]); y := R!.elts[term[2]];
+            info := data[n+1][i]; q := cosetPosition(y^-1);
+            k := info.owner[q]; entry := records[n+1][k];
+            s := info.witness[q]^-1; h := y*s^-1*entry.translate^-1;
+            sign := R!.action(n,i,originalIndex(s));
+            Add(result,[SignInt(term[1])*sign*k,subgroupIndex(h)]);
+        od;
+        return result;
+    end;
+    toOriginal := function(n,word)
+        local result, term, entry;
+        result := [];
+        for term in word do
+            checkCell(n,term[1]); entry := records[n+1][AbsInt(term[1])];
+            Add(result,[SignInt(term[1])*entry.originalCell,
+                originalIndex(elts[term[2]]*entry.translate)]);
+        od;
+        return result;
+    end;
+    boundary := function(n,k)
+        local entry, word;
+        checkCell(n,k);
+        if n=0 then return []; fi;
+        if not IsBound(boundaryCache[n+1][AbsInt(k)]) then
+            entry := records[n+1][AbsInt(k)];
+            word := List(R!.boundary(n,entry.originalCell),term ->
+                [term[1],originalIndex(entry.translate*R!.elts[term[2]])]);
+            boundaryCache[n+1][AbsInt(k)] := fromOriginal(n-1,word);
+        fi;
+        return List(boundaryCache[n+1][AbsInt(k)],
+            term -> [SignInt(k)*term[1],term[2]]);
+    end;
+    stabilizer := function(n,k)
+        checkCell(n,k); return records[n+1][AbsInt(k)].stabilizer;
+    end;
+    action := function(n,k,j)
+        local entry;
+        checkCell(n,k); entry := records[n+1][AbsInt(k)];
+        return R!.action(n,entry.originalCell,
+            originalIndex(entry.translate^-1*elts[j]*entry.translate));
+    end;
+    homotopy := fail;
+    if IsFunction(R!.homotopy) then
+        homotopy := function(n,e)
+            if n=-1 then return fromOriginal(0,R!.homotopy(-1,e)); fi;
+            return fromOriginal(n+1,R!.homotopy(n,toOriginal(n,[e])[1]));
+        end;
+    fi;
+    return Objectify(HapNonFreeResolution,rec(
+        group := K, elts := elts, dimension := dimension, boundary := boundary,
+        stabilizer := stabilizer, action := action, homotopy := homotopy,
+        originalResolution := R, cosetRepresentatives := T,
+        restrictionCells := records, toOriginal := toOriginal,
+        fromOriginal := fromOriginal,
+        properties := [["length",length],
+            ["characteristic",EvaluateProperty(R,"characteristic")],
+            ["reduced",false],["type","resolution"]]
+    ));
+end);
diff --git a/lib/Wunden/simplify.gi b/lib/Wunden/simplify.gi
new file mode 100644
index 00000000..217e0009
--- /dev/null
+++ b/lib/Wunden/simplify.gi
@@ -0,0 +1,123 @@
+InstallGlobalFunction(SimplifyLattice, function(L,kstar,n) #n=dimension of the center of the kstar
+local Changed, b, d, j, m, h, g, g1, g2, dim, face, cell, signs, faceorbit, Stab, topcells, involvedcell1, involvedcell2,
+ind1, ind2, sign1, sign2, incidence1, incidence2, mergecell1, mergecell2, FLgroups, FLcodim1, permFL, FL, cob;
+
+FL:=ForgetG(L,kstar,n);
+permFL:=StructuralCopy(FL);
+FL:=StructuralCopy(FL);
+FLcodim1:=[1..Length(FL[Length(FL)-2])];
+FLgroups:=List([1..Length(FL[Length(FL)-1])], i -> [i]);
+dim:=Length(FL)-1;
+
+repeat
+    Changed:=false;
+    #topcells:=
+    for d in [1 .. Length(FL)-2] do
+        m:=1;
+        # JoinCells can remove an entire dimension during this pass.
+        while IsBound(FL[d]) and m <= Length(FL[d]) do
+            if FL[d][m][1] > 0 and ShouldRemove(FL,d,m) then
+                if d=Length(FL)-2 then #for top dimensional cells, we have to keep track of signs
+                    mergecell1:=Cobound(FL)[d][m][2];
+                    mergecell2:=Cobound(FL)[d][m][3]; #these are the two cells which will be merged
+                    #adjust signs
+                    b:=FLcodim1[m]; #old index of the cell that will be deleted
+                    ind1 := fail;
+                    involvedcell1 := fail;
+                    ind2 := fail;
+                    involvedcell2 := fail; 
+                    for j in FLgroups[mergecell1] do
+                        if b in permFL[dim][AbsoluteValue(j)]{[2..Length(permFL[dim][AbsoluteValue(j)])]} then
+                            #involvedcell1:=permFL[dim][j];
+                            ind1:=AbsoluteValue(j);
+                            sign1:=SignInt(j);
+                            break;
+                        fi;
+                    od;
+                    for j in FLgroups[mergecell2] do
+                        if b in permFL[dim][AbsoluteValue(j)]{[2..Length(permFL[dim][AbsoluteValue(j)])]} then
+                            #involvedcell2:=permFL[dim][j];
+                            ind2:=AbsoluteValue(j);
+                            sign2:=SignInt(j);
+                            break;
+                        fi;
+                    od;
+                    #Print("Merging the cells ", FLgroups[mergecell1], " and ", FLgroups[mergecell2]," along the face ",b,"\n", " the components involved are ",ind1," and ",ind2,"\n");
+                    if ind1 = fail then
+                        Error("Could not find face ", b, " in group ", FLgroups[mergecell1]);
+                    fi;
+                    face:=kstar[dim-1][b];
+                    faceorbit:=face[1];
+                    g:=face[2];
+                    involvedcell1:=kstar[dim][ind1][1];
+                    g1:=kstar[dim][ind1][2];
+                    involvedcell2:=kstar[dim][ind2][1];
+                    g2:=kstar[dim][ind2][2];
+                    for cell in Filtered(L!.boundary(n+dim-1,involvedcell1),p->AbsoluteValue(p[1])=faceorbit) do
+                        h:=L!.elts[cell[2]];
+                        Stab:=kstar[dim-1][b][3];
+                        if Inverse(h)*Inverse(g1)*g in Stab then
+                            incidence1:=SignInt(cell[1]);
+                            break;
+                        fi;
+                    od;
+                    for cell in Filtered(L!.boundary(n+dim-1,involvedcell2),p->AbsoluteValue(p[1])=faceorbit) do
+                        h:=L!.elts[cell[2]];
+                        Stab:=kstar[dim-1][b][3];
+                        if Inverse(h)*Inverse(g2)*g in Stab then
+                            incidence2:=SignInt(cell[1]);
+                            break;
+                        fi;
+                    od;
+                    if sign1*incidence1*sign2*incidence2 = 1 then
+                        FLgroups[mergecell2]:=List(FLgroups[mergecell2], y -> -y);
+                    fi;
+                    Append(FLgroups[mergecell1],FLgroups[mergecell2]);
+                    Remove(FLgroups,mergecell2);
+                    Remove(FLcodim1,m);
+                fi;
+
+                # NEW: if this JoinCells call will delete a cell from the
+                # codim-1 level (FL[Length(FL)-2]) via its d2-side removal,
+                # keep FLcodim1 in sync.
+                if d = Length(FL)-3 then
+                    cob := Cobound(FL)[d][m];
+                    FL:=JoinCells(FL,d,m);
+                    Remove(FLcodim1, cob[3]);
+                else
+                    FL:=JoinCells(FL,d,m);
+                fi;
+                
+                # FL:=JoinCells(FL,d,m);
+                Changed:=true;               
+                m:=1;
+            else
+                m:=m+1;
+            fi;
+
+        od;
+    od;
+until not Changed;
+
+signs := List([1..Length(FLgroups[1])], x -> 0);
+for j in FLgroups[1] do
+    signs[AbsoluteValue(j)] := SignInt(j);
+od;
+return [FL,signs];
+end);
+
+InstallGlobalFunction(HomeomorphicToDisk, function(L,e,n,k)
+local kstar, B, simplified;
+
+B:=false;
+kstar:=kStar(L,e,n,k);
+# A zero cellular boundary can leave no incident cells at this level.
+if kstar[k+1]=[] then return false; fi;
+simplified:=SimplifyLattice(L,kstar,n)[1];
+if Length(simplified)=2 and Length(simplified[1])=1 then
+    B:=true;
+fi;
+return B;
+end);
+############################################
+
diff --git a/lib/Wunden/star.gi b/lib/Wunden/star.gi
new file mode 100644
index 00000000..54861d61
--- /dev/null
+++ b/lib/Wunden/star.gi
@@ -0,0 +1,73 @@
+#Generates the k-star of the e-th cell in dimension n in the equivariant complex L 
+#Output: The i-th entry of the output is a list of cells in dimension k+i of the star. The entries have the form [cell orbit,]
+InstallGlobalFunction(kStar, function(L,e,n,k)
+local star, stardup, dim, f, bnd, cell, Li, g, g1, h, id, cell2, Stabf, Stabcell, K1, K2, l,T,Stabe;
+
+stardup := List([1..k+1], i -> []);
+Stabe:=L!.stabilizer(n,e);
+id := Identity(L!.elts[1]);
+stardup[1]:=[[e,id,Stabe]];
+dim:=1;
+
+for f in [1..L!.dimension(n+dim)] do
+    bnd := List(L!.boundary(n+dim, f),x-> [AbsoluteValue(x[1]),x[2]]);
+    for cell in bnd do
+        Li:= Filtered(stardup[dim],x -> x[1]=cell[1]);
+        g:=L!.elts[cell[2]];
+        for cell2 in Li do
+            Stabf:=L!.stabilizer(n+dim,f);
+            Stabcell:=L!.stabilizer(n+dim-1,cell[1]);
+            K1:=Intersection(Stabf,Stabcell);
+            K2:=RightTransversal(Stabcell,K1);
+            for l in K2 do
+                g1:= Inverse(l)*Inverse(g);
+                Add(stardup[dim+1],[f,g1,Stabf]);
+            od;
+        od;
+    od;
+od;
+
+for dim in [2..k] do
+    for f in [1..L!.dimension(n+dim)] do
+        bnd := List(L!.boundary(n+dim, f),x-> [AbsoluteValue(x[1]),x[2]]);
+        for cell in bnd do
+            Li:= Filtered(stardup[dim],x -> x[1]=cell[1]);
+            g:=L!.elts[cell[2]];
+            for cell2 in Li do
+                h:=cell2[2];
+                Stabf:=L!.stabilizer(n+dim,f);
+                Stabcell:=L!.stabilizer(n+dim-1,cell[1]);
+                K1:=Intersection(Stabf,Stabcell);
+                K2:=RightTransversal(Stabcell,K1);
+                for l in K2 do
+                    g1:= h*Inverse(l)*Inverse(g);
+                    Add(stardup[dim+1],[f,g1,Stabf]);
+                od;
+            od;
+        od;
+    od;
+od;
+
+star := List([1..k+1], i -> []);
+star[1]:=[[e,id,Stabe]];
+for dim in [2..k+1] do
+    for cell in stardup[dim] do
+        T:=true;
+        for cell2 in Filtered(star[dim],x->x[1]=cell[1]) do
+            Stabcell:=L!.stabilizer(n+dim-1,cell[1]);
+            if Inverse(cell[2])*cell2[2] in Stabcell then
+                T:=false;
+                break;
+            fi;
+        od;
+        if T then
+            Add(star[dim],cell);
+        fi;
+    od;
+od;
+
+return star;
+end);
+
+#####################
+

From 528679e1177bdbb5f15c2dc5e3f31f75b96679e4 Mon Sep 17 00:00:00 2001
From: Frederik Wunden <frederik@familie-wunden.de>
Date: Thu, 5 Nov 2026 09:11:31 +0000
Subject: [PATCH 2/2] Load Wunden functions when HAP starts

---
 read.g | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/read.g b/read.g
index 5947d819..3474cd4d 100644
--- a/read.g
+++ b/read.g
@@ -578,3 +578,6 @@ SetInfoLevel(InfoWarning,1); #This is GAP's default level
 
 ##AFTER THOUGHT######################################################
 ReadPackageHap("lib/hap_afterthought.gd");
+
+# Equivariant G-complex simplification.
+ReadPackageHap("lib/Wunden/initWunden.gi");