Problem 246: Tangents to an Ellipse

Project Euler Problem 246 Solution

EulerSolve provides an optimized solution for Project Euler Problem 246, Tangents to an Ellipse, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The ellipse in the problem is the locus of points $Q$ such that $QM+QG=15000$, where the foci are $M=(-2000,1500)$ and $G=(8000,1500)$. For an integer lattice point $P$ outside that ellipse, two tangents can be drawn from $P$ to the ellipse. We must count the lattice points for which the angle between those two tangent rays is greater than $45^\circ$. A naive geometric search would be awkward because every candidate point would seem to require solving for the tangency points. The implementations avoid that completely. They translate the ellipse to its center, derive a closed inequality for the tangent angle, and then count lattice points column by column. Mathematical Approach Let $(X,Y)$ be the original coordinates, and shift to centered coordinates $$x=X-3000,\qquad y=Y-1500.$$ The midpoint of the two foci is $(3000,1500)$, so after translation the foci become $(\pm 5000,0)$. The lattice is unchanged because the shift vector is integral. Standard Form of the Ellipse The focal distance is $2c=10000$, so $c=5000$. The constant distance sum is $2a=15000$, so $a=7500$. Therefore $$b^2=a^2-c^2=7500^2-5000^2=31\,250\,000.$$ The ellipse is thus $$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1,\qquad a^2=56\,250\,000,\qquad b^2=31\,250\,000.$$ All later formulas depend only on these centered quantities....

Detailed mathematical approach

Problem Summary

The circle $c$ has centre $M=(-2000,1500)$ and radius $15000$, and the point $G=(8000,1500)$ lies inside it. The points that are as far from $c$ as from $G$ form an ellipse $e$. From a point $P$ outside $e$ two tangents can be drawn to $e$, and they touch it at $R$ and $S$. We must count the lattice points $P$, the points with two integer coordinates, for which the angle $\angle RPS$ is greater than $45^\circ$. The statement gives no values to check against. ^[1]

The whole solution in one sentence: centre the ellipse, turn the angle condition into exact integer inequalities with the help of the circle on which the two tangents are perpendicular, count the qualifying lattice points of one quadrant column by column, each column being one interval whose ends follow from integer square roots, and use symmetry for the rest of the plane.

Mathematical Approach

The ellipse shrunk by the factor $2500$ serves as a running example whose numbers can be checked by hand.

1. The ellipse

The distance from a point $Q$ to the circle $c$ is $|QM-15000|$, measured along the line through $M$ and $Q$. A point $Q$ on or outside $c$ cannot be as far from $c$ as from $G$: that would need $QM-QG=15000$, but the triangle inequality gives $QM-QG\le MG=10000$. For $Q$ inside $c$ the condition reads $15000-QM=QG$, that is

$$QM+QG=15000.$$

So $e$ is the ellipse with foci $M$ and $G$ whose distances to the two foci add up to $2a=15000$, and $c$ is its circular directrix. ^[2] Its centre is the midpoint $(3000,1500)$ of $MG$, its semi-major axis is $a=7500$, the distance from the centre to a focus is $5000$, and the semi-minor axis $b$ satisfies $b^2=a^2-5000^2=31\,250\,000$. ^[2] From now on we use the centred coordinates $x=X-3000$ and $y=Y-1500$, in which lattice points stay lattice points and the foci are $(\pm f,0)$ with $f=5000$. The distances $d_-$ and $d_+$ of a point to the two foci satisfy $d_-^2-d_+^2=4fx$, so on the ellipse, where $d_-+d_+=2a$, we get $d_--d_+=2fx/a$ and $d_-=a+fx/a$. Squaring this and using $b^2=a^2-f^2$ gives the equation $D(x,y)=0$ with

$$D(x,y)=b^2x^2+a^2y^2-a^2b^2,$$

and since $a\pm fx/a\gt0$ when $|x|\le a$, the steps can be reversed. ^[2] Along every ray from the centre $D$ increases, so $P$ lies outside the ellipse exactly when $D(P)\gt0$. The vertex $(7500,0)$, the point $(10500,1500)$ of the statement, has $D=0$. The running example is the same ellipse shrunk by $2500$: $a=3$, $b^2=5$ and $D=5x^2+9y^2-45$.

2. The two tangents from a point

A line $y=mx+n$ meets the ellipse where

$$(b^2+a^2m^2)\,x^2+2a^2mn\,x+a^2(n^2-b^2)=0.$$

It touches the ellipse when this quadratic has a double root; its discriminant is $4a^2b^2(a^2m^2+b^2-n^2)$, so the condition is $n^2=a^2m^2+b^2$. ^[2] The line passes through $P=(x_0,y_0)$ when $n=y_0-mx_0$, and inserting this gives a quadratic equation for the slopes of the tangents through $P$:

$$(x_0^2-a^2)\,m^2-2x_0y_0\,m+(y_0^2-b^2)=0.$$

Its discriminant is

$$4x_0^2y_0^2-4(x_0^2-a^2)(y_0^2-b^2)=4D(x_0,y_0).$$

So two tangents exist exactly when $P$ lies outside the ellipse. If $x_0^2=a^2$, the equation has degree one and the second tangent is the vertical line $x=x_0$.

In the running example take $P=(0,3)$. The equation is $-9m^2+4=0$, so $m=\pm\frac23$, and the tangents $y=3\pm\frac23x$ touch the ellipse at $R=(-2,\frac53)$ and $S=(2,\frac53)$; indeed $\frac49+\frac{25}{45}=1$. From $\vec{PR}=(-2,-\frac43)$ and $\vec{PS}=(2,-\frac43)$ the dot product gives $\cos\angle RPS=-\frac{5}{13}$, so $\angle RPS\approx112.62^\circ$.

3. A formula for the angle

