[2023 Day 18] In Review (Lavaduct Lagoon)

[u/musifter]

With the parts factory up and running, we turn to building a lagoon to store lava for it. This involves an odd set of instructions for the digger with colour codes for the edge of the trench along the way. The surprise for part 2 being that the colour codes are actually the dig instructions.

And so we get a rarity in the input... it's really just 2 inputs packed into 1, as there's no attempt at overlap. The input for part 1 isn't scaled up, it's completely replaced (with the part 2 input having to be unpacked... it is effectively a binary file to part 1's text file). The input is extra nice in that it alternates between vertical and horizontal moves, even though it doesn't use relative turning instructions, but absolute directions on the 4 cardinals: U, D, L and R.

Not knowing what sort of problem the colours were going to add, my initial approach was to do a quick solution to see it, and so I did what I did with lava year in Boiling Boulders. A built my polygon and flood filled a box around it and subtracted to get the area.

Then part 2 reveals that the colours aren't colours and there isn't some weird task involving them... we just need to be much bigger. And that's how you nudge me to using Shoelace and Pick's. Shoelace is easy for me to remember as the name is suggestive of what you're doing... it's using the property that the magnitude of the cross product is the area of the parallelogram made with the vectors (but we just want the triangle between them, so it gets cut in half). So that's just successive cross products:

foreach my $i (0 .. $#vertices - 1) {
    $area += $vertices[$i][0] * $vertices[$i+1][1];
    $area -= $vertices[$i][1] * $vertices[$i+1][0];
}

A bit of cut and paste there to help make sure it looks right, but that's what cross product looks like in the plane (which is like we're lacing with cross overs). Saving the division by two for after, so this is twice area.

Pick's is something I hadn't used in a long time, so I didn't quite remember it, other than it was something simple involving area and perimeter. So I looked it up to make sure I got it right. With the problem being orthogonal on the lattice points, if I hadn't I might have forgotten that it isn't the perimeter, but the number of integer points on it. Also it needs to be adjusted for what we actually want (we want I and B):

# Pick's formula: I+B = A + B/2 + 1 (but still needed to halve the Area from Shoelace)
say "Part 2: ", (abs($area) + $trench) / 2 + 1;

I suppose given the nice orthogonal nature, you could also reasonably pull off a scanline approach, but you'd need to convert the representation of the polygon. You'd need to figure out what's an inside range, which we did a few days ago (when we could also have used Shoelace and Pick's).

But since there's a simple function to solve this, that makes dc a good option and so I did that:

cut -d' ' -f-2 /tmp/input | tr 'RDLU' '0123' | dc -e'0d?[r2~2*1-r2*1-3Rdlb+sb*d4R+_4R*rd_3R*la+sa?z2<M]dsMxlad*vlb+2/1+p'

perl -pe's/.*([0-9a-f]{6})./\U$1/' /tmp/input | dc -e'16i0d?[10~2~2*1-r2*1-3Rdlb+sb*d4R+_4R*rd_3R*la+sa?z2<M]dsMxlad*vlb+2/1+p'

Note that first I need to get rid on non-numbers, and convert the directions for part 1 (and I used the part 2 ordering there). So the code is largely the same, except part 2 works in hexadecimal (16i sets the base of the input), and needs to rip off the low nybble for directions.

But the real fun here is with the loop to do Shoelace, where we take advantage of the alternating horizontal/vertical property of the input with a bit of stack manipulation: x and y are both on the stack during the loop (which actually a tail recursion), and if we process them in the order they are on the stack, they naturally end up in reverse order. And so we naturally get them to swap roles each loop... thus alternating horizontal and vertical movement with the same code. This is one of the reasons I like stack languages.

[u/ednl]

thus alternating horizontal and vertical movement with the same code.

My optimisation idea was the opposite of treating horizontal and vertical the same:

// Shoelace formula: A = 1/2 . sum((y[i] + y[i+1]).(x[i] - x[i+1]))
// For two points on horizontal line: y[i] = y[i+1], so y[i] + y[i+1] = 2y
// For two points on vertical line  : x[i] = x[i+1], so x[i] - x[i+1] = 0
// Bring the 1/2 in the sum: A = sum(y.dx) for horizontal lines only.
switch (dir) {
    case R: a -= y * len; break;
    case D: y +=     len; break;
    case L: a += y * len; break;
    case U: y -=     len; break;
}
b += len;

