03 - Source Coding: Huffman and Shannon's First Theorem

You want to write the symbols of a source in binary with as few bits as possible, and read the message back without ambiguity. Three questions, in order: which codeword lengths are possible (the Kraft inequality), how far the average length can go down (never below the entropy H; a per-symbol code can stay below H + 1) and which code reaches the minimum (Huffman). In the middle there is a distinction to hold firm: the bound for per-symbol codes and Shannon's source coding theorem are two different statements, the second about long blocks with N tending to infinity. Every listing was run twice with identical output (Python 3.14.4, numpy 2.3.5, fixed seeds, no network).

They presuppose article 01 (entropy) and article 02 (Kullback-Leibler divergence): the lower bound follows from the fact that a divergence is not negative.

Prefix codes and instantaneous decoding

A code assigns a binary word to each symbol; its average length is L = Σ p(x) · |c(x)|.

A code is prefix-free if no word is the beginning of another, and this allows decoding by reading from left to right without looking ahead.

With A = 0, B = 10, C = 11, the sequence 10110 reads 10 | 11 | 0, that is, B, C, A: as soon as the string read coincides with a word, that word is finished, because no other can extend it.

A code is uniquely decodable (UD) if two different sequences of symbols never give the same string. Prefix codes are UD, but the converse is not true, and it is the reason why it matters who states the Kraft inequality. MacKay (§5.2) reports that Kraft proved that, given lengths that satisfy the inequality, a prefix code with those lengths exists, and that the converse, that unique decodability implies the inequality, is due to McMillan (1956): this is why it is also called the Kraft-McMillan inequality.

The listing enumerates all codes with 1 to 4 words of 1 to 3 bits (1470 sets) and for each computes the Kraft sum Σ 2^−l with exact fractions, whether it is prefix-free, and whether it is UD.

For unique decodability it uses the Sardinas-Patterson test, which here is an implementation of mine, not verified against a primary source; as a second independent check it searches by brute force for two sequences of at most 6 words that produce the same string.

from fractions import Fraction
from itertools import combinations, product


def kraft(parole):
    return sum(Fraction(1, 2 ** len(w)) for w in parole)


def senza_prefissi(parole):
    return not any(a != b and b.startswith(a) for a in parole for b in parole)


def univocamente_decodificabile(parole):
    """Sardinas-Patterson test: the code is UD if no 'remainder' becomes a codeword."""
    C = set(parole)
    resti = {b[len(a):] for a in C for b in C if a != b and b.startswith(a)}
    visti = set()
    while resti:
        if resti & C:
            return False
        visti |= resti
        nuovi = set()
        for r in resti:
            for c in C:
                if c.startswith(r) and c != r:
                    nuovi.add(c[len(r):])
                if r.startswith(c) and r != c:
                    nuovi.add(r[len(c):])
        resti = nuovi - visti
    return True


def ambiguo_per_forza_bruta(parole, max_parole=6):
    """Second, independent check: two different sequences of words that give the same string."""
    visti = {}
    for k in range(1, max_parole + 1):
        for seq in product(parole, repeat=k):
            stringa = "".join(seq)
            if stringa in visti and visti[stringa] != seq:
                return True
            visti.setdefault(stringa, seq)
    return False


tutte = ["".join(b) for n in (1, 2, 3) for b in product("01", repeat=n)]   # 14 words
tot = ud = prefix = ud_non_prefix = 0
violazioni_kraft = discordanze = 0
for n in range(1, 5):
    for C in combinations(tutte, n):
        tot += 1
        p, u = senza_prefissi(C), univocamente_decodificabile(C)
        prefix += p
        ud += u
        ud_non_prefix += u and not p
        violazioni_kraft += (p or u) and kraft(C) > 1
        discordanze += (not u) != ambiguo_per_forza_bruta(C)
