Project Euler

Tangents to an Ellipse

Given points M(-2000, 1500) and G(8000, 1500) and a circle c centred at M with radius 15000, the locus of points equidistant from G and from c forms an ellipse e. From an external point P, two tang...

Source sync Apr 19, 2026

Problem #0246

Level Level 15

Solved By 993

Languages C++, Python

Answer 810834388

Length 326 words

geometryalgebrasearch

A definition for an ellipse is:

Given a circle $c$ with centre $M$ and radius $r$ and a point $G$ such that $d(G,M) < r$, the locus of the points that are equidistant from $c$ and $G$ form an ellipse.

The construction of the points of the ellipse is shown below.

Given are the points $M(-2000,1500)$ and $G(8000,1500)$.

Given is also the circle $c$ with centre $M$ and radius $15000$.

The locus of the points that are equidistant from $G$ and $c$ form an ellipse $e$.

From a point $P$ outside $e$ the two tangents $t_1$ and $t_2$ to the ellipse are drawn.

Let the points where $t_1$ and $t_2$ touch the ellipse be $R$ and $S$.

Problem illustration

For how many lattice points $P$ is angle $RPS$ greater than $45$ degrees?

Problem 246: Tangents to an Ellipse

Mathematical Foundation

Theorem 1 (Ellipse from equidistance condition). A point $Q$ is equidistant from the point $G$ and from the circle $c$ (centred at $M$ with radius $r$ ) iff $d(Q, G) = r - d(Q, M)$ , i.e., $d(Q, M) + d(Q, G) = r$ . This is the defining property of an ellipse with foci $M$ and $G$ and semi-major axis $a = r/2$ .

Proof. The distance from $Q$ to the circle $c$ is $|r - d(Q, M)|$ . For $Q$ inside $c$ , this equals $r - d(Q, M)$ . Setting this equal to $d(Q, G)$ gives $d(Q, M) + d(Q, G) = r = 15000$ . By definition, this is an ellipse with foci $M, G$ and major axis $2a = r$ , so $a = 7500$ . $\square$

Lemma 1 (Ellipse parameters). The ellipse $e$ has:

Centre: $O = \left(\frac{-2000 + 8000}{2}, 1500\right) = (3000, 1500)$
Semi-major axis: $a = 7500$
Half-focal distance: $c = \frac{d(M, G)}{2} = \frac{10000}{2} = 5000$
Semi-minor axis: $b = \sqrt{a^2 - c^2} = \sqrt{56250000 - 25000000} = \sqrt{31250000} = 2500\sqrt{5}$

Proof. Direct computation from $d(M, G) = \sqrt{10000^2 + 0^2} = 10000$ , and $b^2 = a^2 - c^2$ . $\square$

Theorem 2 (Isoptic curve). The $\alpha$ -isoptic of the ellipse $\frac{X^2}{a^2} + \frac{Y^2}{b^2} = 1$ (in centred coordinates) is the curve

(X^2 + Y^2 - a^2 - b^2)^2 = 4\cot^2\alpha\,(a^2 Y^2 + b^2 X^2 - a^2 b^2).

For $\alpha = 45°$ , $\cot^2(45°) = 1$ , so the isoptic simplifies to:

f(X, Y) = (X^2 + Y^2 - a^2 - b^2)^2 - 4(a^2 Y^2 + b^2 X^2 - a^2 b^2) = 0.

Proof. The isoptic curve of a conic is a classical result. From a point $(X, Y)$ outside the ellipse, the tangent lines to $\frac{X^2}{a^2} + \frac{Y^2}{b^2} = 1$ have slopes satisfying a quadratic derived from the tangency condition. The angle between them is $\alpha$ iff $\tan\alpha = |m_1 - m_2|/(1 + m_1 m_2)$ , which after algebraic manipulation using the discriminant of the tangent-slope quadratic yields the stated equation. $\square$

Lemma 2 (Counting criterion). A lattice point $P = (x, y)$ in original coordinates satisfies $\angle RPS > 45°$ iff, setting $X = x - 3000$ , $Y = y - 1500$ :

$P$ is strictly outside the ellipse: $\frac{X^2}{a^2} + \frac{Y^2}{b^2} > 1$ , and
$P$ is strictly inside the $45°$ -isoptic: $f(X, Y) < 0$ .

Proof. The angle subtended by the tangent lines decreases as $P$ moves farther from the ellipse. At the isoptic $f = 0$ , the angle equals $45°$ . Inside the isoptic (but outside the ellipse), the angle exceeds $45°$ . $\square$

Editorial

We determine Y range from isoptic extent. We then the isoptic is bounded; find max |Y| by solving f(0, Y) = 0. Finally, solve f(X, Y) = 0 as quadratic in u = X^2.

