Problem 395: Pythagorean Tree

Detailed mathematical approach

Problem Summary

Start from a unit square and call its bottom side the base. On the side opposite the base build a right triangle whose hypotenuse is that side and whose sides are in the ratio \(3:4:5\), with the smaller side on the right with respect to the base. Attach a square to each leg, then repeat the construction on both new squares, each time taking as the base the side touching the triangle. After infinitely many iterations the figure is the Pythagorean tree. At least one rectangle with sides parallel to the largest square encloses it completely, and we want the smallest possible area of such a rectangle, rounded to \(10\) decimal places. ^[1] It is the scalene cousin of the classical right-isosceles Pythagoras tree, which fits inside a \(6\times4\) box. ^[2]

The whole solution in one sentence: the tree is the unit square together with two shrunken, rotated copies of the whole tree, so each of its four extreme coordinates obeys a maximum-recursion that we evaluate by always expanding the most promising copy and bounding everything still unexplored.

Mathematical Approach

1. One step in coordinates

Put the first square at \(S=[0,1]^2\), so its base runs from \((0,0)\) to \((1,0)\) and the opposite side from \((0,1)\) to \((1,1)\). That side is the hypotenuse and has length \(1\), so the legs are \(4/5\) and \(3/5\) long. Let \(P=(x,y)\) be the right-angle vertex, above the side; the longer leg meets the left end \((0,1)\) and the shorter the right end \((1,1)\), so

$$x^2+(y-1)^2=\tfrac{16}{25},\qquad (x-1)^2+(y-1)^2=\tfrac{9}{25}.$$

Subtracting gives \(2x-1=7/25\), so \(x=16/25\), and then \((y-1)^2=144/625\) with \(y\gt1\), so \(P=(16/25,\,37/25)\). The two legs, as vectors \(a=P-(0,1)\) and \(b=(1,1)-P\), are

$$a=\left(\tfrac{16}{25},\tfrac{12}{25}\right),\qquad b=\left(\tfrac{9}{25},-\tfrac{12}{25}\right).$$

In exact arithmetic \(|a|=4/5\), \(|b|=3/5\) and \(a\cdot b=0\), as the right angle at \(P\) demands. Writing \(i(x,y)=(-y,x)\) for the quarter turn, one more exact identity matters later: \(b=-\tfrac34\,i\,a\).

2. The tree is two copies of itself

Identify \((x,y)\) with the complex number \(x+iy\); then \(a=\frac{16+12i}{25}\) and \(b=\frac{9-12i}{25}\). Multiplying by a complex number of modulus \(r\) rotates and scales by \(r\), so the two maps

$$f_1(z)=az+t_1,\qquad f_2(z)=bz+t_2,$$

with \(t_1=i\) and \(t_2=\frac{16+37i}{25}=P\), are similarities with ratios \(4/5\) and \(3/5\). They carry the base of the first square onto the two legs: \(f_1(0)=i\) and \(f_1(1)=t_1+a=P\), while \(f_2(0)=P\) and \(f_2(1)=t_2+b=1+i\). A similarity maps a square to a square, and the construction repeats verbatim inside each child, so the tree \(K\) satisfies

$$K=S\cup f_1(K)\cup f_2(K).$$

We let \(K\) contain the limit points too, which changes nothing about the closed rectangles enclosing it. The equation is an inhomogeneous version of the fixed-point equation of an iterated function system. ^[3] Let \(d_H(X,Y)\) be the Hausdorff distance, the largest distance from a point of either set to the nearest point of the other; the non-empty compact sets form a complete space under it, and for \(T(X)=S\cup f_1(X)\cup f_2(X)\),

$$d_H\bigl(T(X),T(Y)\bigr)\le\tfrac45\,d_H(X,Y),$$

because a union with the fixed \(S\) cannot increase the distance and each \(f_j\) scales it by \(|a|\) or \(|b|\). A contraction on a complete space has one fixed point and its iterates converge to it, so \(K\) is the unique compact solution and the figure after \(n\) iterations, \(K_n=T^n(S)\), tends to it.

3. The smallest rectangle is the bounding box