print(f"codes with 1-4 words of 1-3 bits: {tot}")
print(f"  prefix-free: {prefix}   uniquely decodable: {ud}   UD but with prefixes: {ud_non_prefix}")
print(f"  UD or prefix codes with Kraft sum > 1: {violazioni_kraft}")
print(f"  disagreements between Sardinas-Patterson and brute force (up to 6 words): {discordanze}")
for C in (["0", "10", "11"], ["0", "01", "11"], ["0", "01", "10"]):
    print(C, "Kraft =", kraft(C), "| prefix-free:", senza_prefissi(C), "| UD:", univocamente_decodificabile(C))

# Output:
# codes with 1-4 words of 1-3 bits: 1470
#   prefix-free: 458   uniquely decodable: 686   UD but with prefixes: 228
#   UD or prefix codes with Kraft sum > 1: 0
#   disagreements between Sardinas-Patterson and brute force (up to 6 words): 0
# ['0', '10', '11'] Kraft = 1 | prefix-free: True | UD: True
# ['0', '01', '11'] Kraft = 1 | prefix-free: False | UD: True
# ['0', '01', '10'] Kraft = 1 | prefix-free: False | UD: False

Out of 1470 codes, 458 are prefix-free and 686 UD: 228 are UD even though they have words that begin with others. None of the UD or prefix codes has a Kraft sum greater than 1, and the two methods of deciding UD never disagree.

The last three lines are instructive. {0, 10, 11} is a prefix code. {0, 01, 11} is not, but it is UD (read backwards, the words 0, 10, 11 form a prefix code); to decode it, though, you sometimes have to look ahead, so it is not instantaneous. {0, 01, 10} also has Kraft sum 1, and is not UD: the string 010 reads 0 | 10 or 01 | 0.

The inequality is therefore a condition on the lengths, not on the code: computing the sum is not enough to know that a given code works.

The Kraft inequality, and why it is needed

The inequality says that the lengths l_i of a binary prefix code satisfy Σ 2^−l_i ≤ 1.

The idea can be seen on a complete binary tree of depth L = max l_i, with 2^L leaves: a word of length l "occupies" the 2^L−l leaves below it, and for the code to be prefix-free these leaves do not overlap. Hence Σ 2^L−l_i ≤ 2^L, and dividing by 2^L gives the inequality.

Every short word costs a lot of space (a 1-bit word occupies half the tree), and this is the constraint that prevents giving small lengths to all: three 1-bit words have sum 3/2 > 1 and do not exist.

The practical sense is that the inequality turns a problem about codes into a problem about numbers: every choice of lengths with sum at most 1 can be realised by a prefix code (MacKay, §5.2), and looking for the best code means looking for the integer lengths that minimise L under that constraint.

The bound: entropy, and one bit of margin

Let z = Σ 2^−l_i ≤ 1 and let q_i = 2^−l_i/z, which is a distribution.

From the definition, l_i = log₂(1/q_i) − log₂ z, and substituting in the average length: L = H(p) + D(p||q) + log₂(1/z).

The divergence is non-negative (article 02) and so is log₂(1/z) because z ≤ 1, hence L ≥ H, with equality only if p_i = 2^−l_i for every i: the probabilities must be powers of 1/2. This is MacKay's proof (§5.3), which states it for every uniquely decodable code and uses precisely the Kraft-McMillan inequality.

For the upper bound you choose lengths slightly larger than the ideal ones, l_i = ⌈log₂(1/p_i)⌉.

They satisfy Kraft, because 2^−⌈x⌉ ≤ 2^−x and Σ p_i = 1, and they exceed the ideal length by less than one bit.

From this follows MacKay's theorem 5.1, for per-symbol codes: there exists a prefix code with H(X) ≤ L < H(X) + 1. The "+1" is enormous in proportion when the entropy is small, as we will see shortly.

The Huffman algorithm

MacKay (§5.5) describes the algorithm in two sentences: take the two least probable symbols, which will have the longest codewords, of equal length and differing only in the last bit; merge them into a single symbol and repeat. Each step reduces the alphabet by one, so |A| − 1 steps are needed.

