Qubyte Codes

SPOILER ALERT: If you're doing the 2017 Advent of Code, you may not want to read onward.

Over the holiday period I found a little time to work on the 2017 advent of code challenges. The later ones get quite involved, and I particularly enjoyed day 20. The second part of the task involves a 3D grid containing accelerating particles. Particles which collide at a grid point (occupy that point at the same time) are removed, and we're asked to determine how many particles remain after all collisions have occurred. This one was fun for me because the solution I came up with involved a little mathematics to avoid brute forcing the answer.

At this point in the challenges, writing a parser for an input file is something I'll skip over. The parser gives back an array of objects, each representing a particle, with $p$, $v$, and $a$ fields to represent position, velocity, and acceleration respectively. These quantities are all represented as length three arrays for the $x$, $y$, and $z$ axes. The input looks like:

[
  {
    "p": [-4897, 3080, 2133],
    "v": [-58, -15, -78],
    "a": [17, -7, 0]
  },
  {
    "p": [395, -997, 4914],
    "v": [-30, 66, -69],
    "a": [1,-2,-8]
  },
  {
    "p": [-334, -754, -567],
    "v": [-31, 15, -34],
    "a": [3,1,4]
  },
  // ...
]

The general algorithm I used was to calculate a collision times for each pair of particles, order them by time ascending, and then remove all pairs for each time step when both colliding particles are still in the set of all particles up to that time (earlier collisions can invalidate later collisions since collisions remove colliding particles).

To calculate potential collisions, I used the following:

const collisions = [];

for (let i = 0; i < input.length - 1; i++) {
    for (let j = i + 1; j < input[i].length; j++) {
        const t = collisionTime(input[i], input[j]);

        if (t !== null) {
            if (!collisions[t]) {
                collisions[t] = [];
            }

            collisions[t].push([i, j]);
        }
    }
}

I've used the index of each particle in the input array as its ID for convenience. The snippet grabs all pairs and uses the function collisionTime, which returns either an integer (time of their first collision from $t=0$), or null if they never collide. The collisions array is indexed by time step, and populated with arrays of collisions. I'll come back to collisionTime later, because that's where interesting things happen.

Next I create a set containing the IDs of all particles.

const particles = new Set(Array.from({ length: input.length }, (v, i) => i));

See my Array.from tip for a partial explanation of this, and the MDN article on Set for the rest.

The rest is looping over collision times to discover which collisions are valid (both particles have not yet collided) and remove the particles which collide from the particles set. The size of the set after all collisions have been checked is the number of remaining particles, which is the goal of this problem.

for (const collisionsAtTime of collisions) {
    if (!collisionsAtTime) {
        continue;
    }

    const toDelete = new Set();

    for (const [a, b] of collisionsAtTime) {
        if (particles.has(a) && particles.has(b)) {
            toDelete.add(a);
            toDelete.add(b);
        }
    }

    for (const index of toDelete) {
        particles.delete(index);
    }
}

console.log(particles.size);

I don't care about the time of each set of collisions; I only care about the order. Since collisions is a sparse array (indexed by time step), I do have to check for undefined elements and continue to skip those.

I've modelled all collisions as pairs, so a collision between three or more particles will look like multiple collisions of pairs at the same time. In order to properly count these, particle IDs to be removed are held in the toDelete set and removed after all collisions at a given time step have been checked.

Let's go back to the function collisionTime. For a given pair of particle objects, it determines when their first collision is beginning at $t=0$. Here's a naïve implementation:

// Calculates the Manhatten distance between two particles.
function manhatten(a, b) {
    return Math.abs(a[0] - b[0]) + Math.abs(a[1] - b[1]) + Math.abs(a[2] - b[2]);
}

function vectorAdd(a, b) {
    return [a[0] + b[0], a[1] + b[1], a[2] + b[2]];
}

// Returns a new particle which represents the evolution of a given particle by
// one time step.
function tick(particle) {
    const a = particle.a.slice();
    const v = vectorAdd(particle.v, a);
    const p = vectorAdd(particle.p, v);

    return { p, v, a };
}

function collisionTime(a, b) {
    let nextD = manhatten(a.p, b.p);

    // Trivial case in which particles are initially in the same place.
    if (nextD === 0) {
        return 0;
    }

    let t = 0;
    let d;
    let v;
    let nextV = manhatten(a.v, b.v);
    let particle0 = a;
    let particle1 = b;

    do {
        d = nextD;
        v = nextV;
        particle0 = tick(particle0);
        particle1 = tick(particle1);
        t++;
        nextD = manhatten(particle0.p, particle1.p);

        if (nextD === 0) {
            return t;
        }

        nextV = manhatten(particle0.v, particle1.v);

        // Continue while the particles are moving together, or failing that,
        // the rate at which they're moving apart is slowing (indicating that
        // they will move toward each other in the future).
    } while (nextD < d || nextV < v);

    return null;
}