Let the tangents make the angles $\alpha_1$ and $\alpha_2$ with the $x$-axis, so that their slopes are $m_1=\tan\alpha_1$ and $m_2=\tan\alpha_2$. The angle $\angle RPS$ is $\alpha_1-\alpha_2$ or its supplement, up to sign. Since $\tan(\alpha_1-\alpha_2)=\frac{m_1-m_2}{1+m_1m_2}$ and $\tan(180^\circ-\theta)=-\tan\theta$, ^[3] both possibilities have the same squared tangent:

$$\begin{aligned}\tan^2\angle RPS&=\frac{(m_1-m_2)^2}{(1+m_1m_2)^2}\\&=\frac{(m_1+m_2)^2-4m_1m_2}{(1+m_1m_2)^2}.\end{aligned}$$

By Vieta's formulas ^[4] $m_1+m_2=\frac{2x_0y_0}{x_0^2-a^2}$ and $m_1m_2=\frac{y_0^2-b^2}{x_0^2-a^2}$. By the same computation as for the discriminant above, the numerator becomes $4D(x_0,y_0)/(x_0^2-a^2)^2$, and $1+m_1m_2=s(x_0,y_0)/(x_0^2-a^2)$ with

$$s(x,y)=x^2+y^2-a^2-b^2.$$

When $s\ne0$, the factors $x_0^2-a^2$ cancel, and

$$\tan^2\angle RPS=\frac{4D(x_0,y_0)}{s(x_0,y_0)^2};$$

the case $s=0$ is the right angle of Section 4. If one tangent is vertical, the other has slope $m_2=\frac{y_0^2-b^2}{2x_0y_0}$; for $m_2\ne0$ the angle between them has the squared tangent $1/m_2^2=\frac{4a^2y_0^2}{(y_0^2-b^2)^2}$, which is the formula again, because $D=a^2y_0^2$ and $s=y_0^2-b^2$ when $x_0^2=a^2$. For $P=(0,3)$ in the running example $D=36$ and $s=-5$, so $\tan^2\angle RPS=\frac{144}{25}$, in agreement with $\cos\angle RPS=-\frac5{13}$ and $\sin\angle RPS=\frac{12}{13}$. This is also how the isoptic curves of the ellipse are found: the points that see it under a fixed angle $\alpha$ lie on the quartic $s^2\tan^2\alpha=4D$, together with those that see it under $180^\circ-\alpha$. ^[5]

4. Obtuse or acute: the director circle

The formula gives only the square of the tangent, so it cannot tell $\angle RPS$ from $180^\circ-\angle RPS$; the sign of $s$ can. The tangents are perpendicular exactly when $1+m_1m_2=0$, that is $s=0$; with a vertical tangent, when the other one is horizontal, $y_0^2=b^2$, which again means $s=0$. So $\angle RPS=90^\circ$ exactly on the circle $x^2+y^2=a^2+b^2$, the director circle of the ellipse, ^[6] here of radius $\sqrt{87\,500\,000}\approx9354.1$. Every point of the ellipse is at most $a$ away from the centre, so the director circle lies outside the ellipse, and the points outside the ellipse and inside the director circle form a ring.

The angle $\angle RPS$ changes continuously with $P$, because the points of contact do, and equals $90^\circ$ only on the director circle. The ring is connected: it contains the circle of radius $7501$ about the centre, and the ray from the centre joins every point of the ring to that circle inside the ring. By the intermediate value theorem ^[7] the angle therefore stays on one side of $90^\circ$ in the whole ring, the obtuse side: scaling by $2500$ preserves angles, so the ring point $(0,7500)$ sees the problem's ellipse under the same $112.62^\circ$ as $(0,3)$ sees the running example. The outside of the director circle is connected as well, and from distance $2a=15000$ the ellipse, which fits into the disc of radius $a$, is seen under at most $2\arcsin\frac12=60^\circ$. Hence $\angle RPS\gt90^\circ$ in the ring and $\angle RPS\lt90^\circ$ outside the director circle.

5. The test for a single point

Put $F=4D-s^2$. By Sections 3 and 4 a lattice point $P$ qualifies exactly when $D(P)\gt0$ and either $s(P)\le0$ or $F(P)\gt0$. Indeed, inside or on the director circle, where $s\le0$, the angle is at least $90^\circ$. Outside it, where $s\gt0$, the angle lies between $0^\circ$ and $90^\circ$, where the tangent increases, so $\angle RPS\gt45^\circ$ means $\tan^2\angle RPS\gt1$, that is $4D\gt s^2$; this already forces $D\gt0$.

The case split cannot be skipped. In the ring, $F\gt0$ means $\angle RPS\lt135^\circ$, so the test $F\gt0$ alone would reject the points closest to the ellipse. On the positive $x$-axis of the problem's ellipse, $D=b^2z$ and $s=z-b^2$ with $z=x^2-a^2$, so

$$F(x,0)=4b^2z-(z-b^2)^2=-\bigl(z^2-6b^2z+b^4\bigr),$$

which is positive exactly when

$$(3-2\sqrt2)\,b^2\lt z\lt(3+2\sqrt2)\,b^2,$$

that is for $7850\le x\le15439$. The half-axis therefore splits into four zones:

$7501\le x\le7849$: ring, $F\lt0$, from $177.49^\circ$ at $x=7501$ to $135.02^\circ$ at $x=7849$: counted;
$7850\le x\le9354$: ring, between $90^\circ$ and $135^\circ$: counted;
$9355\le x\le15439$: outside the director circle, $F\gt0$, down to $45.0029^\circ$ at $x=15439$: counted;
$x\ge15440$: $F\lt0$, from $44.9994^\circ$ at $x=15440$ on: not counted.

So the positive $x$-axis contains $15439-7501+1=7939$ qualifying lattice points, and the test $F\gt0$ alone would have lost the $349$ points of the first zone. The same computation on the positive $y$-axis, with $w=y^2-b^2$ and $F(0,y)=-(w^2-6a^2w+a^4)$, gives the qualifying points $5591\le y\le18949$, which are $13359$ points.

6. One quadrant is enough