The listing implements it with a priority queue (heapq): each element is a weight, a tie-break counter and the dictionary of the codewords of the symbols merged so far, and at each merge 0 is prepended to the branch of the first node extracted and 1 to the second.

The counter ordine makes the result deterministic when two weights are equal, and it is the only arbitrary choice: another tie-break gives different words, but the same average length.

The other two functions compute L with exact fractions and the entropy in bits. The block prints nothing: it is the definition used by the following listings.

import heapq
from fractions import Fraction
from itertools import count
from math import log2


def huffman(pesi):
    """pesi: dict symbol -> weight (int, Fraction or float). Returns dict symbol -> codeword."""
    if len(pesi) == 1:
        return {s: "0" for s in pesi}
    ordine = count()                         # deterministic tie-break between equal weights
    heap = [(w, next(ordine), {s: ""}) for s, w in pesi.items()]
    heapq.heapify(heap)
    while len(heap) > 1:
        w1, _, c1 = heapq.heappop(heap)      # the two least probable nodes...
        w2, _, c2 = heapq.heappop(heap)
        fuso = {s: "0" + p for s, p in c1.items()} | {s: "1" + p for s, p in c2.items()}
        heapq.heappush(heap, (w1 + w2, next(ordine), fuso))   # ...become a single node
    return heap[0][2]


def lunghezza_media(pesi, codice):
    tot = sum(pesi.values())
    return sum(Fraction(pesi[s]) * len(codice[s]) for s in pesi) / tot


def entropia(pesi):
    tot = sum(pesi.values())
    return sum(float(w / tot) * log2(tot / w) for w in pesi.values())

The next listing applies the function to two sources: one with probabilities that are powers of 1/2 (0.5, 0.25, 0.125, 0.125) and one with five symbols (0.35, 0.25, 0.20, 0.12, 0.08), chosen for the example.

For each it prints the words, L as an exact fraction, H and the check of H ≤ L < H + 1.

R = Fraction
sorgenti = {
    "powers of 2": {"A": R(1, 2), "B": R(1, 4), "C": R(1, 8), "D": R(1, 8)},
    "five symbols": {"A": R(35, 100), "B": R(25, 100), "C": R(20, 100), "D": R(12, 100), "E": R(8, 100)},
}
for nome, pesi in sorgenti.items():
    codice = huffman(pesi)
    L, H = lunghezza_media(pesi, codice), entropia(pesi)
    print(nome)
    for s in sorted(codice, key=lambda s: (len(codice[s]), s)):
        print(f"   {s} p={float(pesi[s]):.3f}  {codice[s]}")
    print(f"   L = {L} = {float(L):.4f} bit   H = {H:.4f} bit   L - H = {float(L) - H:.4f}   (H <= L < H+1: {H <= L < H + 1})")

# Output:
# powers of 2
#    A p=0.500  0
#    B p=0.250  10
#    C p=0.125  110
#    D p=0.125  111
#    L = 7/4 = 1.7500 bit   H = 1.7500 bit   L - H = 0.0000   (H <= L < H+1: True)
# five symbols
#    A p=0.350  11
#    B p=0.250  10
#    C p=0.200  00
#    D p=0.120  011
#    E p=0.080  010
#    L = 11/5 = 2.2000 bit   H = 2.1531 bit   L - H = 0.0469   (H <= L < H+1: True)

With probabilities that are powers of 1/2 the code reaches the entropy: L = 7/4 = H = 1.75 bits, and the code with lengths 1, 2, 3, 3 is precisely the one the equality case speaks of. With five symbols L = 11/5 = 2.2 bits against H = 2.1531: the gap is 0.0469 bits.

Another tie-break between the two equal weights (0.20 and 0.20) would give different words, but the same average length.

That the code is optimal should not be taken on trust. The next listing verifies it in two ways.

