Problem 513: Integral Median

Project Euler Problem 513 Solution

EulerSolve provides an optimized solution for Project Euler Problem 513, Integral Median, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The task is to count primitive integral-median configurations whose size parameter is at most $n$. A naive scan over geometric candidates would be far too slow. The implementations therefore replace the geometry by an arithmetic parametrization, count the resulting integer lattice points inside rational trapezoids, and only at the end remove the nonprimitive configurations produced by odd rescaling. Mathematical Approach Write $G(n)$ for the number of admissible parameter tuples before primitive filtering, and $C(n)$ for the primitive count required by the problem. The method has two layers: first evaluate $G(n)$ exactly by lattice-point counting, then recover $C(n)$ by subtracting odd dilations of smaller primitive objects. Step 1: Reparametrize the Problem Arithmetically In the parametrization used by the implementations, each admissible configuration is encoded by two positive integer pairs $(\alpha,\beta)$ and $(p,q)$. The admissibility conditions become $$\alpha \gt \sqrt{3}\,\beta,\qquad \beta p \le \alpha q,\qquad (\alpha+\beta)p \gt (\alpha+3\beta)q,\qquad \alpha p-\beta q \le n.$$ There are also parity restrictions. The auxiliary pair $(p,q)$ may lie only in the residue classes $$ (p,q)\equiv (0,1),\ (1,0),\ (1,1)\pmod 2, $$ because the all-even class would represent a nonprimitive encoding....

Detailed mathematical approach

Problem Summary

The task is to count primitive integral-median configurations whose size parameter is at most $n$. A naive scan over geometric candidates would be far too slow. The implementations therefore replace the geometry by an arithmetic parametrization, count the resulting integer lattice points inside rational trapezoids, and only at the end remove the nonprimitive configurations produced by odd rescaling.

Mathematical Approach

Write $G(n)$ for the number of admissible parameter tuples before primitive filtering, and $C(n)$ for the primitive count required by the problem. The method has two layers: first evaluate $G(n)$ exactly by lattice-point counting, then recover $C(n)$ by subtracting odd dilations of smaller primitive objects.

Step 1: Reparametrize the Problem Arithmetically

In the parametrization used by the implementations, each admissible configuration is encoded by two positive integer pairs $(\alpha,\beta)$ and $(p,q)$. The admissibility conditions become

$$\alpha \gt \sqrt{3}\,\beta,\qquad \beta p \le \alpha q,\qquad (\alpha+\beta)p \gt (\alpha+3\beta)q,\qquad \alpha p-\beta q \le n.$$

There are also parity restrictions. The auxiliary pair $(p,q)$ may lie only in the residue classes

$$ (p,q)\equiv (0,1),\ (1,0),\ (1,1)\pmod 2, $$

because the all-even class would represent a nonprimitive encoding. The generating pair $(\alpha,\beta)$ must then follow the compatible even/odd progression for the chosen class, so the implementation advances through these parameters in steps of $2$ instead of testing every integer.

Step 2: For Fixed $(\alpha,\beta)$, the Remaining Points Form a Trapezoid

Fix the generating pair $(\alpha,\beta)$. Rearranging the inequalities gives the allowed range for $p$ as a function of $q$:

$$\frac{\alpha+3\beta}{\alpha+\beta}q \lt p \le \min\left(\frac{\alpha}{\beta}q,\frac{\beta q+n}{\alpha}\right).$$

So for each positive integer $q$, the admissible values of $p$ lie above one open line and below the smaller of two closed lines. The two upper bounds cross at

$$q_{\mathrm{sw}}=\left\lfloor\frac{\beta n}{\alpha^2-\beta^2}\right\rfloor,$$

and the whole region disappears after

$$q_{\max}=\left\lfloor\frac{n(\alpha+\beta)}{\alpha^2+2\alpha\beta-\beta^2}\right\rfloor.$$

Therefore the contribution for fixed $(\alpha,\beta)$ is counted as

$$\text{two closed trapezoids} - \text{one open trapezoid}.$$

This is the geometric core of the solution.

