A circular billiard puzzle by Xavier Durawa: Charlotte's strands reflect inside a circular frame; how many earlier strands does the n-th strand cross on average? A self-contained walk through the geometry, the modular-arithmetic reformulation, the piecewise integration, the antipode trick, and the telescoping sum that lands the exact closed form.
✠
A puzzle by Xavier Durawa for the week of 5/3.
The puzzle
Charlotte is building a web inside a circular frame. She picks a starting point on the boundary uniformly at random, picks a direction pointing into the disk uniformly at random, walks in a straight line, and lays down a strand of silk. When she reaches the boundary she reflects elastically, in exactly the way a billiard ball or a beam of light would, and lays down the next strand. New strands cross many of the earlier ones. On average, how many of the previous strands does the n-th strand cross?
That is the whole question. Getting to a clean answer takes about ten small steps, and we will spell out every one of them. The road runs through some elementary geometry of reflection, a change of variables that collapses two random choices into one, a piecewise integration, a symmetry argument with a satisfying twist, and a telescoping sum. The reward at the end is an exact closed form, simple enough to fit on a single line.
We work in the unit disk {(x,y):x2+y2≤1}. Charlotte’s starting point is on the boundary circle and her initial direction is some unit vector pointing inward.
Step 1: every chord has the same length
Here is the first surprise. Every strand Charlotte lays down has exactly the same length, no matter where the bounces occur. The reflection law arranges this for free.
To see why, look at the very first bounce. Charlotte travels along a chord from A to B, reflects, and continues along BC. We claim ∣AB∣=∣BC∣.
Let O be the centre of the disk. Triangle OAB is isosceles, since the two sides OA and OB are radii of equal length. A standard fact about isosceles triangles is that their two base angles (the angles opposite the two equal sides) are equal. So
∠OAB=∠OBA.
Call this common angle α.
Now we need one fact from elementary plane geometry: at any point of a circle, the radius is perpendicular to the tangent line. Quick justification: the tangent line at B touches the circle only at B, and any other point on the tangent line is outside the circle, i.e. at distance greater than 1 from O. So B is the closest point on the tangent line to O. The closest point on a line to a fixed external point is always the foot of the perpendicular dropped from that point, so OB must meet the tangent at B at a right angle.
Back to the bounce. The chord BA leaves B along a direction making angle α with the radius OB, so it makes angle 90°−α with the tangent. The reflection law says the angle of incidence equals the angle of reflection, both measured against the tangent. So the outgoing chord BC also makes angle 90°−α with the tangent, just on the opposite side. That places BC at angle α from the radius OB on the opposite side. So ∠OBC=α as well.
We have arrived at two isosceles triangles with the same shape. Triangle OAB has OA=OB=1 and base angles α; its apex angle at O is therefore 180°−2α (since angles in a triangle sum to 180°). Triangle OBC has OB=OC=1 and base angle α at B, hence by isosceles symmetry base angle α at C too, and apex angle at O also 180°−2α.
The two triangles each have two unit sides separated by the same included angle, 180°−2α. By side-angle-side, they are congruent. In particular the third sides match: ∣BC∣=∣AB∣.
So every chord has the same length. Equivalently (by the law of cosines applied to OAB, which gives ∣AB∣2=2−2cos∠AOB), every chord subtends the same central angle. Call this constant central angle δ. The boundary points
P0,P1,P2,P3,…
where successive strands meet the circle therefore form an arithmetic progression of angles: if P0 sits at angle θ0, then Pk sits at angle θ0+kδ.
Step 2: the rotation δ is uniform on (0,2π)
The puzzle hands us two pieces of randomness: a uniform starting point and a uniform inward direction. We have just seen that the entire web is determined by the single number δ. So the next question is, what is the distribution of δ?
The starting point itself does no work. The whole problem is rotationally symmetric, so we can place P0 at (1,0) once and forever, and forget about it.
What is left is the random direction. Parametrise inward-pointing unit vectors at (1,0) by an angle ϑ∈(0,π) measured from the upward tangent. The tangent at (1,0) points in the y-direction, so the unit vector
v(ϑ)=(−sinϑ,cosϑ)
points along the upward tangent at ϑ=0, along the inward normal (−1,0) at ϑ=π/2, and along the downward tangent at ϑ=π. The puzzle specifies the direction as uniform, so ϑ is uniform on (0,π). (No Bertrand-style ambiguity here: the random ingredient is the direction at a fixed boundary point, not the chord.)
Charlotte then walks along the line
P(t)=(1,0)+tv(ϑ)=(1−tsinϑ,tcosϑ),t≥0.
She bounces when she meets the boundary again, i.e. when ∥P(t)∥2=1:
The root t=0 is the starting point. The other root, t=2sinϑ, is the next bounce. Substitute it back:
P1=(1−2sin2ϑ,2sinϑcosϑ)=(cos2ϑ,sin2ϑ),
where the last step uses the double-angle identities cos2ϑ=1−2sin2ϑ and sin2ϑ=2sinϑcosϑ.
So P1 sits on the unit circle at angle 2ϑ. The central angle from P0 to P1 is δ=2ϑ.
Since ϑ is uniform on (0,π), the rotation δ=2ϑ is uniform on (0,2π). Setting
β=2πδ,
we have a single random parameter β, uniform on (0,1), that determines the entire web. The whole problem is now parametrised by one number on the unit interval.
Step 3: boundary positions on a circle of length 1
It will pay to forget about angles entirely and measure boundary positions as fractions of the circumference. The boundary becomes the interval [0,1) with the two endpoints glued together. Place P0 at position 0. The bounce point Pk is then at position
kβmod1=the fractional part of kβ.
A note on the notation: for any real x, the symbol xmod1 means the unique y∈[0,1) with x−y∈Z. So 1.7mod1=0.7, 2.5mod1=0.5, −0.3mod1=0.7.
The k-th strand (for k≥1) is the chord whose two boundary endpoints are (k−1)βmod1 and kβmod1. That is the picture we will work with from now on: a sequence of points marching around a length-1 circle in steps of β, with chords joining consecutive points.
Step 4: when do two chords cross inside the disk?
Two chords cross strictly inside the disk if and only if their four endpoints alternate around the boundary. The chord ab separates the open disk into two pieces; a second chord crosses ab iff its two endpoints lie on opposite arcs.
Alternation is invariant under rotating the boundary, so we can shift strand i to sit at {0,β}. For strand j with m=j−i≥1, the endpoints become
{u,u+βmod1},u=mβmod1.
The chord {0,β} cuts the boundary into two arcs: A=(0,β) of length β, and B=(β,1) of length 1−β. Strands i and j cross iff exactly one of {u,u+βmod1} lies in each arc.
Case β<1/2 (top row). Arc A has length β and is the shorter one. Three sub-cases for u:
u∈(0,β): u∈A and u+β∈(β,2β)⊂B. Cross.
u∈(β,1−β): u∈B and u+β∈(2β,1)⊂B. No cross.
u∈(1−β,1): u∈B and u+β−1∈(0,β)=A. Cross.
The cross region is u∈(0,β)∪(1−β,1).
Case β>1/2 (bottom row). Now arc B has length 1−β and is the shorter one. Three sub-cases:
u∈(0,1−β): u∈A and u+β∈(β,1)=B. Cross.
u∈(1−β,β): u∈A and u+β−1∈(0,2β−1)⊂A. No cross.
u∈(β,1): u∈B and u+β−1∈(2β−1,β)⊂A. Cross.
The cross region is u∈(0,1−β)∪(β,1).
Both cases together. In each case, the cross region consists of those u whose circle-distance to the nearer of {0,1} is less than the length of the shorter arc. Writing the circle distance as
d(x)=min(x,1−x),x∈[0,1),
the shorter-arc length is d(β)=min(β,1−β), and the criterion collapses to
strands i and j cross⟺d(mβmod1)<d(β),m=j−i.
The condition has nothing to do with i or j separately, only on the gap m between their indices. That is the simplification we need.
Step 5: the answer is a sum of probabilities
For each m≥1, define
Pm=P(d(mβmod1)<d(β)),β∼Unif(0,1).
By Step 4, Pm is the probability that any two strands whose indices differ by m cross.
Let Nn be the number of previous strands that strand n crosses. For each k=1,…,n−1, strand n versus strand k has m=n−k, so by linearity of expectation (which holds without any independence assumptions),
E[Nn]=∑k=1n−1P(strand n crosses strand k)=∑m=1n−1Pm.
The puzzle has reduced to a single problem: compute Pm.
Step 6: a worked example, P2
Let us compute P2=P(d(2βmod1)<d(β)) by hand. The trick is to chop [0,1] into pieces on which both d(β) and d(2βmod1) are linear in β, and read off the answer piece by piece.
The function d(β) is a tent: it climbs linearly from 0 at β=0 up to 1/2 at β=1/2, then descends to 0 at β=1. Explicitly, d(β)=β on [0,1/2], and d(β)=1−β on [1/2,1].
The function d(2βmod1) is the same shape, but at twice the frequency: two tents stacked side by side. A direct computation on each quarter-interval gives
The probability we want is the total length of the β-interval on which the orange tents (the high-frequency function) are below the green tent (the low-frequency one). Reading off each piece:
Adding the four contributions, P2=0+61+61+0=31.
The same recipe works for general m: the function d(mβmod1) has 2m linear pieces, and on each piece the comparison with d(β) is a single linear inequality whose solution set is an interval (or empty). Machine arithmetic with exact rationals can carry the calculation out to whatever m we want.
Step 7: a table for Pm and a closed form
Carrying out exactly that integration for m=1,2,…,12 in exact rational arithmetic (the code is at the end of this article) gives
P1=0,P2=31,P3=21,P4=157,P5=21,P6=3517,P7=21,P8=6331,P9=21,…
Two patterns leap off the page. Every odd m≥3 gives Pm=1/2 exactly. The even values P2k fit the formula
P2k=4k2−12k2−1.
Quick check: P2=(2−1)/(4−1)=1/3 at k=1; P4=(8−1)/(16−1)=7/15 at k=2; P6=(18−1)/(36−1)=17/35 at k=3. The pattern matches.
Both patterns deserve a proof. The first has a slick one. The second needs more work.
Why Pm=1/2 for every odd m≥3
Here is the nice symmetry. The integral defining Pm is invariant under the substitution β↦β′, where
β′:=(β+21)mod1
is the antipode of β on our length-1 circle. This substitution is a measure-preserving bijection of [0,1) to itself, so any probability over a uniform β is the same as the corresponding probability over a uniform β′.
What does the substitution do to our two circle distances?
First, d(β′). Shifting by 1/2 moves a point to its antipode. The antipode of a point at circle-distance d(β) from 0 sits at circle-distance 21−d(β) from 0. (Check: if β=0.1, then β′=0.6, and d(β′)=1−0.6=0.4=0.5−0.1=0.5−d(β). If β=0.7, then β′=0.2, and d(β′)=0.2=0.5−0.3=0.5−d(β).) So
d(β′)=21−d(β).
Second, d(mβ′mod1). Compute:
mβ′=mβ+2m=mβ+2m−1+21.
For odd m, the term (m−1)/2 is an integer. Reducing modulo 1 kills the integer:
mβ′mod1=(mβ+21)mod1.
So mβ′mod1 is the antipode of mβmod1, and by the same antipode rule,
d(mβ′mod1)=21−d(mβmod1).
Now substitute everything into the inequality d(mβmod1)<d(β). After substitution it becomes
21−d(mβmod1)<21−d(β),
which simplifies to
d(mβmod1)>d(β).
The substitution flipped the direction of the inequality. So
Pm=P(d(mβmod1)<d(β))=P(d(mβmod1)>d(β)).
The equality case d(mβmod1)=d(β) has zero probability (it occurs on a finite set of β-values, which has Lebesgue measure zero), so the two complementary events have probabilities summing to 1. Each must be 1/2.
Why P2k=4k2−12k2−1
The slick substitution above fails for even m, because the term (m−1)/2 is no longer an integer and so does not vanish modulo 1. We have to roll up our sleeves and integrate. The same piecewise-linear approach we used by hand for m=2 goes through for general m=2k; it is just longer.
Halve the interval by symmetry. The substitution β↦1−β is a measure-preserving bijection of [0,1) to itself. It sends β to its mirror image about 1/2. Under it, d(β)=min(β,1−β) is unchanged (since min(1−β,β)=min(β,1−β)), and similarly d(2kβmod1) is unchanged (since d is symmetric: d(1−x)=d(x), and 2k(1−β)mod1=(−2kβ)mod1). So the integral over [0,1] is twice the integral over [0,1/2]. We compute the latter and double.
The structure of d(2kβmod1) on [0,1/2]. On [0,1/2], d(β)=β, the simpler of the two functions. The function h(β):=d(2kβmod1) has period 1/(2k) in β (because adding 1/(2k) to β shifts 2kβ by 1, which is killed by mod1). The interval [0,1/2] contains exactly k full periods, indexed by p=0,1,…,k−1.
The p-th period is the interval [2kp,2kp+1], on which h traces out a tent of height 1/2. Within this period, the tent’s left half is the interval [2kp,4k2p+1], where h rises linearly from 0 to 1/2 with slope 2k:
h(β)=2kβ−pon the left half of the p-th period.
The right half is [4k2p+1,2kp+1], where h falls from 1/2 to 0 with slope −2k:
h(β)=(p+1)−2kβon the right half of the p-th period.
The diagonal y=β has slope 1, while the tent has slope ±2k. The tent crosses y=β once on each half, so the cross region inside the period is two short intervals at the ends.
Integrating piece by piece. We want, on each period, the set of β where h(β)<β.
Left half of period p. Solve 2kβ−p<β:
(2k−1)β<p,i.e.,β<2k−1p.
For p=0 this gives β<0, empty. For p≥1 we intersect with the left-half interval [2kp,4k2p+1]. Note 2kp<2k−1p (denominators), so the lower bound is p/(2k); and one can check 2k−1p≤4k2p+1 holds for p≤k−1 (cross-multiplying gives 4kp≤(2p+1)(2k−1)=4kp+2k−2p−1, equivalent to 2k−2p−1≥0, which is true for p≤k−1). So the upper bound is p/(2k−1), and the contribution from the left half of period p is
2k−1p−2kp=(2k−1)(2k)p⋅(2k−(2k−1))=(2k−1)(2k)p.
Right half of period p. Solve (p+1)−2kβ<β:
(2k+1)β>p+1,i.e.,β>2k+1p+1.
Intersect with the right-half interval [4k2p+1,2kp+1]. By a similar denominator check, 4k2p+1<2k+1p+1 (cross-multiply: (2p+1)(2k+1)=4kp+2k+2p+1 vs. 4k(p+1)=4kp+4k, difference 4k−2k−2p−1=2k−2p−1≥0 for p≤k−1); and 2kp+1≤21 for p≤k−1 (with equality at p=k−1, where the right end of the period is exactly 1/2). So the contribution from the right half of period p is
2kp+1−2k+1p+1=(2k)(2k+1)(p+1)⋅((2k+1)−2k)=(2k)(2k+1)p+1.
Sum the periods. Adding both halves over p=0,1,…,k−1:
measure on [0,21]=p=0∑k−1(2k−1)(2k)p+p=0∑k−1(2k)(2k+1)p+1=(2k−1)(2k)1p=0∑k−1p+(2k)(2k+1)1p=0∑k−1(p+1)=(2k−1)(2k)1⋅2(k−1)k+(2k)(2k+1)1⋅2k(k+1)=4(2k−1)k−1+4(2k+1)k+1.
We used ∑p=0k−1p=(k−1)k/2 and ∑p=0k−1(p+1)=∑q=1kq=k(k+1)/2. Doubling, by the symmetry we opened with,
P2k=2(2k−1)k−1+2(2k+1)k+1.
Putting it on a common denominator 2(2k−1)(2k+1)=2(4k2−1), the numerator becomes
(k−1)(2k+1)+(k+1)(2k−1).
Expand each piece: (k−1)(2k+1)=2k2+k−2k−1=2k2−k−1, and (k+1)(2k−1)=2k2−k+2k−1=2k2+k−1. Adding, (2k2−k−1)+(2k2+k−1)=4k2−2. So
P2k=2(4k2−1)4k2−2=4k2−12k2−1.
Sanity check. At k=1, 1/3, matching P2. At k=2, 7/15, matching P4. At k=3, 17/35, matching P6.
Step 8: summing the Pm‘s
We now want to compute E[Nn]=∑m=1n−1Pm. Split the sum by the parity of m. Among the values m=1,2,…,n−1:
m=1 contributes 0.
The odd values m=3,5,… each contribute 1/2. There are ⌊n/2⌋−1 of them (the odd numbers strictly between 1 and n).
The even values m=2,4,…,2K each contribute P2k=(2k2−1)/(4k2−1). Here K=⌊(n−1)/2⌋, the number of even values among {1,2,…,n−1}.
So
E[Nn]=0+(⌊n/2⌋−1)⋅21+∑k=1K4k2−12k2−1.
The even-m summand has a tidy rewrite:
4k2−12k2−1=21−2(4k2−1)1.
Check: 21−2(4k2−1)1=2(4k2−1)(4k2−1)−1=2(4k2−1)4k2−2=4k2−12k2−1. ✓
Substituting,
E[Nn]=(⌊n/2⌋−1)⋅21+2K−21∑k=1K4k2−11.
The first two terms are just half the count of all m‘s that contribute 1/2 each. The total count of contributors of either parity is n−2 (we have n−1 values of m, of which one, namely m=1, contributes 0). So
(⌊n/2⌋−1)+K=n−2,
which means
(⌊n/2⌋−1)⋅21+2K=2n−2.
A small sanity check at n=5: odd m∈{3} has count 1; even m∈{2,4} has count K=2; total 1+2=3=n−2. ✓
So
E[Nn]=2n−2−21∑k=1K4k2−11,K=⌊(n−1)/2⌋.
The only mystery left is the small tail sum.
Step 9: telescope the sum
The denominator factors:
4k2−1=(2k−1)(2k+1).
Partial fractions then give
(2k−1)(2k+1)1=21(2k−11−2k+11).Quick check. The right side equals 21⋅(2k−1)(2k+1)(2k+1)−(2k−1)=21⋅(2k−1)(2k+1)2=(2k−1)(2k+1)1. ✓
The point of partial fractions here is that the sum now telescopes: each negative term in one summand cancels the positive term in the next. Writing the first few out:
The middle line is where the work happens: −1/3 from the first parenthesis kills the +1/3 from the second, −1/5 kills +1/5, and so on down the line. Only the very first +1 and the very last −1/(2K+1) survive.
Substituting back into Step 8,
E[Nn]=2n−2−21(21−2(2K+1)1)=2n−2−41+4(2K+1)1=42n−5+4(2K+1)1.
Step 10: the exact closed form
The final job is to simplify 2K+1 in terms of n. Recall K=⌊(n−1)/2⌋:
If n is odd, n−1 is even, K=(n−1)/2, so 2K+1=n.
If n is even, n−1 is odd, K=(n−2)/2, so 2K+1=n−1.
Define n∗=n for n odd, and n∗=n−1 for n even. We have proved
E[Nn]=42n−5+4n∗1.
That is the answer. Two things are worth noticing.
First, as n→∞ the residual 1/(4n∗) vanishes and E[Nn]→(2n−5)/4=2n−45. A plausible-looking guess might be that the n-th strand crosses “half of the previous ones,” i.e. (n−1)/2=2n−21. The true answer undershoots that guess by exactly 43 in the limit. The discrepancy comes from the fact that consecutive strands (m=1) share a boundary endpoint and so never cross at all, which systematically depresses the count.
Second, the residual 1/(4n∗) decays like 1/(4n) but does so in pairs of equal height: n=3 and n=4 both give residual 1/12; n=5 and n=6 both give 1/20; and so on. Going from an odd n to the next n+1 does not move n∗ at all, so the residual stays put.
Small-n table
n
E[Nn] exact
decimal
(2n−5)/4
3
1/3
0.3333
0.2500
4
5/6
0.8333
0.7500
5
13/10
1.3000
1.2500
10
34/9
3.7778
3.7500
20
333/38
8.7632
8.7500
30
399/29
13.7586
13.7500
50
1164/49
23.7551
23.7500
100
9653/198
48.7525
48.7500
Monte Carlo verification
Each dot is the average of 80,000 random β draws. For each draw we build the boundary positions {0,β,2β,…,100β}mod1 and, for each strand n, count by direct geometry how many of the previous n−1 strands it crosses. The dots sit on the closed-form curve all the way out to n=100. The inset shows the exact residual E[Nn]−(2n−5)/4=1/(4n∗), which is already about 0.0025 at n=100 and decays as 1/(4n) in pairs of equal height.
Python code
import numpy as npfrom fractions import Fractiondef positions(n, beta): """Boundary positions of P_0, P_1, ..., P_n in [0, 1).""" return (np.arange(n + 1) * beta) % 1.0def to_xy(pos): """Map a position in [0, 1) to its (x, y) on the unit circle.""" a = 2 * np.pi * pos return np.column_stack([np.cos(a), np.sin(a)])def segments_cross(p1, p2, p3, p4, eps=1e-12): """Open segments p1-p2 and p3-p4 cross strictly inside the disk?""" d1, d2 = p2 - p1, p4 - p3 denom = d1[0]*d2[1] - d1[1]*d2[0] if abs(denom) < eps: return False diff = p3 - p1 t = (diff[0]*d2[1] - diff[1]*d2[0]) / denom s = (diff[0]*d1[1] - diff[1]*d1[0]) / denom return eps < t < 1-eps and eps < s < 1-epsdef E_Nn_montecarlo(n_max=100, n_trials=80_000, seed=0): """Estimate E[N_n] for n = 2..n_max by Monte Carlo.""" rng = np.random.default_rng(seed) counts = np.zeros(n_max + 1) for _ in range(n_trials): beta = rng.uniform(0, 1) pts = to_xy(positions(n_max, beta)) for n in range(2, n_max + 1): a, b = pts[n-1], pts[n] counts[n] += sum( segments_cross(a, b, pts[k-1], pts[k]) for k in range(1, n) ) return counts / n_trialsdef exact_E_Nn(n): """Closed form: (2n-5)/4 + 1/(4 n_*) where n_* = n if n is odd, else n - 1.""" n_star = n if n % 2 == 1 else n - 1 return Fraction(2*n - 5, 4) + Fraction(1, 4 * n_star)