First it compares, on 3000 random distributions with 2 to 6 symbols (integer weights, hence exact arithmetic), the Huffman L with the minimum obtained by enumerating all lists of lengths (from 1 to n−1, because a complete binary tree with n leaves is no deeper) that satisfy Kraft.

Then it encodes and decodes all strings up to 6 symbols on the five symbols of the example (19,530) and another 2000 random strings of 300 symbols, and measures the observed average length.

import random
from itertools import product

rng = random.Random(2026)


# For each n, all the lists of lengths (from 1 to n-1) that satisfy Kraft, in integer arithmetic
ammesse = {
    n: [l for l in product(range(1, n), repeat=n)
        if sum(2 ** (n - 1 - x) for x in l) <= 2 ** (n - 1)]
    for n in range(2, 7)
}

# 1) Optimality: comparison with the exhaustive search among ALL the lengths that satisfy Kraft
casi, diversi = 3000, 0
for _ in range(casi):
    n = rng.randint(2, 6)
    pesi = {i: rng.randint(1, 50) for i in range(n)}
    tot = sum(pesi.values())
    L_huffman = lunghezza_media(pesi, huffman(pesi))
    minimo = min(sum(pesi[i] * l[i] for i in range(n)) for l in ammesse[n])
    diversi += L_huffman != Fraction(minimo, tot)
print(f"distributions with L(Huffman) different from the exhaustive minimum: {diversi} out of {casi}")
print("length lists allowed by Kraft, for n=2..6 symbols:", [len(ammesse[n]) for n in range(2, 7)])

# 2) Round trip: all the strings up to 6 symbols, plus 2000 random strings of 300
pesi = {"A": 35, "B": 25, "C": 20, "D": 12, "E": 8}
codice = huffman(pesi)
inverso = {p: s for s, p in codice.items()}


def decodifica(bits):
    out, corrente = [], ""
    for b in bits:
        corrente += b
        if corrente in inverso:              # prefix code: the word is finished as soon as you recognise it
            out.append(inverso[corrente])
            corrente = ""
    assert corrente == ""
    return "".join(out)


controllate = 0
for n in range(1, 7):
    for s in product("ABCDE", repeat=n):
        testo = "".join(s)
        assert decodifica("".join(codice[c] for c in testo)) == testo
        controllate += 1
medie = []
for _ in range(2000):
    testo = "".join(rng.choices("ABCDE", weights=[35, 25, 20, 12, 8], k=300))
    codificato = "".join(codice[c] for c in testo)
    assert decodifica(codificato) == testo
    medie.append(len(codificato) / len(testo))
    controllate += 1
media = sum(medie) / len(medie)
errore = (sum((m - media) ** 2 for m in medie) / (len(medie) - 1)) ** 0.5 / len(medie) ** 0.5
print(f"strings encoded and decoded without errors: {controllate}")
print(f"observed mean length over 600000 symbols: {media:.4f} +/- {errore:.4f} bits (theoretical 2.2000)")

# Output:
# distributions with L(Huffman) different from the exhaustive minimum: 0 out of 3000
# length lists allowed by Kraft, for n=2..6 symbols: [1, 4, 32, 392, 6154]
# strings encoded and decoded without errors: 21530
# observed mean length over 600000 symbols: 2.1987 +/- 0.0005 bits (theoretical 2.2000)

Out of 3000 distributions, 0 differ from the exhaustive minimum; the lists of lengths allowed by Kraft grow fast (1, 4, 32, 392, 6154), which explains why an algorithm is sought and not an enumeration.

The 21,530 strings come back identical after encoding and decoding, and over 600,000 symbols the observed average length is 2.1987 ± 0.0005 bits against the theoretical 2.2000: it is 2.6 standard errors away, more than I expected. With another eight seeds, in a separate trial not reported in the listing, the mean over 600,000 symbols ranges from 2.1993 to 2.2007, so the gap is a fluctuation of this sample and not a defect of the code, which the round trip on 21,530 strings rules out.

The optimality concerns binary per-symbol codes; MacKay (solution of exercise 5.32) describes the modification needed to encode with an alphabet of q symbols.