Pseudocode

Determine Y range from isoptic extent
The isoptic is bounded; find max |Y| by solving f(0, Y) = 0
Solve f(X, Y) = 0 as quadratic in u = X^2
f = (u + Y^2 - a^2 - b^2)^2 - 4(a^2*Y^2 + b^2*u - a^2*b^2) = 0
Expand and solve for u
Count integers X in [-X_iso, X_iso] that are outside ellipse
i.e., |X| > X_ell or X_ell doesn't exist (Y too large for ellipse)

Complexity Analysis

Time: $O(H)$ where $H$ is the vertical extent of the isoptic curve, approximately $2(a + b) \approx 2(7500 + 5590) \approx 26180$ . For each row, solving the quartic (as a quadratic in $X^2$ ) takes $O(1)$ .
Space: $O(1)$ .

Answer

\boxed{810834388}

C++ project_euler/problem_246/solution.cpp

#include <bits/stdc++.h>
using namespace std;

int main() {
    // Problem 246: Tangents to an Ellipse
    // Count lattice points P outside ellipse where angle RPS > 45 degrees
    // Ellipse: center (3000,1500), a=7500, b^2=31250000

    const long long A2 = 56250000LL;
    const long long B2 = 31250000LL;
    const long long S = A2 + B2;
    const long long A2B2 = A2 * B2;
    const int cx = 3000, cy = 1500;

    // f(X,Y) = (X^2+Y^2-S)^2 - 4*(A^2*Y^2+B^2*X^2-A^2*B^2)
    // angle > 45 when: outside ellipse AND (inside director circle OR f < 0)
    // Equivalently: X^2*B2+Y^2*A2 > A2*B2 AND (X^2+Y^2 <= S OR f < 0)

    // For boundary detection, use exact integer arithmetic with __int128.
    // For each Y, find the max |X| where the condition holds by binary search.

    auto isInside = [&](long long X, long long Y2) -> bool {
        // Check: outside ellipse AND (inside director circle OR f < 0)
        long long X2 = X * X;
        __int128 ev = (__int128)X2 * B2 + (__int128)Y2 * A2;
        if (ev <= (__int128)A2B2) return false; // inside or on ellipse

        long long r2 = X2 + Y2;
        if (r2 <= S) return true; // inside director circle

        // Check f < 0
        __int128 diff = r2 - S;
        __int128 lhs = diff * diff;
        __int128 rhs = 4 * ((__int128)A2 * Y2 + (__int128)B2 * X2 - (__int128)A2B2);
        return lhs < rhs;
    };

    long long count = 0;

    // Y range: b ~ 5590, isoptic extends to ~18950 from center
    for (int y = -17500; y <= 20500; y++) {
        long long Y = y - cy;
        long long Y2 = Y * Y;

        // Get approximate iso_hi using floating point
        double p_val = 2.0*Y2 - 2.0*A2 - 6.0*B2;
        double Y4 = (double)Y2 * Y2;
        double q_val = Y4 - (6.0*A2+2.0*B2)*Y2 + (double)A2*A2 + 6.0*A2*B2 + (double)B2*B2;
        double disc = p_val*p_val - 4*q_val;
        if (disc < 0) continue;
        double iso_hi_sq = (-p_val + sqrt(disc)) / 2.0;
        if (iso_hi_sq <= 0) continue;
        double iso_hi_approx = sqrt(iso_hi_sq);

        // Get approximate ell_lo
        long long ell_rhs_num = A2 * (B2 - Y2);
        double ell_lo_approx;
        if (ell_rhs_num <= 0) {
            ell_lo_approx = 0;
        } else {
            ell_lo_approx = sqrt((double)ell_rhs_num / B2);
        }

        if (ell_lo_approx >= iso_hi_approx + 2) continue;

        // Positive X side: find exact range using exact checks
        // Start from approximate boundaries and adjust
        int x_lo = (int)ceil(cx + ell_lo_approx);
        int x_hi = (int)floor(cx + iso_hi_approx) + 1; // +1 for safety

        // Adjust x_lo: find first x where isInside
        // x_lo should be the smallest x >= cx + ell_lo_approx that's inside
        while (x_lo <= x_hi && !isInside(x_lo - cx, Y2)) x_lo++;

        // Adjust x_hi: find last x where isInside
        // Start from approximate and search
        while (x_hi >= x_lo && !isInside(x_hi - cx, Y2)) x_hi--;

        // Verify we haven't missed points above x_hi
        // (unlikely given the +1 safety margin, but check)
        while (x_hi + 1 <= cx + (int)iso_hi_approx + 3 && isInside(x_hi + 1 - cx, Y2)) x_hi++;

        if (x_lo <= x_hi)
            count += x_hi - x_lo + 1;

        // Negative X side (use cx-1 as upper bound when ell_lo=0 to avoid
        // double-counting the x=cx column already covered by positive side)
        int x_lo2 = (int)ceil(cx - iso_hi_approx) - 1; // -1 for safety
        int x_hi2 = (ell_lo_approx < 0.5) ? cx - 1 : (int)floor(cx - ell_lo_approx);

        while (x_lo2 <= x_hi2 && !isInside(x_lo2 - cx, Y2)) x_lo2++;
        while (x_hi2 >= x_lo2 && !isInside(x_hi2 - cx, Y2)) x_hi2--;
        while (x_lo2 - 1 >= cx - (int)iso_hi_approx - 3 && isInside(x_lo2 - 1 - cx, Y2)) x_lo2--;

        if (x_lo2 <= x_hi2)
            count += x_hi2 - x_lo2 + 1;
    }

    cout << count << endl;
    return 0;
}

