Problem 395: Pythagorean Tree

Project Euler Problem 395 Solution

EulerSolve provides an optimized solution for Project Euler Problem 395, Pythagorean Tree, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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.

Problem 395 source code

C++

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <queue>
#include <string>
#include <tuple>
#include <vector>

namespace {

struct Options {
    bool run_checkpoints = true;
};

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

struct Vec2 {
    long double x{0.0L};
    long double y{0.0L};
};

Vec2 operator+(const Vec2& a, const Vec2& b) {
    return Vec2{a.x + b.x, a.y + b.y};
}

Vec2 operator-(const Vec2& a, const Vec2& b) {
    return Vec2{a.x - b.x, a.y - b.y};
}

Vec2 operator*(const Vec2& a, const long double k) {
    return Vec2{a.x * k, a.y * k};
}

long double dot(const Vec2& a, const Vec2& b) {
    return a.x * b.x + a.y * b.y;
}

long double norm(const Vec2& a) {
    return std::sqrt(a.x * a.x + a.y * a.y);
}

Vec2 i_times(const Vec2& z) {
    return Vec2{-z.y, z.x};
}

Vec2 complex_mul(const Vec2& z, const Vec2& c) {
    return Vec2{z.x * c.x - z.y * c.y, z.x * c.y + z.y * c.x};
}

// For c = u + iv, A_c = [[u,-v],[v,u]], so A_c^T d = [u dx + v dy, -v dx + u dy].
Vec2 apply_transpose_mul(const Vec2& d, const Vec2& c) {
    return Vec2{c.x * d.x + c.y * d.y, -c.y * d.x + c.x * d.y};
}

long double support_unit_square(const Vec2& d) {
    return std::max({0.0L, d.x, d.y, d.x + d.y});
}

struct Node {
    long double upper_bound;
    long double offset;
    Vec2 dir;
};

struct NodeCmp {
    bool operator()(const Node& a, const Node& b) const {
        return a.upper_bound < b.upper_bound;
    }
};

class SupportSolver {
   public:
    SupportSolver()
        : a_{16.0L / 25.0L, 12.0L / 25.0L},
          b_{9.0L / 25.0L, -12.0L / 25.0L},
          t1_{0.0L, 1.0L},
          t2_{16.0L / 25.0L, 37.0L / 25.0L} {}

    long double support(const Vec2& direction, const long double tol = 1e-15L) const {
        constexpr long double kRadialBound = 5.0L;  // Global |z| upper bound for the full tree.

        long double best = support_unit_square(direction);
        std::priority_queue<Node, std::vector<Node>, NodeCmp> pq;
        pq.push(Node{norm(direction) * kRadialBound, 0.0L, direction});

        while (!pq.empty() && pq.top().upper_bound > best + tol) {
            const Node cur = pq.top();
            pq.pop();

            best = std::max(best, cur.offset + support_unit_square(cur.dir));

            const Vec2 dir1 = apply_transpose_mul(cur.dir, a_);
            const long double off1 = cur.offset + dot(cur.dir, t1_);
            const long double ub1 = off1 + norm(dir1) * kRadialBound;
            if (ub1 > best + tol) {
                pq.push(Node{ub1, off1, dir1});
            }

            const Vec2 dir2 = apply_transpose_mul(cur.dir, b_);
            const long double off2 = cur.offset + dot(cur.dir, t2_);
            const long double ub2 = off2 + norm(dir2) * kRadialBound;
            if (ub2 > best + tol) {
                pq.push(Node{ub2, off2, dir2});
            }
        }

        return best;
    }

    Vec2 a() const { return a_; }
    Vec2 b() const { return b_; }
    Vec2 t1() const { return t1_; }
    Vec2 t2() const { return t2_; }