Step 3: Convert Each Parity Class to an Ordinary Floor-Sum

For a chosen residue class, write

$$q=2q'+\varepsilon_q,\qquad p=2p'+\varepsilon_p,\qquad \varepsilon_p,\varepsilon_q\in\{0,1\}.$$

After this affine change of variables, each parity-restricted trapezoid turns into an ordinary lattice count with doubled denominator. Closed edges are counted with a usual floor, while open edges are handled by shifting the numerator down by $1$. In that way every residue class is reduced to the same numerical primitive: summing values of a floor of a linear expression.

Step 4: Evaluate the Trapezoids by Euclidean Floor-Sum Reduction

For integers $A,B,M$ and an interval $L \lt x \le U$, define

$$T(A,B,M;L,U)=\sum_{x=L+1}^{U}\left\lfloor\frac{Ax+B}{M}\right\rfloor.$$

This counts lattice points under the rational line $y=(Ax+B)/M$. For a strict upper boundary one instead uses $\lfloor (Ax+B-1)/M\rfloor$. The implementations evaluate these sums exactly by Euclidean-style reduction: first remove the whole-number parts of $A/M$ and $B/M$, then swap slope and denominator on the reduced problem. Once the interval becomes tiny, direct summation finishes the job. Because every admissible region from Step 2 is trapezoidal, $G(n)$ is obtained entirely from a small number of such exact floor-sums.

Step 5: Use a Hyperbola Split for the Large-Parameter Range

Iterating over every possible generating pair would waste work when $\alpha$ is large. The implementations therefore split at

$$R=\left\lfloor\sqrt{\frac{3n}{2}}\right\rfloor.$$

For $\alpha \lt R$, the code fixes $(\alpha,\beta)$ and counts $(p,q)$ directly. For the complementary range it reverses the viewpoint: fix $(p,q)$ and count admissible $(\alpha,\beta)$. Transposing the inequalities from Step 1 gives

$$\beta \le \min\left(\frac{q}{p}\alpha,\frac{p-q}{3q-p}\alpha\right),\qquad \beta \gt \frac{p\alpha-n}{q}.$$

The active upper slope changes exactly when

$$p^2=3q^2,$$

which explains the two-case split in the second major loop. This is a classic hyperbola-style decomposition: one pass is efficient for small generators, the transposed pass is efficient for large generators.

Step 6: Remove Nonprimitive Objects by Odd-Scale Recursion

After the parity normalization, every nonprimitive configuration is an odd dilation of a primitive one. Therefore

$$C(n)=G(n)-\sum_{\substack{d\ge 3\\ d\text{ odd}}} C\!\left(\left\lfloor\frac{n}{d}\right\rfloor\right).$$

The same quotient $\lfloor n/d\rfloor$ appears for many consecutive odd values of $d$, so the implementations group equal quotients into blocks and subtract one cached subproblem multiplied by the number of odd dilations in that block. This is why the primitive filtering stays fast.

Worked Example: Quotient Grouping at $n=10$

For $n=10$, the odd dilations $d\ge 3$ produce

$$\left\lfloor\frac{10}{3}\right\rfloor=3,\qquad \left\lfloor\frac{10}{5}\right\rfloor=2,\qquad \left\lfloor\frac{10}{7}\right\rfloor=\left\lfloor\frac{10}{9}\right\rfloor=1.$$

So the recursion becomes

$$C(10)=G(10)-C(3)-C(2)-2C(1).$$

Two different odd scales collapse to the same subproblem $C(1)$, which is exactly the optimization exploited by quotient grouping. The checkpoint used by the implementations is $C(10)=3$.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. They first compute $R=\lfloor\sqrt{3n/2}\rfloor$ and then evaluate the raw count $G(n)$ across the three admissible parity classes. The first pass iterates over the small generating range and adds the two closed trapezoids minus the open trapezoid from Step 2. The second pass handles the complementary range by fixing the auxiliary pair and counting the transposed trapezoids from Step 5. In both passes, parity is enforced by stepping through the correct congruence classes and by applying affine parity shifts inside the floor-sum evaluator.

