Tucker Osman

AVX-512 Polyline Encoder

2026-01-25

AVX-512 Polyline Encoder

Maybe the fastest and least portable polyline encoder ever written.

Why?

I don't really know why, but instruction set extensions like AVX-512 fascinate me. It's like a reminder of how insanely complicated the x86 instruction set is. At first, nobody wanted to use AVX-512 because it was slow and power-hungry, but now it's slightly less slow and power-hungry, yet its use is still somewhat rare. You'll find AVX-512 used in video encoding, scientific computing, PS3 emulation, and other niche applications. Maybe AVX-512 and similar extensions are interesting to me because they are so niche. Maybe it's because it's an optimization to squeeze every ounce of performance out of a CPU. I really don't know. But I do know that I wanted to try writing some AVX-512 code.

Inspiration

I was inspired to write this after reading this blog post from Cloudflare, which implements a Euclidean distance calculation using AVX. They didn't use AVX-512 for the reasons mentioned above, but this (and the PS3 emulation thing) made me want to try writing some AVX-512 code. I wanted to do something that was a bit more complex than just adding two vectors together, so I thought about what I could do that would be interesting, and I landed on Google's polyline encoding algorithm.

What is Polyline Encoding?

Polyline encoding is a way to encode a series of latitude and longitude coordinates into a compact string format. It's used by Google Maps and other mapping services to efficiently transmit paths and shapes. The algorithm works by converting the coordinates into a series of integers, then encoding those integers into a string using a series of bitwise operations and base64-like encoding. You can read more about it here. You can kind of get the sense that this algorithm lends itself to SIMD optimizations, since it involves a lot of repetitive operations on an array of numbers.

How It Works

Step 1: Memory Allocation

The first thing we need is a temporary buffer. The polyline algorithm involves several intermediate steps before we get the final ASCII string. This buffer, ibuf, will hold the 32-bit integer results from the initial processing stages. We round up the size to the nearest multiple of 32 to ensure we don't write past the end of our buffer when processing the last, possibly incomplete, chunk of data.

uint32_t *ibuf = (uint32_t*)malloc(sizeof(uint32_t)*ROUND_UP(points*2, 32));

Step 2: The Main Processing Loop (Part 1 - Coordinate Deltas)

We'll loop through all the coordinates, processing 16 double-precision values (which is 8 lat/lon pairs) in each iteration.

