Problem 385: Ellipses Inside Triangles

Project Euler Problem 385 Solution

EulerSolve provides an optimized solution for Project Euler Problem 385, Ellipses Inside Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The local C++ validator shows that every admissible triangle can be represented in centered form $$P_1=(x_1,y_1),\qquad P_2=(x_2,y_2),\qquad P_3=(-x_1-x_2,-y_1-y_2),$$ so automatically $$x_1+x_2+x_3=0,\qquad y_1+y_2+y_3=0.$$ The ellipse condition is already converted by the implementation into the Diophantine system $$ (2x_1+x_2)y_1 + (x_1+2x_2)y_2 = 0, $$ $$ x_1^2+x_1x_2+x_2^2 - (y_1^2+y_1y_2+y_2^2) = 39. $$ We must sum the Euclidean areas of all non-degenerate integer triangles whose coordinates lie inside $[-N,N]^2$, for the real target $N=10^9$. The brute-force search in solve_bruteforce_small is only feasible for tiny $N$, so the fast solver rewrites the problem as a finite family of Pell-type recurrences. Mathematical Approach 1. Extract a primitive direction from the $x$-coordinates Write the three $x$-coordinates as $$x_1=du,\qquad x_2=dv,\qquad x_3=-d(u+v),$$ where $d>0$ and $\gcd(u,v)=1$. The linear constraint then becomes $$a y_1 + b y_2 = 0,\qquad a=2u+v,\qquad b=u+2v.$$ All integer solutions of this linear equation have the form $$y_1=\frac{b}{m}k,\qquad y_2=-\frac{a}{m}k,\qquad y_3=-y_1-y_2=\frac{u-v}{m}k,$$ where $m=\gcd(a,b)$ and $k\in\mathbb{Z}$. Because $$2a-b=3u,\qquad 2b-a=3v,$$ every common divisor of $a$ and $b$ must divide $3$. Since $(u,v)$ is primitive, $m\in\{1,3\}$....

Detailed mathematical approach

Problem Summary

Inside every triangle of the plane there is exactly one ellipse of largest area, which the statement grants us. Fix an integer $n\ge1$ and consider the triangles whose three vertices have integer coordinates of absolute value at most $n$ and whose largest inscribed ellipse has its two foci — the points whose distances to any boundary point have a constant sum — at $(\pm\sqrt{13},0)$. Let $A(n)$ be the sum of their areas; we must compute $A(10^9)$. ^[1] For $n=8$ there are two, with vertices $(-4,-3),(-4,3),(8,0)$ and $(4,3),(4,-3),(-8,0)$, each of area $36$, so $A(8)=72$. The statement also gives $A(10)=252$, $A(100)=34632$ and $A(1000)=3529008$.

The whole solution in one sentence: the condition on the foci becomes two polynomial equations in the six coordinates, those collapse into two Pell equations whose solutions grow like powers of $2+\sqrt{3}$, so for $n=10^9$ only $116$ triangles qualify and every one of them can be written down.

Mathematical Approach

1. The largest ellipse inside a triangle

An invertible affine map $p\mapsto Lp+t$ sends ellipses to ellipses and multiplies every area by $|\det L|$, so it carries the largest ellipse inside a triangle to the largest ellipse inside the image triangle. Every triangle is an affine image of the equilateral triangle $E$, so it suffices to find the largest ellipse inside $E$.

The symmetries of $E$ permute the ellipses inside $E$ and preserve areas, so each maps the unique largest one to itself. The rotation $\rho$ by $120^\circ$ about the centre of $E$ therefore maps the maximiser to itself, hence its centre to its centre; as $\rho$ fixes only one point, the maximiser is centred there. A non-circular ellipse has just two axes of symmetry and they are perpendicular, so its rotational symmetries have order at most $2$; since $\rho$ has order $3$, the maximiser is a circle, and the largest such circle is the incircle.

The largest ellipse inside any triangle is therefore the affine image of that incircle. The incircle touches each side of $E$ at its midpoint and affine maps preserve midpoints, so this ellipse is the Steiner inellipse: the unique inscribed ellipse tangent to the sides at their midpoints. Its centre is the centroid and its area is $\pi/(3\sqrt{3})\approx0.6046$ times the triangle's. ^[2]

2. The shape matrix of that ellipse

Place $E$ with centroid at the origin and vertices on the unit circle:

$$\begin{aligned}e_1&=(1,0),\\ e_{2,3}&=\bigl(-\tfrac12,\pm\tfrac{\sqrt3}{2}\bigr).\end{aligned}$$

Its incircle has radius $\tfrac12$, and a direct check gives

$$\sum_{j=1}^{3}e_je_j^{\mathsf T}=\tfrac32 I.$$

Let $T$ have vertices $P_1,P_2,P_3$, centroid $G$ and $V_j=P_j-G$, and let $L$ be the linear map with $Le_j=V_j$, consistent because $\sum e_j=\sum V_j=0$. Then

$$\begin{aligned}\sum_j V_jV_j^{\mathsf T}&=L\Bigl(\sum_j e_je_j^{\mathsf T}\Bigr)L^{\mathsf T}\\ &=\tfrac32LL^{\mathsf T}.\end{aligned}$$

The incircle of $E$ is $\{\tfrac12w:|w|\le1\}$, and the affine map also adds $G$, so its image is

$$\bigl\{G+\tfrac12Lw:|w|\le1\bigr\},$$

an ellipse with centre $G$ and shape matrix

$$Q=\tfrac14LL^{\mathsf T}=\frac16\sum_j V_jV_j^{\mathsf T}.$$

Such an ellipse is the set of points $p$ with $(p-G)^{\mathsf T}Q^{-1}(p-G)\le1$. If $v$ is a unit eigenvector of $Q$, so $Qv=\lambda v$, then $p-G=tv$ turns that into $t^2/\lambda\le1$: the semi-axis along $v$ is $\sqrt\lambda$. A symmetric positive-definite matrix has two perpendicular eigenvector directions and positive eigenvalues, so the eigenvalues of $Q$ are the squared semi-axes; and $w^{\mathsf T}Qw=\tfrac16\sum_j(V_j^{\mathsf T}w)^2$ is positive for $w\ne0$ exactly when the $V_j$ span the plane, that is when $T$ is non-degenerate. For semi-axes $a\ge b$ the two focal distances sum to $2a$ at an end of the major axis and $2\sqrt{b^2+c^2}$ at an end of the minor one, so $c^2=a^2-b^2$ and

$$c^2=\lambda_{\max}-\lambda_{\min}.$$

3. Two equations for the six coordinates

The foci are to be $(\pm\sqrt{13},0)$. Their midpoint is the centre of the ellipse, so $G=0$, which means

$$x_1+x_2+x_3=0,\qquad y_1+y_2+y_3=0.$$

With $G=0$ we have $V_j=P_j$ and

$$6Q=\begin{pmatrix}\sum x_j^2&\sum x_jy_j\\[2pt]\sum x_jy_j&\sum y_j^2\end{pmatrix}.$$