$D$, $s$ and $F$ depend only on $x^2$ and $y^2$, so the qualifying set is symmetric under $x\mapsto-x$ and $y\mapsto-y$. Let $Q$ count the qualifying points with $x\ge0$ and $y\ge0$, $N_x$ those with $y=0,\ x\ge0$ and $N_y$ those with $x=0,\ y\ge0$, and let $N_0$ be $1$ if the origin qualifies and $0$ otherwise. The four closed quadrants contain $Q$ qualifying points each, but an axis point other than the origin lies in two of them and the origin in all four. Removing the double counts, as in the inclusion–exclusion principle, ^[8] gives

$$N=4Q-2N_x-2N_y+N_0:$$

the $2(N_x-N_0)+2(N_y-N_0)$ axis points other than the origin were counted once too often, and the origin three times too often. Here $N_0=0$, since the origin lies inside the ellipse, and Section 5 gave $N_x=7939$ and $N_y=13359$.

7. Every column is one interval

Fix a column $x\ge0$ and write $u=y^2$ and $k=a^2+b^2-x^2$, so that $s=u-k$. Then $F$ becomes a quadratic polynomial in $u$ that opens downwards,

$$G(u)=-u^2+\beta u+\gamma,\qquad\beta=4a^2+2k,$$

with $\gamma=4(b^2x^2-a^2b^2)-k^2$ and discriminant $\Delta=\beta^2+4\gamma$. If $\Delta\gt0$, $G$ is positive exactly between its two roots $u_-\lt u_+$; if $\Delta\le0$, it is never positive. We show that the qualifying $y\ge0$ of a column form a single interval $y_{\min}\le y\le y_{\max}$, possibly empty.

Outside the ellipse. $D$ grows with $y\ge0$, so $D\gt0$ holds exactly from some $y_{\min}$ on.
Columns that meet the director circle, $k\ge0$. The points with $y\ge y_{\min}$ and $y^2\le k$ lie in the ring and all qualify. At $u=k$ we have $s=0$, so $G(k)=4D\gt0$, because the director circle lies outside the ellipse. Hence $\Delta\gt0$ and $u_-\lt k\lt u_+$, and the qualifying points above the circle are exactly those with $k\lt y^2\lt u_+$; they continue the ring part without a gap.
Columns beyond the director circle, $k\lt0$. Here $y_{\min}=0$. If $G(0)\gt0$, then $\Delta\gt0$, $u_-\lt0\lt u_+$ and the column is $0\le y^2\lt u_+$. If $G(0)=F(x,0)\le0$, Section 5 gives $z=x^2-a^2\ge(3+2\sqrt2)\,b^2$, since the other possibility $z\le(3-2\sqrt2)\,b^2$ contradicts $z\gt b^2$, which is $k\lt0$. The vertex of $G$ then lies at $u=\beta/2=2a^2+b^2-z$${}\le2a^2-(2+2\sqrt2)\,b^2$, which is negative because $a^2\lt(1+\sqrt2)\,b^2$: here $a^2=56\,250\,000$ and $(1+\sqrt2)\,b^2\approx75\,444\,174$. So $G$ decreases for $u\ge0$, and the column is empty.

The condition cannot simply be dropped: for the flatter ellipse with $a=10$ and $b^2=19$ the point $(15,0)$ sees the ellipse under only $42.6^\circ$, while $(15,4)$ to $(15,13)$ qualify, and the published method, which counts that column from $y=0$, gives $1180$ points instead of the true $1164$.

8. Where the region ends

A qualifying point outside the director circle satisfies $F\gt0$. Write $X=x^2$ and $Y=y^2$. For fixed $X$, $F$ is a downward parabola in $Y$ with its vertex at $Y=3a^2+b^2-X$${}=200\,000\,000-X$. If $X\ge200\,000\,000$, then $F$ decreases for $Y\ge0$, so $F(x,0)\gt0$ is needed, and by Section 5 this means $X\lt a^2+(3+2\sqrt2)\,b^2$${}\approx238\,388\,348$, so $|x|\le15439$. If $X\lt200\,000\,000$, then $|x|\le14142$ anyway. In the same way, for fixed $Y$ the vertex in $X$ lies at $X=a^2+3b^2-Y$${}=150\,000\,000-Y$; either $|y|\le12247$, or $F(0,y)\gt0$ is needed, which requires $Y\lt b^2+(3+2\sqrt2)\,a^2$${}\approx359\,099\,026$, so $|y|\le18949$. Points of the ring have $|x|,|y|\le9354$. Every qualifying point therefore has $|x|\le15439$ and $|y|\le18949$, and both bounds are attained on the axes. The code avoids irrational numbers: it scans the columns up to $L_x$ and the $y$-axis up to $L_y$, with the larger bounds

$$\begin{aligned}L_x&=\bigl\lfloor\sqrt{a^2+6b^2}\bigr\rfloor+2=15614,\\L_y&=\bigl\lfloor\sqrt{b^2+6a^2}\bigr\rfloor+2=19204,\end{aligned}$$

which are safe because $6\gt3+2\sqrt2\approx5.83$.

9. The ends of a column in exact integers

All bounds of a column are found with the integer square root $\operatorname{isqrt}(n)=\lfloor\sqrt n\rfloor$ ^[9] and a short correction, so that no rounding error can move a point across a boundary.

$y_{\min}$: it is $0$ if $x\gt a$, and $1$ if $x=a$, because the vertex lies on the ellipse. Otherwise we need the least $y$ with $a^2y^2\gt T$, where $T=a^2b^2-b^2x^2$. The start $\operatorname{isqrt}(\lfloor T/a^2\rfloor)$ is too small, and $y$ is increased until $a^2y^2\gt T$.
$y_{\mathrm c}=\operatorname{isqrt}(k)$ is the last $y$ inside or on the director circle when $k\ge0$, and $y_{\mathrm c}=-1$ otherwise; the ring part is $y_{\min}\le y\le y_{\mathrm c}$.
$y_{\max}$: if $\Delta\le0$, there is no outer part. Otherwise the larger root of $G$ is $u_+=(\beta+\sqrt\Delta)/2$, and every integer $u\lt u_+$ satisfies $2u\lt\beta+\lfloor\sqrt\Delta\rfloor+1$, so $u$ is at most

