5 files changed, 999 insertions, 0 deletions

diff --git a/src/arch/arch.h b/src/arch/arch.h
new file mode 100644
index 0000000..d5e0218
--- /dev/null
+++ b/src/arch/arch.h
@@ -0,0 +1,41 @@
+#if defined(CUBE_AVX2)
+
+#include <immintrin.h>
+
+typedef __m256i cube_t;
+
+#if !defined(TEST_H)
+#include "common.h"
+#include "avx2.h"
+#endif
+
+#elif defined(CUBE_NEON)
+
+#include <stdlib.h>
+#include <arm_neon.h>
+
+typedef struct {
+        uint8x16_t corner;
+        uint8x16_t edge;
+} cube_t;
+
+#if !defined(TEST_H)
+#include "common.h"
+#include "neon.h"
+#endif
+
+#else
+
+#include <stdlib.h>
+
+typedef struct {
+        uint8_t corner[8];
+        uint8_t edge[12];
+} cube_t;
+
+#if !defined(TEST_H)
+#include "common.h"
+#include "portable.h"
+#endif
+
+#endif
diff --git a/src/arch/avx2.h b/src/arch/avx2.h
new file mode 100644
index 0000000..59a0d6c
--- /dev/null
+++ b/src/arch/avx2.h
@@ -0,0 +1,319 @@
+#define _co2_avx2 _mm256_set_epi64x(0, 0, 0, INT64_C(0x6060606060606060))
+#define _cocw_avx2 _mm256_set_epi64x(0, 0, 0, INT64_C(0x2020202020202020))
+#define _cp_avx2 _mm256_set_epi64x(0, 0, 0, INT64_C(0x0707070707070707))
+#define _ep_avx2 \
+    _mm256_set_epi64x(INT64_C(0x0F0F0F0F), INT64_C(0x0F0F0F0F0F0F0F0F), 0, 0)
+#define _eo_avx2 \
+    _mm256_set_epi64x(INT64_C(0x10101010), INT64_C(0x1010101010101010), 0, 0)
+
+#define static_cube(c_ufr, c_ubl, c_dfl, c_dbr, c_ufl, c_ubr, c_dfr, c_dbl, \
+    e_uf, e_ub, e_db, e_df, e_ur, e_ul, e_dl, e_dr, e_fr, e_fl, e_bl, e_br) \
+    _mm256_set_epi8(0, 0, 0, 0, e_br, e_bl, e_fl, e_fr, \
+        e_dr, e_dl, e_ul, e_ur, e_df, e_db, e_ub, e_uf, \
+        0, 0, 0, 0, 0, 0, 0, 0, \
+        c_dbl, c_dfr, c_ubr, c_ufl, c_dbr, c_dfl, c_ubl, c_ufr)
+#define zero _mm256_set_epi64x(0, 0, 0, 0)
+#define solved static_cube( \
+    0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
+
+_static void
+pieces(cube_t *cube, uint8_t c[static 8], uint8_t e[static 12])
+{
+        uint8_t aux[32];
+
+        _mm256_storeu_si256((__m256i_u *)aux, *cube);
+        memcpy(c, aux, 8);
+        memcpy(e, aux+16, 12);
+}
+
+_static_inline bool
+equal(cube_t c1, cube_t c2)
+{
+        int32_t mask;
+        __m256i cmp;
+
+        cmp = _mm256_cmpeq_epi8(c1, c2);
+        mask = _mm256_movemask_epi8(cmp);
+
+        return mask == ~0;
+}
+
+_static_inline cube_t
+invertco(cube_t c)
+{
+        cube_t co, shleft, shright, summed, newco, cleanco, ret;
+
+        co = _mm256_and_si256(c, _co2_avx2);
+        shleft = _mm256_slli_epi32(co, 1);
+        shright = _mm256_srli_epi32(co, 1);
+        summed = _mm256_or_si256(shleft, shright);
+        newco = _mm256_and_si256(summed, _co2_avx2);
+        cleanco = _mm256_xor_si256(c, co);
+        ret = _mm256_or_si256(cleanco, newco);
+
+        return ret;
+}
+
+_static_inline cube_t
+compose_epcpeo(cube_t c1, cube_t c2)
+{
+        cube_t b, s, eo2;
+
+        /* Permute and clean unused bits */
+        s = _mm256_shuffle_epi8(c1, c2);
+        b = _mm256_set_epi8(
+                ~0, ~0, ~0, ~0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                ~0, ~0, ~0, ~0, ~0, ~0, ~0, ~0, 0, 0, 0, 0, 0, 0, 0, 0
+        );
+        s = _mm256_andnot_si256(b, s);
+
+        /* Change EO */
+        eo2 = _mm256_and_si256(c2, _eo_avx2);
+        s = _mm256_xor_si256(s, eo2);
+
+        return s;
+}
+
+_static_inline cube_t
+compose_edges(cube_t c1, cube_t c2)
+{
+        return compose_epcpeo(c1, c2);
+}
+
+_static_inline cube_t
+compose_corners(cube_t c1, cube_t c2)
+{
+        /*
+         * We do a full compose. Minor optimizations are possible, like
+         * saving one instruction by not doing EO, but it should not
+         * be significant.
+         */
+        return compose(c1, c2);
+}
+
+_static_inline cube_t
+compose(cube_t c1, cube_t c2)
+{
+        cube_t s, co1, co2, aux, auy1, auy2, auz1, auz2;
+
+        s = compose_epcpeo(c1, c2);
+
+        /* Change CO */
+        co1 = _mm256_and_si256(s, _co2_avx2);
+        co2 = _mm256_and_si256(c2, _co2_avx2);
+        aux = _mm256_add_epi8(co1, co2);
+        auy1 = _mm256_add_epi8(aux, _cocw_avx2);
+        auy2 = _mm256_srli_epi32(auy1, 2);
+        auz1 = _mm256_add_epi8(aux, auy2);
+        auz2 = _mm256_and_si256(auz1, _co2_avx2);
+
+        /* Put together */
+        s = _mm256_andnot_si256(_co2_avx2, s);
+        s = _mm256_or_si256(s, auz2);
+
+        return s;
+}
+
+_static_inline cube_t
+cleanaftershuffle(cube_t c)
+{
+        __m256i b;
+
+        b = _mm256_set_epi8(
+                ~0, ~0, ~0, ~0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                ~0, ~0, ~0, ~0, ~0, ~0, ~0, ~0, 0, 0, 0, 0, 0, 0, 0, 0
+        );
+
+        return _mm256_andnot_si256(b, c);
+}
+
+_static_inline cube_t
+inverse(cube_t c)
+{
+        /* Method taken from Andrew Skalski's vcube[1]. The addition sequence
+         * was generated using [2].
+         * [1] https://github.com/Voltara/vcube
+         * [2] http://wwwhomes.uni-bielefeld.de/achim/addition_chain.html
+         */
+        cube_t v3, vi, vo, vp, ret;
+
+        v3 = _mm256_shuffle_epi8(c, c);
+        v3 = _mm256_shuffle_epi8(v3, c);
+        vi = _mm256_shuffle_epi8(v3, v3);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, v3);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, c);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, v3);
+        vi = _mm256_shuffle_epi8(vi, vi);
+        vi = _mm256_shuffle_epi8(vi, c);
+
+        vo = _mm256_and_si256(c, _mm256_or_si256(_eo_avx2, _co2_avx2));
+        vo = _mm256_shuffle_epi8(vo, vi);
+        vp = _mm256_andnot_si256(_mm256_or_si256(_eo_avx2, _co2_avx2), vi);
+        ret = _mm256_or_si256(vp, vo);
+        ret = cleanaftershuffle(ret);
+        
+        return invertco(ret);
+}
+
+_static_inline int64_t
+coord_co(cube_t c)
+{
+        cube_t co;
+        int64_t mem[4], ret, i, p;
+
+        co = _mm256_and_si256(c, _co2_avx2);
+        _mm256_storeu_si256((__m256i *)mem, co);
+
+        mem[0] >>= 5;
+        for (i = 0, ret = 0, p = 1; i < 7; i++, mem[0] >>= 8, p *= 3)
+                ret += (mem[0] & 3) * p;
+
+        return ret;
+}
+
+_static_inline int64_t
+coord_csep(cube_t c)
+{
+        cube_t cp, shifted;
+        int64_t mask;
+
+        cp = _mm256_and_si256(c, _cp_avx2);
+        shifted = _mm256_slli_epi32(cp, 5);
+        mask = _mm256_movemask_epi8(shifted);
+
+        return mask & 0x7F;
+}
+
+_static_inline int64_t
+coord_cocsep(cube_t c)
+{
+        return (coord_co(c) << 7) + coord_csep(c);
+}
+
+_static_inline int64_t
+coord_eo(cube_t c)
+{
+        cube_t eo, shifted;
+        int64_t mask;
+
+        eo = _mm256_and_si256(c, _eo_avx2);
+        shifted = _mm256_slli_epi32(eo, 3);
+        mask = _mm256_movemask_epi8(shifted);
+
+        return mask >> 17;
+}
+
+_static_inline int64_t
+coord_esep(cube_t c)
+{
+        cube_t ep;
+        int64_t e, mem[4], i, j, jj, k, l, ret1, ret2, bit1, bit2, is1;
+
+        ep = _mm256_and_si256(c, _ep_avx2);
+        _mm256_storeu_si256((__m256i *)mem, ep);
+
+        mem[3] <<= 8;
+        ret1 = ret2 = 0;
+        k = l = 4;
+        for (i = 0, j = 0; i < 12; i++, mem[i/8 + 2] >>= 8) {
+                e = mem[i/8 + 2];
+
+                bit1 = (e & _esepbit1) >> 2;
+                bit2 = (e & _esepbit2) >> 3;
+                is1 = (1 - bit2) * bit1;
+
+                ret1 += bit2 * binomial[11-i][k];
+                k -= bit2;
+
+                jj = j < 8;
+                ret2 += jj * is1 * binomial[7-(j*jj)][l];
+                l -= is1;
+                j += (1-bit2);
+        }
+
+        return ret1 * 70 + ret2;
+}
+
+_static_inline void
+copy_corners(cube_t *dest, cube_t src)
+{
+        *dest = _mm256_blend_epi32(*dest, src, 0x0F);
+}
+
+_static_inline void
+copy_edges(cube_t *dest, cube_t src)
+{
+        *dest = _mm256_blend_epi32(*dest, src, 0xF0);
+}
+
+_static_inline void
+set_eo(cube_t *cube, int64_t eo)
+{
+        int64_t eo12, eotop, eobot;
+        __m256i veo;
+
+        eo12 = (eo << 1) + (_mm_popcnt_u64(eo) % 2);
+        eotop = (eo12 & (1 << 11)) << 17 |
+                (eo12 & (1 << 10)) << 10 |
+                (eo12 & (1 << 9)) << 3 |
+                (eo12 & (1 << 8)) >> 4;
+        eobot = (eo12 & (1 << 7)) << 53 |
+                (eo12 & (1 << 6)) << 46 |
+                (eo12 & (1 << 5)) << 39 |
+                (eo12 & (1 << 4)) << 32 |
+                (eo12 & (1 << 3)) << 25 |
+                (eo12 & (1 << 2)) << 18 |
+                (eo12 & (1 << 1)) << 11 |
+                (eo12 & 1) << 4;
+        veo = _mm256_set_epi64x(eotop, eobot, 0, 0);
+
+        *cube = _mm256_andnot_si256(_eo_avx2, *cube);
+        *cube = _mm256_or_si256(*cube, veo);
+}
+
+_static_inline cube_t
+invcoord_esep(int64_t esep)
+{
+        cube_t eee, ret;
+        int64_t bit1, bit2, i, j, jj, k, l, s, v, w, is1, set1, set2;
+        uint8_t mem[32];
+        uint8_t slice[3] = {0};
+
+        set1 = esep % 70;
+        set2 = esep / 70;
+
+        for (i = 0, j = 0, k = 4, l = 4; i < 12; i++) {
+                v = binomial[11-i][k];
+                jj = j < 8;
+                w = jj * binomial[7-(j*jj)][l];
+                bit2 = set2 >= v;
+                bit1 = set1 >= w;
+                is1 = (1 - bit2) * bit1;
+
+                set2 -= bit2 * v;
+                k -= bit2;
+                set1 -= is1 * w;
+                l -= is1;
+                j += (1-bit2);
+                s = 2*bit2 + (1-bit2)*bit1;
+
+                mem[i+16] = (slice[s]++) | (uint8_t)(s << 2);
+        }
+
+        ret = solved;
+        eee = _mm256_loadu_si256((__m256i_u *)&mem);
+        copy_edges(&ret, eee);
+
+        return ret;
+}
diff --git a/src/arch/common.h b/src/arch/common.h
new file mode 100644
index 0000000..cd2d36e
--- /dev/null
+++ b/src/arch/common.h
@@ -0,0 +1,19 @@
+_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);
+_static_inline cube_t compose_epcpeo(cube_t, cube_t);
+_static_inline cube_t compose_edges(cube_t, cube_t);
+_static_inline cube_t compose_corners(cube_t, cube_t);
+_static_inline cube_t compose(cube_t, cube_t);
+_static_inline cube_t inverse(cube_t);
+
+_static_inline int64_t coord_co(cube_t);
+_static_inline int64_t coord_csep(cube_t);
+_static_inline int64_t coord_cocsep(cube_t);
+_static_inline int64_t coord_eo(cube_t);
+_static_inline int64_t coord_esep(cube_t);
+
+_static_inline void copy_corners(cube_t *, cube_t);
+_static_inline void copy_edges(cube_t *, cube_t);
+_static_inline void set_eo(cube_t *, int64_t);
+_static_inline cube_t invcoord_esep(int64_t);
diff --git a/src/arch/neon.h b/src/arch/neon.h
new file mode 100644
index 0000000..a75f86d
--- /dev/null
+++ b/src/arch/neon.h
@@ -0,0 +1,348 @@
+#define _co2_neon vdupq_n_u8(0x60)
+#define _cocw_neon vdupq_n_u8(0x20)
+#define _cp_neon vdupq_n_u8(0x07)
+#define _ep_neon vcombine_u8(vdupq_n_u8(0x0F), vdupq_n_u8(0x0F))
+#define _eo_neon vcombine_u8(vdupq_n_u8(0x10), vdupq_n_u8(0x10))
+
+// static cube
+#define static_cube(c_ufr, c_ubl, c_dfl, c_dbr, c_ufl, c_ubr, c_dfr, c_dbl, \
+                                        e_uf, e_ub, e_db, e_df, e_ur, e_ul, e_dl, e_dr, e_fr, e_fl, e_bl, e_br) \
+        ((cube_t){ \
+                .corner = {c_ufr, c_ubl, c_dfl, c_dbr, c_ufl, c_ubr, c_dfr, c_dbl, 0, 0, 0, 0, 0, 0, 0, 0}, \
+                .edge = {e_uf, e_ub, e_db, e_df, e_ur, e_ul, e_dl, e_dr, e_fr, e_fl, e_bl, e_br, 0, 0, 0, 0}})
+
+// zero cube
+#define zero \
+        (cube_t) \
+        { \
+                .corner = vdupq_n_u8(0), \
+                .edge = vdupq_n_u8(0) \
+        }
+
+// solved cube
+#define solved static_cube(     \
+        0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
+
+_static void
+pieces(cube_t *cube, uint8_t c[static 8], uint8_t e[static 12])
+{
+        // First 8 bytes of the corner vector are copied from the c array
+        vst1_u8(c, vget_low_u8(cube->corner));
+
+        // 12 bytes of the edge vector are copied from the e array
+        // First 8 bytes
+        vst1_u8(e, vget_low_u8(cube->edge));
+        // Next 4 bytes
+        vst1_lane_u32((uint32_t *)(e + 8), vreinterpret_u32_u8(vget_high_u8(cube->edge)), 0);
+}
+
+_static_inline bool
+equal(cube_t c1, cube_t c2)
+{
+        uint8x16_t cmp_corner, cmp_edge;
+        uint64x2_t cmp_corner_u64, cmp_edge_u64;
+        uint64x2_t cmp_result;
+
+        // compare the corner vectors
+        cmp_corner = vceqq_u8(c1.corner, c2.corner);
+        // compare the edge vectors
+        cmp_edge = vceqq_u8(c1.edge, c2.edge);
+
+        // convert the comparison vectors to 64-bit vectors
+        cmp_corner_u64 = vreinterpretq_u64_u8(cmp_corner);
+        cmp_edge_u64 = vreinterpretq_u64_u8(cmp_edge);
+
+        // combine the comparison vectors
+        cmp_result = vandq_u64(cmp_corner_u64, cmp_edge_u64);
+
+        // check if all the bits are set
+        return vgetq_lane_u64(cmp_result, 0) == ~0ULL && vgetq_lane_u64(cmp_result, 1) == ~0ULL;
+}
+
+_static_inline cube_t
+invertco(cube_t c)
+{
+        cube_t ret;
+        uint8x16_t co, shleft, shright, summed, newco, cleanco;
+
+        co = vandq_u8(c.corner, _co2_neon);
+        shleft = vshlq_n_u8(co, 1);
+        shright = vshrq_n_u8(co, 1);
+        summed = vorrq_u8(shleft, shright);
+        newco = vandq_u8(summed, _co2_neon);
+        cleanco = veorq_u8(c.corner, co);
+        ret.corner = vorrq_u8(cleanco, newco);
+        ret.edge = c.edge;
+        
+        return ret;
+}
+
+_static_inline cube_t
+compose_edges(cube_t c1, cube_t c2)
+{
+        cube_t ret = {0};
+        ret.edge = compose_edges_slim(c1.edge, c2.edge);
+        return ret;
+}
+
+_static_inline cube_t 
+compose_corners(cube_t c1, cube_t c2)
+{
+        cube_t ret = {0};
+        ret.corner = compose_corners_slim(c1.corner, c2.corner);
+        return ret;
+}
+
+_static_inline uint8x16_t 
+compose_edges_slim(uint8x16_t edge1, uint8x16_t edge2)
+{
+        // Masks
+        uint8x16_t p_bits = vdupq_n_u8(_pbits);
+        uint8x16_t eo_bit = vdupq_n_u8(_eobit);
+
+        // Find the index and permutation
+        uint8x16_t p = vandq_u8(edge2, p_bits);
+        uint8x16_t piece1 = vqtbl1q_u8(edge1, p);
+
+        // Calculate the orientation through XOR
+        uint8x16_t orien = vandq_u8(veorq_u8(edge2, piece1), eo_bit);
+
+        // Combine the results
+        uint8x16_t ret = vorrq_u8(vandq_u8(piece1, p_bits), orien);
+
+        // Mask to clear the last 32 bits of the result
+        uint8x16_t mask_last_32 = vsetq_lane_u32(0, vreinterpretq_u32_u8(ret), 3);
+        ret = vreinterpretq_u8_u32(mask_last_32);
+
+        return ret;
+}
+
+_static_inline uint8x16_t 
+compose_corners_slim(uint8x16_t corner1, uint8x16_t corner2)
+{
+        // Masks
+        uint8x16_t p_bits = vdupq_n_u8(_pbits);
+        uint8x16_t cobits = vdupq_n_u8(_cobits);
+        uint8x16_t cobits2 = vdupq_n_u8(_cobits2);
+        uint8x16_t twist_cw = vdupq_n_u8(_ctwist_cw);
+
+        // Find the index and permutation
+        uint8x16_t p = vandq_u8(corner2, p_bits);
+        uint8x16_t piece1 = vqtbl1q_u8(corner1, p);
+
+        // Calculate the orientation
+        uint8x16_t aux = vaddq_u8(vandq_u8(corner2, cobits), vandq_u8(piece1, cobits));
+        uint8x16_t auy = vshrq_n_u8(vaddq_u8(aux, twist_cw), 2);
+        uint8x16_t orien = vandq_u8(vaddq_u8(aux, auy), cobits2);
+
+        // Combine the results
+        uint8x16_t ret = vorrq_u8(vandq_u8(piece1, p_bits), orien);
+
+        // Mask to clear the last 64 bits of the result
+        uint8x16_t mask_last_64 = vsetq_lane_u64(0, vreinterpretq_u64_u8(ret), 1);
+        ret = vreinterpretq_u8_u64(mask_last_64);
+
+        return ret;
+}
+
+_static_inline cube_t 
+compose(cube_t c1, cube_t c2)
+{
+        cube_t ret = {0};
+
+        ret.edge = compose_edges_slim(c1.edge, c2.edge);
+        ret.corner = compose_corners_slim(c1.corner, c2.corner);
+
+        return ret;
+}
+
+_static_inline cube_t 
+inverse(cube_t cube)
+{
+        uint8_t i, piece, orien;
+        cube_t ret;
+
+        // Temp arrays to store the NEON vectors
+        uint8_t edges[16];
+        uint8_t corners[16];
+
+        // Copy the NEON vectors to the arrays
+        vst1q_u8(edges, cube.edge);
+        vst1q_u8(corners, cube.corner);
+
+        uint8_t edge_result[16] = {0};
+        uint8_t corner_result[16] = {0};
+
+        // Process the edges
+        for (i = 0; i < 12; i++)
+        {
+                piece = edges[i];
+                orien = piece & _eobit;
+                edge_result[piece & _pbits] = i | orien;
+        }
+
+        // Process the corners
+        for (i = 0; i < 8; i++)
+        {
+                piece = corners[i];
+                orien = ((piece << 1) | (piece >> 1)) & _cobits2;
+                corner_result[piece & _pbits] = i | orien;
+        }
+
+        // Copy the results back to the NEON vectors
+        ret.edge = vld1q_u8(edge_result);
+        ret.corner = vld1q_u8(corner_result);
+
+        return ret;
+}
+
+_static_inline int64_t 
+coord_co(cube_t c)
+{
+        // Temp array to store the NEON vector
+        uint8_t mem[16];
+        vst1q_u8(mem, c.corner);
+
+        int i, p;
+        int64_t ret;
+
+        for (ret = 0, i = 0, p = 1; i < 7; i++, p *= 3)
+                ret += p * (mem[i] >> _coshift);
+
+        return ret;
+}
+
+_static_inline int64_t
+coord_csep(cube_t c)
+{
+        // Temp array to store the NEON vector
+        uint8_t mem[16];
+        vst1q_u8(mem, c.corner);
+
+        int64_t ret = 0;
+        int i, p;
+        for (ret = 0, i = 0, p = 1; i < 7; i++, p *= 2)
+                ret += p * ((mem[i] & _csepbit) >> 2);
+
+        return ret;
+        return 0;
+}
+
+_static_inline int64_t
+coord_cocsep(cube_t c)
+{
+        return (coord_co(c) << 7) + coord_csep(c);
+}
+
+_static_inline int64_t 
+coord_eo(cube_t c)
+{
+        int64_t ret = 0;
+        int64_t p = 1;
+
+        // Temp array to store the NEON vector
+        uint8_t mem[16];
+        vst1q_u8(mem, c.edge);
+
+        for (int i = 1; i < 12; i++, p *= 2)
+        {
+                ret += p * (mem[i] >> _eoshift);
+        }
+
+        return ret;
+}
+
+_static_inline int64_t 
+coord_esep(cube_t c)
+{
+        int64_t i, j, jj, k, l, ret1, ret2, bit1, bit2, is1;
+
+        // Temp array to store the NEON vector
+        uint8_t mem[16];
+        vst1q_u8(mem, c.edge);
+
+        for (i = 0, j = 0, k = 4, l = 4, ret1 = 0, ret2 = 0; i < 12; i++)
+        {
+                bit1 = (mem[i] & _esepbit1) >> 2;
+                bit2 = (mem[i] & _esepbit2) >> 3;
+                is1 = (1 - bit2) * bit1;
+
+                ret1 += bit2 * binomial[11 - i][k];
+                k -= bit2;
+
+                jj = j < 8;
+                ret2 += jj * is1 * binomial[7 - (j * jj)][l];
+                l -= is1;
+                j += (1 - bit2);
+        }
+
+        return ret1 * 70 + ret2;
+}
+
+_static_inline void
+copy_corners(cube_t *dst, cube_t src)
+{
+        dst->corner = src.corner;
+}
+
+_static_inline void
+copy_edges(cube_t *dst, cube_t src)
+{
+        dst->edge = src.edge;
+}
+
+_static_inline void
+set_eo(cube_t *cube, int64_t eo)
+{
+        // Temp array to store the NEON vector
+        uint8_t mem[16];
+        vst1q_u8(mem, cube->edge);
+        uint8_t i, sum, flip;
+
+        for (sum = 0, i = 1; i < 12; i++, eo >>= 1)
+        {
+                flip = eo % 2;
+                sum += flip;
+                mem[i] = (mem[i] & ~_eobit) | (_eobit * flip);
+        }
+        mem[0] = (mem[0] & ~_eobit) | (_eobit * (sum % 2));
+
+        // Copy the results back to the NEON vector
+        cube->edge = vld1q_u8(mem);
+        return;
+}
+
+_static_inline cube_t
+invcoord_esep(int64_t esep)
+{
+        cube_t ret;
+        int64_t bit1, bit2, i, j, jj, k, l, s, v, w, is1, set1, set2;
+        uint8_t slice[3] = {0};
+
+        ret = solved;
+        uint8_t mem[16];
+        set1 = esep % 70;
+        set2 = esep / 70;
+
+        for (i = 0, j = 0, k = 4, l = 4; i < 12; i++)
+        {
+                v = binomial[11 - i][k];
+                jj = j < 8;
+                w = jj * binomial[7 - (j * jj)][l];
+                bit2 = set2 >= v;
+                bit1 = set1 >= w;
+                is1 = (1 - bit2) * bit1;
+
+                set2 -= bit2 * v;
+                k -= bit2;
+                set1 -= is1 * w;
+                l -= is1;
+                j += (1 - bit2);
+                s = 2 * bit2 + (1 - bit2) * bit1;
+
+                mem[i] = (slice[s]++) | (uint8_t)(s << 2);
+        }
+
+        ret.edge = vld1q_u8(mem);
+        return ret;
+}
diff --git a/src/arch/portable.h b/src/arch/portable.h
new file mode 100644
index 0000000..b8330af
--- /dev/null
+++ b/src/arch/portable.h
@@ -0,0 +1,272 @@
+#define static_cube(c_ufr, c_ubl, c_dfl, c_dbr, c_ufl, c_ubr, c_dfr, c_dbl, \
+    e_uf, e_ub, e_db, e_df, e_ur, e_ul, e_dl, e_dr, e_fr, e_fl, e_bl, e_br) \
+    ((cube_t) { \
+        .corner = { c_ufr, c_ubl, c_dfl, c_dbr, c_ufl, c_ubr, c_dfr, c_dbl }, \
+        .edge = { e_uf, e_ub, e_db, e_df, e_ur, e_ul, \
+                  e_dl, e_dr, e_fr, e_fl, e_bl, e_br } })
+#define zero static_cube( \
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
+#define solved static_cube( \
+    0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
+
+_static void
+pieces(cube_t *cube, uint8_t c[static 8], uint8_t e[static 12])
+{
+        memcpy(c, cube->corner, 8);
+        memcpy(e, cube->edge, 12);
+}
+
+_static_inline bool
+equal(cube_t c1, cube_t c2)
+{
+        uint8_t i;
+        bool ret;
+
+        ret = true;
+        for (i = 0; i < 8; i++)
+                ret = ret && c1.corner[i] == c2.corner[i];
+        for (i = 0; i < 12; i++)
+                ret = ret && c1.edge[i] == c2.edge[i];
+
+        return ret;
+}
+
+_static_inline cube_t
+invertco(cube_t c)
+{
+        uint8_t i, piece, orien;
+        cube_t ret;
+
+        ret = c;
+        for (i = 0; i < 8; i++) {
+                piece = c.corner[i];
+                orien = ((piece << 1) | (piece >> 1)) & _cobits2;
+                ret.corner[i] = (piece & _pbits) | orien;
+        }
+
+        return ret;
+}
+
+_static_inline void
+compose_edges_inplace(cube_t c1, cube_t c2, cube_t *ret)
+{
+        uint8_t i, piece1, piece2, p, orien;
+
+        for (i = 0; i < 12; i++) {
+                piece2 = c2.edge[i];
+                p = piece2 & _pbits;
+                piece1 = c1.edge[p];
+                orien = (piece2 ^ piece1) & _eobit;
+                ret->edge[i] = (piece1 & _pbits) | orien;
+        }
+}
+
+_static_inline void
+compose_corners_inplace(cube_t c1, cube_t c2, cube_t *ret)
+{
+        uint8_t i, piece1, piece2, p, orien, aux, auy;
+
+        for (i = 0; i < 8; i++) {
+                piece2 = c2.corner[i];
+                p = piece2 & _pbits;
+                piece1 = c1.corner[p];
+                aux = (piece2 & _cobits) + (piece1 & _cobits);
+                auy = (aux + _ctwist_cw) >> 2;
+                orien = (aux + auy) & _cobits2;
+                ret->corner[i] = (piece1 & _pbits) | orien;
+        }
+}
+
+_static_inline cube_t
+compose_edges(cube_t c1, cube_t c2)
+{
+        cube_t ret = zero;
+
+        compose_edges_inplace(c1, c2, &ret);
+
+        return ret;
+}
+
+_static_inline cube_t
+compose_corners(cube_t c1, cube_t c2)
+{
+        cube_t ret = zero;
+
+        compose_corners_inplace(c1, c2, &ret);
+
+        return ret;
+}
+
+_static_inline cube_t
+compose(cube_t c1, cube_t c2)
+{
+        cube_t ret = zero;
+
+        compose_edges_inplace(c1, c2, &ret);
+        compose_corners_inplace(c1, c2, &ret);
+
+        return ret;
+}
+
+cube_t
+inverse(cube_t cube)
+{
+        uint8_t i, piece, orien;
+        cube_t ret;
+
+        for (i = 0; i < 12; i++) {
+                piece = cube.edge[i];
+                orien = piece & _eobit;
+                ret.edge[piece & _pbits] = i | orien;
+        }
+
+        for (i = 0; i < 8; i++) {
+                piece = cube.corner[i];
+                orien = ((piece << 1) | (piece >> 1)) & _cobits2;
+                ret.corner[piece & _pbits] = i | orien;
+        }
+
+        return ret;
+}
+
+_static_inline int64_t
+coord_co(cube_t c)
+{
+        int i, p;
+        int64_t ret;
+
+        for (ret = 0, i = 0, p = 1; i < 7; i++, p *= 3)
+                ret += p * (c.corner[i] >> _coshift);
+
+        return ret;
+}
+
+/*
+For corner separation, we consider the axis (a.k.a. tetrad) each
+corner belongs to as 0 or 1 and we translate this sequence into binary.
+Ignoring the last bit, we have a value up to 2^7, but not all values are
+possible. Encoding this as a number from 0 to C(8,4) would save about 40%
+of space, but we are not going to use this coordinate in large tables.
+*/
+_static_inline int64_t
+coord_csep(cube_t c)
+{
+        int i, p;
+        int64_t ret;
+
+        for (ret = 0, i = 0, p = 1; i < 7; i++, p *= 2)
+                ret += p * ((c.corner[i] & _csepbit) >> 2);
+
+        return ret;
+}
+
+_static_inline int64_t
+coord_cocsep(cube_t c)
+{
+        return (coord_co(c) << 7) + coord_csep(c);
+}
+
+_static_inline int64_t
+coord_eo(cube_t c)
+{
+        int i, p;
+        int64_t ret;
+
+        for (ret = 0, i = 1, p = 1; i < 12; i++, p *= 2)
+                ret += p * (c.edge[i] >> _eoshift);
+
+        return ret;
+}
+
+/*
+We encode the edge separation as a number from 0 to C(12,4)*C(8,4).
+It can be seen as the composition of two "subset index" coordinates.
+*/
+_static_inline int64_t
+coord_esep(cube_t c)
+{
+        int64_t i, j, jj, k, l, ret1, ret2, bit1, bit2, is1;
+
+        for (i = 0, j = 0, k = 4, l = 4, ret1 = 0, ret2 = 0; i < 12; i++) {
+                /* Simple version:
+                if (c.edge[i] & _esepbit2) {
+                        ret1 += binomial[11-i][k--];
+                } else {
+                        if (c.edge[i] & _esepbit1)
+                                ret2 += binomial[7-j][l--];
+                        j++;
+                }
+                */
+
+                bit1 = (c.edge[i] & _esepbit1) >> 2;
+                bit2 = (c.edge[i] & _esepbit2) >> 3;
+                is1 = (1 - bit2) * bit1;
+
+                ret1 += bit2 * binomial[11-i][k];
+                k -= bit2;
+
+                jj = j < 8;
+                ret2 += jj * is1 * binomial[7-(j*jj)][l];
+                l -= is1;
+                j += (1-bit2);
+        }
+
+        return ret1 * 70 + ret2;
+}
+
+_static_inline void
+copy_corners(cube_t *dest, cube_t src)
+{
+        memcpy(&dest->corner, src.corner, sizeof(src.corner));
+}
+
+_static_inline void
+copy_edges(cube_t *dest, cube_t src)
+{
+        memcpy(&dest->edge, src.edge, sizeof(src.edge));
+}
+
+_static_inline void
+set_eo(cube_t *cube, int64_t eo)
+{
+        uint8_t i, sum, flip;
+
+        for (sum = 0, i = 1; i < 12; i++, eo >>= 1) {
+                flip = eo % 2;
+                sum += flip;
+                cube->edge[i] = (cube->edge[i] & ~_eobit) | (_eobit * flip);
+        }
+        cube->edge[0] = (cube->edge[0] & ~_eobit) | (_eobit * (sum % 2));
+}
+
+_static_inline cube_t
+invcoord_esep(int64_t esep)
+{
+        cube_t ret;
+        int64_t bit1, bit2, i, j, jj, k, l, s, v, w, is1, set1, set2;
+        uint8_t slice[3] = {0};
+
+        ret = solved;
+        set1 = esep % 70;
+        set2 = esep / 70;
+
+        for (i = 0, j = 0, k = 4, l = 4; i < 12; i++) {
+                v = binomial[11-i][k];
+                jj = j < 8;
+                w = jj * binomial[7-(j*jj)][l];
+                bit2 = set2 >= v;
+                bit1 = set1 >= w;
+                is1 = (1 - bit2) * bit1;
+
+                set2 -= bit2 * v;
+                k -= bit2;
+                set1 -= is1 * w;
+                l -= is1;
+                j += (1-bit2);
+                s = 2*bit2 + (1-bit2)*bit1;
+
+                ret.edge[i] = (slice[s]++) | (uint8_t)(s << 2);
+        }
+
+        return ret;
+}