Because \(K\) is compact, the extreme coordinates \(x_{\min}\), \(x_{\max}\), \(y_{\min}\), \(y_{\max}\) are attained. Every axis-parallel rectangle containing \(K\) contains the four points where they occur, so its width is at least \(x_{\max}-x_{\min}\) and its height at least \(y_{\max}-y_{\min}\) — and \([x_{\min},x_{\max}]\times[y_{\min},y_{\max}]\) is itself such a rectangle. The answer is the product of those two lengths; nothing else about the tree is needed.

4. The tree fits inside a disc of radius 5

Let \(\rho=\max_{z\in K}|z|\), finite because \(K\) is compact. The three pieces of \(K\) give

$$\rho\le\max\left(\sqrt2,\ 1+\tfrac45\rho,\ \tfrac{\sqrt{65}}5+\tfrac35\rho\right),$$

since, writing \(c_1=a\) and \(c_2=b\), we have \(|f_j(z)|\le|c_j||z|+|t_j|\) with \(|t_1|=1\) and \(|t_2|=\sqrt{65}/5\approx1.6125\), while \(\sqrt2\) is the largest modulus on \(S\). If \(\rho\gt5\), all three entries are below \(\rho\): \(\sqrt2\lt5\); \(1+\frac45\rho\lt\rho\) exactly when \(\rho\gt5\); and \(\frac{\sqrt{65}}5+\frac35\rho\lt\rho\) exactly when \(\rho\gt\sqrt{65}/2\approx4.031\). That contradiction forces \(\rho\le5\), and a search for the farthest point gives \(\rho\approx4.3290161979\). Hence for any similarity \(f(z)=cz+t\) built from these maps,

$$f(K)\subseteq\{\,z:|z-t|\le5|c|\,\}.$$

The same two moduli give \(|t_j|+5|c_j|\le5\), since \(1+4=5\) and \(\frac{\sqrt{65}}5+3\approx4.612\); Section 10 needs that.

5. Why enumerating squares cannot work

Because \(0\in S\subseteq K\), Section 2 gives

$$d_H(K_n,K)\le(4/5)^n\,d_H(S,K)\le5\,(4/5)^n,$$

so each extreme coordinate of \(K_n\) is within \(\delta_n=5\,(4/5)^n\) of the true one. As \(K_n\subseteq K\), each side of its box is short by at most \(2\delta_n\), and both sides are at most \(10\); with \(A_n\) and \(A\) for the two areas,

$$0\le A-A_n\le 10\cdot2\delta_n+10\cdot2\delta_n=200\,(4/5)^n.$$

This bound first falls below \(5\cdot10^{-11}\) at \(n=131\). The measured deficit is smaller but no more encouraging: \(0.0904\) at depth \(20\) and \(0.0240\) at depth \(24\), shrinking by only about \(0.72\) per level, which extrapolates to depth \(85\) — some \(2^{85}\approx4\cdot10^{25}\) squares in the last level, where depth \(24\) already has \(2^{24}=16\,777\,216\). The recursion has to be solved, not unrolled.

6. The support function

For a direction \(d\) and a non-empty compact set \(X\), the support function is ^[4]

$$h_X(d)=\max_{z\in X}\langle d,z\rangle.$$

It is the tool we need, because \(x_{\max}=h_K(1,0)\), \(x_{\min}=-h_K(-1,0)\), \(y_{\max}=h_K(0,1)\) and \(y_{\min}=-h_K(0,-1)\). Two rules follow from the definition. The maximum over a union is the larger of the two maxima:

$$h_{X\cup Y}(d)=\max\bigl(h_X(d),h_Y(d)\bigr).$$

For an affine image \(g(z)=Mz+t\) with \(M\) linear, \(\langle d,Mz+t\rangle=\langle d,t\rangle+\langle M^{\mathsf T}d,z\rangle\), and taking the maximum over \(z\in X\) gives

$$h_{g(X)}(d)=\langle d,t\rangle+h_X\bigl(M^{\mathsf T}d\bigr).$$

Note also that \(h_X\) is positively homogeneous, \(h_X(\lambda d)=\lambda h_X(d)\) for \(\lambda\ge0\), and that a linear maximum over \(X\) equals the one over its convex hull — the computation never sees the gaps and overlaps inside the figure.

7. The recursion for the support values

