gh-36884: Method braid of class Link loops when the Link contains loops
    
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes #1234" use "Introduce new method to
calculate 1+1"
-->
<!-- Describe your changes here in detail -->

<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes #12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->

The following link `L` is ambient isotopic to the treefoil knot but has
one additional loop.

![grafik](https://github.com/sagemath/sage/assets/47305845/499d14f9-
fe72-406a-8a23-cd94bd0000f1)

```
sage: L = Link([[1, 7, 2, 6], [3, 1, 4, 8], [5, 5, 6, 4], [7, 3, 8, 2]])
sage: L.regions()`
[[-5], [8, 4, 6, 2], [7, -2], [3, -8], [1, -6, 5, -4], [-1, -3, -7]]
sage: L.seifert_circles()
[[5], [1, 7, 3], [2, 8, 4, 6]]
```

Invoking `L.braid()` does not terminate. The problem already occurs on
the first application of a Vogel move which leads to an extended
`pd_code = [[10, 7, 2, 6], [3, 1, 4, 8], [11, 5, 6, 4], [7, 3, 8, 2],
[12, 9, 5, 1], [11, 9, 12, 10]]`. A diagram of the corresponding link
looks like this:

![grafik](https://github.com/sagemath/sage/assets/47305845/5508593b-
b095-4f87-9019-7cfbfe6b4453)


The labels `C1`, `C2 = [3, 1, 4, 8]`, `D = [11, 9, 12, 10]` and `E =
[12, 9, 5, 1]` correspond to the variables in the code and the numbers
to the edges. In case of `C1`the correct list of edges would be `[5, 11,
6, 4]` but the code produces `[11, 5, 6, 4]`. This is caused by usage of
the `index` method of the `list` class which always returns the first
occurrence of a value in the list, i.e. by the line:

```
C2[C2.index(b)] = newedge + 2
```

To fix it I replace the method by a local function:

```
        def idx(cross, edge):
            r"""
            Return the index of an edge in a crossing taking loops into
account.
            A loop appears as an edge which occurs twice in the
crossing.
            In all cases the second occurrence is the correct one needed
in
            the Vogel algorithm.
            """
            i = cross.index(edge)
            if cross.count(edge) > 1:
                return cross.index(edge, i+1)
            else:
                return i
