Map index along 2D Hilbert to coordinates (x,y) and back, branchless using bit twiddling.

hilbert.js

/*
Map index along 2D Hilbert curve to coordinates (x,y) and back,
branchless using bit twiddling.

Copyright 2012, 2024 Juha Järvi

Permission to use, copy, modify, and/or distribute this software for any
purpose with or without fee is hereby granted.

THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
*/

/** Interleave bits with 0, so input bits are stored in odd bits and even bits are zeroed.
  * Equivalent to Morton code with one coordinate as the input and the other zero.
  * @param {number} t */

function zip(t) {
	t = (t | (t << 8)) & 0x00ff00ff;
	t = (t | (t << 4)) & 0x0f0f0f0f;
	t = (t | (t << 2)) & 0x33333333;
	t = (t | (t << 1)) & 0x55555555;

	return t;
}

/** De-interleave, dropping even bits and packing the remaining bits into the output.
  * @param {number} t */

function unzip(t) {
	t = t & 0x55555555;
	t = (t | (t >>> 1)) & 0x33333333;
	t = (t | (t >>> 2)) & 0x0f0f0f0f;
	t = (t | (t >>> 4)) & 0x00ff00ff;
	t = (t | (t >>> 8)) & 0x0000ffff;

	return t;
}

/** Compute individual parities of two interleaved numbers.
  * If we consider the input a Morton code with x and y bits interleaved,
  * every bit of x and y becomes the parity of that input number
  * ignoring its less significant bits.
  * @param {number} t */

function parity2(t) {
	t ^= t >>> 2;
	t ^= t >>> 4;
	t ^= t >>> 8;
	t ^= t >>> 16;

	return t;
}

/** Sum (mod 3) all corresponding bit pairs in input values, except that:
  * - Output pairs equal to 0 are represented by 3.
  * - Pairs in first input equal to 3 are treated as 0.
  * - Pairs equal to 0 in input are untouched
  *   (nothing will be added to them regardless of 2nd input).
  * Allowing these quirks makes the code faster, calling code works around them.
  * Could be replaced by calling an equivalent without the quirks like this:
  * return ~sum3(~a, ~b);
  * @param {number} a
  * @param {number} b
  * 16 ops. */

function sum3a(a, b) {
	// return ~sum3(~a, ~b); // Slower alternative.

	// Create mask for bit pairs where most significant bit is unset:
	// 10 01 00 becomes:
	// 00 11 11
	let m = (~b >>> 1) & 0x55555555; // 3
	m |= m << 1; // 2

	// Create mask for least significant bits of pairs where at least one bit is set:
	// 10 01 00 becomes:
	// 01 01 00
	const n = (b | (b >>> 1)) & 0x55555555; // 3

	// Leave bit pairs 00 untouched and consider 11 to mean real 00. Then to each bit pair add:
	// 00 if n=00 and m=11 (original pair is 00)
	// 01 if n=01 and m=11 (original pair is 01)
	// 10 if n=01 and m=00 (original pair is 10)
	return (((a >>> 1) ^ a) & n) ^ ((a & n) << 1) ^ (a & m); // 8
}

/** Sum (mod 3) all corresponding bit pairs in input values, except that:
  * - Pairs equal to 3 in input are ignored.
  * Allowing this quirk makes the code faster, calling code avoids it.
  * Could be replaced by an equivalent without the quirk.
  * @param {number} a
  * @param {number} b
  * 13 ops. */

function sum3b(a, b) {
	// return sum3(a, b); // Slower alternative.

	// Bit pair sums (pairs 11 are always filtered out):
	//    00 00 00 01 01 01 10 10 10
	// +  00 01 10 00 01 10 00 01 10
	// =  00 01 10 01 10 00 10 00 01

	// Mark pairs with matching bits set (these will need flipping because 1+1=2 and 2+2=1 mod 3), fill marks towards left:
	//    00 00 00 01 01 01 10 10 10
	//    00 01 10 00 01 10 00 01 10
	// -> 00 00 00 00 11 00 00 00 10
	let m = a & b; // 1
	m |= (m & 0x55555555) << 1; // 3

	// Initial sum calculation by simple OR:
	//    00 00 00 01 01 01 10 10 10
	//    00 01 10 00 01 10 00 01 10
	// -> 00 01 10 01 01 11 10 11 10
	a = a | b; // 1

	// Also mark pairs where sum=3, needs flipping because 3=0 mod 3.
	// -> 00 00 00 00 11 10 00 10 10
	m |= a & ((a & 0x55555555) << 1) // 4;

	// Fill marks towards right:
	// -> 00 00 00 00 11 11 00 11 11
	m |= (m >>> 1) & 0x55555555; // 3

	// Flip marked bit pairs in sum:
	// -> 00 01 10 01 10 00 10 00 01
	return a ^ m; // 1
}

/** Calculate (x, y) coordinates along a 2D Hilbert curve at position t.
  * @param {number} t Index along curve.
  * @param {{x: number, y: number}} xy Object where x and y fields will be set. */

function toHilbert(t, xy) {
	// Mark occurrences of bit patterns 01 and 10 (ignoring 2 least significant bits), store marks in every second bit starting from least significant.
	let count = ((t ^ (t >>> 1)) & 0x55555555) >>> 2;
	// Mark occurrences of bit pattern 11 (ignoring 2 least significant bits), store marks in every second bit starting from second to least significant.
	count |= ((t & (t >>> 1)) & 0x55555555) >>> 1;

	// Marks form two interleaved bit sequences. Compute their parity in parallel giving cumulative counts mod 2 of marked bits.
	// Then flip every second bit of cumulative counts of bit patterns 01 and 10.
	count = parity2(count) ^ 0x44444444;

	// First take every second bit of position along curve starting from least significant, also masked by the cumulative counts of bit patterns 01 and 10.
	let x = count & t;

	// Then XOR them with the other half of bits of position along curve and the cumulative counts of bit patterns 11.
	x ^= (count ^ t) >>> 1;

	// De-interleave x to get rid of junk bits: Throw away every second bit starting from second to least significant and shift remaining bits right to pack all bits on least significant end.
	x = unzip(x);

	// Calculate Y from X by flipping bits left over after de-interleaving original position along curve.
	xy.x = x;
	xy.y = x ^ unzip(t);
}