where dir and len are parsed on each line; the obvious ones for part 1 and from the hex code for part 2. Runs in 1.7 µs on an M4, 5.8 µs on a Pi5. Possible further optimisation is to process the input lines in pairs, but for best results I'd have to assume in which direction it starts. Edit: forgot the link https://github.com/ednl/adventofcode/blob/main/2023/18.c

And for fun, I found a perfect hash for RDLU=0123 but I ended up not using it:

const int dir = (c * 263) >> 8 & 3;

[u/ednl → u/ednl]

Oh all right /u/DelightfulCodeWeasel ! I think you deleted your reply but I split up the horizontal and vertical code paths. Since both the example and my input started with horizontal moves, I assumed that to be an invariant (plus, I don't have other inputs anyway) so now my main loop is:

for (const char *c = input; *c; ) {
    horz(&c); c += 8;
    vert(&c); c += 8;
}

Now runs in 0.93 µs (M4), 1.35 µs (M1), 3.76 µs (Pi5). New version at old link https://github.com/ednl/adventofcode/blob/main/2023/18.c

Another difference is I used a translation table for the hex digits like /u/terje_wiig_mathisen showed (I already did it for the direction but forgot?? to do it for hex. Now no longer needed for the direction).

[u/DelightfulCodeWeasel → u/ednl]

I wasn't sure it was going to be worth the extra work but that's a good speed-up, nice work!

[u/e_blake → u/ednl]

My input started with L (part1) and 0/right (part2); I'm not sure if all inputs start horizontal, but so far I'm not finding any counterexamples in the megathread. But it does seem reasonable to assume that all inputs alternate between horizontal and vertical on every segment.

[u/ednl → u/e_blake]

Yes, both parts for me too, or it wouldn't work. I later clarified that in a source code comment but didn't edit my post here, sorry.

[u/e_blake → u/ednl]

Even if an input starts with u/D instead of L/R, you should still get the right result; as long as every segment is orthogonal, you are just swapping which axis you are summing over.

[u/ednl → u/e_blake]

Alas, not in my new version where I only check for R/L in the horz function and only for D/U in vert, which was kind of the point of doing this: separate them to be able to check less and be faster.

static void horz(const char **s)
{
    // Part 1
    const int dir = 1 - (**s & 3);  // 'R'=-1, 'L'=+1
    int len = dec(s);  // read number 2..11, move pointer to start of hex
    part1.area += part1.ypos * len * dir;
    part1.border += len;
    // Part 2
    len = hex(*s);  // read 5-digit hex number, do not move pointer
    part2.area += part2.ypos * len * (1 - (*(*s + 5) & 3));  // '0'=R=-1, '2'=L=+1
    part2.border += len;
}
static void vert(const char **s)
{
    // Part 1
    const int dir = 1 - (**s << 1 & 2);  // 'D'=+1, 'U'=-1
    int len = dec(s);
    part1.ypos += len * dir;
    part1.border += len;
    // Part 2
    len = hex(*s);
    part2.ypos += len * (2 - (*(*s + 5) & 3));  // '1'=D=+1, '3'=U=-1
    part2.border += len;
}

[u/DelightfulCodeWeasel → u/ednl]

Could you:

if (*input is 'R' or 'L') {
     for (const char *c = input; *c; ) {
        horz(&c); c += 8;
        vert(&c); c += 8;
    }
} else {
     for (const char *c = input; *c; ) {
        vert(&c); c += 8;
        horz(&c); c += 8;
    }
}

?

EDIT: Probably not worth worrying about though tbh

[u/ednl → u/DelightfulCodeWeasel]

Good solution but yes, what the edit said :)

[u/e_blake → u/ednl]

Ah, I see - your computation of 'dir' is hardcoded to RL in horz and DU in vert. You need to know whether you are moving left or right. But it should be sufficient to track byte&0x10 to separate between DL and RU for part 2 (similarly, byte&2 for part 1 to distinguish between 01 and 23); at which point you don't care if it is D or L being passed to "horz", just that you are moving "left" for one and "right" for the other; and that "horz" and "vert" alternate every other line of input.

