From 9ac266c76f39620d8343e46ca41cb09d1534384c Mon Sep 17 00:00:00 2001
From: Sebastiano Tronto <sebastiano@tronto.net>
Date: Sat, 7 Dec 2024 16:38:40 +0100
Subject: Merge the "solver-experiments" branch that I have been working on for
 a few weeks. This include mainly three things:

1. Various tweaks for a total performance gain of around 30%.
2. Take into account symmetries and avoid repeated work. This
   required a re-work of the splitting into tasks before the solve.
3. Add a second fallback table (eoesep). This gives huge performance
   gains for particular scrambles (e.g. superflip).

After merging this commit, remove and re-generate all pruning tables.

Squashed commit of the following:

commit 60f0705d2d69050e6a30581a2810f686d6f69b80
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Sat Dec 7 16:06:48 2024 +0100

    Fix indentation

commit cc5d489a251812b6188c0ba264ac6cb2236f1afe
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Sat Dec 7 15:56:19 2024 +0100

    Updated documentation

commit a3f605dd628546e52564f82b139feb473b0725f3
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Sat Dec 7 14:01:11 2024 +0100

    use eoesep table as second fallback

commit c75e43c9116c64f97e92b0fe038be8a032925a72
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Fri Dec 6 16:13:23 2024 +0100

    First commit for gendata_eoesep

commit fea7688ab8bdc5ae0c3480622e2510a3bcd39248
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Dec 3 17:31:00 2024 +0100

    Add scramble to tool

commit 66866cb71dea4ca8278ecb9e90ff4295771feb42
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Dec 3 17:22:59 2024 +0100

    Added tool to check multiple solutions

commit ef65611c772c3996bddca8d181da3538e0af1674
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Dec 3 17:16:20 2024 +0100

    Write all solutions for symmetric positions

commit e3ded26db7d7d4ae7c0e2488151ed14f581f7b8e
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Nov 26 09:14:54 2024 +0100

    Added symmetry filter (TODO: print excluded solutions)

commit 864c437a9751c58d58562650ca9eba4a9e6ad3eb
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Mon Nov 25 14:51:13 2024 +0100

    Improved task split

commit b88926d36d7ab0c64c5fe3bb954fd15d41267fba
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Fri Nov 22 19:06:41 2024 +0100

    Reworked tasks for multi-threading in view of symmetry filter

commit 26fa653f97df8cd601aecb80eaf89f0a00e9ba9f
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Thu Oct 31 15:37:43 2024 +0100

    Added transform move

commit 19f655ef94d658eaa2fefb5cea3c167a3ec58db6
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Thu Oct 31 09:29:15 2024 +0100

    Clarified doc

commit 3b0fe1e5ef8b628854e30f0f0067300e2763c954
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Thu Oct 31 08:36:58 2024 +0100

    Handle solved cube correctly

commit 57705cbc4982e3abe97a36ed64871738d4f721c0
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Thu Oct 31 08:24:30 2024 +0100

    Close file

commit fc7d462b58bcf3d3a3fbf26c1f3bd04c640a4898
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Oct 29 15:05:48 2024 +0100

    Removed stats tool

commit 359bf7cb49ef405ee76ed662207d47cb2abcc5a9
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Oct 29 15:01:51 2024 +0100

    Updated theory doc

commit 39c315af562bc4ce896f41004388a4c34d475d37
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Oct 29 14:51:40 2024 +0100

    Remove unused constants

commit 57a5d24538aa59a4df9221dad2f99e9f0286bd9d
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Tue Oct 29 10:20:38 2024 +0100

    Add tool to solve scrambles from file

commit 07e2918c216636891b1fa6adecc9756086a901e9
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Mon Oct 28 17:00:00 2024 +0100

    Add make table to tool

commit f5e5266c654eb027a5a35c57cc618555246f5e5e
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Mon Oct 28 09:35:49 2024 +0100

    Remove old solver, other small things

commit a1ec78025b7959dbb845213f7f4e6851ecebc204
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Sun Oct 27 02:00:29 2024 +0200

    Improvements

commit 8eea23dbe888d923e662e24ae969130e2c67b999
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Sat Oct 26 12:24:19 2024 +0200

    Makefile fix

commit 3fc3927beacc78971cefeb42da8d71fe6c015fc1
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Fri Oct 25 18:25:09 2024 +0200

    More performance gains

commit 7b4efa1f9af9722de1ab9ccfc27899825a0d12c4
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date:   Fri Oct 25 15:53:20 2024 +0200

    Alternative solver implementation, small performance gain
---
 src/arch/avx2.h                        |   6 +
 src/arch/common.h                      |  17 +
 src/arch/neon.h                        |  11 +
 src/arch/portable.h                    |  12 +
 src/core/moves.h                       |  54 +--
 src/core/transform.h                   |  40 +++
 src/nissy.c                            |  80 ++---
 src/solvers/h48/coordinate_macros.h    |   1 +
 src/solvers/h48/gendata_cocsep.h       |  34 +-
 src/solvers/h48/gendata_eoesep.h       | 274 +++++++++++++++
 src/solvers/h48/gendata_h48.h          |  62 +++-
 src/solvers/h48/gendata_types_macros.h |  11 +-
 src/solvers/h48/h48.h                  |   3 +-
 src/solvers/h48/solve.h                | 614 ++++++++++++++++++++++++---------
 src/solvers/h48/solve_multithread.h    | 350 -------------------
 src/solvers/h48/stats.h                | 100 ------
 src/utils/constants.h                  |  62 +++-
 17 files changed, 1004 insertions(+), 727 deletions(-)
 create mode 100644 src/solvers/h48/gendata_eoesep.h
 delete mode 100644 src/solvers/h48/solve_multithread.h
 delete mode 100644 src/solvers/h48/stats.h

(limited to 'src')

diff --git a/src/arch/avx2.h b/src/arch/avx2.h
index 2c3371e..b6ff510 100644
--- a/src/arch/avx2.h
+++ b/src/arch/avx2.h
@@ -16,6 +16,12 @@
 #define SOLVED_CUBE STATIC_CUBE( \
     0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
 
+STATIC_INLINE int
+popcount_u32(uint32_t x)
+{
+	return _mm_popcnt_u32(x);
+}
+
 STATIC void
 pieces(cube_t *cube, uint8_t c[static 8], uint8_t e[static 12])
 {
diff --git a/src/arch/common.h b/src/arch/common.h
index 802790c..52e0fc7 100644
--- a/src/arch/common.h
+++ b/src/arch/common.h
@@ -1,3 +1,5 @@
+STATIC_INLINE int popcount_u32(uint32_t);
+
 STATIC void pieces(cube_t *, uint8_t [static 8], uint8_t [static 12]);
 STATIC_INLINE bool equal(cube_t, cube_t);
 STATIC_INLINE cube_t invertco(cube_t);
@@ -18,6 +20,7 @@ STATIC_INLINE void set_eo(cube_t *, int64_t);
 STATIC_INLINE cube_t invcoord_esep(int64_t);
 
 STATIC_INLINE void invcoord_esep_array(int64_t, int64_t, uint8_t[static 12]);
+STATIC_INLINE cube_t invcoord_eoesep(int64_t);
 
 STATIC_INLINE void
 invcoord_esep_array(int64_t set1, int64_t set2, uint8_t mem[static 12])
@@ -44,3 +47,17 @@ invcoord_esep_array(int64_t set1, int64_t set2, uint8_t mem[static 12])
 		mem[i] = (slice[s]++) | (uint8_t)(s << 2);
 	}
 }
+
+STATIC_INLINE cube_t
+invcoord_eoesep(int64_t i)
+{
+	cube_t c;
+	int64_t esep, eo;
+
+	esep = i >> INT64_C(11);
+	eo = i % POW_2_11;
+	c = invcoord_esep(esep);
+	set_eo(&c, eo);
+
+	return c;
+}
diff --git a/src/arch/neon.h b/src/arch/neon.h
index c31c8f9..47dc147 100644
--- a/src/arch/neon.h
+++ b/src/arch/neon.h
@@ -26,6 +26,17 @@ STATIC_INLINE uint8x8_t compose_corners_slim(uint8x8_t, uint8x8_t);
 #define SOLVED_CUBE STATIC_CUBE(	\
 	0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
 
