This is the one hundred and fifth part of the ILP series. For your convenience you can find other parts in the table of contents in Part 1 – Boolean algebra

Today we are going to solve the Sudoku Cube. It’s a variation of the regular Rubik’s cube in which we have digits on the cube and each side must have exactly nine different digits from one to nine. What’s more, digits on a side must all look the same direction.

I took photos of my cube. You can see them below:

The tricky part is what to do. With the regular cube it’s rather clear where each piece should go. The hard part is figuring out how to put it into place. With Sudoku Cube it’s much harder. We don’t actually see where to put each element. The first step is to “solve the grid”. We need to understand where to put elements.

# Solving the grid

First, notice that the number `1`

is symmetrical. We don’t know whether it’s facing “up” or “down”. Therefore, we need to assume it can look both directions. We could actually wonder similarly for number eight (which side is “up”) or six/nine, but for them we can make some assumptions. Eight should be narrower at the top. Nine has a label on it.

Second, we don’t know which directions the sides look into. Also, it may not be obvious from the very beginning, but each center has an orientation. This is different from the regular cube. There are algorithms to do so and we’ll see them later. However, for now we need to assume that any side can look any direction. Sometimes the cube is convenient and the “middle strap” sides look upwards, but we can’t assume that.

Therefore, we need to figure out where pieces go. Once we have that, we can then solve the grid.

We start by defining the permutations in this way:

This looks cryptic, so let’s analyze step by step.

Let’s start with the corners marked with red arrows and numbers. We can see six faces of the cube. Going from left to right, top to bottom, the faces are: up, left, front, right, back, down. We have eight corners of the cube. Therefore, there are eight pieces that we can rotate and put in given places. Places are numbered from oen to eight, and permutations are numbered from zero to two. We are going to encode each piece as a list of numbers and list of rotations.

Let’s take the Picture from the above in which you can see free sides of the cube. Let’s say that the side with numbers `3, 9, 4, 6, 1, 5, 5, 7, 9`

is the front side, and the side with numbers `5, 4, 6, 3, 6, 5, 8, 2, 2`

is the top one. You can see that we have a corner with number `4`

to the right of `9`

and above `5`

. This is the corner that we marked as `1_0`

on the diagram screen. So we can see, that we can encode this corner as having one with the following sequence of numbers: `4, 8, 4`

. Therefore `1_0 = 4, 1_1 = 8, 1_2 = 4`

.

You can see the arrow going up. This is an arbitrary indication which face I consider “upwards” of this given field. Similarly, for other fields. We want to encode the rotation of a piece based on the orientation of the digit. We have four directions: aligned with the arrow encoded as zero, rotated clockwise once by ninety degrees encoded as one, rotated twice encoded as two, rotated three times encoded as three. Therefore, the rotations of the numbers `4, 8, 4`

are `3, 2, 1`

.

We can now take each corner and encode it properly. This should be the encoding:

1 2 3 4 5 6 7 8 9 10 |
var corners = new []{ new Corner(solver, "c1") { Values = new[]{4, 8, 4}, Rotations = new[] {3, 2, 1}}, new Corner(solver, "c2") { Values = new[]{9, 5, 1}, Rotations = new[] {3, 2, 1}}, new Corner(solver, "c3") { Values = new[]{5, 9, 1}, Rotations = new[] {2, 2, 1}}, new Corner(solver, "c4") { Values = new[]{3, 7, 5}, Rotations = new[] {1, 0, 3}}, new Corner(solver, "c5") { Values = new[]{2, 6, 7}, Rotations = new[] {0, 3, 0}}, new Corner(solver, "c6") { Values = new[]{6, 3, 7}, Rotations = new[] {2, 1, 2}}, new Corner(solver, "c7") { Values = new[]{4, 1, 6}, Rotations = new[] {3, 0, 0}}, new Corner(solver, "c8") { Values = new[]{8, 9, 2}, Rotations = new[] {3, 1, 0}} }; |

We can do the same for edges (middles on the sides). We have twelve edges, each edge has two pieces. They are marked with green numbers and arrows. This is how we encode the cube:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
var edges = new []{ new Edge(solver, "e1") { Values = new[]{9, 3}, Rotations = new[] {0, 3}}, new Edge(solver, "e2") { Values = new[]{5, 2}, Rotations = new[] {0, 2}}, new Edge(solver, "e3") { Values = new[]{7, 6}, Rotations = new[] {1, 2}}, new Edge(solver, "e4") { Values = new[]{9, 8}, Rotations = new[] {2, 0}}, new Edge(solver, "e5") { Values = new[]{3, 2}, Rotations = new[] {0, 2}}, new Edge(solver, "e6") { Values = new[]{6, 9}, Rotations = new[] {1, 3}}, new Edge(solver, "e7") { Values = new[]{2, 3}, Rotations = new[] {1, 3}}, new Edge(solver, "e8") { Values = new[]{8, 5}, Rotations = new[] {3, 3}}, new Edge(solver, "e9") { Values = new[]{8, 8}, Rotations = new[] {1, 2}}, new Edge(solver, "e10") { Values = new[]{4, 2}, Rotations = new[] {1, 1}}, new Edge(solver, "e11") { Values = new[]{3, 4}, Rotations = new[] {2, 0}}, new Edge(solver, "e12") { Values = new[]{1, 7}, Rotations = new[] {3, 0}} }; |

