Channel Capacity & Error-Correcting Codes

Before Shannon, engineers believed a simple rule: to get fewer errors, send slower. Repeat each bit more times, and errors drop, but the rate drops to zero with them.

Shannon's A Mathematical Theory of Communication (1948) said this is wrong.

The proof was not constructive. It showed that a random code works, but gave no practical way to decode one. Closing that gap took about 50 years.

Send X with distribution p(x). The receiver sees Y. How much did it learn?

I(X;Y) = H(X) − H(X|Y) = H(Y) − H(Y|X)

• H(X): the doubt before receiving.
• H(X|Y): the doubt left after. Noise shows up here as equivocation.
• The difference is what got through, in bits per channel use.

Note that I(X;Y) depends on both the channel W (fixed by physics) and the input distribution p(x) (our choice).

On BSC(0.1) a uniform input gets 0.531 bits through per use. A skewed input gets less. Capacity is the best we can do by choosing p(x).

def mutual_info(px, W):
    """px: input distribution, W[x][y] = P(y | x). Returns I(X;Y) in bits."""
    ny = len(W[0])
    py = [sum(px[x] * W[x][y] for x in range(len(px))) for y in range(ny)]
    return sum(px[x] * W[x][y] * math.log2(W[x][y] / py[y])
               for x in range(len(px)) for y in range(ny)
               if px[x] > 0 and W[x][y] > 0)

BSC = lambda p: [[1 - p, p], [p, 1 - p]]
BEC = lambda e: [[1 - e, e, 0], [0, e, 1 - e]]      # outputs 0, ?, 1
Z   = lambda p: [[1, 0], [p, 1 - p]]                # 1 flips to 0 w.p. p
print(f"BSC(0.1), uniform input: {mutual_info([0.5, 0.5], BSC(0.1)):.4f}")
print(f"BSC(0.1), px=[0.8,0.2]:  {mutual_info([0.8, 0.2], BSC(0.1)):.4f}")
print(f"BEC(0.3), uniform input: {mutual_info([0.5, 0.5], BEC(0.3)):.4f}")

BSC(0.1), uniform input: 0.5310
BSC(0.1), px=[0.8,0.2]:  0.3578
BEC(0.3), uniform input: 0.7000

• 0 ≤ C ≤ min(log|𝒳|, log|𝒴|).
• I(X;Y) is concave in p(x), so any local max is the global max. The problem is a convex program.
• For symmetric channels the uniform input is optimal.

This is the same alternating-maximization idea as EM.

def capacity(W, iters=2000):
    """Blahut-Arimoto: maximize I(X;Y) over input distributions."""
    nx, ny = len(W), len(W[0])
    px = [1 / nx] * nx
    for _ in range(iters):
        py = [sum(px[x] * W[x][y] for x in range(nx)) for y in range(ny)]
        d = [sum(W[x][y] * math.log2(W[x][y] / py[y]) for y in range(ny) if W[x][y] > 0)
             for x in range(nx)]
        w = [px[x] * 2 ** d[x] for x in range(nx)]
        px = [v / sum(w) for v in w]
    return mutual_info(px, W), px

def z_capacity(q):                     # closed form for the Z-channel
    return math.log2(1 + (1 - q) * q ** (q / (1 - q)))

for name, W, exact in [("BSC(0.11)", BSC(0.11), 1 - h2(0.11)),
                       ("BEC(0.3)", BEC(0.3), 1 - 0.3),
                       ("Z(0.5)", Z(0.5), z_capacity(0.5))]:
    C, px = capacity(W)
    print(f"{name:9}: C={C:.4f}  exact={exact:.4f}  best px=[{px[0]:.3f}, {px[1]:.3f}]")

BSC(0.11): C=0.5001  exact=0.5001  best px=[0.500, 0.500]
BEC(0.3) : C=0.7000  exact=0.7000  best px=[0.500, 0.500]
Z(0.5)   : C=0.3219  exact=0.3219  best px=[0.600, 0.400]