Python project_euler/problem_246/solution.py

import math

def solve():
    """
    Problem 246: Tangents to an Ellipse
    Ellipse: foci M=(-2000,1500), G=(8000,1500), sum of distances = 15000
    Center (3000,1500), a=7500, c=5000, b^2=31250000
    Count lattice points P outside ellipse where angle RPS > 45 degrees.
    Uses isoptic curve theory with exact integer arithmetic at boundaries.
    """
    A2 = 56250000     # a^2 = 7500^2
    B2 = 31250000     # b^2
    S_val = A2 + B2   # a^2 + b^2 = 87500000
    A2B2 = A2 * B2    # 1757812500000000
    cx, cy = 3000, 1500

    def is_valid(X, Y2):
        """Check if lattice point (cx+X, y) with Y^2=Y2 satisfies angle > 45."""
        X2 = X * X
        # Outside ellipse?
        ev = X2 * B2 + Y2 * A2
        if ev <= A2B2:
            return False
        # Inside director circle?
        r2 = X2 + Y2
        if r2 <= S_val:
            return True
        # Inside outer isoptic? f < 0?
        diff = r2 - S_val
        lhs = diff * diff
        rhs = 4 * (A2 * Y2 + B2 * X2 - A2B2)
        return lhs < rhs

    count = 0

    for y in range(-17500, 20501):
        Y = y - cy
        Y2 = Y * Y

        # Approximate isoptic boundary using float
        p_val = 2.0*Y2 - 2.0*A2 - 6.0*B2
        q_val = float(Y2)*Y2 - (6.0*A2+2.0*B2)*Y2 + float(A2)*A2 + 6.0*A2*B2 + float(B2)*B2
        disc = p_val*p_val - 4*q_val
        if disc < 0:
            continue
        iso_hi_sq = (-p_val + math.sqrt(disc)) / 2.0
        if iso_hi_sq <= 0:
            continue
        iso_hi = math.sqrt(iso_hi_sq)

        ell_rhs_num = A2 * (B2 - Y2)
        if ell_rhs_num <= 0:
            ell_lo = 0.0
        else:
            ell_lo = math.sqrt(ell_rhs_num / B2)

        if ell_lo >= iso_hi + 2:
            continue

        # Positive X side with exact boundary checks
        x_lo = math.ceil(cx + ell_lo)
        x_hi = math.floor(cx + iso_hi) + 1

        while x_lo <= x_hi and not is_valid(x_lo - cx, Y2):
            x_lo += 1
        while x_hi >= x_lo and not is_valid(x_hi - cx, Y2):
            x_hi -= 1
        while x_hi + 1 <= cx + int(iso_hi) + 3 and is_valid(x_hi + 1 - cx, Y2):
            x_hi += 1

        if x_lo <= x_hi:
            count += x_hi - x_lo + 1

        # Negative X side (use cx-1 as upper bound when ell_lo=0 to avoid
        # double-counting the x=cx column already covered by positive side)
        x_lo2 = math.ceil(cx - iso_hi) - 1
        x_hi2 = (cx - 1) if ell_lo < 0.5 else math.floor(cx - ell_lo)

        while x_lo2 <= x_hi2 and not is_valid(x_lo2 - cx, Y2):
            x_lo2 += 1
        while x_hi2 >= x_lo2 and not is_valid(x_hi2 - cx, Y2):
            x_hi2 -= 1
        while x_lo2 - 1 >= cx - int(iso_hi) - 3 and is_valid(x_lo2 - 1 - cx, Y2):
            x_lo2 -= 1

        if x_lo2 <= x_hi2:
            count += x_hi2 - x_lo2 + 1

    print(count)

solve()