Thanks for a fascinating puzzle! Yes, it looks like we can do better than a conversion through cartesian coordinates of hexagon centers. It can be done entirely with integer math, though I've included a rational in a matrix below to keep the notation concise.
You're right that both the encoding and decoding processes require loops. Fortunately, because SHM works by orders of magnitude, the loops will never need more than Log_7(n + 1) iterations, where n is the number of tiles in the grid, which is very reasonable.
To keep my notation symmetric, I've opted to use a 3-coordinate grid like the one below (what your example from Red Blob Games calls "Cube Coordinates"). If you prefer, you can omit the z coordinate and replace all instances of .z with -(x + y) to get "Axial Coordinates" instead.
For speed and cleanliness of code, I recommend using lookup tables for encoding and decoding. First, for decoding SHM into cube coordinates, you can build the table with these recurrences:
Rotate(v) = (-v.z, -v.x, -v.y)
Decode[,] = 2D array of ordered triples, [7 x (max_order_of_magnitude + 1)]
Decode[0, k] = (0, 0, 0) for all k = 0...max_order_of_magnitude
Decode[1, 0] = (0, -1, 1)
Decode[d + 1, k] = Rotate(Decode(d, k)) for d = 2...6, all k
Decode[d, k + 1] = Decode(d, k) + 2 * Rotate(Decode(d, k)) for all d, k
Then decoding becomes extremely simple:
point = (0, 0, 0)
order = 0
while(code > 0)
digit = code % 7
point += Decode(digit, order)
code = floor(code/7)
This works because the first row of the table stores the offsets of the cells encoded by the digits 0, 1, 2, 3, 4, 5, 6. The next row stores the offsets for 0, 10, 20, 30, 40, 50, 60, and so on. By summing the offsets contributed by each order of magnitude, you can decode any point within the range of the table.
For encoding, I think the following method will work, but I'll confess I haven't tested it. It uses a much simpler lookup table:
Encode = [0, 5, 1, 6, 3, 4, 2]
and a rotation & scaling matrix:
Untwist = | 5 -4 2 |
| 2 5 -4 |
|-4 2 5 | * 1/21
To encode a cube coordinate into SHM, then...
code = 0
magnitude = 1
while(point != (0, 0, 0))
digitIndex = (((point.x - 2 * point.y) % 7) + 7) % 7
code += magnitude *Encode[digitIndex]
point = Untwist * (point - Decode[digit, 0])
// treating "point" as a column vector, and rounding if necessary
magnitude *= 7
Note that the code generated here (and decoded above) is in standard numeric format (ie. the cell marked "10" in the SHM diagram will be assigned a code of 7 in base 10, or 111 in binary). If you want to display it in base 7, you'll need to use the appropriate radix when converting it to a string.
This works by figuring out the least significant digit of the SHM code at each iteration, (using the fact that these digits repeat in a way that tiles the plane regularly - the LUT is the sequence of digits read from (0, 0, 0) to (6, 0, -6), which repeats at an offset every row). Next it effectively zeroes that digit by subtracting the corresponding offset from our Decode table earlier, and scales the fractal down one level using the Untwist matrix. This acts like a right shift in base-7, in geometric form - it brings the next digit into the least significant position, bringing cell "10" to "1" and so on, where we can use the same trick again to read that digit on the next loop.
I hope that works for you!