$$u_0=\Bigl\lfloor\frac{\beta+\operatorname{isqrt}(\Delta)}{2}\Bigr\rfloor.$$

If $u_0\lt0$, there is no outer part. Otherwise $u$ is decreased while $u\ge0$ and $G(u)\le0$, which gives the largest integer $u\ge0$ with $G(u)\gt0$, or none; then $y=\operatorname{isqrt}(u)$ is decreased while $G(y^2)\le0$, which gives $y_{\max}$. The outer part is $\max(y_{\min},y_{\mathrm c}+1)\le y\le y_{\max}$, empty when the lower end exceeds the upper one.

Over all columns of the problem, the correction of $y_{\min}$ takes exactly one step in each of the $7500$ columns with $x\lt a$, $u$ is decreased $7$ times in total, and $y$ never.

10. Worked example: the ellipse shrunk 2500 times

For $a=3$ and $b^2=5$ we have $D=5x^2+9y^2-45$, $s=x^2+y^2-14$ and $F=4D-s^2$, and the condition of Section 7 holds, since $9\lt5(1+\sqrt2)\approx12.07$. In the column $x=0$, being outside means $y^2\gt5$, so $y_{\min}=3$, and $y_{\mathrm c}=\operatorname{isqrt}(14)=3$. For the upper end, $k=14$, $\beta=64$ and $\gamma=-376$, so $G(u)=-u^2+64u-376$ with $\Delta=2592$ and $\operatorname{isqrt}(\Delta)=50$. Hence $u_0=(64+50)/2=57$, and $G(57)=23\gt0$ needs no correction, so $y_{\max}=\operatorname{isqrt}(57)=7$; indeed $G(49)=359\gt0$ and $G(64)=-376\lt0$. The points with $y=3,\dots,7$ see the ellipse under $112.62^\circ$, $84.26^\circ$, $67.71^\circ$, $56.63^\circ$ and $48.67^\circ$, and $(0,8)$ under only $42.67^\circ$. For the columns $x=0,1,\dots,6$ the interval runs from $y_{\min}=3,3,2,1,0,0,0$ to $y_{\max}=7,7,7,6,5,4,2$; the point $(3,0)$ lies on the ellipse, and no column $x\ge7$ contains a point, since $49\gt9+(3+2\sqrt2)\cdot5\approx38.1$. So $Q=5+5+6+6+6+5+3=36$, $N_x=3$, $N_y=5$ and $N_0=0$, and

$$N=4\cdot36-2\cdot3-2\cdot5=128.$$

How the Code Works

The C++ program ^[10] starts in main, which first calls run_validation unless the option --no-validate is given. That function builds the centred ellipses with $a=r/2$ and $b^2=a^2-(d/2)^2$ for $(r,d)=(6,2)$, $(10,4)$, $(14,6)$, $(18,8)$ and stops with an error message unless count_points returns $154$, $418$, $810$ and $1334$. Then build_ellipse turns the data of the statement into the ellipse of Section 1. It requires a horizontal axis, a centre on the lattice, a radius $r$ above the focal distance $d=|G_x-M_x|=10000$, even $r$ and $d$, and $b^2\gt0$, and stores a2, b2, r2_circle $=a^2+b^2$ and ab2 $=a^2b^2$ in an EllipseData.

count_points obtains $Q$ from count_quadrant, counts $N_x$, $N_y$ and $N_0$ with qualifies, the single-point test of Section 5, along the half-axes up to $L_x$ and $L_y$, and returns $4Q-2N_x-2N_y+N_0$. count_quadrant runs through the columns $x=0,\dots,L_x$; for each it calls min_y_outside for $y_{\min}$, computes $y_{\mathrm c}$ and adds the ring part, then calls max_y_outside for $y_{\max}$ and adds the outer part (Section 9). angle_ineq(x, u) evaluates $G(u)$, that is $F$ with $u=y^2$. max_y_outside forms $k$, $\beta$, $\gamma$ and $\Delta$ in the variables k, b, c and d and carries out Section 9, returning $-1$ when there is no outer part. isqrt_u128 corrects a long double square root with two loops until $r^2\le n\lt(r+1)^2$, so its result is exact.

C++ keeps coordinates and counts in 64-bit integers, forms small products such as x * x in 64 bits and the large geometric expressions in the 128-bit types __int128_t and __uint128_t, a compiler extension. ^[11] Here every intermediate value stays below $1.7\cdot10^{17}$, the largest being $\beta^2=1.6\cdot10^{17}$, so signed 64-bit integers would suffice; the wider types are a safety margin. In the last $175$ columns, $x\ge15440$, the numerator $\beta+\operatorname{isqrt}(\Delta)$ is at most $-24189$; C++ and Java halve it with truncation toward zero and Python rounds down, but in all three the start is negative and max_y_outside returns $-1$ at once. The Python program has the same helper functions, with unbounded integers and math.isqrt; its worker adds the counts of the columns in a range of $x$, and solve applies the axis correction. The Java program keeps coordinates and counts in long and the large expressions in BigInteger, and its isqrtU128 corrects a Math.sqrt guess with exact comparisons. Neither Python nor Java runs the validation cases.

Complexity and Verification

The program examines $L_x+1=15615$ columns, each with a fixed number of integer operations, at most four integer square roots and at most one step of each correction loop (Section 9); the axis counts add $15615+19205$ single-point tests. Counting arithmetic on numbers of this size as constant-time, the work grows linearly when the ellipse is enlarged, while testing every lattice point of the bounding box of the region, here $30879\cdot37899\approx1.2\cdot10^9$ points, would grow with its area. Each column needs a fixed number of integers, so the memory does not grow with the ellipse, and in our runs all three programs finished in well under a second. The method relies on the condition $a^2\lt(1+\sqrt2)\,b^2$ of Section 7, which holds for the problem and all four validation ellipses but is not checked by build_ellipse.