We also have centers. We don’t care about their orientation at this point.

# Formulas

We basically want to do the following:

- For each corner: assign exactly one position to a corner (position from one to eight) and one permutation (from zero to two)
- Do the same for edges (twelve positions and two permutations)
- Make sure that each face has all pieces but center oriented the same way
- Make sure that number one can be oriented upside down (as it looks the same)
- Make sure each face has digits from one to nine

Sounds easy. Let’s see the first implementation.

# First implementation

Let’s start by encoding the grid. First, we’d like to encode the positions from the chart above. For corners, we take the value like `4_1`

and serialize it to one number as `4 * 3 + 1 = 13`

. We do similarly for edges:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
// 0 means empty // Values are 1_0, 1_1, 1_2, 2_0, 2_1, ... // Corners have 3 elements, so: 1_0 = 3 * 1 + 0 = 3 // Edges have 2 elements, so: 1_0 = 2 * 1 + 0 = 2 var diagram = new []{ new[] {0, 0, 0, 16, 17, 26, 0, 0, 0, 0, 0, 0}, new[] {0, 0, 0, 23, 0, 11, 0, 0, 0, 0, 0, 0}, new[] {0, 0, 0, 14, 3, 4, 0, 0, 0, 0, 0, 0}, new[] {17, 22, 13, 12, 2, 3, 5, 10, 25, 24, 16, 15}, new[] {21, 0, 9, 8, 0, 4, 5, 0, 14, 15, 0, 20}, new[] {19, 24, 11, 9, 6, 6, 7, 12, 23, 21, 18, 18}, new[] {0, 0, 0, 10, 7, 8, 0, 0, 0, 0, 0, 0}, new[] {0, 0, 0, 25, 0, 13, 0, 0, 0, 0, 0, 0}, new[] {0, 0, 0, 20, 19, 22, 0, 0, 0, 0, 0, 0}, }; |

Now, we encode the orientation of each field (arrows). We assume `-1`

means empty field:

1 2 3 4 5 6 7 8 9 10 11 |
var fieldsOrientations = new []{ new[] {-1, -1, -1, 3, 0, 0, -1, -1, -1, -1, -1, -1}, new[] {-1, -1, -1, 3, 0, 1, -1, -1, -1, -1, -1, -1}, new[] {-1, -1, -1, 2, 2, 1, -1, -1, -1, -1, -1, -1}, new[] {3, 0, 0, 3, 0, 0, 3, 0, 0, 3, 0, 0}, new[] {3, 0, 1, 3, 0, 1, 3, 0, 1, 3, 0, 1}, new[] {2, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 1}, new[] {-1, -1, -1, 3, 0, 0, -1, -1, -1, -1, -1, -1}, new[] {-1, -1, -1, 3, 0, 1, -1, -1, -1, -1, -1, -1}, new[] {-1, -1, -1, 2, 2, 1, -1, -1, -1, -1, -1, -1} }; |

Now, let’s make sure that each piece is put in a different position:

1 2 3 4 5 |
// Make corners in distinct positions solver.Set<AllDifferent>(corners[0].FinalPosition, corners.Skip(1).Select(c => c.FinalPosition).ToArray()); // Make edges in distinct positions solver.Set<AllDifferent>(edges[0].FinalPosition, edges.Skip(1).Select(e => e.FinalPosition).ToArray()); |

Let’s now start the magic. The idea is as follows: we ask the solver to assign positions and rotations to the pieces. We then iterate over all possibilities, check if a given possibility was selected by the solver, and then calculate the outcome. We then add the requirements to make the grid solved. Let’s begin with the corners.

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 |
Func<int, int, Tuple<IVariable, IVariable>> setupSingleCorner = (x, y) => { int number = diagram[x][y]; int position = number / 3 - 1; int permutation = number % 3; IVariable piecePermutation = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(2)); IVariable finalPermutation = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(2)); IVariable pieceRotationValue = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(3)); IVariable finalRotation = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(3)); IVariable cornerValue = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(9)).Set<GreaterOrEqual>(solver.FromConstant(1)); foreach(var corner in corners){ var matchingCorner = solver.Operation<IsEqual>(corner.FinalPosition, solver.FromConstant(position)); solver.Operation<MaterialImplication>( matchingCorner, solver.Operation<IsEqual>(corner.FinalPermutation, piecePermutation) ).Set<Equal>(solver.FromConstant(1)); for(int possiblePermutation=0;possiblePermutation<3;++possiblePermutation){ var matchingPermutation = matchingCorner.Operation<Conjunction>(finalPermutation.Operation<IsEqual>(solver.FromConstant(possiblePermutation))); solver.Operation<MaterialImplication>( matchingPermutation, solver.Operation<IsEqual>(pieceRotationValue, solver.FromConstant(corner.Rotations[possiblePermutation])) ).Set<Equal>(solver.FromConstant(1)); solver.Operation<MaterialImplication>( matchingPermutation, solver.Operation<IsEqual>(cornerValue, solver.FromConstant(corner.Values[possiblePermutation])) ).Set<Equal>(solver.FromConstant(1)); } } solver.Set<Equal>( finalPermutation, solver.Operation<Remainder>( solver.Operation<Addition>(solver.FromConstant(permutation), piecePermutation), solver.FromConstant(3) ) ); solver.Set<Equal>( finalRotation, solver.Operation<Remainder>( solver.Operation<Addition>(solver.FromConstant(fieldsOrientations[x][y]), pieceRotationValue), solver.FromConstant(4) ) ); return Tuple.Create(cornerValue, finalRotation); }; |

