Problem 904: Pythagorean Angle

Project Euler Problem 904 Solution

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

Problem Summary The problem asks for $$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$ with $N=45000$ and $L=10^{10}$. For one target angle $\alpha$ measured in degrees, we consider integer right triangles coming from Pythagorean triples. Each admissible triangle determines a special angle $\theta$, and $f(\alpha,L)$ is the perimeter of the scaled triangle whose angle is closest to $\alpha$, subject to the hypotenuse bound $c\le L$ after scaling. The implementation does not brute-force all triples. Instead it converts the geometry into a one-variable function of the ratio $u=m/n$, inverts that function numerically, and then searches only for nearby rational values that can really occur in primitive Pythagorean triples. Mathematical Approach The entire solution rests on two facts. First, every primitive integer right triangle has Euclid parameters $m>n>0$. Second, the relevant angle depends only on the quotient $u=m/n$, so the continuous optimization is one-dimensional before it is converted back to integer data. Primitive triples, scaling, and the quantity to maximize Every primitive Pythagorean triple can be written as $$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$ with $$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$ Any non-primitive solution is just a scaled copy $(ka,kb,kc)$ of such a primitive triangle....

Detailed mathematical approach

Problem Summary

The problem asks for

$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$

with $N=45000$ and $L=10^{10}$. For one target angle $\alpha$ measured in degrees, we consider integer right triangles coming from Pythagorean triples. Each admissible triangle determines a special angle $\theta$, and $f(\alpha,L)$ is the perimeter of the scaled triangle whose angle is closest to $\alpha$, subject to the hypotenuse bound $c\le L$ after scaling.

The implementation does not brute-force all triples. Instead it converts the geometry into a one-variable function of the ratio $u=m/n$, inverts that function numerically, and then searches only for nearby rational values that can really occur in primitive Pythagorean triples.

Mathematical Approach

The entire solution rests on two facts. First, every primitive integer right triangle has Euclid parameters $m>n>0$. Second, the relevant angle depends only on the quotient $u=m/n$, so the continuous optimization is one-dimensional before it is converted back to integer data.

Primitive triples, scaling, and the quantity to maximize

Every primitive Pythagorean triple can be written as

$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$

with

$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$

Any non-primitive solution is just a scaled copy $(ka,kb,kc)$ of such a primitive triangle. Scaling does not change the associated angle, so once a primitive triple is chosen, the best allowed scale is simply

$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor.$$

The value returned by $f(\alpha,L)$ is therefore

$$k(a+b+c).$$

So the mathematical task is: among primitive pairs $(m,n)$, find the one whose angle is closest to $\alpha$, then multiply its perimeter by the largest admissible scale.

The angle as a function of a single ratio

The problem-specific angle is encoded through

$$\tan\theta=\frac{3ab}{2c^2}.$$

Substituting Euclid's formulas and writing

$$u=\frac{m}{n}>1$$

gives

$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}.$$

This is the key collapse: the angle no longer depends on the absolute size of $m$ and $n$, only on their ratio.

For a target angle $\alpha$, the code computes

$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right).$$

Using the tangent-difference identity, the angular error can be measured by

$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|.$$

Because $\tan x$ is strictly increasing on the small angle range relevant here, minimizing the true angle difference is equivalent to minimizing $\Delta$.

The shape of $t(u)$ and why two branches appear

Differentiating yields

$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}.$$

For $u>1$, there is a single critical point, obtained from $u^4-6u^2+1=0$:

$$u_{\max}=1+\sqrt{2}.$$

At that point,

$$t(u_{\max})=\frac34,$$

so $t(u)$ increases on $(1,1+\sqrt2)$ and decreases on $(1+\sqrt2,\infty)$. Therefore a target value $t_\alpha<3/4$ has two inverse images:

$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$

Those are the two continuous ratios around which the search must be organized.

There is also a useful problem-specific observation: in the actual outer sum, $\alpha=\sqrt[3]{i}$ with $1\le i\le45000$, so

$$\alpha\le\sqrt[3]{45000}\approx 35.57^\circ\lt\arctan\!\left(\frac34\right)\approx 36.87^\circ.$$

Hence every angle used in the final computation lies below the peak. The production run always has two branches, even though the implementations still handle the peak case for completeness.

From continuous roots to admissible fractions

