kuniga.me > NP-Incompleteness > Asymmetric Numeral Systems

Asymmetric Numeral Systems

20 Aug 2026

Jarosław Duda is a Polish professor at the Jagiellonian University in Kraków. He developed a family of entropy coding methods called asymmetric numeral systems (ANS), mainly used in data compression. He wanted these to remain patent-free but has had mixed success.

Google worked with Duda around 2014 in a paper Mixed boolean-token ANS coefficient coding and tried to patent a coder for video, but Duda pushed back and Google abandoned the attempt. In 2019 Microsoft was able to patent a variant and since then other patents have been granted internationally.

A variant of ANS, known as FSE, is used by Meta’s compression library called zstd. I wanted to learn more about it and decided to study this algorithm first.

Setup

We’ll develop the Python code as we follow along, both to make the explanations more precise and to define common boilerplate gradually so each code snippet stays shorter. Since we’ll be developing and comparing encoders, we define an abstract class:

class Encoder(ABC):
    @abstractmethod
    def desc(self) -> str:
        "Human-friendly description of this encoder"

    @abstractmethod
    def encode(self, input: str) -> bytes:
        "Gets an input, encodes to bytes"

    @abstractmethod
    def decode(self, encoded: bytes) -> str:
        "Decodes bytes to the string"

A dummy encoder is simply converting strings to bytes using UTF-8. Not very useful but it gives us a working implementation:

class UTFEncoder(Encoder):
    def desc(self) -> str:
        return "UTF-8"

    def encode(self, input: str) -> bytes:
        return input.encode('utf-8')

    def decode(self, encoded: bytes) -> str:
        return encoded.decode("utf-8")

A simple way to evaluate the encoders is to encode an input, decode it and verify it matches the original input since we’re assuming lossless encoders. We can also compute the compression rate as (1 - encoded/decoded) and the time it takes to encode and decode.

encoders = [
    UTFEncoder(),
]

decoded_sz = len(data)
for encoder in encoders:
    start = time.perf_counter()

    encoded = encoder.encode(data)
    encoded_sz = len(encoded)
    rate = (1 - encoded_sz/decoded_sz) * 100

    output = encoder.decode(encoded)
    assert data == output

    elapsed = time.perf_counter() - start

    rows.append([
        encoder.desc(),
        elapsed,
        rate,
    ])

print(tabulate(
    rows,
    headers=["Encoder", "Time (s)", "Compression Rate"],
    tablefmt="rounded_grid",
    floatfmt=".3f",
))

Huffman Coding

We wrote about Huffman coding before. The idea is: given a probability distribution of symbols, we construct the Huffman tree which is a binary tree where symbols are at the leaves. Then we build a map from each symbol to a binary code representing the path to the corresponding leaf (0 is left, 1 is right).

We can then encode an input string character by character as a byte string. This encoding has the property of being prefix-free, which allows us to decode the byte string greedily by using the 0s and 1s to traverse the Huffman tree.

The Huffman encoding minimizes this function:

where $p_i$ is the probability of symbol $i$ and $l_i$ is the length of the path from the root to the leaf containing $i$. We then talked about entropy which is defined as:

and that while the Huffman encoding doesn’t minimize it, it’s within $H + 1$.

Intuitively the gap is that $l_i$ is an integer, while $\log(p_i)$ is a real value. This is where Asymmetric Numeral Systems come in. Before we go there, let’s cover a more familiar concept but with an unfamiliar name: symmetric numeral systems.

We won’t cover the code for Huffman again, but we can assume a class HuffmanEncoder exists extending Encoder.

Symmetric Numeral Systems

Consider how we convert a binary string $b_n b_{n-1} b_{n-2} \cdots b_0$ (most significant bit first) to decimal. In Python it can look like:

def btoi(bits):
    r = 0
    for b in bits:
        r = r * 2 + int(b)
    return r

Note how we use one bit per symbol 0 or 1 and note how the integer roughly doubles each time. Another way to implement this is via:

def encode_bit(x, b):
    match b:
        case '0':
            return 2 * x
        case '1':
            return 2 * x + 1

def btoi(bits):
    r = 0
    for b in bits:
        r = encode_bit(r, b)
    return r

One way to interpret encode_bit() is as an explicit map where 0 gets the even integers and 1 gets the odd ones. Then we create a mapping from the natural numbers to their respective subset.

    0  1  2  3  4  5 ...
0:  0  2  4  6  8 10 ...
1:  1  3  5  7  9 11 ...

Another way to see this is as a partition of the natural numbers into two subsets and the association of one to 0 and another to 1. As we know, this function is reversible so we can convert an integer to a string of bits (if we know how many bits we started with).

A different way to see this is as a reversible operation is to look at which partition of the naturals the integer we’re decoding is at. If it’s an even number, it means it was last mapped by a 0, whereas if it’s odd, it was mapped by a 1. This interpretation will be more useful soon.

We can generalize this encoding for an arbitrary base, and thus we can also encode an input that used any alphabet. For example, if we have 3 symbols and the input ACBAABCBA, we can encode it as an integer by converting this as if it was a string of digits in base 3.

Asymmetric Numeral Systems

Suppose now that 0s are 3x more likely than 1s to appear in the input. Intuitively we want to use fewer bits to encode 0 than 1. We can now distribute the image non-uniformly: for every 4 numbers we assign the first 3 to 0 and the fourth to 1:

    0  1  2  3  4  5 ...
0:  0  1  2  4  5  6 ...
1:  3  7 11 15 19 23 ...

In the original example, we said each input maps to roughly 2x its value. Now for 1 the input grows at 4x rate, but for 0 it’s more like 4/3x, which is closer to their distribution!

We can represent this distribution in encode_bit() via:

