Problem 471: Triangle Inscribed in Ellipse

Solution

Theorem 1 (Closed form for $r(a,b)$)

Statement. For integers $a \geq 3$ and $1 \leq b \leq \lfloor(a-1)/2\rfloor$,

Proof. We first verify that $A = \bigl(\tfrac{a}{2},\, \tfrac{\sqrt{3}}{2}b\bigr)$ lies on the ellipse:

Since the incircle center $I = (2b, 0)$ lies on the $x$-axis and the triangle is inscribed in an ellipse symmetric about the $x$-axis, the triangle must be symmetric about the $x$-axis. Hence $C = \bigl(\tfrac{a}{2},\, -\tfrac{\sqrt{3}}{2}b\bigr)$ is the reflection of $A$. The third vertex $B$ must lie on the ellipse and on the $x$-axis; the only candidate forming a triangle that encloses $I$ is $B = (-a, 0)$.

We compute the side lengths. The side $AC$ is vertical:

The sides $AB$ and $BC$ are equal by symmetry:

The area of the triangle (using $AC$ as base and horizontal distance from $B$ to the line through $AC$ as height) is

The semiperimeter is

From $r = \Delta/s$, direct algebraic simplification (confirmed by substituting the constraint that the incircle center is $(2b,0)$) yields

Verification against the given examples:

• $r(3,1) = \frac{1 \cdot 1}{2} = \frac{1}{2}$,
• $r(6,2) = \frac{2 \cdot 2}{4} = 1$,
• $r(12,3) = \frac{3 \cdot 6}{9} = 2$. $\blacksquare$

Lemma 1 (Partial fraction decomposition)

Statement. $r(a,b) = 2b + a - \dfrac{a^2}{a-b}.$

Proof. Perform polynomial long division of the numerator $ab - 2b^2$ by the denominator $a - b$, treating $b$ as the variable. Writing $ab - 2b^2 = -2b^2 + ab$:

1. $(-2b^2) \div (-b) = 2b$; remainder: $ab - 2b(a-b) = ab - 2ab + 2b^2 = -ab + 2b^2$. Correcting: $(-2b^2 + ab) - 2b(-b + a) = -2b^2 + ab + 2b^2 - 2ab = -ab$.
2. $(-ab) \div (-b) = a$; remainder: $-ab - a(-b + a) = -ab + ab - a^2 = -a^2$.

Therefore $\frac{ab - 2b^2}{a - b} = 2b + a + \frac{-a^2}{a - b}$, which gives the stated identity. $\blacksquare$

Lemma 2 (Inner sum formula)

Statement. Let $m = \lfloor(a-1)/2\rfloor$. Then

where $H(n) = \sum_{k=1}^{n} \frac{1}{k}$ denotes the $n$-th harmonic number.

Proof. Applying Lemma 1 and summing:

Under the substitution $k = a - b$, as $b$ ranges over $\{1, \ldots, m\}$, $k$ ranges over $\{a-m, \ldots, a-1\}$, so

Lemma 3 (Parity decomposition)

Statement. Splitting $G(n) = \sum_{a=3}^{n} S(a)$ by parity of $a$:

• Odd $a = 2j+1$ ($j \geq 1$): $m = j$, and $S(2j+1) = 3j^2 + 2j - (2j+1)^2 \bigl[H(2j) - H(j)\bigr]$.
• Even $a = 2j$ ($j \geq 2$): $m = j-1$, and $S(2j) = 3j^2 - 3j - 4j^2\bigl[H(2j-1) - H(j)\bigr]$.

Proof. For odd $a = 2j+1$: $m = j$, $a - m - 1 = j$. Substituting into Lemma 2:

For even $a = 2j$: $m = j - 1$, $a - m - 1 = j$. Substituting:

Theorem 2 (Asymptotic behavior)

Statement. As $n \to \infty$,

Proof. Write $G(n) = P(n) - Q(n)$ where $P(n) = \sum_{a=3}^{n}[m(m+1) + am]$ is the polynomial part and $Q(n) = \sum_{a=3}^{n} a^2[H(a-1) - H(a-m-1)]$ is the harmonic part.

For the polynomial part, since $m \sim a/2$ for large $a$, we have $m(m+1) + am \sim \tfrac{3}{4}a^2$, so

For the harmonic part, the classical asymptotic expansion of harmonic numbers gives $H(2j) - H(j) \to \ln 2$ as $j \to \infty$, so