The $x$-axis is an axis of the ellipse exactly when $Q$ is diagonal, that is when $\sum x_jy_j=0$; the eigenvalues are then $Q_{11}$ and $Q_{22}$, and the major axis is the $x$-axis with $c^2=13$ exactly when $Q_{11}-Q_{22}=13$. The geometric condition is therefore equivalent to

$$\begin{aligned}\sum_j x_jy_j&=0,\\ \sum_j \bigl(x_j^2-y_j^2\bigr)&=78.\end{aligned}$$

Both are symmetric in the vertices, as they must be. Eliminating $P_3=-P_1-P_2$ and writing

$$q(s,t)=s^2+st+t^2,$$

one finds $\sum x_j^2=2q(x_1,x_2)$, the same for the $y$-part, and the linear form below for $\sum x_jy_j$. The system becomes

$$(2x_1+x_2)y_1+(x_1+2x_2)y_2=0,$$

$$q(x_1,x_2)-q(y_1,y_2)=39.$$

Check the triangle $(4,3),(4,-3),(-8,0)$ of the statement: $\sum x_jy_j=12-12+0=0$ and $\sum(x_j^2-y_j^2)=7+7+64=78$. Its shape matrix is $Q=\operatorname{diag}(16,3)$, so the semi-axes are $4$ and $\sqrt{3}$ and indeed $16-3=13$.

4. The same conditions from Marden's theorem

Read the vertices as complex numbers $z_j=x_j+iy_j$ and put $p(z)=(z-z_1)(z-z_2)(z-z_3)$. Marden's theorem, found by Jörg Siebeck in 1864, says that the foci of the Steiner inellipse are the zeros of $p'$. ^[3] Demanding $p'(z)=3z^2-39$ makes the coefficient of $z^2$ in $p$ vanish and that of $z$ equal $-39$, that is $\sum z_j=0$ and $\sum_{i\lt j}z_iz_j=-39$, hence

$$\sum_j z_j^2=-2\sum_{i\lt j}z_iz_j=78.$$

Its real part is $\sum(x_j^2-y_j^2)=78$ and its imaginary part $2\sum x_jy_j=0$: the equations of Section 3. Both routes agree.

5. A primitive direction in the $x$-coordinates

Since $4q(s,t)=(2s+t)^2+3t^2$, the form $q$ is positive definite, so $q(y_1,y_2)\ge0$ and therefore $q(x_1,x_2)=39+q(y_1,y_2)\ge39$. In particular $(x_1,x_2)\ne(0,0)$, so we may set

$$d=\gcd(x_1,x_2)\gt0,\quad x_1=du,\quad x_2=dv,$$

with $\gcd(u,v)=1$. Abbreviate $a=2u+v$ and $b=u+2v$. From $2a-b=3u$ and $2b-a=3v$ we get $(a,b)\ne(0,0)$, and every common divisor of $a,b$ divides $3u$ and $3v$, hence divides $3$: so $m=\gcd(a,b)$ is $1$ or $3$. Dividing the linear equation by $d$ leaves $ay_1+by_2=0$, with integer solutions exactly

$$y_1=\frac bm\,k,\quad y_2=-\frac am\,k,\quad k\in\mathbb Z,$$

and then $y_3=-y_1-y_2=\frac{u-v}{m}k$, because $a-b=u-v$.

6. The quadratic equation in $d$ and $k$

Write $q_0=q(u,v)$, so $q(x_1,x_2)=d^2q_0$. Expanding gives the identity

$$q(b,-a)=a^2-ab+b^2=3q_0,$$

so $q(y_1,y_2)=3q_0k^2/m^2$ and the quadratic equation becomes

$$q_0\Bigl(d^2-\frac{3k^2}{m^2}\Bigr)=39.$$

Which value does $m$ take? Because $q(u,v)=(u-v)^2+3uv$, we have $q_0\equiv(u-v)^2\pmod 3$, so $3\mid q_0$ exactly when $u\equiv v\pmod3$; and $3\mid a$ means $v\equiv-2u\equiv u\pmod 3$, which also forces $3\mid b$. Hence

$$m=3\iff 3\mid q_0,\qquad m=1\iff 3\nmid q_0.$$

Suppose $m=1$. Then $q_0(d^2-3k^2)=39$ with $3\nmid q_0$, so $q_0\in\{1,13\}$ and $d^2-3k^2\in\{39,3\}$. Reducing modulo $3$ gives $3\mid d$; writing $d=3d'$ and dividing by $3$ turns the two cases into $3d'^2-k^2=13$ and $3d'^2-k^2=1$, and modulo $3$ both demand $k^2\equiv2\pmod3$, which no square satisfies. So $m=1$ is impossible and always $m=3$, leaving

$$q_0\,(3d^2-k^2)=117.$$

7. Only two norms survive

The last equation forces $q_0\mid117=3^2\cdot13$ with $3\mid q_0$, so $q_0\in\{3,9,39,117\}$. A norm $q(u,v)$ with $\gcd(u,v)=1$ is never divisible by $9$. Indeed, Section 6 gives $u\equiv v\pmod3$; writing $v=u+3t$ turns it into

$$q(u,v)=3u^2+9ut+9t^2,$$

and $3\nmid u$, since $3\mid u$ would force $3\mid v$ against $\gcd(u,v)=1$. Hence $u^2\equiv1$ and $q_0\equiv3\pmod 9$, ruling out $9$ and $117$. ^[4] Both survivors occur, since $q(1,1)=3$ and $q(2,5)=39$.

From $4q_0=(2v+u)^2+3u^2\ge3u^2$ we get $|u|,|v|\le2\sqrt{q_0/3}$, so a scan over $|u|,|v|\le7$ is exhaustive. It finds $6$ primitive pairs of norm $3$, namely

$$(1,1),\ (1,-2),\ (-2,1)$$

and their negatives, and $12$ of norm $39$, namely

$$(2,5),\ (5,2),\ (2,-7),\ (-7,2),\ (5,-7),\ (-7,5)$$

and their negatives. The whole infinite search has shrunk to $18$ directions.

8. Two Pell equations

The two norms give $3d^2-k^2=39$ and $3d^2-k^2=3$. In both cases $3\mid k^2$, hence $3\mid k$; put $k=3\kappa$ and divide by $3$:

$$d^2-3\kappa^2=13\quad(q_0=3),$$

$$d^2-3\kappa^2=1\quad(q_0=39).$$

The second is Pell's equation for $\sqrt{3}$, with fundamental solution $(d,\kappa)=(2,1)$. ^[5] Multiplying $d+\kappa\sqrt3$ by the unit $2+\sqrt3$ preserves $d^2-3\kappa^2$, giving the step

$$(d,\kappa)\mapsto(2d+3\kappa,\ d+2\kappa),$$

or, in the variables $d$ and $k=3\kappa$ used by the programs,

$$(d,k)\mapsto(k+2d,\ 2k+3d).$$