Multiplication by \(c=u+iv\) has the matrix \(\begin{pmatrix}u&-v\\v&u\end{pmatrix}\), whose transpose is \(\begin{pmatrix}u&v\\-v&u\end{pmatrix}\): that is multiplication by the conjugate \(\overline c\). So \(M^{\mathsf T}d=\overline c\,d\) and, in particular, \(|M^{\mathsf T}d|=|c|\,|d|\). Applying the two rules of Section 6 to \(K=S\cup f_1(K)\cup f_2(K)\) gives a single scalar equation for \(h=h_K\):

$$\begin{aligned}h(d)=\max\bigl(&h_S(d),\\&\langle d,t_1\rangle+h(\overline a\,d),\\&\langle d,t_2\rangle+h(\overline b\,d)\bigr).\end{aligned}$$

The maximum of a linear function over a square is attained at a corner, and \(S\) has the corners \(0\), \(1\), \(i\) and \(1+i\), so

$$h_S(d)=\max\left(0,\ d_x,\ d_y,\ d_x+d_y\right).$$

8. One level worked out by hand

Replace the inner \(h\) in Section 7 by \(h_S\) and the recursion returns the extremes of the figure after one iteration. Take \(d=(0,1)\). Here \(h_S(d)=1\); the left branch has \(\langle d,t_1\rangle=1\) and \(\overline a\,d=\left(\tfrac{12}{25},\tfrac{16}{25}\right)\) with \(h_S=\tfrac{28}{25}\), giving \(\tfrac{53}{25}\); the right branch has \(\langle d,t_2\rangle=\tfrac{37}{25}\) and \(\overline b\,d=\left(-\tfrac{12}{25},\tfrac9{25}\right)\) with \(h_S=\tfrac9{25}\), giving \(\tfrac{46}{25}\). The maximum is \(\tfrac{53}{25}=2.12\). The other three directions come out as

• \(d=(1,0)\): \(\max\left(1,\,0+\tfrac{16}{25},\,\tfrac{16}{25}+\tfrac{21}{25}\right)=\tfrac{37}{25}\);
• \(d=(-1,0)\): \(\max\left(0,\,0+\tfrac{12}{25},\,-\tfrac{16}{25}+0\right)=\tfrac{12}{25}\);
• \(d=(0,-1)\): \(\max\left(0,\,-1+0,\,-\tfrac{37}{25}+\tfrac{12}{25}\right)=0\).

So the depth-\(1\) figure has the exact bounding box \(\left[-\tfrac{12}{25},\tfrac{37}{25}\right]\times\left[0,\tfrac{53}{25}\right]\) of area \(\tfrac{2597}{625}=4.1552\), which is what the C++ checkpoint verifies. Two, three and four levels deep the exact areas are \(8.20207616\), \(13.0035093504\) and \(16.63059765362688\) — still nowhere near the limit, as Section 5 predicts.

9. Unrolling the recursion into words

For a finite word \(w=w_1w_2\ldots w_n\) over \(\{1,2\}\) put \(f_w=f_{w_1}\circ\cdots\circ f_{w_n}\), again a similarity, say \(f_w(z)=c_wz+t_w\) with \(c_w\) the product of the corresponding factors, so \(|c_w|=(4/5)^k(3/5)^{n-k}\) for \(k\) the number of \(1\)s in \(w\). As \(f_w(S)\) is the square of the word \(w\), expanding Section 7 repeatedly turns the recursion into one maximum over all words,

$$h(d)=\sup_w\bigl(s_w+h_S(d_w)\bigr),\quad s_w=\langle d,t_w\rangle,\ \ d_w=\overline{c_w}\,d.$$

Appending one symbol is cheap: from \(t_{wj}=f_w(t_j)=c_wt_j+t_w\) we get

$$s_{wj}=s_w+\langle d_w,t_j\rangle,\qquad d_{wj}=\overline{c_j}\,d_w,$$

so a node stores only the scalar \(s_w\) and the vector \(d_w\), never the word. Section 4 bounds everything below \(w\), since \(f_w(K)\) lies in the disc of radius \(5|c_w|\) about \(t_w\):

$$h_{f_w(K)}(d)\le s_w+5\,|d_w|.$$

The square \(f_w(S)\) alone already gives the lower bound \(s_w+h_S(d_w)\). Since \(s_w\) changes as the word grows, the upper bound need not fall — along the extremal branch it tends to the support value. What falls geometrically is the gap between the two bounds, at most \(5|d_w|\) because \(h_S(d_w)\ge0\), and each symbol multiplies \(|d_w|\) by \(4/5\) or \(3/5\).