A primitive triple corresponds to a reduced fraction

$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$

with size constraint

$$p^2+q^2\le L.$$

If $p/q$ is close to one of the continuous roots $u_\ast$, then $p\approx u_\ast q$, so the bound above suggests the natural denominator scale

$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$

The implementations turn this into the Farey order

$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor.$$

They first find the Farey neighbors of order $N_u$ that bracket $u_\ast$. Then they walk outward in the Farey sequence until they encounter the nearest fractions on the lower and upper sides that satisfy all primitive-triple conditions. Those nearby admissible fractions are the concrete integer candidates that get evaluated.

This is why the solution is fast: it never scans the whole disk $m^2+n^2\le L$. It solves the continuous inverse problem first and only then inspects a tiny local set of rational approximations around the relevant roots.

Worked example: $\alpha=30^\circ$ and $L=100$

For the small checkpoint angle,

$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0.57735.$$

One excellent admissible fraction is

$$\frac{m}{n}=\frac92.$$

It generates

$$a=9^2-2^2=77,\qquad b=2\cdot 9\cdot 2=36,\qquad c=9^2+2^2=85.$$

Since $c=85\le100$, the largest allowed scale is $k=\lfloor 100/85\rfloor=1$. The associated tangent value is

$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0.57550,$$

so the error proxy is

$$\Delta=\frac{|0.57550-0.57735|}{1+0.57550\cdot0.57735}\approx 0.0013875.$$

This candidate beats the nearby admissible alternatives, so

$$f(30,100)=77+36+85=198.$$

The implementations use this value as a validation checkpoint in the C++ version, and the same mathematics appears in all three languages.

How the Code Works

Inverting the target angle

For each $\alpha$, the C++, Python, and Java implementations compute $t_\alpha$. When $t_\alpha<3/4$, they use bisection once on the increasing branch $(1,1+\sqrt2)$ and once on the decreasing branch $(1+\sqrt2,\infty)$, producing two continuous targets $u_-$ and $u_+$. If the target ever reached or exceeded the peak value $3/4$, the search would collapse to the neighborhood of $u=1+\sqrt2$.

Building and scoring integer candidates

For each continuous target $u_\ast$, the implementation computes the corresponding Farey order $N_u$, constructs the bracketing neighbors, and then walks outward to the closest admissible fractions on both sides. Every valid fraction is turned back into a primitive triple $(a,b,c)$, scaled by $k=\lfloor L/c\rfloor$, and scored using the same quantity $\Delta$.

If two candidates are numerically tied in angular error, the implementation prefers the one with larger scaled area. Since the scaled area is

$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2},$$

it is enough to compare $k^2ab$, which avoids any division.

Summing over all cube-root angles

After the single-angle routine is in place, the final answer is obtained by evaluating it at

$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$

and summing the resulting perimeters. The C++ and Java implementations split this outer loop across available processor threads and add partial sums at the end. The Python implementation performs the same computation serially.

Complexity Analysis

For one angle, the inverse stage uses a fixed number of bisection iterations, and the local rational search advances through only nearby Farey neighbors. In the literal implementations, even the outward walk is guarded by a fixed iteration cap, so the coded per-angle cost is bounded by a constant. For the full problem this means the runtime scales essentially linearly with $N$.

Memory usage is $O(1)$ per worker for the single-angle search. The threaded implementations need $O(T)$ additional storage for $T$ partial sums, which is negligible compared with the arithmetic work.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=904
Pythagorean triples: Wikipedia - Pythagorean triple
Euclid's formula: MathWorld - Pythagorean Triple
Farey sequences: Wikipedia - Farey sequence
Continued fractions: Wikipedia - Continued fraction
Tangent addition and subtraction identities: Wikipedia - List of trigonometric identities

Problem 904 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

struct Fraction {
    long long p;
    long long q;
};

struct Candidate {
    bool valid = false;
    long double diff = 0.0L;
    unsigned __int128 area = 0;
    uint64_t sum = 0;
};

static long double tan_theta_from_u(long double u) {
    long double u2 = u * u;
    return 3.0L * u * (u2 - 1.0L) / ((u2 + 1.0L) * (u2 + 1.0L));
}

