# Entropy Coding

In these lab notes I derive an arithmetic coder that is optimally efficient using 64 bit registers.

## Range multiplication

Given two sub-intervals of $\left[{0, 1}\right)$:

\begin{aligned} B &= \left[{b_0 , b_1 }\right) \\ S &= \left[{s_0 , s_1 }\right) \end{aligned}

and we want to compute the product of $S$ into $B$:

$S \cdot B = \left[{bₒ + s_0 (b_1 - b_0 ), bₒ + s_1 (b_1 - b_0 )}\right)$

## 64 bit intervals

We want to represent non-empty subintervals of $\left[{0,1}\right)$ with 64 bit integers, $\mathbb{B} ^{64}$.

Take $\left({\mathtt{b},\mathtt{r}}\right) \in \left({\mathbb{B} ^{64}}\right)^2$ with the following map:

$\left({\mathtt{b},\mathtt{r}}\right) \mapsto \left[{\frac{\mathtt{b}}{2^{64}}, \frac{\mathtt{b}+ \mathtt{r}+ 1}{2^{64}}}\right)$

Some examples of this map:

\begin{aligned} \left({0,0}\right) &\mapsto \left[{0, 2^- ^{64}}\right) \\ \left({0,2^{64}-1}\right) &\mapsto \left[{0, 1}\right) \\ \left({2^{64}-1, 0}\right) &\mapsto \left[{1 - 2^- ^{64}, 1}\right) \\ \left({2^{64}-1, 2^{64}-1}\right) &\mapsto \left[{1 - 2^- ^{64}, 2 - 2^- ^{64}}\right) \\ \end{aligned}

Where the last example is invalid. For validity we require

$\mathtt{b}+ \mathtt{r}+ 1 \leq 2^{64}$

A reverse map can be created as:

$\left[{l, h}\right) \mapsto \left({\left\lceil {l \cdot 2^{64}}\right\rceil , \left\lfloor {h \cdot 2^{64}}\right\rfloor - \left\lceil {l \cdot 2^{64}}\right\rceil - 1}\right)$

In order for $\mathtt{b}$ to be in $\mathbb{B} ^{64}$ we require $0 \leq \left\lceil {l \cdot 2^{64}}\right\rceil < 2^{64}$. We already have $0 \leq l < 1$. To satisfy the upper bound we additionally require:

$l \leq 1 - 2^- ^{64}$

In order for $\mathtt{r}$ to be in $\mathbb{B} ^{64}$ we require $0 \leq \left\lfloor {h \cdot 2^{64}}\right\rfloor - \left\lceil {l \cdot 2^{64}}\right\rceil - 1 < 2^{64}$. The upper bound is already satisfied. To satisfy the lower bound we require:

\begin{aligned} 0 &\leq \left\lfloor {h \cdot 2^{64}}\right\rfloor - \left\lceil {l \cdot 2^{64}}\right\rceil - 1 \\ &< \left({h - l}\right) \cdot 2^{64} - 1 - 1 \\ l + 2^- ^{63} &< h \end{aligned}

Since $h \leq 1$ this also implies the stronger bound $l < 1 - 2^- ^{63}$ on $l$. For the reverse map to work we need $l + 2^- ^{63} < h \leq 1$.

Going full circle results in:

$\left[{l, h}\right) \mapsto \left({\mathtt{b},\mathtt{r}}\right) \mapsto \left[{\frac{\left\lceil {l \cdot 2^{64}}\right\rceil }{2^{64}}, \frac{\left\lfloor {h \cdot 2^{64}}\right\rfloor }{2^{64}}}\right)$

It can be easily shown that the later interval is always a subinterval of the former. It is also easily seen that this operation is idempotent.

## 64 bit interval multiply

\begin{aligned} \left({\mathtt{b}_B,\mathtt{r}_B}\right) &\mapsto \left[{\frac{\mathtt{b}_B}{2^{64}}, \frac{\mathtt{b}_B + \mathtt{r}_B + 1}{2^{64}}}\right) \\ \left({\mathtt{b}_S,\mathtt{r}_S}\right) &\mapsto \left[{\frac{\mathtt{b}_S}{2^{64}}, \frac{\mathtt{b}_S + \mathtt{r}_S + 1}{2^{64}}}\right) \end{aligned}

$\left[{ \frac{\mathtt{b}_B}{2^{64}} + \frac{\mathtt{b}_S}{2^{64}} \left({\frac{\mathtt{b}_B + \mathtt{r}_B + 1}{2^{64}} - \frac{\mathtt{b}_B}{2^{64}}}\right), \frac{\mathtt{b}_B}{2^{64}} + \frac{\mathtt{b}_S + \mathtt{r}_S + 1}{2^{64}} \left({\frac{\mathtt{b}_B + \mathtt{r}_B + 1}{2^{64}} - \frac{\mathtt{b}_B}{2^{64}}}\right) }\right)$

