From charlesreid1

Permutation Counts

For counting the number of permutations of a Rubik's Cube, see Rubiks Cube/Numbers

Permutation Algebra

In order to create a system for talking about and dealing with permutations, we follow Volume 3 of Knuth's The Art of Computer Programming (see AOCP).

Representation of a Rubik's Cube

The mathematical. representation of a Rubik's Cube is discussed on the Rubiks Cube/Tuple page, but in short, that page concludes that we should represent the Rubik's Cube by using a particular permutation of an n-tuple, where n is the number of faces.

A 3x3 Rubik's Cube state is represented by a particular permutation of the integers 1 to 36. The solved state is (1 2 3 4 ... 35 36).

A 4x4 Rubik's Cube state is represented by a particular permutation of the integers 1 to 96. The solved state is (1 2 3 4 ... 95 96).

Writing a Permutation

Knuth introduces the top-bottom notation for permutations, in which a particular permutation of a sequence of n integers is written by first writing each element of the sequence in increasing order on the top row, then writing the occurrence of each element in the order it occurs on the bottom row:


a = \bigl(\begin{smallmatrix}
  1 & 2 & 3 & \cdots & n-1 & n \\
  2 & 3 & 4 & \cdots &  n  & 1
\end{smallmatrix}\bigr)

Intercalation Product

We now define an intercalation operation with two permutations; this is basically an operation that interleaves two permutations. We will see why this is useful in a moment.

Suppose we have two permutations, \alpha and \beta, of four objects \{a, b, c, d\}, each occurring multiple times:


\alpha = \bigl(\begin{smallmatrix}
  a & a & b & c & d \\
  c & a & d & a & b 
\end{smallmatrix}\bigr)


\beta = \bigl(\begin{smallmatrix}
  a & b & d & d & d \\
  b & d & d & a & d 
\end{smallmatrix}\bigr)

Then we can define the intercalation product \alpha \top \beta as the elements of these permutations combined in a way that interleaves elements of both, in a way that groups all elements by the letter on the top row, but sorts within those letters according to the original order in \alpha and \beta. For our example:


\alpha \top \beta = \bigl(\begin{smallmatrix}
  a & a & b & c & d \\
  c & a & d & a & b 
\end{smallmatrix}\bigr) \top \bigl(\begin{smallmatrix}
  a & b & d & d & d \\
  b & d & d & a & d 
\end{smallmatrix}\bigr) = 
\bigl(\begin{smallmatrix}
  a & a & a & b & b & c & d & d & d & d \\
  c & a & b & d & d & a & b & d & a & d
\end{smallmatrix}\bigr)

This is basically an interleaving operation. All top-bottom pairs with a at the top are grouped together; all from \alpha come first, and all from \beta come second.

The first two a items in \alpha \top \beta come from \alpha, the third a item comes from \beta.

Etc.

Properties of Intercalation

Before using the intercalation product, let's define a few properties :

  • If \alpha \top \pi = \beta \top \pi or \pi \top \alpha = \pi \top \beta this implies \alpha = \beta
  • Identity element exists such that \epsilon \top \alpha = \alpha \top \epsilon = \alpha
  • Commutative property only holds if \alpha independent of \beta; then \alpha \top \beta = \beta \top \alpha. This property does not hold in general.

Cycles and Intercalation

Cycles are elements that are swapped in some prescribed way. For example, suppose we have a sequence,


\bigl(\begin{smallmatrix}
  x_1 & x_2 & \dots & x_n
\end{smallmatrix}\bigr)

and further suppose that we shift this sequence to the left by one element, and write in the two line notation:


\bigl(\begin{smallmatrix}
    x_1 & x_2 & \dots & x_{n-1} & x_n \\
    x_2 & x_3 & \dots & x_{n} & x_1 
\end{smallmatrix}\bigr)

We can observe that the product of disjoint cycles is the same as their intercalation.

Now we have a full "product" permutation.

Factoring Permutations

An Example

Suppose we have a permutation:


\pi = \bigl(\begin{smallmatrix}
    a & a & b & b & b & b & b & c & c & c & d & d & d & d & d \\
    d & b & c & b & c & a & c & d & a & d & d & b & b & b & d
\end{smallmatrix}\bigr)

