An involution is a permutation that equals its inverse; thought of another way, it's a matching on an n-vertex complete graph, where a position in the permutation is mapped to its partner in the matching (or itself if it is unmatched). Here are the 26 involutions on 5 elements (26 is a telephone number, the numbers that count involutions):

In this sequence, each two consecutive involutions differ from each other in a simple way: two consecutive elements mapped to themselves may become mapped to each other or vice versa, one of the elements of a pair of items mapped to each other may shift position by one, or two consecutive elements that belong to two different pairs may swap partners. These types of changes generate an $$(n-1)$$-regular multigraph on the set of involutions, but it's not a simple graph because sometimes it's possible to get the same change in two different ways.

The overall structure of the sequence is that it starts with all of the involutions that keep the last item fixed (recursively generated as the involutions on one fewer item), and then as in the Steinhausâ€“Johnsonâ€“Trotter algorithm it sweeps the partner of the last item back and forth over the remaining items, at the end of each sweep performing a step in the recursive sequence for two fewer items. Because of the low ratio of recursive steps to sweep steps, the time to generate the sequence is constant per involution. This algorithm demonstrates the existence of a Hamiltonian cycle in the multigraph of involutions and their changes.

My implementation of all this is now towards the end of PADS/Permutations.py. I'm very glad that it was easy to do some thorough unit testing for this one; without that, I would have had a difficult time convincing myself both that the algorithm is correct and that my implementation is correct, because the details are a bit messy.