10. Best-first branch and bound

The two bounds of Section 9 give a classical branch-and-bound search. ^[5] Start with \(\beta=h_S(d)\) and a max-heap holding one node, the empty word, with offset \(0\), direction \(d\) and bound \(5|d|\). Then look at the node on top, the one with the largest bound: stop if that bound is at most \(\beta+\varepsilon\); otherwise remove it, replace \(\beta\) by the larger of \(\beta\) and \(s_w+h_S(d_w)\), and push its children, each only if its own bound still exceeds the updated \(\beta+\varepsilon\).

In exact arithmetic the result is a lower bound short of the truth by at most \(\varepsilon\). On one side \(\beta\) is always the exact support value of some square of the tree, so \(\beta\le h(d)\). On the other, every word is expanded or lies below a node still in the heap or below a discarded one; a discarded node had bound at most \(\beta+\varepsilon\) then, \(\beta\) only grows, and at the end every node left in the heap has bound at most \(\beta+\varepsilon\) too. Hence \(h(d)\le\beta+\varepsilon\).

It also terminates. Expanding a node raises \(\beta\) to at least \(s_w\), because \(h_S\ge0\), so its child \(wj\) is pushed only if

$$s_w+\langle d_w,t_j\rangle+5|c_j|\,|d_w|\gt s_w+\varepsilon,$$

and as \(\langle d_w,t_j\rangle\le|t_j|\,|d_w|\) and \(|t_j|+5|c_j|\le5\) by Section 4, this forces \(5|d_w|\gt\varepsilon\). With \(|d_w|\le|d|\,(4/5)^n\) at depth \(n\), no node deeper than \(\lceil\log(5|d|/\varepsilon)/\log(5/4)\rceil\) exists — \(162\) for \(|d|=1\) and \(\varepsilon=10^{-15}\) — so only finitely many words arise.

11. Where the extremes sit

All three options in Section 7 are positively homogeneous in \(d\), so which branch wins depends only on the direction of \(d\), not its length. Since \(b=-\frac34\,i\,a\) conjugates to \(\overline b=\frac34\,i\,\overline a\), the two steps turn a direction by the same angle up to a quarter turn, and the four cardinal directions differ by quarter turns too, so the four searches can synchronise — and at \(120\) digits with \(\varepsilon=10^{-90}\) they do: the optimal words for \(+x\), \(-x\), \(+y\) and \(-y\) begin \(2212\), \(1112\), \(1212\) and \(1111\), after which all four stand at the unit direction \(\left(-\tfrac{336}{625},\tfrac{527}{625}\right)\) and then agree throughout.

That common itinerary looks periodic but is not. Writing \(\alpha=\arg a\) gives \(\cos\alpha=4/5\), and by Niven's theorem a cosine at a rational multiple of \(\pi\) can only be \(0\), \(\pm\tfrac12\) or \(\pm1\), ^[6] so \(\alpha/\pi\) is irrational; after \(n\) steps a direction that started at angle \(\theta\) stands at \(\theta-n\alpha+m\cdot90^\circ\) for an integer \(m\), and no direction is ever revisited. Yet \(83\alpha\) misses a multiple of \(90^\circ\) by only \(0.2^\circ\), the closest such near-return under \(200\) steps, and the tail does repeat with period \(83\) for its first \(470\) symbols before breaking. A long near-period is not a period: the extremes are not fixed points of one similarity, and evaluating the recursion numerically is the right tool. To four decimals, \(x_{\max}\approx3.1307\) at \((3.1307,1.5455)\), \(x_{\min}\approx-3.2353\) at \((-3.2353,1.9247)\), \(y_{\max}\approx4.3209\) at \((-0.0873,4.3209)\) and \(y_{\min}\approx-0.1161\) at \((-1.6610,-0.1161)\): the tree dips below the base line, so the rectangle is not bounded by \(y=0\).

How the Code Works

All three programs run the same search. ^[10] Each has Vec2, dot and norm, the transpose helper apply_transpose_mul(d, c) returning \(\overline c\,d\) as derived in Section 7 — the C++ file carries that derivation as a comment — and support_unit_square for \(h_S(d)=\max(0,d_x,d_y,d_x+d_y)\); Java spells the last two applyTransposeMul and supportUnitSquare. Only C++ also has i_times (the quarter turn \(iz\)) and complex_mul, needed by its checkpoint generator.

