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 $\text{delim}[0,1)$:

$$\begin{array}{rl}B& =\text{delim}[b_0,b_1)\\ S& =\text{delim}[s_0,s_1)\end{array}$$

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

$$S\cdot B=\text{delim}[b_o+s_0(b_1-b_0),b_o+s_1(b_1-b_0))$$

64 bit intervals

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

Take $\text{p}\mathtt{b},\mathtt{r}\in \text{p}{{\mathbb{B}}^{64}}^2$ with the following map:

$$\text{p}\mathtt{b},\mathtt{r}\mapsto \text{delim}[\frac{\mathtt{b}}{2^{64}},\frac{\mathtt{b}+\mathtt{r}+1}{2^{64}})$$

Some examples of this map:

\begin{aligned} \p{0,0} &↦ \delim[{0, 2⁻^{64}}) \\ \p{0,2^{64}-1} &↦ \delim[{0, 1}) \\ \p{2^{64}-1, 0} &↦ \delim[{1 - 2⁻^{64}, 1}) \\ \p{2^{64}-1, 2^{64}-1} &↦ \delim[{1 - 2⁻^{64}, 2 - 2⁻^{64}}) \\ \end{aligned}

Where the last example is invalid. For validity we require

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

A reverse map can be created as:

$$\text{delim}[l,h)\mapsto \text{p}\text{ceil}l\cdot 2^{64},\text{floor}h\cdot 2^{64}-\text{ceil}l\cdot 2^{64}-1$$

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

l ≤ 1 - 2⁻^{64}

In order for $\mathtt{r}$ to be in ${\mathbb{B}}^{64}$ we require $0\le \text{floor}h\cdot 2^{64}-\text{ceil}l\cdot 2^{64}-1<2^{64}$. The upper bound is already satisfied. To satisfy the lower bound we require:

$$\begin{array}{rl}0& \le \text{floor}h\cdot 2^{64}-\text{ceil}l\cdot 2^{64}-1\\ & <\text{p}h-l\cdot 2^{64}-1-1\\ l+2^{-63}& <h\end{array}$$

Since $h\le 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\le 1$.

Going full circle results in:

$$\text{delim}[l,h)\mapsto \text{p}\mathtt{b},\mathtt{r}\mapsto \text{delim}[\frac{\text{ceil}l\cdot 2^{64}}{2^{64}},\frac{\text{floor}h\cdot 2^{64}}{2^{64}})$$

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{array}{rl}\text{p}{\mathtt{b}}_B,{\mathtt{r}}_B& \mapsto \text{delim}[\frac{{\mathtt{b}}_B}{2^{64}},\frac{{\mathtt{b}}_B+{\mathtt{r}}_B+1}{2^{64}})\\ \text{p}{\mathtt{b}}_S,{\mathtt{r}}_S& \mapsto \text{delim}[\frac{{\mathtt{b}}_S}{2^{64}},\frac{{\mathtt{b}}_S+{\mathtt{r}}_S+1}{2^{64}})\end{array}$$

$$\text{delim}[\frac{{\mathtt{b}}_B}{2^{64}}+\frac{{\mathtt{b}}_S}{2^{64}}\text{p}\frac{{\mathtt{b}}_B+{\mathtt{r}}_B+1}{2^{64}}-\frac{{\mathtt{b}}_B}{2^{64}},\frac{{\mathtt{b}}_B}{2^{64}}+\frac{{\mathtt{b}}_S+{\mathtt{r}}_S+1}{2^{64}}\text{p}\frac{{\mathtt{b}}_B+{\mathtt{r}}_B+1}{2^{64}}-\frac{{\mathtt{b}}_B}{2^{64}})$$

$$\text{delim}[\frac{{\mathtt{b}}_B}{2^{64}}+\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{128}},\frac{{\mathtt{b}}_B}{2^{64}}+\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{128}})$$

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

$$\begin{array}{rl}\frac{{\mathtt{b}}_B}{2^{64}}+\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{128}}+2^{-63}& <\frac{{\mathtt{b}}_B}{2^{64}}+\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{128}}\\ {\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1+2^{64}& <\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1\end{array}$$

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

$$\begin{array}{rl}\frac{{\mathtt{b}}_B}{2^{64}}+\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{128}}& \le 1\\ {\mathtt{b}}_B\cdot 2^{64}+\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1& \le 2^{128}\end{array}$$

Mapping this back to $\text{p}{{\mathbb{B}}^{64}}^2$ using the reverse map:

$$\begin{array}{rl}{\mathtt{b}}_B^{\prime }& =\text{ceil}l\cdot 2^{64}\\ & =\text{ceil}\text{p}\frac{{\mathtt{b}}_B}{2^{64}}+\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{128}}\cdot 2^{64}\\ & ={\mathtt{b}}_B+\text{ceil}\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{64}}\end{array}$$

Since we have ${\mathtt{b}}_S<2^{64}$ and ${\mathtt{r}}_B+1\le 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{array}{rl}{\mathtt{r}}_B^{\prime }& =\text{floor}h\cdot 2^{64}-\text{ceil}l\cdot 2^{64}-1\\ & =\text{floor}\text{p}\frac{{\mathtt{b}}_B}{2^{64}}+\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{128}}\cdot 2^{64}-{\mathtt{b}}_B+\text{ceil}\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{64}}-1\\ & =\text{floor}\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{64}}-\text{ceil}\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{64}}-1\end{array}$$

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

