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 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.

Exterior Points and the First Column Bound

Two real tangents exist only when the point lies strictly outside the ellipse. In centered coordinates this means

$$\frac{x^2}{a^2}+\frac{y^2}{b^2} \gt 1,$$

or equivalently

$$b^2x^2+a^2y^2 \gt a^2b^2.$$

For a fixed nonnegative $x$, this gives the lower edge of the admissible vertical column. If $x^2\gt a^2$, then $y=0$ is already outside the ellipse. If $x^2\le a^2$, then we need

$$y^2 \gt \frac{a^2b^2-b^2x^2}{a^2},$$

so the first valid integer $y$ is the smallest integer strictly above that threshold.

Deriving the Tangent-Angle Formula

Take an exterior point $P=(x_0,y_0)$. A line through $P$ with slope $m$ has equation $y=mx+c$, where $c=y_0-mx_0$. Such a line is tangent to the ellipse exactly when

$$c^2=a^2m^2+b^2.$$

Substituting $c=y_0-mx_0$ gives a quadratic equation for the two tangent slopes:

$$\left(x_0^2-a^2\right)m^2-2x_0y_0m+\left(y_0^2-b^2\right)=0.$$

If the roots are $m_1$ and $m_2$, then

$$m_1+m_2=\frac{2x_0y_0}{x_0^2-a^2},\qquad m_1m_2=\frac{y_0^2-b^2}{x_0^2-a^2}.$$

The angle $\theta$ between the two tangent rays satisfies

$$\tan^2\theta=\frac{(m_1-m_2)^2}{(1+m_1m_2)^2}=\frac{4\left(b^2x_0^2+a^2y_0^2-a^2b^2\right)}{\left(x_0^2+y_0^2-a^2-b^2\right)^2}.$$

This identity is the central geometric formula used in all three implementations.

The Director Circle and the $45^\circ$ Test

The denominator vanishes on

$$x^2+y^2=a^2+b^2,$$

the director circle of the ellipse. On this circle the two tangents are perpendicular. That gives a clean case split.

If a point is outside the ellipse but satisfies $x^2+y^2\le a^2+b^2$, then the tangent angle is at least $90^\circ$, so it automatically exceeds $45^\circ$.

If $x^2+y^2\gt a^2+b^2$, then $0\lt\theta\le 90^\circ$, so the problem condition is equivalent to $\tan\theta\gt 1$. Using the previous formula yields

$$4\left(b^2x^2+a^2y^2-a^2b^2\right)-\left(a^2+b^2-x^2-y^2\right)^2 \gt 0.$$

This polynomial inequality is exactly the test performed by the implementations.

Reducing the Count to One Quadrant

The ellipse, the exterior condition, and the angle inequality all depend only on $x^2$ and $y^2$. Therefore the valid set is symmetric under reflections across both coordinate axes. It is enough to count the centered lattice points with $x\ge 0$ and $y\ge 0$, then correct for axis overcounting.

Counting a Fixed $x$-Column

For one nonnegative integer $x$, the admissible $y$-values form a contiguous interval beginning at the first point outside the ellipse.

The first portion of that interval is automatic: every integer

$$y_{\min}(x)\le y\le \left\lfloor\sqrt{a^2+b^2-x^2}\right\rfloor$$

qualifies whenever it lies between the ellipse and the director circle.

Beyond the director circle, set $u=y^2$ and write

$$k=a^2+b^2-x^2.$$

Then the angle condition becomes

$$-u^2+(4a^2+2k)u+4\left(b^2x^2-a^2b^2\right)-k^2 \gt 0.$$

This is a downward-opening quadratic in $u$, so the valid values are exactly those below its larger root. Hence

$$u_{\max}(x)=\frac{4a^2+2k+\sqrt{(4a^2+2k)^2+4\left(4\left(b^2x^2-a^2b^2\right)-k^2\right)}}{2},$$

and the largest admissible integer $y$ is $\lfloor\sqrt{u_{\max}(x)}\rfloor$, adjusted by a tiny integer correction because the inequality is strict.

Worked Example: Why the Search Is Finite

Take the $x$-axis, so $y=0$. Being outside the ellipse means $x^2\gt a^2$. Set

$$z=x^2-a^2.$$

The $45^\circ$ inequality becomes

$$4b^2z \gt (b^2-z)^2,$$

which simplifies to

$$z^2-6b^2z+b^4 \lt 0.$$

Therefore

$$b^2(3-2\sqrt{2}) \lt z \lt b^2(3+2\sqrt{2}).$$

In particular, every qualifying axis point satisfies $z\lt 6b^2$, so

$$x^2 \lt a^2+6b^2.$$

The same argument on the $y$-axis gives $y^2\lt b^2+6a^2$. This is why the code only scans finitely many columns and rows, using those rounded-up bounds as safe limits.

How the Code Works

Geometric Precomputation

The C++, Python, and Java implementations first reduce the problem data to $a^2$, $b^2$, $a^2+b^2$, and $a^2b^2$. Once those are known, every geometric decision becomes an integer comparison in centered coordinates.

Column-by-Column Counting

For each nonnegative integer $x$, the implementation computes the first integer $y$ outside the ellipse. It then counts every point up to the director circle immediately, because those points satisfy the angle condition automatically. If the column extends beyond the director circle, it solves the quadratic inequality in $u=y^2$ to obtain the last admissible $y$ in that column and adds the remaining points.

No tangent points are ever constructed explicitly. The algorithm works entirely with the derived inequalities, which is why it stays fast and numerically stable.

Exact Integer Arithmetic and Symmetry Repair

The implementations use integer square roots and exact integer comparisons, then make short upward or downward corrections when a strict inequality sits just past the square-root estimate. That avoids off-by-one errors near the boundary.

If $N_{Q1}$ is the number of qualifying centered lattice points with $x\ge 0$ and $y\ge 0$, then the full answer is reconstructed as

$$N=4N_{Q1}-2N_x-2N_y+N_0,$$

where $N_x$ and $N_y$ count the qualifying points on the centered axes and $N_0$ is the origin correction. The C++ and Java implementations parallelize the column sweep, while the Python implementation applies the same mathematics with a serial or process-based split depending on the runtime environment.

Complexity Analysis

Let $X_{\max}$ be the safe horizontal scan bound, which is on the order of $\sqrt{a^2+6b^2}$. The main loop processes each integer $x$ from $0$ to $X_{\max}$, and each column requires only constant-time arithmetic plus a few final boundary corrections. The running time is therefore essentially linear in the scanned width.

Memory usage is $O(1)$ beyond a handful of counters and geometric constants, or $O(T)$ if one includes thread-local accumulators for $T$ workers. The important point is that the algorithm never searches over tangent points on the ellipse; it counts lattice points directly from closed formulas.

Footnotes and References

Problem page: https://projecteuler.net/problem=246
Ellipse: Wikipedia - Ellipse
Director circle: Wikipedia - Director circle
Tangent: Wikipedia - Tangent
Quadratic equation: Wikipedia - Quadratic equation

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());
    }
}