After $G(n)$ is known, the implementation computes the primitive total $C(n)$ with memoized recursion. Cached values avoid repeated work, and the subtraction over odd dilations is grouped by equal quotients $\lfloor n/d\rfloor$, so the same smaller argument is solved only once. The three language versions are mathematically identical.

Complexity Analysis

The geometric part of one distinct subproblem is organized as a hyperbola decomposition. The direct and transposed passes together create about $O(n)$ trapezoid evaluations for an argument of size $n$, and each trapezoid evaluation is reduced by Euclidean floor-sum transformations, giving roughly logarithmic work per call. The primitive-filter recursion is much cheaper than subtracting every odd scale independently because memoization and quotient grouping collapse many scales to the same smaller argument. The overall method is comfortably subquadratic in practice and far faster than direct enumeration of candidate triangles.

Footnotes and References

Problem page: https://projecteuler.net/problem=513
Triangle median: Wikipedia — Median (geometry)
Apollonius's theorem: Wikipedia — Apollonius's theorem
Lattice-point counting: Wikipedia — Lattice point
Möbius inversion and primitive counting: Wikipedia — Möbius inversion formula

Problem 513 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <unordered_map>

namespace {

using i64 = long long;

constexpr bool kClosed = true;
constexpr bool kOpen = false;

struct Options {
    int n = 100000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) return false;
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) return false;

    unsigned long long parsed = 0ULL;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') return false;
        const unsigned long long d = static_cast<unsigned long long>(ch - '0');
        if (parsed > (std::numeric_limits<unsigned long long>::max() - d) / 10ULL) return false;
        parsed = parsed * 10ULL + d;
    }
    if (parsed > static_cast<unsigned long long>(std::numeric_limits<int>::max())) return false;
    value = static_cast<int>(parsed);
    return true;
}

bool parse_arguments(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 (parse_int_after_prefix(arg, "--n=", options.n)) continue;
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1;
}

inline i64 floor_div(i64 a, i64 b) {
    assert(b > 0);
    if (a >= 0) return a / b;
    return -((-a + b - 1) / b);
}

i64 isqrt_floor(i64 n) {
    i64 r = static_cast<i64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1) * (r + 1) <= n) ++r;
    while (r * r > n) --r;
    return r;
}

i64 points_in_trapezoid(i64 slope,
                        i64 intercept,
                        i64 denominator,
                        i64 lower_domain,
                        i64 upper_domain,
                        bool boundary) {
    i64 result = 0;
    while (true) {
        assert(denominator > 0);
        if (std::llabs(upper_domain - lower_domain) <= 8) {
            i64 s = 0;
            const i64 adjustment = boundary ? 0 : 1;
            if (upper_domain > lower_domain) {
                for (i64 x = lower_domain + 1; x <= upper_domain; ++x) {
                    s += floor_div(slope * x + intercept - adjustment, denominator);
                }
            } else {
                for (i64 x = upper_domain + 1; x <= lower_domain; ++x) {
                    s += floor_div(slope * x + intercept - adjustment, denominator);
                }
                s = -s;
            }
            result += s;
            break;
        }

        result += (upper_domain - lower_domain) * floor_div(intercept, denominator);
        intercept = ((intercept % denominator) + denominator) % denominator;
        result += ((upper_domain - lower_domain) * (upper_domain + lower_domain + 1) / 2) *
                  floor_div(slope, denominator);
        slope = ((slope % denominator) + denominator) % denominator;

        if (slope != 0) {
            const i64 upper_value = floor_div(slope * upper_domain + intercept, denominator);
            const i64 lower_value = floor_div(slope * lower_domain + intercept, denominator);
            result += upper_domain * upper_value - lower_domain * lower_value;
            lower_domain = upper_value;
            upper_domain = lower_value;
            std::swap(slope, denominator);
            intercept = -intercept;
            boundary = !boundary;
        } else {
            if (intercept == 0 && !boundary) result -= upper_domain - lower_domain;
            break;
        }
    }
    return result;
}

