May 5, 2025
This post is the second of a two-part series that describes my computer-assisted solution to a physical puzzle I got. You will probably want to read the beginning of part I to understand what the puzzle is about, but other than that this post is pretty self-contained :)
Last time we defined a list allValidPieces of all the ways that a piece may
be fitted in the box. Its type is the nested list [[V3 Int]], because we
represent a single piece as a list of the voxels that compose it. Looks like
we're going to deal with lots of nested lists... Let's try to prevent confusion
by defining a type synonym:
type Shape = [V3 Int]
Shape will typically represent a puzzle piece, but it can represent any chunk
of voxels. We'll use it below to represent a set of pieces lumped together.
Ok, moving on. Let's recall what we are after: we need to find a way to fit 25
pieces in the box simultaneously. In other words, we're looking for a
subcollection of allValidPieces that has 25 elements, such that no two
elements interpenetrate one another. Which begs the question: what does it mean,
formally, for two shapes to interpenetrate one another?
Recall that we represented a piece as the list of voxels that compose it. Two shapes, then, are in conflict when they have at least one voxel in common — in other words, when there is any voxel from the first that is also an element of the second:
conflict :: Shape -> Shape -> Bool
conflict shape1 shape2 = any (\voxel -> elem voxel shape2) shape1
At this point, I'd say we're finally done modeling the problem and can start looking for a solution!
25 pieces is a lot, we can't just repeatedly pick lists of 25 pieces until we
come across one where pieces happen to not conflict with one another. We need to
build our solution step by step, and for that, we need a way to find out what
pieces are compatible with a collection of pieces already in place. Out of
performance concerns, we won't encode the collection of pieces already in place
using a list of Shapes, but we'll lump them all together in one single
Shape composed of the voxels that are already occupied by some piece:
newCompatiblePieces :: Shape -- ^ Voxels to avoid
-> [Shape] -- ^ All pieces that do not conflict with the input shape
newCompatiblePieces occupiedVoxels =
filter
(\candidatePiece -> not $ conflict candidatePiece occupiedVoxels)
allValidPieces
Note that in the above function, we are looking for all pieces that do not conflict with the input, but the pieces in the output list will generally conflict with one another. That's not a problem: the output list should not be seen as a list of pieces to add, but as a list of options which we can choose from.
Of course, some of the potentially compatible pieces we find, while they do not directly conflict with any piece already in place, will lead to a dead-end nevertheless: there is no solution compatible with both the pieces in place and the new piece. That's actually kind of the point of our puzzle!: we can't solve it just by adding piece after piece as long as they fit.
With an imperative apprach, we'd deal with this problem using explicit backtracking, keeping track, at each step, of the pieces we've tried, and going back one step whenever we don't find what we're looking for. This is hard to think about! In Haskell, we have a cheat code, and that's called:
List MonadI like to think of monads and do notation as a way for the programmer to pretend
they have access to a pure value when they don't actually. In the IO monad,
writing userInput <- getLine is a way to pretend you have a value of type
String, even though you need to wait for getLine to be executed at runtime
for the value to actually be defined. In the Maybe monad, writing
result <- someComputation is a way to pretend that someComputation yielded
a normal value, which you can access under the name result, even though
someComputation might be Nothing — in which case the whole do block becomes
Nothing.
Enter the List monad (also known as the [] monad, but that's harder to
pronounce). If someList is a list of type [a], you can write
value <- someList and proceed as though value were an actual value of type
a, even though there may be zero or more than one elements in the list.
What mechanism makes that possible? When you write value1 <- someList1, the
universe splits into several parallel subuniverses, and everything you write
under that line in the do block is computed/executed as many times as there are
elements in someList1. If you write value2 <- someList2 below that, each
of those parallel universes splits into parallel universes again, one for each
element of someList2. If you're not careful, the number of branches can explode
quickly!
With that in mind, we can implement our function that recursively looks for
subsolutions of the puzzle: lists of n compatible pieces (for some input n
between 0 and 25) that do not conflict with a collection of pieces (again,
lumped together into a single Shape) assumed to already be in place:
subsolutions :: Int -- ^ How many pieces should be in the subsolution we're looking for?
-> Shape -- ^ Unavailable voxels, already occupied by some piece
-> [[Shape]] -- ^ List of subsolutions, each of which is itself a list of pieces
subsolutions 0 _ = [[]]
subsolutions n occupiedVoxels = do
newPiece <- filter
(\piece -> not (conflict occupiedVoxels piece))
allValidPieces
let updatedOccupiedVoxels = newPiece <> occupiedVoxels
otherPieces <- subsolutions (n - 1) updatedOccupiedVoxels
return $ newPiece : otherPieces
The base case is n = 0: there is exactly one solution with zero pieces, that's the list with no piece in it. Now look at what happens in the other branch:
<-), that is, pieces
that do not conflict with our list of occupied voxels;<- symbol, we split the universe into as many parallel
universes as there are candidate pieces; in each subuniverse, newPiece
will be defined to be the corresponding candidate piece;Now you might think "gosh, that's a lot of universes to split into, how is my computer ever going to deal with those?", and you'd be right. There are two reasons for hope, though:
Anyway, the final list of solutions (there may be more than one!) to our problem is then:
allSolutions :: [[Shape]]
allSolutions = subsolutions 25 []
and we get an arbitrary solution by taking the first element of the list with
head.
Alas, my computer hangs and never returns a value. It seems that we were too careless splitting the universe many times, and the computer can't handle the resulting combinatorial explosion. Honestly, this isn't too surprising. Our approach is as bruteforce as it gets — which is actually good for a first approach: as Donald Knuth puts it, "premature optimization is the root of all evil (or at least most of it) in programming".
Ok, but what now? I guess we'll need to dig into optimization literature, or maybe turn the problem into a SAT formula and fire up a SAT solver... Just kidding! We're still going to bruteforce our problem, we just need to be smarter about it.
We don't need to look far, in the way we built our final list of solutions, to
find sources of redundancy — (sub)solutions that are equivalent to one another,
but are nevertheless registered separately. For two given compatible piece
dispositions piece1 and piece2, for instance, [piece1, piece2] and
[piece2, piece1] count as two subsolutions, even though they obviously
represent the same situation. The search will be performed once for each, and
if it leads nowhere, well... too bad, you will have wasted twice as much time.
And this will be true for any combination of compatible pieces, at any recursive
call. No wonder the algorithm is inefficient!
So, how are we to reduce the search space? Well, assume that you have a
subsolution and you're looking for a piece to add to it in order to get closer
to the full solution. Now take any free voxel freeVoxel, that is, any voxel
that is not part of any of the subsolution's pieces. Since we ultimately want
to fill the 5x5x5 cube, we know for sure that the full solution will feature a
piece that contains freeVoxel. So we can restrict our search to only pieces
that do!
Let's first implement a function that picks any free voxel in the box. We just
define a list composed of all of the box's voxels, filter out those that are
contained in one of the pieces in place, and take the first element of the
filtered list using head:
pickFreeVoxel :: Shape -> V3 Int
pickFreeVoxel alreadyOccupiedVoxels =
head $
filter (\voxel -> not $ elem voxel alreadyOccupiedVoxels) $
[V3 x y z | x <- [1..5], y <- [1..5], z <- [1..5]]
It isn't hard, now, to modify our subsolutions function so that it implements
our "smarter" bruteforce algorithm:
subsolutionsSmart :: Int -- ^ How many pieces should be in the subsolution we're looking for?
-> Shape -- ^ Unavailable voxels, already occupied by some piece
-> [[Shape]]
subsolutionsSmart 0 _ = [[]]
subsolutionsSmart n occupiedVoxels = do
let freeVoxel = pickFreeVoxel occupiedVoxels
newPiece <- filter
(\piece -> elem freeVoxel piece &&
not (conflict occupiedVoxels piece)
)
allValidPieces
let updatedOccupiedVoxels = newPiece <> occupiedVoxels
otherPieces <- subsolutionsSmart (n - 1) updatedOccupiedVoxels
return $ newPiece : otherPieces
The only change is that each time we're looking for a new piece, we pick a
freeVoxel and filter out pieces that don't contain it. Note that while
newPiece is on the left of a left arrow <-, because we want to be
considering all pieces in the list, freeVoxel is defined using a normal
let assignment. That's because, even though any free voxel will do, we only
want to consider one.
And finally, to get the full list of solutions, we do exactly the same as above:
allSolutionsSmart :: [[Shape]]
allSolutionsSmart = subsolutionsSmart 25 []
We only want one solution, so we print the first element (or head) of that
list, and we measure the time it takes:
main :: IO ()
main = do
startTime <- getCurrentTime
print $ head allSolutionsSmart
endTime <- getCurrentTime
putStrLn $ "Found in " <>
show (realToFrac (diffUTCTime endTime startTime) :: Double) <>
" seconds."
And my computer does find a solution, in a little more than 30 minutes!
I don't have a moral to this story, it was just a fun case study. But here are three suggestions of ways the program could be improved.
Optimization. On the one hand, the fact that a bruteforcing approach was able to find us a solution at all is pretty nice. On the other hand, 30 minutes is still quite a long time! The algorithm could probably be optimized, for instance by favoring subsolutions that are more "compact", for some definition of "compact"...
If you run the program at home, the output you'll get is a frankly illegible
list of lists of vector values that look like V3 2 4 1. I was able to test my
solution in practice by manually drawing the corresponding pieces, and then packing the pieces
according to my drawing, but I won't lie, the process was pretty tedious. So
another possible improvement could be to find a way to pretty-print the solution(s)
the algorithm finds. It's not obvious how to do so, especially with text output!
Finally, you might want to get all possible solutions, and for that, to get rid of solutions that are rotated and/or mirrored versions of other solutions. Better yet, you could try to avoid generating them in the first place (in the same way that we found a trick to avoid generating solutions that were the same up to a rearrangement of the piece order).
Feel free to tell me if you make progress on either front!