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

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