# The H48 optimal solver

This document contains information on the H48 Rubik's Cube optimal solver.
The implementation of the solver is still in progress. This document
partly describes ideas that have not been implemented yet, and it will
be updated to reflect the actual implementation.

I highly encourage the reader to check out Jaap Scherphuis'
[Computer Puzzling page](https://www.jaapsch.net/puzzles/compcube.htm)
before reading this document.

Other backround material that will be referenced throughout the text includes:

* [Cube coordinates](https://sebastiano.tronto.net/speedcubing/coordinates)
  by Sebastiano Tronto
* [Nxopt](https://github.com/rokicki/cube20src/blob/master/nxopt.md)
  by Tomas Rokicki
* [The Mathematics behind Cube Explorer](https://kociemba.org/cube.htm)
  by Herber Kociemba

Initially I intended to give minimum information here, citing the
above resources when needed, but I ended up rewriting even some of the
basic things in this document.

## Optimal solvers

### IDA*

The basic idea behind the solver is an iterative-deepening
depth-first search with pruning tables, also known as IDA*. You
can find a detailed explanation of this idea in
[Jaap's page](https://www.jaapsch.net/puzzles/compcube.htm#tree).

To summarize, the IDA* algorithm finds the shortest solution(s) for a
scrambled puzzle by using multiple depth-first searches at increasing
depth. A breadth-first search is not applicable in this scenario because
of the large amount of memory required.  At each depth, the puzzle's state
is evaluated to obtain an estimated lower bound for how many moves are
required to solve the puzzle; this is done by employing large pruning
tables, as well as other techniques. If the lower bound exceeds the
number of moves required to reach the current target depth, the branch is
pruned. Otherwise, every possible move is tried, except those that would
"cancel" with the preceding moves, obtaining new positions at depth N+1.

### Pruning tables and coordinates

A pruning table associates to a cube position a value that is a
lower bound for the number of moves it takes to solve that position.
To acces this value, the cube must be turned into an index for the
given table. This is done using a
[*coordinate*](https://sebastiano.tronto.net/speedcubing/coordinates),
by which we mean a function from the set of (admissible, solvable) cube
positions to the set of non-negative integers. The maximum value N of a
coordinate must be less than the size of the pruning table; preferably
there are no "gaps", that is the coordinate reaches all values between
0 and N and the pruning table has size N+1.

Coordinates are often derived from the cosets of a *target subgroup*.
As an example, consider the pruning table that has one entry for each
possible corner position, and the stored values denote the length of
an optimal corner solution for that position. In this situation, the
target subgroup is the set of all cube positions with solved corners
(for the Mathematician: this is indeed a subgroup of the group of all
cube positions). This group consists of `12! * 2^10` positions.  The set
of cosets of this subgroup can be identified with the set of corner
configurations, and it has size `8! * 3^7`.  By looking at corners,
from each cube position we can compute a number from `0` to `8! * 3^7 - 1`
and use it as an index for the pruning table.

The pruning table can be filled in many ways. A common way is starting
with the coordinate of the solved cube (often 0) and visiting the set
of all coordinates with a breadth-first search. To do this, one needs
to apply a cube move to a coordinate; this can either be done directly
(for example using transistion tables, as explained in Jaap's page) or
by converting the coordinate back into a valid cube position, applying
the move and then computing the coordinate of the new position.

The largest a pruning table is, the more information it contains,
the faster the solvers that uses it is. It is therefore convenient to
compress the tables as much as possible. This can be achieved by using
4 bits per position to store the lower bound, or as little as 2 bits
per position with some caveats (see Jaap's page or Rokicki's nxopt
description for two possible ways to achieve this).  Tables that use
only 1 bit per position are possible, but the only information they
can provide is wether the given position requires more or less than a
certain number of moves to be solved.

### Symmetries

Another, extremely effective way of reducing the size of a pruning table
is using symmetries. In the example above of the corner table, one can
notice that all corner positions that are one quarter-turn away from being
solved are essentially the same: they can be turned one into the other
by rotating the cube and optionally applying a mirroring transformation.
Indeed, every position belongs to a set of up to 48 positions that are
"essentially the same" in this way.  Keeping a separate entry in the
pruning table for each of these positions is a waste of memory.

There are ways to reduce a coordinate by symmetry, grouping every
transformation-equivalent position into the same index. You can
find a description in my
[cube coordinates page](https://sebastiano.tronto.net/speedcubing/coordinates).
By doing so, the size of the pruning table is reduced without loss of
information.

Some well-known solvers (Cube Explorer, nxopt) do not take advantage
of the full group of 48 symmetries of the cube, but they use only the
16 symmetries that fix the UD axis. However, they make up for this by
doing 3 pruning table lookups, one for each axis of the cube.

Some cube positions are self-symmetric: when applying certain
transformations to them, they remain the same. For example, the
cube obtained by applying U2 to a solved cube is invariant with
respect to the mirroring along the RL-plane. This fact has a couple
of consequences:

* Reducing by a group of N symmetries does not reduce the size of
  the pruning table by a factor of N, but by a slightly smaller
  factor. If the size of the coordinate is large, this fact is
  harmless.
* When combining symmetry-reduced coordinates with other coordinates,
  one has to be extra careful with handling self-symmetric coordinates.
  An informal explanation of this phenomenon can be found in
  the [docs/nasty-symmetries.md document](./nasty-symmetries.md).

## The H48 solver

### The target subgroup: H48 coordinates

The H48 solver uses a target group that is invariant under all 48
symmetries. This group is defined as follows:

* Each corner is placed in the *corner tetrad* it belongs to, and it
  is oriented. Here by corner tetrad we mean one of the two sets of
  corners {UFR, UBL, DFL, DBR} and {UFL, UBR, DFR, DBL}. For a corner
  that is placed in its own tetrad, the orientation does not depend on
  the reference axis chosen, and an oriented corner can be defined
  has having the U or D sticker facing U or D.
* Each edge is placed in the *slice* it belongs to. The three edge slices
  are M = {UF, UB, DB, DF}, S = {UR, UL, DL, DR} and E = {FR, FL, BL, BR}.

Other options are available for corners: for example, we could
impose that the corner permutation is even or that corners are in
*half-turn reduction* state, that is they are in the group generated
by {U2, D2, R2, L2, F2, B2}. We settled on the corner group described
above because it is easy to calculate and it gives enough options
for pruning table sizes, depending on how edge orientation is
considered (see the next section).

We call the coordinate obtained from the target subgroup described
above the **h0** coordinate. This coordinate has `C(8,4) * 3^7 *
C(12,8) * C(8,4) = 5.304.568.500` possible values.  If the corner
part of the coordinate is reduced by symmetry, it consists of only
3393 elements, giving a total of `3393 * C(12,8) * C(8,4) =
117.567.450` values. This is not very large for a pruning table,
as with 2 bits per entry it would occupy less than 30MB. However,
it is possible to take edge orientation (partially) into account
to produce larger tables.

### Edge orientation

Like for corners, the orientation of an edge is well-defined independently
of any axis when said edge is placed in the slice it belongs to. We may
modify the target subgroup defined in the previous section by imposing
that the edges are all oriented. This yields a new coordinate, which we
call **h11**, whose full size with symmetry-reduced corners is `3393 *
C(12,8) * C(8,4) * 2^11 = 240.778.137.600`. A pruning table based on
this coordinate with 2 bits per entry would occupy around 60GB, which
is a little too much for most personal computers.

One can wonder if it is possible to use a coordinate that considers
the orientation of only *some* of the edges, which we may call **h1**
to **h10**. Such coordinates do exist, but they are not invariant under
the full symmetry group: indeed an edge whose orientation we keep track
of could be moved to any of the untracked edges' positions by one of
the symmetries, making the whole coordinate ill-defined.

It is however possible to compute the symmetry-reduced pruning tables
for these coordinates. One way to construct them is by taking the **h11**
pruning table and "forgetting" about some of the edge orientation values,
collapsing 2 (or a power thereof) values to one by taking the minimum.
It is also possible to compute these tables directly, as explained in
the **Pruning table computation** section below.

### Coordinate computation for pruning value estimation

In order to access the pruning value for a given cube position, we first
need to compute its coordinate value. Let's use for simplicity the **h0**
coordinate, but everything will be valid for any other coordinate of
the type described in the previous section.

As described in the symmetric-composite coordinate section of
[my coordinates page](https://sebastiano.tronto.net/speedcubing/coordinates),
we first need to compute the part of the coordinate where the symmetry
is applied, that is the corner separation + orientation (also called
*cocsep*). The value of the coordinate depends on the equivalence
class of the corner configuration, and is memorized in a table called
`cocsepdata`. To avoid lengthy computations, the index for a cube position
in this table is determined by the corner orientation (between `0` and
`3^7-1`) and a binary representation of the corner separation part (an
8-bit number with 4 zero digits and 4 one digits, where zeros denote
corners in one tetrad and 1 denotes corners in the other; the last bit
can be ignored, as it can easily be deduced from the other 7). This is
slightly less space-efficient than computing the actual corner-separation
coordinate as a value between `0` and `C(8,4)`, but it is much faster.

We can now access the value in the `cocsepdata` table corresponding to
our cube position. This table contains the value of the symmetry-reduced
coordinate and the so-called *ttrep* (transformation to representative)
necessary to compute the full **h0** coordinate. For convenience,
this table also stores a preliminary pruning value based on the corner
coordinate.  If this value is large enough, we may skip the computation of
the **h0** coordinate and any further estimate. These three values can be
stored in a single 32 bit integer, so that the table uses less than 12MB.

### Getting the pruning value

Once we have computed the full h0 coordinate, we can access the correct
entry in the full pruning table. As mentioned above, the pruning table
can be one of three kinds:

* 4 bits per entry, or `k4`: In this case the pruning value (between 0
  and 15) can be simply read off the table.
* 2 bits per entry, or `k2`: Tables of this kind work as described by
  Rokicki in the
  [nxopt document](https://github.com/rokicki/cube20src/blob/master/nxopt.md).
  In this case the pruning table also has a *base value*, that determines
  the offset to be added to each entry (each entry can only be 0, 1, 2 or 3).
  If the base value is `b`, a pruning value of 1, 2 or 3 can be used directly
  as a lower bound of b+1, b+2 and b+3 respectively. However, a value of 0
  could mean that the actual lower bound is anything between 0 and b, so we
  cannot take b as a lower bound. Instead we have to use a pruning value from
  another table - see the section "Fallback tables" below.
* 1 bit per entry, or `k1`: With one bit per entry, the only information we
  can get from the pruning table is wether or not the current position
  requires more or fewer moves than a fixed base value b. This can still be
  valuable if most positions are more or less equally split between two
  pruning values.
  (Work in progress - `k1` tables not available in the code yet)

### Fallback tables

When a pruning table does not store the exact lower bound value, for
example in the case of `k2` table as described above, we need to refine
our estimate using another table, which we call *fallback table*.
In the current implementation, we actually use two different tables:

* **h0** table: for larger `k2` tables we get a fallback value from the full
  table for the **h0** coordinate.
* edges-only table: to improve the pruning value for certain specific
  scrambles, namely whose with solved corners such as the superflip,
  we employ a second table that takes into account the edges of a full
  **h11** coordinate. This table is small (around 1MB).

More tables could be used to refine the fallback estimate, but each
additional table leads to longer lookup times, especially if it is
too large to fit in cache.

### Estimation refinements

After computing the pruning value, there are a number of different tricks
that can be used to improve the estimation.

#### Inverse estimate

A cube position and its inverse will, in general, give a different
pruning value. If the normal estimate is not enough to prune the branch,
the inverse of the position can be computed and a new estimate can be
obtained from there.

#### Reducing the branching factor - tight bounds

*This trick and the next are and adaptation of a similar technique
introduced by Rokicki in nxopt.*

We can take advantage of the fact that the **h0** (and **h11**) coordinate
is invariant under the subgroup `<U2, D2, R2, L2, F2, B2>`. Suppose we
are looking for a solution of length `m`, we are at depth `d` and the
pruning value `p` *for the inverse position* is a strict bound, that is
`d+p = m`.  In this situation, from the inverse point of view the solution
cannot end with a 180° move: if this was the case, since these moves are
contained in the target group for **h0** (or **h11**), it must be that
we were in the target subgroup as well *before the last move*, i.e. we
have found a solution for the **h0** coordinate of length `p-1`. But this
is in contradiction with the fact that the inverse pruning value is `p`.

From all of this, we conclude that trying any 180° move as next move is
useless (because it would be the last move from the inverse position).
We can therefore reduce the branching factor significantly and only try
quarter-turn moves.

#### Reducing the branching factor - switching to inverse

We can expand on the previous trick by using a technique similar to NISS.

Suppose that we have a strict bound for the *normal position*. As above,
we can deduce that the last move cannot be a 180° move, but this does not
tell us anything about the possibilities for the next move. However, if
we *switch to the inverse scramble* we can take advantage of this fact
as described above. For doing this, we need to replace the cube with its
inverse, and keep track of the moves done from now on so that we can
invert them at the end to construct the final solution.

When this technique is used, it is also possible to avoid some table
lookups: if the last move applied to the cube is a 180° move *on the
inverse position*, then the coordinate on the normal position has not
changed. Thus if we keep track of the last computed pruning value,
we can reuse it and avoid an expensive table lookup.

### Other optimizations

Other possible (low-level) optimizations include:

* **Avoid inverse computation**: computing the inverse of a cube position
  is expensive. We can avoid doing that (for the inverse pruning value
  estimate) if we bring along both the normal and the inverse cube during
  the search, and we use *premoves* to apply moves to the inverse scramble.
* **Multi-threading (multiple scrambles)**: It is easy to parallelize this
  algorithm when solving multiple cubes at once, by firing up multiple
  instances of the solver. It is important to make sure that the same
  (read-only) pruning table is used for all instances, to avoid expensive
  memory duplication.
* **Multi-threading (single scramble)**: It is also possible to parallelize
  the search for a single scramble. For example, we can generate 18 different
  cubes, one for each possible starting move, and solve each of them in a
  separate thread. Some coordination between threads is necessary to stop
  the search when the desired number of solutions has been found.

## Pruning table computation

### 4 bits tables for h0 and h11

Computing the pruning table for a "real" h48 coordinate, that **h0**
or **h11**, using 4 bits per entry is quite simple. The method currently
implemented works as follows.

First, we set the value of the solved cube's coordinate to 0 and every
other entry to 15, the largest 4-bit integer. Then we iteratively scan
through the table and compute the coordinates at depth n+1 from those at
depth n as follows: for each coordinate at depth n, we compute a valid
representative for it, we apply each of the possible 18 moves to it, and
we set their value to n+1. Once the table is filled or we have reached
depth 15, we stop (there are no **h0** coordinates at depth 16 or more,
but I currently don't know if this is the case for **h11**).

Unfortunately what I described above is an oversimplification: one
also must take into account the case where the corner coordinate is
self-symmetric, and handle it accordingly by making sure that every
symmetric variation of the same coordinate has its value set at the
right iteration.

When this algorithm reaches the last few iterations, the coordinates at
depth n are many more than those whose depth is still unknown. Therefore
it is convenient to look at the unset coordinate and check if they have
any neighbor whose depth is i; if so, we can set such a coordinate to i+1.
Thanks to [torchlight](https://github.com/torchlight) for suggesting
this optimization.

Finally, this algorithm can be easily parallelized by dividing the set
of coordinates into separate sections, but one must take into account
that a coordinate and its neighbors are usually not in the same section.

This method is currently implemented only for **h0**, since the 4 bits
table for **h11** would require around 115GB of RAM to compute, and I
don't have this much memory.

### 2 bits tables with for h0 and h11

(Work in progress - currently there is no specialized routine for
computing 2 bits table for "real" coordinates; instead, an optimized
version of the generic method explained below is used)

### A generic method for intermediate coordinates (from h1 to h10)

(Work in progress - this method is quite slow and it may be replaced
by a better algorithm in the future)

Computing the pruning tables for intermediate coordinates is not as
simple.  The reason is that these coordinates are not invariant under
the 48 transformations of the cube, but we are treating them as such.
Take for example the case of **h10**. Each **h10** coordinate is
represented by either of two **h11** coordinates: the one where the 11th
edge is correctly oriented and the one where this edge is flipped (of
course, the orientation of the 12th edge can be deduced from the parity
of the orientation of the other 11). We can take any cube whose **h11**
coordinate is one of these two to represent our **h10** coordinate, but
the set of neighbors will change depending on which representative we
picked. In other words: if I take one of the two cubes and make a move,
the position reached is not necessarily obtainable by applying a move
to the other cube. This means that the same algorithm that we described
for the "real coordinate" case cannot be applied here.

It is still possible to do a brute-force depth-first seach to fill
these tables.  To make it reasonably fast, we apply a couple of
optimizations. First of all, we restrict ourselves to computing a 2 bits
table, so we can stop the search early.

The second optimization we employ consists in pre-computing all possible
**h11** coordinates at a fixed depth, and storing their depth in a hash
map indexed by their **h11** coordinate value. This way we do not have
to brute-force our way through from depth 0, and it make this method
considerably faster. Unfortunately, since the number of coordinates
at a certain depth increases exponentially with the depth, we are for
now limited at storing the coordinates at depth 8. (Work in progress -
we may experiment with depth 9 in the future)

This method works for the "real" coordinates **h0** and **h11** as well.
Moreover, in this case one can optimize it further by avoiding to repeat
the search from a coordinate that has already been visited. (Work
in progress - this method will be replaced in the future by a more
efficient one)