Besides the self-test of run_validation, the following was verified for this article:

The four reference counts and the $128$ of Section 10 agree with two brute forces over a large square: the exact test of Section 5, and a geometric computation that locates $R$ and $S$ on the parametrisation $(a\cos t,b\sin t)$ and measures $\angle RPS$ from $\vec{PR}$ and $\vec{PS}$ with $60$-digit arithmetic, without the angle formula.
For the problem's ellipse, the exact test on all of the roughly $2\cdot10^9$ lattice points with $|x|\le20000$ and $|y|\le25000$, without symmetry or column formulas, reproduces the programs' total and every column count; every column is one interval, with extreme coordinates $15439$ and $18949$.
The angle from the actual tangent points, in double precision, classifies every lattice point with $|x|\le16000$ and $|y|\le19500$ like the exact test. No outside lattice point has an angle within $2\cdot10^{-9}$ radians of $45^\circ$, none lies on the director circle, and only the two vertices lie on the ellipse.
The C++, Python and Java programs print the same number.

Footnotes and References

The references below supply the background facts used above. The single-point test, the interval property of the columns, the extent of the region and all counts are derived in this article.

Project Euler 246 — Tangents to an Ellipse. The official statement defines the ellipse as the locus of the points equidistant from a circle c (centre M, radius r) and a point G with d(G, M) < r, gives M(−2000, 1500), G(8000, 1500) and r = 15000, and asks for how many lattice points P outside the ellipse the angle RPS between the points of contact R and S of the two tangents is greater than 45 degrees.
Wikipedia — Ellipse. An ellipse is the set of points whose distances to two foci add up to 2a, and also the set of points as far from a circular directrix, the circle of radius 2a about one focus, as from the other focus. Its semi-axes satisfy b² = a² − c², with c the distance from the centre to a focus, its standard equation is x²/a² + y²/b² = 1, and its tangents of slope m are y = mx ± √(m²a² + b²).
Wikipedia — List of trigonometric identities. Among the angle sum and difference identities is tan(α − β) = (tan α − tan β) / (1 + tan α tan β), and among the reflections tan(π − θ) = −tan θ.
Wikipedia — Vieta's formulas. The roots r₁ and r₂ of ax² + bx + c satisfy r₁ + r₂ = −b/a and r₁r₂ = c/a.
Wikipedia — Orthoptic (geometry). An isoptic is the set of points where two tangents of a curve meet at a fixed angle. The α-isoptics of the ellipse x²/a² + y²/b² = 1 are the two parts of the degree-4 curve (x² + y² − a² − b²)² tan²α = 4(a²y² + b²x² − a²b²); the proof substitutes the two tangent slopes into tan²α = ((m₁ − m₂)/(1 + m₁m₂))². The orthoptic, where the tangents are perpendicular, is the director circle.
Wikipedia — Director circle. The director circle of an ellipse consists of the points where two perpendicular tangents of the ellipse cross; it has the centre of the ellipse and the radius √(a² + b²), and from each of its points the ellipse spans a right angle.
Wikipedia — Intermediate value theorem. A continuous function on [a, b] takes every value between f(a) and f(b); more generally, continuous maps send connected sets to connected sets, the form used in Section 4.
Wikipedia — Inclusion–exclusion principle. For finite sets |A ∪ B| = |A| + |B| − |A ∩ B|, and in general the size of a union is the alternating sum of the sizes of all intersections.
Wikipedia — Integer square root. The integer square root of n ≥ 0 is the greatest integer m with m² ≤ n, that is ⌊√n⌋; the article describes exact algorithms, one with a final verification and correction step.
C++ — Euler246.cpp. The linked immutable C++ revision contains min_y_outside, angle_ineq, max_y_outside, qualifies, count_quadrant, count_points, build_ellipse and run_validation with its four reference counts. It is the implementation record for this article.
GCC manual — 128-bit Integers. GCC provides __int128 and unsigned __int128 on targets with an integer mode wide enough to hold 128 bits.

Problem 246 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