This implementation takes a pair of particles, and evolves both while either they are approaching each other, or accelerating toward each other. If the two particles are at zero distance during the evolution, they have collided. It works, but it's slow.

All I really need is the location of each particle at a given time. Since the particles move under a very simple set of rules, I should be able to predict where a particle will be at any time without needing to resort to loops. If I can predict the location of two particles, I can calculate the relative location of one particle with respect to another, and ultimately I can solve an equation for which that location is the same (zero distance).

First I have to predict the location of a particle though. We are told that at each time, the acceleration vector is added to the velocity vector, and then the updated velocity vector is added to the position vector. The position of a particle is thus (in one dimension, where zeros denote initial values):

$$p_t = p_0 + v_0 t + a\sum_{n=0}^t n$$

That sum can be swapped out for an expression (high schoolers will know this but I, with the full weight of two degrees in physics, had forgotten and had to look it up). This is equivalent:

$$p_t = p_0 + v_0 t + a\frac{(t + 1)t}{2}$$

Armed with this equation, I need to calculate the difference in position between two particles:

$$\Delta p_t = \Delta p_0 + \Delta v_0 t + \Delta a\frac{(t + 1)t}{2}$$

Where the delta signifies a difference (the above is the equation of one particle subtracted from the equation for a another particle). The solutions for $t$ I want are when the left hand side (difference in position at time step $t$) is zero. To avoid a division and retain integers as long as possible, I multiply both sides by two and rewrite slightly to collect powers of $t$.

$$0 = \Delta a t^2 + \left(2\Delta v_0 + \Delta a\right) t + 2\Delta p_0$$

This is a quadratic equation! When the acceleration is zero this collapses down to a linear equation with the solution:

$$t = -\frac{\Delta p_0}{\Delta v_0}$$

When acceleration is non-zero, I have to solve a quadratic equation:

$$t = \frac{-b\pm\sqrt{b^2-4ac}}{2a}$$

where

$$\begin{align} a &= \Delta a \\ b &= 2\Delta v_0 + \Delta a \\ c &= 2\Delta p_0 \end{align}$$

I'm interested in solutions which are natural numbers (real positive integers). To get natural number solutions to quadratic equations, I wrote these functions:

function isNaturalNumber(num) {
    return num >>> 0 === num;
}

// Provides only solutions which are natural numbers.
function solveQuadratic(a, b, c) {
  if (a === 0) {
    const solution = -c / b;

    return isNaturalNumber(solution) ? [solution] : [];
  }

  const rootDiscriminant = Math.sqrt(b * b - 4 * a * c);

  return [
    (-b + rootDiscriminant) / (2 * a),
    (-b - rootDiscriminant) / (2 * a)
  ].filter(isNaturalNumber);
}

The function isNaturalNumber uses the right bit shift operator. This is a trick to get an integer out of a number. Positive numbers will be floored, negative numbers will overflow (become large positive integers), and NaN will become zero. Only positive integers will pass the strict equality check. solveQuadratic itself encodes the solution to both the linear case (acceleration of zero) and the quadratic case. JavaScript is on our side here with how it represents numbers.

This solution is in one dimension. For collisions in three dimensions, this equation must be solved for each, and only when all three dimensions have a solution which is the same can I say that the two particles are projected to collide. At last, here is the improved implementation of collisionTime:

function collisionTime(a, b) {
    let commonSolutions;

    // Calculate the solutions for each axis and initialise or whittle down
    // commonSolutions.
    for (let i = 0; i < 3; i++) {
        const dp = a.p[i] - b.p[i];
        const dv = a.v[i] - b.v[i];
        const da = a.a[i] - b.a[i];

        const axisSolutions = solveQuadratic(da, 2 * dv + da, 2 * dp);

        // Most of the time there are no solutions, so this is a good
        // optimisation.
        if (axisSolutions.length === 0) {
          return null;
        }

        if (!commonSolutions) {
            // Initialise the set with solutions for the first axis.
            commonSolutions = new Set(axisSolutions);
        } else {
            // Remove solutions for other axes which this axis does't have.
            for (const solution of commonSolutions) {
                if (!axisSolutions.includes(solution)) {
                    commonSolutions.delete(solution);
                }
            }
        }

        // If the set of common solutions is ever empty, there will be no
        // collision.
        if (commonSolutions.size === 0) {
            return null;
        }
    }

    // I only want the earliest collision.
    return Math.min(...commonSolutions);
}

This function first calculates solutions for the x-axis. Solutions which are not in common with both the y-axis and z-axis are subsequently removed. If the set of common solutions is ever empty, the function returns null. If any axis has no solutions I can take a shortcut and immediately return null. If the set has remaining entries, these represent times when the particles will occupy the same coordinates in all three dimensions (a collision). The lowest (earliest) is plucked out as the solution!

On my machine, the naïve implementation takes ~3.2 seconds with the input I'm given. With the improved implementation, it takes less than 0.1 seconds. The full solution can be found in this gist.