SupportSolver stores the constants of Section 2 as a_, b_, t1_ and t2_, and its support(direction, tol) is the search of Section 10 with kRadialBound = 5 from Section 4. A Node holds one upper bound and the pair \((s_w,d_w)\) of Section 9, and NodeCmp compares the bounds so that std::priority_queue is a max-heap on them. ^[7] The loop condition pq.top().upper_bound > best + tol is the stopping rule itself; inside, best rises to cur.offset + support_unit_square(cur.dir) and each child is built by the update formulas of Section 9, as off1 = cur.offset + dot(cur.dir, t1_) and dir1 = apply_transpose_mul(cur.dir, a_), then pushed only if its bound beats best + tol.

finite_depth_bbox(depth) is the naive generator of Section 5, kept only for the checkpoint: a square is an AffineSquare{p, z}, a base point with a base vector, whose corners are p, p+z, p+iz, p+z+iz, and whose children follow \(f_1\) and \(f_2\). run_checkpoints compares its depth-\(1\) output with \(\left[-\tfrac{12}{25},\tfrac{37}{25}\right]\times\left[0,\tfrac{53}{25}\right]\), then checks the recursion of Section 7 in the four cardinal directions — both through nearly_equal, whose tolerance is relative: the difference must not exceed the tolerance, \(10^{-12}\) and \(2\cdot10^{-13}\) respectively, times the larger value. main exits with an error on an unknown argument or a failed checkpoint; otherwise solve_area evaluates the four support values, and their product is cast to double and printed with setprecision(10).

The ports keep that structure in solve and differ only in how each language expresses a priority queue and in its arithmetic. Python's heapq is a min-heap, so the key is the negated bound and a node is the flat tuple (-ub, offset, dir.x, dir.y) of plain floats, which ties compare without error. Java's PriorityQueue keeps the least element at the head, so Node.compareTo reverses the test as Double.compare(o.upperBound, this.upperBound), and Locale.US pins the decimal point. ^[8] C++ accumulates in long double, whose width is platform-dependent, ^[9] while the ports use double precision. One behavioural difference: C++ leaves the loop as soon as the top of the heap cannot beat best + tol, whereas the ports test the popped node and continue, discarding the few hundred nodes still queued. The result is identical, and only C++ carries the checkpoints.

Complexity and Verification

A query that expands \(N\) nodes pushes at most two children each time, so the heap costs \(O(N\log N)\) time and \(O(N)\) memory; the area needs four queries and the C++ checkpoints add twelve more, three per cardinal direction. Section 10 bounds the depth by \(O(\log(1/\varepsilon))\), but a binary tree of that depth could still hold exponentially many nodes: the shrinking gap of Section 9 only guarantees that a branch is abandoned once its square falls more than \(5|d_w|\) short of \(\beta\), so the surviving set being a thin tube around the optimal word is an observation, not a consequence of the bound. At \(\varepsilon=10^{-15}\) we measured \(396\), \(390\), \(386\) and \(385\) expansions for \(+x\), \(-x\), \(+y\) and \(-y\) — \(1557\) in all — with at most about \(294\) nodes in the heap and a deepest word of \(107\) symbols. A few hundred records of four numbers is negligible memory, and on one ordinary machine every program finished in milliseconds, most of it process start-up.

Two errors must stay below the last printed decimal. The truncation \(\varepsilon\) leaves each side short by at most \(2\varepsilon\) in exact arithmetic, so the area is short by at most \(2\varepsilon(W+H)+4\varepsilon^2\), with \(2(W+H)\approx21.6\); measured against the \(\varepsilon=10^{-15}\) run, \(\varepsilon=10^{-12}\) moves the area by \(1.7\cdot10^{-11}\) and \(\varepsilon=10^{-14}\) by \(1.4\cdot10^{-13}\). Floating-point rounding, which \(\varepsilon\) does not control, is the second: against a \(50\)-digit reference the double-precision support values deviate by up to \(5.1\cdot10^{-15}\), in two of the four directions even exceeding the true support, so the exact lower-bound property does not survive rounding, and the area deviates by about \(1.7\cdot10^{-14}\). Ten decimals need a total error below \(5\cdot10^{-11}\), and the exact area lies about \(3.0\cdot10^{-11}\) from the nearest rounding boundary, so both sources together are three orders of magnitude too small to change the printed value — whereas \(\varepsilon=10^{-12}\) alone would be of the same order as that margin. The \(50\)-digit reference rounds to the same ten decimals that all three programs print.

