Problem 460: An Ant on the Move

Project Euler Problem 460 Solution

EulerSolve provides an optimized solution for Project Euler Problem 460, An Ant on the Move, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The ant starts at $A(0,1)$ and must reach $B(d,1)$, where $d\in 2\mathbb{Z}_{>0}$. Every move is a straight segment between lattice points, always going forward in $x$. If a segment starts at height $y_0$ and ends at height $y_1$, then its speed is $$v=\begin{cases} y_0, & y_0=y_1,\\ \dfrac{y_1-y_0}{\ln y_1-\ln y_0}, & y_0\ne y_1. \end{cases}$$ The objective is to minimize the total travel time; that minimum is denoted by $F(d)$. The published checkpoints are $F(4)=2.960516287$, $F(10)=4.668187834$, and $F(100)=9.217221972$. The real target is $F(10000)$, so a brute-force search over all lattice paths is far too large. Mathematical Approach The implementations solve the problem on the left half of the picture first and then mirror that half-path. Write $$w=\frac d2.$$ The task is therefore to understand good paths from $(0,1)$ to the vertical line $x=w$, and then double the best half-cost. Step 1: Rewrite Every Segment Cost For one segment from $(x_0,y_0)$ to $(x_1,y_1)$, let $$\Delta x=x_1-x_0,\qquad \Delta y=y_1-y_0,\qquad L=\sqrt{(\Delta x)^2+(\Delta y)^2}.$$ The travel time is $L/v$. It is convenient to write it with a reciprocal-speed factor $$\psi(y_0,y_1)=\begin{cases} \dfrac1{y_0}, & y_0=y_1,\\[4pt] \dfrac{\ln y_1-\ln y_0}{y_1-y_0}, & y_0\ne y_1....

Detailed mathematical approach

Problem Summary

The ant starts at $A(0,1)$ and must reach $B(d,1)$, where $d\in 2\mathbb{Z}_{>0}$. Every move is a straight segment between lattice points, always going forward in $x$. If a segment starts at height $y_0$ and ends at height $y_1$, then its speed is

$$v=\begin{cases} y_0, & y_0=y_1,\\ \dfrac{y_1-y_0}{\ln y_1-\ln y_0}, & y_0\ne y_1. \end{cases}$$

The objective is to minimize the total travel time; that minimum is denoted by $F(d)$. The published checkpoints are $F(4)=2.960516287$, $F(10)=4.668187834$, and $F(100)=9.217221972$. The real target is $F(10000)$, so a brute-force search over all lattice paths is far too large.

Mathematical Approach

The implementations solve the problem on the left half of the picture first and then mirror that half-path. Write

$$w=\frac d2.$$

The task is therefore to understand good paths from $(0,1)$ to the vertical line $x=w$, and then double the best half-cost.

Step 1: Rewrite Every Segment Cost

For one segment from $(x_0,y_0)$ to $(x_1,y_1)$, let

$$\Delta x=x_1-x_0,\qquad \Delta y=y_1-y_0,\qquad L=\sqrt{(\Delta x)^2+(\Delta y)^2}.$$

The travel time is $L/v$. It is convenient to write it with a reciprocal-speed factor

$$\psi(y_0,y_1)=\begin{cases} \dfrac1{y_0}, & y_0=y_1,\\[4pt] \dfrac{\ln y_1-\ln y_0}{y_1-y_0}, & y_0\ne y_1. \end{cases}$$

Then every segment contributes

$$\tau\bigl((x_0,y_0)\to(x_1,y_1)\bigr)=L\,\psi(y_0,y_1).$$

When $y_1\to y_0$, the second branch tends to $1/y_0$, so the horizontal and non-horizontal cases fit together continuously. For $y_0\ne y_1$, the speed is the logarithmic mean of $y_0$ and $y_1$, and $\psi$ is its reciprocal.

Step 2: Use Symmetry and Work on a Half-Path

The endpoints $A(0,1)$ and $B(d,1)$ are symmetric with respect to the line $x=w$. Reflecting a path in that line preserves Euclidean lengths and also preserves the pair of endpoint heights of every reflected segment, so it preserves travel time as well.