A descent shows that this one step generates every solution. Let $d\gt0$, $\kappa\ge0$ and $d^2-3\kappa^2=N$ with $N\in\{1,13\}$. The inverse step sends $(d,\kappa)$ to $(2d-3\kappa,\,2\kappa-d)$, again a solution; here $2d-3\kappa\gt0$ because $4d^2=12\kappa^2+4N\gt9\kappa^2$, and $2\kappa-d\lt\kappa$ because $d\gt\kappa$. Also $2\kappa-d\ge0$ as soon as $4\kappa^2\ge3\kappa^2+N$, that is $\kappa^2\ge N$. So from any solution with $\kappa^2\ge N$ we step down to a strictly smaller non-negative $\kappa$, until $\kappa^2\lt N$. For $N=1$ only $\kappa=0,d=1$ remains; for $N=13$ only $\kappa\le3$, of which just $\kappa=1,2$ make $3\kappa^2+13$ a square, giving $d=4$ and $d=5$.

Back in $(d,k)$, the admissible pairs form three chains, two for $3d^2-k^2=39$ ^[6] and one for $3d^2-k^2=3$ ^[7]:

$(4,3),(11,18),(40,69),$ $(149,258),(556,963),\ldots$
$(5,6),(16,27),(59,102),$ $(220,381),(821,1422),\ldots$
$(1,0),(2,3),(7,12),$ $(26,45),(97,168),\ldots$

The first member $(1,0)$ of the last chain has $k=0$, so all three $y_j$ vanish and the points lie on the $x$-axis; it is the only solution with $k=0$ and must be discarded. Every other pair gives a genuine triangle.

9. Which pairs fit in the box

With $m=3$ and $k=3\kappa$ the vertices of the triangle belonging to a direction $(u,v)$ and a pair $(d,k)$ are

$$(x_1,x_2,x_3)=d\,(u,\ v,\ -u-v),$$

$$(y_1,y_2,y_3)=\kappa\,(b,\ -a,\ u-v).$$

Hence $\max_j|x_j|=dX$ with $X=\max(|u|,|v|,|u+v|)$ and $\max_j|y_j|=|\kappa|\,Y$ with $Y=\max(|a|,|b|,|u-v|)$. All six directions of norm $3$ have $X=2,Y=3$, and all twelve of norm $39$ have $X=7,Y=12$. So the triangle lies in $[-n,n]^2$ exactly when

$$q_0=3:\quad 2d\le n,\quad k\le n,$$

$$q_0=39:\quad 7d\le n,\quad 4k\le n.$$

Along every chain $k^2=3d^2-c$ with $c\gt0$, so $k\lt\sqrt3\,d$; therefore $k\lt2d$ and $4k\lt4\sqrt3\,d\lt7d$. In both families the bound on $d$ is the binding one, and the bound on $k$ follows.

10. The area, and a worked example

The shoelace formula ^[8] gives twice the area as the absolute value of

$$x_1(y_2-y_3)+x_2(y_3-y_1)+x_3(y_1-y_2).$$

Substituting the coordinates above and using $a+u-v=3u$, $b-u+v=3v$ and $a+b=3(u+v)$, the summands become $-3d\kappa u^2$, $-3d\kappa v^2$ and $-3d\kappa(u+v)^2$, with sum $-6d\kappa\,q_0$. So the area is

$$3q_0\,d\,|\kappa|=q_0\,d\,|k|.$$

Take $q_0=3$, $(u,v)=(1,1)$, so $a=b=3$ and $u-v=0$, with $(d,k)=(4,3)$, that is $\kappa=1$. The vertices are $(4,3)$, $(4,-3)$, $(-8,0)$ and the area is $3\cdot4\cdot3=36$: one of the two triangles of $A(8)=72$, the direction $(-1,-1)$ giving the other.

Take $q_0=39$, $(u,v)=(2,5)$, so $a=9$, $b=12$, $u-v=-3$, with $(d,k)=(2,3)$ and $\kappa=1$. The vertices are $(4,12)$, $(10,-9)$, $(-14,-3)$ and the area is $39\cdot2\cdot3=234$; the shoelace formula gives $468/2=234$. As a check, $\sum x_jy_j=48-90+42=0$ and $\sum(x_j^2-y_j^2)=-128+19+187=78$.

11. From pairs to triangles

Fix a pair $(d,k)$ with $k\gt0$. Its $6$ directions of norm $3$, or $12$ of norm $39$, with the two signs of $k$, give $12$ or $24$ labelled triples, all distinct: the direction fixes $(x_1,x_2)=(du,dv)$, and flipping the sign of $k$ changes $(y_1,y_2)$ because $(a,b)\ne(0,0)$.

Conversely, relabelling changes none of $d$, $q_0$ and $|k|$: the number $d=\gcd(x_1,x_2,x_3)$ is symmetric in the vertices, and then $\sum x_j^2=2d^2q_0$ and $\sum y_j^2=\tfrac23q_0k^2$ fix $q_0$ and $|k|$. So all $6$ labellings of a triangle built from $(d,k)$ are among those triples, which therefore split into classes of $6$: each pair gives exactly $2$ triangles when $q_0=3$ and $4$ when $q_0=39$, all of area $q_0dk$.

Geometrically they are images of one another under $p\mapsto\bar p$, the reflection in the $x$-axis, and $p\mapsto-p$. For $q_0=3$ the directions $(1,1),(1,-2),(-2,1)$ all produce

$$\{(d,3\kappa),(d,-3\kappa),(-2d,0)\},$$

which is symmetric about the $x$-axis, while the other three produce its point reflection. Writing $S_3(n)$ for the pairs of the first two chains with $2d\le n$ and $S_{39}(n)$ for those of the third with $7d\le n$ and $k\gt0$:

$$\begin{aligned}A(n)=\;&6\sum_{(d,k)\in S_3(n)}dk\\ +\;&156\sum_{(d,k)\in S_{39}(n)}dk.\end{aligned}$$

For $n=10$ only $(4,3)$ and $(5,6)$ qualify and $S_{39}$ is empty, so $A(10)=6(12+30)=252$. For $n=100$ the set $S_3$ holds $(4,3),(11,18),(40,69)$ and $(5,6),(16,27)$, with $\sum dk=3432$, while $S_{39}$ holds $(2,3)$ and $(7,12)$, with $\sum dk=90$; so $A(100)=6\cdot3432+156\cdot90=34632$. Both match the statement, as does

$$A(1000)=6\cdot131712+156\cdot17556=3529008.$$

At $n=10^9$ the two sets have $30$ and $14$ elements, so exactly $2\cdot30+4\cdot14=116$ triangles contribute.

How the Code Works

All three programs run the same enumeration. The C++ entry point ^[9] parses --n= and --skip-checkpoints, runs run_checkpoints by default, calls solve_fast and prints the total, with the triangle count on standard error; Python and Java call solve() with the bound fixed at $10^9$ and carry no checkpoints. The names below are the C++ ones, which Python shares and Java camel-cases.