$\left[{ \frac{\mathtt{b}_B}{2^{64}} + \frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{128}}, \frac{\mathtt{b}_B}{2^{64}} + \frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{128}} }\right)$

The first bound in $l + 2^- ^{63} < h \leq 1$ gives:

\begin{aligned} \frac{\mathtt{b}_B}{2^{64}} + \frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{128}} + 2^- ^{63} &< \frac{\mathtt{b}_B}{2^{64}} + \frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{128}} \\ \mathtt{b}_S\left({\mathtt{r}_B + 1}\right) + 2^{64} &< \left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right) \end{aligned}

The second bound in $l + 2^- ^{63} < h \leq 1$ gives:

\begin{aligned} \frac{\mathtt{b}_B}{2^{64}} + \frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{128}} & \leq 1 \\ \mathtt{b}_B \cdot 2^{64} + \left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right) & \leq 2^{128} \end{aligned}

Mapping this back to $\left({\mathbb{B} ^{64}}\right)^2$ using the reverse map:

\begin{aligned} \mathtt{b}_B' &= \left\lceil {l \cdot 2^{64}}\right\rceil \\ &= \left\lceil {\left({ \frac{\mathtt{b}_B}{2^{64}} + \frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{128}} }\right)\cdot 2^{64}}\right\rceil \\ &= \mathtt{b}_B + \left\lceil {\frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{64}} }\right\rceil \\ \end{aligned}

Since we have $\mathtt{b}_S < 2^{64}$ and $\mathtt{r}_B + 1 \leq 2^{64}$ we can implement this using a 64 bit multiply with 128 bit result, and an 128 bit and 64 bit add as $\mathtt{b}_S \cdot \mathtt{r}_B + \mathtt{b}_S$, or in code:

uint64 h, l;
std::tie(h, l) = mul128(b_S, r_B);
std::tie(h, l) = add128(h, l, b_S);
const uint64 t = h + (l > 0 ? 1 : 0);
b_B += t;

Where the final term implements the ceiling operator.

\begin{aligned} \mathtt{r}_B' &= \left\lfloor {h \cdot 2^{64}}\right\rfloor - \left\lceil {l \cdot 2^{64}}\right\rceil - 1 \\ &= \left\lfloor {\left({\frac{\mathtt{b}_B}{2^{64}} + \frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{128}}}\right)\cdot 2^{64}}\right\rfloor - \mathtt{b}_B + \left\lceil {\frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rceil - 1 \\ &= \left\lfloor {\frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rfloor - \left\lceil {\frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rceil - 1 \\ \end{aligned}

To obtain a valid result we need $0 \leq \mathtt{r}_B' < 2^{64}$. The lower bound goes like:

\begin{aligned} 0 &\leq \left\lfloor {\frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rfloor - \left\lceil {\frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rceil - 1 \\ 0 &\leq \frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{64}} - \frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{64}} - 1 \\ 1 &< \frac{\left({\mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{64}} \\ 2^{64} &\leq \left({\mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right) \\ \mathtt{r}_S &\geq \frac{2^{64}}{\mathtt{r}_B + 1} - 1 \\ \mathtt{r}_S &\geq 1 \end{aligned}

Where in the last step I have used $\mathtt{r}_B \geq 2^{63}$, which will be enforced later, in normalization. The upper bound $\mathtt{r}_B' < 2^{64}$ goes like:

\begin{aligned} 2^{64} &> \left\lfloor {\frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rfloor - \left\lceil {\frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rceil - 1 \\ 2^{64} &> \frac{\left({\mathtt{b}_S + \mathtt{r}_S + 1}\right)\left({\mathtt{r}_B + 1}\right)}{2^{64}} - \left\lceil {\frac{\mathtt{b}_S\left({\mathtt{r}_B + 1}\right)}{2^{64}}}\right\rceil - 1 \\ \end{aligned}

The second term we already have as t. For the first term we need another tricky multiply. In this case we have $\mathtt{b}_S + \mathtt{r}_S + 1 \leq 2^{64}$ and $\mathtt{r}_B + 1 \leq 2^{64}$, so the intermediate result can actually be $2^{64}$. The $-1$ will make sure the final result will be in $\mathbb{B} ^{64}$. We must be careful, but modular arithmetic will handle the overflow fine. Let’s rewrite the multiplication in values $<2^{64}$:

$\left({\mathtt{b}_S + \mathtt{r}_S}\right)\cdot \mathtt{r}_B + \left({\mathtt{b}_S + \mathtt{r}_S}\right) + \mathtt{r}_B + 1$

Now we can calculate the final interval:

std::tie(h, l) = mul128(b_S + r_S, r_B);
std::tie(h, l) = add128(h, l, b_S + r_S);
std::tie(h, l) = add128(h, l, r_B);
std::tie(h, l) = add128(h, l, 1);
r_B = h - t - 1;

## Normalization

Given an interval $\left[{l,h}\right)$ we want to extract the prefix bits that won’t change anymore. To determine which bits won’t change it is useful to look at $h-l$:

\begin{aligned} l &= \mathtt{0.0001101101000111111110110010001001}\ldots \\ h &= \mathtt{0.0001101101001000000011000010000111}\ldots \\ h - l &= \mathtt{0.} \underbrace{\mathtt{000000000000}}_{\mathrm{settled}} \underbrace{\mathtt{0000000}}_{\mathrm{outstanding}} \underbrace{\mathtt{100001111111110}\ldots }_{\mathrm{active}} \end{aligned}

The settled bits can be written directly to the output. The outstanding bits can still change because of a carry, but are otherwise settled. To normalize this interval we output the first 12 bits, note that there are 7 bits outstanding, and rescale the interval by $2^{12}^+ ^7$. There is one edge case we can consider:

\begin{aligned} h - l &= \mathtt{0.0000000000000000000100000000000000}\ldots \\ &= \mathtt{0.00000000000000000000} \underbrace{\mathtt{11111111111111}\ldots }_{\mathrm{active}} \end{aligned}

Writing the number in two different but equal ways can result in a different number of leading zeros. We take the one that results in the largest number of leading zeros. So in general we want to rescale by $2^n$ with $n$ given by:

$n = \left\lfloor {-\log _2 \left({h - l}\right)}\right\rfloor - 1$

The interval is normalized if $n=0$ and thus iff $½ < h - l$.

\begin{aligned} l' = l\cdot 2^n &= \mathtt{1101101000111111.110110010001001}\ldots \\ h' = h\cdot 2^n &= \mathtt{1101101001000000.011000010000111}\ldots \\ \end{aligned}

Here we can be in one of two situations, either the integral parts are the same, or $h$’s is one larger that $l$’s. In this case $h$’s is in fact one larger. We can subtract the integral part of $l$:

\begin{aligned} l'' = l'-\left\lfloor {l'}\right\rfloor &= \mathtt{0.110110010001001}\ldots \\ h'' = h'-\left\lfloor {l'}\right\rfloor &= \mathtt{1.011000010000111}\ldots \\ \end{aligned}

We have now scaled our $\left[{l,h}\right)$ interval to a subinterval of $\left[{0,2}\right)$, but we also have

$0 \leq l < 1$

$½ < h-l \leq 1$

In summary, an interval is normalized when it is a subinterval of $\left[{0,2}\right)$ and the above inequalities hold.

After normalization we end up in one of two cases:

• $\left[{h,l}\right) ⊆ \left[{0,1}\right)$
• $\left[{h,l}\right) ⊈ \left[{0,1}\right)$

Our situation is as above, there are no bits outstanding. This is when we have a carry outstanding and $h > 1$. After further narrowing of the interval we will eventually end up in the first case and flush the carry buffer, or we will end up in a third case:

• $\left[{h,l}\right) ⊆ \left[{1,2}\right)$

In this third case we should add the carry and subtract $1$ from both $h$ and $l$. We are then effectively back in the first case.

## 64 bit normalization

$\left[{\frac{\mathtt{b}}{2^{64}}, \frac{\mathtt{b}+ \mathtt{r}+ 1}{2^{64}}}\right)$

For this to be normalized in $\left[{0, 2}\right)$ we must have

And the normalization condition $½ < h-l \leq 1$ reduces to $2^{63} \leq \mathtt{r}< 2^{64}$. This essentially means that $\mathtt{r}\in \mathbb{B} ^{64}$ and $\mathtt{r}$ must have the most significant bit set.

\begin{aligned} n &= \left\lfloor {-\log _2 \left({\frac{\mathtt{r}+ 1}{2^{64}}}\right)}\right\rfloor - 1 \\ &= 63 - \left\lceil {\log _2 \left({\mathtt{r}+ 1}\right)}\right\rceil \\ \end{aligned}

The ceiling and the $+1$ cancel, except when $\mathtt{r}=0$. What remains is essentially a count leading zeroes operation:

const uint n = r == 0 ? 63 : count_leading_zeros(r);

\begin{aligned} l'' &= l\cdot 2^n - \left\lfloor {l\cdot 2^n }\right\rfloor = \left[{\frac{\mathtt{b}}{2^{64}^- ^n }}\right]_1 \\ h'' &= h\cdot 2^n - \left\lfloor {l\cdot 2^n }\right\rfloor = \frac{\mathtt{b}+ \mathtt{r}+ 1}{2^{64}^- ^n } - \left\lfloor {\frac{\mathtt{b}}{2^{64}^- ^n }}\right\rfloor \\ \end{aligned}

Where $\left[{x}\right]_1 = x - \left\lfloor {x}\right\rfloor$ is the fractional part operator.

\begin{aligned} \mathtt{b}'' &= \left\lceil {l'' \cdot 2^{64}}\right\rceil \\ &= \left\lceil {\left[{\frac{\mathtt{b}}{2^{64}^- ^n }}\right]_1 \cdot 2^{64}}\right\rceil \\ &= \left\lceil {\left[{\frac{\mathtt{b}}{2^{64}^- ^n } \cdot 2^{64}}\right]_{2^{64}}}\right\rceil \\ &= \left[{\mathtt{b}\cdot 2^n }\right]_{2^{64}} \end{aligned}

Where $\left[{x}\right]_{2^{64}} = x \;\mathrm{mod}\; 2^{64}$ is the modular operator. Since $\mathtt{b}\cdot 2^n$ is strictly integer, the ceiling has no effect and the operation reduces to a simple right shift:

b <<= n;

\begin{aligned} \mathtt{r}'' &= \left\lfloor {h'' \cdot 2^{64}}\right\rfloor - \left\lceil {l'' \cdot 2^{64}}\right\rceil - 1 \\ &= \left\lfloor {h'' \cdot 2^{64}}\right\rfloor - \mathtt{b}'' - 1 \\ &= \left\lfloor {\left({\frac{\mathtt{b}+ \mathtt{r}+ 1}{2^{64}^- ^n } - \left\lfloor {\frac{\mathtt{b}}{2^{64}^- ^n }}\right\rfloor }\right) \cdot 2^{64}}\right\rfloor - \mathtt{b}'' - 1 \\ &= \left\lfloor {\left({\mathtt{b}+ \mathtt{r}+ 1}\right)\cdot 2^n - \left\lfloor {\frac{\mathtt{b}}{2^{64}^- ^n }}\right\rfloor \cdot 2^{64} }\right\rfloor - \mathtt{b}'' - 1 \\ &= \left\lfloor {\left({\frac{\mathtt{b}\cdot 2^n }{2^{64}} - \left\lfloor {\frac{\mathtt{b}\cdot 2^n }{2^{64}}}\right\rfloor }\right)\cdot 2^{64} + \left({\mathtt{r}+ 1}\right)\cdot 2^n }\right\rfloor - \mathtt{b}'' - 1 \\ &= \left\lfloor {\left[{\mathtt{b}\cdot 2^n }\right]_{2^{64}} + \left({\mathtt{r}+ 1}\right)\cdot 2^n }\right\rfloor - \mathtt{b}'' - 1 \\ &= \left\lfloor {\mathtt{b}'' + \left({\mathtt{r}+ 1}\right)\cdot 2^n }\right\rfloor - \mathtt{b}'' - 1 \\ &= \mathtt{r}\cdot 2^n + 2^n - 1 \\ \end{aligned}

This is essentially shifting $\mathtt{r}$ to the right $n$ places while shifting in ones. We don’t need to worry about overflow here because $\mathtt{r}$ is strictly less than $2^{64}^- ^n$.

r <<= n;
r += (1UL << n) - 1;

## Appendix: 128 Bit arithmetic

std::pair<std::uint64_t, std::uint64_t> mul128_emu(std::uint64_t a, std::uint64_t b)
{
using std::uint64_t;
const uint64_t u1 = (a & 0xffffffff);
const uint64_t v1 = (b & 0xffffffff);
uint64_t t = (u1 * v1);
const uint64_t w3 = (t & 0xffffffff);
uint64_t k = (t >> 32);

a >>= 32;
t = (a * v1) + k;
k = (t & 0xffffffff);
uint64_t w1 = (t >> 32);

b >>= 32;
t = (u1 * b) + k;
k = (t >> 32);

const uint64_t h = (a * b) + w1 + k;
const uint64_t l = (t << 32) + w3;
return std::make_pair(h, l);
}

std::pair<uint64, uint64> add128(uint64 h, uint64 l, uint64 n)
{
l += n;
h += (l < n) ? 1 : 0;
return std::make_pair(h, l);
}