/** Calculate index along a 2D Hilbert curve at given (x, y) coordinates.
  * @param {number} x
  * @param {number} y */

function fromHilbert(x, y) {
	// Every second bit of result comes directly from x^y interleaved.
	let h = x ^ y;
	// Flip the curve along x=y diagonal, interchanging x and y.
	x ^= h;

	// Interleave.
	h = zip(h);
	x = zip(x);

	// j contains all adjacent bit pairs of x^y.
	// x^y:	ab cd ef gh
	// h:   0a 0b 0c 0d 0e 0f 0g 0h
	// j:   0a ab bc cd de ef fg gh
	let j = h ^ (h >>> 1);

	// Magic.
	let r = (h & x) ^ (((~h & ~x) & 0x55555555) << 1);

	// Clear bit pairs 11 in j (effectively set bit pairs to remainder modulo 3).
	// 11 10 01 00 becomes:
	// 00 10 01 00
	let m = j & ((j & 0x55555555) << 1);
	j ^= m | (m >>> 1);

	j >>>= 2;

	// Handle curve order 1.
	r ^= (r >>> 2);

	// Handle higher orders.
	r ^= sum3a(r >>> 4, j);
	j = sum3b(j, j >>> 4);
	r ^= sum3a(r >>> 8, j);
	j = sum3b(j, j >>> 8);
	r ^= sum3a(r >>> 16, j);

	// For 64-bit support and curve orders above 16,
	// double the length of all bit masks and add:
	// j = sum3b(j, j >>> 16);
	// r ^= sum3a(r >>> 32, j);

	// Magic.
	return (h ^ (
		(r & (h << 1)) ^
		((r & ~h & 0x55555555) << 1) ^
		(x << 1)
	)) >>> 0;
}

// Supplementary unused functions, might be nice to optimize and replace the current quirky ones!

/** Sum (mod 3) all corresponding bit pairs in input values.
  * @param {number} a
  * @param {number} b
  * 26 ops. */

function sum3(a, b) {
	let m, j2;

	// If both bits of second input are set, copy first input to output at the end.
	m = b & (b >>> 1) & 0x55555555; // 3
	// If both bits of first input are set, set them in output at the end.
	j2 = a & (a >>> 1) & 0x55555555; // 3
	j2 |= (j2 << 1) | ((m | (m << 1)) & a); // 6

	// Bit pair sums (pairs 11 are always filtered out):
	//    00 00 00 01 01 01 10 10 10
	// +  00 01 10 00 01 10 00 01 10
	// =  00 01 10 01 10 00 10 00 01

	// Mark pairs with matching bits set (these will need flipping because 1+1=2 and 2+2=1 mod 3), fill marks towards left:
	//    00 00 00 01 01 01 10 10 10
	//    00 01 10 00 01 10 00 01 10
	// -> 00 00 00 00 11 00 00 00 10
	m = a & b; // 1
	m |= (m & 0x55555555) << 1; // 3

	// Initial sum calculation by simple OR:
	//    00 00 00 01 01 01 10 10 10
	//    00 01 10 00 01 10 00 01 10
	// -> 00 01 10 01 01 11 10 11 10
	a = a | b; // 1

	// Also mark pairs where sum=3, needs flipping because 3=0 mod 3.
	// -> 00 00 00 00 11 10 00 10 10
	m |= a & ((a & 0x55555555) << 1); // 4

	// Fill marks towards right:
	// -> 00 00 00 00 11 11 00 11 11
	m |= (m >>> 1) & 0x55555555; // 3

	// Flip marked bit pairs in sum:
	// -> 00 01 10 01 10 00 10 00 01
	// Then set some bits in output if one of inputs had both bits set.
	return a ^ m ^ j2; // 2
}

/** Sum (mod 3) all corresponding bit pairs in input values.
  * Works identically to sum3.
  * @param {number} a
  * @param {number} b
  * 28 ops. */

function sum3_alternative(a, b) {
	let c, m, n;

	// a 11 11 11 11 10 10 10 10 01 01 01 01 00 00 00 00
	// b 11 10 01 00 11 10 01 00 11 10 01 00 11 10 01 00
	// c 11 11 11 11 10 01 00 10 01 00 10 01 00 10 01 00
	//   ^^       ^^ ^^       ^^ ^^       ^^ ^^       ^^
	// if(b == 00 || b == 11) c = a;
	m = (b ^ (b >>> 1)) & 0x55555555; // 3
	m = m | (m << 1); // 2
	c = a & ~m; // 2

	//      ^^ ^^                               ^^ ^^
	// else {
	//     if(a == 00 || a == 11) c = a | b;
	n = (a ^ (a >>> 1)) & 0x55555555; // 3
	n = ~(n | (n << 1)) & m; // 4
	c |= (a | b) & n; // 3

	//                  ^^             ^^
	//     if(a == b) c |= ~b;
	// }
	n = a ^ b; // 1
	n = (n | (n >>> 1)) & 0x55555555; // 3
	n = ~(n | (n << 1)) & m; // 4

	return c | (~b & n); // 3
}