FastMortonToHilbert explained

fast_hilbert_to_morton.cpp

uint32_t fastMortonToHilbert3_32(const uint32_t morton, const unsigned iteration) noexcept
{
    // one specific permutation table for one specific case of morton <-> hilbert has very interesting properties
    // it consists of 3 permutations (rot 1, xor 3, (rot 1, xor 6)) as well as their inverse
    // rotations are all to the right
    // the first column of this permutation table encodes the xor component (how convenient)
    // the rotations simply have to be calculated but are always 0, 1, or 2 as a rotation of 3 is equiv. to 0

    constexpr uint32_t permTableXors = 0x5356'0300;  // first digit of the 8 permutations, also encodes XOR component
    constexpr uint32_t permTableRots = 0x2021'2021;  // right-shifts of each permutation
    constexpr uint32_t mask = 0x9249'2492;

    uint32_t hilbert = morton;

    uint32_t block = iteration * 3;
    uint32_t hcode = (hilbert >> block) & 0b111;
    uint32_t totalRot = 0, totalXor = 0;
    while (block) {
        block -= 3;
        hcode <<= 2;                  // leftshift previous octal digit so that we can use it in lookup tables
                                      // when rightshifting in the next steps, this now means that we implicitly
                                      // rightshift by 4 * (original hcode)
                                      // leftshifting by 1 here would not have been enough because digits would overlap

        uint32_t permRot = (permTableRots >> hcode) & 0b11;
        // compute the new total rotation with the help of a magic constant
        // this is a shortcut for (totalRot + permShift) % 3
        totalRot = (0b1001 >> (4 - totalRot - permRot)) & 0b11;

        // we must also apply our rotation to our total xor
        totalXor = rightRotOctal(totalXor, permRot);
        totalXor ^= permTableXors >> hcode;
        totalXor &= 0b111;               // convert back to octal digit

        uint32_t mcode = (morton >> block) & 0b111;

        hcode = mcode;
        hcode = rightRotOctal(hcode, totalRot);
        hcode ^= totalXor;
        hcode &= 0b111;               // convert back to octal digit

        hilbert ^= (mcode ^ hcode) << block; // replace current digit
    }

    // these last steps encode morton digits as hilbert digits in parallel in a snake-curve pattern
    // consider the bit-order zyx
    hilbert ^= (hilbert >> 1) & mask; // for each digit flip y if z is set
    hilbert ^= (hilbert & mask) >> 1; // for each digit flip x if y is set
    return hilbert;
}