I(X;Y) = H(Y) − H(Y|X)
    = H(Y) − ∑_x p(x) H(Y | X=x)
    = H(Y) − ∑_x p(x) h(p)
    = H(Y) − h(p)
    ≤ 1 − h(p).

Given X, the output is just "flipped or not", so H(Y|X) = h(p) whatever we send. Y is binary, so H(Y) ≤ 1. A uniform input makes Y uniform and reaches it. ∎

for p in [0.0, 0.01, 0.05, 0.11, 0.2, 0.5]:
    print(f"p={p:<5} C = 1 - h(p) = {1 - h2(p):.4f}")

p=0.0   C = 1 - h(p) = 1.0000
p=0.01  C = 1 - h(p) = 0.9192
p=0.05  C = 1 - h(p) = 0.7136
p=0.11  C = 1 - h(p) = 0.5001
p=0.2   C = 1 - h(p) = 0.2781
p=0.5   C = 1 - h(p) = 0.0000

• At p = 1/2 the output is independent of the input: nothing gets through.
• p = 1 is perfect: just flip every bit back.
• At 11% bit flips you can still send half a bit per use, reliably.

Question. A link flips 11% of its bits at random. How many real data bits can each channel use carry?

1. Step 1. Binary entropy of the noise:

   h(p) = −p log₂ p − (1−p) log₂(1−p)

2. Step 2. First term: log₂ 0.11 ≈ −3.184, so −0.11 × −3.184 ≈ 0.350.
3. Step 3. Second term: log₂ 0.89 ≈ −0.168, so −0.89 × −0.168 ≈ 0.150.
4. Step 4. Add: h(0.11) ≈ 0.500 bits of doubt per received bit.
5. Step 5. C = 1 − h(p) ≈ 0.500 bits per use.

Binary erasure channel

Let E be the event "erased". It is independent of X and can be read off Y.

H(X|Y) = P(E)·H(X|E) + P(not E)·0
           = ε H(X)
I(X;Y) = H(X) − εH(X) = (1−ε) H(X) ≤ 1 − ε.

Equality at uniform input. ∎

Intuition: a fraction ε of bits vanish. Even a genie who tells the sender which bits will be lost could not do better than 1 − ε. Shannon says you can match the genie without feedback.

Compare: BEC(0.11) has capacity 0.89. BSC(0.11) has 0.5. Knowing where errors are is worth a lot.

Z-channel

A 1 turns into 0 with probability q. Let α = P(X=1). Then

I = h(α(1−q)) − α h(q).

Set the derivative in α to zero and simplify:

C = log₂(1 + (1−q) q^q/(1−q))

For q = 0.5: C = log₂ 1.25 ≈ 0.3219, reached at α = 0.4, matching Blahut–Arimoto's [0.600, 0.400].

The best input is not uniform: we send the fragile symbol 1 less often. For asymmetric channels, the input distribution is part of the design.

Two surprises in one line:

• Rate and reliability are not at war. Fix any R < C. Longer blocks give lower error, at the same rate.
• The operational limit equals an information quantity. "Best achievable rate" is a question about all codes. max I(X;Y) is a formula. They are equal.

Stronger forms:

• Error exponent: the error falls like 2^−nE(R) with E(R) > 0 for R < C.
• Strong converse (Wolfowitz 1957): above C the error goes to 1, not just away from 0.

Capacity is a hard wall with a gentle slope in front of it. Engineering since 1948 has been about getting close to the wall at practical cost.