namespace {

using i64 = long long;
using u64 = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;

constexpr i64 kMx = -2000;
constexpr i64 kMy = 1500;
constexpr i64 kGx = 8000;
constexpr i64 kGy = 1500;
constexpr i64 kRadius = 15000;

u64 isqrt_u128(u128 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while (static_cast<u128>(r + 1) * (r + 1) <= x) ++r;
    while (static_cast<u128>(r) * r > x) --r;
    return r;
}

struct EllipseData {
    i64 a2 = 0;
    i64 b2 = 0;
    i64 r2_circle = 0;
    i128 ab2 = 0;
};

i64 min_y_outside(i64 x, const EllipseData& e) {
    u128 lhs = static_cast<u128>(e.b2) * x * x;
    u128 rhs = static_cast<u128>(e.a2) * e.b2;
    if (lhs > rhs) return 0;
    if (lhs == rhs) return 1;
    u128 rem = rhs - lhs;
    u128 q = rem / static_cast<u128>(e.a2);
    i64 y = static_cast<i64>(isqrt_u128(q));
    while (static_cast<u128>(e.a2) * y * y <= rem) ++y;
    return y;
}

i128 angle_ineq(i64 x, i128 u, const EllipseData& e) {
    // Positive when the acute angle between tangents is greater than 45 degrees.
    i128 x2 = static_cast<i128>(x) * x;
    i128 term1 = 4 * (static_cast<i128>(e.b2) * x2 +
                      static_cast<i128>(e.a2) * u - e.ab2);
    i128 t = static_cast<i128>(e.a2) + e.b2 - x2 - u;
    return term1 - t * t;
}

i64 max_y_outside(i64 x, const EllipseData& e) {
    i128 x2 = static_cast<i128>(x) * x;
    i128 k = static_cast<i128>(e.a2) + e.b2 - x2;
    i128 b = 4 * static_cast<i128>(e.a2) + 2 * k;
    i128 c = 4 * (static_cast<i128>(e.b2) * x2 - e.ab2) - k * k;
    i128 d = b * b + 4 * c;
    if (d <= 0) return -1;
    u128 sqrt_d = isqrt_u128(static_cast<u128>(d));
    i128 u_max = (b + static_cast<i128>(sqrt_d)) / 2;
    if (u_max < 0) return -1;
    while (u_max >= 0 && angle_ineq(x, u_max, e) <= 0) --u_max;
    if (u_max < 0) return -1;
    i64 y = static_cast<i64>(isqrt_u128(static_cast<u128>(u_max)));
    while (y >= 0 && angle_ineq(x, static_cast<i128>(y) * y, e) <= 0) --y;
    return y;
}

bool qualifies(i64 x, i64 y, const EllipseData& e) {
    i128 x2 = static_cast<i128>(x) * x;
    i128 y2 = static_cast<i128>(y) * y;
    if (static_cast<i128>(e.b2) * x2 + static_cast<i128>(e.a2) * y2 <= e.ab2) {
        return false;
    }
    i128 r2 = x2 + y2;
    if (r2 <= e.r2_circle) return true;
    return angle_ineq(x, y2, e) > 0;
}

i64 count_quadrant(const EllipseData& e, int threads) {
    i64 x_limit = static_cast<i64>(isqrt_u128(static_cast<u128>(e.a2) +
                                             static_cast<u128>(6) * e.b2)) +
                  2;
    int use_threads = std::max(1, threads);
    use_threads = std::min<i64>(use_threads, x_limit + 1);
    std::vector<i64> locals(use_threads, 0);

    auto worker = [&](int tid, i64 start, i64 end) {
        i64 local = 0;
        for (i64 x = start; x < end; ++x) {
            i64 y_min = min_y_outside(x, e);
            i64 y_circle = -1;
            i64 x2 = x * x;
            if (x2 <= e.r2_circle) {
                // Inside the director circle (r^2 <= a^2 + b^2), the tangent angle is obtuse.
                y_circle = static_cast<i64>(isqrt_u128(e.r2_circle - x2));
            }
            if (y_circle >= y_min) {
                local += y_circle - y_min + 1;
            }
            i64 y_start = std::max(y_min, y_circle + 1);
            i64 y_max = max_y_outside(x, e);
            if (y_max >= y_start) {
                local += y_max - y_start + 1;
            }
        }
        locals[tid] = local;
    };

    i64 total_x = x_limit + 1;
    i64 chunk = (total_x + use_threads - 1) / use_threads;
    std::vector<std::thread> pool;
    pool.reserve(use_threads);
    for (int t = 0; t < use_threads; ++t) {
        i64 start = t * chunk;
        i64 end = std::min(start + chunk, total_x);
        if (start >= end) break;
        pool.emplace_back(worker, t, start, end);
    }
    for (auto& th : pool) th.join();

    i64 sum = 0;
    for (i64 v : locals) sum += v;
    return sum;
}

i64 count_points(const EllipseData& e, int threads) {
    i64 count_q1 = count_quadrant(e, threads);

    i64 x_limit = static_cast<i64>(isqrt_u128(static_cast<u128>(e.a2) +
                                             static_cast<u128>(6) * e.b2)) +
                  2;
    i64 y_limit = static_cast<i64>(isqrt_u128(static_cast<u128>(e.b2) +
                                             static_cast<u128>(6) * e.a2)) +
                  2;

    i64 count_x_axis = 0;
    for (i64 x = 0; x <= x_limit; ++x) {
        if (qualifies(x, 0, e)) ++count_x_axis;
    }
    i64 count_y_axis = 0;
    for (i64 y = 0; y <= y_limit; ++y) {
        if (qualifies(0, y, e)) ++count_y_axis;
    }
    i64 origin = qualifies(0, 0, e) ? 1 : 0;

    return 4 * count_q1 - 2 * count_x_axis - 2 * count_y_axis + origin;
}

EllipseData build_ellipse(i64 mx, i64 my, i64 gx, i64 gy, i64 r) {
    if (my != gy) {
        std::cerr << "Ellipse axis is not horizontal.\n";
        std::exit(1);
    }
    if ((mx + gx) % 2 != 0 || (my + gy) % 2 != 0) {
        std::cerr << "Ellipse center is not on lattice.\n";
        std::exit(1);
    }
    i64 d = std::llabs(gx - mx);
    if (r <= d) {
        std::cerr << "Radius must exceed focal distance.\n";
        std::exit(1);
    }
    if ((r % 2) != 0 || (d % 2) != 0) {
        std::cerr << "Radius and focal distance must be even.\n";
        std::exit(1);
    }
    i64 a = r / 2;
    i64 c = d / 2;
    i64 a2 = a * a;
    i64 b2 = a2 - c * c;
    if (b2 <= 0) {
        std::cerr << "Invalid ellipse parameters.\n";
        std::exit(1);
    }
    EllipseData e;
    e.a2 = a2;
    e.b2 = b2;
    e.r2_circle = a2 + b2;
    e.ab2 = static_cast<i128>(a2) * b2;
    return e;
}

void run_validation() {
    struct Check {
        i64 r;
        i64 d;
        i64 expected;
    };
    const Check checks[] = {
        {6, 2, 154},
        {10, 4, 418},
        {14, 6, 810},
        {18, 8, 1334},
    };
    for (const auto& chk : checks) {
        i64 a = chk.r / 2;
        i64 c = chk.d / 2;
        EllipseData e;
        e.a2 = a * a;
        e.b2 = e.a2 - c * c;
        e.r2_circle = e.a2 + e.b2;
        e.ab2 = static_cast<i128>(e.a2) * e.b2;
        i64 got = count_points(e, 1);
        if (got != chk.expected) {
            std::cerr << "Validation failed for r=" << chk.r << " d=" << chk.d
                      << ": got " << got << ", expected " << chk.expected << '\n';
            std::exit(1);
        }
    }
}

void usage(const char* argv0) {
    std::cerr << "Usage: " << argv0 << " [-t THREADS] [--no-validate]\n";
}

}  // namespace

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

    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;
    bool validate = true;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "-t" && i + 1 < argc) {
            threads = std::max(1, std::stoi(argv[++i]));
        } else if (arg == "--no-validate") {
            validate = false;
        } else if (arg == "--validate") {
            validate = true;
        } else {
            usage(argv[0]);
            return 1;
        }
    }

    if (validate) run_validation();

    EllipseData ellipse = build_ellipse(kMx, kMy, kGx, kGy, kRadius);
    i64 answer = count_points(ellipse, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

import math
import multiprocessing

def isqrt_u128(x):
    if x < 0: return 0
    return math.isqrt(x)

class EllipseData:
    def __init__(self, a2, b2):
        self.a2 = a2
        self.b2 = b2
        self.r2_circle = a2 + b2
        self.ab2 = a2 * b2

def min_y_outside(x, e):
    lhs = e.b2 * x * x
    rhs = e.a2 * e.b2
    if lhs > rhs: return 0
    if lhs == rhs: return 1
    rem = rhs - lhs
    q = rem // e.a2
    y = isqrt_u128(q)
    while e.a2 * y * y <= rem:
        y += 1
    return y

def angle_ineq(x, u, e):
    x2 = x * x
    term1 = 4 * (e.b2 * x2 + e.a2 * u - e.ab2)
    t = e.a2 + e.b2 - x2 - u
    return term1 - t * t

def max_y_outside(x, e):
    x2 = x * x
    k = e.a2 + e.b2 - x2
    b = 4 * e.a2 + 2 * k
    c = 4 * (e.b2 * x2 - e.ab2) - k * k
    d = b * b + 4 * c
    if d <= 0: return -1
    
    sqrt_d = isqrt_u128(d)
    u_max = (b + sqrt_d) // 2
    if u_max < 0: return -1
    
    while u_max >= 0 and angle_ineq(x, u_max, e) <= 0:
        u_max -= 1
        
    if u_max < 0: return -1
    y = isqrt_u128(u_max)
    while y >= 0 and angle_ineq(x, y * y, e) <= 0:
        y -= 1
        
    return y

def qualifies(x, y, e):
    x2 = x * x
    y2 = y * y
    if e.b2 * x2 + e.a2 * y2 <= e.ab2:
        return False
    r2 = x2 + y2
    if r2 <= e.r2_circle:
        return True
    return angle_ineq(x, y2, e) > 0

def worker(args):
    start, end, e = args
    local = 0
    for x in range(start, end):
        y_min = min_y_outside(x, e)
        y_circle = -1
        x2 = x * x
        if x2 <= e.r2_circle:
            y_circle = isqrt_u128(e.r2_circle - x2)
            
        if y_circle >= y_min:
            local += y_circle - y_min + 1
            
        y_start = max(y_min, y_circle + 1)
        y_max = max_y_outside(x, e)
        
        if y_max >= y_start:
            local += y_max - y_start + 1
            
    return local

def solve():
    kMx = -2000
    kMy = 1500
    kGx = 8000
    kGy = 1500
    kRadius = 15000
    
    d = abs(kGx - kMx)
    a = kRadius // 2
    c = d // 2
    a2 = a * a
    b2 = a2 - c * c
    
    e = EllipseData(a2, b2)
    
    x_limit = isqrt_u128(e.a2 + 6 * e.b2) + 2
    y_limit = isqrt_u128(e.b2 + 6 * e.a2) + 2
    
    threads = max(1, multiprocessing.cpu_count())
    total_x = x_limit + 1
    chunk = (total_x + threads - 1) // threads
    
    tasks = []
    for t in range(threads):
        start = t * chunk
        end = min(start + chunk, total_x)
        if start >= end: break
        tasks.append((start, end, e))
        
    count_q1 = 0
    if threads > 1 and len(tasks) > 1:
        with multiprocessing.Pool(threads) as pool:
            results = pool.map(worker, tasks)
            count_q1 = sum(results)
    else:
        for t in tasks:
            count_q1 += worker(t)
            
    count_x_axis = sum(1 for x in range(x_limit + 1) if qualifies(x, 0, e))
    count_y_axis = sum(1 for y in range(y_limit + 1) if qualifies(0, y, e))
    origin = 1 if qualifies(0, 0, e) else 0
    
    ans = 4 * count_q1 - 2 * count_x_axis - 2 * count_y_axis + origin
    return str(ans)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.*;

public class Euler246 {

    static long isqrtU128(BigInteger x) {
        if (x.compareTo(BigInteger.ZERO) <= 0)
            return 0;
        // BigInteger has sqrt() in newer Java, but let's just use double as an initial
        // guess, then refine.
        long r = (long) Math.sqrt(x.doubleValue());

        while (BigInteger.valueOf(r + 1).multiply(BigInteger.valueOf(r + 1)).compareTo(x) <= 0) {
            r++;
        }
        while (BigInteger.valueOf(r).multiply(BigInteger.valueOf(r)).compareTo(x) > 0) {
            r--;
        }
        return r;
    }

    static class EllipseData {
        long a2, b2;
        long r2Circle;
        BigInteger ab2;

        EllipseData(long a2, long b2) {
            this.a2 = a2;
            this.b2 = b2;
            this.r2Circle = a2 + b2;
            this.ab2 = BigInteger.valueOf(a2).multiply(BigInteger.valueOf(b2));
        }
    }

    static long minYOutside(long x, EllipseData e) {
        BigInteger lhs = BigInteger.valueOf(e.b2).multiply(BigInteger.valueOf(x)).multiply(BigInteger.valueOf(x));
        BigInteger rhs = e.ab2;
        if (lhs.compareTo(rhs) > 0)
            return 0;
        if (lhs.equals(rhs))
            return 1;

        BigInteger rem = rhs.subtract(lhs);
        BigInteger q = rem.divide(BigInteger.valueOf(e.a2));
        long y = isqrtU128(q);

        while (BigInteger.valueOf(e.a2).multiply(BigInteger.valueOf(y)).multiply(BigInteger.valueOf(y))
                .compareTo(rem) <= 0) {
            y++;
        }
        return y;
    }

    static BigInteger angleIneq(long x, BigInteger u, EllipseData e) {
        BigInteger x2 = BigInteger.valueOf(x).multiply(BigInteger.valueOf(x));

        BigInteger term1Inner = BigInteger.valueOf(e.b2).multiply(x2)
                .add(BigInteger.valueOf(e.a2).multiply(u))
                .subtract(e.ab2);

        BigInteger term1 = term1Inner.multiply(BigInteger.valueOf(4));

        BigInteger t = BigInteger.valueOf(e.a2).add(BigInteger.valueOf(e.b2))
                .subtract(x2).subtract(u);

        return term1.subtract(t.multiply(t));
    }

    static long maxYOutside(long x, EllipseData e) {
        BigInteger x2 = BigInteger.valueOf(x).multiply(BigInteger.valueOf(x));

        BigInteger k = BigInteger.valueOf(e.a2).add(BigInteger.valueOf(e.b2)).subtract(x2);
        BigInteger b = BigInteger.valueOf(4 * e.a2).add(k.multiply(BigInteger.valueOf(2)));

        BigInteger cInner = BigInteger.valueOf(e.b2).multiply(x2).subtract(e.ab2);
        BigInteger c = cInner.multiply(BigInteger.valueOf(4)).subtract(k.multiply(k));

        BigInteger d = b.multiply(b).add(c.multiply(BigInteger.valueOf(4)));

        if (d.compareTo(BigInteger.ZERO) <= 0)
            return -1;

        long sqrtD = isqrtU128(d);
        BigInteger uMax = b.add(BigInteger.valueOf(sqrtD)).divide(BigInteger.valueOf(2));

        if (uMax.compareTo(BigInteger.ZERO) < 0)
            return -1;

        while (uMax.compareTo(BigInteger.ZERO) >= 0 && angleIneq(x, uMax, e).compareTo(BigInteger.ZERO) <= 0) {
            uMax = uMax.subtract(BigInteger.ONE);
        }

        if (uMax.compareTo(BigInteger.ZERO) < 0)
            return -1;

        long y = isqrtU128(uMax);
        while (y >= 0 && angleIneq(x, BigInteger.valueOf(y).multiply(BigInteger.valueOf(y)), e)
                .compareTo(BigInteger.ZERO) <= 0) {
            y--;
        }

        return y;
    }

    static boolean qualifies(long x, long y, EllipseData e) {
        BigInteger x2 = BigInteger.valueOf(x).multiply(BigInteger.valueOf(x));
        BigInteger y2 = BigInteger.valueOf(y).multiply(BigInteger.valueOf(y));

        BigInteger lhs = BigInteger.valueOf(e.b2).multiply(x2).add(BigInteger.valueOf(e.a2).multiply(y2));
        if (lhs.compareTo(e.ab2) <= 0)
            return false;

        BigInteger r2 = x2.add(y2);
        if (r2.compareTo(BigInteger.valueOf(e.r2Circle)) <= 0)
            return true;

        return angleIneq(x, y2, e).compareTo(BigInteger.ZERO) > 0;
    }

    public static String solve() {
        long kMx = -2000;
        long kGx = 8000;
        long kRadius = 15000;

        long d = Math.abs(kGx - kMx);
        long a = kRadius / 2;
        long c = d / 2;
        long a2 = a * a;
        long b2 = a2 - c * c;

        EllipseData e = new EllipseData(a2, b2);

        long xLimit = isqrtU128(BigInteger.valueOf(e.a2).add(BigInteger.valueOf(6 * e.b2))) + 2;
        long yLimit = isqrtU128(BigInteger.valueOf(e.b2).add(BigInteger.valueOf(6 * e.a2))) + 2;

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        long totalX = xLimit + 1;
        long chunk = (totalX + threads - 1) / threads;

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int t = 0; t < threads; ++t) {
            final long start = t * chunk;
            final long end = Math.min(start + chunk, totalX);
            if (start >= end)
                break;

            futures.add(executor.submit(() -> {
                long local = 0;
                for (long x = start; x < end; ++x) {
                    long yMin = minYOutside(x, e);
                    long yCircle = -1;
                    long x2 = x * x;
                    if (x2 <= e.r2Circle) {
                        yCircle = isqrtU128(BigInteger.valueOf(e.r2Circle).subtract(BigInteger.valueOf(x2)));
                    }
                    if (yCircle >= yMin) {
                        local += yCircle - yMin + 1;
                    }

                    long yStart = Math.max(yMin, yCircle + 1);
                    long yMax = maxYOutside(x, e);
                    if (yMax >= yStart) {
                        local += yMax - yStart + 1;
                    }
                }
                return local;
            }));
        }

        long countQ1 = 0;
        for (Future<Long> f : futures) {
            try {
                countQ1 += f.get();
            } catch (Exception ex) {
            }
        }
        executor.shutdown();

        long countXAxis = 0;
        for (long x = 0; x <= xLimit; ++x) {
            if (qualifies(x, 0, e))
                countXAxis++;
        }

        long countYAxis = 0;
        for (long y = 0; y <= yLimit; ++y) {
            if (qualifies(0, y, e))
                countYAxis++;
        }

        long origin = qualifies(0, 0, e) ? 1 : 0;

        long ans = 4 * countQ1 - 2 * countXAxis - 2 * countYAxis + origin;
        return String.valueOf(ans);
    }

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