What was verified for this article, all by separate code:

• The geometry in exact rationals — the apex \((16/25,37/25)\), the leg lengths \(4/5\) and \(3/5\), \(a\cdot b=0\), \(b=-\tfrac34\,i\,a\) — and, by an independent enumeration of all square corners, the depth-\(2\) box \(\left[-\tfrac{684}{625},\tfrac{1204}{625}\right]\times\left[0,\tfrac{1697}{625}\right]\), which the depth-limited recursion of Section 8 reproduces exactly, as it does the depth-\(1\), \(3\) and \(4\) boxes.
• Both C++ checkpoints, recomputed outside the program: residuals below \(2\cdot10^{-15}\) for the recursion of Section 7.
• A second, provably valid upper bound in place of the disc of radius \(5\): the exact depth-\(10\) box widened on every side by \(5\,(4/5)^{10}\approx0.5369\), which by Section 5 contains \(K\), its support taken over the four corners. It returns the same four values — three bit for bit, the fourth within \(3\cdot10^{-16}\) — in \(2000\) expansions instead of \(1557\).
• An enumeration of square corners through depth \(22\), which gives lower bounds, together with \(s_w+5|d_w|\) over its \(2^{22}\) last-level words, which bounds the unexplored tails. With no pruning anywhere they bracket every support value, for example \(3.1292829\le h(1,0)\le3.1310500\).
• The radius of Section 4: a search for the farthest point of the tree gives \(\rho\approx4.3290161979\), below \(5\).

Footnotes and References

The references below supply the background facts used above; the coordinates, the fixed-point equation, the radius bound \(5\), the support recursion, the correctness of the search and the synchronisation of the optimal words are derived here.

1. Project Euler 395 — Pythagorean Tree. The official statement attaches a right triangle with sides in a 3-4-5 ratio to the side of a unit square opposite its base, the smaller side on the right, puts a square on each leg and repeats for ever; it asks for the smallest area of an enclosing rectangle with sides parallel to the largest square, to 10 decimal places.
2. Wikipedia — Pythagoras tree (fractal). It describes the classical right-isosceles version, in which every iteration adds squares of total area equal to the first square's, the whole tree fits in a box of 6L by 4L for a largest square of side L, and squares begin to overlap at the fifth iteration.
3. Wikipedia — Iterated function system. A finite set of contractions on a complete metric space has, as Hutchinson proved in 1981, a unique non-empty compact fixed set under the operator sending a set to the closure of the union of its images; Section 2 argues the same way for the operator that also keeps the first square.
4. Wikipedia — Support function. It defines the support function as the supremum of the inner product with a direction and states that it is positively homogeneous, that a translation adds the inner product with the shift vector, and that a closed convex set is determined by it; the union and linear-image rules are proved in Section 6.
5. Wikipedia — Branch and bound. The method bounds the best value attainable in each part of a split search space and discards a part whose bound cannot contain an optimum; a priority queue ordered by that bound gives the best-first variant of Section 10.
6. Wikipedia — Niven's theorem. At a rational multiple of pi, a sine or cosine can be rational only if its value is 0, plus or minus one half, or plus or minus one — the fact Section 11 applies to an angle of cosine 4/5.
7. cppreference — std::priority_queue. The adaptor gives constant-time access to the largest element by default, with logarithmic insertion and extraction, and accepts a comparator that changes the ordering; NodeCmp is that comparator here.
8. Python — heapq documents a min-heap whose smallest element is always heap[0], and recommends a tuple led by the priority, hence the negated bound. Java — PriorityQueue documents that the head is the least element under the natural ordering or a supplied comparator, hence the reversed compareTo.
9. cppreference — fundamental types. long double need not match any IEEE-754 format: many x86 implementations use the 80-bit x87 format, some platforms binary128, and some make it identical to double, so its width is not something a solution may rely on.
10. C++ — Euler395.cpp. The linked immutable revision contains SupportSolver, apply_transpose_mul, finite_depth_bbox, kRadialBound and the run_checkpoints harness that --skip-checkpoints disables. It is the implementation record for this article; the cross-checks above are additional validation, not tests packaged in that file.