The C++, Python, and Java implementations therefore search for a left half-path from $(0,1)$ to some midpoint $(w,Y)$, then reflect it to obtain a full path from $A$ to $B$. Under that model,

$$F(d)=2\min_{1\le Y\le w} D(w,Y),$$

where $D(x,y)$ denotes the minimum time to reach $(x,y)$ on the left half. The searched half-path is restricted to nondecreasing heights, so the reflected second half automatically supplies the matching descent back to $y=1$.

Step 3: Dynamic Programming Recurrence

Set the base state

$$D(0,1)=0.$$

For every lattice point $(x,y)$ with $0\le x\le w$ and $1\le y\le w$, the half-path cost is obtained from an earlier point $(x_0,y_0)$ with $x_0\le x$ and $y_0\le y$:

$$D(x,y)=\min_{\substack{0\le x_0\le x\\1\le y_0\le y}}\left(D(x_0,y_0)+\sqrt{(x-x_0)^2+(y-y_0)^2}\,\psi(y_0,y)\right).$$

This is a shortest-path computation on an acyclic state graph: $x$ never decreases, and on the searched half-path the height never decreases either. Any state on the boundary $x=w$ is already a complete half-solution and yields the full candidate time $2D(w,y)$.

Step 4: Worked Example for $d=4$

For $d=4$, we have $w=2$. A natural half-path is

$$ (0,1)\to(1,2)\to(2,2). $$

The first segment has length $\sqrt2$ and reciprocal speed

$$\psi(1,2)=\frac{\ln 2-\ln 1}{2-1}=\ln 2.$$

So its time is $\sqrt2\ln 2$. The second segment is horizontal at height $2$, length $1$, hence time $1/2$. The half-cost is therefore

$$D(2,2)=\sqrt2\ln 2+\frac12.$$

Reflecting this half-path gives the full path

$$ (0,1)\to(1,2)\to(2,2)\to(3,2)\to(4,1), $$

and thus

$$F(4)=2\left(\sqrt2\ln 2+\frac12\right)=2.9605162869\ldots,$$

which matches the stated checkpoint $2.960516287$. This is a useful small-scale example of the general recurrence.

Step 5: Why the Full DP Is Still Too Large

If we allowed every state $(x,y)$ to examine every predecessor $(x_0,y_0)$, then there would be $O(w^2)$ states and up to $O(w^2)$ predecessors per state. For $d=10000$, we have $w=5000$, so a literal implementation of the full recurrence would be far too slow.

The implementations therefore keep the same segment-cost formula and the same half-path DP idea, but prune the search to a narrow region where near-optimal states are expected to lie.

Step 6: Restrict the Search to a Reference Arc

The reference curve used by the implementations is the quarter-circle

$$x_{\mathrm{arc}}(y)=w-\sqrt{w^2-y^2},\qquad 1\le y\le w,$$

which is the left branch of the circle centered at $(w,0)$ with radius $w$.

Only source states near that circle are expanded, namely those satisfying

$$\left|\sqrt{(w-x_0)^2+y_0^2}-w\right|\le \varepsilon.$$

For each target height $y$, destination $x$-coordinates are restricted to a narrow horizontal band:

$$\max\bigl(x_0,\ x_{\mathrm{arc}}(y)-W\bigr)\le x\le \min\bigl(w,\ x_{\mathrm{arc}}(y)+1\bigr).$$

In the released implementations, the default values are $W=30$ and $\varepsilon=0.5$. This is not the unrestricted theoretical DP; it is a numerically guided pruning strategy. However, it reproduces all published checkpoints and is exactly the algorithm used by the C++, Python, and Java implementations.

How the Code Works

The C++, Python, and Java implementations allocate a table of best half-path times for all lattice states with $0\le x\le w$ and $1\le y\le w$. Every entry starts at infinity except the initial state $(0,1)$, whose cost is $0$.

They then precompute two geometric ingredients: $\ln y$ for each height $y$, and the reference-arc position $x_{\mathrm{arc}}(y)$. Precomputing logarithms is important because the reciprocal-speed factor is evaluated repeatedly for many pairs of heights.