Where Huffman loses, and the theorem that surpasses it

A per-symbol code cannot use fewer than one bit per symbol, even when the symbol is almost certain. For a source with probabilities (0.99, 0.01) the entropy is 0.0808 bits, but Huffman assigns 1 bit to each of the two symbols and the average length is 1: twelve times the entropy.

MacKay (§5.6) expresses it by saying that Huffman codes, although presented as "optimal", have many practical defects: they are optimal among per-symbol codes, but in practice we do not want a per-symbol code.

The remedy is to encode blocks of N symbols as if they were a single symbol. The listing tests it for N from 1 to 12: the probability of a block depends only on the number of rare symbols, and using the integer weight 99^N−k (the true probability is the weight divided by 100^N) the average length per symbol is exact.

Next to it, it prints the bound H + 1/N, which follows from theorem 5.1 applied to the block source, whose entropy is N times that of the single symbol.

from itertools import product
from math import log2

H = 0.99 * log2(1 / 0.99) + 0.01 * log2(1 / 0.01)
print(f"entropy of one symbol with p = (0.99, 0.01): H = {H:.4f} bit")
print("  N    words   L/N (bits per symbol)   H + 1/N")
for N in (1, 2, 4, 6, 8, 10, 12):
    # blocks of N symbols: the probability depends only on the number k of rare symbols
    # integer weight 99^(N-k): the true probability is weight / 100^N, so the calculation is exact
    pesi = {b: 99 ** (N - sum(b)) for b in product((0, 1), repeat=N)}
    L = lunghezza_media(pesi, huffman(pesi))
    print(f"{N:3d}  {len(pesi):6d}   {float(L) / N:20.4f}   {H + 1 / N:8.4f}")

# Output:
# entropy of one symbol with p = (0.99, 0.01): H = 0.0808 bit
#   N    words   L/N (bits per symbol)   H + 1/N
#   1       2                 1.0000     1.0808
#   2       4                 0.5150     0.5808
#   4      16                 0.2727     0.3308
#   6      64                 0.1953     0.2475
#   8     256                 0.1572     0.2058
#  10    1024                 0.1362     0.1808
#  12    4096                 0.1222     0.1641

The length per symbol falls from 1 to 0.1222 with blocks of 12 symbols, always below the bound H + 1/N, and tends to the entropy, but slowly: at N = 12 it is still 51% above 0.0808, and the code has 4096 words.

This is the operational form of the source coding theorem. MacKay (theorem 4.1, p. 78) states it as follows: N independent and identically distributed random variables, each with entropy H(X), can be compressed into more than N H(X) bits with negligible risk of information loss as N → ∞; conversely, if they are compressed into fewer than N H(X) bits it is virtually certain that information will be lost.

What it guarantees, and what it does not. It guarantees that H is the boundary: you can reach above it, not below.

It guarantees nothing for a finite N, it does not construct the code, and it presupposes that the distribution is known and that the source is i.i.d.; Shannon, in his 1948 paper (theorem 3, on typical sequences), works with ergodic sources, which are more general. Nor does it say that the cost in memory and time of the code is acceptable: the code with blocks of 12 has 4096 words, with blocks of 30 it would have 2³⁰ possible words, over a billion.

Arithmetic coding: the whole sequence as a number

Arithmetic coding abolishes the obligation of an integer number of bits per symbol. MacKay (ch. 6) attributes its invention to Elias, Rissanen and Pasco, and its practical realisation to Witten and others (1987). The idea is to assign to each sequence a sub-interval of [0, 1) whose width is the probability of the sequence: each symbol narrows the interval in proportion to its own probability.

Then a binary fraction is chosen inside the interval: if the width is P, ⌈log₂(1/P)⌉ bits are enough.