For each corner from the diagram, we deserialize its position and permutation (lines 2-4). Next, we need to prepare some variables (lines 6-10) and then decipher the solution.

We iterate over each possible piece (line 12) and see if the corner was assigned to the field we analyze now (line 13). If that’s the case, then we copy the value of the assigned permutation to the temporary variable (lines 15-18).

Next, we do some magic. We would like to check all possible permutations to see which one was selected. To do that, we need to calculate what we call here `finalPermutation`

. However, this is some kind of a cyclic dependency. First, we copy the permutation assigned by the solver and called `corner.FinalPermutation`

to the variable `piecePermutation`

(lines 15-18). Next, we need to add the permutation of the field to that value, calculate the modulo three, and assign that to the variable `finalPermutation`

(lines 35-41). However, we then (or technically earlier) use the variable `finalPermutation`

in line 21 to answer whether the current permutation we analyze is the one that was selected. Once we have that, we can copy the rotation of the corner (lines 23-26) and the value of the corner (lines 28-32).

Next, we do another time travel and we calculate the final orientation of the field by adding the orientation of the corner and the orientation (direction of the arrow) of the field. That’s in lines 43-49.

Finally, we return both the value and the final orientation (line 51).

Okay, we can now iterate over all corners of a given face:

1 2 3 4 5 6 7 8 9 10 |
Func<int, int, Tuple<IVariable, IVariable>[]> setupCorners = (x, y) => { var valuesAndRotations = new []{ setupSingleCorner(x, y), setupSingleCorner(x, y+2), setupSingleCorner(x+2, y+2), setupSingleCorner(x+2, y) }; return valuesAndRotations; }; |

We need to do a very similar code with the edges:

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 |
Func<int, int, Tuple<IVariable, IVariable>> setupSingleEdge = (x, y) => { int number = diagram[x][y]; int position = number / 2 - 1; int permutation = number % 2; IVariable piecePermutation = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(1)); IVariable finalPermutation = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(1)); IVariable pieceRotationValue = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(3)); IVariable finalRotation = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(3)); IVariable edgeValue = solver.CreateAnonymous(Domain.PositiveOrZeroInteger).Set<LessOrEqual>(solver.FromConstant(9)).Set<GreaterOrEqual>(solver.FromConstant(1)); calculatedRotations[x][y] = finalRotation; foreach(var edge in edges){ var matchingEdge = solver.Operation<IsEqual>(edge.FinalPosition, solver.FromConstant(position)); solver.Operation<MaterialImplication>( matchingEdge, solver.Operation<IsEqual>(edge.FinalPermutation, piecePermutation) ).Set<Equal>(solver.FromConstant(1)); for(int possiblePermutation=0;possiblePermutation<2;++possiblePermutation){ var matchingPermutation = matchingEdge.Operation<Conjunction>(finalPermutation.Operation<IsEqual>(solver.FromConstant(possiblePermutation))); solver.Operation<MaterialImplication>( matchingPermutation, solver.Operation<IsEqual>(pieceRotationValue, solver.FromConstant(edge.Rotations[possiblePermutation])) ).Set<Equal>(solver.FromConstant(1)); solver.Operation<MaterialImplication>( matchingPermutation, solver.Operation<IsEqual>(edgeValue, solver.FromConstant(edge.Values[possiblePermutation])) ).Set<Equal>(solver.FromConstant(1)); } } solver.Set<Equal>( finalPermutation, solver.Operation<Remainder>( solver.Operation<Addition>(solver.FromConstant(permutation), piecePermutation), solver.FromConstant(2) ) ); solver.Set<Equal>( finalRotation, solver.Operation<Remainder>( solver.Operation<Addition>(solver.FromConstant(fieldsOrientations[x][y]), pieceRotationValue), solver.FromConstant(4) ) ); return Tuple.Create(edgeValue, finalRotation); }; |

We iterate over all edges of a face:

1 2 3 4 5 6 7 8 9 10 |
Func<int, int, Tuple<IVariable, IVariable>[]> setupEdges = (x, y) => { var valuesAndRotations = new []{ setupSingleEdge(x, y+1), setupSingleEdge(x+1, y), setupSingleEdge(x+1, y+2), setupSingleEdge(x+2, y+1) }; return valuesAndRotations; }; |

And we also get the value of the center:

1 2 3 |
Func<int, int, IVariable> setupCenter = (x, y) => { return solver.FromConstant(int.Parse("" + centers[x][y])); }; |