   private:
    Vec2 a_;
    Vec2 b_;
    Vec2 t1_;
    Vec2 t2_;
};

struct AffineSquare {
    Vec2 p;  // base start
    Vec2 z;  // base vector
};

std::tuple<long double, long double, long double, long double> finite_depth_bbox(const int depth) {
    const Vec2 a{16.0L / 25.0L, 12.0L / 25.0L};
    const Vec2 b{9.0L / 25.0L, -12.0L / 25.0L};
    std::vector<AffineSquare> current;
    current.push_back(AffineSquare{Vec2{0.0L, 0.0L}, Vec2{1.0L, 0.0L}});

    long double xmin = std::numeric_limits<long double>::infinity();
    long double xmax = -std::numeric_limits<long double>::infinity();
    long double ymin = std::numeric_limits<long double>::infinity();
    long double ymax = -std::numeric_limits<long double>::infinity();

    for (int level = 0; level <= depth; ++level) {
        for (const AffineSquare& sq : current) {
            const Vec2 iz = i_times(sq.z);
            const Vec2 c0 = sq.p;
            const Vec2 c1 = sq.p + sq.z;
            const Vec2 c2 = sq.p + iz;
            const Vec2 c3 = sq.p + sq.z + iz;

            for (const Vec2& c : {c0, c1, c2, c3}) {
                xmin = std::min(xmin, c.x);
                xmax = std::max(xmax, c.x);
                ymin = std::min(ymin, c.y);
                ymax = std::max(ymax, c.y);
            }
        }

        if (level == depth) {
            break;
        }

        std::vector<AffineSquare> next;
        next.reserve(current.size() * 2);
        for (const AffineSquare& sq : current) {
            const Vec2 iz = i_times(sq.z);
            const Vec2 left_p = sq.p + iz;
            const Vec2 left_z = complex_mul(sq.z, a);

            const Vec2 right_p = sq.p + iz + left_z;
            const Vec2 right_z = complex_mul(sq.z, b);

            next.push_back(AffineSquare{left_p, left_z});
            next.push_back(AffineSquare{right_p, right_z});
        }
        current.swap(next);
    }

    return {xmin, xmax, ymin, ymax};
}

bool nearly_equal(const long double a, const long double b, const long double tol = 1e-12L) {
    const long double scale = std::max(std::fabsl(a), std::fabsl(b));
    if (scale == 0.0L) {
        return true;
    }
    return std::fabsl(a - b) <= tol * scale;
}

bool run_checkpoints() {
    const auto [xmin1, xmax1, ymin1, ymax1] = finite_depth_bbox(1);
    if (!nearly_equal(xmin1, -12.0L / 25.0L) ||
        !nearly_equal(xmax1, 37.0L / 25.0L) ||
        !nearly_equal(ymin1, 0.0L) ||
        !nearly_equal(ymax1, 53.0L / 25.0L)) {
        std::cerr << "Checkpoint failed: finite depth-1 geometry\n";
        return false;
    }

    const SupportSolver solver;
    const std::vector<Vec2> dirs{{1.0L, 0.0L}, {-1.0L, 0.0L}, {0.0L, 1.0L}, {0.0L, -1.0L}};
    for (const Vec2& d : dirs) {
        const long double h = solver.support(d);
        const long double rhs = std::max({
            support_unit_square(d),
            dot(d, solver.t1()) + solver.support(apply_transpose_mul(d, solver.a())),
            dot(d, solver.t2()) + solver.support(apply_transpose_mul(d, solver.b())),
        });
        if (!nearly_equal(h, rhs, 2e-13L)) {
            std::cerr << "Checkpoint failed: support fixed-point residual\n";
            return false;
        }
    }

    return true;
}

long double solve_area() {
    const SupportSolver solver;
    const long double xmax = solver.support(Vec2{1.0L, 0.0L});
    const long double xmin = -solver.support(Vec2{-1.0L, 0.0L});
    const long double ymax = solver.support(Vec2{0.0L, 1.0L});
    const long double ymin = -solver.support(Vec2{0.0L, -1.0L});
    return (xmax - xmin) * (ymax - ymin);
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << std::fixed << std::setprecision(10) << static_cast<double>(solve_area()) << '\n';
    return 0;
}

Python

import heapq
import math

class Vec2:
    def __init__(self, x, y):
        self.x = float(x)
        self.y = float(y)

def norm(v):
    return math.sqrt(v.x * v.x + v.y * v.y)

def apply_transpose_mul(d, c):
    return Vec2(c.x * d.x + c.y * d.y, -c.y * d.x + c.x * d.y)

def dot(a, b):
    return a.x * b.x + a.y * b.y

def support_unit_square(d):
    return max(0.0, d.x, d.y, d.x + d.y)

class SupportSolver:
    def __init__(self):
        self.a = Vec2(16.0 / 25.0, 12.0 / 25.0)
        self.b = Vec2(9.0 / 25.0, -12.0 / 25.0)
        self.t1 = Vec2(0.0, 1.0)
        self.t2 = Vec2(16.0 / 25.0, 37.0 / 25.0)

    def support(self, direction, tol=1e-15):
        kRadialBound = 5.0
        best = support_unit_square(direction)
        
        pq = []
        heapq.heappush(pq, (-norm(direction) * kRadialBound, 0.0, direction.x, direction.y))
        