Next, for each horizontal position $x_0$, the implementation precomputes the heights $y_0$ that lie inside the thin annulus around the reference circle. Those lattice points are treated as the expandable source states.

When one source state is expanded, the implementation sweeps through all target heights $y\ge y_0$. For each such $y$, it evaluates the reciprocal-speed factor once, computes the admissible horizontal band around the reference arc, and relaxes all destinations $(x,y)$ inside that band.

Whenever a relaxation reaches the boundary $x=w$, the best half-answer seen so far is updated. After all source states have been processed, the final result returned is exactly twice that best half-cost, formatted to nine decimal places.

Complexity Analysis

Let $w=d/2$, and let $W$ be the horizontal band width. The unpruned monotone DP has $O(w^2)$ states and $O(w^2)$ predecessor checks per state, so its naive time complexity is $O(w^4)$, with $O(w^2)$ memory.

In the implemented version, the annulus around the reference circle has constant thickness, so it contains only $O(w)$ relevant source states. Each expanded source considers $O(wW)$ destination cells, because it scans all target heights and only a band of width $O(W)$ in the horizontal direction. That gives a practical running cost of roughly $O(w^2W)$, while memory remains $O(w^2)$.

For the actual target $d=10000$, we have $w=5000$ and the default $W=30$, which is why the pruned method is feasible whereas the unrestricted DP is not.

Footnotes and References

Problem page: https://projecteuler.net/problem=460
Logarithmic mean: Wikipedia — Logarithmic mean
Dynamic programming: Wikipedia — Dynamic programming
Shortest path problem: Wikipedia — Shortest path problem
Circle: Wikipedia — Circle

Problem 460 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <vector>

namespace {

struct Options {
    int d = 10'000;
    int window = 30;
    double radial_eps = 0.5;
    bool run_checkpoints = true;
};

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

bool parse_double_after_prefix(const std::string& arg, const std::string& prefix, double& out) {
    if (arg.rfind(prefix, 0U) != 0U) return false;
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) return false;
    try {
        out = std::stod(tail);
    } catch (...) {
        return false;
    }
    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, "--d=", options.d)) continue;
        if (parse_int_after_prefix(arg, "--window=", options.window)) continue;
        if (parse_double_after_prefix(arg, "--radial-eps=", options.radial_eps)) continue;

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.d >= 2 && (options.d % 2 == 0) && options.window >= 1 && options.radial_eps > 0.0;
}