1. Random codebook. Pick p(x) that achieves C. Draw M = 2^nR codewords, each symbol i.i.d. from p(x).
2. Joint typicality decoding. On receiving yⁿ, output the unique w with (xⁿ(w), yⁿ) jointly typical. If none or many, declare error.
3. Right codeword fails? By the joint AEP, the true pair is typical with probability → 1.
4. A wrong codeword passes? A wrong xⁿ is independent of yⁿ. It looks jointly typical with probability about 2^−nI(X;Y).
5. Union bound over the 2^nR wrong words: P(err) ≤ δ + 2^nR 2^{−n(C−3ε)} → 0 if R < C − 3ε.
6. Averaging. The average random code is good, so some fixed code is good. Throw away the worst half of its codewords to make the max error small too.

Sphere-packing picture: each codeword's noisy outputs form a cloud. We can fit about 2^nI non-overlapping clouds into the output space, so about nI bits per block.

def random_code_error(n, R, p, trials, rng):
    M = 2 ** round(n * R)
    code = [rng.getrandbits(n) for _ in range(M)]
    errors = 0
    for _ in range(trials):
        i = rng.randrange(M)
        noise = sum(1 << b for b in range(n) if rng.random() < p)
        y = code[i] ^ noise
        best = min(range(M), key=lambda j: bin(code[j] ^ y).count("1"))
        errors += best != i
    return errors / trials

rng = random.Random(42)
p, R = 0.05, 0.25                          # C = 1 - h(0.05) = 0.714
print(f"C = {1 - h2(p):.3f}, R = {R}")
for n in [8, 16, 24, 32, 40]:
    print(f"n={n:>2}  M={2 ** round(n * R):>4}  block error={random_code_error(n, R, p, 2000, rng):.4f}")

C = 0.714, R = 0.25
n= 8  M=   4  block error=0.0125
n=16  M=  16  block error=0.0105
n=24  M=  64  block error=0.0045
n=32  M= 256  block error=0.0000
n=40  M=1024  block error=0.0000

No structure at all: codewords are random bit strings. The decoder picks the nearest codeword (maximum likelihood on a BSC with p < 1/2).

The rate stays fixed at 0.25, yet the block error falls as n grows. 0 errors in 2000 trials at n = 32.

The history of coding theory is the search for codes that are random enough to near capacity but structured enough to decode fast.

Why: let E flag an error. Knowing Ŵ, first say whether E happened (≤ 1 bit). If it did, W is one of M−1 others (≤ log M bits, weighted by P_e).

Putting it together

Take W uniform, so H(W) = nR.

nR = H(W) = H(W|Ŵ) + I(W;Ŵ)
    ≤ 1 + P_e nR + nC.

Divide by nR:

P_e ≥ 1 − C/R − 1/(nR).

If R > C, the error stays above 1 − C/R > 0 forever. ∎

def converse_bound(R, C, n):
    return 1 - C / R - 1 / (n * R)
C = 1 - h2(0.11)
for R in [0.6, 0.75, 0.9]:
    print(f"R={R}: P_e >= {converse_bound(R, C, 1000):.3f}  (n=1000, C={C:.3f})")

R=0.6: P_e >= 0.165  (n=1000, C=0.500)
R=0.75: P_e >= 0.332  (n=1000, C=0.500)
R=0.9: P_e >= 0.443  (n=1000, C=0.500)