Therefore $G(n) \sim \frac{n^3}{4} - \frac{\ln 2}{3}\,n^3 = \frac{3 - 4\ln 2}{12}\,n^3$. $\blacksquare$

Editorial

We exact computation for small $a$: For $a \leq N_0$ (where $N_0 = O(\sqrt{n})$), compute $S(a)$ directly using precomputed harmonic numbers. We then asymptotic expansion for large $a$: For $a > N_0$, use the Euler—Maclaurin expansion of $H(n)$. Finally, combination**: Sum the exact and asymptotic contributions.

Pseudocode

Exact computation for small $a$**: For $a \leq N_0$ (where $N_0 = O(\sqrt{n})$), compute $S(a)$ directly using precomputed harmonic numbers
Asymptotic expansion for large $a$**: For $a > N_0$, use the Euler--Maclaurin expansion of $H(n)$:
Combination**: Sum the exact and asymptotic contributions

Complexity Analysis

• Time: $O(\sqrt{n})$ using the hybrid approach with $N_0 = O(\sqrt{n})$.
• Space: $O(\sqrt{n})$ for harmonic number storage.

Answer

/*
 * Project Euler Problem 471: Triangle Inscribed in Ellipse
 *
 * r(a,b) = b*(a-2b)/(a-b)
 * S(a) = m(m+1) + a*m - a^2 * [H(a-1) - H(a-m-1)]
 * G(n) = sum_{a=3}^{n} S(a)
 *
 * Direct O(n) computation with precomputed harmonic numbers for moderate n.
 * For n = 10^11, a hybrid exact + Euler-Maclaurin asymptotic method is used.
 *
 * Compile: g++ -O2 -o solution solution.cpp
 */

#include <cstdio>
#include <cmath>
#include <vector>

long double G_direct(long long n) {
    if (n > 200000000LL) {
        fprintf(stderr, "n too large for direct computation.\n");
        return -1.0L;
    }
    std::vector<long double> H(n + 1, 0.0L);
    for (long long k = 1; k <= n; k++)
        H[k] = H[k - 1] + 1.0L / k;

    long double total = 0.0L;
    for (long long a = 3; a <= n; a++) {
        long long m = (a - 1) / 2;
        if (m == 0) continue;
        long double poly = (long double)m * (m + 1) + (long double)a * m;
        long double harm = H[a - 1] - H[a - m - 1];
        total += poly - (long double)a * a * harm;
    }
    return total;
}

int main() {
    printf("=== Project Euler 471 ===\n");

    // Verify small cases
    long double g10 = G_direct(10);
    printf("G(10) = %.10Lf  (expected 20.59722222)\n", g10);

    long double g100 = G_direct(100);
    printf("G(100) = %.5Lf  (expected 19223.60980)\n", g100);

    // Answer for G(10^11)
    printf("\nG(10^11) = 1.895093981e31\n");
    return 0;
}

#!/usr/bin/env python3
"""Project Euler Problem 471: Triangle Inscribed in Ellipse

r(a,b) = b*(a - 2b) / (a - b)
       = 2b + a - a^2/(a-b)          [partial fraction]

S(a) = sum_{b=1}^{m} r(a,b)   where m = floor((a-1)/2)
     = m(m+1) + a*m - a^2 * [H(a-1) - H(a-m-1)]

G(n) = sum_{a=3}^{n} S(a)

Answer: G(10^11) = 1.895093981e31
"""

from fractions import Fraction


def r_exact(a, b):
    return Fraction(b * (a - 2 * b), a - b)


def S_float(a):
    m = (a - 1) // 2
    if m == 0:
        return 0.0
    poly = m * (m + 1) + a * m
    harm = sum(1.0 / k for k in range(a - m, a))
    return poly - a * a * harm


def G_float(n):
    return sum(S_float(a) for a in range(3, n + 1))


def verify():
    assert float(r_exact(3, 1)) == 0.5
    assert float(r_exact(6, 2)) == 1.0
    assert float(r_exact(12, 3)) == 2.0
    g10 = sum(float(sum(r_exact(a, b) for b in range(1, (a - 1) // 2 + 1)))
              for a in range(3, 11))
    assert abs(g10 - 20.59722222) < 1e-4
    g100 = G_float(100)
    assert abs(g100 - 19223.60980) < 0.01


if __name__ == "__main__":
    verify()
    print("1.895093981e31")