The listing implements the encoder and the decoder with exact fractions, hence without rounding errors. intervallo builds the table of cumulative probabilities and narrows low and width symbol by symbol; codifica looks for the smallest k such that 2^−k ≤ width and takes the smallest multiple of 2^−k that is not less than low, which by construction lies in the interval. The decoder must know the number N of symbols (or an end symbol is needed).

I test it on all sequences up to 10 symbols and on 300 random sequences of 1000, with P(A) = 0.99 and P(B) = 0.01.

import random
from math import ceil, log2


def intervallo(seq, prob):
    """Interval [low, low+width) of [0,1) associated with the sequence, with exact fractions."""
    cum, c = {}, Fraction(0)
    for s in sorted(prob):
        cum[s], c = c, c + prob[s]
    low, width = Fraction(0), Fraction(1)
    for s in seq:
        low, width = low + width * cum[s], width * prob[s]
    return low, width, cum


def codifica(seq, prob):
    low, width, _ = intervallo(seq, prob)
    k = 0
    while Fraction(1, 2 ** k) > width:        # smallest k with 2^-k <= width
        k += 1
    m = -((-low.numerator * 2 ** k) // low.denominator)       # ceil(low * 2^k)
    return format(m, f"0{k}b") if k else ""


def decodifica(bits, n, prob):
    x = Fraction(int(bits, 2), 2 ** len(bits)) if bits else Fraction(0)
    _, _, cum = intervallo("", prob)
    low, width, out = Fraction(0), Fraction(1), []
    for _ in range(n):
        for s in sorted(prob):
            a = low + width * cum[s]
            if a <= x < a + width * prob[s]:
                out.append(s)
                low, width = a, width * prob[s]
                break
    return "".join(out)


prob = {"A": Fraction(99, 100), "B": Fraction(1, 100)}
rng = random.Random(5)
# All the sequences up to 10 symbols and 300 random sequences of 1000: exact round trip
for n in range(1, 11):
    for i in range(2 ** n):
        s = "".join("AB"[(i >> j) & 1] for j in range(n))
        assert decodifica(codifica(s, prob), n, prob) == s
lunghezze, ideali = [], []
for _ in range(300):
    s = "".join(rng.choices("AB", weights=[99, 1], k=1000))
    bits = codifica(s, prob)
    assert decodifica(bits, 1000, prob) == s
    _, width, _ = intervallo(s, prob)
    lunghezze.append(len(bits))
    ideali.append(ceil(-log2(float(width))) if width > 0 else None)
H = 0.99 * log2(1 / 0.99) + 0.01 * log2(1 / 0.01)
media = sum(lunghezze) / len(lunghezze)
print("exact round trips on 2046 short sequences and 300 of 1000 symbols")
errore = (sum((x - media) ** 2 for x in lunghezze) / (len(lunghezze) - 1)) ** 0.5 / len(lunghezze) ** 0.5
print(f"mean length of the arithmetic code: {media:.1f} +/- {errore:.1f} bits per 1000 symbols (N·H = {1000 * H:.1f})")
print(f"largest gap from ceil(-log2 P(sequence)): {max(abs(a - b) for a, b in zip(lunghezze, ideali))} bit")
print("Huffman symbol by symbol: 1000 bits per 1000 symbols")

# Output:
# exact round trips on 2046 short sequences and 300 of 1000 symbols
# mean length of the arithmetic code: 81.5 +/- 1.3 bits per 1000 symbols (N·H = 80.8)
# largest gap from ceil(-log2 P(sequence)): 0 bit
# Huffman symbol by symbol: 1000 bits per 1000 symbols

All sequences come back identical, and the code length always coincides with ⌈−log₂ P(sequence)⌉ (largest gap 0 bits).

On average it is 81.5 ± 1.3 bits per 1000 symbols, against N H = 80.8 and the 1000 bits of Huffman per symbol: the gap of 0.7 bits is within the error and of the same order as the rounding up (at most 1 bit). The listing uses unlimited-precision rationals, which grow with the message: real encoders use finite-precision integer arithmetic, and how much they lose I have not measured. As with Huffman, the decoder must know the same probabilistic model.

Real compressors, on data built here

The comparison with zlib, bz2 and lzma from the standard library uses two sources built by me, of 200,000 symbols, one byte per symbol. The first has independent symbols with probabilities (1/2, 1/4, 1/8, 1/8): per-symbol Huffman reaches the entropy. The second is a two-state Markov chain that stays in the current state with probability 0.95: the entropy of a symbol taken alone is 1 bit, but the entropy conditioned on the previous symbol, H(X_n|X_n−1) = 0.95 log₂(1/0.95) + 0.05 log₂(1/0.05) = 0.2864 bits, is much smaller (it is the quantity of article 01 applied to a source with memory).

The listing also prints a line for zlib limited to Huffman coding alone (Z_HUFFMAN_ONLY, without matches).

import bz2
import lzma
import zlib

import numpy as np

rng = np.random.default_rng(3)
N = 200_000


def solo_huffman(d):
    c = zlib.compressobj(9, zlib.DEFLATED, 15, 9, zlib.Z_HUFFMAN_ONLY)   # no LZ77 matches
    return c.compress(d) + c.flush()


def bit_per_simbolo(dati):
    return {nome: 8 * len(f(dati)) / N for nome, f in (
        ("zlib -9", lambda d: zlib.compress(d, 9)),
        ("zlib Huffman only", solo_huffman),
        ("bz2 -9", lambda d: bz2.compress(d, 9)),
        ("lzma -9", lambda d: lzma.compress(d, preset=9)),
    )}


# Source 1: independent symbols, p = (1/2, 1/4, 1/8, 1/8), one byte per symbol
x = rng.choice(4, size=N, p=[0.5, 0.25, 0.125, 0.125])
dati = (x + 65).astype(np.uint8).tobytes()                                    # 'A'..'D'
freq = np.bincount(x) / N
H1 = float(-(freq * np.log2(freq)).sum())
L1 = float(lunghezza_media(dict(enumerate(np.bincount(x).tolist())), huffman(dict(enumerate(np.bincount(x).tolist())))))
print(f"source 1 (i.i.d.): empirical entropy {H1:.4f} bits/symbol; Huffman: {L1:.4f}")
for nome, v in bit_per_simbolo(dati).items():
    print(f"   {nome}: {v:.4f} bits/symbol")
# DEFLATE adds an end-of-block code to the symbols (about one every 32768 symbols): I add it to the counts
con_eob = dict(enumerate(np.bincount(x).tolist())) | {"EOB": 6}
parole = huffman(con_eob)
print(f"   Huffman on the same counts + 6 end-of-block: {float(lunghezza_media({k: v for k, v in con_eob.items() if k != 'EOB'}, parole)):.4f} bits/symbol")

# Source 2: two-state Markov chain, it stays in the current state with probability 0.95
resta = rng.random(N) < 0.95
stato = np.zeros(N, dtype=np.int64)
for i in range(1, N):
    stato[i] = stato[i - 1] if resta[i] else 1 - stato[i - 1]
dati = (stato + 48).astype(np.uint8).tobytes()                                # '0' and '1'
marginale = np.bincount(stato) / N
H_marg = float(-(marginale * np.log2(marginale)).sum())
tasso = 0.95 * np.log2(1 / 0.95) + 0.05 * np.log2(1 / 0.05)   # H(X_n | X_{n-1})
print(f"source 2 (Markov): marginal entropy {H_marg:.4f}; "
      f"conditional entropy H(X_n|X_(n-1)) = {tasso:.4f} bits/symbol")
conti = dict(enumerate(np.bincount(stato).tolist()))
print(f"   Huffman symbol by symbol: {float(lunghezza_media(conti, huffman(conti))):.4f} bits/symbol")
con_eob = conti | {"EOB": 6}
print(f"   Huffman on the same counts + 6 end-of-block: {float(lunghezza_media(conti, huffman(con_eob))):.4f} bits/symbol")
for nome, v in bit_per_simbolo(dati).items():
    print(f"   {nome}: {v:.4f} bits/symbol")

# Output:
# source 1 (i.i.d.): empirical entropy 1.7510 bits/symbol; Huffman: 1.7510
#    zlib -9: 2.1342 bits/symbol
#    zlib Huffman only: 1.8789 bits/symbol
#    bz2 -9: 2.0758 bits/symbol
#    lzma -9: 2.0128 bits/symbol
#    Huffman on the same counts + 6 end-of-block: 1.8759 bits/symbol
# source 2 (Markov): marginal entropy 1.0000; conditional entropy H(X_n|X_(n-1)) = 0.2864 bits/symbol
#    Huffman symbol by symbol: 1.0000 bits/symbol
#    Huffman on the same counts + 6 end-of-block: 1.4986 bits/symbol
#    zlib -9: 0.3930 bits/symbol
#    zlib Huffman only: 1.4930 bits/symbol
#    bz2 -9: 0.3496 bits/symbol
#    lzma -9: 0.4083 bits/symbol

On the memoryless source no general-purpose compressor goes below the entropy, and none reaches it: lzma 2.01, bz2 2.08, zlib 2.13 bits per symbol against 1.751.

Interesting is zlib: with only the LZ77 matches turned off it falls to 1.879, better than with everything on, because on an alphabet of four letters repetitions are random and exploiting them costs more than it saves. Why 1.879 and not 1.751? I hypothesise the DEFLATE end-of-block code (about one every 32,768 symbols, my estimate, which I have not verified in the zlib source): it is one more rare symbol, and Huffman lengthens one of the words by a bit.

With 6 end-of-block symbols the calculation gives 1.876, and for the Markov source 1.4986 against a measured 1.4930: with only two symbols, one more costs half a bit.

On the source with memory the opposite happens: per-symbol Huffman cannot do better than 1 bit, while the compressors reach 0.35-0.41 bits per symbol (bz2 0.3496, zlib 0.3930, lzma 0.4083), a little more than a third, despite not knowing the model.

They do not reach the 0.2864 bits of the conditional entropy. It is the experimental confirmation that the entropy of a symbol taken alone is the right measure only for memoryless sources, and that a compressor that exploits repetitions actually exploits the dependence between symbols.

As for DEFLATE, the specification (RFC 1951) defines it as the combination of the LZ77 algorithm and Huffman coding: the zlib and gzip format therefore contains Huffman as its last stage, but it is not "just Huffman".

Typical mistakes

• Believing that H is reachable by every per-symbol code. It is so only if the probabilities are powers of 1/2; otherwise the average length lies between H and H + 1.
• Using Kraft to conclude that a code works. {0, 01, 10} has sum 1 and is not uniquely decodable.
• Confusing the per-symbol bound with Shannon's theorem. The former holds for every source and gives H + 1; the latter concerns blocks, with N tending to infinity.
• Taking for granted that Huffman is the best possible. It is the best among per-symbol codes: with (0.99, 0.01) it gives 1 bit against 0.0808.
• Expecting a general-purpose compressor to approach the entropy of the symbol. On the memoryless source it stays 0.26-0.38 bits above; on the one with memory it goes far below.

The five mistakes stem from one and the same confusion: talking about "compression" without saying with respect to which model, to which class of codes and with which block length. Every comparison with the entropy makes sense only if these three elements are declared.

What you can do now

You can verify whether a set of lengths can be realised by a prefix code and you know that the Kraft sum says nothing about a given code. You can prove that no per-symbol code lies below H and that one exists below H + 1, build the Huffman code, test it by round trip and compare it with the exhaustive search.

You also know where it stops: for very concentrated sources you need blocks or arithmetic coding, and Shannon's theorem says where the boundary is, not how to get there. The next step, in article 04, moves the question from compressing to transmitting over a channel that makes errors, and in article 05 error-correcting codes are built. Lossy compression is in article 06, and the index is on the series page.