Convert 64-bit Gray code to integer

Question

This algorithm encodes a 64 bit integer. As input I have a 64 bit integer a, for example:

(a is a single 64 bit number,which I split in 8 bytes for readability) a=11110110,10101110,00010000,01100110,11101011,11000001,11001011,01100011

To calculate b=f(a), I begin with the rightmost digit (the less significant digit) and write '0'. Then, I work my way leftward. Each time I encounter a '1' in a, I switch between writing '0' and '1'. For example, if a=1011, then b=0010.

Here is an example of 64 bit a and b=f(a)

a=11110110,10101110,00010000,01100110,11101011,11000001,11001011,01100011
b=10100100,11001011,11100000,01000100,10110010,10000001,01110010,01000010

This is dumb code that calculates it for a list of bytes:

# Slow function to memoize f(numbers_to_memoize) up to 8 bits
# Returns the encoded list
def encode_reversed_gray_code(numbers: list) -> list:
    # Convert to binary
    numbers_bin = np.vectorize(np.binary_repr)(numbers, width=bit_length)

    # Convert each string to a list of digits
    numbers_bin_list = [list(x) for x in numbers_bin]

    # encode reversed gray code
    gray_code = []
    for _list in numbers_bin_list:
        reversed_number = _list[::-1]

        next_digit_is_1 = False
        encoded_list_reversed = ["0"]
        for x in reversed_number:
            if x == "1":
                next_digit_is_1 = not next_digit_is_1
            encoded_list_reversed.append(["0", "1"][next_digit_is_1])

        gray_code.append(int("".join(encoded_list_reversed[-2::-1]), 2))

    # print("".join(new_list[::-1]))
    return np.asarray(gray_code).astype(np.uint8)

That code is slow, but runs only once. I use it to accelerate the calculation for 64 bit.

I need a lot more speed, so I thought to memoize the result in a table, but the table would have 2⁶⁴ uint64, which is too much memory.

So, I thought to memoize only up to 8 bits, which only requires 256 bytes of memory. And build the 64 bit code as a paste of 8 bytes encoded. This will allow me to use the memoization tables for 8-bit integers to encode larger integers.

With the above code, I made an 8-bit table named gray_code_starting_in_0

The 8-bit table is a convenient approach, because most calculations will involve much smaller numbers than 64 bit. So, the small table will likely be useful on its own, requiring only one or two lookups, and I guess will fit easily into the cpu cache, meanwhile an 16 bit will probably not.

The caveat is that although the encoding of the entire 64-bit uint starts with 0, the individual bytes may not. A new byte encoding does not necessarily start with 0, but with the most significant bit of the last byte XORed. If the last byte was byte_a, and its encoding was code(byte_a), the next byte starts with (byte_a ^ code(byte_a)) >> 7.

So I made a second table memoizing the same encoding, but starting with "1" instead of "0": gray_code_starting_in_0 = not gray_code_starting_in_0

this is the complete code:

import numpy as np


bit_length = 8
numbers_to_memoize = np.arange(0, 2**bit_length).astype(np.uint8)


# Slow function to memoize f(numbers_to_memoize) up to 8 bits
# Returns the encoded list
def encode_reversed_gray_code(numbers: list) -> (list, list):
    # Convert to binary
    numbers_bin = np.vectorize(np.binary_repr)(numbers, width=bit_length)

    # Convert each string to a list of digits
    numbers_bin_list = [list(x) for x in numbers_bin]

    # encode reversed gray code
    gray_code = []
    for _list in numbers_bin_list:
        reversed_number = _list[::-1]

        next_digit_is_1 = False
        encoded_list_reversed = ["0"]
        for x in reversed_number:
            if x == "1":
                next_digit_is_1 = not next_digit_is_1
            encoded_list_reversed.append(["0", "1"][next_digit_is_1])

        gray_code.append(int("".join(encoded_list_reversed[-2::-1]), 2))

    # print("".join(new_list[::-1]))
    return np.asarray(gray_code).astype(np.uint8)


gray_code_starting_in_0 = encode_reversed_gray_code(numbers_to_memoize)


# the bitwise_not, calculates the code as if it started with 1 instead of 0
gray_code_starting_in_1 = np.vectorize(np.bitwise_not)(gray_code_starting_in_0).astype(
    np.uint8
)

# Dummy variable to be overwritten, to avoid creating a new one each time a code is caculated
chunks_encoded = np.asarray([0] * 8, dtype=np.uint8)
# Weights to reconvert the array of bytes into a 64 bit uint
chunk_weights = 2 ** np.arange(0, 8**2, 8, dtype=np.uint64)