double solve_fast(const int d, const int window, const double radial_eps) {
    const int w = d / 2;
    const int h = w;
    const double inf = std::numeric_limits<double>::infinity();

    std::vector<std::vector<double>> g(static_cast<std::size_t>(h + 1),
                                       std::vector<double>(static_cast<std::size_t>(w + 1), inf));
    g[1][0] = 0.0;

    std::vector<double> log_y(static_cast<std::size_t>(h + 1), 0.0);
    for (int y = 1; y <= h; ++y) {
        log_y[static_cast<std::size_t>(y)] = std::log(static_cast<double>(y));
    }

    std::vector<int> x_ref(static_cast<std::size_t>(h + 1), 0);
    for (int y = 1; y <= h; ++y) {
        const int inner = w * w - y * y;
        x_ref[static_cast<std::size_t>(y)] = w - static_cast<int>(std::sqrt(static_cast<double>(inner)));
    }

    std::vector<std::vector<int>> source_y(static_cast<std::size_t>(w));
    const double outer2 = (static_cast<double>(w) + radial_eps) * (static_cast<double>(w) + radial_eps);
    const double inner_r = std::max(0.0, static_cast<double>(w) - radial_eps);
    const double inner2 = inner_r * inner_r;
    for (int x0 = 0; x0 < w; ++x0) {
        const double dx = static_cast<double>(w - x0);
        const double dx2 = dx * dx;
        const double low = inner2 - dx2;
        const double high = outer2 - dx2;
        if (high < 1.0) continue;

        const int y_lo = std::max(1, static_cast<int>(std::ceil(std::sqrt(std::max(0.0, low)))));
        const int y_hi = std::min(h, static_cast<int>(std::floor(std::sqrt(std::max(0.0, high)))));
        if (y_lo > y_hi) continue;

        auto& ys = source_y[static_cast<std::size_t>(x0)];
        ys.reserve(static_cast<std::size_t>(y_hi - y_lo + 1));
        for (int y0 = y_lo; y0 <= y_hi; ++y0) {
            const double r = std::sqrt(dx2 + static_cast<double>(y0) * static_cast<double>(y0));
            if (std::fabs(r - static_cast<double>(w)) <= radial_eps + 1e-12) {
                ys.push_back(y0);
            }
        }
    }

    double best_half = inf;

    for (int x0 = 0; x0 < w; ++x0) {
        const auto& ys = source_y[static_cast<std::size_t>(x0)];
        for (const int y0 : ys) {
            const double base = g[static_cast<std::size_t>(y0)][static_cast<std::size_t>(x0)];
            if (!std::isfinite(base)) continue;

            for (int y = y0; y <= h; ++y) {
                const int dy = y - y0;
                const double inv_v = (dy == 0)
                                         ? 1.0 / static_cast<double>(y0)
                                         : (log_y[static_cast<std::size_t>(y)] -
                                            log_y[static_cast<std::size_t>(y0)]) /
                                               static_cast<double>(dy);

                const int xr = x_ref[static_cast<std::size_t>(y)];
                const int lx = std::max(x0, xr - window);
                const int rx = std::min(w, xr + 1);

                for (int x = lx; x <= rx; ++x) {
                    const int dx = x - x0;
                    const double len = std::sqrt(static_cast<double>(dx * dx + dy * dy)) * inv_v;
                    double& dst = g[static_cast<std::size_t>(y)][static_cast<std::size_t>(x)];
                    const double cand = base + len;
                    if (cand < dst) {
                        dst = cand;
                        if (x == w && cand < best_half) best_half = cand;
                    }
                }
            }
        }
    }

    return 2.0 * best_half;
}

bool close_to(const double a, const double b, const double eps) {
    return std::fabs(a - b) <= eps;
}

bool run_checkpoints(const int window, const double radial_eps) {
    if (!close_to(solve_fast(4, window, radial_eps), 2.960516287, 5e-10)) {
        std::cerr << "Checkpoint failed: F(4)\n";
        return false;
    }
    if (!close_to(solve_fast(10, window, radial_eps), 4.668187834, 5e-10)) {
        std::cerr << "Checkpoint failed: F(10)\n";
        return false;
    }
    if (!close_to(solve_fast(100, window, radial_eps), 9.217221972, 5e-10)) {
        std::cerr << "Checkpoint failed: F(100)\n";
        return false;
    }
    return true;
}

}  // namespace

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

    if (options.run_checkpoints && !run_checkpoints(options.window, options.radial_eps)) {
        return 2;
    }

    const double ans = solve_fast(options.d, options.window, options.radial_eps);
    std::cout << std::fixed << std::setprecision(9) << ans << '\n';
    return 0;
}

Python

import math

def solve():
    d = 10000
    window = 30
    radial_eps = 0.5
    w = d // 2
    h = w
    INF = float('inf')

    g = [[INF] * (w + 1) for _ in range(h + 1)]
    g[1][0] = 0.0

    log_y = [0.0] * (h + 1)
    for y in range(1, h + 1):
        log_y[y] = math.log(y)

    x_ref = [0] * (h + 1)
    for y in range(1, h + 1):
        inner = w * w - y * y
        x_ref[y] = w - int(math.sqrt(max(0, inner)))

    # Precompute source_y for each x0
    source_y = [[] for _ in range(w)]
    outer2 = (w + radial_eps) ** 2
    inner_r = max(0.0, w - radial_eps)
    inner2 = inner_r ** 2
    for x0 in range(w):
        dx = w - x0
        dx2 = dx * dx
        low = inner2 - dx2
        high = outer2 - dx2
        if high < 1.0: continue
        y_lo = max(1, math.ceil(math.sqrt(max(0.0, low))))
        y_hi = min(h, int(math.sqrt(max(0.0, high))))
        for y0 in range(y_lo, y_hi + 1):
            r = math.sqrt(dx2 + y0 * y0)
            if abs(r - w) <= radial_eps + 1e-12:
                source_y[x0].append(y0)

    best_half = INF
    for x0 in range(w):
        for y0 in source_y[x0]:
            base = g[y0][x0]
            if not math.isfinite(base): continue
            for y in range(y0, h + 1):
                dy = y - y0
                inv_v = (1.0 / y0) if dy == 0 else (log_y[y] - log_y[y0]) / dy
                xr = x_ref[y]
                lx = max(x0, xr - window)
                rx = min(w, xr + 1)
                for x in range(lx, rx + 1):
                    dxi = x - x0
                    length = math.sqrt(dxi * dxi + dy * dy) * inv_v
                    cand = base + length
                    if cand < g[y][x]:
                        g[y][x] = cand
                        if x == w and cand < best_half:
                            best_half = cand

    return f'{2.0 * best_half:.9f}'

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

