Notebook

Calculating ballistic trajectories¶

https://adventofcode.com/2021/day/17

We are tasked with solving an optimisation problem; given a target area, what is the highest a ballistic missile will reach when fired from the submarine.

Given the target area with bounds $x_{min}..x_{max}, y_{min}..y_{max}$ and assuming that the trench is always below the submarine (negative values for $y$), you can already figure out the minimum and maximum velocities:

The minimum $x$ velocity, $vx_{min}$, must reach the nearest side of the target $x_{min}$, and because it decreases monotically with every step essentially is a triangle number. That means $x_{min} <= \frac {vx_{min} (vx_{min} - 1)} {2}$, or $vx_{min} = \lceil \frac 1 2 (-1 + \sqrt {8x_{min} + 1}) \rceil$
$vx_{max}$ should simply not overshoot in a single step, so $vx_{max} = x_{max}$
$vy_{min}$ must hit the bottom of the target area in one step, so $vy_{min} = y_{max}$
The trajectory upwards is exactly mirrored (symetrical) because y changes by -1 each step. If you shoot your probe upwards with initial velocity $y0$, by the time the probe reaches $y=0$ it'll have the same velocity in the opposite direction minus one, $-y0 - 1$. It can't overshoot the target area, to not overshoot the area you can't go beyond the lower bound so the velocity can't be greater than the (absolute) lower bound minus 1: $vy_{max} = |y_{min}| - 1$

The maximum height $h_{max}$ for any given $vy$ is another triangle number, so $h_{max} = \frac {vy (vy + 1)} 2$; and you can reach that maximum height by finding the maximum upwards velocity $vy_{max}$, a value we can calculate without having to know $vx$. So all we need to know is $vy_{max}$, which is simply $|y_{min}| - 1$.

In [1]:

from __future__ import annotations

import re
from dataclasses import dataclass

INPUTLINE = re.compile(
    r"target area: x=(?P<xmin>-?\d+)\.\.(?P<xmax>-?\d+), "
    r"y=(?P<ymin>-?\d+)\.\.(?P<ymax>-?\d+)"
).fullmatch


@dataclass
class TargetArea:
    xmin: int
    xmax: int
    ymin: int
    ymax: int

    @classmethod
    def from_line(cls, line: str) -> TargetArea:
        args = {k: int(v) for k, v in INPUTLINE(line).groupdict().items()}
        return cls(**args)

    @property
    def max_height(self) -> int:
        vy_max = abs(self.ymin) - 1
        return (vy_max * (vy_max + 1)) // 2


assert TargetArea.from_line("target area: x=20..30, y=-10..-5").max_height == 45

In [2]:

import aocd

target = aocd.get_data(day=17, year=2021)
print("Part 1:", TargetArea.from_line(target).max_height)

Part 1: 7875

Part 2, find all possible velocities¶

Now we need to find all possible velocities. We already know the bounds of $vx$ and $vy$, which limits how many integer velocities we need to consider. We can't simply give all combinations of velocities, however; some will overshoot the area in one or the other direction at the exact time $t$ that the downward or sideways movement would be within the target area.

E.g. from the puzzle example, the bounds are $vx=7..30, vy=-5..9$, but shooting the probe with initial velocity (17,-4) doesn't work, even though those two values fall within the velocity bounds. That's because $vx = 17$ will always miss the target bounds. Similarily, there will be values for $vy$ that will always miss the target, and for some values of $vx, vy$ there is no time $t$ where both $x$ and $y$ are within the target area.

Instead, we calculate what bounds of $t$ can be set for a given $vy$ to still hit the area, then verify that $vx$ will hit the bounds in the same range of $t$. That's easy to verify for the limited ranges involved.

The simplest case is if $vy$ is negative; the position $y$ for a given time $t$ is:

$ y = \frac {t (2vy - (t - 1))} {2} $

which can be rearranged as a quadradic equation:

$ -t^2 + (2vy + 1)t - 2y = 0 $

Given that we know the bounds of $y$ we can find bounds for $t$ for any given $vy$, by plugging in those bounds into the positive root of the quadratic equation:

$ t = \frac {1 + 2 vy + \sqrt {1 + 4 vy + 4 vy^2 - 8 y}} {2} $

The bounds for $t$ need to be rounded; to reach $y_{max}$ or beyond, we need to round $t$ up, to not overshoot $y_{min}$ you need to round $t$ down:

$ t_{min} = \lceil \frac {1 + 2 vy + \sqrt {1 + 4 vy + 4 vy^2 - 8 y_{max}}} {2} \rceil $ $ t_{max} = \lfloor \frac {1 + 2 vy + \sqrt {1 + 4 vy + 4 vy^2 - 8 y_{min}}} {2} \rfloor $

If $t_{min} > t_{max}$ then the probe can't reach the target area, so we can skip this $vy$ value.

For $y = 0$, just add 1 to the $t_{min}..t_{max}$ range calculated for $vy = -1$, because that's what the velocity will be at $t = 1$. For positive values of $vy$, you reach $y = 0$ after $t = 2(vy + 1)$, after which it's the same solution as for $vy = -vy - 1$ because that's what the velocity will be at that point. These two relationships mean we can start our search at $vy = y_{min}$ and loop up to $vy = -1$ and take along the $vy >= 0$ solutions from there, where the zero or larger value $pvy = -vy - 1$ and the delta for $t$, $dt = 2(-vy) - 1$.

Once we have determined bounds for $t$, we can calculate all possible $vx$ values that will still fall within the $x_{min}..x_{max}$ range at those times. For a given $vx$, if $t >= vx$ the value comes to rest inside the target area at $t = vx$ at $x = \frac {vx (vx + 1)} {2}$, or, for values where $t < vx$, at a distance $x = \frac {vx (vx - 1)} {2} - \frac {(vx - t) (vx - 1 + 1)} {2}$

I found it easiest to just loop over the product the $t_{min}..t_{max}$ and $vx_{min}..vx_{max}$ ranges, and if the distance falls within the target area increment the counter, and if you overshoot the area, break out of the loop early. I also needed to account for $vx$ values that fall inside the target area for multiple values of $t$.

I may revisit this part later as I'm sure there is a closed-form solution to calculate the bounds of $vx$ for a given $t_{min}..t_{max}$ without looping. You can then just take the length of that range and add that to a counter, rather than generate all possible velocity pairs as tuples.

In [3]:

import math
from typing import Iterable


def _pos_quad(y: int, vy: int) -> float:
    return (1 + (2 * vy) + math.sqrt(1 + (4 * vy) + (4 * vy * vy) - (8 * y))) / 2


def _t_bounds(vy: int, y_min: int, y_max: int) -> tuple[int, int]:
    return math.ceil(_pos_quad(y_max, vy)), int(_pos_quad(y_min, vy))


class ExhaustiveTargetArea(TargetArea):
    def __iter__(self) -> Iterable[tuple[int, int]]:
        xmin, xmax = self.xmin, self.xmax
        vxmin, vxmax = math.ceil((-1 + math.sqrt(8 * xmin + 1)) / 2), xmax
        vymin, vymax = self.ymin, abs(self.ymin) - 1

        def _find_vx(tmin: int, tmax: int) -> Iterable[int]:
            seen = set()
            for t in range(tmin, tmax + 1):
                for vx in range(vxmin, vxmax + 1):
                    if vx in seen:
                        continue
                    x = vx * (vx + 1) // 2
                    if t < vx:
                        x -= (vx - t) * (vx - t + 1) // 2
                    if x > xmax:
                        # this vx overshoots at t, all further will too
                        break
                    elif x >= xmin:
                        seen.add(vx)
                        yield vx

        for vy in range(vymin, 0):
            tmin, tmax = _t_bounds(vy, self.ymin, self.ymax)
            if tmin > tmax:
                # at no point will the probe hit the target
                continue
            yield from ((vx, vy) for vx in _find_vx(tmin, tmax))
            # the positive vy equal to -vy - 1 reaches the same negative
            # velocity once it returns to y=0 at t + 2(-vy) - 1, so for all
            # negative vy where -vy - 1 <= vymax.
            pvy, dt = -vy - 1, 2 * -vy - 1
            if pvy > vymax:
                continue
            yield from ((vx, pvy) for vx in _find_vx(tmin + dt, tmax + dt))

    def __len__(self) -> int:
        return sum(1 for _ in self)


assert len(ExhaustiveTargetArea.from_line("target area: x=20..30, y=-10..-5")) == 112

In [4]:

print("Part 2:", len(ExhaustiveTargetArea.from_line(target)))

Part 2: 2321