Kazakhstan FAIO 2024 · Correction IOAI Training

D. Flight Speed 2: Data Lost

Cette page donne toute la méthode. Lisez-la une fois que vous avez cherché.

Métrique: exact coefficients on hidden tests
Évaluée sur: self-generated tests (official tests are hidden)
Notre score: 10000 random checks vs numpy solve

Nous avons résolu l'épreuve nous-mêmes : téléchargé les données, mené les expériences et mesuré le score ci-dessous. L'énoncé et les données appartiennent à leurs organisateurs. Cette correction est rédigée en anglais.

The problem in brief

A ship's distance over time follows a parabola $y = ax^2 + bx + c$ , and three points of the graph were recorded before the sensor failed. Your program reads six numbers and must print $a$ , $b$ and $c$ . The twist is in how the points are encoded: the first and third points are given directly, but the second is described relative to the others, $x_2$ to the left of the first point and $y_2$ above the third one. So the actual second point is $(x_1 - x_2,\ y_3 + y_2)$ .

How you are scored

Hidden tests compare the three printed numbers with the true coefficients (the tolerance is not stated, so we print 10 decimals). The README gives no sample, so we check the program on 10 000 random inputs against a generic linear solver.

The key idea

Three points determine a parabola, through a linear system. Each point $(u_i, v_i)$ gives one equation that is linear in the unknowns:

a u_i^2 + b u_i + c = v_i, \qquad i = 1, 2, 3 .

This is a $3\times3$ Vandermonde system. It has exactly one solution when the three $u_i$ are distinct, because a nonzero polynomial of degree at most 2 cannot vanish at three different points. (If two points shared the same $x$ , there would be no function through them at all.)

Rather than inverting a matrix, we can read the coefficients off divided differences, the discrete cousins of derivatives. Write $P_1, P_2, P_3$ for the three points:

s_{12} = \frac{v_2 - v_1}{u_2 - u_1}, \qquad s_{23} = \frac{v_3 - v_2}{u_3 - u_2},