We can now setup a single face:

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 |
Action<int, int> setupSingleFace = (x, y) => { var valuesAndRotations = setupCorners(x, y).Concat(setupEdges(x, y)).ToArray(); var valuesAndRotationsAndOnes = valuesAndRotations.Select(t => { var anotherPossibleRotation = solver.Operation<Condition>( t.Item1.Operation<IsEqual>(solver.FromConstant(1)), t.Item2.Operation<Addition>(solver.FromConstant(2)).Operation<Remainder>(solver.FromConstant(4)), t.Item2 ); return Tuple.Create(t.Item1, t.Item2, anotherPossibleRotation); }).ToArray(); // Make pairwise rotations equal = make pieces (but center) look the same direction for(int id=1; id < valuesAndRotationsAndOnes.Length; ++id){ solver.Set<Equal>( solver.Operation<Disjunction>( solver.Operation<IsEqual>(valuesAndRotationsAndOnes[id].Item2, valuesAndRotationsAndOnes[id-1].Item2), solver.Operation<IsEqual>(valuesAndRotationsAndOnes[id].Item2, valuesAndRotationsAndOnes[id-1].Item3), solver.Operation<IsEqual>(valuesAndRotationsAndOnes[id].Item3, valuesAndRotationsAndOnes[id-1].Item2), solver.Operation<IsEqual>(valuesAndRotationsAndOnes[id].Item3, valuesAndRotationsAndOnes[id-1].Item3) ), solver.FromConstant(1) ); } var center = setupCenter(x + 1, y+1); // Make sudoku numbers solver.Set<AllDifferent>(center, valuesAndRotationsAndOnes.Select(t => t.Item1).ToArray()); }; |

We take all the values and rotations of corners and edges (line 2). We then need to take all the ones and calculate additional rotation. We check if the value is one (line 6) and if that’s the case, then we add another rotation (upside-down), otherwise we return the same rotation (lines 7-8).

Finally, we take all the sides and make sure that they all have the same rotations. So for each pair of pieces, either their rotations are equal, or their rotation and alternative rotation are equal, or their alternative rotations are equal (lines 15-25).

Finally, we add the center (line 27) and make sure that all numbers are different (line 30) which means that we have the sudoku requirement met.

We can now configure all faces:

1 2 3 4 5 6 7 8 9 10 |
Action setupFaces = () => { setupSingleFace(0, 3); setupSingleFace(3, 0); setupSingleFace(3, 3); setupSingleFace(3, 6); setupSingleFace(3, 9); setupSingleFace(6, 3); }; setupFaces(); |

Let’s now add the goal, solve the problem, and print the result:

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 |
solver.AddGoal("goal", solver.FromConstant(0)); solver.Solve(); var template = new []{ " --- ", " | | ", " | | ", " | | ", " --- --- --- ---", "| | | | |", "| | | | |", "| | | | |", " --- --- --- --- ", " | | ", " | | ", " | | ", " --- " }; var numbersGrid = template.Select(x => x.ToArray()).ToArray(); var rotationsGrid = template.Select(x => x.ToArray()).ToArray(); Action<int, int, int, int> printSingleCorner = (x, y, destinationX, destinationY) => { int number = diagram[x][y]; int position = number / 3 - 1; int permutation = number % 3; var matchingCorner = corners.First(c => (int)Math.Round(solver.GetValue(c.FinalPosition)) == position); var piecePermutation = (int)Math.Round(solver.GetValue(matchingCorner.FinalPermutation)); var finalPermutation = (permutation + piecePermutation) % 3; numbersGrid[destinationX][destinationY] = matchingCorner.Values[finalPermutation].ToString()[0]; rotationsGrid[destinationX][destinationY] = matchingCorner.Rotations[finalPermutation].ToString()[0]; }; Action<int, int, int, int> printCorners = (x, y, destinationX, destinationY) => { printSingleCorner(x, y, destinationX, destinationY); printSingleCorner(x, y+2, destinationX, destinationY + 2); printSingleCorner(x+2, y+2, destinationX + 2, destinationY + 2); printSingleCorner(x+2, y, destinationX + 2, destinationY); }; Action<int, int, int, int> printSingleEdge = (x, y, destinationX, destinationY) => { int number = diagram[x][y]; int position = number / 2 - 1; int permutation = number % 2; var matchingEdge = edges.First(e => (int)Math.Round(solver.GetValue(e.FinalPosition)) == position); var piecePermutation = (int)Math.Round(solver.GetValue(matchingEdge.FinalPermutation)); var finalPermutation = (permutation + piecePermutation) % 2; numbersGrid[destinationX][destinationY] = matchingEdge.Values[finalPermutation].ToString()[0]; rotationsGrid[destinationX][destinationY] = matchingEdge.Rotations[finalPermutation].ToString()[0]; }; Action<int, int, int, int> printEdges = (x, y, destinationX, destinationY) => { printSingleEdge(x, y+1, destinationX, destinationY + 1); printSingleEdge(x+1, y, destinationX + 1, destinationY); printSingleEdge(x+1, y+2, destinationX + 1, destinationY + 2); printSingleEdge(x+2, y+1, destinationX + 2, destinationY + 1); }; Action<int, int, int, int> printCenters = (x, y, destinationX, destinationY) => { numbersGrid[destinationX+1][destinationY+1] = centers[x+1][y+1]; }; Action<int, int, int, int> printSingleFace = (x, y, destinationX, destinationY) => { printCorners(x, y, destinationX, destinationY); printEdges(x, y, destinationX, destinationY); printCenters(x, y, destinationX, destinationY); }; Action printFaces = () => { printSingleFace(0, 3, 1, 5); printSingleFace(3, 0, 5, 1); printSingleFace(3, 3, 5, 5); printSingleFace(3, 6, 5, 9); printSingleFace(3, 9, 5, 13); printSingleFace(6, 3, 9, 5); }; printFaces(); for(int row = 0; row < template.Length; ++ row){ Console.Write(new String(numbersGrid[row])); Console.Write(" "); Console.WriteLine(new String(rotationsGrid[row])); } |

This works. In theory. I tried running that with NEOS for eight hours and it didn’t solve the problem, unfortunately. It was simply too complex.

Let’s see another solution.

# Second implementation

In this solution we’ll apply many optimizations. We’ll store the positions as bitmasks instead of integers to use boolean flags (which are much faster). We won’t calculate the obtained direction of the face, but we’ll enforce which face must be selected. Finally, we won’t check the options and see if they were selected, but rather we will calculate what was selected end enforce the options. Let’s start.

First, we create a bitmask for directions of each face. We have six faces (sides) and four directions for each:

1 2 3 |
// face0_up, face0_right, face0_down, face0_left, face1_up ... // front = 0, right = 1, back = 2, left = 3, up = 4, down = 5 var facesDirections = Enumerable.Range(0, 6 * 4).Select(d => solver.CreateAnonymous(Domain.BinaryInteger)).ToArray(); |

Next, we encode a bitmask what value was assigned to each corner side. We have six faces, four corners for each face, and nine numbers to choose from:

1 2 3 |
// face0_corner_top_right_value1, face0_corner_top_right_value2, ..., face0_corner_top_right_value9, face0_corner_bottom_right_value1, // ..., face0_corner_bottom_right_value9, face0_corner_bottom_left_value1, ..., face0_corner_top_left_value_9, face1_corner_top_right_value1... var cornersNumbers = Enumerable.Range(0, 6*4*9).Select(n => solver.CreateAnonymous(Domain.BinaryInteger)).ToArray(); |

We do the similar for edges:

1 2 3 4 |
// face0_top_value1, face_0_top_value2, ..., face0_top_value9, face0_right_value1 // ..., face0_right_value9, face0_bottom_value1, ..., face0_bottom_value9, face0_left_value1, ... // face0_left_value9, face1_top_value1... var edgeNumbers = Enumerable.Range(0, 6*4*9).Select(n => solver.CreateAnonymous(Domain.BinaryInteger)).ToArray(); |

Now, we can use these bitmasks to constrain the solution. We make sure that digits on a face look all the same direction. So, for each face we take the direction flags and make sure that only one flag is selected:

1 2 3 4 |
// Each face has exactly one direction for(var faceId = 0; faceId < 6; ++faceId){ solver.Operation<Addition>(facesDirections.Skip(faceId * 4).Take(4).ToArray()).Set<Equal>(solver.FromConstant(1)); } |

We then make sure each corner has exactly one number. We have six faces and four corners on each face:

1 2 3 4 |
// Each corner has exactly one number for(var cornerId = 0; cornerId < 6 * 4; ++cornerId){ solver.Operation<Addition>(cornersNumbers.Skip(cornerId * 9).Take(9).ToArray()).Set<Equal>(solver.FromConstant(1)); } |

We now need to see how corners are stored. Previously, `corner.FinalPosition`

was an integer. Now it’s a long bitmasks of positions and permutations. We have eight possible positions and three permutations:

1 2 3 4 5 6 7 8 9 10 11 12 |
public class Corner$id { public int[] Values {get;set;} public int[] Rotations {get;set;} public IVariable[] FinalPosition {get; set;} public string Prefix {get;set;} public Corner(IMilpManager solver, string prefix){ // position1_permutation0, position1_permutation1, position1_permutation2, position2_permutation0, ... FinalPosition = Enumerable.Range(0, 8*3).Select(d => solver.CreateAnonymous(Domain.BinaryInteger)).ToArray(); Prefix = prefix; } } |

Next, we make sure that each corner piece has exactly one position and one rotation selected:

1 2 3 4 |
// Each corner is in one place foreach(var corner in corners){ solver.Operation<Addition>(corner.FinalPosition).Set<Equal>(solver.FromConstant(1)); } |

Next we need to make sure that each corner of the cube has exactly one piece assigned. We have eight corners, for each corner we take bit flags of each piece, and make sure that only one was selected:

1 2 3 4 5 6 |
// Exactly one piece in given corner for(var cornerId = 0; cornerId < 8; ++cornerId){ solver.Operation<Addition>( corners.SelectMany(corner => corner.FinalPosition.Skip(cornerId * 3).Take(3)).ToArray() ).Set<Equal>(solver.FromConstant(1)); } |

We do the same for edges:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
// Each edge has exactly one number for(var edgeId = 0; edgeId < 6 * 4; ++edgeId){ solver.Operation<Addition>(edgeNumbers.Skip(edgeId * 9).Take(9).ToArray()).Set<Equal>(solver.FromConstant(1)); } // Each edge is in one place foreach(var edge in edges){ solver.Operation<Addition>(edge.FinalPosition).Set<Equal>(solver.FromConstant(1)); } // Exactly one piece in given edge for(var edgeId = 0; edgeId < 12; ++edgeId){ solver.Operation<Addition>( edges.SelectMany(edge => edge.FinalPosition.Skip(edgeId * 2).Take(2)).ToArray() ).Set<Equal>(solver.FromConstant(1)); } |

Now we need to make sure that each side has sudoku numbers:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
var faceToCenters = new Dictionary<int, Tuple<int, int>>{ {0, Tuple.Create(4, 4)}, {1, Tuple.Create(4, 7)}, {2, Tuple.Create(4, 10)}, {3, Tuple.Create(4, 1)}, {4, Tuple.Create(1, 4)}, {5, Tuple.Create(7, 4)} }; // Sudoku numbers on each face for(var faceId = 0; faceId < 6; ++faceId){ var centerPosition = faceToCenters[faceId]; var centerNumber = int.Parse("" + centers[centerPosition.Item1][centerPosition.Item2]); var cornerDigits = Enumerable.Range(0, 4).Select(p => cornersNumbers.Skip(faceId * 4 * 9).Skip(p * 9).Take(9).ToArray()).ToArray(); var edgeDigits = Enumerable.Range(0, 4).Select(p => edgeNumbers.Skip(faceId * 4 * 9).Skip(p * 9).Take(9).ToArray()).ToArray(); var centerDigit = new [] { Enumerable.Range(0, 9).Select(i => i == centerNumber - 1 ? solver.FromConstant(1) : solver.FromConstant(0)).ToArray() }; var vectors = cornerDigits.Concat(edgeDigits).Concat(centerDigit).ToArray(); for(var digit = 0; digit < 9; ++ digit){ solver.Operation<Addition>(vectors.Select(v => v[digit]).ToArray()).Set<Equal>(solver.FromConstant(1)); } } |

We need two additional mappings from corners and edges to faces:

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 |
var cornersToFaces = new Dictionary<Tuple<int, int>, int> { {Tuple.Create(0, 3), 4}, {Tuple.Create(0, 5), 4}, {Tuple.Create(2, 3), 4}, {Tuple.Create(2, 5), 4}, {Tuple.Create(3, 0), 3}, {Tuple.Create(3, 2), 3}, {Tuple.Create(5, 0), 3}, {Tuple.Create(5, 2), 3}, {Tuple.Create(3, 3), 0}, {Tuple.Create(3, 5), 0}, {Tuple.Create(5, 3), 0}, {Tuple.Create(5, 5), 0}, {Tuple.Create(3, 6), 1}, {Tuple.Create(3, 8), 1}, {Tuple.Create(5, 6), 1}, {Tuple.Create(5, 8), 1}, {Tuple.Create(3, 9), 2}, {Tuple.Create(3, 11), 2}, {Tuple.Create(5, 9), 2}, {Tuple.Create(5, 11), 2}, {Tuple.Create(6, 3), 5}, {Tuple.Create(6, 5), 5}, {Tuple.Create(8, 3), 5}, {Tuple.Create(8, 5), 5}, }; var edgesToFaces = new Dictionary<Tuple<int, int>, int> { {Tuple.Create(0, 4), 4}, {Tuple.Create(1, 3), 4}, {Tuple.Create(1, 5), 4}, {Tuple.Create(2, 4), 4}, {Tuple.Create(3, 1), 3}, {Tuple.Create(4, 0), 3}, {Tuple.Create(4, 2), 3}, {Tuple.Create(5, 1), 3}, {Tuple.Create(3, 4), 0}, {Tuple.Create(4, 3), 0}, {Tuple.Create(4, 5), 0}, {Tuple.Create(5, 4), 0}, {Tuple.Create(3, 7), 1}, {Tuple.Create(4, 6), 1}, {Tuple.Create(4, 8), 1}, {Tuple.Create(5, 7), 1}, {Tuple.Create(3, 10), 2}, {Tuple.Create(4, 9), 2}, {Tuple.Create(4, 11), 2}, {Tuple.Create(5, 10), 2}, {Tuple.Create(6, 4), 5}, {Tuple.Create(7, 3), 5}, {Tuple.Create(7, 5), 5}, {Tuple.Create(8, 4), 5}, }; |

We also need to know where each corner is (whether it’s right top, right bottom, left bottom, left top). Similarly for edges (top, right, bottom, left):

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 |
var cornersToDiagonalEnds = new Dictionary<Tuple<int, int>, int> { {Tuple.Create(0, 3), 3}, {Tuple.Create(0, 5), 0}, {Tuple.Create(2, 3), 2}, {Tuple.Create(2, 5), 1}, {Tuple.Create(3, 0), 3}, {Tuple.Create(3, 2), 0}, {Tuple.Create(5, 0), 2}, {Tuple.Create(5, 2), 1}, {Tuple.Create(3, 3), 3}, {Tuple.Create(3, 5), 0}, {Tuple.Create(5, 3), 2}, {Tuple.Create(5, 5), 1}, {Tuple.Create(3, 6), 3}, {Tuple.Create(3, 8), 0}, {Tuple.Create(5, 6), 2}, {Tuple.Create(5, 8), 1}, {Tuple.Create(3, 9), 3}, {Tuple.Create(3, 11), 0}, {Tuple.Create(5, 9), 2}, {Tuple.Create(5, 11), 1}, {Tuple.Create(6, 3), 3}, {Tuple.Create(6, 5), 0}, {Tuple.Create(8, 3), 2}, {Tuple.Create(8, 5), 1}, }; var edgesToAxisEnds = new Dictionary<Tuple<int, int>, int> { {Tuple.Create(0, 4), 0}, {Tuple.Create(1, 3), 3}, {Tuple.Create(1, 5), 1}, {Tuple.Create(2, 4), 2}, {Tuple.Create(3, 1), 0}, {Tuple.Create(4, 0), 3}, {Tuple.Create(4, 2), 1}, {Tuple.Create(5, 1), 2}, {Tuple.Create(3, 4), 0}, {Tuple.Create(4, 3), 3}, {Tuple.Create(4, 5), 1}, {Tuple.Create(5, 4), 2}, {Tuple.Create(3, 7), 0}, {Tuple.Create(4, 6), 3}, {Tuple.Create(4, 8), 1}, {Tuple.Create(5, 7), 2}, {Tuple.Create(3, 10), 0}, {Tuple.Create(4, 9), 3}, {Tuple.Create(4, 11), 1}, {Tuple.Create(5, 10), 2}, {Tuple.Create(6, 4), 0}, {Tuple.Create(7, 3), 3}, {Tuple.Create(7, 5), 1}, {Tuple.Create(8, 4), 2}, }; |

With all that prepared we can finally assert the corners:

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 |
Action<int, int> setupSingleCorner = (x, y) => { int number = diagram[x][y]; int position = number / 3 - 1; int permutation = number % 3; var face = cornersToFaces[Tuple.Create(x, y)]; var diagonalEnd = cornersToDiagonalEnds[Tuple.Create(x, y)]; foreach(var corner in corners){ for(var piecePermutation = 0; piecePermutation < 3; ++piecePermutation){ var wasPicked = corner.FinalPosition.Skip(position * 3).Skip(piecePermutation).First(); var finalPermutation = (permutation + piecePermutation) % 3; var finalRotation = (corner.Rotations[finalPermutation] + fieldsOrientations[x][y]) % 4; var value = corner.Values[finalPermutation]; if(value == 1){ solver.Operation<MaterialImplication>( wasPicked, solver.Operation<Disjunction>( facesDirections.Skip(face * 4).Take(4).ToArray()[finalRotation], facesDirections.Skip(face * 4).Take(4).ToArray()[(finalRotation + 2)%4] ) ).Set<Equal>(solver.FromConstant(1)); }else{ solver.Operation<MaterialImplication>( wasPicked, facesDirections.Skip(face * 4).Take(4).ToArray()[finalRotation] ).Set<Equal>(solver.FromConstant(1)); } solver.Operation<MaterialImplication>( wasPicked, cornersNumbers.Skip(face * 4 * 9).Skip(diagonalEnd * 9).Take(9).ToArray()[value - 1] ).Set<Equal>(solver.FromConstant(1)); } } }; |

We take some metadata about the piece we are in (lines 2-6). Next, we iterate over all corners and permutations. However, this time instead of asking “is this corner in that permutation assigned to this field”, we simply take what was assigned and enforce the proper bitflags.

We first take the bit indicating whether a given piece was assigned to this corner and this permutation (line 10). Next, we calculate the final permutation in this field (line 11) and the final rotation including directions of the arrows (line 12). We also take the value of the piece (line 13). Notice that this time we don’t extract these values from the ILP variables. We precalculate them outside of the model which is much faster.

We can now do two simple constraints. First, if this digit is one (line 15) then we enforce this constraint: if this piece was picked (line 17), then direction of the face needs to be either the same as the direction of the piece (line 19) or the opposite (line 20). If the digit is not one, then we enforce the constraint, that if the piece was picked, than the boolean flag indicating the direction of the face must be selected (line 26).

The crucial difference is that we don’t calculate the direction within the model now. We do the causation outside of the model. We encode something like “I don’t know what piece was selected, but if the piece I’m currently looking at was selected, then this bit flag of the face direction must be selected too”.

The next constraint we add is for the numbers: if the piece was picked (line 31) then the number assigned to the corner is selected (line 32). Again, we don’t care about the actual value if it was selected or not. We just make sure that if it was selected, then some other bit flag is set.

We can now setup all corners:

1 2 3 4 5 6 |
Action<int, int> setupCorners = (x, y) => { setupSingleCorner(x, y); setupSingleCorner(x, y+2); setupSingleCorner(x+2, y+2); setupSingleCorner(x+2, y); }; |

We do very similar stuff for edges:

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 |
Action<int, int> setupSingleEdge = (x, y) => { int number = diagram[x][y]; int position = number / 2 - 1; int permutation = number % 2; var face = edgesToFaces[Tuple.Create(x, y)]; var axisEnd = edgesToAxisEnds[Tuple.Create(x, y)]; foreach(var edge in edges){ for(var piecePermutation = 0; piecePermutation < 2; ++piecePermutation){ var wasPicked = edge.FinalPosition.Skip(position * 2).Skip(piecePermutation).First(); var finalPermutation = (permutation + piecePermutation) % 2; var finalRotation = (edge.Rotations[finalPermutation] + fieldsOrientations[x][y]) % 4; var value = edge.Values[finalPermutation]; if(value == 1){ solver.Operation<MaterialImplication>( wasPicked, solver.Operation<Disjunction>( facesDirections.Skip(face * 4).Take(4).ToArray()[finalRotation], facesDirections.Skip(face * 4).Take(4).ToArray()[(finalRotation + 2)%4] ) ).Set<Equal>(solver.FromConstant(1)); }else{ solver.Operation<MaterialImplication>( wasPicked, facesDirections.Skip(face * 4).Take(4).ToArray()[finalRotation] ).Set<Equal>(solver.FromConstant(1)); } solver.Operation<MaterialImplication>( wasPicked, edgeNumbers.Skip(face * 4 * 9).Skip(axisEnd * 9).Take(9).ToArray()[value - 1] ).Set<Equal>(solver.FromConstant(1)); } } }; Action<int, int> setupEdges = (x, y) => { setupSingleEdge(x, y+1); setupSingleEdge(x+1, y); setupSingleEdge(x+1, y+2); setupSingleEdge(x+2, y+1); }; |

Finally, we encode faces:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
Action<int, int> setupSingleFace = (x, y) => { setupCorners(x, y); setupEdges(x, y); }; Action setupFaces = () => { setupSingleFace(0, 3); setupSingleFace(3, 0); setupSingleFace(3, 3); setupSingleFace(3, 6); setupSingleFace(3, 9); setupSingleFace(6, 3); }; setupFaces(); |

We run it and this is the result we get:

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 |
Tried aggregator 2 times. MIP Presolve eliminated 7285 rows and 3504 columns. MIP Presolve modified 1149 coefficients. Aggregator did 3596 substitutions. Reduced MIP has 979 rows, 844 columns, and 4885 nonzeros. Reduced MIP has 844 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.03 sec. (15.56 ticks) Found incumbent of value 0.000000 after 0.03 sec. (24.88 ticks) Root node processing (before b&c): Real time = 0.03 sec. (25.19 ticks) Parallel b&c, 12 threads: Real time = 0.00 sec. (0.00 ticks) Sync time (average) = 0.00 sec. Wait time (average) = 0.00 sec. ------------ Total (root+branch&cut) = 0.03 sec. (25.19 ticks) --- --- |239| |333| |568| |3 3| |147| |333| --- --- --- --- --- --- --- --- |829|523|681|529| |222|111|111|111| |157|618|974|347| |2 2|1 1|1 1|1 1| |634|497|532|681| |222|111|111|111| --- --- --- --- --- --- --- --- |863| |111| |219| |1 1| |457| |111| --- --- |

It took 30 milliseconds to find the answer. First implementation couldn’t do that in 8 hours. Not to mention that the first implementation ran on Gurobi in NEOS, and the second one was running with CPLEX on my local machine.

# Reading the solution

Let’s read the solution. We can see that when we hold the cube in a way that the top side has center 6, front side has center 1, and bottom side has center 1, then the top side looks to the left, left side looks down, front and right and back and down sides look to the right. We can also see all the numbers assigned, so we know the “colors” of the pieces. We can now apply the regular algorithms for solving the cube.

The last remaining piece is how to rotate centers. There are two algorithms for that:

- Top center clockwise and front center counter-clockwise: f b’ l r’ u d’ f’ u’ d l’ r f’ b u
- Top center 180 degrees: u r l uu r’ l’ u r l uu r’ l’

We can now easily solve the cube:

# Lessons learned

There are some important tricks we applied to optimize the solution:

- We stored numbers as bitflags instead of storing them as integers. This is a common trick. If you know the range of the variables, then don’t use integers. Bitflags are much faster.
- We precalcualted results outside of the model. Instead of creating a variable with a value of the face rotation, we calculated the rotation outside of the model, and then added a constraint (basically, an “if” condition in the ILP)
- We didn’t compare numbers. Generally, don’t do that if it’s possible. Instead, just create a boolean flag indicating the result of comparison and use it later onw. The difference is subtle but improves the performance a lot.
- We didn’t use “if else” conditions. They look easy, but they are really the performance killer.

With all these improvements in place, we managed to solve the cube in milliseconds.