def rep_error(n, p):
    """Exact bit error of an n-fold repetition code (n odd) with majority vote."""
    return sum(math.comb(n, k) * p**k * (1 - p)**(n - k) for k in range(n // 2 + 1, n + 1))

def simulate_rep(n, p, trials, rng):
    errors = 0
    for _ in range(trials):
        flips = sum(rng.random() < p for _ in range(n))
        errors += flips > n // 2
    return errors / trials

rng = random.Random(1948)
p = 0.1
for n in [1, 3, 5, 7, 9, 11]:
    print(f"n={n:>2} rate={1/n:.3f}  exact={rep_error(n, p):.2e}  "
          f"simulated={simulate_rep(n, p, 200_000, rng):.2e}")

n= 1 rate=1.000  exact=1.00e-01  simulated=9.92e-02
n= 3 rate=0.333  exact=2.80e-02  simulated=2.77e-02
n= 5 rate=0.200  exact=8.56e-03  simulated=8.64e-03
n= 7 rate=0.143  exact=2.73e-03  simulated=2.68e-03
n= 9 rate=0.111  exact=8.91e-04  simulated=1.00e-03
n=11 rate=0.091  exact=2.96e-04  simulated=2.70e-04

Each step down in error costs a lot of rate. To get error → 0 the rate must go → 0.

Shannon says rate 0.5 with error 10⁻⁹ is possible on this channel. Repetition at error 10^−3.5 already runs at rate 0.09.

Proof. ≤ d−1 flips cannot turn one codeword into another. If balls of radius t around codewords overlapped, the triangle inequality gives two codewords within 2t < d. ∎

def hamming(a, b):
    return sum(x != y for x, y in zip(a, b))

def min_distance(code):
    return min(hamming(a, b) for a, b in itertools.combinations(code, 2))

rep3 = ["000", "111"]
parity = [a + str(a.count("1") % 2) for a in map("".join, itertools.product("01", repeat=3))]
print(hamming("1011101", "1001001"), min_distance(rep3), min_distance(parity))
for name, d in [("rep3", 3), ("parity", 2), ("Hamming(7,4)", 3)]:
    print(f"{name:13} d={d}: detects {d - 1}, corrects {(d - 1) // 2}")

2 3 2
rep3          d=3: detects 2, corrects 1
parity        d=2: detects 1, corrects 0
Hamming(7,4)  d=3: detects 2, corrects 1

• Generator matrix G (k × n): encode message m as c = mG.
• Parity-check matrix H ((n−k) × n): c is a codeword iff Hc^T = 0.
• Systematic form: G = [I_k | P], H = [P^T | I_n−k]. The message appears in the codeword as is.

If e is a single 1 at position i, then He^T is column i of H.

So: if all columns of H are nonzero and distinct, the syndrome of a single error names its position. That is the whole idea of Hamming codes.

Syndrome decoding needs a table of 2^n−k syndromes, not 2^k codewords. For Hamming(7,4): 8 entries.

Use all 7 nonzero 3-bit columns in H. Then every single error has its own syndrome. With 3 parity bits we protect 4 data bits: rate 4/7 ≈ 0.571, d = 3.

G = [ I₄ | P ] =
1 0 0 0 | 1 1 0
0 1 0 0 | 1 0 1
0 0 1 0 | 0 1 1
0 0 0 1 | 1 1 1

H = [ P^T | I₃ ] =
1 1 0 1 | 1 0 0
1 0 1 1 | 0 1 0
0 1 1 1 | 0 0 1

Parity bits: p₁ = d₁+d₂+d₄, p₂ = d₁+d₃+d₄, p₃ = d₂+d₃+d₄.

One flipped bit breaks exactly the circles that contain it. Each bit sits in a different set of circles, so the broken circles point to it. The pattern of broken circles is the syndrome.

Message 1011 encodes to 1011010 (checked on the next slide).

def encode(msg):                       # msg: 4 bits -> 7-bit codeword
    return [sum(m * g for m, g in zip(msg, col)) % 2 for col in zip(*G)]

def syndrome(r):
    return tuple(sum(h * x for h, x in zip(row, r)) % 2 for row in Hm)

# each single-bit error at position i has syndrome = column i of H
FIX = {tuple(col): i for i, col in enumerate(zip(*Hm))}

def decode(r):
    r = list(r)
    s = syndrome(r)
    if any(s):
        r[FIX[s]] ^= 1                 # flip the bit the syndrome points to
    return r[:4]                       # systematic: data is the first 4 bits

assert all(syndrome(encode(m)) == (0, 0, 0) for m in itertools.product([0, 1], repeat=4))
fixed = 0
for m in itertools.product([0, 1], repeat=4):
    c = encode(m)
    for i in range(7):
        r = c.copy(); r[i] ^= 1
        fixed += decode(r) == list(m)
print(f"corrected {fixed} of {16 * 7} single-bit errors")
print("codeword for 1011:", "".join(map(str, encode([1, 0, 1, 1]))))
codewords = ["".join(map(str, encode(m))) for m in itertools.product([0, 1], repeat=4)]
print("min distance:", min_distance(codewords))

corrected 112 of 112 single-bit errors
codeword for 1011: 1011010
min distance: 3

G and Hm are the matrices from the last slide, as lists of lists.

• Every one of the 16 messages × 7 error positions is fixed: 112 of 112.
• Decoding is one small matrix multiply plus one table lookup. No search.

wrong = 0
for m in itertools.product([0, 1], repeat=4):
    c = encode(m)
    for i, j in itertools.combinations(range(7), 2):
        r = c.copy(); r[i] ^= 1; r[j] ^= 1
        wrong += decode(r) != list(m)
print(f"double errors decoded wrongly: {wrong} of {16 * 21}")

double errors decoded wrongly: 336 of 336

Two errors always fool it: the syndrome points at a third bit, which makes 3 errors. A d = 3 code corrects 1 error, never 2. Adding an overall parity bit (the [8,4,4] code) turns this into "correct 1, detect 2", the SECDED used in ECC RAM.

Encode the message 0110

d₁ d₂ d₃ d₄ = 0 1 1 0. Add mod 2 (1 + 1 = 0):

• p₁ = d₁+d₂+d₄ = 0+1+0 = 1
• p₂ = d₁+d₃+d₄ = 0+1+0 = 1
• p₃ = d₂+d₃+d₄ = 1+1+0 = 0

Sent codeword: 0110 110.

The channel flips bit 3

Received: 0100 110. The decoder does not know which bit flipped.

Decode: compute the syndrome

Each row of H rechecks one circle. Add the received bits where the row has a 1:

• Row 1 (bits 1, 2, 4, 5): 0+1+0+1 = 0
• Row 2 (bits 1, 3, 4, 6): 0+0+0+1 = 1
• Row 3 (bits 2, 3, 4, 7): 1+0+0+0 = 1

Syndrome 011 matches column 3 of H. So flip bit 3 back and read 0110. Fixed.

def block_error(p, n, t):
    """P(more than t errors in n bits)."""
    return 1 - sum(math.comb(n, k) * p**k * (1 - p)**(n - k) for k in range(t + 1))
for p in [0.01, 0.05, 0.1]:
    print(f"p={p}: uncoded 4 bits fail {block_error(p, 4, 0):.4f}   "
          f"Hamming(7,4) fails {block_error(p, 7, 1):.4f}")

p=0.01: uncoded 4 bits fail 0.0394   Hamming(7,4) fails 0.0020
p=0.05: uncoded 4 bits fail 0.1855   Hamming(7,4) fails 0.0444
p=0.1: uncoded 4 bits fail 0.3439   Hamming(7,4) fails 0.1497

At p = 0.01, a 4-bit block fails 20× less often for a 75% bandwidth cost. Compare rep3: rate 1/3 for a similar gain per bit.

Hamming(7,4) gets close to rep3's reliability at 1.7× the rate. But both are far below capacity 0.919.

For a noisy channel, approaching the asymptotic capacity with vanishing error generally requires growing block length. Special noiseless or zero-error cases need separate treatment.

Proof. Radius-t balls around codewords are disjoint. Each holds ∑ C(n,i)(q−1)ⁱ words. They all fit in the qⁿ words of the space. ∎

Tietäväinen–van Lint (1973): over prime-power alphabets, every nontrivial perfect code has the parameters of a Hamming code or of one of the two Golay codes, [23,12,7] binary and [11,6,5] ternary.

"Trivial" perfect codes also exist: the whole space, a single word, and odd-length binary repetition codes. That is why rep(5,1) shows up as PERFECT below.

def hamming_bound_ok(n, k, d, q=2):
    t = (d - 1) // 2
    ball = sum(math.comb(n, i) * (q - 1) ** i for i in range(t + 1))
    return q ** k * ball, q ** n

for name, n, k, d in [("Hamming(7,4)", 7, 4, 3), ("Hamming(15,11)", 15, 11, 3),
                      ("Golay(23,12)", 23, 12, 7), ("rep(5,1)", 5, 1, 5),
                      ("parity(8,7)", 8, 7, 2), ("BCH(15,7)", 15, 7, 5)]:
    used, total = hamming_bound_ok(n, k, d)
    tag = "PERFECT" if used == total else f"{used / total:.3f} of space"
    print(f"{name:15} 2^k*|ball|={used:>8}  2^n={total:>8}  {tag}")

Hamming(7,4)    2^k*|ball|=     128  2^n=     128  PERFECT
Hamming(15,11)  2^k*|ball|=   32768  2^n=   32768  PERFECT
Golay(23,12)    2^k*|ball|= 8388608  2^n= 8388608  PERFECT
rep(5,1)        2^k*|ball|=      32  2^n=      32  PERFECT
parity(8,7)     2^k*|ball|=     128  2^n=     256  0.500 of space
BCH(15,7)       2^k*|ball|=   15488  2^n=   32768  0.473 of space

Golay: 2¹² (1 + 23 + 253 + 1771) = 2¹² · 2¹¹ = 2²³. A small miracle. Voyager 1 and 2 used the (extended) Golay code for color images of Jupiter and Saturn.

Proof. Delete the last d − 1 positions of every codeword. Two codewords differ in at least d places, so they still differ. We now have q^k distinct words of length n − d + 1. So q^k ≤ q^n−d+1. ∎

Binary MDS codes are trivial (repetition, parity, full space). You need a large alphabet, like bytes, for good MDS codes. That is what Reed–Solomon uses.

for name, n, k, d in [("Hamming(7,4)", 7, 4, 3), ("Golay(23,12)", 23, 12, 7),
                      ("RS(255,223)", 255, 223, 33), ("RS(7,3) over GF(8)", 7, 3, 5)]:
    print(f"{name:19} d={d:>2}  n-k+1={n - k + 1:>2}  {'MDS' if d == n - k + 1 else ''}".rstrip())

Hamming(7,4)        d= 3  n-k+1= 4
Golay(23,12)        d= 7  n-k+1=12
RS(255,223)         d=33  n-k+1=33  MDS
RS(7,3) over GF(8)  d= 5  n-k+1= 5  MDS

RS(255,223) with 32 check bytes corrects 16 byte errors or 32 erasures per block. It is the CCSDS deep-space standard.

• Two different polynomials of degree < k agree on at most k − 1 points. So codewords differ in at least n − k + 1 places: MDS.
• Erasures: any k surviving values fix f by Lagrange interpolation.
• Errors: Berlekamp–Massey or Berlekamp–Welch decode up to (n−k)/2 wrong symbols in polynomial time.
• Real systems use GF(2⁸), so a symbol is a byte. A burst of bad bits hits only a few bytes.

This prime-field example demonstrates Reed-Solomon evaluation and erasure recovery by interpolation. It does not implement an unknown-location error decoder such as Berlekamp-Massey.

P = 929                                    # prime field used by PDF417 barcodes

def rs_encode(msg, n):
    """Codeword = values of the message polynomial at x = 0..n-1 (mod P)."""
    return [sum(c * pow(x, i, P) for i, c in enumerate(msg)) % P for x in range(n)]

def rs_recover(points, k):
    """Lagrange interpolation from any k (x, y) pairs -> message coefficients."""
    coeffs = [0] * k
    for j, (xj, yj) in enumerate(points[:k]):
        basis, denom = [1], 1              # prod (x - xm) for m != j
        for m, (xm, _) in enumerate(points[:k]):
            if m != j:
                basis = [(a - xm * b) % P for a, b in zip([0] + basis, basis + [0])]
                denom = denom * (xj - xm) % P
        scale = yj * pow(denom, P - 2, P) % P
        coeffs = [(c + scale * b) % P for c, b in zip(coeffs, basis)]
    return coeffs

msg = [72, 105, 33]                        # "Hi!" as numbers, k = 3
cw = rs_encode(msg, 7)                     # n = 7, can lose any 4
print("codeword:", cw)
survivors = [(x, cw[x]) for x in (1, 4, 6)]    # 4 of 7 symbols erased
print("recovered:", rs_recover(survivors, 3))

codeword: [72, 210, 414, 684, 91, 493, 32]
recovered: [72, 105, 33]

Historical teaching excerpts may omit support definitions. No companion test file is required for the complete lab below: it includes its own implementation, input data, and checks.

• An error pattern E(x) is missed only if G(x) divides E(x).
• With r = 32: every single-bit error and every burst of length ≤ 32 is caught. Random corruption slips by with probability about 2⁻³².
• Used in Ethernet, ZIP, PNG, SATA. Detect, then ask again (ARQ) or drop.

Use a keyed MAC, digital signature, or a hash whose expected value is obtained through a trusted channel. An attacker can recompute an unauthenticated hash or CRC after changing data.

def crc32(data: bytes) -> int:
    crc = 0xFFFFFFFF
    for byte in data:
        crc ^= byte
        for _ in range(8):
            crc = (crc >> 1) ^ (0xEDB88320 if crc & 1 else 0)
    return crc ^ 0xFFFFFFFF

msg = b"The quick brown fox"
print(hex(crc32(msg)), hex(zlib.crc32(msg)))
bad = bytearray(msg); bad[4] ^= 0b00000100          # flip one bit
print(hex(crc32(bytes(bad))), crc32(bytes(bad)) != crc32(msg))
parity_bit = sum(bin(b).count("1") for b in msg) % 2
print("parity bit:", parity_bit)

0xb74574de 0xb74574de
0x3d3911bd True
parity bit: 1

Our bitwise CRC-32 matches zlib.crc32. 0xEDB88320 is the IEEE 802.3 polynomial with its bits reversed (the "reflected" form).

Historical teaching excerpts may omit support definitions. No companion test file is required for the complete lab below: it includes its own implementation, input data, and checks.

def polarize(eps, levels):
    z = [eps]
    for _ in range(levels):
        z = [w for e in z for w in (2 * e - e * e, e * e)]
    return z

for levels in [1, 4, 8, 12, 16]:
    z = polarize(0.5, levels)
    good = sum(e < 1e-3 for e in z) / len(z)
    bad = sum(e > 1 - 1e-3 for e in z) / len(z)
    print(f"N=2^{levels:<2}  good={good:.3f}  bad={bad:.3f}  middle={1 - good - bad:.3f}")

N=2^1   good=0.000  bad=0.000  middle=1.000
N=2^4   good=0.062  bad=0.062  middle=0.875
N=2^8   good=0.266  bad=0.266  middle=0.469
N=2^12  good=0.385  bad=0.385  middle=0.229
N=2^16  good=0.446  bad=0.446  middle=0.108

On BEC(0.5), erasure rates split into near 0 and near 1. The mean stays 0.5 (the test asserts it), so the good fraction tends to 1 − ε = C.

Sketch. I(X;Y) = h(Y) − h(Z). Variance of Y is at most P+N, and the Gaussian has the most differential entropy for a given variance: h = ½ log(2πeσ²). So I ≤ ½ log(2πe(P+N)) − ½ log(2πeN). A Gaussian input reaches it. ∎

SNR in dB is 10 log₁₀(P/N). Each +3 dB of SNR (at high SNR) buys about one more bit per sample pair, or B more bits/s.

def shannon_hartley(B_hz, snr_db):
    return B_hz * math.log2(1 + 10 ** (snr_db / 10))

for name, B, snr in [("phone line 3.1 kHz", 3.1e3, 35), ("Wi-Fi 20 MHz", 20e6, 25),
                     ("Wi-Fi 6 160 MHz", 160e6, 35), ("5G NR 100 MHz", 100e6, 20),
                     ("5G mmWave 400 MHz", 400e6, 15)]:
    print(f"{name:19} SNR={snr:>2} dB  C={shannon_hartley(B, snr) / 1e6:9.3f} Mbit/s")
print(f"Shannon limit Eb/N0 = ln 2 = {10 * math.log10(math.log(2)):.2f} dB")

phone line 3.1 kHz  SNR=35 dB  C=    0.036 Mbit/s
Wi-Fi 20 MHz        SNR=25 dB  C=  166.188 Mbit/s
Wi-Fi 6 160 MHz     SNR=35 dB  C= 1860.353 Mbit/s
5G NR 100 MHz       SNR=20 dB  C=  665.821 Mbit/s
5G mmWave 400 MHz   SNR=15 dB  C= 2011.123 Mbit/s
Shannon limit Eb/N0 = ln 2 = -1.59 dB

• Phone line: 36 kbit/s. V.34 modems hit 33.6 kbit/s, right at the wall.
• Wi-Fi 6, one stream, 160 MHz, 1024-QAM rate 5/6: 1201 Mbit/s, about 65% of this bound.
• Real links pass these numbers per antenna with MIMO: k spatial streams give up to k parallel channels.

Noise power grows with bandwidth: N = N₀B. So

C(B) = B log₂(1 + P / (N₀B)).

for B in [1e6, 10e6, 100e6, 1e9]:
    S_over_N0 = 1e8                        # fixed received power / noise density
    print(f"B={B/1e6:>6.0f} MHz  C={B * math.log2(1 + S_over_N0 / B) / 1e6:7.2f} Mbit/s")
print(f"limit B->inf: {1e8 / math.log(2) / 1e6:.2f} Mbit/s")

B=     1 MHz  C=   6.66 Mbit/s
B=    10 MHz  C=  34.59 Mbit/s
B=   100 MHz  C= 100.00 Mbit/s
B=  1000 MHz  C= 137.50 Mbit/s
limit B->inf: 144.27 Mbit/s

More bandwidth helps, but with a fixed power the gains flatten. As B → ∞, C → (P/N₀) log₂ e.

Below −1.59 dB no code, however long, can work. This number is the yardstick for every modern code: "how many dB from the Shannon limit?"

Power-limited links use low-rate codes and simple modulation. Bandwidth-limited links pack many bits per symbol (256- or 1024-QAM) and need high SNR.

1. What is the capacity of a binary erasure channel that erases 25% of bits?
2. What is the capacity of a BSC with p = 0.5? And with p = 1?
3. A code has minimum distance 5. How many errors can it correct? How many can it detect?
4. What is the Hamming distance between 1011010 and 0110110?
5. A 1 MHz channel has SNR = 15 (linear). What is its capacity?
6. You want rate 0.6 over a BSC with p = 0.11. Can any code make errors vanish?

Try each one before you look to the right.

Noise sets a speed limit, not a quality limit.

• Below C, long random-like codes make errors vanish. Above it, nothing can.
• Simple codes (repetition, Hamming) show the idea but sit far from C.
• Reed–Solomon owns bursts and erasures: QR codes, CDs, RAID, deep space.
• LDPC, turbo and polar codes reach within a fraction of a dB of Shannon's 1948 bound. Your phone uses them every second.

The big picture: source coding removes redundancy down to H. Channel coding adds it back, in a smart form, up to C. If H < C, the message gets through.