i64 points_in_trapezoid_mod2(i64 slope,
                             i64 intercept,
                             i64 denominator,
                             i64 lower_domain,
                             i64 upper_domain,
                             bool boundary,
                             i64 x_residue,
                             i64 y_residue) {
    if ((y_residue & 1LL) != 0) intercept += denominator;
    if ((x_residue & 1LL) != 0) {
        intercept -= slope;
        ++lower_domain;
        ++upper_domain;
    }
    return points_in_trapezoid(2 * slope, intercept, 2 * denominator,
                               floor_div(lower_domain, 2), floor_div(upper_domain, 2), boundary);
}

i64 f_value(i64 n) {
    i64 result = 0;
    const i64 three_halves_n = n + n / 2;
    const i64 root = isqrt_floor(three_halves_n);
    const int ij[3][2] = {{0, 1}, {1, 0}, {1, 1}};

    for (int tc = 0; tc < 3; ++tc) {
        const i64 i = ij[tc][0];
        const i64 j = ij[tc][1];

        i64 max_t = 1;
        for (i64 s = 2; s < root; ++s) {
            if (3 * (max_t + 1) * (max_t + 1) <= s * s) ++max_t;
            if (i == j || (s & 1LL) == 0) {
                for (i64 t = ((s - 1) & 1LL) + 1; t <= max_t; t += 2) {
                    const i64 v_mid = t * n / ((s - t) * (s + t));
                    const i64 v_max = n * (s + t) / (s * s + 2 * s * t - t * t);
                    result += points_in_trapezoid_mod2(s, 0, t, 0, v_mid, kClosed, j, i);
                    result += points_in_trapezoid_mod2(t, n, s, v_mid, v_max, kClosed, j, i);
                    result -= points_in_trapezoid_mod2(s + 3 * t, 0, s + t, 0, v_max, kOpen, j, i);
                }
            }
        }

        const i64 u_max = three_halves_n / root;
        for (i64 u = 1 + ((i + 1) & 1LL); u <= u_max; u += 2) {
            const i64 v_max_outer = std::min(n / 2, u - 1);
            for (i64 v = 1 + ((j + 1) & 1LL); v <= v_max_outer; v += 2) {
                const i64 mid_s = (v + n) / u;
                const int residue_count = (i == j ? 2 : 1);
                for (int sr = 0; sr < residue_count; ++sr) {
                    const i64 s_residue = (i == j ? sr : 0);

                    i64 min_s = 0, max_s = -1;
                    i64 slope0 = 0, slope1 = 1;
                    if (u * u < 3 * v * v) {
                        min_s = root;
                        max_s = n * (3 * v - u) / (2 * u * v + v * v - u * u);
                        slope0 = u - v;
                        slope1 = 3 * v - u;
                    } else {
                        min_s = std::max(root, floor_div(u + v - 1, v));
                        max_s = n * u / ((u - v) * (u + v));
                        slope0 = v;
                        slope1 = u;
                    }

                    if (max_s >= min_s) {
                        result += points_in_trapezoid_mod2(slope0, 0, slope1,
                                                           min_s - 1, max_s, kClosed,
                                                           s_residue, s_residue);
                        if (mid_s < max_s) {
                            result -= points_in_trapezoid_mod2(u, -n, v,
                                                               std::max(mid_s, min_s - 1), max_s,
                                                               kOpen, s_residue, s_residue);
                        }
                    }
                }
            }
        }
    }
    return result;
}

std::unordered_map<i64, i64> memo_F;

i64 F(i64 n) {
    if (n <= 0) return 0;
    const auto it = memo_F.find(n);
    if (it != memo_F.end()) return it->second;

    i64 result = f_value(n);
    i64 k = 3;
    i64 n_over_k = n / k;
    while (k <= n_over_k) {
        result -= F(n_over_k);
        k += 2;
        n_over_k = n / k;
    }

    i64 min_k = n / (n_over_k + 1);
    while (n_over_k) {
        const i64 max_k = n / n_over_k;
        const i64 left = (min_k + 1) + (min_k & 1LL);
        const i64 right = max_k - ((max_k + 1) & 1LL);
        i64 count = 0;
        if (right >= left) count = (right - left) / 2 + 1;
        result -= F(n_over_k) * count;
        --n_over_k;
        min_k = max_k;
    }

    memo_F.emplace(n, result);
    return result;
}