static Fraction farey_prev(const Fraction& lo, const Fraction& hi, long long N) {
    long long k = (N + hi.q) / lo.q;
    return {k * lo.p - hi.p, k * lo.q - hi.q};
}

static Fraction farey_next(const Fraction& lo, const Fraction& hi, long long N) {
    long long k = (N + lo.q) / hi.q;
    return {k * hi.p - lo.p, k * hi.q - lo.q};
}

static void farey_bounds(long double x, long long N, Fraction& lo, Fraction& hi) {
    std::vector<long long> a;
    long double y = x;
    for (int i = 0; i < 64; ++i) {
        long long ai = static_cast<long long>(std::floor(y));
        a.push_back(ai);
        long double frac = y - static_cast<long double>(ai);
        if (std::fabsl(frac) < 1e-18L) {
            break;
        }
        y = 1.0L / frac;
    }

    std::vector<long long> ps;
    std::vector<long long> qs;
    long long p_minus2 = 0;
    long long p_minus1 = 1;
    long long q_minus2 = 1;
    long long q_minus1 = 0;

    for (long long ai : a) {
        long long p = ai * p_minus1 + p_minus2;
        long long q = ai * q_minus1 + q_minus2;
        ps.push_back(p);
        qs.push_back(q);
        if (q > N) {
            break;
        }
        p_minus2 = p_minus1;
        p_minus1 = p;
        q_minus2 = q_minus1;
        q_minus1 = q;
    }

    int idx = 0;
    while (idx < static_cast<int>(qs.size()) && qs[idx] <= N) {
        ++idx;
    }

    if (idx == 0) {
        lo = {0, 1};
        hi = {1, 0};
        return;
    }

    if (idx == static_cast<int>(qs.size())) {
        lo = {ps.back(), qs.back()};
        hi = lo;
        return;
    }

    long long p_prev = ps[idx - 1];
    long long q_prev = qs[idx - 1];
    long long p_prevprev = (idx >= 2) ? ps[idx - 2] : 1;
    long long q_prevprev = (idx >= 2) ? qs[idx - 2] : 0;

    long long t = (N - q_prevprev) / q_prev;
    long long p_cand = p_prevprev + t * p_prev;
    long long q_cand = q_prevprev + t * q_prev;

    if (static_cast<long double>(p_prev) / static_cast<long double>(q_prev) < x) {
        lo = {p_prev, q_prev};
        hi = {p_cand, q_cand};
    } else {
        lo = {p_cand, q_cand};
        hi = {p_prev, q_prev};
    }
}

static bool is_valid_fraction(const Fraction& f, long long L) {
    if (f.q <= 0 || f.p <= 0) {
        return false;
    }
    if (f.p <= f.q) {
        return false;
    }
    if (((f.p + f.q) & 1LL) == 0) {
        return false;
    }
    if (std::gcd(f.p, f.q) != 1) {
        return false;
    }
    unsigned long long m = static_cast<unsigned long long>(f.p);
    unsigned long long n = static_cast<unsigned long long>(f.q);
    unsigned long long m2 = m * m;
    unsigned long long n2 = n * n;
    return m2 + n2 <= static_cast<unsigned long long>(L);
}

static Fraction find_valid_lower(Fraction lo, Fraction hi, long long N, long long L) {
    if (lo.p <= lo.q) {
        return {0, 0};
    }
    for (int steps = 0; steps < 100000; ++steps) {
        if (is_valid_fraction(lo, L)) {
            return lo;
        }
        Fraction prev = farey_prev(lo, hi, N);
        hi = lo;
        lo = prev;
        if (lo.q == 0) {
            break;
        }
    }
    return {0, 0};
}

static Fraction find_valid_upper(Fraction lo, Fraction hi, long long N, long long L) {
    for (int steps = 0; steps < 100000; ++steps) {
        if (is_valid_fraction(hi, L)) {
            return hi;
        }
        Fraction next = farey_next(lo, hi, N);
        lo = hi;
        hi = next;
        if (hi.q == 0) {
            break;
        }
    }
    return {0, 0};
}