FF = 0xFF

# code to split a 64-bit integer into bytes
import struct

pack = struct.pack
unpack = struct.unpack
# Usage:
# 8_byte_chunks_of_n = unpack("8B", pack("Q", n))

from math import ceil


def fast_64_bit_encoding(n: int) -> int:
    # split n in byte sized chunks
    n = int(n)

    # split 64 bit integer into 8 bit chunks
    number_of_chunks = ceil(n.bit_length() / 8)
    chunks = unpack("8B", pack("Q", n))[:number_of_chunks]

    next_starting_bit_is_1 = False
    for index, chunk in enumerate(chunks):
        # print(",".join(np.vectorize(np.binary_repr)(chunks[::-1], width=bit_length)))
        # print(
        #     ",".join(
        #         np.vectorize(np.binary_repr)(chunks_encoded[::-1], width=bit_length)
        #     )
        # )
        # print()
        if next_starting_bit_is_1:
            chunks_encoded[index] = gray_code_starting_in_1[chunk]
        else:
            chunks_encoded[index] = gray_code_starting_in_0[chunk]

        next_starting_bit_is_1 = (chunk ^ chunks_encoded[index]) >> 7 == 1

        # # assert that next starting bit is correct
        # assert (
        #     ((chunk ^ chunks_encoded[index]) & int("10000000", 2)) > 0
        # ) == next_starting_bit_is_1, "".join(
        #     [
        #         f""" chunk {np.binary_repr(chunk, width=bit_length)}\n""",
        #         f"""code:  {np.binary_repr(chunks_encoded[index], width=bit_length)}\n""",
        #         f"""next bit: { next_starting_bit_is_1}""",
        #     ]
        # )

    return np.dot(chunks_encoded, chunk_weights)


# Random uint64 to test the encoding
a = int("1" + "".join(np.random.choice(["0", "1"], size=63)), 2)
answer = f"{np.binary_repr(fast_64_bit_encoding(a))}"

# split "a" and "answer" in bytes
bin_a = bin(a)[2:]
a_split = [bin_a[i : i + 8] for i in range(0, len(bin_a), 8)]
answer = [answer[i : i + 8] for i in range(0, len(answer), 8)]

print("a=" + ",".join(a_split))
print("b=" + ",".join(answer))

Note: the OEIS sequence A006068 calculates the same thing, but from left to right. That webpage has slow code (but simpler) for that calculation:

def oeis_A006068(n: np.uint64):
    s = 1
    while True:
        ns = n >> s
        if ns == 0:
            break
        n = n ^ ns
        s <<= 1
    return n

Answer 1

While it's not your code, I'll point out that the routine from OEIS can be significantly shortened by using an assignment expression, available since Python 3.8:

def oeis_A006068(n: int):
    s = 1
    while ns := n >> s:
        n ^= ns
        s <<= 1
    return n

Here's an unrolled version of it for 64-bit arguments:

def oeis_A006068_64bit(n: np.uint64):
    n ^= n >> 1
    n ^= n >> 2
    n ^= n >> 4
    n ^= n >> 8
    n ^= n >> 16
    n ^= n >> 32
    return n

Either of these will probably be significantly faster than your approach, since they involve relatively few basic operations (bytecode instructions).

---

Converting to the Gray-code representation is simpler than converting back. In your code, you do in effect

for g in range(256):
    table[g] = from_gray(g)

you could simplify the construction of the table by instead doing this:

for i in range(256):
    table[to_gray(i)] = i

or more precisely

def make_lookup_table():
    table = np.empty(1 << bit_length, dtype=np.uint8)

    for i in range(1 << bit_length):
        table[i ^ (i >> 1)] = i

    return table


gray_code_starting_in_0 = make_lookup_table()

(That doesn't produce the same output as your code, but I think your code is buggy: for one thing, it only produces even numbers.)

---

numbers_bin_list = [list(x) for x in numbers_bin]

It's unnecessary to convert the strings to lists. Iterating over the strings directly will return the same characters.

for _list in numbers_bin_list:

Leading underscores in Python are usually reserved for variables that are private in some sense. I would call the variable list_ or lst or (probably better) a descriptive name like number_bin, or even a short loop-variable name like n.

# Dummy variable to be overwritten, to avoid creating a new one each time a code is caculated
chunks_encoded = np.asarray([0] * 8, dtype=np.uint8)

Using a global variable as scratch space is generally a bad idea since it makes the code not thread safe.

FF = 0xFF

I don't think this is a good idea. Either use 0xFF directly or give the constant a descriptive name like BYTE_MASK. (In any case, it isn't used.)