for(int i = 0; i < points*2; i+=16) {

Handling Edges with Masking

What if we have, say, 19 points instead of a neat multiple of 16? The last loop iteration would try to read past our input data. To prevent this, we use "masks". A mask is like a stencil; it tells the CPU which data lanes in the register to use and which to ignore. Here, loadm1 and loadm2 are calculated to specify exactly how many valid data points are left to load. If we're not at the end, the mask is 0xFF (binary 11111111), meaning "load all 8 doubles". If we are near the end, the mask will have fewer 1s, preventing out-of-bounds reads. We load two 512-bit registers, in1 and in2, with our coordinate data, safely applying the masks.

  __mmask8 loadm1 = ((i+8-points*2) < 0) ? 0xFF : (0xFF >> (i+8-points*2)) & 0xFF;
  __mmask8 loadm2 = ((i+16-points*2) < 0) ? 0xFF : (0xFF >> (i+16-points*2)) & 0xFF;
  __m512d in1 = _mm512_maskz_loadu_pd(loadm1, a+i);
  __m512d in2 = _mm512_maskz_loadu_pd(loadm2, a+i+8);

Calculating Deltas in Parallel

The polyline algorithm doesn't encode absolute coordinates; it encodes the difference from the previous point. To calculate this, we need to subtract the previous point from the current one (e.g., point[N] - point[N-1]). In our register, the data looks like [lat1, lon1, lat2, lon2, ...]. To subtract, we need to shuffle the data so we can perform the operation [lat2-lat1, lon2-lon1, lat3-lat2, ...].

The _mm512_permutexvar_pd instruction is a powerful shuffle. We use it to create sub1 and sub2 which hold the preceding coordinate for each point in the register. For the very first point in a chunk, we need the last point from the previous chunk. This is handled by the if(i != 0) block, which carefully injects a[i-1] and a[i-2] into the calculation. For the very first chunk, we subtract zero. Finally, we perform the subtraction to get the deltas.

  __m512d sub1 = _mm512_maskz_permutexvar_pd(0xFC, _mm512_set_epi64(5, 4, 3, 2, 1, 0, 0, 0), in1);
  __m512d sub2 = _mm512_maskz_permutexvar_pd(0xFC, _mm512_set_epi64(5, 4, 3, 2, 1, 0, 0, 0), in2);
  if(i != 0)
      sub1 = _mm512_maskz_add_pd(loadm1, sub1, _mm512_set_pd(0, 0, 0, 0, 0, 0, a[i-1], a[i-2]));
  else
      sub1 = _mm512_maskz_add_pd(loadm1, sub1, _mm512_setzero_pd());

  sub2 = _mm512_maskz_add_pd(loadm2, sub2, _mm512_set_pd(0, 0, 0, 0, 0, 0, a[i+7], a[i+6]));

  in1 = _mm512_sub_pd(in1, sub1);
  in2 = _mm512_sub_pd(in2, sub2);

Step 3: Polyline Encoding Logic

Now we follow the core steps of the polyline algorithm, but applied to 16 values at once!

  1. Scale and Round: Multiply by 1e5 and convert to 32-bit integers. This is done with _mm512_mul_pd followed by _mm512_cvtpd_epi32. We then combine our two 256-bit integer results back into a single 512-bit register, mul.
  __m256i mul1 = _mm512_cvtpd_epi32(_mm512_mul_pd(in1, _mm512_set1_pd(1e5)));
  __m256i mul2 = _mm512_cvtpd_epi32(_mm512_mul_pd(in2, _mm512_set1_pd(1e5)));

  __m512i mul = _mm512_inserti32x8(_mm512_setzero_si512(), mul1, 0);
  mul = _mm512_inserti32x8(mul, mul2, 1);
  1. Left Shift: Shift each integer left by one bit. _mm512_slli_epi32.
  __m512i out = _mm512_slli_epi32(mul, 1);
  1. Invert if Negative: If the original number was negative, we need to invert the shifted value (a bitwise NOT). We first find all negative numbers by comparing mul with zero, which creates a mask (_mm512_cmp_epi32_mask). Then, we use that mask to selectively apply a bitwise XOR with -1 (which is equivalent to a NOT) on only the values that were originally negative.
  __mmask16 mask = _mm512_cmp_epi32_mask(mul, _mm512_setzero_si512(), _MM_CMPINT_LT);
  out = _mm512_mask_xor_epi32(out, mask, out, _mm512_set1_epi32(-1));

Finally, we store the processed 32-bit integers into our intermediate buffer ibuf.

  _mm512_storeu_epi32(ibuf + i, out);
}

Step 4: The Second Loop - Converting to 5-bit Chunks

The next stage of the algorithm is to break each integer into a series of 5-bit chunks. Each chunk is then transmitted with a 1-bit continuation flag indicating if there are more chunks for the same number. This is where things get really interesting with bit manipulation instructions. This loop processes 8 integers from ibuf at a time.

uint32_t out_idx = 0;
for(int i = 0; i < points*2; i+=8) {

Spreading Bits with PDEP

This was a large speed boost for me. Initially, I was using bextr to extract five Bits at a time in a loop, but that was slow. Using 8 _pdep_u64 calls in parallel is way faster. Over two times faster, actually. We use it to take the bits from our integer and "deposit" them into a new location defined by a mask. Our mask, 0x1F1F1F1F1F1F, essentially says "take the first 5 bits of the source and place them in the first byte, the next 5 bits in the second byte, and so on". This instantly converts our 32-bit integer into a series of bytes, each containing a 5-bit chunk. The __builtin_bswap64 is used to fix the byte order (endianness) after the deposit. All 8 of these operations are done to prepare the values for a single AVX-512 register x.

  __m512i x = _mm512_set_epi64(__builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+0], 0x1F1F1F1F1F1F)), __builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+1], 0x1F1F1F1F1F1F)),
                               __builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+2], 0x1F1F1F1F1F1F)), __builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+3], 0x1F1F1F1F1F1F)),
                               __builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+4], 0x1F1F1F1F1F1F)), __builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+5], 0x1F1F1F1F1F1F)),
                               __builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+6], 0x1F1F1F1F1F1F)), __builtin_bswap64(_pdep_u64((uint64_t)ibuf[i+7], 0x1F1F1F1F1F1F)));