static Candidate evaluate_fraction(const Fraction& f, long double t_alpha, long long L) {
    Candidate cand;
    if (f.p == 0 || f.q == 0) {
        return cand;
    }
    unsigned long long m = static_cast<unsigned long long>(f.p);
    unsigned long long n = static_cast<unsigned long long>(f.q);
    unsigned long long m2 = m * m;
    unsigned long long n2 = n * n;
    unsigned long long a = m2 - n2;
    unsigned long long b = 2ULL * m * n;
    unsigned long long c = m2 + n2;
    unsigned long long k = static_cast<unsigned long long>(L) / c;
    if (k == 0) {
        return cand;
    }

    long double t = 3.0L * static_cast<long double>(a) * static_cast<long double>(b)
        / (2.0L * static_cast<long double>(c) * static_cast<long double>(c));
    long double diff = std::fabsl(t - t_alpha) / (1.0L + t * t_alpha);

    cand.valid = true;
    cand.diff = diff;
    cand.sum = static_cast<uint64_t>(k) * static_cast<uint64_t>(a + b + c);
    cand.area = static_cast<unsigned __int128>(k) * static_cast<unsigned __int128>(k)
        * static_cast<unsigned __int128>(a) * static_cast<unsigned __int128>(b);
    return cand;
}

static Candidate best_candidate_for_u(long double u_target, long double t_alpha, long long L) {
    Candidate best;
    if (u_target <= 1.0L) {
        return best;
    }

    long double R = sqrtl(static_cast<long double>(L));
    long double denom = sqrtl(1.0L + u_target * u_target);
    long long N = static_cast<long long>(std::floor(R / denom));
    if (N < 1) {
        N = 1;
    }

    Fraction lo{0, 1};
    Fraction hi{1, 0};
    farey_bounds(u_target, N, lo, hi);

    Fraction lo_valid = find_valid_lower(lo, hi, N, L);
    Fraction hi_valid = find_valid_upper(lo, hi, N, L);

    if (lo_valid.q != 0) {
        best = evaluate_fraction(lo_valid, t_alpha, L);
    }
    if (hi_valid.q != 0) {
        Candidate cand = evaluate_fraction(hi_valid, t_alpha, L);
        if (!best.valid) {
            best = cand;
        } else if (cand.valid) {
            const long double eps = 1e-21L;
            if (cand.diff + eps < best.diff) {
                best = cand;
            } else if (std::fabsl(cand.diff - best.diff) <= eps) {
                if (cand.area > best.area) {
                    best = cand;
                }
            }
        }
    }

    return best;
}

static uint64_t compute_f(long double alpha_deg, long long L) {
    static const long double kPi = acosl(-1.0L);
    static const long double kUMin = 1.0L + sqrtl(2.0L);
    static const long double kTMax = 0.75L;

    long double alpha_rad = alpha_deg * kPi / 180.0L;
    long double t_alpha = tanl(alpha_rad);

    long double u_low = 1.0L;
    long double u_high = kUMin;
    if (t_alpha < kTMax) {
        long double low = 1.0L;
        long double high = kUMin;
        for (int iter = 0; iter < 100; ++iter) {
            long double mid = 0.5L * (low + high);
            if (tan_theta_from_u(mid) < t_alpha) {
                low = mid;
            } else {
                high = mid;
            }
        }
        u_low = 0.5L * (low + high);

        low = kUMin;
        high = std::max(kUMin + 1.0L, 3.0L / t_alpha + 2.0L);
        while (tan_theta_from_u(high) > t_alpha) {
            high *= 2.0L;
        }
        for (int iter = 0; iter < 100; ++iter) {
            long double mid = 0.5L * (low + high);
            if (tan_theta_from_u(mid) > t_alpha) {
                low = mid;
            } else {
                high = mid;
            }
        }
        u_high = 0.5L * (low + high);
    }

    Candidate best = best_candidate_for_u(u_low, t_alpha, L);
    Candidate alt = best_candidate_for_u(u_high, t_alpha, L);
    if (!best.valid || (alt.valid && (alt.diff < best.diff - 1e-21L
        || (std::fabsl(alt.diff - best.diff) <= 1e-21L && alt.area > best.area)))) {
        best = alt;
    }

    return best.sum;
}

