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

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