import struct

pack = struct.pack
unpack = struct.unpack

You can just say from struct import pack, unpack instead. Also, imports should generally be at the top of the file.

chunks = unpack("8B", pack("Q", n))[:number_of_chunks]

The unpack isn't necessary; you can slice and iterate over the return value of pack.

The byte order of pack("Q", ...) is platform-dependent. If you want little-endian output you should say pack("<Q", ...).

There's a built-in method of int that would work better here: n.to_bytes(number_of_chunks, 'little'). That way you needn't import the struct module at all, and it will work if number_of_chunks > 8.

return np.dot(chunks_encoded, chunk_weights)

Similarly, this could be return int.from_bytes(chunks_encoded, 'little').

Answer 2

I assume that your first example is depicted in human-legible (i.e. big-endian) order. The right-most digit is indeed then the least-significant digit. When you call binary_repr, that converts into a big-endian string, which is not particularly useful for this application. Instead you should be looking at something like unpackbits with bitorder='little', which will produce a little-endian array. You should be working with little-endian arrays (and not lists), based on your description of processing order.

I need a lot more speed

How much? How slow is your current method? Where is it slow? You haven't profiled, and you really need to.

For true vectorisation it isn't enough to call np.vectorize; the entire inner loop needs to be rewritten. The inner loop should assume that both the input and output are little-endian uint8 arrays of 64 bits, like:

def encode_reversed_gray_code(numbers: np.ndarray) -> np.ndarray:
    increasing = numbers[:-1].cumsum(dtype=np.uint8)
    output = np.zeros(shape=numbers.size, dtype=np.uint8)
    np.bitwise_and(increasing, 1, out=output[1:], dtype=np.uint8)
    return output

An XOR version is also possible:

def encode_reversed_gray_code(numbers: np.ndarray) -> np.ndarray:
    output = np.zeros(shape=numbers.size, dtype=np.uint8)
    np.bitwise_xor.accumulate(numbers[:-1], out=output[1:], dtype=np.uint8)
    return output

with demo code

a_in = np.array(
    (
        1,1,0,0,0,1,1,0,
        1,1,0,1,0,0,1,1,
        1,0,0,0,0,0,1,1,
        1,1,0,1,0,1,1,1,
        0,1,1,0,0,1,1,0,
        0,0,0,0,1,0,0,0,
        0,1,1,1,0,1,0,1,
        0,1,1,0,1,1,1,1,
    ),
    dtype=np.uint8,
)
a_out = encode_reversed_gray_code(a_in)
print(a_in.reshape((-1, 8)))
print(a_out.reshape((-1, 8)))

and output

[[1 1 0 0 0 1 1 0]
 [1 1 0 1 0 0 1 1]
 [1 0 0 0 0 0 1 1]
 [1 1 0 1 0 1 1 1]
 [0 1 1 0 0 1 1 0]
 [0 0 0 0 1 0 0 0]
 [0 1 1 1 0 1 0 1]
 [0 1 1 0 1 1 1 1]]
[[0 1 0 0 0 0 1 0]
 [0 1 0 0 1 1 1 0]
 [1 0 0 0 0 0 0 1]
 [0 1 0 0 1 1 0 1]
 [0 0 1 0 0 0 1 0]
 [0 0 0 0 0 1 1 1]
 [1 1 0 1 0 0 1 1]
 [0 0 1 0 0 1 0 1]]

This can still be memoized. Depending on a few things, it would be easy to increase the memoisation width to 16, which will still fit in memory (assuming a typical desktop computer). Though this is perhaps splitting hairs, what you've shown is not really traditional memoisation. Memoisation is the process of computing only the requested output and caching it for reuse. You've pre-computed a complete LUT of all possible inputs and outputs. That's maybe fine, and may or may not be justifiable (depending on the speed of the newly vectorised gray operation). Construction of the LUT can itself be vectorised:

t_lut = np.uint16
lut_bytes_in = (
    np.arange(1 + np.iinfo(t_lut).max, dtype=t_lut)
    .view(dtype=np.uint8)
    .reshape((-1, t_lut().nbytes))
)
lut_bits_in = np.unpackbits(lut_bytes_in, axis=1, bitorder='little')
lut_bits_out = np.zeros_like(lut_bits_in)
np.bitwise_xor.accumulate(
    array=lut_bits_in[:, :-1],
    out=lut_bits_out[:, 1:],
    axis=1,
    dtype=np.uint8,
)