static uint64_t compute_F(int N, long long L) {
    unsigned int threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 1;
    }

    std::vector<uint64_t> partial(threads, 0);
    std::vector<std::thread> pool;
    pool.reserve(threads);

    for (unsigned int t = 0; t < threads; ++t) {
        pool.emplace_back([t, threads, N, L, &partial]() {
            uint64_t local = 0;
            for (int i = 1 + static_cast<int>(t); i <= N; i += static_cast<int>(threads)) {
                long double alpha = cbrtl(static_cast<long double>(i));
                local += compute_f(alpha, L);
            }
            partial[t] = local;
        });
    }

    for (auto& th : pool) {
        th.join();
    }

    uint64_t total = 0;
    for (uint64_t v : partial) {
        total += v;
    }
    return total;
}

static bool run_validation() {
    bool ok = true;
    if (compute_f(30.0L, 100LL) != 198ULL) {
        std::cerr << "Validation failed: f(30, 10^2) != 198\n";
        ok = false;
    }
    if (compute_f(10.0L, 1000000LL) != 1600158ULL) {
        std::cerr << "Validation failed: f(10, 10^6) != 1600158\n";
        ok = false;
    }
    if (compute_F(10, 1000000LL) != 16684370ULL) {
        std::cerr << "Validation failed: F(10, 10^6) != 16684370\n";
        ok = false;
    }
    if (ok) {
        std::cerr << "Validation checkpoints passed.\n";
    }
    return ok;
}

int main() {
    if (!run_validation()) {
        return 1;
    }

    uint64_t result = compute_F(45000, 10000000000LL);
    std::cout << result << "\n";
    return 0;
}

Python

import math
from typing import List

class Fraction:
    def __init__(self, p, q):
        self.p = p
        self.q = q

def tan_theta_from_u(u):
    u2 = u * u
    return 3.0 * u * (u2 - 1.0) / ((u2 + 1.0) * (u2 + 1.0))

def farey_prev(lo, hi, N):
    k = (N + hi.q) // lo.q
    return Fraction(k * lo.p - hi.p, k * lo.q - hi.q)

def farey_next(lo, hi, N):
    k = (N + lo.q) // hi.q
    return Fraction(k * hi.p - lo.p, k * hi.q - lo.q)

def farey_bounds(x, N):
    a = []
    y = x
    for _ in range(64):
        ai = int(math.floor(y))
        a.append(ai)
        frac = y - ai
        if abs(frac) < 1e-15:
            break
        y = 1.0 / frac

    ps = []
    qs = []
    p_minus2 = 0
    p_minus1 = 1
    q_minus2 = 1
    q_minus1 = 0

    for ai in a:
        p = ai * p_minus1 + p_minus2
        q = ai * q_minus1 + q_minus2
        ps.append(p)
        qs.append(q)
        if q > N:
            break
        p_minus2 = p_minus1
        p_minus1 = p
        q_minus2 = q_minus1
        q_minus1 = q

    idx = 0
    while idx < len(qs) and qs[idx] <= N:
        idx += 1

    if idx == 0:
        return Fraction(0, 1), Fraction(1, 0)

    if idx == len(qs):
        return Fraction(ps[-1], qs[-1]), Fraction(ps[-1], qs[-1])

    p_prev = ps[idx - 1]
    q_prev = qs[idx - 1]
    p_prevprev = ps[idx - 2] if idx >= 2 else 1
    q_prevprev = qs[idx - 2] if idx >= 2 else 0

    t = (N - q_prevprev) // q_prev
    p_cand = p_prevprev + t * p_prev
    q_cand = q_prevprev + t * q_prev

    if p_prev / q_prev < x:
        return Fraction(p_prev, q_prev), Fraction(p_cand, q_cand)
    else:
        return Fraction(p_cand, q_cand), Fraction(p_prev, q_prev)

def is_valid_fraction(f, L):
    if f.q <= 0 or f.p <= 0:
        return False
    if f.p <= f.q:
        return False
    if ((f.p + f.q) & 1) == 0:
        return False
    if math.gcd(f.p, f.q) != 1:
        return False
    return f.p * f.p + f.q * f.q <= L

def find_valid_lower(lo, hi, N, L):
    if lo.p <= lo.q:
        return Fraction(0, 0)
    for _ in range(100000):
        if is_valid_fraction(lo, L):
            return lo
        prev = farey_prev(lo, hi, N)
        hi = lo
        lo = prev
        if lo.q == 0:
            break
    return Fraction(0, 0)