build_directions(n) implements Sections 7 and 9: for each $q_0\in\{3,39\}$ it sets pell_n to -117/q0, that is $-39$ and $-3$, the right-hand sides written as $k^2-3d^2$; it scans $|u|,|v|$ up to int(sqrt(q0))+2, keeps the pairs with $\gcd(u,v)=1$, $q(u,v)=q_0$ and $m=\gcd(2u+v,u+2v)=3$, and computes the local scales x_scale $=X$ and y_scale $=Y/3$. A kept record holds $q_0,u,v,m$, the Pell constant and the bounds d_max $=\lfloor n/X\rfloor$ and k_max $=\lfloor 3n/Y\rfloor$, not the scales; a record with a zero bound is dropped, so the list has all $18$ directions once $n\ge7$ and only the six of norm $3$ for $2\le n\le6$. The scan radii $3$ and $8$ exceed the true bounds $2$ and $7$, so both scans are exhaustive.

generate_pell_pairs(pell_n, d_max, k_max) is Section 8. Its seeds are hard-coded — $(d,k)=(4,3)$ and $(5,6)$ for $-39$, $(2,3)$ for $-3$ — so the degenerate $(1,0)$ is never produced. While both coordinates are inside d_max and k_max it appends $(d,k)$, computes next_k = 2*k + 3*d and next_d = k + 2*d from the old values and replaces both at once, so $(4,3)$ becomes $(11,18)$; both grow, so stopping at the first violation loses nothing. Each version then sorts and deduplicates its list. C++ and Java guard the next step against overflow — safe at the target size, though not for an arbitrary --n=, since 3*d inside the guard can itself overflow.

build_triangle_key(dir, d, k, n) rebuilds the vertices as in Section 9 — d*u, d*v, -d*(u+v) and (b/m)*k, -(a/m)*k, ((u-v)/m)*k — rejects a triple that leaves the box or has two equal vertices, and returns the sorted triple with the area q0*abs(d*k) of Section 10. Its divisibility tests on a, b and u-v are redundant, since $m=\gcd(a,b)=3$ divides all three.

solve_fast(n) takes the largest d_max and k_max of each family, builds the two pair lists once, and then constructs a triangle for every direction, every pair inside that direction's bounds and both signs of $k$, using the sorted vertex triple as a hash-map key. That is the deduplication of Section 11: each triangle is constructed six times, once per labelling. At $n=10^9$ that makes $6\cdot30+12\cdot14=348$ direction-and-pair combinations, $696$ candidates and $116$ distinct keys. C++ and Python store a key only when absent, Java rewrites the equal area, and all three stop if one key ever carries two different areas — a live check on Section 10.

solve_bruteforce_small(n), in the C++ file only, is an independent cubic-time search: for every $(x_1,y_1,x_2)$ it solves the linear equation for $y_2$, treating the vanishing coefficient $x_1+2x_2$ separately, then tests the quadratic equation, the box, distinctness and the shoelace area. It builds $x_3,y_3$ as $-x_1-x_2,-y_1-y_2$, so it checks Sections 5 to 11 and not the geometry of Sections 1 to 3. run_checkpoints compares the fast method against the four given values of $A(n)$ and against this brute force at $n=20$.

The three languages differ mainly in integer width, and a crude bound shows why. Expanding the recurrence, consecutive products satisfy $d'k'/(dk)=7+6d/k+2k/d$, and $6d/k+2k/d\ge4\sqrt3$ because $(\sqrt{6d/k}-\sqrt{2k/d})^2\ge0$, so the ratio exceeds $13$; a chain's sum is thus below $\tfrac{13}{12}$ of its last product, which is below $\sqrt3(n/X)^2$. Two chains of weight $6$ with $X=2$ and one of weight $156$ with $X=7$ give $A(n)\lt\sqrt3\,(13/4+169/49)\,n^2\lt12n^2$. At $n=10^9$ that is $1.2\cdot10^{19}$, above the signed 64-bit maximum $2^{63}-1\approx9.22\cdot10^{18}$, which is why the C++ accumulates in unsigned __int128 while Python's integers are unbounded. The Java port sums into a long unchecked; the true total stays in range at the target size, but not from $n=1\,328\,254\,382$ on, where Java would overflow silently. Single areas are safe, since $3dk\lt1.5n^2$ and $39dk\lt1.4n^2$.

Complexity and Verification

Building the directions scans $7^2+17^2=338$ pairs $(u,v)$: constant work. Each chain multiplies $d$ by roughly $2+\sqrt3\approx3.732$ per step, so the number $M$ of pairs is $O(\log n)$; at $n=10^9$ each chain of norm $3$ holds $15$ and the chain of norm $39$ holds $14$, so $M=44$. Sorting the lists costs $O(M\log M)$ comparisons and the enumeration $O(M)$ hash-map operations of expected constant cost, so with unit-cost arithmetic the total time is $O(\log n\log\log n)$ and the storage $O(\log n)$ records, the final map holding $116$ entries. All three programs therefore finish in a small fraction of a second and print the same value; in the C++ run with checkpoints the cubic brute force at $n=20$ dominates the time.

The following checks were run for this article, independently of the packaged ones. Maximising the area of an inscribed ellipse numerically — over its centre and its symmetric $2\times2$ shape matrix, under the three tangency inequalities — returns for $(4,3),(4,-3),(-8,0)$ the semi-axes $4$ and $\sqrt3$ with the foci $(\pm\sqrt{13},0)$, for $(5,6),(5,-6),(-10,0)$ the semi-axes $5$ and $2\sqrt3$ with the same foci, and for an equilateral triangle the incircle; the area ratio is $\pi/(3\sqrt3)$ throughout, as Section 1 predicts. The foci computed from $Q$ also agree with the zeros of $p'$, confirming Marden's theorem numerically.

A fully independent enumeration over every unordered triple of lattice points in $[-n,n]^2$, with no assumption about the centroid and the foci computed from $Q$, returns exactly the two triangles of the statement for $n=8$, then $A(10)=252$, $A(12)=252$ and $A(20)=1188$. A brute force over the two Diophantine equations agrees with the fast method for every $n\le60$. The closed formula of Section 11 agrees with building and deduplicating the triangles of the eighteen directions for every $n\le400$ and for $n=1000,2000,5000,12345$ and $99999$; that comparison shares the Pell parametrisation, so it tests the weights $6$ and $156$ and the duplicate handling, not the geometry. A search over $0\lt d\le200\,000$ with $k\gt0$ finds no solution of the two Pell equations outside the three chains, checking the seeds in that range; the descent of Section 8 is what proves them complete. Finally $A(n)\lt12n^2$ holds in every tested case, the largest ratio observed being $A(679)/679^2\approx7.65$.

Footnotes and References

The references below supply the background facts used above. The reduction of the focus condition to two symmetric equations, the exclusion of $m=1$ and of the norms $9$ and $117$, the descent behind the three chains, the counts of $2$ and $4$ triangles per pair and the closed formula for $A(n)$ are derived explicitly here.