+/* TODO: optimize this (use intrinsics?) */
+STATIC_INLINE int
+{
+	int ret;
+
+	for (ret = 0; x != 0; x >>= 1)
+		ret += x & 1;
+
+	return ret;
+}
+
 STATIC void
 pieces(cube_t *cube, uint8_t c[static 8], uint8_t e[static 12])
 {
diff --git a/src/arch/portable.h b/src/arch/portable.h
index e9a8dd0..b8c5211 100644
--- a/src/arch/portable.h
+++ b/src/arch/portable.h
@@ -9,6 +9,18 @@
 #define SOLVED_CUBE STATIC_CUBE( \
     0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
 
+/* TODO: optimize this (use bit tricks?) */
+STATIC_INLINE int
+popcount_u32(uint32_t x)
+{
+	int ret;
+
+	for (ret = 0; x != 0; x >>= 1)
+		ret += x & 1;
+
+	return ret;
+}
+
 STATIC void
 pieces(cube_t *cube, uint8_t c[static 8], uint8_t e[static 12])
 {
diff --git a/src/core/moves.h b/src/core/moves.h
index e3a9bc1..7d4972e 100644
--- a/src/core/moves.h
+++ b/src/core/moves.h
@@ -2,8 +2,8 @@
 #define PREMOVE(M, c) compose(MOVE_CUBE_ ## M, c)
 
 STATIC_INLINE bool allowednextmove(uint8_t *, uint8_t);
+STATIC_INLINE uint32_t allowednextmove_mask(uint8_t *, uint8_t);
 
-STATIC_INLINE uint8_t inverse_trans(uint8_t);
 STATIC_INLINE uint8_t movebase(uint8_t);
 STATIC_INLINE uint8_t moveaxis(uint8_t);
 STATIC_INLINE uint32_t disable_moves(uint32_t, uint8_t);
@@ -38,26 +38,38 @@ STATIC cube_t applymoves(cube_t, const char *);
 STATIC bool
 allowednextmove(uint8_t *moves, uint8_t n)
 {
-	uint8_t base[3], axis[3];
+	return n == 0 ? true :
+	    allowednextmove_mask(moves, n-1) & (1 << moves[n-1]);
+}
+
+STATIC uint32_t
+allowednextmove_mask(uint8_t *moves, uint8_t n)
+{
+	uint32_t result;
+	uint8_t base1, base2, axis1, axis2;
+
+	result = MM_ALLMOVES;
 
-	if (n < 2)
-		return true;
+	if (n == 0)
+		return result;
 
-	base[0] = movebase(moves[n-1]);
-	axis[0] = moveaxis(moves[n-1]);
-	base[1] = movebase(moves[n-2]);
-	axis[1] = moveaxis(moves[n-2]);
+	base1 = movebase(moves[n-1]);
+	axis1 = moveaxis(moves[n-1]);
+	result = disable_moves(result, base1 * 3);
 
-	if (base[0] == base[1] || (axis[0] == axis[1] && base[0] < base[1]))
-		return false;
+	if (base1 % 2)
+		result = disable_moves(result, (base1 - 1) * 3);
 
-	if (n == 2)
-		return true;
+	if (n == 1)
+	return result;
 
-	base[2] = movebase(moves[n-3]);
-	axis[2] = moveaxis(moves[n-3]);
+	base2 = movebase(moves[n-2]);
+	axis2 = moveaxis(moves[n-2]);
 
-	return axis[1] != axis[2] || base[0] != base[2];
+	if(axis1 == axis2)
+		result = disable_moves(result, base2 * 3);
+
+	return result;
 }
 
 STATIC_INLINE uint32_t 
@@ -66,12 +78,6 @@ disable_moves(uint32_t current_result, uint8_t base_index)
 	return current_result & ~(7 << base_index);
 }
 
-STATIC_INLINE uint8_t
-inverse_trans(uint8_t t)
-{
-	return inverse_trans_table[t];
-}
-
 STATIC_INLINE uint8_t
 movebase(uint8_t move)
 {
@@ -84,6 +90,12 @@ moveaxis(uint8_t move)
 	return move / 6;
 }
 
+STATIC_INLINE uint8_t
+moveopposite(uint8_t move)
+{
+	return movebase(move) == 2 * moveaxis(move) ? move + 3 : move - 3;
+}
+
 STATIC cube_t
 move(cube_t c, uint8_t m)
 {
diff --git a/src/core/transform.h b/src/core/transform.h
index 73636fe..f956718 100644
--- a/src/core/transform.h
+++ b/src/core/transform.h
@@ -21,6 +21,9 @@ STATIC cube_t transform_edges(cube_t, uint8_t);
 STATIC cube_t transform_corners(cube_t, uint8_t);
 STATIC cube_t transform(cube_t, uint8_t);
 STATIC cube_t applytrans(cube_t, const char *);
+STATIC_INLINE uint8_t inverse_trans(uint8_t);
+STATIC uint8_t transform_move(uint8_t, uint8_t);
+STATIC uint64_t symmetry_mask(cube_t);
 
 STATIC cube_t
 transform_edges(cube_t c, uint8_t t)
@@ -354,3 +357,40 @@ applytrans(cube_t cube, const char *buf)
 
 	return transform(cube, t);
 }
+
+STATIC_INLINE uint8_t
+inverse_trans(uint8_t t)
+{
+	return inverse_trans_table[t];
+}
+
+STATIC uint8_t
+transform_move(uint8_t m, uint8_t t)
+{
+	uint8_t a, base, modifier;
+
+	a = moveaxis(m);
+	base = trans_move_table[t][a];
+	if (movebase(m) != 2 * a)
+		base = moveopposite(base);
+
+	modifier = m % 3;
+	if (t >= TRANS_UFm)
+		modifier = 2 - modifier;
+
+	return base + modifier;
+}
+
+STATIC uint64_t
+symmetry_mask(cube_t cube)
+{
+	uint64_t t, ret;
+	cube_t transformed;
+
+	for (t = 0, ret = 0; t < 48; t++) {
+		transformed = transform(cube, t);
+		ret |= ((uint64_t)equal(cube, transformed)) << t;
+	}
+
+	return ret;
+}
diff --git a/src/nissy.c b/src/nissy.c
index c09c99e..f933cab 100644
--- a/src/nissy.c
+++ b/src/nissy.c
@@ -12,12 +12,15 @@
 #include "core/core.h"
 #include "solvers/solvers.h"
 
-int parse_h48_solver(const char *, uint8_t [static 1], uint8_t [static 1]);
+long long parse_h48_solver(
+    const char *, uint8_t [static 1], uint8_t [static 1]);
 STATIC bool checkdata(const char *, const tableinfo_t *);
 STATIC bool distribution_equal(const uint64_t [static INFO_DISTRIBUTION_LEN],
     const uint64_t [static INFO_DISTRIBUTION_LEN], uint8_t);
 STATIC long long write_result(cube_t, char [static NISSY_SIZE_B32]);
 STATIC size_t my_strnlen(const char *, size_t); 
+STATIC long long nissy_gendata_unsafe(
+    const char *, unsigned long long, char *);
 
 #define GETCUBE_OPTIONS(S, F) { .option = S, .fix = F }
 struct {
@@ -28,19 +31,13 @@ struct {
 	GETCUBE_OPTIONS(NULL, NULL)
 };
 
-int
+long long
 parse_h48_solver(const char *buf, uint8_t h[static 1], uint8_t k[static 1])
 {
 	const char *fullbuf = buf;
 
 	buf += 3;
 
-	if (!strcmp(buf, "stats")) {
-		*h = 0;
-		*k = 4;
-		return 0;
-	}
-
 	if (*buf != 'h')
 		goto parse_h48_solver_error;
 	buf++;
@@ -62,9 +59,9 @@ parse_h48_solver(const char *buf, uint8_t h[static 1], uint8_t k[static 1])
 parse_h48_solver_error:
 	*h = 0;
 	*k = 0;
-	LOG("Error parsing solver: must be in \"h48h*k*\" format"
-	    " or \"h48stats\", but got %s\n", fullbuf);
-	return -1;
+	LOG("Error parsing solver: must be in \"h48h*k*\" format,"
+	    " but got %s\n", fullbuf);
+	return NISSY_ERROR_INVALID_SOLVER;
 }
 
 STATIC bool
@@ -359,20 +356,6 @@ nissy_getcube(
 	return write_result(c, result);
 }
 
-long long
-nissy_datasize(
-	const char *solver
-)
-{
-	if (solver == NULL) {
-		LOG("Error: 'solver' argument is NULL\n");
-		return NISSY_ERROR_NULL_POINTER;
-	}
-
-	/* gendata() handles a NULL *data as a "dryrun" request */
-	return nissy_gendata(solver, 0, NULL);
-}
-
 long long
 nissy_datainfo(
         uint64_t data_size,
@@ -422,6 +405,15 @@ nissy_datainfo(
 	return NISSY_OK;
 }
 
+long long
+nissy_datasize(
+	const char *solver
+)
+{
+	/* gendata() handles a NULL *data as a "dryrun" request */
+	return nissy_gendata_unsafe(solver, 0, NULL);
+}
+
 long long
 nissy_gendata(
 	const char *solver,
@@ -429,7 +421,17 @@ nissy_gendata(
 	char data[data_size]
 )
 {
-	int p;
+	return nissy_gendata_unsafe(solver, data_size, data);
+}
+
+STATIC long long
+nissy_gendata_unsafe(
+	const char *solver,
+	unsigned long long data_size,
+	char *data
+)
+{
+	long long parse_ret;
 	gendata_h48_arg_t arg;
 
 	if (solver == NULL) {
@@ -440,10 +442,10 @@ nissy_gendata(
 	arg.buf_size = data_size;
 	arg.buf = data;
 	if (!strncmp(solver, "h48", 3)) {
-		p = parse_h48_solver(solver, &arg.h, &arg.k);
+		parse_ret = parse_h48_solver(solver, &arg.h, &arg.k);
 		arg.maxdepth = 20;
-		if (p != 0)
-			return NISSY_ERROR_UNKNOWN;
+		if (parse_ret != NISSY_OK)
+			return parse_ret;
 		return gendata_h48(&arg);
 	} else {
 		LOG("gendata: unknown solver %s\n", solver);
@@ -495,7 +497,7 @@ nissy_solve(
 )
 {
 	cube_t c;
-	int p;
+	long long parse_ret;
 	uint8_t h, k;
 
 	if (solver == NULL) {
@@ -521,20 +523,12 @@ nissy_solve(
 	}
 
 	if (!strncmp(solver, "h48", 3)) {
-		if (!strcmp(solver, "h48stats"))
-			return solve_h48stats(c, maxmoves, data, sols);
-
-		p = parse_h48_solver(solver, &h, &k);
-		if (p != 0) {
-			LOG("solve: unknown solver %s\n", solver);
-			return NISSY_ERROR_INVALID_SOLVER;
-		} else {
-			return THREADS > 1 ?
-			    solve_h48_multithread(c, minmoves, maxmoves,
-			        maxsols, data_size, data, sols_size, sols, stats) :
-			    solve_h48(c, minmoves, maxmoves, maxsols,
+		parse_ret = parse_h48_solver(solver, &h, &k);
+		if (parse_ret == NISSY_OK)
+			return solve_h48(c, minmoves, maxmoves, maxsols,
 			        data_size, data, sols_size, sols, stats);
-		}
+		else
+			return parse_ret;
 	} else {
 		LOG("solve: unknown solver '%s'\n", solver);
 		return NISSY_ERROR_INVALID_SOLVER;
diff --git a/src/solvers/h48/coordinate_macros.h b/src/solvers/h48/coordinate_macros.h
index 5b67177..04462d6 100644
--- a/src/solvers/h48/coordinate_macros.h
+++ b/src/solvers/h48/coordinate_macros.h
@@ -2,5 +2,6 @@
 
 #define COCLASS_MASK (UINT32_C(0xFFFF) << UINT32_C(16))
 #define COCLASS(x)   (((x) & COCLASS_MASK) >> UINT32_C(16))
+#define ECLASS(x)    COCLASS(x)
 #define TTREP_MASK   (UINT32_C(0xFF) << UINT32_C(8))
 #define TTREP(x)     (((x) & TTREP_MASK) >> UINT32_C(8))
diff --git a/src/solvers/h48/gendata_cocsep.h b/src/solvers/h48/gendata_cocsep.h
index b46bcf3..291ae48 100644
--- a/src/solvers/h48/gendata_cocsep.h
+++ b/src/solvers/h48/gendata_cocsep.h
@@ -1,27 +1,16 @@
-STATIC_INLINE bool get_visited(const uint8_t *, int64_t);
-STATIC_INLINE void set_visited(uint8_t *, int64_t);
+STATIC_INLINE bool gendata_cocsep_get_visited(const uint8_t *, int64_t);
+STATIC_INLINE void gendata_cocsep_set_visited(uint8_t *, int64_t);
 
 STATIC size_t gendata_cocsep(
-    char [static COCSEP_FULLSIZE+INFOSIZE], uint64_t *, cube_t *);
+    char [static COCSEP_FULLSIZE], uint64_t *, cube_t *);
 STATIC uint32_t gendata_cocsep_dfs(cocsep_dfs_arg_t *);
 STATIC void getdistribution_cocsep(const uint32_t *, uint64_t [static 21]);
 
 STATIC_INLINE int8_t get_h48_cdata(cube_t, const uint32_t *, uint32_t *);
 
-/*
-Each element of the cocsep table is a uint32_t used as follows:
-  - Lowest 8-bit block: pruning value
-  - Second-lowest 8-bit block: "ttrep" (transformation to representative)
-  - Top 16-bit block: symcoord value
-After the data as described above, more auxiliary information is appended:
-  - A uint32_t representing the number of symmetry classes
-  - A uint32_t representing the highest value of the pruning table
-  - One uint32_t for each "line" of the pruning table, representing the number
-    of positions having that pruning value.
-*/
 STATIC size_t
 gendata_cocsep(
-	char buf[static COCSEP_FULLSIZE+INFOSIZE],
+	char buf[static COCSEP_FULLSIZE],
 	uint64_t *selfsim,
 	cube_t *rep
 )
@@ -36,7 +25,7 @@ gendata_cocsep(
 		goto gendata_cocsep_return_size;
 
 	memset(buf, 0xFF, COCSEP_FULLSIZE);
-	buf32 = (uint32_t *)((char *)buf + INFOSIZE);
+	buf32 = (uint32_t *)(buf + INFOSIZE);
 	if (selfsim != NULL)
 		memset(selfsim, 0, sizeof(uint64_t) * COCSEP_CLASSES);
 
@@ -69,7 +58,7 @@ gendata_cocsep(
 		info.distribution[i] = cc;
 	}
 
-	writetableinfo(&info, COCSEP_FULLSIZE+INFOSIZE, buf);
+	writetableinfo(&info, COCSEP_FULLSIZE, buf);
 
 	DBG_ASSERT(n == COCSEP_CLASSES, 0,
 	    "cocsep: computed %" PRIu16 " symmetry classes, "
@@ -101,9 +90,10 @@ gendata_cocsep_dfs(cocsep_dfs_arg_t *arg)
 
 	i = coord_cocsep(arg->cube);
 	olddepth = (uint8_t)(arg->buf32[i] & 0xFF);
-	if (olddepth < arg->depth || get_visited(arg->visited, i))
+	if (olddepth < arg->depth ||
+	    gendata_cocsep_get_visited(arg->visited, i))
 		return 0;
-	set_visited(arg->visited, i);
+	gendata_cocsep_set_visited(arg->visited, i);
 
 	if (arg->depth == arg->maxdepth) {
 		if ((arg->buf32[i] & 0xFF) != 0xFF)
@@ -118,7 +108,7 @@ gendata_cocsep_dfs(cocsep_dfs_arg_t *arg)
 				arg->selfsim[*arg->n] |= UINT64_C(1) << t;
 			if (COCLASS(arg->buf32[j]) != UINT32_C(0xFFFF))
 				continue;
-			set_visited(arg->visited, j);
+			gendata_cocsep_set_visited(arg->visited, j);
 			tinv = inverse_trans(t);
 			olddepth = arg->buf32[j] & 0xFF;
 			cc += olddepth == 0xFF;
@@ -155,13 +145,13 @@ getdistribution_cocsep(const uint32_t *table, uint64_t distr[static 21])
 }
 
 STATIC_INLINE bool
-get_visited(const uint8_t *a, int64_t i)
+gendata_cocsep_get_visited(const uint8_t *a, int64_t i)
 {
 	return a[VISITED_IND(i)] & VISITED_MASK(i);
 }
 
 STATIC_INLINE void
-set_visited(uint8_t *a, int64_t i)
+gendata_cocsep_set_visited(uint8_t *a, int64_t i)
 {
 	a[VISITED_IND(i)] |= VISITED_MASK(i);
 }
diff --git a/src/solvers/h48/gendata_eoesep.h b/src/solvers/h48/gendata_eoesep.h
new file mode 100644
index 0000000..d3b73ea
--- /dev/null
+++ b/src/solvers/h48/gendata_eoesep.h
@@ -0,0 +1,274 @@
+STATIC int64_t coord_eoesep_sym(cube_t, const uint32_t [static ESEP_MAX]);
+STATIC size_t gendata_esep_classes(
+    uint32_t [static ESEP_MAX], uint16_t [static ESEP_CLASSES]);
+STATIC size_t gendata_eoesep(char [static EOESEP_FULLSIZE], uint8_t);
+STATIC uint32_t gendata_eoesep_bfs(uint8_t, uint8_t [static EOESEP_BUF],
+    uint32_t [static ESEP_MAX], uint16_t [static ESEP_CLASSES]);
+STATIC uint32_t gendata_eoesep_fromnew(uint8_t, uint8_t [static EOESEP_BUF],
+    uint32_t [static ESEP_MAX], uint16_t [static ESEP_CLASSES]);
+STATIC uint32_t gendata_eoesep_fromdone(uint8_t, uint8_t [static EOESEP_BUF],
+    uint32_t [static ESEP_MAX], uint16_t [static ESEP_CLASSES]);
+STATIC uint32_t gendata_eoesep_marksim(int64_t, uint8_t,
+    uint8_t [static EOESEP_BUF], uint32_t [static ESEP_MAX]);
+STATIC bool gendata_eoesep_next(cube_t, uint8_t,
+    uint8_t [static EOESEP_BUF], uint32_t [static ESEP_MAX]);
+STATIC uint8_t get_eoesep_pval(const uint8_t *, int64_t);
+STATIC uint8_t get_eoesep_pval_cube(const void *, cube_t);
+STATIC void set_eoesep_pval(uint8_t *, int64_t, uint8_t);
+
+STATIC int64_t
+coord_eoesep_sym(cube_t c, const uint32_t esep_classes[static ESEP_MAX])
+{
+	uint8_t ttrep;
+	uint32_t edata, class;
+	int64_t esep, eo;
+
+	esep = coord_esep(c);
+	edata = esep_classes[esep];
+	class = ECLASS(edata);
+	ttrep = TTREP(edata);
+	eo = coord_eo(transform(c, ttrep));
+
+	return (class << UINT32_C(11)) + eo;
+}
+
+STATIC size_t
+gendata_esep_classes(
+	uint32_t esep_classes[static ESEP_MAX],
+	uint16_t rep[static ESEP_CLASSES]
+)
+{
+	bool visited[ESEP_MAX];
+	uint8_t t;
+	uint32_t class, cl, ti;
+	int64_t i, j;
+	cube_t c;
+
+	memset(visited, 0, ESEP_MAX * sizeof(bool));
+	class = 0;
+	for (i = 0; i < ESEP_MAX; i++) {
+		if (visited[i])
+			continue;
+		c = invcoord_esep(i);
+		for (t = 0; t < 48; t++) {
+			j = coord_esep(transform(c, t));
+			cl = class << UINT32_C(16);
+			ti = inverse_trans(t) << UINT32_C(8);
+			esep_classes[j] = cl | ti;
+			visited[j] = true;
+		}
+		rep[class] = i;
+		class++;
+	}
+
+	return class;
+}
+
+STATIC size_t
+gendata_eoesep(char buf[static EOESEP_FULLSIZE], uint8_t maxdepth)
+{
+	uint8_t *buf8, d;
+	uint16_t rep[ESEP_CLASSES];
+	uint32_t *esep_classes, done, level;
+	int64_t coord;
+	tableinfo_t info;
+
+	if (buf == NULL)
+		goto gendata_eoesep_return_size;
+
+	LOG("Computing eoesep data\n");
+	memset(buf, 0xFF, EOESEP_FULLSIZE);
+	esep_classes = (uint32_t *)(buf + INFOSIZE);
+	buf8 = (uint8_t *)(buf + INFOSIZE + 4*ESEP_MAX);
+	gendata_esep_classes(esep_classes, rep);
+
+	info = (tableinfo_t) {
+		.solver = "eoesep data for h48",
+		.type = TABLETYPE_SPECIAL,
+		.infosize = INFOSIZE,
+		.fullsize = EOESEP_FULLSIZE,
+		.hash = 0,
+		.entries = EOESEP_TABLESIZE,
+		.classes = ESEP_CLASSES,
+		.bits = 4,
+		.base = 0,
+		.maxvalue = 11,
+		.next = 0
+	};
+
+	coord = 0; /* Assumed coordinate of solved cube */
+	set_eoesep_pval(buf8, coord, 0);
+	done = 1;
+	info.distribution[0] = 1;
+	for (d = 1; d <= maxdepth && done < EOESEP_TABLESIZE; d++) {
+		level = gendata_eoesep_bfs(d, buf8, esep_classes, rep);
+		done += level;
+		info.distribution[d] = level;
+	}
+
+	writetableinfo(&info, EOESEP_FULLSIZE, buf);
+
+	LOG("eoesep data computed\n");
+
+gendata_eoesep_return_size:
+	return EOESEP_FULLSIZE;
+}
+
+STATIC uint32_t
+gendata_eoesep_bfs(
+	uint8_t d,
+	uint8_t buf8[EOESEP_BUF],
+	uint32_t esep_classes[static ESEP_MAX],
+	uint16_t rep[static ESEP_CLASSES]
+)
+{
+	if (d < 9)
+		return gendata_eoesep_fromdone(d, buf8, esep_classes, rep);
+	else
+		return gendata_eoesep_fromnew(d, buf8, esep_classes, rep);
+}
+
+STATIC uint32_t
+gendata_eoesep_fromdone(
+	uint8_t d,
+	uint8_t buf8[EOESEP_BUF],
+	uint32_t esep_classes[static ESEP_MAX],
+	uint16_t rep[static ESEP_CLASSES]
+)
+{
+	uint8_t pval;
+	int64_t i, esep, eo, coord, done;
+
+	done = 0;
+	for (i = 0; i < (int64_t)ESEP_CLASSES; i++) {
+		esep = rep[i];
+		for (eo = 0; eo < POW_2_11; eo++) {
+			coord = (i << INT64_C(11)) + eo;
+			pval = get_eoesep_pval(buf8, coord);
+			if (pval != d-1)
+				continue;
+
+			coord = (esep << INT64_C(11)) + eo;
+			done += gendata_eoesep_marksim(
+			    coord, d, buf8, esep_classes);
+		}
+	}
+
+	return done;
+}
+
+STATIC uint32_t
+gendata_eoesep_fromnew(
+	uint8_t d,
+	uint8_t buf8[EOESEP_BUF],
+	uint32_t esep_classes[static ESEP_MAX],
+	uint16_t rep[static ESEP_CLASSES]
+)
+{
+	uint8_t pval;
+	int64_t i, esep, eo, coord, done;
+	cube_t c;
+
+	done = 0;
+	for (i = 0; i < (int64_t)ESEP_CLASSES; i++) {
+		esep = rep[i];
+		for (eo = 0; eo < POW_2_11; eo++) {
+			coord = (i << INT64_C(11)) + eo;
+			pval = get_eoesep_pval(buf8, coord);
+			if (pval != 15)
+				continue;
+
+			c = invcoord_eoesep((esep << INT64_C(11)) + eo);
+			if (gendata_eoesep_next(c, d, buf8, esep_classes)) {
+				set_eoesep_pval(buf8, coord, d);
+				done++;
+			}
+		}
+	}
+
+	return done;
+}
+
+STATIC uint32_t
+gendata_eoesep_marksim(
+	int64_t i,
+	uint8_t d,
+	uint8_t buf8[static EOESEP_BUF],
+	uint32_t esep_classes[static ESEP_MAX]
+)
+{
+	uint8_t t, m, pval;
+	cube_t c, moved, transformed;
+	uint32_t done;
+	int64_t coord;
+
+	done = 0;
+	c = invcoord_eoesep(i);
+	for (m = 0; m < 18; m++) {
+		moved = move(c, m);
+		for (t = 0; t < 48; t++) {
+			transformed = transform(moved, t);
+			coord = coord_eoesep_sym(transformed, esep_classes);
+			pval = get_eoesep_pval(buf8, coord);
+			if (pval > d) {
+				set_eoesep_pval(buf8, coord, d);
+				done++;
+			}
+		}
+	}
+
+	return done;
+}
+
+STATIC bool
+gendata_eoesep_next(
+	cube_t c,
+	uint8_t d,
+	uint8_t buf8[static EOESEP_BUF],
+	uint32_t esep_classes[static ESEP_MAX]
+)
+{
+	uint8_t m, t, pval;
+	int64_t coord;
+	cube_t moved, transformed;
+
+	for (t = 0; t < 48; t++) {
+		transformed = transform(c, t);
+		for (m = 0; m < 18; m++) {
+			moved = move(transformed, m);
+			coord = coord_eoesep_sym(moved, esep_classes);
+			pval = get_eoesep_pval(buf8, coord);
+			if (pval == d-1)
+				return true;
+		}
+	}
+
+	return false;
+}
+
+STATIC uint8_t
+get_eoesep_pval(const uint8_t *table, int64_t i)
+{
+	return (table[EOESEP_INDEX(i)] & EOESEP_MASK(i)) >> EOESEP_SHIFT(i);
+}
+
+STATIC uint8_t
+get_eoesep_pval_cube(const void *data, cube_t c)
+{
+	int64_t coord;
+	const uint8_t *table;
+	const uint32_t *esep_classes;
+
+	esep_classes = (const uint32_t *)data;
+	table = (const uint8_t *)data + 4*ESEP_MAX;
+	coord = coord_eoesep_sym(c, esep_classes);
+
+	return get_eoesep_pval(table, coord);
+}
+
+STATIC void
+set_eoesep_pval(uint8_t *table, int64_t i, uint8_t val)
+{
+	table[EOESEP_INDEX(i)] = (table[EOESEP_INDEX(i)] & (~EOESEP_MASK(i)))
+	    | (val << EOESEP_SHIFT(i));
+}
diff --git a/src/solvers/h48/gendata_h48.h b/src/solvers/h48/gendata_h48.h
index 55c4f34..ae0b732 100644
--- a/src/solvers/h48/gendata_h48.h
+++ b/src/solvers/h48/gendata_h48.h
@@ -61,13 +61,13 @@ gendata_h48short(gendata_h48short_arg_t *arg)
 	return arg->map->n;
 }
 
-/* Generic function that dispatches to the data generators */
 STATIC int64_t
 gendata_h48(gendata_h48_arg_t *arg)
 {
-	uint64_t size, cocsepsize, h48size, fallbacksize;
+	uint64_t size, cocsepsize, h48size, fallbacksize, fallback2size, of;
+	long long r;
 	void *cocsepdata_offset;
-	tableinfo_t cocsepinfo, h48info;
+	tableinfo_t cocsepinfo, h48info, fallbackinfo;
 	gendata_h48_arg_t arg_h0k4;
 
 	if (arg == NULL) {
@@ -78,7 +78,8 @@ gendata_h48(gendata_h48_arg_t *arg)
 	cocsepsize = COCSEP_FULLSIZE;
 	h48size = INFOSIZE + H48_TABLESIZE(arg->h, arg->k);
 	fallbacksize = arg->k == 2 ? INFOSIZE + H48_TABLESIZE(0, 4) : 0;
-	size = cocsepsize + h48size + fallbacksize;
+	fallback2size = EOESEP_FULLSIZE;
+	size = cocsepsize + h48size + fallbacksize + fallback2size;
 
 	if (arg->buf == NULL)
 		return size; /* Dry-run */
@@ -110,18 +111,22 @@ gendata_h48(gendata_h48_arg_t *arg)
 		return NISSY_ERROR_INVALID_SOLVER;
 	}
 
-	if (readtableinfo(arg->buf_size, arg->buf, &cocsepinfo) != NISSY_OK) {
+	r = readtableinfo(arg->buf_size, arg->buf, &cocsepinfo);
+	if (r != NISSY_OK) {
 		LOG("gendata_h48: could not read info for cocsep table\n");
 		return NISSY_ERROR_UNKNOWN;
 	}
 
 	cocsepinfo.next = cocsepsize;
-	if (writetableinfo(&cocsepinfo, arg->buf_size, arg->buf) != NISSY_OK) {
+	r = writetableinfo(&cocsepinfo, arg->buf_size, arg->buf);
+	if (r != NISSY_OK) {
 		LOG("gendata_h48: could not write info for cocsep table"
 		    " with updated 'next' value\n");
 		return NISSY_ERROR_UNKNOWN;
 	}
 
+	/* Add h0k4 fallback table */
+
 	if (arg->k == 2) {
 		arg_h0k4 = *arg;
 		arg_h0k4.h = 0;
@@ -134,16 +139,42 @@ gendata_h48(gendata_h48_arg_t *arg)
 
 		gendata_h48h0k4(&arg_h0k4);
 
-		if (readtableinfo_n(arg->buf_size, arg->buf, 2, &h48info)
-		    != NISSY_OK) {
-			LOG("gendata_h48: could not read info for h48 table\n");
+	}
+
+	/* Add eoesep fallback table */
+
+	gendata_eoesep((char *)arg->buf + (size - fallback2size), 20);
+
+	/* Update tableinfo with correct next values */
+
+	r = readtableinfo_n(arg->buf_size, arg->buf, 2, &h48info);
+	if (r != NISSY_OK) {
+		LOG("gendata_h48: could not read info for h48 table\n");
+		return NISSY_ERROR_UNKNOWN;
+	}
+	h48info.next = h48size;
+	r = writetableinfo(&h48info,
+	     arg->buf_size - cocsepsize, (char *)arg->buf + cocsepsize);
+	if (r != NISSY_OK) {
+		LOG("gendata_h48: could not write info for h48 table\n");
+		return NISSY_ERROR_UNKNOWN;
+	}
+
+	if (arg->k == 2) {
+		r = readtableinfo_n(arg->buf_size, arg->buf, 3, &fallbackinfo);
+		if (r != NISSY_OK) {
+			LOG("gendata_h48: could not read info for h48 "
+			    "fallback table\n");
 			return NISSY_ERROR_UNKNOWN;
 		}
 
-		h48info.next = h48size;
-		if (writetableinfo(&h48info, arg->buf_size - cocsepsize,
-		    (char *)arg->buf + cocsepsize) != NISSY_OK) {
-			LOG("gendata_h48: could not write info for h48 table\n");
+		of = cocsepsize + h48size;
+		fallbackinfo.next = fallbacksize;
+		r = writetableinfo(&fallbackinfo,
+		    arg->buf_size - of, (char *)arg->buf + of);
+		if (r != NISSY_OK) {
+			LOG("gendata_h48: could not write info for h48 "
+			    "fallback table\n");
 			return NISSY_ERROR_UNKNOWN;
 		}
 	}
@@ -302,9 +333,8 @@ gendata_h48k2(gendata_h48_arg_t *arg)
 	 *
 	 * The following values for the base have been hand-picked. I first
 	 * performed some statistics on the frequency of these values, but
-	 * they turned out to be unreliable. I have not figured out why yet.
-	 * In the end I resorted to generating the same table with multiple
-	 * base value and see what was best.
+	 * they turned out to be unreliable. In the end I generated the same
+	 * table with multiple base value and see what was best.
 	 *
 	 * A curious case is h3, which has this distribution for base 8:
 	 *   [0] = 6686828
diff --git a/src/solvers/h48/gendata_types_macros.h b/src/solvers/h48/gendata_types_macros.h
index 73a46ff..d7c88af 100644
--- a/src/solvers/h48/gendata_types_macros.h
+++ b/src/solvers/h48/gendata_types_macros.h
@@ -3,13 +3,22 @@
 #define COCSEP_VISITEDSIZE  DIV_ROUND_UP(COCSEP_TABLESIZE, (size_t)8)
 #define COCSEP_FULLSIZE     (INFOSIZE + (size_t)4 * COCSEP_TABLESIZE)
 
+#define ESEP_MAX            (COMB_12_4 * COMB_8_4)
+#define ESEP_CLASSES        ((size_t)782)
+#define EOESEP_TABLESIZE    (ESEP_CLASSES << (size_t)11)
+#define EOESEP_BUF          DIV_ROUND_UP(EOESEP_TABLESIZE, 2)
+#define EOESEP_FULLSIZE     (INFOSIZE + EOESEP_BUF + (size_t)4 * ESEP_MAX)
+#define EOESEP_INDEX(i)     ((i)/2)
+#define EOESEP_SHIFT(i)     (UINT8_C(4) * (uint8_t)((i) % 2))
+#define EOESEP_MASK(i)      (UINT8_C(0xF) << EOESEP_SHIFT(i))
+
 #define VISITED_IND(i)      ((uint32_t)(i) / UINT32_C(8))
 #define VISITED_MASK(i)     (UINT32_C(1) << ((uint32_t)(i) % UINT32_C(8)))
 
 #define CBOUND_MASK         UINT32_C(0xFF)
 #define CBOUND(x)           ((x) & CBOUND_MASK)
 
-#define H48_COORDMAX_NOEO   ((int64_t)(COCSEP_CLASSES * COMB_12_4 * COMB_8_4))
+#define H48_COORDMAX_NOEO   ((int64_t)(COCSEP_CLASSES * ESEP_MAX))
 #define H48_COORDMAX(h)     (H48_COORDMAX_NOEO << (int64_t)(h))
 #define H48_DIV(k)          ((size_t)8 / (size_t)(k))
 #define H48_TABLESIZE(h, k) DIV_ROUND_UP((size_t)H48_COORDMAX((h)), H48_DIV(k))
diff --git a/src/solvers/h48/h48.h b/src/solvers/h48/h48.h
index 3cf02b7..d5a67ef 100644
--- a/src/solvers/h48/h48.h
+++ b/src/solvers/h48/h48.h
@@ -2,7 +2,6 @@
 #include "map.h"
 #include "gendata_types_macros.h"
 #include "gendata_cocsep.h"
+#include "gendata_eoesep.h"
 #include "gendata_h48.h"
-#include "stats.h"
 #include "solve.h"
-#include "solve_multithread.h"
diff --git a/src/solvers/h48/solve.h b/src/solvers/h48/solve.h
index d506596..8666007 100644
--- a/src/solvers/h48/solve.h
+++ b/src/solvers/h48/solve.h
@@ -1,9 +1,26 @@
+#define STARTING_MOVES 3
+#define STARTING_CUBES 3240 /* Number of 3-move sequences */
+
+typedef struct {
+	cube_t cube;
+	uint8_t moves[STARTING_MOVES];
+	uint64_t symmask0;
+} solve_h48_task_t;
+
 typedef struct {
+	cube_t start_cube;
+	uint64_t symmask0;
 	cube_t cube;
 	cube_t inverse;
-	int8_t nmoves;
 	int8_t depth;
+	int8_t nmoves;
 	uint8_t moves[MAXLEN];
+	int8_t npremoves;
+	uint8_t premoves[MAXLEN];
+	int8_t lb_normal;
+	int8_t lb_inverse;
+	bool use_lb_normal;
+	bool use_lb_inverse;
 	_Atomic int64_t *nsols;
 	int64_t maxsolutions;
 	uint8_t h;
@@ -11,209 +28,385 @@ typedef struct {
 	uint8_t base;
 	const uint32_t *cocsepdata;
 	const uint8_t *h48data;
-	const uint8_t *h48data_fallback;
+	const uint8_t *h48data_fallback_h0k4;
+	const void *h48data_fallback_eoesep;
 	uint64_t solutions_size;
-	char **nextsol;
-	uint8_t nissbranch;
-	int8_t npremoves;
-	uint8_t premoves[MAXLEN];
-	long long nodes_visited;
-	long long table_fallbacks;
-} dfsarg_solveh48_t;
+	uint64_t *solutions_used;
+	char **solutions;
+	uint32_t movemask_normal;
+	uint32_t movemask_inverse;
+	int64_t nodes_visited;
+	int64_t table_fallbacks;
+	int64_t table_lookups;
+	int ntasks;
+	solve_h48_task_t *tasks;
+	int thread_id;
+	pthread_mutex_t *solutions_mutex;
+} dfsarg_solve_h48_t;
 
-STATIC uint32_t allowednextmove_h48(uint8_t *, uint8_t, uint8_t);
-
-STATIC void solve_h48_appendsolution(dfsarg_solveh48_t *);
-STATIC_INLINE bool solve_h48_stop(dfsarg_solveh48_t *);
-STATIC int64_t solve_h48_dfs(dfsarg_solveh48_t *);
-STATIC int64_t solve_h48(cube_t, int8_t, int8_t, int8_t, uint64_t,
+typedef struct {
+	cube_t cube;
+	int8_t nmoves;
+	uint8_t moves[STARTING_MOVES];
+	int8_t minmoves;
+	int8_t maxmoves;
+} dfsarg_solve_h48_maketasks_t;
+
+STATIC int64_t solve_h48_appendsolution(dfsarg_solve_h48_t *);
+STATIC bool solve_h48_appendmoves(dfsarg_solve_h48_t *, int8_t,
+    uint8_t *, uint8_t);
+STATIC bool solve_h48_appendchar(dfsarg_solve_h48_t *, char);
+STATIC_INLINE bool solve_h48_stop(dfsarg_solve_h48_t *);
+STATIC int64_t solve_h48_maketasks(
+    dfsarg_solve_h48_t *, dfsarg_solve_h48_maketasks_t *,
+    solve_h48_task_t [static STARTING_CUBES], int *);
+STATIC void *solve_h48_runthread(void *);
+STATIC int64_t solve_h48_dfs(dfsarg_solve_h48_t *);
+STATIC int64_t solve_h48(cube_t, int8_t, int8_t, uint64_t, uint64_t,
     const void *, uint64_t, char *, long long [static NISSY_SIZE_SOLVE_STATS]);
 
-STATIC uint32_t
-allowednextmove_h48(uint8_t *moves, uint8_t n, uint8_t h48branch)
+STATIC int64_t
+solve_h48_appendsolution(dfsarg_solve_h48_t *arg)
 {
-	uint32_t result = MM_ALLMOVES;
-	if (h48branch & MM_NORMALBRANCH)
-		result &= MM_NOHALFTURNS;
-	if (n < 1)
-		return result;
+	uint8_t t;
+	int64_t ret;
+	uint64_t solstart;
+
+	if (*arg->nsols >= arg->maxsolutions)
+		return 0;
+
+	solstart = *arg->solutions_used;
+	invertmoves(arg->premoves, arg->npremoves, arg->moves + arg->nmoves);
 
-	uint8_t base1 = movebase(moves[n-1]);
-	uint8_t axis1 = moveaxis(moves[n-1]);
+	/* Do not append the solution in case premoves cancel with normal */
+	if (arg->npremoves > 0 && !allowednextmove(arg->moves, arg->nmoves+1))
+		return 0;
+	if (arg->npremoves > 1 && !allowednextmove(arg->moves, arg->nmoves+2))
+		return 0;
+
+	for (t = 0, ret = 0; t < 48 && *arg->nsols < arg->maxsolutions; t++) {
+		if (!(arg->symmask0 & (UINT64_C(1) << (uint64_t)t)))
+			continue;
 
-	result = disable_moves(result, base1 * 3);
-	if (base1 % 2)
-		result = disable_moves(result, (base1 - 1) * 3);
+		if (!solve_h48_appendmoves(arg, arg->nmoves + arg->npremoves,
+		    arg->moves, t))
+			goto solve_h48_appendsolution_error;
 
-	if (n == 1)
-		return result;
+		LOG("Solution found: %s\n", *arg->solutions + solstart);
 
-	uint8_t base2 = movebase(moves[n-2]);
-	uint8_t axis2 = moveaxis(moves[n-2]);
+		if (!solve_h48_appendchar(arg, '\n'))
+			goto solve_h48_appendsolution_error;
+		(*arg->nsols)++;
+		ret++;
+	}
 
-	if(axis1 == axis2)
-		result = disable_moves(result, base2 * 3);
+	return ret;
 
-	return result;
+solve_h48_appendsolution_error:
+	LOG("Could not append solution to buffer: size too small\n");
+	return NISSY_ERROR_BUFFER_SIZE;
 }
 
-STATIC void
-solve_h48_appendsolution(dfsarg_solveh48_t *arg)
+STATIC bool
+solve_h48_appendmoves(
+	dfsarg_solve_h48_t *arg,
+	int8_t n,
+	uint8_t *moves,
+	uint8_t t
+)
 {
+	int i;
 	int64_t strl;
-	uint8_t invertedpremoves[MAXLEN];
-	char *solution = *arg->nextsol; 
+	uint8_t mm[MAXLEN];
 
-	strl = writemoves(
-	    arg->moves, arg->nmoves, arg->solutions_size, *arg->nextsol);
+	for (i = 0; i < n; i++)
+		mm[i] = transform_move(moves[i], t);
 
-	if (strl < 0)
-		goto solve_h48_appendsolution_error;
-	*arg->nextsol += strl-1; 
-	arg->solutions_size -= strl-1;
+	strl = writemoves(mm, n, arg->solutions_size - *arg->solutions_used,
+	    *arg->solutions + *arg->solutions_used);
 
-	if (arg->npremoves) {
-		**arg->nextsol = ' ';
-		(*arg->nextsol)++;
-		arg->solutions_size--;
+	if (strl < 0)
+		return false;
 
-		invertmoves(arg->premoves, arg->npremoves, invertedpremoves);
-		strl = writemoves(invertedpremoves,
-		    arg->npremoves, arg->solutions_size, *arg->nextsol);
+	*arg->solutions_used += MAX(0, strl-1);
+	return true;
+}
 
-		if (strl < 0)
-			goto solve_h48_appendsolution_error;
-		*arg->nextsol += strl-1;
-		arg->solutions_size -= strl-1;
-	}
-	LOG("Solution found: %s\n", solution);
+STATIC bool
+solve_h48_appendchar(dfsarg_solve_h48_t *arg, char c)
+{
+	if (arg->solutions_size <= *arg->solutions_used)
+		return false;
 
-	**arg->nextsol = '\n';
-	(*arg->nextsol)++;
-	arg->solutions_size--;
-	(*arg->nsols)++;
+	*(*arg->solutions + *arg->solutions_used) = c;
+	(*arg->solutions_used)++;
 
-solve_h48_appendsolution_error:
-	/* We could add some logging, but writemoves() already does */
-	return;
+	return true;
 }
 
 STATIC_INLINE bool
-solve_h48_stop(dfsarg_solveh48_t *arg)
+solve_h48_stop(dfsarg_solve_h48_t *arg)
 {
 	uint32_t data, data_inv;
-	int8_t cbound, cbound_inv, h48bound, h48bound_inv;
-	int64_t coord, coord_inv;
+	int64_t coord;
+	int8_t target, nh;
+	uint8_t pval_cocsep, pval_eoesep;
 
+	target = arg->depth - arg->nmoves - arg->npremoves;
+	if (target <= 0 || *arg->nsols == arg->maxsolutions)
+		return true;
+
+	arg->movemask_normal = arg->movemask_inverse = MM_ALLMOVES;
 	arg->nodes_visited++;
 
-	arg->nissbranch = MM_NORMAL;
-	cbound = get_h48_cdata(arg->cube, arg->cocsepdata, &data);
-	if (cbound + arg->nmoves + arg->npremoves > arg->depth)
+	/* Preliminary probing using last computed bound, if possible */
+
+	if ((arg->use_lb_normal && arg->lb_normal > target) ||
+	    (arg->use_lb_inverse && arg->lb_inverse > target))
 		return true;
 
-	cbound_inv = get_h48_cdata(arg->inverse, arg->cocsepdata, &data_inv);
-	if (cbound_inv + arg->nmoves + arg->npremoves > arg->depth)
+	/* Preliminary corner probing */
+
+	if (get_h48_cdata(arg->cube, arg->cocsepdata, &data) > target ||
+	    get_h48_cdata(arg->inverse, arg->cocsepdata, &data_inv) > target)
 		return true;
 
-	coord = coord_h48_edges(arg->cube, COCLASS(data), TTREP(data), arg->h);
-	h48bound = get_h48_pval(arg->h48data, coord, arg->k);
+	/* Inverse probing */
 
-	/* If the h48 bound is > 0, we add the base value.    */
-	/* Otherwise, we use the fallback h0k4 value instead. */
+	if (!arg->use_lb_inverse) {
+		coord = coord_h48_edges(
+		    arg->inverse, COCLASS(data_inv), TTREP(data_inv), arg->h);
+		arg->lb_inverse = get_h48_pval(arg->h48data, coord, arg->k);
+		arg->table_lookups++;
 
-	if (arg->k == 2) {
-		if (h48bound == 0) {
+		if (arg->k == 2 && arg->lb_inverse == 0) {
 			arg->table_fallbacks++;
-			h48bound = get_h48_pval(
-			    arg->h48data_fallback, coord >> arg->h, 4);
+
+			pval_cocsep = get_h48_pval(
+			    arg->h48data_fallback_h0k4, coord >> arg->h, 4);
+			pval_eoesep = get_eoesep_pval_cube(
+			    arg->h48data_fallback_eoesep, arg->inverse);
+			arg->lb_inverse = MAX(pval_cocsep, pval_eoesep);
 		} else {
-			h48bound += arg->base;
+			arg->lb_inverse += arg->base;
 		}
+
+		arg->use_lb_inverse = true;
 	}
-	if (h48bound + arg->nmoves + arg->npremoves > arg->depth)
+
+	if (arg->lb_inverse > target)
 		return true;
-	if (h48bound + arg->nmoves + arg->npremoves == arg->depth)
-		arg->nissbranch = MM_INVERSEBRANCH;
-
-	coord_inv = coord_h48_edges(
-	    arg->inverse, COCLASS(data_inv), TTREP(data_inv), arg->h);
-	h48bound_inv = get_h48_pval(arg->h48data, coord_inv, arg->k);
-	if (arg->k == 2) {
-		if (h48bound_inv == 0) {
+	nh = arg->lb_inverse == target;
+	arg->movemask_normal = nh * MM_NOHALFTURNS + (1-nh) * MM_ALLMOVES;
+
+	/* Normal probing */
+
+	if (!arg->use_lb_normal) {
+		coord = coord_h48_edges(
+		    arg->cube, COCLASS(data), TTREP(data), arg->h);
+		arg->lb_normal = get_h48_pval(arg->h48data, coord, arg->k);
+		arg->table_lookups++;
+
+		if (arg->k == 2 && arg->lb_normal == 0) {
 			arg->table_fallbacks++;
-			h48bound_inv = get_h48_pval(
-			    arg->h48data_fallback, coord_inv >> arg->h, 4);
+
+			pval_cocsep = get_h48_pval(
+			    arg->h48data_fallback_h0k4, coord >> arg->h, 4);
+			pval_eoesep = get_eoesep_pval_cube(
+			    arg->h48data_fallback_eoesep, arg->cube);
+			arg->lb_normal = MAX(pval_cocsep, pval_eoesep);
 		} else {
-			h48bound_inv += arg->base;
+			arg->lb_normal += arg->base;
 		}
+
+		arg->use_lb_normal = true;
 	}
-	if (h48bound_inv + arg->nmoves + arg->npremoves > arg->depth)
+
+	if (arg->lb_normal > target)
 		return true;
-	if (h48bound_inv + arg->nmoves + arg->npremoves == arg->depth)
-		arg->nissbranch = MM_NORMALBRANCH;
+	nh = arg->lb_normal == target;
+	arg->movemask_inverse = nh * MM_NOHALFTURNS + (1-nh) * MM_ALLMOVES;
 
 	return false;
 }
 
 STATIC int64_t
-solve_h48_dfs(dfsarg_solveh48_t *arg)
+solve_h48_dfs(dfsarg_solve_h48_t *arg)
 {
-	dfsarg_solveh48_t nextarg;
-	int64_t ret;
-	uint8_t m;
-
-	if (*arg->nsols == arg->maxsolutions)
-		return 0;
-
-	if (solve_h48_stop(arg))
-		return 0;
+	int64_t ret, n;
+	uint8_t m, lbn, lbi;
+	uint32_t mm_normal, mm_inverse;
+	bool ulbi, ulbn;
+	cube_t backup_cube, backup_inverse;
 
 	if (issolved(arg->cube)) {
 		if (arg->nmoves + arg->npremoves != arg->depth)
 			return 0;
-		solve_h48_appendsolution(arg);
-		return 1;
+		pthread_mutex_lock(arg->solutions_mutex);
+		ret = solve_h48_appendsolution(arg);
+		pthread_mutex_unlock(arg->solutions_mutex);
+		return ret;
 	}
 
-	nextarg = *arg;
+	if (solve_h48_stop(arg))
+		return 0;
+
+	backup_cube = arg->cube;
+	backup_inverse = arg->inverse;
+	lbn = arg->lb_normal;
+	lbi = arg->lb_inverse;
+	ulbn = arg->use_lb_normal;
+	ulbi = arg->use_lb_inverse;
+
 	ret = 0;
-	uint32_t allowed;
-	if(arg->nissbranch & MM_INVERSE) {
-		allowed = allowednextmove_h48(arg->premoves, arg->npremoves, arg->nissbranch);
+	mm_normal = allowednextmove_mask(arg->moves, arg->nmoves) &
+	    arg->movemask_normal;
+	mm_inverse = allowednextmove_mask(arg->premoves, arg->npremoves) &
+	    arg->movemask_inverse;
+	if (popcount_u32(mm_normal) <= popcount_u32(mm_inverse)) {
+		arg->nmoves++;
 		for (m = 0; m < 18; m++) {
-			if(allowed & (1 << m)) {
-				nextarg.npremoves = arg->npremoves + 1;
-				nextarg.premoves[arg->npremoves] = m;
-				nextarg.inverse = move(arg->inverse, m);
-				nextarg.cube = premove(arg->cube, m);
-				ret += solve_h48_dfs(&nextarg);
-			}
+			if (!(mm_normal & (1 << m)))
+				continue;
+			arg->moves[arg->nmoves-1] = m;
+			arg->cube = move(backup_cube, m);
+			arg->inverse = premove(backup_inverse, m);
+			arg->lb_inverse = lbi;
+			arg->use_lb_normal = false;
+			arg->use_lb_inverse = ulbi && m % 3 == 1;
+			n = solve_h48_dfs(arg);
+			if (n < 0)
+				return n;
+			ret += n;
 		}
+		arg->nmoves--;
 	} else {
-		allowed = allowednextmove_h48(arg->moves, arg->nmoves, arg->nissbranch);
+		arg->npremoves++;
 		for (m = 0; m < 18; m++) {
-			if (allowed & (1 << m)) {
-				nextarg.nmoves = arg->nmoves + 1;
-				nextarg.moves[arg->nmoves] = m;
-				nextarg.cube = move(arg->cube, m);
-				nextarg.inverse = premove(arg->inverse, m);
-				ret += solve_h48_dfs(&nextarg);
-			}
+			if(!(mm_inverse & (1 << m)))
+				continue;
+			arg->premoves[arg->npremoves-1] = m;
+			arg->inverse = move(backup_inverse, m);
+			arg->cube = premove(backup_cube, m);
+			arg->lb_normal = lbn;
+			arg->use_lb_inverse = false;
+			arg->use_lb_normal = ulbn && m % 3 == 1;
+			n = solve_h48_dfs(arg);
+			if (n < 0)
+				return n;
+			ret += n;
 		}
+		arg->npremoves--;
 	}
 
-	arg->nodes_visited = nextarg.nodes_visited;
-	arg->table_fallbacks = nextarg.table_fallbacks;
+	arg->cube = backup_cube;
+	arg->inverse = backup_inverse;
+
 	return ret;
 }
 
+STATIC void *
+solve_h48_runthread(void *arg)
+{
+	int i, j;
+	solve_h48_task_t task;
+	dfsarg_solve_h48_t * dfsarg;
+	cube_t cube;
+
+	dfsarg = (dfsarg_solve_h48_t *)arg;
+	cube = dfsarg->start_cube;
+
+	for (i = dfsarg->thread_id; i < dfsarg->ntasks; i += THREADS) {
+		task = dfsarg->tasks[i];
+		memcpy(dfsarg->moves, task.moves, STARTING_MOVES);
+		dfsarg->cube = cube;
+		for (j = 0; j < STARTING_MOVES; j++)
+			dfsarg->cube = move(
+			    dfsarg->cube, dfsarg->moves[j]);
+		dfsarg->inverse = inverse(dfsarg->cube);
+		dfsarg->nmoves = STARTING_MOVES;
+		dfsarg->npremoves = 0;
+		dfsarg->lb_normal = 0;
+		dfsarg->lb_inverse = 0;
+		dfsarg->use_lb_normal = false;
+		dfsarg->use_lb_inverse = false;
+		dfsarg->movemask_normal = MM_ALLMOVES;
+		dfsarg->movemask_inverse = MM_ALLMOVES;
+
+		solve_h48_dfs(dfsarg);
+	}
+
+	return NULL;
+}
+
+STATIC int64_t
+solve_h48_maketasks(
+	dfsarg_solve_h48_t *solve_arg,
+	dfsarg_solve_h48_maketasks_t *maketasks_arg,
+	solve_h48_task_t tasks[static STARTING_CUBES],
+	int *ntasks
+)
+{
+	int r;
+	int64_t appret;
+	uint8_t m, t;
+	uint32_t mm;
+	cube_t backup_cube;
+
+	if (issolved(maketasks_arg->cube)) {
+		if (maketasks_arg->nmoves > maketasks_arg->maxmoves ||
+		    maketasks_arg->nmoves < maketasks_arg->minmoves ||
+		    *solve_arg->nsols >= solve_arg->maxsolutions)
+			return NISSY_OK;
+		memcpy(solve_arg->moves,
+		    maketasks_arg->moves, maketasks_arg->nmoves);
+		solve_arg->nmoves = maketasks_arg->nmoves;
+		appret = solve_h48_appendsolution(solve_arg);
+		return appret < 0 ? appret : NISSY_OK;
+	}
+
+	if (maketasks_arg->nmoves == STARTING_MOVES) {
+		tasks[*ntasks].cube = maketasks_arg->cube;
+		memcpy(tasks[*ntasks].moves,
+		    maketasks_arg->moves, STARTING_MOVES);
+		(*ntasks)++;
+		return NISSY_OK;
+	}
+
+	mm = allowednextmove_mask(maketasks_arg->moves, maketasks_arg->nmoves);
+
+	maketasks_arg->nmoves++;
+	backup_cube = maketasks_arg->cube;
+	for (m = 0; m < 18; m++) {
+		if (!(mm & (1 << m)))
+			continue;
+		maketasks_arg->moves[maketasks_arg->nmoves-1] = m;
+		maketasks_arg->cube = move(backup_cube, m);
+		r = solve_h48_maketasks(
+		    solve_arg, maketasks_arg, tasks, ntasks);
+		if (r < 0)
+			return r;
+
+		/* Avoid symmetry-equivalent moves from the starting cube */
+		if (maketasks_arg->nmoves == 1)
+			for (t = 0; t < 48; t++)
+				if (solve_arg->symmask0 &
+				    (UINT64_C(1) << (uint64_t)t))
+					mm &= ~(UINT32_C(1) <<
+					    (uint32_t)transform_move(m, t));
+	}
+	maketasks_arg->nmoves--;
+	maketasks_arg->cube = backup_cube;
+
+	return NISSY_OK;
+}
+
 STATIC int64_t
 solve_h48(
 	cube_t cube,
 	int8_t minmoves,
 	int8_t maxmoves,
-	int8_t maxsolutions,
+	uint64_t maxsolutions,
 	uint64_t data_size,
 	const void *data,
 	uint64_t solutions_size,
@@ -221,61 +414,150 @@ solve_h48(
 	long long stats[static NISSY_SIZE_SOLVE_STATS]
 )
 {
+	int i, ntasks, eoesep_table_index;
+	int8_t d;
 	_Atomic int64_t nsols;
-	dfsarg_solveh48_t arg;
-	tableinfo_t info, fbinfo;
+	dfsarg_solve_h48_t arg[THREADS];
+	solve_h48_task_t tasks[STARTING_CUBES];
+	dfsarg_solve_h48_maketasks_t maketasks_arg;
+	long double fallback_rate, lookups_per_node;
+	uint64_t solutions_used, symmask, offset;
+	int64_t nodes_visited, table_lookups, table_fallbacks;
+	tableinfo_t info, fbinfo, fbinfo2;
+	const uint32_t *cocsepdata;
+	const uint8_t *fallback, *h48data;
+	const void *fallback2;
+	pthread_t thread[THREADS];
+	pthread_mutex_t solutions_mutex;
 
 	if(readtableinfo_n(data_size, data, 2, &info) != NISSY_OK)
 		goto solve_h48_error_data;
 
-	arg = (dfsarg_solveh48_t) {
-		.cube = cube,
-		.inverse = inverse(cube),
-		.nsols = &nsols,
-		.maxsolutions = maxsolutions,
-		.h = info.h48h,
-		.k = info.bits,
-		.base = info.base,
-		.cocsepdata = (uint32_t *)((char *)data + INFOSIZE),
-		.h48data = (uint8_t *)data + COCSEP_FULLSIZE + INFOSIZE,
-		.solutions_size = solutions_size,
-		.nextsol = &solutions,
-		.nodes_visited = 0,
-		.table_fallbacks = 0
-	};
+	cocsepdata = (uint32_t *)((char *)data + INFOSIZE);
+	h48data = (uint8_t *)data + COCSEP_FULLSIZE + INFOSIZE;
 
+	/* Read fallback table(s) */
+	fallback = NULL;
+	if (readtableinfo_n(data_size, data, 3, &fbinfo) != NISSY_OK)
+		goto solve_h48_error_data;
+	offset = info.next;
+	eoesep_table_index = 3;
 	if (info.bits == 2) {
-		if (readtableinfo_n(data_size, data, 3, &fbinfo) != NISSY_OK)
-			goto solve_h48_error_data;
 		/* We only support h0k4 as fallback table */
 		if (fbinfo.h48h != 0 || fbinfo.bits != 4)
 			goto solve_h48_error_data;
-		arg.h48data_fallback = arg.h48data + info.next;
-	} else {
-		arg.h48data_fallback = NULL;
+		fallback = h48data + offset;
+		offset += fbinfo.next;
+		eoesep_table_index++;
+	}
+
+	if (readtableinfo_n(data_size, data, eoesep_table_index, &fbinfo2)
+	    != NISSY_OK)
+		goto solve_h48_error_data;
+
+	/* Some heuristic check to see that it is eoesep */
+	if (fbinfo2.bits != 4 || fbinfo2.type != TABLETYPE_SPECIAL)
+		goto solve_h48_error_data;
+	fallback2 = h48data + offset;
+
+	symmask = symmetry_mask(cube);
+	for (i = 0; i < THREADS; i++) {
+		arg[i] = (dfsarg_solve_h48_t) {
+			.start_cube = cube,
+			.cube = cube,
+			.symmask0 = symmask,
+			.nsols = &nsols,
+			.maxsolutions = maxsolutions,
+			.h = info.h48h,
+			.k = info.bits,
+			.base = info.base,
+			.cocsepdata = cocsepdata,
+			.h48data = h48data,
+			.h48data_fallback_h0k4 = fallback,
+			.h48data_fallback_eoesep = fallback2,
+			.solutions_size = solutions_size,
+			.solutions_used = &solutions_used,
+			.solutions = &solutions,
+			.nodes_visited = 0,
+			.table_fallbacks = 0,
+			.table_lookups = 0,
+			.thread_id = i,
+			.solutions_mutex = &solutions_mutex,
+		};
+
 	}
 
 	nsols = 0;
-	for (arg.depth = minmoves;
-	     arg.depth <= maxmoves && nsols < maxsolutions;
-	     arg.depth++)
-	{
-		LOG("Found %" PRId64 " solutions, searching at depth %"
-		    PRId8 "\n", nsols, arg.depth);
-		arg.nmoves = 0;
-		arg.npremoves = 0;
-		solve_h48_dfs(&arg);
+	solutions_used = 0;
+
+	pthread_mutex_init(&solutions_mutex, NULL);
+
+	maketasks_arg = (dfsarg_solve_h48_maketasks_t) {
+		.cube = cube,
+		.nmoves = 0,
+		.minmoves = minmoves,
+		.maxmoves = maxmoves,
+	};
+	ntasks = 0;
+	solve_h48_maketasks(&arg[0], &maketasks_arg, tasks, &ntasks);
+	if (ntasks < 0)
+		goto solve_h48_error_solutions_buffer;
+	if (*arg[0].nsols >= (int64_t)maxsolutions)
+		goto solve_h48_done;
+
+	for (i = 0; i < THREADS; i++) {
+		arg[i].ntasks = ntasks;
+		arg[i].tasks = tasks;
+	}
+
+	LOG("Prepared %d tasks\n", ntasks);
+
+	for (
+	    d = MAX(minmoves, STARTING_MOVES + 1);
+	    d <= maxmoves && nsols < (int64_t)maxsolutions;
+	    d++
+	) {
+		if (d >= 10)
+			LOG("Found %" PRId64 " solutions, searching at depth %"
+			    PRId8 "\n", nsols, d);
+		for (i = 0; i < THREADS; i++) {
+			arg[i].depth = d;
+			pthread_create(
+			    &thread[i], NULL, solve_h48_runthread, &arg[i]);
+		}
+		for (i = 0; i < THREADS; i++)
+			pthread_join(thread[i], NULL);
+	}
+
+solve_h48_done:
+	if (!solve_h48_appendchar(&arg[0], '\0'))
+		goto solve_h48_error_solutions_buffer;
+
+	nodes_visited = table_lookups = table_fallbacks = 0;
+	for (i = 0; i < THREADS; i++) {
+		nodes_visited += arg[i].nodes_visited;
+		table_fallbacks += arg[i].table_fallbacks;
+		table_lookups += arg[i].table_lookups;
 	}
-	**arg.nextsol = '\0';
 
-	stats[0] = arg.nodes_visited;
-	stats[1] = arg.table_fallbacks;
-	LOG("Nodes visited: %lld\nTable fallbacks: %lld\n",
-	    arg.nodes_visited, arg.table_fallbacks);
+	stats[0] = nodes_visited;
+	stats[1] = table_lookups;
+	stats[2] = table_fallbacks;
+	lookups_per_node = table_lookups / (long double)nodes_visited;
+	fallback_rate = nodes_visited == 0 ? 0.0 :
+	    (table_fallbacks * 100) / (long double)table_lookups;
+	LOG("Nodes visited: %" PRId64 "\n", nodes_visited);
+	LOG("Lookups: %" PRId64 " (%.3Lf per node)\n",
+	    table_lookups, lookups_per_node);
+	LOG("Table fallbacks: %" PRId64 " (%.3Lf%%)\n",
+	    table_fallbacks, fallback_rate);
 
 	return nsols;
 
 solve_h48_error_data:
 	LOG("solve_h48: error reading table\n");
 	return NISSY_ERROR_DATA;
+
+solve_h48_error_solutions_buffer:
+	return NISSY_ERROR_BUFFER_SIZE;
 }
diff --git a/src/solvers/h48/solve_multithread.h b/src/solvers/h48/solve_multithread.h
deleted file mode 100644
index 72f519e..0000000
--- a/src/solvers/h48/solve_multithread.h
+++ /dev/null
@@ -1,350 +0,0 @@
-#define MAX_QUEUE_SIZE 244
-#define BFS_DEPTH 2
-
-typedef struct {
-	dfsarg_solveh48_t tasks[MAX_QUEUE_SIZE];
-	int front;
-	int rear;
-	int tasks_count;
-	int active;
-	pthread_mutex_t mutex;
-	pthread_cond_t cond;
-	pthread_cond_t active_cond;
-	atomic_bool terminate;
-	_Atomic long long nodes_visited_global;
-	_Atomic long long table_fallbacks_global;
-} task_queue_t;
-
-STATIC void solve_h48_appendsolution_thread(dfsarg_solveh48_t *, task_queue_t *);
-STATIC void init_queue(task_queue_t *);
-STATIC void submit_task(task_queue_t *, dfsarg_solveh48_t);
-STATIC void copy_queue(task_queue_t *, task_queue_t *, int, _Atomic int64_t *);
-STATIC void *start_thread(void *);
-STATIC int64_t solve_h48_bfs(dfsarg_solveh48_t *, task_queue_t *, int8_t);
-STATIC int64_t solve_h48_single(dfsarg_solveh48_t *, task_queue_t *);
-STATIC int64_t solve_h48_multithread(cube_t, int8_t, int8_t, int8_t, uint64_t,
-    const void *, uint64_t, char *, long long [static NISSY_SIZE_SOLVE_STATS]);
-
-STATIC void
-solve_h48_appendsolution_thread(dfsarg_solveh48_t *arg, task_queue_t *tq)
-{
-	pthread_mutex_lock(&tq->mutex);
-	int64_t strl = 0;
-	uint8_t invertedpremoves[MAXLEN];
-	char *solution = *arg->nextsol;
-
-	strl = writemoves(
-	    arg->moves, arg->nmoves, arg->solutions_size, *arg->nextsol);
-
-	if (strl < 0)
-		goto solve_h48_appendsolution_thread_error;
-	*arg->nextsol += strl-1;
-	arg->solutions_size -= strl-1;
-
-	if (arg->npremoves)
-	{
-		**arg->nextsol = ' ';
-		(*arg->nextsol)++;
-		arg->solutions_size--;
-
-		invertmoves(arg->premoves, arg->npremoves, invertedpremoves);
-		strl = writemoves(invertedpremoves,
-		    arg->npremoves, arg->solutions_size, *arg->nextsol);
-
-		if (strl < 0)
-			goto solve_h48_appendsolution_thread_error;
-		*arg->nextsol += strl-1;
-		arg->solutions_size -= strl-1;
-	}
-	LOG("Solution found: %s\n", solution);
-
-	**arg->nextsol = '\n';
-	(*arg->nextsol)++;
-	arg->solutions_size--;
-	(*arg->nsols)++;
-
-solve_h48_appendsolution_thread_error:
-	/* We could add some logging, but writemoves() already does */
-	pthread_mutex_unlock(&tq->mutex);
-}
-
-STATIC void
-init_queue(task_queue_t *queue)
-{
-	queue->front = 0;
-	queue->rear = 0;
-	queue->tasks_count = 0;
-	queue->active = 0;
-	queue->terminate = ATOMIC_VAR_INIT(false);
-	pthread_mutex_init(&queue->mutex, NULL);
-	pthread_cond_init(&queue->cond, NULL);
-	pthread_cond_init(&queue->active_cond, NULL);
-}
-
-STATIC void
-submit_task(task_queue_t *queue, dfsarg_solveh48_t task)
-{
-	pthread_mutex_lock(&queue->mutex);
-	queue->tasks[queue->rear] = task;
-	queue->rear = (queue->rear + 1) % MAX_QUEUE_SIZE;
-	queue->tasks_count++;
-	pthread_cond_broadcast(&queue->cond);
-	pthread_mutex_unlock(&queue->mutex);
-}
-
-STATIC void
-copy_queue(task_queue_t *src, task_queue_t *dest, int depth, _Atomic int64_t *nsols)
-{
-	pthread_mutex_lock(&dest->mutex);
-	for (int i = src->front; i != src->rear; i = (i + 1) % MAX_QUEUE_SIZE)
-	{
-			dest->tasks[i] = src->tasks[i];
-			dest->tasks[i].depth = depth;
-	}
-	dest->front = src->front;
-	dest->rear = src->rear;
-	dest->tasks_count = src->tasks_count;
-	pthread_cond_broadcast(&dest->cond);
-	pthread_mutex_unlock(&dest->mutex);
-}
-
-STATIC void *
-start_thread(void *arg)
-{
-	task_queue_t *queue = (task_queue_t *)arg;
-	while (1) {
-		pthread_mutex_lock(&queue->mutex);
-		while (queue->tasks_count == 0 && !queue->terminate) {
-			pthread_cond_wait(&queue->cond, &queue->mutex);
-		}
-		if (queue->tasks_count == 0 && queue->terminate) {
-			pthread_mutex_unlock(&queue->mutex);
-			break;
-		}
-
-		if (queue->tasks_count > 0) {
-			dfsarg_solveh48_t task = queue->tasks[queue->front];
-			queue->front = (queue->front + 1) % MAX_QUEUE_SIZE;
-			queue->tasks_count--;
-			queue->active++;
-			pthread_mutex_unlock(&queue->mutex);
-			
-			solve_h48_single(&task, queue);
-			queue->nodes_visited_global += task.nodes_visited;
-			queue->table_fallbacks_global += task.table_fallbacks;
-
-			pthread_mutex_lock(&queue->mutex);
-			queue->active--;
-
-			if(queue->tasks_count == 0 && queue->active == 0)
-				pthread_cond_signal(&queue->active_cond);
-		}
-		pthread_mutex_unlock(&queue->mutex);
-	}
-	return NULL;
-}
-
-STATIC int64_t
-solve_h48_bfs(dfsarg_solveh48_t *arg_zero, task_queue_t *tq, int8_t maxmoves)
-{
-	dfsarg_solveh48_t queue[MAX_QUEUE_SIZE];
-	int front = 0, rear = 0;
-	dfsarg_solveh48_t nextarg;
-	int depth = 0;
-	int nodes_at_current_depth = 1;
-	int nodes_at_next_depth = 0;
-	queue[rear++] = *arg_zero;
-
-	dfsarg_solveh48_t task_pool[MAX_QUEUE_SIZE];
-
-	while (front < rear){
-		dfsarg_solveh48_t arg = queue[front++];
-		nodes_at_current_depth--;
-
-		if (*arg.nsols == arg.maxsolutions)
-			return 1;
-
-		if (issolved(arg.cube)){
-			if (arg.nmoves + arg.npremoves >= arg.depth && arg.nmoves + arg.npremoves <= maxmoves)
-			solve_h48_appendsolution(&arg);
-			continue;
-		}
-
-		arg.nissbranch = MM_NORMAL;
-		uint32_t allowed = allowednextmove_h48(arg.moves, arg.nmoves, arg.nissbranch);
-
-		for (uint8_t m = 0; m < 18; m++){
-			if (allowed & (1 << m)){
-				nextarg = arg;
-				nextarg.nmoves = arg.nmoves + 1;
-				nextarg.moves[arg.nmoves] = m;
-				nextarg.cube = move(arg.cube, m);
-				nextarg.inverse = premove(arg.inverse, m);
-
-				if (nextarg.nmoves == BFS_DEPTH){
-					dfsarg_solveh48_t *task = &task_pool[rear % MAX_QUEUE_SIZE];
-					*task = nextarg;
-					submit_task(tq, *task);
-				} else {
-					queue[rear++] = nextarg;
-					nodes_at_next_depth++;
-				}
-			}
-		}
-		if (nodes_at_current_depth == 0){
-			nodes_at_current_depth = nodes_at_next_depth;
-			nodes_at_next_depth = 0;
-			LOG("Found %" PRId64 " solutions, searching at depth %" PRId8 "\n", *nextarg.nsols, depth++);
-		}
-		if (depth == BFS_DEPTH) return 0;
-	}
-	return 1;
-}
-
-STATIC int64_t
-solve_h48_single(dfsarg_solveh48_t *arg, task_queue_t *tq)
-{
-	dfsarg_solveh48_t nextarg;
-	int64_t ret;
-	uint8_t m;
-
-	if (*arg->nsols == arg->maxsolutions)
-		return 0;
-
-	if (solve_h48_stop(arg))
-		return 0;
-
-	if (issolved(arg->cube)){
-		if (arg->nmoves + arg->npremoves != arg->depth)
-			return 0;
-		solve_h48_appendsolution_thread(arg, tq);
-		return 1;
-	}
-
-	nextarg = *arg;
-	ret = 0;
-	uint32_t allowed;
-	if (arg->nissbranch & MM_INVERSE){
-		allowed = allowednextmove_h48(arg->premoves, arg->npremoves, arg->nissbranch);
-		for (m = 0; m < 18; m++){
-			if (allowed & (1 << m)){
-				nextarg.npremoves = arg->npremoves + 1;
-				nextarg.premoves[arg->npremoves] = m;
-				nextarg.inverse = move(arg->inverse, m);
-				nextarg.cube = premove(arg->cube, m);
-				ret += solve_h48_single(&nextarg, tq);
-			}
-		}
-	} else {
-		allowed = allowednextmove_h48(arg->moves, arg->nmoves, arg->nissbranch);
-		for (m = 0; m < 18; m++){
-			if (allowed & (1 << m)){
-				nextarg.nmoves = arg->nmoves + 1;
-				nextarg.moves[arg->nmoves] = m;
-				nextarg.cube = move(arg->cube, m);
-				nextarg.inverse = premove(arg->inverse, m);
-				ret += solve_h48_single(&nextarg, tq);
-			}
-		}
-	}
-
-	arg->nodes_visited = nextarg.nodes_visited;
-	arg->table_fallbacks = nextarg.table_fallbacks;
-	return ret;
-}
-
-STATIC int64_t
-solve_h48_multithread(
-	cube_t cube,
-	int8_t minmoves,
-	int8_t maxmoves,
-	int8_t maxsolutions,
-	uint64_t data_size,
-	const void *data,
-	uint64_t solutions_size,
-	char *solutions,
-	long long stats[static NISSY_SIZE_SOLVE_STATS]
-)
-{
-	_Atomic int64_t nsols = 0;
-	int p_depth = 0;
-	dfsarg_solveh48_t arg;
-	tableinfo_t info, fbinfo;
-	pthread_t threads[THREADS];
-
-	if (readtableinfo_n(data_size, data, 2, &info) != NISSY_OK)
-		goto solve_h48_multithread_error_data;
-
-	arg = (dfsarg_solveh48_t){
-		.cube = cube,
-		.inverse = inverse(cube),
-		.nsols = &nsols,
-		.depth = minmoves,
-		.maxsolutions = maxsolutions,
-		.h = info.h48h,
-		.k = info.bits,
-		.base = info.base,
-		.cocsepdata = (uint32_t *)((char *)data + INFOSIZE),
-		.h48data = (uint8_t *)data + COCSEP_FULLSIZE + INFOSIZE,
-		.solutions_size = solutions_size,
-		.nextsol = &solutions,
-		.nodes_visited = 0,
-		.table_fallbacks = 0
-	};
-
-	if (info.bits == 2) {
-		if (readtableinfo_n(data_size, data, 3, &fbinfo) != NISSY_OK)
-			goto solve_h48_multithread_error_data;
-		/* We only support h0k4 as fallback table */
-		if (fbinfo.h48h != 0 || fbinfo.bits != 4)
-			goto solve_h48_multithread_error_data;
-		arg.h48data_fallback = arg.h48data + info.next;
-	} else {
-		arg.h48data_fallback = NULL;
-	}
-
-	task_queue_t q;
-	init_queue(&q);
-	if (solve_h48_bfs(&arg, &q, maxmoves))
-		return nsols;
-
-	task_queue_t nq;
-	init_queue(&nq);
-
-	nq.nodes_visited_global = nq.table_fallbacks_global = 0;
-	for (int i = 0; i < THREADS; i++) {
-		pthread_create(&threads[i], NULL, &start_thread, &nq);
-	}
-
-	nsols = 0;
-	for (p_depth = minmoves > BFS_DEPTH ? minmoves : BFS_DEPTH;
-		 p_depth <= maxmoves && nsols < maxsolutions;
-		 p_depth++)
-	{
-		LOG("Found %" PRId64 " solutions, "
-		    "searching at depth %" PRId8 "\n", nsols, p_depth);
-		copy_queue(&q, &nq, p_depth, &nsols);
-		
-		pthread_mutex_lock(&nq.mutex);
-		while (nq.active > 0 || nq.tasks_count > 0)
-			pthread_cond_wait(&nq.active_cond, &nq.mutex);
-		pthread_mutex_unlock(&nq.mutex);
-	}
-
-	atomic_store(&nq.terminate, true);
-	pthread_cond_broadcast(&nq.cond);
-
-	for (int i = 0; i < THREADS; i++) {
-		pthread_join(threads[i], NULL);
-	}
-	**arg.nextsol = '\0';
-
-	stats[0] = nq.nodes_visited_global;
-	stats[1] = nq.table_fallbacks_global;
-	LOG("Nodes visited: %lld\nTable fallbacks: %lld\n",
-	    nq.nodes_visited_global, nq.table_fallbacks_global);
-
-	return nsols;
-
-solve_h48_multithread_error_data: LOG("solve_h48: error reading table\n"); return NISSY_ERROR_DATA;
-}
diff --git a/src/solvers/h48/stats.h b/src/solvers/h48/stats.h
deleted file mode 100644
index 2a7e887..0000000
--- a/src/solvers/h48/stats.h
+++ /dev/null
@@ -1,100 +0,0 @@
-/*
-The h48stats solver computes how many moves it takes to solve to
-each of the 12 h48 coordinates, one for each value of h from 0 to 11.
-The solutions array is filled with the length of the solutions. The
-solutions array is therefore not a printable string.
-*/
-
-typedef struct {
-	cube_t cube;
-	int8_t nmoves;
-	int8_t depth;
-	uint8_t moves[MAXLEN];
-	const uint32_t *cocsepdata;
-	const uint8_t *h48data;
-	char *s;
-} dfsarg_solveh48stats_t;
-
-STATIC int64_t solve_h48stats_dfs(dfsarg_solveh48stats_t *);
-STATIC int64_t solve_h48stats(cube_t, int8_t, const void *, char [static 12]);
-
-STATIC int64_t
-solve_h48stats_dfs(dfsarg_solveh48stats_t *arg)
-{
-	const int64_t limit = 11;
-
-	int8_t bound, u;
-	uint8_t m;
-	uint32_t d;
-	int64_t coord, h;
-	dfsarg_solveh48stats_t nextarg;
-
-	/* Check cocsep lower bound (corners only) */
-	bound = get_h48_cdata(arg->cube, arg->cocsepdata, &d);
-	if (bound + arg->nmoves > arg->depth)
-		return 0;
-
-	/* Check h48 lower bound for h=0 (esep, but no eo) */
-	coord = coord_h48_edges(arg->cube, COCLASS(d), TTREP(d), 0);
-	bound = get_h48_pval(arg->h48data, coord, 4);
-	if (bound + arg->nmoves > arg->depth)
-		return 0;
-	
-	/* Update all other values, if solved */
-	coord = coord_h48_edges(arg->cube, COCLASS(d), TTREP(d), 11);
-	for (h = 0; h <= limit; h++) {
-		u = coord >> (11-h) == 0 && arg->s[h] == 99;
-		arg->s[h] = u * arg->nmoves + (1-u) * arg->s[h];
-	}
-
-	if (arg->s[limit] != 99)
-		return 0;
-
-	nextarg = *arg;
-	nextarg.nmoves = arg->nmoves + 1;
-	for (m = 0; m < 18; m++) {
-		nextarg.moves[arg->nmoves] = m;
-		if (!allowednextmove(nextarg.moves, nextarg.nmoves)) {
-			/* If a move is not allowed, neither are its 180
-			 * and 270 degree variations */
-			m += 2;
-			continue;
-		}
-		nextarg.cube = move(arg->cube, m);
-		solve_h48stats_dfs(&nextarg);
-	}
-
-	return 0;
-}
-
-STATIC int64_t
-solve_h48stats(
-	cube_t cube,
-	int8_t maxmoves,
-	const void *data,
-	char solutions[static 12]
-)
-{
-	int i;
-	dfsarg_solveh48stats_t arg;
-
-	arg = (dfsarg_solveh48stats_t) {
-		.cube = cube,
-		.cocsepdata = (uint32_t *)((char *)data + INFOSIZE),
-		.h48data = (uint8_t *)data + COCSEP_FULLSIZE + INFOSIZE,
-		.s = solutions
-	};
-
-	for (i = 0; i < 12; i++)
-		solutions[i] = (char)99;
-
-	for (arg.depth = 0;
-	     arg.depth <= maxmoves && solutions[11] == 99;
-	     arg.depth++)
-	{
-		arg.nmoves = 0;
-		solve_h48stats_dfs(&arg);
-	}
-
-	return 0;
-}
diff --git a/src/utils/constants.h b/src/utils/constants.h
index ab1b9d8..dd2b92a 100644
--- a/src/utils/constants.h
+++ b/src/utils/constants.h
@@ -93,12 +93,11 @@ STATIC int64_t binomial[12][12] = {
 #define TRANS_BDm UINT8_C(46)
 #define TRANS_BLm UINT8_C(47)
 
-#define MM_NORMAL         UINT8_C(0x00)
-#define MM_INVERSE        UINT8_C(0x01)
-#define MM_INVERSEBRANCH  UINT8_C(0x03)
-#define MM_NORMALBRANCH   UINT8_C(0x02)
-#define MM_ALLMOVES       UINT32_C(0x3FFFF)
-#define MM_NOHALFTURNS    UINT32_C(0x2DB6D)
+#define NMOVES  (1+MOVE_B3)
+#define NTRANS  (1+TRANS_BLr)
+
+#define MM_ALLMOVES    UINT32_C(0x3FFFF)
+#define MM_NOHALFTURNS UINT32_C(0x2DB6D)
 
 #define CORNER_UFR      UINT8_C(0)
 #define CORNER_UBL      UINT8_C(1)
@@ -296,3 +295,54 @@ static uint8_t inverse_trans_table[48] = {
 	[TRANS_BLr] = TRANS_RDr,
 	[TRANS_BLm] = TRANS_LDm,
 };
+
+static uint8_t trans_move_table[48][3] = {
+	[TRANS_UFr] = { MOVE_U, MOVE_R, MOVE_F },
+	[TRANS_UFm] = { MOVE_U, MOVE_L, MOVE_F },
+	[TRANS_ULr] = { MOVE_U, MOVE_F, MOVE_L },
+	[TRANS_ULm] = { MOVE_U, MOVE_F, MOVE_R },
+	[TRANS_UBr] = { MOVE_U, MOVE_L, MOVE_B },
+	[TRANS_UBm] = { MOVE_U, MOVE_R, MOVE_B },
+	[TRANS_URr] = { MOVE_U, MOVE_B, MOVE_R },
+	[TRANS_URm] = { MOVE_U, MOVE_B, MOVE_L },
+	[TRANS_DFr] = { MOVE_D, MOVE_L, MOVE_F },
+	[TRANS_DFm] = { MOVE_D, MOVE_R, MOVE_F },
+	[TRANS_DLr] = { MOVE_D, MOVE_B, MOVE_L },
+	[TRANS_DLm] = { MOVE_D, MOVE_B, MOVE_R },
+	[TRANS_DBr] = { MOVE_D, MOVE_R, MOVE_B },
+	[TRANS_DBm] = { MOVE_D, MOVE_L, MOVE_B },
+	[TRANS_DRr] = { MOVE_D, MOVE_F, MOVE_R },
+	[TRANS_DRm] = { MOVE_D, MOVE_F, MOVE_L },
+	[TRANS_RUr] = { MOVE_R, MOVE_F, MOVE_U },
+	[TRANS_RUm] = { MOVE_L, MOVE_F, MOVE_U },
+	[TRANS_RFr] = { MOVE_R, MOVE_D, MOVE_F },
+	[TRANS_RFm] = { MOVE_L, MOVE_D, MOVE_F },
+	[TRANS_RDr] = { MOVE_R, MOVE_B, MOVE_D },
+	[TRANS_RDm] = { MOVE_L, MOVE_B, MOVE_D },
+	[TRANS_RBr] = { MOVE_R, MOVE_U, MOVE_B },
+	[TRANS_RBm] = { MOVE_L, MOVE_U, MOVE_B },
+	[TRANS_LUr] = { MOVE_L, MOVE_B, MOVE_U },
+	[TRANS_LUm] = { MOVE_R, MOVE_B, MOVE_U },
+	[TRANS_LFr] = { MOVE_L, MOVE_U, MOVE_F },
+	[TRANS_LFm] = { MOVE_R, MOVE_U, MOVE_F },
+	[TRANS_LDr] = { MOVE_L, MOVE_F, MOVE_D },
+	[TRANS_LDm] = { MOVE_R, MOVE_F, MOVE_D },
+	[TRANS_LBr] = { MOVE_L, MOVE_D, MOVE_B },
+	[TRANS_LBm] = { MOVE_R, MOVE_D, MOVE_B },
+	[TRANS_FUr] = { MOVE_F, MOVE_L, MOVE_U },
+	[TRANS_FUm] = { MOVE_F, MOVE_R, MOVE_U },
+	[TRANS_FRr] = { MOVE_F, MOVE_U, MOVE_R },
+	[TRANS_FRm] = { MOVE_F, MOVE_U, MOVE_L },
+	[TRANS_FDr] = { MOVE_F, MOVE_R, MOVE_D },
+	[TRANS_FDm] = { MOVE_F, MOVE_L, MOVE_D },
+	[TRANS_FLr] = { MOVE_F, MOVE_D, MOVE_L },
+	[TRANS_FLm] = { MOVE_F, MOVE_D, MOVE_R },
+	[TRANS_BUr] = { MOVE_B, MOVE_R, MOVE_U },
+	[TRANS_BUm] = { MOVE_B, MOVE_L, MOVE_U },
+	[TRANS_BRr] = { MOVE_B, MOVE_D, MOVE_R },
+	[TRANS_BRm] = { MOVE_B, MOVE_D, MOVE_L },
+	[TRANS_BDr] = { MOVE_B, MOVE_L, MOVE_D },
+	[TRANS_BDm] = { MOVE_B, MOVE_R, MOVE_D },
+	[TRANS_BLr] = { MOVE_B, MOVE_U, MOVE_L },
+	[TRANS_BLm] = { MOVE_B, MOVE_U, MOVE_R },
+};
-- 
cgit v1.3