bool run_checkpoints() {
    if (F(10) != 3) {
        std::cerr << "Validation failed: F(10)\n";
        return false;
    }
    if (F(50) != 165) {
        std::cerr << "Validation failed: F(50)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    memo_F.reserve(1 << 15);
    if (options.run_checkpoints && !run_checkpoints()) {
        return 1;
    }

    std::cout << F(options.n) << '\n';
    return 0;
}

Python

import sys
import math

sys.setrecursionlimit(2000)

memo_F = {}

def floor_div(a, b):
    return a // b

def isqrt_floor(n):
    r = int(math.sqrt(n))
    while (r + 1) ** 2 <= n:
        r += 1
    while r * r > n:
        r -= 1
    return r

def points_in_trapezoid(slope, intercept, denominator, lower_domain, upper_domain, boundary):
    result = 0
    while True:
        if abs(upper_domain - lower_domain) <= 8:
            s = 0
            adjustment = 0 if boundary else 1
            if upper_domain > lower_domain:
                for x in range(lower_domain + 1, upper_domain + 1):
                    s += floor_div(slope * x + intercept - adjustment, denominator)
            else:
                for x in range(upper_domain + 1, lower_domain + 1):
                    s += floor_div(slope * x + intercept - adjustment, denominator)
                s = -s
            result += s
            break

        result += (upper_domain - lower_domain) * floor_div(intercept, denominator)
        intercept %= denominator
        
        result += ((upper_domain - lower_domain) * (upper_domain + lower_domain + 1) // 2) * floor_div(slope, denominator)
        slope %= denominator

        if slope != 0:
            upper_value = floor_div(slope * upper_domain + intercept, denominator)
            lower_value = floor_div(slope * lower_domain + intercept, denominator)
            result += upper_domain * upper_value - lower_domain * lower_value
            lower_domain, upper_domain = upper_value, lower_value
            slope, denominator = denominator, slope
            intercept = -intercept
            boundary = not boundary
        else:
            if intercept == 0 and not boundary:
                result -= upper_domain - lower_domain
            break
            
    return result

def points_in_trapezoid_mod2(slope, intercept, denominator, lower_domain, upper_domain, boundary, x_residue, y_residue):
    if (y_residue & 1) != 0:
        intercept += denominator
    if (x_residue & 1) != 0:
        intercept -= slope
        lower_domain += 1
        upper_domain += 1
        
    return points_in_trapezoid(2 * slope, intercept, 2 * denominator,
                               floor_div(lower_domain, 2), floor_div(upper_domain, 2), boundary)

def f_value(n):
    result = 0
    three_halves_n = n + n // 2
    root = isqrt_floor(three_halves_n)
    ij = [(0, 1), (1, 0), (1, 1)]
    
    for i, j in ij:
        max_t = 1
        for s in range(2, root):
            if 3 * (max_t + 1) ** 2 <= s * s:
                max_t += 1
            if i == j or (s & 1) == 0:
                start_t = ((s - 1) & 1) + 1
                for t in range(start_t, max_t + 1, 2):
                    v_mid = t * n // ((s - t) * (s + t))
                    v_max = n * (s + t) // (s * s + 2 * s * t - t * t)
                    result += points_in_trapezoid_mod2(s, 0, t, 0, v_mid, True, j, i)
                    result += points_in_trapezoid_mod2(t, n, s, v_mid, v_max, True, j, i)
                    result -= points_in_trapezoid_mod2(s + 3 * t, 0, s + t, 0, v_max, False, j, i)
                    
        u_max = three_halves_n // root
        start_u = 1 + ((i + 1) & 1)
        for u in range(start_u, u_max + 1, 2):
            v_max_outer = min(n // 2, u - 1)
            start_v = 1 + ((j + 1) & 1)
            for v in range(start_v, v_max_outer + 1, 2):
                mid_s = (v + n) // u
                residue_count = 2 if i == j else 1
                
                for sr in range(residue_count):
                    s_residue = sr if i == j else 0
                    
                    if u * u < 3 * v * v:
                        min_s = root
                        max_s = n * (3 * v - u) // (2 * u * v + v * v - u * u)
                        slope0 = u - v
                        slope1 = 3 * v - u
                    else:
                        min_s = max(root, floor_div(u + v - 1, v))
                        max_s = n * u // ((u - v) * (u + v))
                        slope0 = v
                        slope1 = u
                        
                    if max_s >= min_s:
                        result += points_in_trapezoid_mod2(slope0, 0, slope1,
                                                           min_s - 1, max_s, True,
                                                           s_residue, s_residue)
                        if mid_s < max_s:
                            result -= points_in_trapezoid_mod2(u, -n, v,
                                                               max(mid_s, min_s - 1), max_s,
                                                               False, s_residue, s_residue)
    return result

def F(n):
    if n <= 0: return 0
    if n in memo_F: return memo_F[n]
    
    result = f_value(n)
    k = 3
    n_over_k = n // k
    while k <= n_over_k:
        result -= F(n_over_k)
        k += 2
        n_over_k = n // k
        
    min_k = n // (n_over_k + 1) if n_over_k + 1 > 0 else n
        
    while n_over_k > 0:
        max_k = n // n_over_k
        left = (min_k + 1) + (min_k & 1)
        right = max_k - ((max_k + 1) & 1)
        count = 0
        if right >= left:
            count = (right - left) // 2 + 1
        result -= F(n_over_k) * count
        n_over_k -= 1
        min_k = max_k
        
    memo_F[n] = result
    return result

def solve():
    n = 100000
    ans = F(n)
    return str(ans)

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

Java

import java.util.HashMap;
import java.util.Map;

public class Euler513 {

    static Map<Long, Long> memoF = new HashMap<>();

    static long floorDiv(long a, long b) {
        return Math.floorDiv(a, b);
    }

    static long isqrtFloor(long n) {
        long r = (long) Math.sqrt((double) n);
        while ((r + 1) * (r + 1) <= n)
            r++;
        while (r * r > n)
            r--;
        return r;
    }

    static long pointsInTrapezoid(long slope, long intercept, long denominator, long lowerDomain, long upperDomain,
            boolean boundary) {
        long result = 0;
        while (true) {
            if (Math.abs(upperDomain - lowerDomain) <= 8) {
                long s = 0;
                long adjustment = boundary ? 0 : 1;
                if (upperDomain > lowerDomain) {
                    for (long x = lowerDomain + 1; x <= upperDomain; x++) {
                        s += floorDiv(slope * x + intercept - adjustment, denominator);
                    }
                } else {
                    for (long x = upperDomain + 1; x <= lowerDomain; x++) {
                        s += floorDiv(slope * x + intercept - adjustment, denominator);
                    }
                    s = -s;
                }
                result += s;
                break;
            }

            result += (upperDomain - lowerDomain) * floorDiv(intercept, denominator);
            intercept = ((intercept % denominator) + denominator) % denominator;

            result += ((upperDomain - lowerDomain) * (upperDomain + lowerDomain + 1) / 2)
                    * floorDiv(slope, denominator);
            slope = ((slope % denominator) + denominator) % denominator;

            if (slope != 0) {
                long upperValue = floorDiv(slope * upperDomain + intercept, denominator);
                long lowerValue = floorDiv(slope * lowerDomain + intercept, denominator);
                result += upperDomain * upperValue - lowerDomain * lowerValue;
                lowerDomain = upperValue;
                upperDomain = lowerValue;
                long tmp = slope;
                slope = denominator;
                denominator = tmp;
                intercept = -intercept;
                boundary = !boundary;
            } else {
                if (intercept == 0 && !boundary)
                    result -= upperDomain - lowerDomain;
                break;
            }
        }
        return result;
    }

    static long pointsInTrapezoidMod2(long slope, long intercept, long denominator, long lowerDomain, long upperDomain,
            boolean boundary, long xResidue, long yResidue) {
        if ((yResidue & 1) != 0)
            intercept += denominator;
        if ((xResidue & 1) != 0) {
            intercept -= slope;
            lowerDomain++;
            upperDomain++;
        }
        return pointsInTrapezoid(2 * slope, intercept, 2 * denominator, floorDiv(lowerDomain, 2),
                floorDiv(upperDomain, 2), boundary);
    }

    static long fValue(long n) {
        long result = 0;
        long threeHalvesN = n + n / 2;
        long root = isqrtFloor(threeHalvesN);
        int[][] ij = { { 0, 1 }, { 1, 0 }, { 1, 1 } };

        for (int[] pair : ij) {
            long i = pair[0];
            long j = pair[1];

            long maxT = 1;
            for (long s = 2; s < root; s++) {
                if (3 * (maxT + 1) * (maxT + 1) <= s * s)
                    maxT++;
                if (i == j || (s & 1) == 0) {
                    for (long t = ((s - 1) & 1) + 1; t <= maxT; t += 2) {
                        long vMid = t * n / ((s - t) * (s + t));
                        long vMax = n * (s + t) / (s * s + 2 * s * t - t * t);
                        result += pointsInTrapezoidMod2(s, 0, t, 0, vMid, true, j, i);
                        result += pointsInTrapezoidMod2(t, n, s, vMid, vMax, true, j, i);
                        result -= pointsInTrapezoidMod2(s + 3 * t, 0, s + t, 0, vMax, false, j, i);
                    }
                }
            }

            long uMax = threeHalvesN / root;
            for (long u = 1 + ((i + 1) & 1); u <= uMax; u += 2) {
                long vMaxOuter = Math.min(n / 2, u - 1);
                for (long v = 1 + ((j + 1) & 1); v <= vMaxOuter; v += 2) {
                    long midS = (v + n) / u;
                    int residueCount = (i == j ? 2 : 1);
                    for (int sr = 0; sr < residueCount; sr++) {
                        long sResidue = (i == j ? sr : 0);

                        long minS = 0, maxS = -1;
                        long slope0 = 0, slope1 = 1;
                        if (u * u < 3 * v * v) {
                            minS = root;
                            maxS = n * (3 * v - u) / (2 * u * v + v * v - u * u);
                            slope0 = u - v;
                            slope1 = 3 * v - u;
                        } else {
                            minS = Math.max(root, floorDiv(u + v - 1, v));
                            maxS = n * u / ((u - v) * (u + v));
                            slope0 = v;
                            slope1 = u;
                        }

                        if (maxS >= minS) {
                            result += pointsInTrapezoidMod2(slope0, 0, slope1, minS - 1, maxS, true, sResidue,
                                    sResidue);
                            if (midS < maxS) {
                                result -= pointsInTrapezoidMod2(u, -n, v, Math.max(midS, minS - 1), maxS, false,
                                        sResidue, sResidue);
                            }
                        }
                    }
                }
            }
        }
        return result;
    }

    static long F(long n) {
        if (n <= 0)
            return 0;
        if (memoF.containsKey(n))
            return memoF.get(n);

        long result = fValue(n);
        long k = 3;
        long nOverK = n / k;
        while (k <= nOverK) {
            result -= F(nOverK);
            k += 2;
            nOverK = n / k;
        }

        long minK = (nOverK + 1 > 0) ? n / (nOverK + 1) : n;
        while (nOverK > 0) {
            long maxK = n / nOverK;
            long left = (minK + 1) + (minK & 1);
            long right = maxK - ((maxK + 1) & 1);
            long count = 0;
            if (right >= left)
                count = (right - left) / 2 + 1;
            result -= F(nOverK) * count;
            nOverK--;
            minK = maxK;
        }

        memoF.put(n, result);
        return result;
    }

    public static void main(String[] args) {
        long n = 100000;
        System.out.println(F(n));
    }
}