the slopes of two chords. For a parabola, the slope of the chord between $u$ and $u'$ is exactly $a(u + u') + b$ (the derivative $2ax + b$ taken at the chord's midpoint). Subtracting the two chord slopes eliminates $b$ :

s_{23} - s_{12} = a\,(u_3 - u_1) \quad\Longrightarrow\quad a = \frac{s_{23} - s_{12}}{u_3 - u_1}.

Then $b = s_{12} - a(u_1 + u_2)$ and $c = v_1 - a u_1^2 - b u_1$ . The second divided difference $a$ measures how fast the slope changes: it is half the (constant) second derivative of the trajectory, that is, half the acceleration.

The three points, with P2 rebuilt from its relative description, and the parabola through them

The full solution

x1, y1, x2, y2, x3, y3 = map(float, input().split())
(u1, v1), (u2, v2), (u3, v3) = (x1, y1), (x1 - x2, y3 + y2), (x3, y3)
s12 = (v2 - v1) / (u2 - u1)
s23 = (v3 - v2) / (u3 - u2)
a = (s23 - s12) / (u3 - u1)
b = s12 - a * (u1 + u2)
c = v1 - a * u1 * u1 - b * u1

For example, the input 6 50 4 30 9 40 describes the points $(6, 50)$ , $(2, 70)$ and $(9, 40)$ , and the program prints 0.2380952381 -6.9047619048 82.8571428571.

Results

check	result
10 000 random inputs in $(-1000, 1000)$ : coefficients vs `numpy.linalg.solve`	identical (relative tolerance $10^{-6}$ )
worst relative residual $\lvert a u^2 + bu + c - v\rvert / (1 + \lvert v\rvert)$	$2.2\cdot 10^{-11}$
reading the second point as $(x_2, y_2)$ instead, input `1 2 3 4 5 6`	$(0, 1, 1)$ instead of $(0.524, -2.143, 3.619)$ : wrong

What to remember

"Passes through these points" is a linear system in the coefficients, even though the curve is not a line: the model is linear in its parameters. This is the same fact that makes polynomial regression a linear regression.
$n$ points with distinct $x$ determine a unique polynomial of degree $n-1$ ; divided differences compute it without a matrix inverse.
The second divided difference is a discrete second derivative: for a trajectory, $2a$ is the acceleration.
Rebuild the data exactly as described before modelling: here the whole difficulty is the relative encoding of the second point.

Code complet (solution.py)

"""FAIO 2024 qualification, problem D: Flight Speed 2 -- the parabola through three points.

stdin : x1 y1 x2 y2 x3 y3   (the 2nd point is (x1 - x2, y3 + y2), NOT (x2, y2))
stdout: a b c  with  y = a x^2 + b x + c  through the three points

Method: Newton's divided differences (no matrix inverse needed), exact for a quadratic.
Run: python solution.py < input.txt ;  python solution.py --check ;  python solution.py --figure
"""
import sys


def points(x1, y1, x2, y2, x3, y3):
    return (x1, y1), (x1 - x2, y3 + y2), (x3, y3)


def parabola(p, q, r):
    (u1, v1), (u2, v2), (u3, v3) = p, q, r
    s12 = (v2 - v1) / (u2 - u1)               # slope of the chord P1P2
    s23 = (v3 - v2) / (u3 - u2)               # slope of the chord P2P3
    a = (s23 - s12) / (u3 - u1)               # second divided difference
    b = s12 - a * (u1 + u2)
    c = v1 - a * u1 * u1 - b * u1
    return a, b, c


def fmt(v):
    s = f"{v:.10f}".rstrip("0").rstrip(".")
    return "0" if s in ("-0", "") else s


def solve():
    vals = list(map(float, sys.stdin.read().split()[:6]))
    print(" ".join(fmt(v) for v in parabola(*points(*vals))))


def check():
    import random
    import numpy as np
    rng = random.Random(0)
    worst = 0.0; n = 0
    while n < 10000:
        vals = [rng.uniform(-999, 999) for _ in range(6)]
        pts = points(*vals)
        xs = [p[0] for p in pts]
        if min(abs(xs[i] - xs[j]) for i in range(3) for j in range(i)) < 1e-3:
            continue                           # two points on one vertical: no parabola
        a, b, c = parabola(*pts)
        worst = max(worst, max(abs(a * x * x + b * x + c - y) / (1 + abs(y)) for x, y in pts))
        A = np.vander(xs, 3)
        ref = np.linalg.solve(A, [p[1] for p in pts])
        assert np.allclose([a, b, c], ref, rtol=1e-6, atol=1e-6), (vals, (a, b, c), ref)
        n += 1
    print(f"10000 random cases: matches numpy.linalg.solve; worst relative residual {worst:.2e}")
    # the classic trap: using (x2, y2) as the second point
    vals = [1, 2, 3, 4, 5, 6]
    print("example 1 2 3 4 5 6 ->", parabola(*points(*vals)), "| wrong reading (x2, y2) ->",
          parabola((1, 2), (3, 4), (5, 6)))


def figure():
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    import numpy as np
    vals = (6, 50, 4, 30, 9, 40)
    p, q, r = points(*vals)
    a, b, c = parabola(p, q, r)
    xs = np.linspace(0, 11, 200)
    fig, ax = plt.subplots(figsize=(8, 4.6))
    ax.plot(xs, a * xs ** 2 + b * xs + c, color="#2171b5", lw=2, label=f"y = {a:.3g} x² + {b:.3g} x + {c:.3g}")
    for (x, y), nm in zip((p, q, r), ("P1 = (x1, y1)", "P2 = (x1 - x2, y3 + y2)", "P3 = (x3, y3)")):
        ax.plot(x, y, "o", color="crimson", ms=8); ax.annotate(nm, (x, y), xytext=(6, 6), textcoords="offset points")
    ax.annotate("", xy=(q[0], p[1] + 4), xytext=(p[0], p[1] + 4), arrowprops=dict(arrowstyle="->"))
    ax.text((p[0] + q[0]) / 2, p[1] + 5, "x2 to the left", ha="center", fontsize=9)
    ax.annotate("", xy=(r[0] + 0.4, q[1]), xytext=(r[0] + 0.4, r[1]), arrowprops=dict(arrowstyle="->"))
    ax.text(r[0] + 0.55, (q[1] + r[1]) / 2, "y2 above P3", fontsize=9)
    ax.set_xlabel("time x"); ax.set_ylabel("distance y"); ax.legend(loc="lower right")
    ax.set_title("Example input 6 50 4 30 9 40: rebuild P2, then fit the parabola")
    fig.tight_layout(); fig.savefig("figures/parabola.png", dpi=150, facecolor="white")


if __name__ == "__main__":
    if "--check" in sys.argv:
        check()
    elif "--figure" in sys.argv:
        figure()
    else:
        solve()