def find_valid_upper(lo, hi, N, L):
    for _ in range(100000):
        if is_valid_fraction(hi, L):
            return hi
        nxt = farey_next(lo, hi, N)
        lo = hi
        hi = nxt
        if hi.q == 0:
            break
    return Fraction(0, 0)

def evaluate_fraction(f, t_alpha, L):
    if f.p == 0 or f.q == 0:
        return None
    m = f.p
    n = f.q
    m2 = m * m
    n2 = n * n
    a = m2 - n2
    b = 2 * m * n
    c = m2 + n2
    k = L // c
    if k == 0:
        return None

    t = 3.0 * a * b / (2.0 * c * c)
    diff = abs(t - t_alpha) / (1.0 + t * t_alpha)

    sum_val = k * (a + b + c)
    area = k * k * a * b
    return {'diff': diff, 'sum': sum_val, 'area': area}

def best_candidate_for_u(u_target, t_alpha, L):
    if u_target <= 1.0:
        return None

    R = math.sqrt(L)
    denom = math.sqrt(1.0 + u_target * u_target)
    N = int(math.floor(R / denom))
    if N < 1:
        N = 1

    lo, hi = farey_bounds(u_target, N)
    lo_valid = find_valid_lower(lo, hi, N, L)
    hi_valid = find_valid_upper(lo, hi, N, L)

    best = None
    if lo_valid.q != 0:
        best = evaluate_fraction(lo_valid, t_alpha, L)

    if hi_valid.q != 0:
        cand = evaluate_fraction(hi_valid, t_alpha, L)
        if best is None:
            best = cand
        elif cand is not None:
            eps = 1e-15
            if cand['diff'] + eps < best['diff']:
                best = cand
            elif abs(cand['diff'] - best['diff']) <= eps:
                if cand['area'] > best['area']:
                    best = cand

    return best

def compute_f(alpha_deg, L):
    kPi = math.acos(-1.0)
    kUMin = 1.0 + math.sqrt(2.0)
    kTMax = 0.75

    alpha_rad = alpha_deg * kPi / 180.0
    t_alpha = math.tan(alpha_rad)

    u_low = 1.0
    u_high = kUMin
    if t_alpha < kTMax:
        low = 1.0
        high = kUMin
        for _ in range(100):
            mid = 0.5 * (low + high)
            if tan_theta_from_u(mid) < t_alpha:
                low = mid
            else:
                high = mid
        u_low = 0.5 * (low + high)

        low = kUMin
        high = max(kUMin + 1.0, 3.0 / t_alpha + 2.0)
        while tan_theta_from_u(high) > t_alpha:
            high *= 2.0
        for _ in range(100):
            mid = 0.5 * (low + high)
            if tan_theta_from_u(mid) > t_alpha:
                low = mid
            else:
                high = mid
        u_high = 0.5 * (low + high)

    best = best_candidate_for_u(u_low, t_alpha, L)
    alt = best_candidate_for_u(u_high, t_alpha, L)

    if best is None:
        best = alt
    elif alt is not None:
        eps = 1e-15
        if alt['diff'] + eps < best['diff']:
            best = alt
        elif abs(alt['diff'] - best['diff']) <= eps:
            if alt['area'] > best['area']:
                best = alt

    return best['sum'] if best else 0

def solve():
    N = 45000
    L = 10000000000
    
    total = 0
    for i in range(1, N + 1):
        alpha = i ** (1/3)
        total += compute_f(alpha, L)
        
    return str(total)

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

Java

public class Euler904 {

    static class Fraction {
        long p;
        long q;

        Fraction(long p, long q) {
            this.p = p;
            this.q = q;
        }
    }

    static class Candidate {
        boolean valid = false;
        double diff = 0.0;
        long sum = 0;
        java.math.BigInteger area = java.math.BigInteger.ZERO;
    }

    static double tanThetaFromU(double u) {
        double u2 = u * u;
        return 3.0 * u * (u2 - 1.0) / ((u2 + 1.0) * (u2 + 1.0));
    }

    static Fraction fareyPrev(Fraction lo, Fraction hi, long N) {
        long k = (N + hi.q) / lo.q;
        return new Fraction(k * lo.p - hi.p, k * lo.q - hi.q);
    }