Java

public class Euler460 {

    static double solveFast(int d, int window, double radialEps) {
        int w = d / 2;
        int h = w;
        double inf = Double.POSITIVE_INFINITY;

        double[][] g = new double[h + 1][w + 1];
        for (int i = 0; i <= h; i++) {
            for (int j = 0; j <= w; j++) {
                g[i][j] = inf;
            }
        }
        g[1][0] = 0.0;

        double[] logY = new double[h + 1];
        for (int y = 1; y <= h; y++) {
            logY[y] = Math.log(y);
        }

        int[] xRef = new int[h + 1];
        for (int y = 1; y <= h; y++) {
            int inner = w * w - y * y;
            xRef[y] = w - (int) Math.sqrt(inner);
        }

        int[][] sourceY = new int[w][];
        double outer2 = (w + radialEps) * (w + radialEps);
        double innerR = Math.max(0.0, w - radialEps);
        double inner2 = innerR * innerR;

        for (int x0 = 0; x0 < w; x0++) {
            double dx = w - x0;
            double dx2 = dx * dx;
            double low = inner2 - dx2;
            double high = outer2 - dx2;
            if (high < 1.0) {
                sourceY[x0] = new int[0];
                continue;
            }

            int yLo = Math.max(1, (int) Math.ceil(Math.sqrt(Math.max(0.0, low))));
            int yHi = Math.min(h, (int) Math.floor(Math.sqrt(Math.max(0.0, high))));
            if (yLo > yHi) {
                sourceY[x0] = new int[0];
                continue;
            }

            int[] temp = new int[yHi - yLo + 1];
            int count = 0;
            for (int y0 = yLo; y0 <= yHi; y0++) {
                double r = Math.sqrt(dx2 + (double) y0 * y0);
                if (Math.abs(r - w) <= radialEps + 1e-12) {
                    temp[count++] = y0;
                }
            }
            int[] ys = new int[count];
            System.arraycopy(temp, 0, ys, 0, count);
            sourceY[x0] = ys;
        }

        double bestHalf = inf;

        for (int x0 = 0; x0 < w; x0++) {
            int[] ys = sourceY[x0];
            for (int y0 : ys) {
                double base = g[y0][x0];
                if (!Double.isFinite(base))
                    continue;

                for (int y = y0; y <= h; y++) {
                    int dy = y - y0;
                    double invV = (dy == 0) ? 1.0 / y0 : (logY[y] - logY[y0]) / dy;

                    int xr = xRef[y];
                    int lx = Math.max(x0, xr - window);
                    int rx = Math.min(w, xr + 1);

                    double[] gy = g[y];
                    for (int x = lx; x <= rx; x++) {
                        int dx = x - x0;
                        double len = Math.sqrt(dx * dx + dy * dy) * invV;
                        double cand = base + len;
                        if (cand < gy[x]) {
                            gy[x] = cand;
                            if (x == w && cand < bestHalf) {
                                bestHalf = cand;
                            }
                        }
                    }
                }
            }
        }

        return 2.0 * bestHalf;
    }

    public static String solve() {
        double ans = solveFast(10000, 30, 0.5);
        return String.format(java.util.Locale.US, "%.9f", ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}