def encode_bit(x, b):
    match b:
        case '0':
            return 4 * (x // 3) + (x % 3)
        case '1':
            return 4 * x + 3

This numeral system, in which different digits have different weights, is called Asymmetric Numeral Systems or ANS. As in the symmetric case, we can think of a general base $n$ system in which digit/symbol $i$ occurs with probability $p_i$.

We find a suitable denominator and then the rational that is the closest approximation to the probability. For example, if we have 3 symbols (A, B, C) with probabilities $0.5$, $0.3$ and $0.2$ respectively and we choose a denominator $M = 16$, we can have $A = 8 / 16$, $B = 5 / 16$ and $C = 3 / 16$.

We can then distribute a chunk of 16 integers between them: the first 8 goes to A, the next 5 to B and the last 3 to C:

0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
A  A  A  A  A  A  A  A  B  B  B  B  B  C  C  C

Let $M$ be the chunk size, $f_i$ how many slots of $M$ we assign to symbol $i$, and $o_i$ the offset in the chunk where $i$ starts. So in the example above we have: $f_A = 8, f_B = 5, f_C = 3$, $o_A = 0, o_B = 8, o_C = 13$.

To find to which integer an input $x$ maps for a symbol $i$, we can first find the right chunk by doing $\lfloor{x / f_i\rfloor}$. So for $A$, 0-7 will map to the 1st chunk, 8-15 to the 2nd, etc. Then we compute the offset of that chunk by multiplying by $M$. Now we add the offset of $i$ within that chunk, $o_i$ and finally the offset of the $k$-th copy of $i$ in that segment, $k = x \pmod {f_i}$. The formula is a bit convoluted but you can try with a few examples to convince yourself.

Now to the implementation. We’ll first define the class SimpleANSEncoder extending Encoder, which computes $f$ (as freq) and $o$ (as offsets):

FrequencyTable = Sequence[tuple[str, int]]

class SimpleANSEncoder(Encoder):

    def __init__(self, raw_freq: FrequencyTable) -> None:
        self.M = 2**12

        # least frequent first
        sorted_freq = sorted(raw_freq, key=lambda item: item[1])
        total = sum(cnt for _, cnt in sorted_freq)

        freq_norm = [
            (s, max(int(cnt/total * self.M), 1)) for s, cnt in sorted_freq
        ]

        self.offsets = {}
        acc = 0
        for s, f in freq_norm:
            self.offsets[s] = acc
            acc += f

        # tweak last symbol freq so it adds up to M
        last = freq_norm[-1]
        freq_norm[-1] = (last[0], last[1] + self.M - acc)

        self.freq = dict(freq_norm)

The code is a bit complicated because we want to make sure that sum(prob) == M, while we’re rounding floats to integers. A trick is to have the most frequent symbol have M - the sum of the probabilities of the least frequent items.

For incorporating a symbol into the integer x we follow the formula above:

def encode_symbol(self, x: int, symbol: str) -> int:
    f = self.freq[symbol]
    off = self.offsets[symbol]
    return self.M * (x // f) + (x % f) + off

The encode() consists in building the giant integer and encoding it as a byte string. We can use BytesIO to read/write to this byte string:

def encode_int(cursor: BytesIO, x: int) -> bytes:
    cursor.write(x.to_bytes(4, byteorder="big"))

def encode(self, input: str) -> bytes:
    r = 0
    for c in input:
        r = self.encode_symbol(r, c)

    cursor = BytesIO()
    encode_int(cursor, len(input))
    encode_bigint(cursor, r)
    return cursor.getvalue()

We also include the length of the original text because otherwise the decoder doesn’t know how many symbols to generate back. For example, in the binary to decimal encoding, leading zeros would be lost during encoding, so 0001 gets encoded as 1 and decoded as 1.

The decoding is only a bit more complicated. Finding the symbol is easy: we just need to mod $M$, since this gives the offset within a chunk and then we can do a binary search on off to find to which offset it belongs. Alternatively we store an array of size $M$ with the symbols repeated, which we can compute upfront, in __init__(). In our example:

class SimpleANSEncoder(Encoder):
    def __init__(self, raw_freq: FrequencyTable) -> None:
        # ...
        self.lookup = "".join(s * self.freq[s] for s, _ in sorted_freq)

    def get_symbol(self, x):
        return self.lookup[x % self.M]

The decode_symbol() is roughly the inverse of encode_symbol once we have the symbol:

def decode_symbol(self, x: int) -> tuple[str, int]:
    symbol = self.get_symbol(x)
    p = self.freq[symbol]
    off = self.offsets[symbol]
    return symbol, p * (x // self.M) + (x % self.M) - off

And so is the decode():

def decode_int(cursor: BytesIO) -> int:
    return int.from_bytes(cursor.read(4), byteorder="big")

def decode(self, encoded: bytes) -> str:
    cursor = BytesIO(encoded)
    sz = decode_int(cursor)
    value = decode_bigint(cursor)

    output = ''
    for _ in range(sz):
        symbol, value = self.decode_symbol(value)
        output += symbol
    return output[::-1]

Efficiency

The rational approximation probability of a symbol is $q_i = f_i / M$. The rough growth of x when going through encode() is roughly dominated by $M / f_i$, so $1 / q_i$. The approximate number of bits in this factor is $\log_2 1/q_i = - \log_2 q_i$. So the expected cost per symbol is a weighted sum:

Where $L$ is the cost of a random variable representing a symbol. This expected value is also equal to the cross-entropy which indicates how many bits on average we need if the symbols come from a distribution $P$, but we encode it using a distribution $Q$ and is denoted by $H(P, Q)$.

We can compare this with the true entropy:

which we’ve seen in Huffman coding is the theoretical optimal. So the closer we can approximate the rational $q_i$ to $p_i$ the better, so in theory we could use a gigantic value for $M$ but in practice working with large values is prohibitive.

Speaking of large integers, encoding an input as an integer has the major downside of the integer growing exponentially with the size of the input.

Renormalization

One way to avoid working with arbitrarily large integers is to move the lowest bits to a bit stream whenever it grows too big. It will be part of the encoded output but is not going to participate in future multiplications. The idea is to have the invariant $M \le x \lt 2M$. If multiplying $x$ by a factor would tip it over $2M$ we first reduce it.

We define a class RANSEncoder extending from SimpleANSEncoder since they share a lot of the methods. For encoding, encode_symbol() is the same, but now we need to emit bits periodically. We use bitarray for that:

def encode_bitarray(cursor: BytesIO, x: bitarray) -> None:
    encode_int(cursor, len(x))
    cursor.write(x.tobytes())

class RANSEncoder(SimpleANSEncoder):

    # ...

    def encode(self, input: str) -> bytes:
        stream = bitarray()
        x = self.M
        for b in input:
            while (y := self.encode_symbol(x, b)) >= 2 * self.M:
                stream.append(x & 1)
                x >>= 1
            x = y

        cursor = BytesIO()
        encode_int(cursor, len(input))
        encode_bitarray(cursor, stream)
        encode_bigint(cursor, x)
        return cursor.getvalue()

Note how we start $x$ as $M$ as opposed to $0$ as we did before, so that we maintain the invariant $M \le x \lt 2M$ inside the loop.

For decoding, we can also maintain the invariant and consume bits from the stream to restore it:

def decode_bitarray(cursor: BytesIO):
    x = bitarray()
    bit_cnt = decode_int(cursor.read(4))
    byte_cnt = (bit_cnt + 7) // 8
    x.frombytes(cursor.read(byte_cnt))
    x = x[:bit_cnt]  # remove padding
    return x

class RANSEncoder(SimpleANSEncoder):

    # ...

    def decode(self, encoded: bytes) -> str:
        cursor = BytesIO(encoded)
        input_size = decode_int(cursor.read(4))
        stream = decode_bitarray(cursor)

        out = ''
        x = decode_int(cursor.read())
        for _ in range(input_size):
            s, x = self.decode_symbol(x)
            while x < self.M:
                bit = stream.pop()
                x = (x << 1) | bit
            out += s

        return out[::-1]

Not only does this allow us to operate with a small integer range, the set of values we pass for x to encode_symbol() is bounded between [M, 2M[, which is very useful as we’ll see next.

Tabular ANS

Recall that in our 3-symbol example with distinct weights we laid them out like this:

0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
A  A  A  A  A  A  A  A  B  B  B  B  B  C  C  C

First the A’s, then the B’s and then the C’s. The problem with this approach is that C always gets penalized. For example, if the input is, say 1, and the symbol is C, it will be mapped to 14, while B gets mapped to 9 and A gets mapped to 1. So C grows a lot faster.

It will be more balanced if we interleave the symbols more evenly within one batch, for example:

0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
A  B  A  C  A  B  A  B  A  C  A  B  A  B  A  C

Now for input 1, A maps to 2, B to 5 and C to 9. The problem with this arbitrary order is that it becomes much harder to devise a formula that maps integers to their corresponding output. That’s where the renormalization helps: if we have a fixed range of values of input ([M, 2M[), we can write a table with M entries that map them exactly to the order above! This variant is called the Tabular ANS.

Let’s create the TRANSEncoder also extending from SimpleANSEncoder. The main difference on the initialization is that we’ll “shuffle” the lookup array and then pre-compute the encode and decode table.

class TRANSEncoder(SimpleANSEncoder):
    def __init__(self, raw_freq: FrequencyTable) -> None:
        super().__init__(raw_freq)

        step = (self.M >> 1) + (self.M >> 3) + 3
        lookup: list[str] = [''] * self.M

        pos = 0
        for symbol, cnt in self.freq.items():
            for _ in range(cnt):
                lookup[pos] = symbol
                pos = (pos + step) % self.M

        self.lookup = lookup
        self.precompute_encode()
        self.precompute_decode()

The method precompute_encode() is analogous to RANSEncoder’s encode() but it uses the positions computed above. It returns a table that for each M <= x < 2*M and symbol s, stores a pair (x, bits), where x is the new value of x after incorporating s and the bits that will go into the stream.

def precompute_encode(self) -> None:
    positions = defaultdict(list)
    for p, s in enumerate(self.lookup):
        positions[s].append(p)

    def encode_symbol(x: int, s: str) -> int:
        slots = positions[s]
        f = len(slots)
        return self.M * (x // f) + slots[x % f]

    self.encode_table: list[dict[str, EncodeEntry]] = [
        {} for _ in range(self.M)
    ]
    for x0 in range(self.M, 2 * self.M):
        for s in self.freq:
            x = x0

            slots = positions[s]
            f = len(slots)
            bits = bitarray()
            while (y := encode_symbol(x, s)) >= 2 * self.M:
                bits.append(x & 1)
                x >>= 1

            self.encode_table[x0 - self.M][s] = (
                y,
                bits,
            )

With this table, encoding becomes very simple (omitting the writes):

def encode(self, input: str) -> bytes:
    stream = bitarray()
    x = self.M
    for b in input:
        x, bits = self.encode_table[x - self.M][b]
        stream += bits
    # ...

The decode table only depends on M <= x < 2*M and it maps to a pair: the previous state and how many bits it needs to consume from the stream to grow x back to the range.

def precompute_decode(self) -> None:
    cnt = defaultdict(int)
    decode_table = []
    for s in self.lookup:
        rank = cnt[s]
        cnt[s] += 1

        prev_x = self.freq[s] + rank
        bit_cnt = 0
        while (prev_x << bit_cnt) < self.M:
            bit_cnt += 1

        prev_x <<= bit_cnt

        decode_table.append((
            prev_x,
            bit_cnt
        ))
    self.decode_table  = decode_table

Decoding becomes:

def decode(self, encoded: bytes) -> str:
    cursor = BytesIO(encoded)
    input_size = decode_int(cursor)
    stream = decode_bitarray(cursor)

    out = ''
    x = decode_bigint(cursor)
    for _ in range(input_size):
        off = x - self.M
        s = self.get_symbol(off)
        base, bit_cnt = self.decode_table[off]

        suffix = 0
        for _ in range(bit_cnt):
            suffix = (suffix << 1) | stream.pop()

        x = base | suffix
        out += s

    return out[::-1]

Experiments

I used Carroll’s Alice in Wonderland as a training dataset to estimate the character frequencies, then encoded Shakespeare’s As You Like It (125k characters):

As expected, the no-op UTF-8 is pretty efficient! But it doesn’t compress anything. Huffman is very fast and is only slightly worse compression than ANS. Working with gigantic integers is very slow.

tANS was slower than rANS mostly because the overhead of constructing the tables wasn’t enough to pay off for 125k characters. To verify that, I also tried encoding a bigger text ($10^8$ characters):

Here the overhead paid off for tANS. It also achieved a better compression rate with the more spread out symbols in the chunk.

Conclusion

I found it hard to grok ANS at first until realizing it is, in a way, a generalization of the binary to decimal algorithm!

I found it difficult to implement these in Python, and got a lot of “off-by-one” errors. Luckily, testing a lossless encoder is relatively easy, we just need to test the result is the same as the input.

Related Posts

In The Cardinality of Complex Numbers we talk about bijective mapping between $\mathbb{C}$ and $\mathbb{R}$, which one can see as a lossless encoding. This also vaguely reminds me of the Gödel number system which maps mathematical expressions to natural numbers!

One of the most surprising coincidences is that in the recent post about Folly F14 Map we discussed the idea of “shuffling” elements across chunks much in the same way TRANSEncoder shuffles the symbols over the array. Both choose the size to be a power of 2 and the stride or step to be an odd number.

Another connection is that Folly F14 Map is a hash table and the first hashing algorithm I learned and possibly the only one I ever implemented is related to Symmetric Numeral System: we treat the key string as a number in base, say 128 (for ASCII), and then convert it to an decimal modulo a prime number.

Like Huffman, ANS is a lossless compression, but in many cases a lossy compression is a tradeoff we can take, such as in Linear Predictive Coding and T-Digest.

Finally, a weaker connection is with Numerical Representations as inspiration for Data Structures, because it’s also based on number systems.

Source code | RSS feed | Tags: python