$$\begin{array}{rl}0& \le \text{floor}\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{64}}-\text{ceil}\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{64}}-1\\ 0& \le \frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{64}}-\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{64}}-1\\ 1& <\frac{\text{p}{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{64}}\\ 2^{64}& \le \text{p}{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1\\ {\mathtt{r}}_S& \ge \frac{2^{64}}{{\mathtt{r}}_B+1}-1\\ {\mathtt{r}}_S& \ge 1\end{array}$$

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

$$\begin{array}{rl}{2}^{64}& >\text{floor}\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{64}}-\text{ceil}\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{64}}-1\\ 2^{64}& >\frac{\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+1\text{p}{\mathtt{r}}_B+1}{2^{64}}-\text{ceil}\frac{{\mathtt{b}}_S\text{p}{\mathtt{r}}_B+1}{2^{64}}-1\end{array}$$

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\le 2^{64}$ and ${\mathtt{r}}_B+1\le 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}$:

$$\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S\cdot {\mathtt{r}}_B+\text{p}{\mathtt{b}}_S+{\mathtt{r}}_S+{\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 $\text{delim}[l,h)$ 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{array}{rl}l& =0.0001101101000111111110110010001001\dots \\ h& =0.0001101101001000000011000010000111\dots \\ h-l& =0.\underset{\mathrm{settled}}{\underbrace{000000000000}}\underset{\mathrm{outstanding}}{\underbrace{0000000}}\underset{\mathrm{active}}{\underbrace{100001111111110\dots }}\end{array}$$

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{array}{rl}h-l& =0.0000000000000000000100000000000000\dots \\ & =0.00000000000000000000\underset{\mathrm{active}}{\underbrace{11111111111111\dots }}\end{array}$$

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=\text{floor}-{\log }_2\text{p}h-l-1$$

The interval is normalized if $n=0$ and thus iff $\frac{1}{2}<h-l$.

$$\begin{array}{rl}{l}^{\prime }=l\cdot 2^n& =1101101000111111.110110010001001\dots \\ h^{\prime }=h\cdot 2^n& =1101101001000000.011000010000111\dots \end{array}$$

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{array}{rl}{l}^{\prime \prime }=l^{\prime }-\text{floor}{l}^{\prime }& =0.110110010001001\dots \\ h^{\prime \prime }=h^{\prime }-\text{floor}{l}^{\prime }& =1.011000010000111\dots \end{array}$$

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

$$0\le l<1$$

$$½<h-l\le 1$$

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

After normalization we end up in one of two cases:

• $\text{delim}[h,l)\subseteq \text{delim}[0,1)$
• $\text{delim}[h,l)⊈\text{delim}[0,1)$

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:

• $\text{delim}[h,l)\subseteq \text{delim}[1,2)$

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

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

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

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

$$\begin{array}{rl}n& =\text{floor}-{\log }_2\text{p}\frac{\mathtt{r}+1}{2^{64}}-1\\ & =63-\text{ceil}{\log }_2\text{p}\mathtt{r}+1\end{array}$$

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·2^n - \floor{l·2^n} = \mod{\frac{\b}{2^{64}⁻^n}}_1 \\ h'' &= h·2^n - \floor{l·2^n} = \frac{\b + \r + 1}{2^{64}⁻^n} - \floor{\frac{\b}{2^{64}⁻^n}}\\ \end{aligned}

Where $\mathrm{mod}{x}_1=x-\text{floor}x$ is the fractional part operator.

\begin{aligned} \b'' &= \ceil{l'' · 2^{64}} \\ &= \ceil{\mod{\frac{\b}{2^{64}⁻^n}}_1 · 2^{64}} \\ &= \ceil{\mod{\frac{\b}{2^{64}⁻^n}_1 · 2^{64}}_{2^{64}}} \\ &= \mod{\b · 2^n}_{2^{64}} \end{aligned}

Where $\mathrm{mod}{x}_{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} \r'' &= \floor{h'' · 2^{64}} - \ceil{l'' · 2^{64}} - 1 \\ &= \floor{h'' · 2^{64}} - \b'' - 1 \\ &= \floor{\p{\frac{\b + \r + 1}{2^{64}⁻^n} - \floor{\frac{\b}{2^{64}⁻^n}}} · 2^{64}} - \b'' - 1 \\ &= \floor{\p{\b + \r + 1}·2^n - \floor{\frac{\b}{2^{64}⁻^n}}· 2^{64} } - \b'' - 1 \\ \end{aligned}

$$\begin{array}{rl}& =\text{floor}\text{p}\frac{\mathtt{b}\cdot 2^n}{2^{64}}-\text{floor}\frac{\mathtt{b}\cdot 2^n}{2^{64}}\cdot 2^{64}+\text{p}\mathtt{r}+1\cdot 2^n-{\mathtt{b}}^{\prime \prime }-1\\ & =\text{floor}\mathrm{mod}\mathtt{b}\cdot 2_{2^{64}}^n+\text{p}\mathtt{r}+1\cdot 2^n-{\mathtt{b}}^{\prime \prime }-1\\ & =\text{floor}{\mathtt{b}}^{\prime \prime }+\text{p}\mathtt{r}+1\cdot 2^n-{\mathtt{b}}^{\prime \prime }-1\\ & =\mathtt{r}\cdot 2^n+2^n-1\end{array}$$

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

pair<uint64_t, uint64_t> mul128_emu(uint64_t a, uint64_t b)
{
	using 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 make_pair(h, l);
}

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