and suppose that we want to "factor" the permutation - that is, to write the permutation as the intercalation product of independent, disjoint cycles \pi = \alpha \top \beta \top \dots \top \gamma.

Assembling Factors

We can extract each factor one at a time using the following algorithm:

Start by supposing that the first factor \alpha contains the first symbol a. If this assumption is true, then \alpha must map a to the same thing that the final permutation maps a to, namely, the first column of \pi, which is the combination \begin{smallmatrix}
    a \\
    d 
\end{smallmatrix}
. This means that the first entry in \alpha must turn a into d.

Now we suppose that \alpha contains d, as a consequence of the prior step. Let's find the leftmost d in the top line, and see what symbol it will map to. It maps to b. This means that the first entry in \alpha must turn d into d and should therefore contain \begin{smallmatrix}
    d \\
    d
\end{smallmatrix}

Now we start with the outcome of the previous step, which is another d. Since we already used the first d-d pairing, we look for the second leftmost d, which is part of the combination/mapping \begin{smallmatrix}
    d \\
    b
\end{smallmatrix}
. This implies the mapping d-b should be in \alpha, and we use the mapping outcome b as the starting point for the next step.

If we keep doing this, eventually we will get:


\alpha = \bigl(\begin{smallmatrix}
    a & d & d & b & c & d & b & b & c \\
    d & d & b & c & d & b & b & c & a
\end{smallmatrix}\bigr)

Why Does This Work

Let's pause for a moment and see what's happening. What we're doing is following a thread between the top and bottom rows of the permutation; this thread tells us how elements are being moved around to create permutations.

(A simpler but easier way to see this is by comparing two permutations of 1 2 3 4 5 6: consider the permutation 2 1 3 4 6 5, versus the permutation 2 4 5 6 1 3. The first permutation swaps positions 0 and 1, and positions 4 and 5, independently; the second permutation mixes all positions together.)

We are assembling \alpha piece by piece, by pulling out pairs from the top and bottom row of \pi and putting them into \alpha. At some point we will come back to the starting point, the symbol a, and we will be finished finding the first factor \alpha, which is a disjoint cycle.

By starting from the top row and following where it leads in the bottom row, and continuing until we return to the original starting element in the top row, we can carve up the permutation into groups of pieces exchanged with one another and not with any other pieces, or groups of pieces that don't move.

Forming All Factors

We then continue the process of assembling factors from what is left of \pi. (Note that if \pi is prime, every element of \pi will appear in \alpha and there will be no further products.) Eventually we will have a number of factors,


\pi = \alpha \top \beta \top \dots \top \gamma

The result of applying this procedure, of skipping between the top and bottom rows, is a set of 1 or more permutation factors. If we continue with the example \pi above, we'll eventually get:

Original:


\pi = \bigl( \begin{smallmatrix}
    a & a & b & b & b & b & b & c & c & c & d & d & d & d & d \\
    d & b & c & b & c & a & c & d & a & d & d & b & b & b & d
\end{smallmatrix}\bigr)

First factor pulled out:


\pi = \bigl( \begin{smallmatrix}
    a & d & d & b & c & d & b & b & c \\
    d & d & b & c & d & b & b & c & a
\end{smallmatrix} \bigr)
\top
\bigl( \begin{smallmatrix}
    a & b & b & c & d & d \\
    b & a & c & d & b & d
\end{smallmatrix} \bigr)

Second factor pulled out:


\pi = \bigl( \begin{smallmatrix}
    a & d & d & b & c & d & b & b & c \\
    d & d & b & c & d & b & b & c & a
\end{smallmatrix} \bigr)
\top
\bigl( \begin{smallmatrix}
    a & b \\
    b & a 
\end{smallmatrix} \bigr)
\top 
\bigl( \begin{smallmatrix}
    b & c & d & d \\
    c & d & b & d
\end{smallmatrix} \bigr)

and the last factor is a piece that is not moved in this permutation:


\pi = \bigl( \begin{smallmatrix}
    a & d & d & b & c & d & b & b & c \\
    d & d & b & c & d & b & b & c & a
\end{smallmatrix} \bigr)
\top
\bigl( \begin{smallmatrix}
    a & b \\
    b & a 
\end{smallmatrix} \bigr)
\top 
\bigl( \begin{smallmatrix}
    b & c & d  \\
    c & d & b
\end{smallmatrix} \bigr)
\top 
\bigl( \begin{smallmatrix}
    d \\
    d
\end{smallmatrix} \bigr)