Project Euler 385 — Ellipses Inside Triangles. The official statement notes that every triangle contains a unique ellipse of largest area, defines A(n) as the sum of the areas of the triangles with integer vertices of absolute value at most n whose largest inscribed ellipse has foci (±√13, 0), lists the two triangles behind A(8) = 72, and asks for A(10⁹) after giving A(10), A(100) and A(1000).
Wikipedia — Steiner inellipse. The article defines it as the unique inscribed ellipse tangent to the sides at their midpoints, states that it has the greatest area of all inellipses, that its area is π/(3√3) times the triangle's and that its centre is the centroid, and reduces its proofs to an equilateral triangle — the argument of Section 1.
Wikipedia — Marden's theorem. For a cubic polynomial p with non-collinear zeros, the foci of the Steiner inellipse of the triangle they form are the zeros of p′. The article credits the discovery to Jörg Siebeck in 1864 and reports that Dan Kalman titled his paper "Marden's theorem" after first reading the result in Marden's book. Section 4 uses it as a second derivation.
OEIS A034017 — numbers primitively represented by x² + xy + y². Its positive terms are the numbers represented with gcd(x, y) = 1: 1, 3, 7, 13, 19, 21, 31, 37, 39, … (the listed initial 0 is not such a representation). Its formula states that each is a product of primes of the form 3a + 1 times at most one factor 3, so none is divisible by 9 — the general form of the congruence proved in Section 7.
Wikipedia — Pell's equation. The article treats x² − ny² = 1, states that all solutions follow from the powers of the fundamental solution, and lists (2, 1) for n = 3. For x² − ny² = N it notes that the solutions fall into finitely many classes, each arising from a class representative times a solution of x² − ny² = 1 — the structure of Section 8.
OEIS A077236 and OEIS A077235. Both satisfy a(n) = 4a(n−1) − a(n−2), the first starting 4, 11, 40, 149, 556 and the second 5, 16, 59, 220, 821, and both state that a(n)² − 3b(n)² = 13 for a companion b. These are the two chains of d-values for the norm 3, with k = 3b(n).
OEIS A001075. Defined by a(0) = 1, a(1) = 2, a(n) = 4a(n−1) − a(n−2), it begins 1, 2, 7, 26, 97, 362, 1351, and the entry states that these are the x-values of x² − 3y² = 1, the y-values being A001353. This is the chain of d-values for the norm 39, its first term the degenerate pair discarded in Section 8.
Wikipedia — Shoelace formula. For three vertices the determinant form gives the signed area as half of (x₁y₂ − x₂y₁) + (x₂y₃ − x₃y₂) + (x₃y₁ − x₁y₃); Section 10 uses the equivalent arrangement x₁(y₂ − y₃) + x₂(y₃ − y₁) + x₃(y₁ − y₂).
C++ — Euler385.cpp. The linked immutable revision contains build_directions, generate_pell_pairs, build_triangle_key, solve_fast, solve_bruteforce_small and run_checkpoints, and is the implementation record for this article; the independent checks of the verification section are not part of it.

Problem 385 source code

C++

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr i64 kDefaultN = 1'000'000'000LL;

struct Options {
    i64 n = kDefaultN;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

struct Point {
    i64 x = 0;
    i64 y = 0;

    bool operator<(const Point& other) const {
        if (x != other.x) {
            return x < other.x;
        }
        return y < other.y;
    }

    bool operator==(const Point& other) const {
        return x == other.x && y == other.y;
    }
};

struct TriangleKey {
    std::array<Point, 3> vertices{};

    bool operator==(const TriangleKey& other) const {
        return vertices == other.vertices;
    }
};

struct TriangleKeyHash {
    std::size_t operator()(const TriangleKey& key) const {
        std::size_t h = 0ULL;
        for (const Point& p : key.vertices) {
            const std::size_t hx = std::hash<i64>{}(p.x);
            const std::size_t hy = std::hash<i64>{}(p.y);
            h ^= hx + 0x9e3779b97f4a7c15ULL + (h << 6U) + (h >> 2U);
            h ^= hy + 0x9e3779b97f4a7c15ULL + (h << 6U) + (h >> 2U);
        }
        return h;
    }
};

struct Direction {
    int q0 = 0;
    int u = 0;
    int v = 0;
    int m = 3;
    int pell_n = 0;  // k^2 - 3 d^2 = pell_n.
    i64 d_max = 0;
    i64 k_max = 0;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
            if (parsed_u64 > static_cast<u64>(std::numeric_limits<i64>::max())) {
                std::cerr << "--n is too large for this implementation.\n";
                return false;
            }
            options.n = static_cast<i64>(parsed_u64);
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    constexpr std::size_t kMinUnitsForParallel = 8ULL;