    static Fraction fareyNext(Fraction lo, Fraction hi, long N) {
        long k = (N + lo.q) / hi.q;
        return new Fraction(k * hi.p - lo.p, k * hi.q - lo.q);
    }

    static Fraction[] fareyBounds(double x, long N) {
        java.util.List<Long> a = new java.util.ArrayList<>();
        double y = x;
        for (int i = 0; i < 64; ++i) {
            long ai = (long) Math.floor(y);
            a.add(ai);
            double frac = y - ai;
            if (Math.abs(frac) < 1e-15) {
                break;
            }
            y = 1.0 / frac;
        }

        java.util.List<Long> ps = new java.util.ArrayList<>();
        java.util.List<Long> qs = new java.util.ArrayList<>();
        long pMinus2 = 0;
        long pMinus1 = 1;
        long qMinus2 = 1;
        long qMinus1 = 0;

        for (long ai : a) {
            long p = ai * pMinus1 + pMinus2;
            long q = ai * qMinus1 + qMinus2;
            ps.add(p);
            qs.add(q);
            if (q > N) {
                break;
            }
            pMinus2 = pMinus1;
            pMinus1 = p;
            qMinus2 = qMinus1;
            qMinus1 = q;
        }

        int idx = 0;
        while (idx < qs.size() && qs.get(idx) <= N) {
            ++idx;
        }

        if (idx == 0) {
            return new Fraction[] { new Fraction(0, 1), new Fraction(1, 0) };
        }

        if (idx == qs.size()) {
            return new Fraction[] { new Fraction(ps.get(idx - 1), qs.get(idx - 1)),
                    new Fraction(ps.get(idx - 1), qs.get(idx - 1)) };
        }

        long pPrev = ps.get(idx - 1);
        long qPrev = qs.get(idx - 1);
        long pPrevPrev = (idx >= 2) ? ps.get(idx - 2) : 1;
        long qPrevPrev = (idx >= 2) ? qs.get(idx - 2) : 0;

        long t = (N - qPrevPrev) / qPrev;
        long pCand = pPrevPrev + t * pPrev;
        long qCand = qPrevPrev + t * qPrev;

        if ((double) pPrev / qPrev < x) {
            return new Fraction[] { new Fraction(pPrev, qPrev), new Fraction(pCand, qCand) };
        } else {
            return new Fraction[] { new Fraction(pCand, qCand), new Fraction(pPrev, qPrev) };
        }
    }

    static long gcd(long a, long b) {
        return b == 0 ? a : gcd(b, a % b);
    }

    static boolean isValidFraction(Fraction f, long L) {
        if (f.q <= 0 || f.p <= 0)
            return false;
        if (f.p <= f.q)
            return false;
        if (((f.p + f.q) & 1L) == 0)
            return false;
        if (gcd(f.p, f.q) != 1)
            return false;
        long m = f.p;
        long n = f.q;
        if (m > Math.sqrt(L))
            return false;
        return m * m + n * n <= L;
    }

    static Fraction findValidLower(Fraction lo, Fraction hi, long N, long L) {
        if (lo.p <= lo.q)
            return new Fraction(0, 0);
        for (int steps = 0; steps < 100000; ++steps) {
            if (isValidFraction(lo, L))
                return lo;
            Fraction prev = fareyPrev(lo, hi, N);
            hi = lo;
            lo = prev;
            if (lo.q == 0)
                break;
        }
        return new Fraction(0, 0);
    }

    static Fraction findValidUpper(Fraction lo, Fraction hi, long N, long L) {
        for (int steps = 0; steps < 100000; ++steps) {
            if (isValidFraction(hi, L))
                return hi;
            Fraction next = fareyNext(lo, hi, N);
            lo = hi;
            hi = next;
            if (hi.q == 0)
                break;
        }
        return new Fraction(0, 0);
    }

    static Candidate evaluateFraction(Fraction f, double tAlpha, long L) {
        Candidate cand = new Candidate();
        if (f.p == 0 || f.q == 0)
            return cand;
        long m = f.p;
        long n = f.q;
        long m2 = m * m;
        long n2 = n * n;
        long a = m2 - n2;
        long b = 2L * m * n;
        long c = m2 + n2;
        long k = L / c;
        if (k == 0)
            return cand;

        double t = 3.0 * a * b / (2.0 * c * c);
        double diff = Math.abs(t - tAlpha) / (1.0 + t * tAlpha);

        cand.valid = true;
        cand.diff = diff;
        cand.sum = k * (a + b + c);
        java.math.BigInteger bK = java.math.BigInteger.valueOf(k);
        java.math.BigInteger bA = java.math.BigInteger.valueOf(a);
        java.math.BigInteger bB = java.math.BigInteger.valueOf(b);
        cand.area = bK.multiply(bK).multiply(bA).multiply(bB);
        return cand;
    }