Thus the permutation \pi can be expressed as the intercalation of four independent cycles.

To relate this back to the Rubik's Cube, we can start with a sequence of interest, like U R D D B, and write the tuple representing the outcome of this sequence when it is applied to the cube. In this way we represent a move sequence as a tuple or as a permutation.

Next, we factor this permutation the way we factored \pi, into the intercalation product of independent cycles. These are groups of pieces being swapped each time the cycle is applied.

Now if one group of pieces being swapped each time the cycle is applied is of size 4 and another group is of size 3 and yet another is of size 20, then we know that the overall sequence order (the number of times the sequence should be applied to a solved cube to return it back to the solved state) is LCM(3,4,20) = 60.

This procedure illustrates what Knuth refers to as Theorem A.

(Note: had we initially assumed \alpha contained b instead of a, we would end up starting by pulling out a different factor, but we would ultimately end up with the same set of four factors.)

Application to Rubiks Cube

Sequence R

Let's look at an example of applying this to a 4x4 Rubiks Revenge cube. We will consider the sequence "R" (turning the right face clockwise once).

Start with the solved state:

(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)

Now, when we apply the move R, we get:

(1 2 3 36 5 6 7 40 9 10 11 44 13 14 15 48 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 84 37 38 39 88 41 42 43 92 45 46 47 96 61 57 53 49 62 58 54 50 63 59 55 51 64 60 56 52 16 66 67 68 12 70 71 72 8 74 75 76 4 78 79 80 81 82 83 77 85 86 87 73 89 90 91 69 93 94 95 65)

When we carry out the Agorithm A procedure and factor this permutation, we get many factors that are independent (meaning they do not move when we carry out the sequence). We have six groups of four faces that cycle through one another (one corner piece with three faces, two double edges with two faces, and one square piece with one face, for a total of six faces):

Factor sizes: {1, 4}
Factors:
[36, 84, 77, 4]
[40, 88, 73, 8]
[44, 92, 69, 12]
[48, 96, 65, 16]
[61, 64, 52, 49]
[57, 63, 56, 50]
[53, 62, 60, 51]
[58, 59, 55, 54]
Independent Pieces: [1, 2, 3, 5, 6, 7, 9, 10, 11, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 41, 42, 43, 45, 46, 47, 66, 67, 68, 70, 71, 72, 74, 75, 76, 78, 79, 80, 81, 82, 83, 85, 86, 87, 89, 90, 91, 93, 94, 95]
Least common multiple: 4

The largest set of faces that are exchanged is 4, and the smallest is 1. No other groups of faces being exchanged have any other sizes. This means that if we apply the sequence 4 times, each of those groups of faces being interchanged will have returned to their original state.

This tells us what we already knew: that if we apply the sequence "R", it rotates groups of pieces in a sequence of 4 moves each, so overall the order of this permutation is 4 - if we apply the sequence R to a solved 4x4 Rubik's Revenge cube 4 times, the cube will return to the solved state.

To formalize this, if we have cycles with arbitrary lengths, we must apply the sequence a number of times equal to the least common multiple of each factor's size. (For example, if we had a cycle of length 3 above, the cycle order would have been 12 - because the sequence must be applied 12 times before the 4-cycle face exchanges "sync up" with the 3-cycle face exchanges.)

Sequence U R U' R'

Let's look at a slightly longer sequence: U R U' R'.

We start with the n-tuple representing the solved state (top row):

(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)

and the n-tuple representing this permutation (bottom row):

(1 2 3 77 5 6 7 73 9 10 11 69 16 12 8 20 17 18 19 36 21 22 23 24 25 26 27 28 29 30 31 32 49 50 51 33 37 38 39 40 41 42 43 44 45 46 47 48 13 56 60 64 53 54 55 34 57 58 59 35 61 62 63 4 96 66 67 68 14 70 71 72 15 74 75 76 65 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 52)

Now, if we factor out the n-tuple, we get some factors (independent cycles) of length 3 and length 6:

Factor sizes: {1, 3, 6}
Factors:
[77, 65, 96, 52, 64, 4]
[73, 15, 8]
[69, 14, 12]
[16, 20, 36, 33, 49, 13]
[50, 56, 34]
[51, 60, 35]
Independent Pieces: [1, 2, 3, 5, 6, 7, 9, 10, 11, 17, 18, 19, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 53, 54, 55, 57, 58, 59, 61, 62, 63, 66, 67, 68, 70, 71, 72, 74, 75, 76, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95]
Least common multiple: 6

The least common multiple of 1, 3, and 6 is 6, meaning the order of the sequence U R U' R' is 6. If we apply the sequence U R U' R' six times to a solved cube, it will return to the solved state.

Sequence U R

Now we'll examine the longer sequence U R, which has an order of 105.

We start with the n-tuple representing the solved state (top row):

(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)

and the n-tuple representing this permutation (bottom row):

(13 9 5 1 14 10 6 2 15 11 7 3 48 44 40 36 33 34 35 84 21 22 23 24 25 26 27 28 29 30 31 32 61 57 53 49 37 38 39 88 41 42 43 92 45 46 47 96 16 66 67 68 62 58 54 50 63 59 55 51 64 60 56 52 17 18 19 20 12 70 71 72 8 74 75 76 4 78 79 80 81 82 83 77 85 86 87 73 89 90 91 69 93 94 95 65)

and this factors into:

Factor sizes: {1, 3, 4, 7, 15}
Factors:
[13, 48, 96, 65, 17, 33, 61, 64, 52, 68, 20, 84, 77, 4, 1]
[9, 15, 40, 88, 73, 8, 2]
[5, 14, 44, 92, 69, 12, 3]
[10, 11, 7, 6]
[36, 49, 16]
[34, 57, 63, 56, 50, 66, 18]
[35, 53, 62, 60, 51, 67, 19]
[58, 59, 55, 54]
Independent Pieces: [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 37, 38, 39, 41, 42, 43, 45, 46, 47, 70, 71, 72, 74, 75, 76, 78, 79, 80, 81, 82, 83, 85, 86, 87, 89, 90, 91, 93, 94, 95]
Least common multiple: 420

However, the adventurous cuber will find, when actually carrying out this move sequence, that the order is in fact 105, and not 420.

So what gives? Why is the predicted order 4 times larger than the one we actually see?

This move sequence illustrates a special case that we did not see earlier: by indexing each of the four center faces with different integers, we have made the implicit assumption that the square faces are not interchangeable. This was necessary for corner and double-edge pieces, because orientation was important - double edge blue faces are not interchangeable, and the double edge with a blue face linked to a red face is different and distinct from a double edge with a blue face linked to an orange face.

However, square faces are not linked to any other faces, so distinguishing between the NW and SE faces do not matter. But rather than complicate matters by changing our tuple representation, we can simply account for the fact that the squares are interchangeable.

We can factor the n-tuple and get a list of independent cycles. The factors have orders of 1, 3, 4, 7, and 15. So nominally, the sequence U R would have an order of LCM(1,3,4,7,15) = 420. However, empirically, we see that the square faces can be interchanged, so the cycle of four square faces does not contribute to the overall order of the sequence.

Thus, the order of the sequence U R is actually LCM(1,3,7,15) = 105.

We must look for cycles that contain only the squares, and eliminate them from the count, when computing the order of a sequence.

Factors:
[13, 48, 96, 65, 17, 33, 61, 64, 52, 68, 20, 84, 77, 4, 1]
[9, 15, 40, 88, 73, 8, 2]
[5, 14, 44, 92, 69, 12, 3]
[10, 11, 7, 6]
[36, 49, 16]
[34, 57, 63, 56, 50, 66, 18]
[35, 53, 62, 60, 51, 67, 19]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[37]
[38]
[39]
[41]
[42]
[43]
[45]
[46]
[47]
[58, 59, 55, 54]
[70]
[71]
[72]
[74]
[75]
[76]
[78]
[79]
[80]
[81]
[82]
[83]
[85]
[86]
[87]
[89]
[90]
[91]
[93]
[94]
[95]

Code

To implement a program to turn a sequence of moves into a permutation, and factor the permutation into its component cycles, we utilized the NxNxN rubiks solver here: https://github.com/dwalton76/rubiks-cube-NxNxN-solver

We also made modifications to implement some functionality not contained in the original cube object. Forked version is here: https://charlesreid1.com:3000/charlesreid1/rubiks-cube-nnn-solver