        while pq:
            neg_ub, offset, dx, dy = heapq.heappop(pq)
            ub = -neg_ub
            
            if ub <= best + tol:
                continue
                
            cur_dir = Vec2(dx, dy)
            best = max(best, offset + support_unit_square(cur_dir))
            
            dir1 = apply_transpose_mul(cur_dir, self.a)
            off1 = offset + dot(cur_dir, self.t1)
            ub1 = off1 + norm(dir1) * kRadialBound
            if ub1 > best + tol:
                heapq.heappush(pq, (-ub1, off1, dir1.x, dir1.y))
                
            dir2 = apply_transpose_mul(cur_dir, self.b)
            off2 = offset + dot(cur_dir, self.t2)
            ub2 = off2 + norm(dir2) * kRadialBound
            if ub2 > best + tol:
                heapq.heappush(pq, (-ub2, off2, dir2.x, dir2.y))
                
        return best

def solve():
    solver = SupportSolver()
    xmax = solver.support(Vec2(1.0, 0.0))
    xmin = -solver.support(Vec2(-1.0, 0.0))
    ymax = solver.support(Vec2(0.0, 1.0))
    ymin = -solver.support(Vec2(0.0, -1.0))
    
    ans = (xmax - xmin) * (ymax - ymin)
    return "{:.10f}".format(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.PriorityQueue;

public class Euler395 {
    static class Vec2 {
        double x, y;

        Vec2(double x, double y) {
            this.x = x;
            this.y = y;
        }
    }

    static double norm(Vec2 v) {
        return Math.sqrt(v.x * v.x + v.y * v.y);
    }

    static Vec2 applyTransposeMul(Vec2 d, Vec2 c) {
        return new Vec2(c.x * d.x + c.y * d.y, -c.y * d.x + c.x * d.y);
    }

    static double dot(Vec2 a, Vec2 b) {
        return a.x * b.x + a.y * b.y;
    }

    static double supportUnitSquare(Vec2 d) {
        return Math.max(0.0, Math.max(d.x, Math.max(d.y, d.x + d.y)));
    }

    static class Node implements Comparable<Node> {
        double upperBound;
        double offset;
        Vec2 dir;

        Node(double upperBound, double offset, Vec2 dir) {
            this.upperBound = upperBound;
            this.offset = offset;
            this.dir = dir;
        }

        @Override
        public int compareTo(Node o) {
            return Double.compare(o.upperBound, this.upperBound);
        }
    }

    static class SupportSolver {
        Vec2 a = new Vec2(16.0 / 25.0, 12.0 / 25.0);
        Vec2 b = new Vec2(9.0 / 25.0, -12.0 / 25.0);
        Vec2 t1 = new Vec2(0.0, 1.0);
        Vec2 t2 = new Vec2(16.0 / 25.0, 37.0 / 25.0);

        double support(Vec2 direction, double tol) {
            double kRadialBound = 5.0;
            double best = supportUnitSquare(direction);

            PriorityQueue<Node> pq = new PriorityQueue<>();
            pq.offer(new Node(norm(direction) * kRadialBound, 0.0, direction));

            while (!pq.isEmpty()) {
                Node cur = pq.poll();
                if (cur.upperBound <= best + tol)
                    continue;

                best = Math.max(best, cur.offset + supportUnitSquare(cur.dir));

                Vec2 dir1 = applyTransposeMul(cur.dir, a);
                double off1 = cur.offset + dot(cur.dir, t1);
                double ub1 = off1 + norm(dir1) * kRadialBound;
                if (ub1 > best + tol) {
                    pq.offer(new Node(ub1, off1, dir1));
                }

                Vec2 dir2 = applyTransposeMul(cur.dir, b);
                double off2 = cur.offset + dot(cur.dir, t2);
                double ub2 = off2 + norm(dir2) * kRadialBound;
                if (ub2 > best + tol) {
                    pq.offer(new Node(ub2, off2, dir2));
                }
            }
            return best;
        }
    }

    static String solve() {
        SupportSolver solver = new SupportSolver();
        double xmax = solver.support(new Vec2(1.0, 0.0), 1e-15);
        double xmin = -solver.support(new Vec2(-1.0, 0.0), 1e-15);
        double ymax = solver.support(new Vec2(0.0, 1.0), 1e-15);
        double ymin = -solver.support(new Vec2(0.0, -1.0), 1e-15);

        double area = (xmax - xmin) * (ymax - ymin);
        return String.format(java.util.Locale.US, "%.10f", area);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}