Finding the End of Each Sequence

We need to know where the last 5-bit chunk is for each original number, because it should not have the continuation bit set. The other chunks should. The next block of code is a very clever, if complex, sequence of SSE instructions (operating on 128-bit values) to create a mask1 that identifies the last valid (non-zero) chunk for each number. This avoids a slow, branching loop. It's a highly specialized bit-twiddling algorithm to essentially "fill in" the gaps in a mask.

I would be lying if I said I fully understood every step of this sequence. It was adapted from a StackOverflow answer that I can't seem to find again. The gist is that we create a mask of all non-zero bytes, then use a series of shifts and minimums to propagate the information about which bytes are non-zero. The final result is a mask where only the last byte of each sequence of non-zero bytes is set.

  __mmask64 mask1 = _mm512_cmpgt_epi8_mask(x, _mm512_setzero_si512());
  const __m128i nibble_mask = _mm_set1_epi8(0x0F);
  __m128i v = _mm_set_epi64x(mask1, mask1);
  __m128i t;

  t = _mm_and_si128(nibble_mask, v);
  v = _mm_and_si128(_mm_srli_epi16(v, 4), nibble_mask);
  t = _mm_shuffle_epi8(_mm_set_epi64((__m64)0x0001000200010003, (__m64)0x0001000200010008), t);
  v = _mm_shuffle_epi8(_mm_set_epi64((__m64)0x0405040604050407, (__m64)0x0405040604050408), v);
  v = _mm_min_epu8(v, t);
  v = _mm_maskz_unpackhi_epi8(0x5555, v, v);

  __m128i mask = _mm_sllv_epi16(_mm_set1_epi16(-1), v);
  __m128i ormaskv = _mm_slli_epi16(mask, 1);

  mask = _mm_xor_epi32(mask, _mm_set1_epi32(-1));
  ormaskv = _mm_xor_epi32(ormaskv, _mm_set1_epi32(-1));
  mask = _mm_packus_epi16(mask, _mm_setzero_si128());
  ormaskv = _mm_packus_epi16(ormaskv, _mm_setzero_si128());
  mask = _mm_xor_epi32(mask, _mm_set1_epi32(-1));
  ormaskv = _mm_xor_epi32(ormaskv, _mm_set1_epi32(-1));

  mask1 = _mm_extract_epi64(mask, 0);
  __mmask64 ormask = _mm_extract_epi64(ormaskv, 0);

Setting the Continuation Bits

Now we apply the continuation bit. We OR every byte with 0x20 (binary 00100000) to set the 6th bit. Then, we use our calculated ormask to undo this operation for the last chunk of each number, effectively leaving it without the continuation bit.

  __m512i y = _mm512_or_epi64(x, _mm512_set1_epi8(0x20));
  x = _mm512_mask_add_epi8(x, ormask, y, _mm512_setzero_si512());

Add 63

The final step before output is to add 63 to every chunk. We use mask1, which we derived earlier, to ensure we only add 63 to the bytes that are part of our numbers, and not to any zero-padding bytes.

  x = _mm512_maskz_add_epi8(mask1, x, _mm512_set1_epi8(63));

Compressing and Storing the Result

Our register x now contains the final byte values, but they are interspersed with zero-padding (e.g., [B1, B2, 0, 0, B3, B4, B5, 0, ...]). We need to compact these bytes together.

First, we create a mask (mask2) of all non-zero bytes. Then, _mm512_maskz_compress_epi8 does exactly what it sounds like: it takes the bytes indicated by mask2 and packs them tightly to the beginning of the register. The number of bytes we wrote is simply the population count (popcnt) of the mask. We store the compressed block of bytes into our final outbuf and increment our output index by the number of bytes written.

  __mmask64 mask2 = _mm512_cmpgt_epi8_mask(x, _mm512_setzero_si512());
  x = _mm512_maskz_compress_epi8(mask2, x);
  _mm512_storeu_epi64(outbuf + out_idx, x);
  out_idx += _mm_popcnt_u64(mask2);
}

And that's it! After the loop finishes, outbuf contains the complete, compact polyline string. We just need to free the temporary buffer we allocated at the start. By using AVX-512, we've transformed a serial, one-by-one process into a massively parallel operation.

free(ibuf);

← Back to Sidequests