Integrating T-Shirt

Recently, my son spent a week at an Olympiad Maths training camp at Oxford University, delving into a treasure trove of challenging and enjoyable problems. At the end of the adventure, each kid got a T-shirt with a scary, complicated function printed on it. The moment I saw that shirt, I knew I wasn’t just going to admire it — I was going to write code to integrate it!

So where do I start?

Probably from the beginning: What the hell is this function?

At first glance, the formula on the T-shirt looks intimidating: integrals, expectations, limits, indicators, and a mysterious “mod 2”. But behind it lies a surprisingly elegant piece of probabilistic geometry.

Let’s unpack what’s really going on.

The setup: throwing random disks

Start with the unit square, randomly place disk centers inside the square and assign each disk a random radius. The key detail is that the radii shrink as grows, so the total expected disk area stays constant.

Get Michael Leznik’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

Here’s how that looks in code:

 1. import numpy as np
 2.  
 3. n = 500
 4. rng = np.random.default_rng()
 5.  
 6. # Disk centers uniformly distributed in [0,1]^2
 7. centers = rng.uniform(0, 1, size=(n, 2))
 8.  
 9. # Radii scaled so total expected area does not explode
10. r_max = np.sqrt(3 / (n * np.pi))
11. radii = rng.uniform(0, r_max, size=n)

As n increases, we get more disks — but each one is smaller.

Counting coverage — but only modulo 2

Pick a point x in the square and ask: How many disks cover this point? Instead of keeping the exact count, we only care whether the number is odd or even. Odd counts as 1; even counts as 0. This is what the “mod 2” in the formula means. To approximate this across the whole square, we sample many random points:

1. # Monte Carlo points used to approximate the integral
2. n_points = 200_000
3. points = rng.uniform(0, 1, size=(n_points, 2))

Now we track parity — not total coverage:

 1. # 0 = even coverage, 1 = odd coverage
 2. parity = np.zeros(n_points, dtype=np.int8)
 3.  
 4. for (cx, cy), r in zip(centers, radii):
 5.     dx = points[:, 0] - cx
 6.     dy = points[:, 1] - cy
 7.  
 8.     inside = (dx*dx + dy*dy) <= r*r
 9.  
10.     # XOR implements addition modulo 2
11.     parity ^= inside

Each disk flips the parity of the points it covers.

From parity to area

The integral in the formula measures the area of points with odd coverage. Because the square has area 1, this is simply the fraction of sampled points whose parity is 1:

1. odd_coverage_area = parity.mean()
2.

This single number is the value of the integral for one random disk configuration.

Why do we take an expectation?

Every random configuration gives a slightly different odd-coverage area.
To get something stable, we repeat the experiment and average.

1. def estimate_expected_area(n, trials=20):
2.     values = []
3.     for _ in range(trials):
4.         values.append(run_single_simulation(n))
5.     return np.mean(values)

This corresponds exactly to the expectation in the formula.

The limit

The limit outside the expectation is crucial.

As n grows:

• Disks become tiny
• Boundary effects vanish
• The geometry becomes locally Poisson

Numerically, we just increase n and observe convergence:

1. for n in [100, 200, 500, 1000, 2000]:
2.     estimate = estimate_expected_area(n)
3.     print(f"n={n:4d} → expected odd area ≈ {estimate:.4f}")

After running this code an unreasonably large but impressive 100,000 times, the entire expression collapses to ~0.475

1. import math
2.  
3. limit_value = (1 - math.exp(-3)) / 2
4. print(limit_value)

I suspect that some of the readers might integrate this function analytically, so please be my guests and drop your solution here.