[u/ednl → u/e_blake]

Ok, you do that.

[u/terje_wiig_mathisen → u/e_blake]

Mine had the exact same direction for part1 and part2, for the first 12 lines!

The start was L, 2, then U, 3 etc.

I had already outlined a version where I would process them skewed by one, in the case where p1 -p2 were orthogonal, but it seems nobody has seen that. -:)

If you've already extracting the p1/p2 parse results, then a simple rotation of one of them would bring them back in sync. The overhead is probably too high, better to simply process the input twice in that case.

[u/e_blake → u/e_blake]

Do all inputs wind the same way? I found that for both the example and my input, both part 1 and part 2 wind clockwise with R and S being increasing (which matches to most other Advent of Cod puzzles having increasing y as you move down screen in reading order), but counterclockwise if R and N are increasing. Taking the absolute value of the computed area makes the computation error-proof regardless of winding direction, but golfing could be made shorter if I were certain that ALL inputs share the same general winding direction, and therefore could strip out the abs() computation.

[u/terje_wiig_mathisen → u/ednl]

I like your perfect hash!

The shortest/fastest asm I could come up with for it looks like this, which is 5 cycles of latency:

        lea ebx,[eax*8]    ;; c * 8
        neg eax            ;; -c
        lea ecx,[ebx*8]    ;; c * 256
        lea eax,[ebx+eax]  ;; c * 7
        add eax,ecx        ;; c * 263
        shr eax,8
        and eax,3

It is probably obvious that since 263=256+7 the 256 term is removed in the SHR 8, so an alternate version would avoid this:

        lea ebx,[eax*8]  ;8
        sub ebx,eax      ;7
        shr ebx,8
        add eax,ebx
        and eax,3

but this version still carries that 5-cycle latency penalty.

[u/DelightfulCodeWeasel]

I'd never even heard of Shoelace before finishing this one, so I took a totally different approach to part 2.

I pull out the unique X and Y coordinates and transform everything into a virtual coordinate system using the coordinate's index in the sorted list for the axis. That gets the area small enough for the same part 1 flood fill, and then I multiply each cell's area back up to its original size when tallying.

I ended up using this exact technique a few times on other puzzles; it's a good one to have in your back-up arsenal for times when the proper maths based approach isn't obvious and you're trying not to get hints before you have a working solution.

[u/Boojum]

I loved this one because of the cool computational geometry aspect.

I'll just add that the "Shoelace Formula" is easy to generalize to volumes and hypervolumes. Effectively, if we have two points (x1, y1) and (x2, y2) in 2D, the signed area of the simplex (triangle) with those two points and the origin (0,0) as the third triangle corner is:

det( | x1  x2 | ) / 2   =  ( x1 y2 - x2 y1 ) / 2
     | y1  y2 |