    if (!allow_multithreading || workload_units < 2ULL || workload_units < kMinUnitsForParallel) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

i64 abs_i64(const i64 x) {
    return x >= 0 ? x : -x;
}

u64 abs_to_u64(const i64 x) {
    return static_cast<u64>(x >= 0 ? x : -x);
}

std::string to_string_u128(u128 value) {
    if (value == 0U) {
        return "0";
    }

    std::string digits;
    while (value > 0U) {
        const unsigned digit = static_cast<unsigned>(value % 10U);
        digits.push_back(static_cast<char>('0' + digit));
        value /= 10U;
    }
    std::reverse(digits.begin(), digits.end());
    return digits;
}

std::vector<std::pair<i64, i64>> generate_pell_pairs(const int pell_n,
                                                      const i64 d_max,
                                                      const i64 k_max) {
    // Problem-specific seeds for k^2 - 3 d^2 = -39 and = -3.
    std::vector<std::pair<i64, i64>> seeds;
    if (pell_n == -39) {
        seeds = {{3, 4}, {6, 5}};
    } else if (pell_n == -3) {
        seeds = {{3, 2}};
    } else {
        return {};
    }

    std::vector<std::pair<i64, i64>> all;
    for (const auto [k0, d0] : seeds) {
        i64 k = k0;
        i64 d = d0;

        while (d <= d_max && k <= k_max) {
            all.emplace_back(d, k);

            if (k > (std::numeric_limits<i64>::max() - 3LL * d) / 2LL ||
                d > (std::numeric_limits<i64>::max() - k) / 2LL) {
                break;
            }

            const i64 next_k = 2LL * k + 3LL * d;
            const i64 next_d = k + 2LL * d;
            k = next_k;
            d = next_d;
        }
    }

    std::sort(all.begin(), all.end());
    all.erase(std::unique(all.begin(), all.end()), all.end());
    return all;
}

std::vector<Direction> build_directions(const i64 n) {
    std::vector<Direction> directions;

    for (const int q0 : {3, 39}) {
        const int pell_n = -117 / q0;  // Derived from q0(3 d^2 - k^2) = 117.
        const int limit = static_cast<int>(std::sqrt(static_cast<long double>(q0))) + 2;

        for (int u = -limit; u <= limit; ++u) {
            for (int v = -limit; v <= limit; ++v) {
                if (u == 0 && v == 0) {
                    continue;
                }
                if (std::gcd(abs_i64(u), abs_i64(v)) != 1) {
                    continue;
                }
                if (u * u + u * v + v * v != q0) {
                    continue;
                }

                const int a = 2 * u + v;
                const int b = u + 2 * v;
                const int m = static_cast<int>(std::gcd(abs_i64(a), abs_i64(b)));
                if (m != 3) {
                    continue;
                }

                const i64 x_scale = std::max({abs_i64(u), abs_i64(v), abs_i64(u + v)});
                const i64 y_scale_raw = std::max({abs_i64(a), abs_i64(b), abs_i64(u - v)});
                if (x_scale == 0 || y_scale_raw % m != 0) {
                    continue;
                }
                const i64 y_scale = y_scale_raw / m;
                if (y_scale == 0) {
                    continue;
                }

                Direction dir;
                dir.q0 = q0;
                dir.u = u;
                dir.v = v;
                dir.m = m;
                dir.pell_n = pell_n;
                dir.d_max = n / x_scale;
                dir.k_max = n / y_scale;
                if (dir.d_max > 0 && dir.k_max > 0) {
                    directions.push_back(dir);
                }
            }
        }
    }

    std::sort(directions.begin(),
              directions.end(),
              [](const Direction& lhs, const Direction& rhs) {
                  if (lhs.q0 != rhs.q0) {
                      return lhs.q0 < rhs.q0;
                  }
                  if (lhs.u != rhs.u) {
                      return lhs.u < rhs.u;
                  }
                  return lhs.v < rhs.v;
              });
    return directions;
}

bool build_triangle_key(const Direction& dir,
                        const i64 d,
                        const i64 k,
                        const i64 n,
                        TriangleKey& out_key,
                        u64& out_area) {
    const i64 u = dir.u;
    const i64 v = dir.v;
    const i64 m = dir.m;

    const i64 a = 2LL * u + v;
    const i64 b = u + 2LL * v;
    if (a % m != 0LL || b % m != 0LL || (u - v) % m != 0LL) {
        return false;
    }

    const i64 x1 = d * u;
    const i64 x2 = d * v;
    const i64 x3 = -d * (u + v);

    const i64 y1 = (b / m) * k;
    const i64 y2 = -(a / m) * k;
    const i64 y3 = ((u - v) / m) * k;

    if (abs_i64(x1) > n || abs_i64(x2) > n || abs_i64(x3) > n || abs_i64(y1) > n ||
        abs_i64(y2) > n || abs_i64(y3) > n) {
        return false;
    }

    std::array<Point, 3> points = {{{x1, y1}, {x2, y2}, {x3, y3}}};
    std::sort(points.begin(), points.end());

    if (points[0] == points[1] || points[1] == points[2]) {
        return false;
    }

    out_key.vertices = points;
    out_area = static_cast<u64>(dir.q0) * abs_to_u64(d * k);
    return true;
}

u128 sum_triangle_areas_from_map(const std::unordered_map<TriangleKey, u64, TriangleKeyHash>& map) {
    u128 total = 0U;
    for (const auto& entry : map) {
        total += static_cast<u128>(entry.second);
    }
    return total;
}

u128 solve_fast(const i64 n,
                const bool allow_multithreading,
                const unsigned requested_threads,
                std::size_t* out_triangle_count = nullptr) {
    const std::vector<Direction> directions = build_directions(n);
    if (directions.empty()) {
        if (out_triangle_count != nullptr) {
            *out_triangle_count = 0ULL;
        }
        return 0U;
    }

    i64 max_d_n39 = 0;
    i64 max_k_n39 = 0;
    i64 max_d_n3 = 0;
    i64 max_k_n3 = 0;

    for (const Direction& dir : directions) {
        if (dir.pell_n == -39) {
            max_d_n39 = std::max(max_d_n39, dir.d_max);
            max_k_n39 = std::max(max_k_n39, dir.k_max);
        } else if (dir.pell_n == -3) {
            max_d_n3 = std::max(max_d_n3, dir.d_max);
            max_k_n3 = std::max(max_k_n3, dir.k_max);
        }
    }

    const std::vector<std::pair<i64, i64>> pairs_n39 = generate_pell_pairs(-39, max_d_n39, max_k_n39);
    const std::vector<std::pair<i64, i64>> pairs_n3 = generate_pell_pairs(-3, max_d_n3, max_k_n3);

    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, directions.size());
    std::vector<std::unordered_map<TriangleKey, u64, TriangleKeyHash>> partial_maps(
        static_cast<std::size_t>(threads));
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    for (unsigned t = 0U; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            auto& local = partial_maps[static_cast<std::size_t>(t)];

            for (std::size_t i = static_cast<std::size_t>(t); i < directions.size();
                 i += static_cast<std::size_t>(threads)) {
                const Direction& dir = directions[i];
                const auto& pairs = (dir.pell_n == -39) ? pairs_n39 : pairs_n3;

                for (const auto [d, k_abs] : pairs) {
                    if (d > dir.d_max || k_abs > dir.k_max) {
                        continue;
                    }

                    for (const i64 sign : {-1LL, 1LL}) {
                        TriangleKey key;
                        u64 area = 0ULL;
                        if (!build_triangle_key(dir, d, sign * k_abs, n, key, area)) {
                            continue;
                        }

                        const auto it = local.find(key);
                        if (it == local.end()) {
                            local.emplace(key, area);
                        } else if (it->second != area) {
                            std::cerr << "Internal inconsistency: same triangle got two areas.\n";
                            std::abort();
                        }
                    }
                }
            }
        });
    }

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

    std::unordered_map<TriangleKey, u64, TriangleKeyHash> merged;
    for (auto& local : partial_maps) {
        for (auto& kv : local) {
            const auto it = merged.find(kv.first);
            if (it == merged.end()) {
                merged.emplace(std::move(kv));
            } else if (it->second != kv.second) {
                std::cerr << "Internal inconsistency during merge.\n";
                std::abort();
            }
        }
    }

    if (out_triangle_count != nullptr) {
        *out_triangle_count = merged.size();
    }
    return sum_triangle_areas_from_map(merged);
}

u128 solve_bruteforce_small(const int n) {
    std::unordered_map<TriangleKey, u64, TriangleKeyHash> triangles;

    for (int x1 = -n; x1 <= n; ++x1) {
        for (int y1 = -n; y1 <= n; ++y1) {
            for (int x2 = -n; x2 <= n; ++x2) {
                const i64 coeff = static_cast<i64>(x1) + 2LL * static_cast<i64>(x2);
                const i64 rhs = -(2LL * static_cast<i64>(x1) + static_cast<i64>(x2)) *
                                static_cast<i64>(y1);

                const auto check_candidate = [&](const i64 y2) {
                    if (y2 < -n || y2 > n) {
                        return;
                    }

                    const i64 re = static_cast<i64>(x1) * static_cast<i64>(x1) -
                                   static_cast<i64>(y1) * static_cast<i64>(y1) +
                                   static_cast<i64>(x1) * static_cast<i64>(x2) -
                                   static_cast<i64>(y1) * y2 +
                                   static_cast<i64>(x2) * static_cast<i64>(x2) - y2 * y2;
                    if (re != 39LL) {
                        return;
                    }

                    const i64 x3 = -static_cast<i64>(x1) - static_cast<i64>(x2);
                    const i64 y3 = -static_cast<i64>(y1) - y2;
                    if (abs_i64(x3) > n || abs_i64(y3) > n) {
                        return;
                    }

                    std::array<Point, 3> points = {{{x1, y1}, {x2, y2}, {x3, y3}}};
                    std::sort(points.begin(), points.end());
                    if (points[0] == points[1] || points[1] == points[2]) {
                        return;
                    }

                    const i64 area2 = abs_i64((points[1].x - points[0].x) * (points[2].y - points[0].y) -
                                              (points[1].y - points[0].y) * (points[2].x - points[0].x));
                    if (area2 == 0 || (area2 & 1LL) != 0LL) {
                        return;
                    }

                    TriangleKey key;
                    key.vertices = points;
                    triangles[key] = static_cast<u64>(area2 / 2LL);
                };

                if (coeff == 0LL) {
                    if (rhs != 0LL) {
                        continue;
                    }
                    for (int y2 = -n; y2 <= n; ++y2) {
                        check_candidate(static_cast<i64>(y2));
                    }
                } else {
                    if (rhs % coeff != 0LL) {
                        continue;
                    }
                    check_candidate(rhs / coeff);
                }
            }
        }
    }

    return sum_triangle_areas_from_map(triangles);
}