    static Candidate bestCandidateForU(double uTarget, double tAlpha, long L) {
        Candidate best = new Candidate();
        if (uTarget <= 1.0)
            return best;

        double R = Math.sqrt(L);
        double denom = Math.sqrt(1.0 + uTarget * uTarget);
        long N = (long) Math.floor(R / denom);
        if (N < 1)
            N = 1;

        Fraction[] bounds = fareyBounds(uTarget, N);
        Fraction lo = bounds[0];
        Fraction hi = bounds[1];

        Fraction loValid = findValidLower(lo, hi, N, L);
        Fraction hiValid = findValidUpper(lo, hi, N, L);

        if (loValid.q != 0) {
            best = evaluateFraction(loValid, tAlpha, L);
        }
        if (hiValid.q != 0) {
            Candidate cand = evaluateFraction(hiValid, tAlpha, L);
            if (!best.valid) {
                best = cand;
            } else if (cand.valid) {
                double eps = 1e-15;
                if (cand.diff + eps < best.diff) {
                    best = cand;
                } else if (Math.abs(cand.diff - best.diff) <= eps) {
                    if (cand.area.compareTo(best.area) > 0) {
                        best = cand;
                    }
                }
            }
        }

        return best;
    }

    static long computeF(double alphaDeg, long L) {
        double kPi = Math.acos(-1.0);
        double kUMin = 1.0 + Math.sqrt(2.0);
        double kTMax = 0.75;

        double alphaRad = alphaDeg * kPi / 180.0;
        double tAlpha = Math.tan(alphaRad);

        double uLow = 1.0;
        double uHigh = kUMin;
        if (tAlpha < kTMax) {
            double low = 1.0;
            double high = kUMin;
            for (int iter = 0; iter < 100; ++iter) {
                double mid = 0.5 * (low + high);
                if (tanThetaFromU(mid) < tAlpha) {
                    low = mid;
                } else {
                    high = mid;
                }
            }
            uLow = 0.5 * (low + high);

            low = kUMin;
            high = Math.max(kUMin + 1.0, 3.0 / tAlpha + 2.0);
            while (tanThetaFromU(high) > tAlpha) {
                high *= 2.0;
            }
            for (int iter = 0; iter < 100; ++iter) {
                double mid = 0.5 * (low + high);
                if (tanThetaFromU(mid) > tAlpha) {
                    low = mid;
                } else {
                    high = mid;
                }
            }
            uHigh = 0.5 * (low + high);
        }

        Candidate best = bestCandidateForU(uLow, tAlpha, L);
        Candidate alt = bestCandidateForU(uHigh, tAlpha, L);

        if (!best.valid) {
            best = alt;
        } else if (alt.valid) {
            double eps = 1e-15;
            if (alt.diff + eps < best.diff) {
                best = alt;
            } else if (Math.abs(alt.diff - best.diff) <= eps) {
                if (alt.area.compareTo(best.area) > 0) {
                    best = alt;
                }
            }
        }

        return best.sum;
    }

    public static String solve() {
        int N = 45000;
        long L = 10000000000L;

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        long[] partial = new long[threads];
        Thread[] pool = new Thread[threads];

        for (int t = 0; t < threads; ++t) {
            final int tIdx = t;
            pool[t] = new Thread(() -> {
                long local = 0;
                for (int i = 1 + tIdx; i <= N; i += threads) {
                    double alpha = Math.cbrt(i);
                    local += computeF(alpha, L);
                }
                partial[tIdx] = local;
            });
            pool[t].start();
        }

        long total = 0;
        try {
            for (int t = 0; t < threads; t++) {
                if (pool[t] != null)
                    pool[t].join();
                total += partial[t];
            }
        } catch (InterruptedException e) {
            e.printStackTrace();
        }

        return Long.toString(total);
    }

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