Detailed mathematical approach

Problem Summary

Problem 395 asks for the area of the axis-aligned bounding box that contains the entire infinite Pythagorean tree. The local C++, Python, and Java solutions do not enumerate squares to a fixed depth. Instead they model the tree as a self-similar set \(K \subset \mathbb{R}^2\), answer four support-function extremal queries, and finally obtain the numerical area \(28.2453753155\).

Mathematical Approach

Self-Similar Geometry

In the C++ checkpoint code, a square is represented by a start point \(p\) and a base vector \(z\). Its four corners are

$p,\qquad p+z,\qquad p+i z,\qquad p+z+i z,$

where \(i(x,y)=(-y,x)\) is a quarter-turn. The two child squares are obtained by multiplying the base vector by the complex-like constants

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

These have lengths \(|a|=\frac{4}{5}\) and \(|b|=\frac{3}{5}\), so each recursive step both rotates and shrinks. For the root square \(S=[0,1]^2\), the child translations used by the solver are

$t_1=(0,1),\qquad t_2=t_1+a=\left(\frac{16}{25},\frac{37}{25}\right).$

If \(M_c\) denotes the linear map corresponding to complex multiplication by \(c=(u,v)\), then

$M_c(x,y)=(ux-vy,\ vx+uy).$

With \(A=M_a\) and \(B=M_b\), the entire infinite tree is the unique compact set satisfying the affine fixed-point equation

$K=S\cup (t_1+A K)\cup (t_2+B K).$

Support-Function Fixed Point

For any direction \(d\in\mathbb{R}^2\), define the support function

$h(d)=\sup_{x\in K} d\cdot x.$

Two standard identities turn the set equation into a scalar recurrence:

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

$h_{t+M X}(d)=d\cdot t+h_X(M^{\mathsf{T}}d).$

Therefore the exact support value satisfies

$h(d)=\max\left(h_S(d),\ d\cdot t_1+h(A^{\mathsf{T}}d),\ d\cdot t_2+h(B^{\mathsf{T}}d)\right).$

The helper support_unit_square is exact because the unit square has vertices \((0,0)\), \((1,0)\), \((0,1)\), and \((1,1)\), so

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

The transpose action implemented by apply_transpose_mul is

$M_c^{\mathsf{T}}(d_x,d_y)=(u d_x+v d_y,\ -v d_x+u d_y).$

Global Radius Bound and Branch-and-Bound

The heap search needs an upper bound for any unfinished subtree. The solver uses the global radial bound \(R=5\), and this follows directly from the same self-similar equation. If every point of \(K\) satisfies \(\|x\|\le R\), then the root square contributes at most \(\sqrt{2}\), the left subtree contributes at most \(1+\frac{4R}{5}\), and the right subtree contributes at most \(\frac{\sqrt{65}}{5}+\frac{3R}{5}\). So it is enough that

$R\ge \max\left(\sqrt{2},\ 1+\frac{4R}{5},\ \frac{\sqrt{65}}{5}+\frac{3R}{5}\right).$

This reduces to \(R\ge 5\) and \(R\ge \frac{\sqrt{65}}{2}\), hence \(R=5\) is safe. For any affine copy \(u+M K\), we then obtain the upper bound

$\sup_{x\in u+M K} d\cdot x \le d\cdot u+\|M^{\mathsf{T}}d\|R.$

After a finite child word \(w\), the program stores precisely the pair \(\text{offset}=d\cdot u_w\) and \(\text{dir}=M_w^{\mathsf{T}}d\). The exact contribution is \(\text{offset}+h(\text{dir})\), the lower bound is \(\text{offset}+h_S(\text{dir})\), and the upper bound is \(\text{offset}+\|\text{dir}\|R\). The priority queue always expands the node with the largest upper bound first. Nodes with upper bound at most best + tol are discarded, and because each child multiplies the direction norm by \(\frac{4}{5}\) or \(\frac{3}{5}\), promising bounds decay geometrically.

Bounding Box Extraction

Once support values are available, the axis-aligned box comes from the four coordinate directions:

$x_{\max}=h(1,0),\qquad x_{\min}=-h(-1,0),$

$y_{\max}=h(0,1),\qquad y_{\min}=-h(0,-1).$

Evaluating the local solver gives

$x_{\min}\approx -3.235295198730329,\qquad x_{\max}\approx 3.130653549285805,$

$y_{\min}\approx -0.116075304973770,\qquad y_{\max}\approx 4.320871399047746,$

so the final area is

$A=(x_{\max}-x_{\min})(y_{\max}-y_{\min})\approx 28.245375315480086.$

Checkpoint Used by the C++ Version

Before solving the infinite problem, the C++ file verifies the first generation exactly. The left child has corners \((0,1)\), \((\frac{16}{25},\frac{37}{25})\), \((-\frac{12}{25},\frac{41}{25})\), and \((\frac{4}{25},\frac{53}{25})\). The right child has corners \((\frac{16}{25},\frac{37}{25})\), \((1,1)\), \((\frac{28}{25},\frac{46}{25})\), and \((\frac{37}{25},\frac{34}{25})\). Together with the root square, this gives the exact depth-1 box

$\left[-\frac{12}{25},\frac{37}{25}\right]\times\left[0,\frac{53}{25}\right],$

which is exactly what finite_depth_bbox(1) checks.

How the Code Works

SupportSolver stores the constants \(a\), \(b\), \(t_1\), and \(t_2\). The function apply_transpose_mul computes \(M_c^{\mathsf{T}}d\), support_unit_square evaluates the base square exactly, and the heap stores triples \((\text{upper\_bound},\text{offset},\text{dir})\). The C++ version uses long double, checks the depth-1 geometry, and also verifies the support fixed-point identity in the four cardinal directions. The Python and Java versions implement the same best-first search with the same constants, but without the extra checkpoint harness.

Complexity Analysis

If one support query expands \(N_d\) nodes, the heap operations cost \(O(N_d\log N_d)\) time and \(O(N_d)\) memory. The final area requires four such queries, so the overall cost is the same order up to a factor of four. This is very different from naive depth-\(n\) generation, which touches \(2^n\) squares and still only approximates the infinite tree. Here the search is adaptive: entire subtrees are pruned as soon as their best possible upper bound cannot beat the current optimum.

Further Reading

1. Problem page: https://projecteuler.net/problem=395
2. Pythagorean tree: Wikipedia — Pythagoras tree
3. Iterated function system: Wikipedia — Iterated function system
4. Support function: Wikipedia — Support function
5. Branch and bound: Wikipedia — Branch and bound