bool run_checkpoints(const bool allow_multithreading, const unsigned requested_threads) {
    struct Checkpoint {
        int n;
        u64 expected;
    };

    const std::vector<Checkpoint> known = {
        {8, 72ULL},
        {10, 252ULL},
        {100, 34'632ULL},
        {1'000, 3'529'008ULL},
    };

    for (const Checkpoint& cp : known) {
        const u128 got = solve_fast(cp.n, allow_multithreading, requested_threads);
        if (got != static_cast<u128>(cp.expected)) {
            std::cerr << "Checkpoint failed for A(" << cp.n << "): expected " << cp.expected
                      << ", got " << to_string_u128(got) << '\n';
            return false;
        }
    }

    const u128 brute20 = solve_bruteforce_small(20);
    const u128 fast20 = solve_fast(20, false, 1U);
    if (brute20 != fast20) {
        std::cerr << "Bruteforce cross-check failed for n=20: brute=" << to_string_u128(brute20)
                  << ", fast=" << to_string_u128(fast20) << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.n < 0) {
        std::cerr << "--n must be non-negative.\n";
        return 1;
    }

    if (options.run_checkpoints &&
        !run_checkpoints(options.allow_multithreading, options.requested_threads)) {
        return 1;
    }

    std::size_t triangle_count = 0ULL;
    const u128 answer = solve_fast(options.n,
                                   options.allow_multithreading,
                                   options.requested_threads,
                                   &triangle_count);

    std::cout << to_string_u128(answer) << '\n';
    std::cerr << "Triangles counted: " << triangle_count << '\n';
    return 0;
}

Python

import math

def generate_pell_pairs(pell_n, d_max, k_max):
    if pell_n == -39:
        seeds = [(3, 4), (6, 5)]
    elif pell_n == -3:
        seeds = [(3, 2)]
    else:
        return []
        
    all_pairs = []
    for k0, d0 in seeds:
        k, d = k0, d0
        while d <= d_max and k <= k_max:
            all_pairs.append((d, k))
            next_k = 2 * k + 3 * d
            next_d = k + 2 * d
            k, d = next_k, next_d
    return sorted(list(set(all_pairs)))

def build_directions(n):
    directions = []
    for q0 in (3, 39):
        pell_n = -117 // q0
        limit = int(math.sqrt(q0)) + 2
        for u in range(-limit, limit + 1):
            for v in range(-limit, limit + 1):
                if u == 0 and v == 0: continue
                if math.gcd(abs(u), abs(v)) != 1: continue
                if u * u + u * v + v * v != q0: continue
                
                a = 2 * u + v
                b = u + 2 * v
                m = math.gcd(abs(a), abs(b))
                if m != 3: continue
                
                x_scale = max(abs(u), abs(v), abs(u + v))
                y_scale_raw = max(abs(a), abs(b), abs(u - v))
                if x_scale == 0 or y_scale_raw % m != 0: continue
                
                y_scale = y_scale_raw // m
                if y_scale == 0: continue
                
                d_max = n // x_scale
                k_max = n // y_scale
                
                if d_max > 0 and k_max > 0:
                    directions.append({
                        'q0': q0, 'u': u, 'v': v, 'm': m,
                        'pell_n': pell_n, 'd_max': d_max, 'k_max': k_max
                    })
                    
    directions.sort(key=lambda d: (d['q0'], d['u'], d['v']))
    return directions

def build_triangle_key(dir_obj, d, k, n):
    u, v, m = dir_obj['u'], dir_obj['v'], dir_obj['m']
    a = 2 * u + v
    b = u + 2 * v
    if a % m != 0 or b % m != 0 or (u - v) % m != 0:
        return None
        
    x1 = d * u
    x2 = d * v
    x3 = -d * (u + v)
    
    y1 = (b // m) * k
    y2 = -(a // m) * k
    y3 = ((u - v) // m) * k
    
    if max(abs(x1), abs(x2), abs(x3), abs(y1), abs(y2), abs(y3)) > n:
        return None
        
    points = sorted([(x1, y1), (x2, y2), (x3, y3)])
    if points[0] == points[1] or points[1] == points[2]:
        return None
        
    area = dir_obj['q0'] * abs(d * k)
    return tuple(points), area

def solve():
    n = 1000000000
    directions = build_directions(n)
    if not directions: return "0"
    
    max_d_n39 = max_k_n39 = max_d_n3 = max_k_n3 = 0
    for d in directions:
        if d['pell_n'] == -39:
            max_d_n39 = max(max_d_n39, d['d_max'])
            max_k_n39 = max(max_k_n39, d['k_max'])
        elif d['pell_n'] == -3:
            max_d_n3 = max(max_d_n3, d['d_max'])
            max_k_n3 = max(max_k_n3, d['k_max'])
            
    pairs_n39 = generate_pell_pairs(-39, max_d_n39, max_k_n39)
    pairs_n3 = generate_pell_pairs(-3, max_d_n3, max_k_n3)
    
    triangles = {}
    for d_obj in directions:
        pairs = pairs_n39 if d_obj['pell_n'] == -39 else pairs_n3
        for d, k_abs in pairs:
            if d > d_obj['d_max'] or k_abs > d_obj['k_max']:
                continue
            for sign in (-1, 1):
                res = build_triangle_key(d_obj, d, sign * k_abs, n)
                if not res: continue
                key, area = res
                if key not in triangles:
                    triangles[key] = area
                elif triangles[key] != area:
                    raise Exception("Inconsistency: same triangle got two areas.")
                    
    total_area = sum(triangles.values())
    return str(total_area)

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

Java

import java.util.*;

public class Euler385 {
    static class Point implements Comparable<Point> {
        long x, y;

        Point(long x, long y) {
            this.x = x;
            this.y = y;
        }

        @Override
        public int compareTo(Point o) {
            if (this.x != o.x)
                return Long.compare(this.x, o.x);
            return Long.compare(this.y, o.y);
        }

        @Override
        public boolean equals(Object o) {
            if (!(o instanceof Point))
                return false;
            Point p = (Point) o;
            return this.x == p.x && this.y == p.y;
        }

        @Override
        public int hashCode() {
            int hx = Long.hashCode(x), hy = Long.hashCode(y);
            return hx ^ (hy + 0x9e3779b9 + (hx << 6) + (hx >> 2));
        }
    }

    static class TriangleKey {
        Point[] pts;

        TriangleKey(Point p1, Point p2, Point p3) {
            pts = new Point[] { p1, p2, p3 };
            Arrays.sort(pts);
        }

        @Override
        public boolean equals(Object o) {
            if (!(o instanceof TriangleKey))
                return false;
            TriangleKey tk = (TriangleKey) o;
            return pts[0].equals(tk.pts[0]) && pts[1].equals(tk.pts[1]) && pts[2].equals(tk.pts[2]);
        }

        @Override
        public int hashCode() {
            int h = 0;
            for (Point p : pts) {
                h ^= p.hashCode() + 0x9e3779b9 + (h << 6) + (h >> 2);
            }
            return h;
        }
    }

    static class Pair implements Comparable<Pair> {
        long d, k;

        Pair(long d, long k) {
            this.d = d;
            this.k = k;
        }

        @Override
        public int compareTo(Pair o) {
            if (this.d != o.d)
                return Long.compare(this.d, o.d);
            return Long.compare(this.k, o.k);
        }

        @Override
        public boolean equals(Object o) {
            if (!(o instanceof Pair))
                return false;
            Pair p = (Pair) o;
            return this.d == p.d && this.k == p.k;
        }

        @Override
        public int hashCode() {
            return Long.hashCode(d) ^ Long.hashCode(k);
        }
    }

    static class Direction {
        int q0, u, v, m, pell_n;
        long d_max, k_max;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = a % b;
            a = b;
            b = t;
        }
        return Math.abs(a);
    }

    static List<Pair> generatePellPairs(int pellN, long dMax, long kMax) {
        List<Pair> seeds = new ArrayList<>();
        if (pellN == -39) {
            seeds.add(new Pair(4, 3)); // (d, k) form => k0=3, d0=4
            seeds.add(new Pair(5, 6)); // k0=6, d0=5
        } else if (pellN == -3) {
            seeds.add(new Pair(2, 3)); // k0=3, d0=2
        } else {
            return new ArrayList<>();
        }

        List<Pair> all = new ArrayList<>();
        for (Pair s : seeds) {
            long k = s.k, d = s.d;
            while (d <= dMax && k <= kMax) {
                all.add(new Pair(d, k));
                if (k > (Long.MAX_VALUE - 3 * d) / 2 || d > (Long.MAX_VALUE - k) / 2)
                    break;
                long nK = 2 * k + 3 * d;
                long nD = k + 2 * d;
                k = nK;
                d = nD;
            }
        }

        Set<Pair> set = new HashSet<>(all);
        all = new ArrayList<>(set);
        Collections.sort(all);
        return all;
    }

    static List<Direction> buildDirections(long n) {
        List<Direction> dirs = new ArrayList<>();
        int[] q0s = { 3, 39 };
        for (int q0 : q0s) {
            int pellN = -117 / q0;
            int limit = (int) Math.sqrt(q0) + 2;
            for (int u = -limit; u <= limit; u++) {
                for (int v = -limit; v <= limit; v++) {
                    if (u == 0 && v == 0)
                        continue;
                    if (gcd(u, v) != 1)
                        continue;
                    if (u * u + u * v + v * v != q0)
                        continue;

                    int a = 2 * u + v;
                    int b = u + 2 * v;
                    int m = (int) gcd(a, b);
                    if (m != 3)
                        continue;

                    long xScale = Math.max(Math.max(Math.abs(u), Math.abs(v)), Math.abs(u + v));
                    long yScaleRaw = Math.max(Math.max(Math.abs(a), Math.abs(b)), Math.abs(u - v));
                    if (xScale == 0 || yScaleRaw % m != 0)
                        continue;

                    long yScale = yScaleRaw / m;
                    if (yScale == 0)
                        continue;

                    long dMax = n / xScale;
                    long kMax = n / yScale;

                    if (dMax > 0 && kMax > 0) {
                        Direction d = new Direction();
                        d.q0 = q0;
                        d.u = u;
                        d.v = v;
                        d.m = m;
                        d.pell_n = pellN;
                        d.d_max = dMax;
                        d.k_max = kMax;
                        dirs.add(d);
                    }
                }
            }
        }
        return dirs;
    }

    static class Result {
        TriangleKey tk;
        long area;
    }

    static Result buildTriangleKey(Direction dir, long d, long k, long n) {
        long u = dir.u, v = dir.v, m = dir.m;
        long a = 2 * u + v;
        long b = u + 2 * v;

        if (a % m != 0 || b % m != 0 || (u - v) % m != 0)
            return null;

        long x1 = d * u;
        long x2 = d * v;
        long x3 = -d * (u + v);

        long y1 = (b / m) * k;
        long y2 = -(a / m) * k;
        long y3 = ((u - v) / m) * k;

        if (Math.abs(x1) > n || Math.abs(x2) > n || Math.abs(x3) > n ||
                Math.abs(y1) > n || Math.abs(y2) > n || Math.abs(y3) > n) {
            return null;
        }

        Point p1 = new Point(x1, y1), p2 = new Point(x2, y2), p3 = new Point(x3, y3);
        TriangleKey tk = new TriangleKey(p1, p2, p3);

        if (tk.pts[0].equals(tk.pts[1]) || tk.pts[1].equals(tk.pts[2]))
            return null;

        Result r = new Result();
        r.tk = tk;
        r.area = dir.q0 * Math.abs(d * k);
        return r;
    }

    static String solve() {
        long n = 1000000000L;
        List<Direction> dirs = buildDirections(n);
        if (dirs.isEmpty())
            return "0";

        long maxDN39 = 0, maxKN39 = 0, maxDN3 = 0, maxKN3 = 0;
        for (Direction d : dirs) {
            if (d.pell_n == -39) {
                maxDN39 = Math.max(maxDN39, d.d_max);
                maxKN39 = Math.max(maxKN39, d.k_max);
            } else if (d.pell_n == -3) {
                maxDN3 = Math.max(maxDN3, d.d_max);
                maxKN3 = Math.max(maxKN3, d.k_max);
            }
        }

        List<Pair> pairsN39 = generatePellPairs(-39, maxDN39, maxKN39);
        List<Pair> pairsN3 = generatePellPairs(-3, maxDN3, maxKN3);

        Map<TriangleKey, Long> map = new HashMap<>();

        for (Direction dir : dirs) {
            List<Pair> pairs = (dir.pell_n == -39) ? pairsN39 : pairsN3;
            for (Pair p : pairs) {
                if (p.d > dir.d_max || p.k > dir.k_max)
                    continue;
                for (long sign : new long[] { -1, 1 }) {
                    Result res = buildTriangleKey(dir, p.d, sign * p.k, n);
                    if (res == null)
                        continue;

                    if (map.containsKey(res.tk) && map.get(res.tk) != res.area) {
                        throw new RuntimeException("Inconsistency in triangle area");
                    }
                    map.put(res.tk, res.area);
                }
            }
        }

        long total = 0;
        for (long a : map.values())
            total += a;
        return Long.toString(total);
    }

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