```


> In all cases the second occurrence is the correct one needed in the
Vogel algorithm.

The according argumentation goes as follows:

If the loop is clockwise oriented the edge has a positive sign in the
region, else a negative sign. In the first case the outgoing side of the
edge always occurs after the incoming one in the PD-code (edge 1 in the
first picture below and edge 2 in the second). In the second case the
incoming side of the edge always occurs after the outgoing one in the
PD-code (edge 2 in the first picture below and edge 1 in the second).

![grafik](https://github.com/sagemath/sage/assets/47305845/224324c4-
d392-4944-8406-586838b164ea) ![grafik](https://github.com/sagemath/sage/
assets/47305845/21548500-777f-4a79-99d7-0e59c3bfd8a1)

Thus, this corresponds to the fact that in the first case `a` and `b`
leave the crossings `C1` and `C2` (in the original link) whereas they
arrive there in the second case.

Another way to fix the issue would be to remove all loops from the
diagram before starting the Vogel algorithm. But this would harm regular
isotopy which might violate the expectation of the user. Thus, I
implement this as an optional feature. This gives the user the
possibility to obtain a braid of smaller index if he is fine with
ambient isotopy.


### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [ ] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: #36884
Reported by: Sebastian Oehms
Reviewer(s): Sebastian Oehms, Travis Scrimshaw

Author: Release Manager

src/sage/knots/link.py

@@ -610,10 +610,17 @@ def __ne__(self, other):
         """
         return not self.__eq__(other)
 
-    def braid(self):
+    def braid(self, remove_loops=False):
         r"""
         Return a braid representation of ``self``.
 
+        INPUT:
+
+        - ``remove_loops`` -- boolean (default: ``False``). If set to ``True``
+          loops will be removed first. This can reduce the number of strands
+          needed for an ambient isotopic braid closure. However, this can lead
+          to a loss of the regular isotopy.
+
         OUTPUT: an element in the braid group
 
         .. WARNING::
@@ -634,6 +641,14 @@ def braid(self):
             sage: L.braid()
             (s0*s1^-1)^2*s1^-1
 
+        using ``remove_loops=True``::
+
+            sage: L = Link([[2, 7, 1, 1], [7, 3, 9, 2], [4, 11, 3, 9], [11, 5, 5, 4]])
+            sage: L.braid()
+            s0*s1^-1*s2*s3^-1
+            sage: L.braid(remove_loops=True)
+            1
+
         TESTS::
 
             sage: L = Link([])
@@ -648,7 +663,20 @@ def braid(self):
             sage: A = Link([[[1, 2, -2, -1, -3, -4, 4, 3]], [1, 1, 1, 1]])
             sage: A.braid()
             s0*s1*s2*s3
+
+        Check that :issue:`36884` is solved::
+
+            sage: L = Link([[1, 7, 2, 6], [3, 1, 4, 8], [5, 5, 6, 4], [7, 3, 8, 2]])
+            sage: L.braid()
+            s0^3*s1*s0*s1^-1
+            sage: L.braid(remove_loops=True)
+            s^3
         """
+        if remove_loops:
+            L = self.remove_loops()
+            if L != self:
+                return L.braid(remove_loops=remove_loops)
+
         if self._braid is not None:
             return self._braid
 
@@ -657,8 +685,8 @@ def braid(self):
         if len(comp) > 1:
             L1 = Link(comp[0])
             L2 = Link(flatten(comp[1:], max_level=1))
-            b1 = L1.braid()
-            b2 = L2.braid()
+            b1 = L1.braid(remove_loops=remove_loops)
+            b2 = L2.braid(remove_loops=remove_loops)
             n1 = b1.parent().strands()
             n2 = b2.parent().strands()
             t1 = list(b1.Tietze())
@@ -667,13 +695,26 @@ def braid(self):
             self._braid = B(t1 + t2)
             return self._braid
 
-        # look for possible Vogel moves, perform them and call recursively to the modified link
         pd_code = self.pd_code()
         if not pd_code:
             B = BraidGroup(2)
             self._braid = B.one()
             return self._braid
 
+        # look for possible Vogel moves, perform them and call recursively to the modified link
+        def idx(cross, edge):
+            r"""
+            Return the index of an edge in a crossing taking loops into account.
+            A loop appears as an edge which occurs twice in the crossing.
+            In all cases the second occurrence is the correct one needed in
+            the Vogel algorithm (see :issue:`36884`).
+            """
+            i = cross.index(edge)
+            if cross.count(edge) > 1:
+                return cross.index(edge, i+1)
+            else:
+                return i
+
         seifert_circles = self.seifert_circles()
         newedge = max(flatten(pd_code)) + 1
         for region in self.regions():
@@ -702,12 +743,12 @@ def braid(self):
                             # C1   C2   existing crossings
                             # -------------------------------------------------
                             C1 = newPD[newPD.index(heads[a])]
-                            C1[C1.index(a)] = newedge + 1
+                            C1[idx(C1, a)] = newedge + 1
                             C2 = newPD[newPD.index(tails[b])]
-                            C2[C2.index(b)] = newedge + 2
+                            C2[idx(C2, b)] = newedge + 2
                             newPD.append([newedge + 3, newedge, b, a]) # D
                             newPD.append([newedge + 2, newedge, newedge + 3, newedge + 1]) # E
-                            self._braid = Link(newPD).braid()
+                            self._braid = Link(newPD).braid(remove_loops=remove_loops)
                             return self._braid
                         else:
                             # -------------------------------------------------
@@ -723,12 +764,12 @@ def braid(self):
                             #  /   \
                             # -------------------------------------------------
                             C1 = newPD[newPD.index(heads[-a])]
-                            C1[C1.index(-a)] = newedge + 1
+                            C1[idx(C1, -a)] = newedge + 1
                             C2 = newPD[newPD.index(tails[-b])]
-                            C2[C2.index(-b)] = newedge + 2
+                            C2[idx(C2, -b)] = newedge + 2
                             newPD.append([newedge + 2, newedge + 1, newedge + 3, newedge]) # D
                             newPD.append([newedge + 3, -a, -b, newedge]) # E
-                            self._braid = Link(newPD).braid()
+                            self._braid = Link(newPD).braid(remove_loops=remove_loops)
                             return self._braid
 
         # We are in the case where no Vogel moves are necessary.
@@ -2362,6 +2403,50 @@ def regions(self):
             regions.append(region)
         return regions
 
+    def remove_loops(self):
+        r"""
+        Return an ambient isotopic link in which all loops are removed.
+
+        EXAMPLES::
+
+            sage: b = BraidGroup(4)((3, 2, -1, -1))
+          sage: L = Link(b)
+          sage: L.remove_loops()
+          Link with 2 components represented by 2 crossings
+          sage: K4 = Link([[1, 7, 2, 6], [3, 1, 4, 8], [5, 5, 6, 4], [7, 3, 8, 2]])
+          sage: K3 = K4.remove_loops()
+          sage: K3.pd_code()
+          [[1, 7, 2, 4], [3, 1, 4, 8], [7, 3, 8, 2]]
+          sage: U = Link([[1, 2, 2, 1]])
+          sage: U.remove_loops()
+          Link with 1 component represented by 0 crossings
+        """
+        pd = self.pd_code()
+        new_pd = []
+        loop_crossings = []
+        for cr in pd:
+            if len(set(cr)) == 4:
+                new_pd.append(list(cr))
+            else:
+                loop_crossings.append(cr)
+        if not loop_crossings:
+            return self
+        if not new_pd:
+            # trivial knot
+            return type(self)([])
+        new_edges = flatten(new_pd)
+        for cr in loop_crossings:
+            rem = set([e for e in cr if e in new_edges])
+            if len(rem) == 2:
+                # put remaining edges together
+                a, b = sorted(rem)
+                for ncr in new_pd:
+                    if b in ncr:
+                        ncr[ncr.index(b)] = a
+                        break
+        res = type(self)(new_pd)
+        return res.remove_loops()
+
     @cached_method
     def mirror_image(self):
         r"""