We have a main driver program to perform the factoring here: https://charlesreid1.com:3000/charlesreid1/rubiks-cube-cycles

Let's run through the essential components of how we do this in code.

Basic Cube Operations

Some of the basic cube operations that we need to be able to do:

  • Turn a representation of the cube into integers
  • Get the index of the center squares on each face
  • Obtain the index shifts that occur when a particular move sequence is applied (e.g., the rotation U will move face 13 to face 1, face 1 to face 4, etc.

The latter part was the extension that I needed to make to the NxNxN rubiks solver.

Face Maps for Moves

To construct a face map (a mapping of which "before" face index would become which "after" face index), we did not need to implement any new functionality - we simply based it on the existing method to apply a rotation, which already computes which indices to swap. But instead of swapping faces, we simply save the index pairs and return them.

Unfortunately, the original code is a bit messy and inelegant, so this method is too, but it was much less complicated than the alternative.

    def rotation_map(self, action):
        """
        This returns a set of tuples
        that correspond to the movements
        of each piece on the cube
        for the move "action".
        """
        results = []

        if action[-1] in ("'", "`"):
            reverse = True
            action = action[0:-1]
        else:
            reverse = False

        # 2Uw, Uw and 2U all mean rotate the top 2 U rows
        # 3Uw and 3U mean rotate the top 3 U rows
        if len(action) >= 2 and action[0].isdigit() and action[1].isdigit():
            rows_to_rotate = int(action[0:2])
            action = action[2:]
        elif action[0].isdigit():
            rows_to_rotate = int(action[0])
            action = action[1:]
        else:
            # Uw also means rotate the top 2 U rows
            if 'w' in action:
                rows_to_rotate = 2
            else:
                rows_to_rotate = 1

        # We've accounted for this so remove it
        if 'w' in action:
            action = action.replace('w', '')

        # The digit at the end indicates how many quarter turns to do
        if action[-1].isdigit():
            quarter_turns = int(action[-1])
            action = action[0:-1]
        else:
            quarter_turns = 1

        side_name = action

        # Skip moves x, y, z

        side = self.sides[side_name]
        min_pos = side.min_pos
        max_pos = side.max_pos

        # Rotation of the squares of the face itself
        for turn in range(quarter_turns):

            # Number the squares of the faces
            # in a 2D array, then re-use all the
            # existing transforms
            oldface = []
            counter = min_pos
            for ii in range(self.size):
                facerow = []
                for jj in range(self.size):
                    facerow.append(counter)
                    counter += 1
                oldface.append(facerow)

            # This is not what we want, because it returns
            # the face colors, not the face index numbers:
            #oldface = side.get_face_as_2d_list()

            if reverse:
                newface = rotate_counter_clockwise(oldface)
            else:
                newface = rotate_clockwise(oldface)

            oldface = compress_2d_list(oldface)
            newface = compress_2d_list(newface)

            for (j,k) in zip(oldface,newface):
                results.append((j,k))



        # Again skip rotation of entire cube with x, y, z


        if side_name == "U":

            for turn in range(quarter_turns):

                # rotate the connecting row(s) of the surrounding sides
                for row in range(rows_to_rotate):
                    left_first_square = self.squares_per_side + 1 + (row * self.size)
                    left_last_square = left_first_square + self.size - 1

                    front_first_square = (self.squares_per_side * 2) + 1 + (row * self.size)
                    front_last_square = front_first_square + self.size - 1

                    right_first_square = (self.squares_per_side * 3) + 1 + (row * self.size)
                    right_last_square = right_first_square + self.size - 1

                    back_first_square = (self.squares_per_side * 4) + 1 + (row * self.size)
                    back_last_square = back_first_square + self.size - 1

                    if reverse:
                        for square_index in range(left_first_square, left_last_square + 1):
                            old_index = square_index
                            new_index = square_index + (3 * self.squares_per_side)
                            results.append((old_index,new_index))

                        for square_index in range(front_first_square, front_last_square + 1):
                            old_index = square_index
                            new_index = square_index - self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(right_first_square, right_last_square + 1):
                            old_index = square_index
                            new_index = square_index - self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(back_first_square, back_last_square + 1):
                            old_index = square_index
                            new_index = square_index - self.squares_per_side
                            results.append((old_index,new_index))

                    else:
                        for square_index in range(left_first_square, left_last_square + 1):
                            old_index = square_index
                            new_index = square_index + self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(front_first_square, front_last_square + 1):
                            old_index = square_index
                            new_index = square_index + self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(right_first_square, right_last_square + 1):
                            old_index = square_index
                            new_index = square_index + self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(back_first_square, back_last_square + 1):
                            old_index = square_index
                            new_index = square_index - (3 * self.squares_per_side)
                            results.append((old_index,new_index))

        elif side_name == "L":

            for turn in range(quarter_turns):

                # rotate the connecting row(s) of the surrounding sides
                for row in range(rows_to_rotate):

                    top_first_square = 1 + row
                    top_last_square = top_first_square + ((self.size - 1) * self.size)

                    front_first_square = (self.squares_per_side * 2) + 1 + row
                    front_last_square = front_first_square + ((self.size - 1) * self.size)

                    down_first_square = (self.squares_per_side * 5) + 1 + row
                    down_last_square = down_first_square + ((self.size - 1) * self.size)

                    back_first_square = (self.squares_per_side * 4) + self.size - row
                    back_last_square = back_first_square + ((self.size - 1) * self.size)

                    top_squares = []
                    for square_index in range(top_first_square, top_last_square + 1, self.size):
                        top_squares.append(square_index)

                    front_squares = []
                    for square_index in range(front_first_square, front_last_square + 1, self.size):
                        front_squares.append(square_index)

                    down_squares = []
                    for square_index in range(down_first_square, down_last_square + 1, self.size):
                        down_squares.append(square_index)

                    back_squares = []
                    for square_index in range(back_first_square, back_last_square + 1, self.size):
                        back_squares.append(square_index)

                    if reverse:
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = front_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(front_first_square, front_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

                        back_squares = list(reversed(back_squares))
                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = back_squares[index]
                            results.append((old_index,new_index))

                        top_squares = list(reversed(top_squares))
                        for (index, square_index) in enumerate(range(back_first_square, back_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

                    else:
                        back_squares = list(reversed(back_squares))
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = back_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(front_first_square, front_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = front_squares[index]
                            results.append((old_index,new_index))

                        down_squares = list(reversed(down_squares))
                        for (index, square_index) in enumerate(range(back_first_square, back_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

        elif side_name == "F":

            for turn in range(quarter_turns):

                # rotate the connecting row(s) of the surrounding sides
                for row in range(rows_to_rotate):
                    top_first_square = (self.squares_per_side - self.size) + 1 - (row * self.size)
                    top_last_square = top_first_square + self.size - 1

                    left_first_square = self.squares_per_side + self.size - row
                    left_last_square = left_first_square + ((self.size - 1) * self.size)

                    down_first_square = (self.squares_per_side * 5) + 1 + (row * self.size)
                    down_last_square = down_first_square + self.size - 1

                    right_first_square = (self.squares_per_side * 3) + 1 + row
                    right_last_square = right_first_square + ((self.size - 1) * self.size)

                    #log.info("top first %d, last %d" % (top_first_square, top_last_square))
                    #log.info("left first %d, last %d" % (left_first_square, left_last_square))
                    #log.info("down first %d, last %d" % (down_first_square, down_last_square))
                    #log.info("right first %d, last %d" % (right_first_square, right_last_square))

                    top_squares = []
                    for square_index in range(top_first_square, top_last_square + 1):
                        top_squares.append(square_index)

                    left_squares = []
                    for square_index in range(left_first_square, left_last_square + 1, self.size):
                        left_squares.append(square_index)

                    down_squares = []
                    for square_index in range(down_first_square, down_last_square + 1):
                        down_squares.append(square_index)

                    right_squares = []
                    for square_index in range(right_first_square, right_last_square + 1, self.size):
                        right_squares.append(square_index)

                    if reverse:
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1)):
                            old_index = square_index
                            new_index = right_squares[index]
                            results.append((old_index,new_index))

                        top_squares = list(reversed(top_squares))
                        for (index, square_index) in enumerate(range(left_first_square, left_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1)):
                            old_index = square_index
                            new_index = left_squares[index]
                            results.append((old_index,new_index))

                        down_squares = list(reversed(down_squares))
                        for (index, square_index) in enumerate(range(right_first_square, right_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

                    else:
                        left_squares = list(reversed(left_squares))
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1)):
                            old_index = square_index
                            new_index = left_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(left_first_square, left_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

                        right_squares = list(reversed(right_squares))
                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1)):
                            old_index = square_index
                            new_index = right_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(right_first_square, right_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

        elif side_name == "R":

            for turn in range(quarter_turns):

                # rotate the connecting row(s) of the surrounding sides
                for row in range(rows_to_rotate):

                    top_first_square = self.size - row
                    top_last_square = self.squares_per_side

                    front_first_square = (self.squares_per_side * 2) + self.size - row
                    front_last_square = front_first_square + ((self.size - 1) * self.size)

                    down_first_square = (self.squares_per_side * 5) + self.size - row
                    down_last_square = down_first_square + ((self.size - 1) * self.size)

                    back_first_square = (self.squares_per_side * 4) + 1 + row
                    back_last_square = back_first_square + ((self.size - 1) * self.size)

                    #log.info("top first %d, last %d" % (top_first_square, top_last_square))
                    #log.info("front first %d, last %d" % (front_first_square, front_last_square))
                    #log.info("down first %d, last %d" % (down_first_square, down_last_square))
                    #log.info("back first %d, last %d" % (back_first_square, back_last_square))

                    top_squares = []
                    for square_index in range(top_first_square, top_last_square + 1, self.size):
                        top_squares.append(square_index)

                    front_squares = []
                    for square_index in range(front_first_square, front_last_square + 1, self.size):
                        front_squares.append(square_index)

                    down_squares = []
                    for square_index in range(down_first_square, down_last_square + 1, self.size):
                        down_squares.append(square_index)

                    back_squares = []
                    for square_index in range(back_first_square, back_last_square + 1, self.size):
                        back_squares.append(square_index)

                    if reverse:
                        back_squares = list(reversed(back_squares))
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = back_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(front_first_square, front_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = front_squares[index]
                            results.append((old_index,new_index))

                        down_squares = list(reversed(down_squares))
                        for (index, square_index) in enumerate(range(back_first_square, back_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

                    else:
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = front_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(front_first_square, front_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

                        back_squares = list(reversed(back_squares))
                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = back_squares[index]
                            results.append((old_index,new_index))

                        top_squares = list(reversed(top_squares))
                        for (index, square_index) in enumerate(range(back_first_square, back_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

        elif side_name == "B":

            for turn in range(quarter_turns):

                # rotate the connecting row(s) of the surrounding sides
                for row in range(rows_to_rotate):
                    top_first_square = 1 + (row * self.size)
                    top_last_square = top_first_square + self.size - 1

                    left_first_square = self.squares_per_side + 1 + row
                    left_last_square = left_first_square + ((self.size - 1) * self.size)

                    down_first_square = (self.squares_per_side * 6)  - self.size + 1 - (row * self.size)
                    down_last_square = down_first_square + self.size - 1

                    right_first_square = (self.squares_per_side * 3) + self.size - row
                    right_last_square = right_first_square + ((self.size - 1) * self.size)

                    #log.info("top first %d, last %d" % (top_first_square, top_last_square))
                    #log.info("left first %d, last %d" % (left_first_square, left_last_square))
                    #log.info("down first %d, last %d" % (down_first_square, down_last_square))
                    #log.info("right first %d, last %d" % (right_first_square, right_last_square))

                    top_squares = []
                    for square_index in range(top_first_square, top_last_square + 1):
                        top_squares.append(square_index)

                    left_squares = []
                    for square_index in range(left_first_square, left_last_square + 1, self.size):
                        left_squares.append(square_index)

                    down_squares = []
                    for square_index in range(down_first_square, down_last_square + 1):
                        down_squares.append(square_index)

                    right_squares = []
                    for square_index in range(right_first_square, right_last_square + 1, self.size):
                        right_squares.append(square_index)

                    if reverse:
                        left_squares = list(reversed(left_squares))
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1)):
                            old_index = square_index
                            new_index = left_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(left_first_square, left_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

                        right_squares = list(reversed(right_squares))
                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1)):
                            old_index = square_index
                            new_index = right_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(right_first_square, right_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

                    else:
                        for (index, square_index) in enumerate(range(top_first_square, top_last_square + 1)):
                            old_index = square_index
                            new_index = right_squares[index]
                            results.append((old_index,new_index))

                        top_squares = list(reversed(top_squares))
                        for (index, square_index) in enumerate(range(left_first_square, left_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = top_squares[index]
                            results.append((old_index,new_index))

                        for (index, square_index) in enumerate(range(down_first_square, down_last_square + 1)):
                            old_index = square_index
                            new_index = left_squares[index]
                            results.append((old_index,new_index))

                        down_squares = list(reversed(down_squares))
                        for (index, square_index) in enumerate(range(right_first_square, right_last_square + 1, self.size)):
                            old_index = square_index
                            new_index = down_squares[index]
                            results.append((old_index,new_index))

        elif side_name == "D":

            for turn in range(quarter_turns):

                # rotate the connecting row(s) of the surrounding sides
                for row in range(rows_to_rotate):
                    left_first_square = (self.squares_per_side * 2) - self.size + 1 - (row * self.size)
                    left_last_square = left_first_square + self.size - 1

                    front_first_square = (self.squares_per_side * 3) - self.size + 1 - (row * self.size)
                    front_last_square = front_first_square + self.size - 1

                    right_first_square = (self.squares_per_side * 4) - self.size + 1 - (row * self.size)
                    right_last_square = right_first_square + self.size - 1

                    back_first_square = (self.squares_per_side * 5) - self.size + 1 - (row * self.size)
                    back_last_square = back_first_square + self.size - 1

                    #log.info("left first %d, last %d" % (left_first_square, left_last_square))
                    #log.info("front first %d, last %d" % (front_first_square, front_last_square))
                    #log.info("right first %d, last %d" % (right_first_square, right_last_square))
                    #log.info("back first %d, last %d" % (back_first_square, back_last_square))

                    if reverse:
                        for square_index in range(left_first_square, left_last_square + 1):
                            old_index = square_index
                            new_index = square_index + self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(front_first_square, front_last_square + 1):
                            old_index = square_index
                            new_index = square_index + self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(right_first_square, right_last_square + 1):
                            old_index = square_index
                            new_index = square_index + self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(back_first_square, back_last_square + 1):
                            old_index = square_index
                            new_index = square_index - (3 * self.squares_per_side)
                            results.append((old_index,new_index))

                    else:
                        for square_index in range(left_first_square, left_last_square + 1):
                            old_index = square_index
                            new_index = square_index + (3 * self.squares_per_side)
                            results.append((old_index,new_index))

                        for square_index in range(front_first_square, front_last_square + 1):
                            old_index = square_index
                            new_index = square_index - self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(right_first_square, right_last_square + 1):
                            old_index = square_index
                            new_index = square_index - self.squares_per_side
                            results.append((old_index,new_index))

                        for square_index in range(back_first_square, back_last_square + 1):
                            old_index = square_index
                            new_index = square_index - self.squares_per_side
                            results.append((old_index,new_index))

        return results

Algorithm A

The most challenging part of carrying out this operation is the implementation of Algorithm A. While Knuth describes it in his AOCP Volume 3 in a straightfoward way, there are numerous details and subtleties and questions that come up that are not addressed in AOCP.

To do this, we implemented a bit vector that kept track of which top-bottom pairs had already been separated into a factor, and kept pulling out factors until every item in the original tuple had been visited.

Here is the algorithm implemented in Python:

def factor_permutation(perm_top,perm_bot):
    """
    Factor a permutation into its lowest terms
    """
    MAX = 96
    # Need a way to also mark them as used... bit vector
    used_vector = [0,]*len(perm_top)

    i = 0
    start = perm_top[0]
    used_vector[0] = 1

    factors = []

    # If we still have values to pick out:
    while(0 in used_vector):

        factor = []

        while(True):
            used_vector[i] = 1
            leader = perm_top[i]
            follower = perm_bot[i]

            i = perm_top.index(follower)
            while(used_vector[i]==1):
                i += 1
                if(i>=MAX):
                    break

            if(i>=MAX):
                break
            elif(follower==start):
                break
            else:
                factor.append(follower)

        # add start to end
        factor.append(start)

        factors.append(factor)
        try:
            #import pdb; pdb.set_trace()
            i = used_vector.index(0)
            start = perm_top[i]
        except ValueError:
            break

    factorsize = set()
    check = 0
    for factor in factors:
        factorsize.add(len(factor))
        check += len(factor)
    return factors



Flags