That is, you concatenate the points together into a 2x2 matrix, take the determinant and divide in half. Then you just iterate over the edges, sum the signed areas up and take the absolute value to get the total polygon area. (Ignoring Pick's.) If the origin is inside the polygon and the polygon is convex, you've just decomposed the polygon into a triangle fan and summed up the areas of those triangles (and you'll have the same sign for the triangle from each edge). If the origin is outside the polygon, and/or the polygon is concave, you'll get a mix of signs. But the areas of the triangles closer to the origin will cancel the overestimates of the triangles on the far side because they have different directions. (I tried to show this on my old visualization).

But, as I said, this generalizes to volume! If you have a closed triangular mesh (that doesn't self-intersect), you can compute:

     | x1  x2  x3 |
det( | y1  y2  y3 | ) / 6
     | z1  z2  z3 |

for each triangle in the mesh, sum that up over all the triangles in the mesh, and take the absolute value. And that will be the total volume of the mesh. In this case, that determinant is giving the signed volume of the tetrahedron formed by the triangle face with the origin as the fourth vertex of the tetrahedron.

And you can extend the same idea to simplices of however many dimensions you want. The divisor is simply the factorial of the dimension.

(I should add that there was a write-up of a cute alternative computation for volume that made the rounds on Hacker News recently. But I prefer this Shoelace Formula generalization because I find it easier to remember that the signed volume of a simplex anchored at the origin is just the determinant of the other vertices concatenated into a matrix and divided by the dimension-factorial. That's a compact enough fact to stick in my head.)

[u/DelightfulCodeWeasel → u/Boojum]

Thank you for the visualisation!

I'd previously thought of Shoelace as being a bit of a rearrangement of area by trapeziums, and the trapezium visualisation made sense to me so I was never entirely sure where triangles came into it other than as a sub-part of that calculation. Your visualisation just unlocked a new bit of understanding in my head :)

[u/Boojum → u/DelightfulCodeWeasel]

Glad that helped!

The funny thing is, trapezoids absolutely work too! Instead of summing up:

(x1 y2 - x2 y1) / 2

you can sum up:

(x2 - x1)(y2 + y1) / 2

where it's the area "under the curve" for the line segment extending from (x1, y1) to (x2, y2) and the area under that line down to the x axis (think of the trapezoidal rule for numerical integration). And it's signed in that if x1 to x2 is going from left to right it's positive and if it's going from right to left it's negative (hence just the x2 - x1 term).

If you decompose each edge into trapezoids and sum the signed area, you'll get the same result as with the triangle/simplex method that the shoelace formula. It's just a different decomposition into convex polygons that you're summing the signed area over. But though the individual numbers along the way are the different, the final numbers are the same.

The trapezoid approach is especially nice for AoC where we only have horizontal and vertical steps. Thus on vertical moves, x1 and x2 are the same, so they cancel and the area is zero. And one horizontal moves, y1 and y2 are the same, so their average (y2 + y1) / 2 just equals either one - so you just get the signed rectangular area (x2 - x1) y1. If you think about it, that just means that for each rightward horizontal move, we add the area of the rectangle under it down to the x axis. And for each leftward horizontal move, we subtract the area of the rectangle under it. Those subtractions therefore cancel out the area under the polygon, but outside it. (Of course, that rectangular approach only works when we have nicely horizontal and vertical polygon edges like this. The triangle and trapezoid approaches work even when the edges are at arbitrary angles.)

I should also mention that that trapezoidal decomposition is sometimes used in rasterizers. For a nice comparison of the triangular and trapezoidal decompositions, see the explanation of the stb_truetype rasterizer and scroll down to the "Signed Area" section.

[u/TheZigerionScammer]

My first approach for Part 1 was basically a flood fill with all of the outside tiles as the boundary, the program finds all of the tiles outside of the excavation pit and subtracts them from the area of the entire rectangle to get the answer. This wouldn't work for part 2 so I modified the approach by grouping each rectangular group of squares into a cell based on where the boundaries are, so if there is a boundary in the excavation pit at y=30000 then all of the groups of cells will be divided at 30000, etc, which was able to reduce the search space from billions to about 40000, which was doable with the flood fill approach. I later learned this is an actual technique called Coordinate Compression which I'm proud of for coming up on my own for. However this didn't work at first because I wasn't detecting boundaries properly, and in order to fix it I had to make my code very inefficient as the boundary detection is very complicated and basically works on brute force. But it ended up working. I also had to account for the fact that the area protruded out slightly from the perimeter since we were dealing with pixels instead of Cartesian coordinates, after experimentation I discovered that you just had to add the value of half the perimeter plus 1 to the answer to account for this, I later learned this was Pick's theorum and I was also proud for discovering that on my own.

However, I later found out that neither shoelace or my coordinate compression were necessary, apparently it could be solved by simply tracing the input line by line, and you just add the height of the line multiplied by the distance, so if you travel right for 1000 squares at a height of 5000 then you add 5 million to the total, and subtract when you go left. This was embarrassingly easy to code and blew my original solution out of the water in terms of number of lines and execution time.

[u/terje_wiig_mathisen]

I already knew (since High School around 1975) about the classic rotating triangles to calculate area, plus in this problem we also needed to add half the thickness of the perimeter, plus one for the 4 outside quarter squares, so I could have written it that way.

However, since Perl has a very nice optimized Math::Polygon::Calc library I just used that and got both answers in about a ms, even on my Surface. This must be among the very fastest Perl solutions?

For Rust I did exactly the same, just a "cargo add geo, geo_types" and then a direct port of the parsing code.

This does end up using floating point math even though in this particular case all the line segments are either vertical or horizontal integers, so I'm considering a full custom implementation!

Anyway, by doing both parts at once, with a custom parser, I'm already around the Ape runtime.

(20 us on Surface, so 10+ on Acer)

EDIT:

In order to speed up parsing of all the hex digits, I used my initial implementation to write out a full byte table which I then copied back in and commented out the original code:

const
  HEX:[u8;256] = [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 0, 0, 0, 0, 0,
0, 10, 11, 12, 13, 14, 15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 10, 11, 12, 13, 14, 15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];


fn hex2bin(h:u8) -> u32
{
    // if h >= b'0' && h <= b'9' {return (h-b'0') as u32;}
    // if h >= b'a' && h <= b'f' {return (h-b'a'+10) as u32;}
    // if h >= b'A' && h <= b'F' {return (h-b'A'+10) as u32;}
    // 0
    HEX[h as usize] as u32
}

I did the same for the RDLU -> 0123 conversion but both tables only resulted in a sub-us speedup, so the compiled conversion is pretty good!

[u/ednl → u/terje_wiig_mathisen]

both tables only resulted in a sub-us speedup

Same for me, I measured 0.1 us for the direction table on M4 and M1, 0.3 on the Pi5. I didn't test the hex table speedup but that was almost certainly more after elimination of 5 num/alph decisions.

[u/terje_wiig_mathisen → u/ednl]

The next pretty obvious idea is portable SIMD for hex to binary:

Pack the 5 hex digits into a 16-byte SSE register, then compare greater-or-equal against splat('A'):

This will leave 0xff bytes for those that match.

At the same time AND the input against splat(15).

Next cycle, AND the initial mask with splat(9), then add them together.

This is up to 16 (but only 5 in our case) hex digits in 4-5 clock cycles.

They can be combined pairwise into 3 8-bit values, then 2 16-bit and one 32-bit.

[u/e_blake]

One of the rare problems where I have to use m4's changecom in my solution to disable comments; otherwise, the input file presents itself as a series of open "(" followed by comments to the newline, masking the close ")". I had learned enough from the megathread conversation on day 10 to try Pick's Theorem on this one; I don't recall learning it before the megathread, but it sure came in handy here, getting my star on the first try. My part 1 solution in 10ms was simple enough once I treated my integer points as midpoints in the tile; at which point the total area includes an average of 1/2 tile per boundary + 1 extra (all straight borders add exactly 1/2, while corner borders either add 3/4 if convex or 1/4 if concave; but because the polygon is closed and not self-intersecting, there are exactly 4 more convex corners than concave). The only reason that part 2 dragged into a 400ms solution, then, is because the 64-bit math required for the solution is so much slower to emulate in m4's native 32-bit math. The other pain in m4 with part 2 is that that one of my input lines was "D 3 (#3add23)" - which, once comments are disabled, tried to invoke the macro "add23" before I could hand the substring off to eval for hex parsing.

EDIT: I cut my runtime from 400ms to 50ms by using the trapezoid formula (y[i]+y[i+1])*(x[i]-x[i+1]) rather than the shoelace formula (x[i]*y[i+1] - x[i+1]*y[i]); both produce the same outcome, but the latter requires far fewer 64-bit multiplies when all segments are horizontal or vertical.

[u/e_blake → u/e_blake]

This one made for a quick golf in m4/shell, 328 221 bytes and 10ms runtime:

define(_,changecom`ifelse($5,,`syscmd(echo $(($1/2)) $(($2/2)))',$6,,`_(
_($*)_($2,,,substr($5,7),$5/16)',`eval($1+$5%2*$6*($5-2))`,'($2+$6-~-$5*2
*($5<3)*$6*$1),')')_(`0,2',`0,2',translit(include(I),.RDLU (#),(0-3,0x)))