Project Euler Problem 525

Solution

Answer: 44.69921807

Let the ellipse be

$$\frac{x^2}{a^2}+\frac{(y-b)^2}{b^2}=1$$

so in body coordinates (centered at the ellipse center) the boundary is

$$r(t)=(a\cos t,; b\sin t), \qquad -\frac{\pi}{2}\le t\le \frac{3\pi}{2}.$$

At $t=-\pi/2$, the bottom point touches the ground.

1. Mathematical analysis

When the ellipse rolls without slipping, the contact point changes continuously around the boundary.

Step 1: Contact geometry

For the point $r(t)$, the tangent vector is

$$r'(t)=(-a\sin t,; b\cos t),$$

whose magnitude is

$$q(t)=\sqrt{a^2\sin^2 t+b^2\cos^2 t}.$$

This is the infinitesimal arc-length factor:

$$ds=q(t),dt.$$

The outward normal of the ellipse at parameter $t$ is proportional to

$$\left(\frac{\cos t}{a},\frac{\sin t}{b}\right).$$

Rolling on a horizontal line means this normal points vertically downward in world coordinates. Thus the ellipse rotation angle $\phi(t)$ satisfies

$$\phi(t) = -\frac{\pi}{2} -\operatorname{atan2} !\left(\frac{\sin t}{b},\frac{\cos t}{a}\right).$$

To make the rotation continuous over one turn, we unwrap this angle so that the ellipse rotates through $2\pi$.

Step 2: Center coordinates

Because there is no slipping, the ground-contact point advances by the boundary arc length:

$$x_{\text{contact}}(t) = \int q(t),dt.$$

If the rotated contact point in world coordinates is

$$R_{\phi(t)}r(t) = (x_r(t),y_r(t)),$$

then the ellipse center is

$$X(t)=x_{\text{contact}}(t)-x_r(t), \qquad Y(t)=-y_r(t).$$

The path length is therefore

$$C(a,b) = \int \sqrt{ \left(\frac{dX}{dt}\right)^2 + \left(\frac{dY}{dt}\right)^2 },dt.$$

Numerical integration over one full turn gives the desired value.

As a check, evaluating $C(2,4)$ reproduces the problem statement:

$$C(2,4)\approx 21.38816906.$$

2. Python implementation

import mpmath as mp

mp.mp.dps = 50  # high precision

def C(a, b):
    # Rotation angle of ellipse so tangent is horizontal
    def phi_raw(t):
        return -mp.pi / 2 - mp.atan2(mp.sin(t) / b,
                                      mp.cos(t) / a)

    # Unwrap angle to make one continuous full turn
    def phi(t):
        p = phi_raw(t)
        if t > -mp.pi / 2:
            while p > 0:
                p -= 2 * mp.pi
        return p

    # Arc-length factor of ellipse boundary
    def q(t):
        return mp.sqrt(a*a * mp.sin(t)**2 +
                       b*b * mp.cos(t)**2)

    # Rotated boundary point in world coordinates
    def x_rel(t):
        p = phi(t)
        x = a * mp.cos(t)
        y = b * mp.sin(t)
        return x * mp.cos(p) - y * mp.sin(p)

    def y_rel(t):
        p = phi(t)
        x = a * mp.cos(t)
        y = b * mp.sin(t)
        return x * mp.sin(p) + y * mp.cos(p)

    # Speed of ellipse center
    def integrand(t):
        dx = q(t) - mp.diff(x_rel, t)
        dy = -mp.diff(y_rel, t)
        return mp.sqrt(dx*dx + dy*dy)

    return mp.quad(integrand,
                   [-mp.pi/2, 3*mp.pi/2])

# Verify sample
print(C(2, 4))

# Required quantity
ans = C(1, 4) + C(3, 4)
print(ans)
print(f"{ans:.8f}")

3. Code walkthrough

phi_raw(t) computes the rotation angle required so the tangent at the contact point is horizontal.
phi(t) unwraps the angle to ensure a continuous full revolution.
q(t) computes the differential arc length $ds/dt$ of the ellipse boundary.
x_rel(t), y_rel(t) compute the contact point location after rotation.
The center position comes from:
advancing by boundary arc length (no slipping),
subtracting the rotated contact-point offset.
integrand(t) computes the instantaneous speed of the center.
Numerical integration over one full traversal of the boundary yields $C(a,b)$.

Sample verification

The program gives:

$$C(2,4)=21.3881690605786\ldots$$

which matches the provided value

$$21.38816906.$$

4. Final verification

We obtain

$$C(1,4)+C(3,4) = 44.69921806756434\ldots$$

Rounded to 8 digits after the decimal point:

Answer: 44.69921807