diff --git a/spatial/curve/doc.go b/spatial/curve/doc.go
new file mode 100644
index 0000000..847c4cf
--- /dev/null
+++ b/spatial/curve/doc.go
@@ -0,0 +1,38 @@
+// Copyright ©2024 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package curve defines space filling curves. A space filling curve is a curve
+// whose range contains the entirety of a finite k-dimensional space. Space
+// filling curves can be used to map between a linear space and a 2D, 3D, or 4D
+// space.
+//
+// # Hilbert Curves
+//
+// The Hilbert curve is a continuous, space filling curve first described by
+// David Hilbert. The implementation of Hilbert2D is based on example code from
+// the [Hilbert curve (Wikipedia)]. The implementation of Hilbert3D and
+// Hilbert4D are extrapolated from Hilbert2D.
+//
+// Technically, a Hilbert curve is a continuous, fractal, space-filling curve of
+// infinite length, constructed as the limit of a series of piecewise linear
+// curves. We refer to the kth piecewise linear curve as the k-order Hilbert
+// curve. The first-order 2D Hilbert curve looks like ⊐. The k-order
+// n-dimensional curve can be constructed from 2ⁿ copies of the (k-1) order
+// curve. These (finite) Hilbert curves are iterative/recursive, and the order
+// refers to the iteration number.
+//
+// The runtime of Space and Curve scales as O(n∙k) where n is the dimension and
+// k is the order. The length of the curve is 2^(n∙k).
+//
+// # Limitations
+//
+// An n-dimensional, k-order Hilbert curve will not be fully usable if 2ⁿᵏ
+// overflows int (which is dependent on architecture). Len will overflow if n∙k
+// ≥ bits.UintSize-1. Curve will overflow if n∙k > bits.UintSize-1 for some
+// values of v. Space will not overflow, but it cannot be called with values
+// that do not fit in a signed integer, thus only a subset of the curve can be
+// utilized.
+//
+// [Hilbert curve (Wikipedia)]: https://en.wikipedia.org/w/index.php?title=Hilbert_curve&oldid=1011599190
+package curve
diff --git a/spatial/curve/hilbert.go b/spatial/curve/hilbert.go
new file mode 100644
index 0000000..9edecfe
--- /dev/null
+++ b/spatial/curve/hilbert.go
@@ -0,0 +1,334 @@
+// Copyright ©2024 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package curve
+
+import (
+	"errors"
+	"fmt"
+)
+
+// ErrUnderflow is returned by Hilbert curve constructors when the power is less
+// than 1.
+var ErrUnderflow = errors.New("order is less than 1")
+
+// ErrOverflow is returned (wrapped) by Hilbert curve constructors when the
+// power would cause Len and Pos to overflow.
+var ErrOverflow = errors.New("overflow int")
+
+// The size of an int. Taken from src/math/const.go.
+const intSize = 32 << (^uint(0) >> 63) // 32 or 64
+
+// Hilbert2D is a 2-dimensional Hilbert curve.
+type Hilbert2D struct{ order int }
+
+// NewHilbert2D constructs a [Hilbert2D] of the given order. NewHilbert2D
+// returns [ErrOverflow] (wrapped) if the order would cause Len and Pos to
+// overflow.
+func NewHilbert2D(order int) (Hilbert2D, error) {
+	v := Hilbert2D{order: order}
+
+	// The order must be greater than zero.
+	if order < 1 {
+		return v, ErrUnderflow
+	}
+
+	// The product of the order and number of dimensions must not exceed or
+	// equal the size of an int.
+	if order*2 >= intSize {
+		return v, fmt.Errorf("a 2-dimensional, %d-order Hilbert curve will %w", order, ErrOverflow)
+	}
+
+	return v, nil
+}
+
+// Dims returns the spatial dimensions of the curve, which is {2ᵏ, 2ᵏ}, where k
+// is the order.
+func (h Hilbert2D) Dims() []int { return []int{1 << h.order, 1 << h.order} }
+
+// Len returns the length of the curve, which is 2ⁿᵏ, where n is the dimension
+// (2) and k is the order.
+//
+// Len will overflow if order is ≥ 16 on architectures where [int] is 32-bits,
+// or ≥ 32 on architectures where [int] is 64-bits.
+func (h Hilbert2D) Len() int { return 1 << (2 * h.order) }
+
+func (h Hilbert2D) rot(n int, v []int, d int) {
+	switch d {
+	case 0:
+		swap{0, 1}.do(n, v)
+	case 3:
+		flip{0, 1}.do(n, v)
+	}
+}
+
+// Pos returns the linear position of the 3-spatial coordinate v along the
+// curve. Pos modifies v.
+//
+// Pos will overflow if order is ≥ 16 on architectures where [int] is 32-bits,
+// or ≥ 32 on architectures where [int] is 64-bits.
+func (h Hilbert2D) Pos(v []int) int {
+	var d int
+	const dims = 2
+	for n := h.order - 1; n >= 0; n-- {
+		rx := (v[0] >> n) & 1
+		ry := (v[1] >> n) & 1
+		rd := ry<<1 | (ry ^ rx)
+		d += rd << (dims * n)
+		h.rot(h.order, v, rd)
+	}
+	return d
+}
+
+// Coord writes the spatial coordinates of pos to dst and returns it. If dst is
+// nil, Coord allocates a new slice of length 2; otherwise Coord is
+// allocation-free.
+//
+// Coord panics if dst is not nil and len(dst) is not equal to 2.
+func (h Hilbert2D) Coord(dst []int, pos int) []int {
+	if dst == nil {
+		dst = make([]int, 2)
+	} else if len(dst) != 2 {
+		panic("len(dst) must equal 2")
+	}
+	for n := 0; n < h.order; n++ {
+		e := pos & 3
+		h.rot(n, dst[:], e)
+
+		ry := e >> 1
+		rx := (e>>0 ^ e>>1) & 1
+		dst[0] += rx << n
+		dst[1] += ry << n
+		pos >>= 2
+	}
+	return dst
+}
+
+// Hilbert3D is a 3-dimensional Hilbert curve.
+type Hilbert3D struct{ order int }
+
+// NewHilbert3D constructs a [Hilbert3D] of the given order. NewHilbert3D
+// returns [ErrOverflow] (wrapped) if the order would cause Len and Pos to
+// overflow.
+func NewHilbert3D(order int) (Hilbert3D, error) {
+	v := Hilbert3D{order: order}
+
+	// The order must be greater than zero.
+	if order < 1 {
+		return v, ErrUnderflow
+	}
+
+	// The product of the order and number of dimensions must not exceed or
+	// equal the size of an int.
+	if order*3 >= intSize {
+		return v, fmt.Errorf("a 3-dimensional, %d-order Hilbert curve will %w", order, ErrOverflow)
+	}
+
+	return v, nil
+}
+
+// Dims returns the spatial dimensions of the curve, which is {2ᵏ, 2ᵏ, 2ᵏ}, where
+// k is the order.
+func (h Hilbert3D) Dims() []int { return []int{1 << h.order, 1 << h.order, 1 << h.order} }
+
+// Len returns the length of the curve, which is 2ⁿᵏ, where n is the dimension
+// (3) and k is the order.
+//
+// Len will overflow if order is ≥ 11 on architectures where [int] is 32-bits,
+// or ≥ 21 on architectures where [int] is 64-bits.
+func (h Hilbert3D) Len() int { return 1 << (3 * h.order) }
+
+func (h Hilbert3D) rot(reverse bool, n int, v []int, d int) {
+	switch d {
+	case 0:
+		do2(reverse, n, v, swap{1, 2}, swap{0, 2})
+	case 1, 2:
+		do2(reverse, n, v, swap{0, 2}, swap{1, 2})
+	case 3, 4:
+		invert{0, 1}.do(n, v)
+	case 5, 6:
+		do2(reverse, n, v, flip{0, 2}, flip{1, 2})
+	case 7:
+		do2(reverse, n, v, flip{1, 2}, flip{0, 2})
+	}
+}
+
+// Pos returns the linear position of the 4-spatial coordinate v along the
+// curve. Pos modifies v.
+//
+// Pos will overflow if order is ≥ 11 on architectures where [int] is 32-bits,
+// or ≥ 21 on architectures where [int] is 64-bits.
+func (h Hilbert3D) Pos(v []int) int {
+	var d int
+	const dims = 3
+	for n := h.order - 1; n >= 0; n-- {
+		rx := (v[0] >> n) & 1
+		ry := (v[1] >> n) & 1
+		rz := (v[2] >> n) & 1
+		rd := rz<<2 | (rz^ry)<<1 | (rz ^ ry ^ rx)
+		d += rd << (dims * n)
+		h.rot(false, h.order, v, rd)
+	}
+	return d
+}
+
+// Coord writes the spatial coordinates of pos to dst and returns it. If dst is
+// nil, Coord allocates a new slice of length 3; otherwise Coord is
+// allocation-free.
+//
+// Coord panics if dst is not nil and len(dst) is not equal to 3.
+func (h Hilbert3D) Coord(dst []int, pos int) []int {
+	if dst == nil {
+		dst = make([]int, 3)
+	} else if len(dst) != 3 {
+		panic("len(dst) must equal 3")
+	}
+	for n := 0; n < h.order; n++ {
+		e := pos & 7
+		h.rot(true, n, dst[:], e)
+
+		rz := e >> 2
+		ry := (e>>1 ^ e>>2) & 1
+		rx := (e>>0 ^ e>>1) & 1
+		dst[0] += rx << n
+		dst[1] += ry << n
+		dst[2] += rz << n
+		pos >>= 3
+	}
+	return dst
+}
+
+// Hilbert4D is a 4-dimensional Hilbert curve.
+type Hilbert4D struct{ order int }
+
+// NewHilbert4D constructs a [Hilbert4D] of the given order. NewHilbert4D
+// returns [ErrOverflow] (wrapped) if the order would cause Len and Pos to
+// overflow.
+func NewHilbert4D(order int) (Hilbert4D, error) {
+	v := Hilbert4D{order: order}
+
+	// The order must be greater than zero.
+	if order < 1 {
+		return v, ErrUnderflow
+	}
+
+	// The product of the order and number of dimensions must not exceed or
+	// equal the size of an int.
+	if order*4 >= intSize {
+		return v, fmt.Errorf("a 4-dimensional, %d-order Hilbert curve will %w", order, ErrOverflow)
+	}
+
+	return v, nil
+}
+
+// Dims returns the spatial dimensions of the curve, which is {2ᵏ, 2ᵏ, 2ᵏ, 2ᵏ},
+// where k is the order.
+func (h Hilbert4D) Dims() []int { return []int{1 << h.order, 1 << h.order, 1 << h.order, 1 << h.order} }
+
+// Len returns the length of the curve, which is 2ⁿᵏ, where n is the dimension
+// (4) and k is the order.
+//
+// Len will overflow if order is ≥ 8 on architectures where [int] is 32-bits, or
+// ≥ 16 on architectures where [int] is 64-bits.
+func (h Hilbert4D) Len() int { return 1 << (4 * h.order) }
+
+func (h Hilbert4D) rot(reverse bool, n int, v []int, d int) {
+	switch d {
+	case 0:
+		do2(reverse, n, v, swap{1, 3}, swap{0, 3})
+	case 1, 2:
+		do2(reverse, n, v, swap{0, 3}, swap{1, 3})
+	case 3, 4:
+		do2(reverse, n, v, flip{0, 1}, swap{2, 3})
+	case 5, 6:
+		do2(reverse, n, v, flip{1, 2}, swap{2, 3})
+	case 7, 8:
+		invert{0, 2}.do(n, v)
+	case 9, 10:
+		do2(reverse, n, v, flip{1, 2}, flip{2, 3})
+	case 11, 12:
+		do2(reverse, n, v, flip{0, 1}, flip{2, 3})
+	case 13, 14:
+		do2(reverse, n, v, flip{0, 3}, flip{1, 3})
+	case 15:
+		do2(reverse, n, v, flip{1, 3}, flip{0, 3})
+	}
+}
+
+// Pos returns the linear position of the 2-spatial coordinate v along the
+// curve. Pos modifies v.
+//
+// Pos will overflow if order is ≥ 8 on architectures where [int] is 32-bits, or
+// ≥ 16 on architectures where [int] is 64-bits.
+func (h Hilbert4D) Pos(v []int) int {
+	var d int
+	const dims = 4
+	for n := h.order - 1; n >= 0; n-- {
+		var e int
+		for i := dims - 1; i >= 0; i-- {
+			v := v[i] >> n & 1
+			e = e<<1 | (e^v)&1
+		}
+
+		d += e << (dims * n)
+		h.rot(false, h.order, v, e)
+	}
+	return d
+}
+
+// Coord writes the spatial coordinates of pos to dst and returns it. If dst is
+// nil, Coord allocates a new slice of length 4; otherwise Coord is
+// allocation-free.
+//
+// Coord panics if dst is not nil and len(dst) is not equal to 4.
+func (h Hilbert4D) Coord(dst []int, pos int) []int {
+	if dst == nil {
+		dst = make([]int, 4)
+	} else if len(dst) != 4 {
+		panic("len(dst) must equal 4")
+	}
+	N := 4
+	for n := 0; n < h.order; n++ {
+		e := pos & (1<<N - 1)
+		h.rot(true, n, dst[:], e)
+
+		for i, e := 0, e; i < N; i++ {
+			dst[i] += (e ^ e>>1) & 1 << n
+			e >>= 1
+		}
+		pos >>= N
+	}
+	return dst
+}
+
+type op interface{ do(int, []int) }
+
+// invert I and J
+type invert struct{ i, j int }
+
+func (c invert) do(n int, v []int) { v[c.i], v[c.j] = v[c.i]^(1<<n-1), v[c.j]^(1<<n-1) }
+
+// swap I and J
+type swap struct{ i, j int }
+
+func (c swap) do(n int, v []int) { v[c.i], v[c.j] = v[c.j], v[c.i] }
+
+// swap and invert I and J
+type flip struct{ i, j int }
+
+func (c flip) do(n int, v []int) { v[c.i], v[c.j] = v[c.j]^(1<<n-1), v[c.i]^(1<<n-1) }
+
+// do2 executes the given operations, optionally in reverse.
+//
+// Generic specialization reduces allocation (because it can eliminate interface
+// value boxing) and improves performance
+func do2[A, B op](reverse bool, n int, v []int, a A, b B) {
+	if reverse {
+		b.do(n, v)
+		a.do(n, v)
+	} else {
+		a.do(n, v)
+		b.do(n, v)
+	}
+}
diff --git a/spatial/curve/hilbert_test.go b/spatial/curve/hilbert_test.go
new file mode 100644
index 0000000..b68b14e
--- /dev/null
+++ b/spatial/curve/hilbert_test.go
@@ -0,0 +1,457 @@
+// Copyright ©2024 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package curve
+
+import (
+	"errors"
+	"fmt"
+	"reflect"
+	"slices"
+	"strings"
+	"testing"
+
+	"golang.org/x/exp/rand"
+)
+
+func ExampleHilbert2D_Pos() {
+	h := Hilbert2D{order: 3}
+
+	for y := 0; y < 1<<h.order; y++ {
+		for x := 0; x < 1<<h.order; x++ {
+			if x > 0 {
+				fmt.Print("  ")
+			}
+			fmt.Printf("%02x", h.Pos([]int{x, y}))
+		}
+		fmt.Println()
+	}
+	// Output:
+	// 00  01  0e  0f  10  13  14  15
+	// 03  02  0d  0c  11  12  17  16
+	// 04  07  08  0b  1e  1d  18  19
+	// 05  06  09  0a  1f  1c  1b  1a
+	// 3a  39  36  35  20  23  24  25
+	// 3b  38  37  34  21  22  27  26
+	// 3c  3d  32  33  2e  2d  28  29
+	// 3f  3e  31  30  2f  2c  2b  2a
+}
+
+func ExampleHilbert3D_Pos() {
+	h := Hilbert3D{order: 2}
+
+	for z := 0; z < 1<<h.order; z++ {
+		for y := 0; y < 1<<h.order; y++ {
+			for x := 0; x < 1<<h.order; x++ {
+				if x > 0 {
+					fmt.Print("  ")
+				}
+				fmt.Printf("%02x", h.Pos([]int{x, y, z}))
+			}
+			fmt.Println()
+		}
+		fmt.Println()
+	}
+	// Output:
+	// 00  07  08  0b
+	// 01  06  0f  0c
+	// 1a  1b  10  13
+	// 19  18  17  14
+	//
+	// 03  04  09  0a
+	// 02  05  0e  0d
+	// 1d  1c  11  12
+	// 1e  1f  16  15
+	//
+	// 3c  3b  36  35
+	// 3d  3a  31  32
+	// 22  23  2e  2d
+	// 21  20  29  2a
+	//
+	// 3f  38  37  34
+	// 3e  39  30  33
+	// 25  24  2f  2c
+	// 26  27  28  2b
+}
+
+func ExampleHilbert4D_Pos() {
+	h := Hilbert4D{order: 2}
+	N := 1 << h.order
+	for z := 0; z < N; z++ {
+		if z > 0 {
+			s := strings.Repeat("═", N*4-2)
+			s = s + strings.Repeat("═╬═"+s, N-1)
+			fmt.Println(s)
+		}
+		for y := 0; y < N; y++ {
+			for w := 0; w < N; w++ {
+				if w > 0 {
+					fmt.Print(" ║ ")
+				}
+				for x := 0; x < N; x++ {
+					if x > 0 {
+						fmt.Print("  ")
+					}
+					fmt.Printf("%02x", h.Pos([]int{x, y, z, w}))
+				}
+			}
+			fmt.Println()
+		}
+	}
+	// Output:
+	// 00  0f  10  13 ║ 03  0c  11  12 ║ fc  f3  ee  ed ║ ff  f0  ef  ec
+	// 01  0e  1f  1c ║ 02  0d  1e  1d ║ fd  f2  e1  e2 ║ fe  f1  e0  e3
+	// 32  31  20  23 ║ 35  36  21  22 ║ ca  c9  de  dd ║ cd  ce  df  dc
+	// 33  30  2f  2c ║ 34  37  2e  2d ║ cb  c8  d1  d2 ║ cc  cf  d0  d3
+	// ═══════════════╬════════════════╬════════════════╬═══════════════
+	// 07  08  17  14 ║ 04  0b  16  15 ║ fb  f4  e9  ea ║ f8  f7  e8  eb
+	// 06  09  18  1b ║ 05  0a  19  1a ║ fa  f5  e6  e5 ║ f9  f6  e7  e4
+	// 3d  3e  27  24 ║ 3a  39  26  25 ║ c5  c6  d9  da ║ c2  c1  d8  db
+	// 3c  3f  28  2b ║ 3b  38  29  2a ║ c4  c7  d6  d5 ║ c3  c0  d7  d4
+	// ═══════════════╬════════════════╬════════════════╬═══════════════
+	// 76  77  6c  6d ║ 79  78  6b  6a ║ 86  87  94  95 ║ 89  88  93  92
+	// 75  74  63  62 ║ 7a  7b  64  65 ║ 85  84  9b  9a ║ 8a  8b  9c  9d
+	// 42  41  5c  5d ║ 45  46  5b  5a ║ ba  b9  a4  a5 ║ bd  be  a3  a2
+	// 43  40  53  52 ║ 44  47  54  55 ║ bb  b8  ab  aa ║ bc  bf  ac  ad
+	// ═══════════════╬════════════════╬════════════════╬═══════════════
+	// 71  70  6f  6e ║ 7e  7f  68  69 ║ 81  80  97  96 ║ 8e  8f  90  91
+	// 72  73  60  61 ║ 7d  7c  67  66 ║ 82  83  98  99 ║ 8d  8c  9f  9e
+	// 4d  4e  5f  5e ║ 4a  49  58  59 ║ b5  b6  a7  a6 ║ b2  b1  a0  a1
+	// 4c  4f  50  51 ║ 4b  48  57  56 ║ b4  b7  a8  a9 ║ b3  b0  af  ae
+}
+
+func TestConstructors(t *testing.T) {
+	const intSize = 32 << (^uint(0) >> 63) // 32 or 64
+
+	t.Run("2D/Ok", func(t *testing.T) {
+		_, err := NewHilbert2D(intSize/2 - 1)
+		noError(t, err)
+	})
+
+	t.Run("3D/Ok", func(t *testing.T) {
+		_, err := NewHilbert3D(intSize / 3)
+		noError(t, err)
+	})
+
+	t.Run("4D/Ok", func(t *testing.T) {
+		_, err := NewHilbert4D(intSize/4 - 1)
+		noError(t, err)
+	})
+
+	t.Run("2D/Underflow", func(t *testing.T) {
+		_, err := NewHilbert2D(0)
+		errorIs(t, err, ErrUnderflow)
+	})
+
+	t.Run("3D/Underflow", func(t *testing.T) {
+		_, err := NewHilbert3D(0)
+		errorIs(t, err, ErrUnderflow)
+	})
+
+	t.Run("4D/Underflow", func(t *testing.T) {
+		_, err := NewHilbert4D(0)
+		errorIs(t, err, ErrUnderflow)
+	})
+
+	t.Run("2D/Overflow", func(t *testing.T) {
+		_, err := NewHilbert2D(intSize / 2)
+		errorIs(t, err, ErrOverflow)
+	})
+
+	t.Run("3D/Overflow", func(t *testing.T) {
+		_, err := NewHilbert3D(intSize/3 + 1)
+		errorIs(t, err, ErrOverflow)
+	})
+
+	t.Run("4D/Overflow", func(t *testing.T) {
+		_, err := NewHilbert4D(intSize / 4)
+		errorIs(t, err, ErrOverflow)
+	})
+}
+
+func TestHilbert2D(t *testing.T) {
+	for ord := 1; ord <= 4; ord++ {
+		t.Run(fmt.Sprintf("Order/%d", ord), func(t *testing.T) {
+			testCurve(t, Hilbert2D{order: ord})
+		})
+	}
+
+	cases := map[int][]int{
+		1: {
+			0, 1,
+			3, 2,
+		},
+		2: {
+			0x0, 0x3, 0x4, 0x5,
+			0x1, 0x2, 0x7, 0x6,
+			0xE, 0xD, 0x8, 0x9,
+			0xF, 0xC, 0xB, 0xA,
+		},
+		3: {
+			0x00, 0x01, 0x0E, 0x0F, 0x10, 0x13, 0x14, 0x15,
+			0x03, 0x02, 0x0D, 0x0C, 0x11, 0x12, 0x17, 0x16,
+			0x04, 0x07, 0x08, 0x0B, 0x1E, 0x1D, 0x18, 0x19,
+			0x05, 0x06, 0x09, 0x0A, 0x1F, 0x1C, 0x1B, 0x1A,
+			0x3A, 0x39, 0x36, 0x35, 0x20, 0x23, 0x24, 0x25,
+			0x3B, 0x38, 0x37, 0x34, 0x21, 0x22, 0x27, 0x26,
+			0x3C, 0x3D, 0x32, 0x33, 0x2E, 0x2D, 0x28, 0x29,
+			0x3F, 0x3E, 0x31, 0x30, 0x2F, 0x2C, 0x2B, 0x2A,
+		},
+	}
+
+	for order, expected := range cases {
+		t.Run(fmt.Sprintf("Case/%d", order), func(t *testing.T) {
+			testCurveCase(t, Hilbert2D{order: order}, order, expected)
+		})
+	}
+}
+
+func TestHilbert3D(t *testing.T) {
+	for ord := 1; ord <= 4; ord++ {
+		t.Run(fmt.Sprintf("Order/%d", ord), func(t *testing.T) {
+			testCurve(t, Hilbert3D{order: ord})
+		})
+	}
+
+	cases := map[int][]int{
+		1: {
+			0, 1,
+			3, 2,
+
+			7, 6,
+			4, 5,
+		},
+		2: {
+			0x00, 0x07, 0x08, 0x0B,
+			0x01, 0x06, 0x0F, 0x0C,
+			0x1A, 0x1B, 0x10, 0x13,
+			0x19, 0x18, 0x17, 0x14,
+
+			0x03, 0x04, 0x09, 0x0A,
+			0x02, 0x05, 0x0E, 0x0D,
+			0x1D, 0x1C, 0x11, 0x12,
+			0x1E, 0x1F, 0x16, 0x15,
+
+			0x3C, 0x3B, 0x36, 0x35,
+			0x3D, 0x3A, 0x31, 0x32,
+			0x22, 0x23, 0x2E, 0x2D,
+			0x21, 0x20, 0x29, 0x2A,
+
+			0x3F, 0x38, 0x37, 0x34,
+			0x3E, 0x39, 0x30, 0x33,
+			0x25, 0x24, 0x2F, 0x2C,
+			0x26, 0x27, 0x28, 0x2B,
+		},
+	}
+
+	for order, expected := range cases {
+		t.Run(fmt.Sprintf("Case/%d", order), func(t *testing.T) {
+			testCurveCase(t, Hilbert3D{order: order}, order, expected)
+		})
+	}
+}
+
+func TestHilbert4D(t *testing.T) {
+	for ord := 1; ord <= 4; ord++ {
+		t.Run(fmt.Sprintf("Order %d", ord), func(t *testing.T) {
+			testCurve(t, Hilbert4D{order: ord})
+		})
+	}
+}
+
+func BenchmarkHilbert(b *testing.B) {
+	const O = 10
+	for N := 2; N <= 4; N++ {
+		b.Run(fmt.Sprintf("%dD/Pos", N), func(b *testing.B) {
+			for ord := 1; ord <= O; ord++ {
+				b.Run(fmt.Sprintf("Order %d", ord), func(b *testing.B) {
+					h := newCurve(ord, N)
+					v := make([]int, N)
+					for i := range v {
+						v[i] = rand.Intn(1 << ord)
+					}
+					u := make([]int, N)
+					for n := 0; n < b.N; n++ {
+						copy(u, v)
+						h.Pos(u)
+					}
+				})
+			}
+		})
+
+		b.Run(fmt.Sprintf("%dD/Coord", N), func(b *testing.B) {
+			for ord := 1; ord <= O; ord++ {
+				b.Run(fmt.Sprintf("Order %d", ord), func(b *testing.B) {
+					h := newCurve(ord, N)
+					d := rand.Intn(1 << (ord * N))
+					v := make([]int, N)
+					for n := 0; n < b.N; n++ {
+						h.Coord(v, d)
+					}
+				})
+			}
+		})
+	}
+}
+
+type curve interface {
+	Dims() []int
+	Len() int
+	Pos(v []int) int
+	Coord(dst []int, pos int) []int
+}
+
+func newCurve(order, dim int) curve {
+	switch dim {
+	case 2:
+		return Hilbert2D{order: order}
+	case 3:
+		return Hilbert3D{order: order}
+	case 4:
+		return Hilbert4D{order: order}
+	}
+	panic("invalid dimension")
+}
+
+// testCurve verifies that Pos and Coord (of C) are inverses of each other and
+// that the spatial coordinates V and U - corresponding to linear coordinates D
+// and D+1 - are exactly one unit (euclidean) distant from each other.
+func testCurve(t *testing.T, c curve) {
+	t.Helper()
+
+	// Stop if the error count reaches 10
+	var errc int
+	fail := func() {
+		if errc < 10 {
+			errc++
+			t.Fail()
+		} else {
+			t.FailNow()
+		}
+	}
+
+	// Map between linear and spatial coordinates, and verify that Pos and Coord
+	// are inverses of each other
+	m := map[int][]int{}
+	curveRange(c, func(v []int) {
+		d := c.Pos(slices.Clone(v))
+		u := c.Coord(nil, d)
+		if !reflect.DeepEqual(v, u) {
+			t.Logf("Space is not the inverse of Curve for d=%d %v", d, v)
+			fail()
+		}
+
+		m[d] = slices.Clone(v)
+	})
+
+	D := 1
+	for _, v := range c.Dims() {
+		D *= v
+	}
+
+	// For each possible pairs of linear coordinates D and D+1, verify that the
+	// corresponding spatial coordinates V and U are exactly one unit apart
+	// (euclidean distance)
+	for d := 0; d < D-1; d++ {
+		v, u := m[d], m[d+1]
+		if !adjacent(v, u) {
+			t.Logf("points %x and %x are not adjacent\n    %v -> %v", d, d+1, v, u)
+			fail()
+		}
+	}
+}
+
+// curveRange ranges over the n-dimensional coordinate space of the curve,
+// calling fn on each element of the space.
+func curveRange(c curve, fn func([]int)) {
+	size := c.Dims()
+	dimRange(len(size), size, make([]int, len(size)), fn)
+}
+
+// dimRange ranges over the coordinate space defined by size, calling fn on each
+// element of the space. Call dimRange with dim = len(size).
+func dimRange(dim int, size []int, v []int, fn func([]int)) {
+	if dim == 0 {
+		fn(v)
+		return
+	}
+
+	for i := 0; i < size[dim-1]; i++ {
+		v[dim-1] = i
+		dimRange(dim-1, size, v, fn)
+	}
+}
+
+// adjacent returns true if the euclidean distance between v and u is
+// exactly one. v and u must be the same length.
+//
+// In other words, given d(i) = abs(v[i] - u[i]), adjacent returns false if
+// d(i) > 1 for any i or if d(i) == 1 for more than a single i, or if d(i)
+// == 0 for all i.
+func adjacent(v, u []int) bool {
+	n := 0
+	for i := range v {
+		x := v[i] - u[i]
+		if x == 0 {
+			continue
+		}
+		if x < -1 || 1 < x {
+			return false
+		}
+		n++
+	}
+
+	return n == 1
+}
+
+// testCurveCase verifies that the curve produces the expected sequence of
+// values.
+func testCurveCase(t *testing.T, c curve, order int, expected []int) {
+	t.Helper()
+
+	dim := len(c.Dims())
+	actual := make([]int, len(expected))
+	for i := range expected {
+		v := coord(i, order, dim)
+		actual[i] = c.Pos(slices.Clone(v))
+	}
+	if !reflect.DeepEqual(expected, actual) {
+		t.Fatalf("unexpected result:\ngot:  %v\nwant: %v", actual, expected)
+	}
+
+	for i, d := range expected {
+		v := coord(i, order, dim)
+		if !reflect.DeepEqual(v, c.Coord(nil, d)) {
+			t.Fatalf("[%d] expected %v for d = %d", i, v, d)
+		}
+	}
+}
+
+// coord returns the nth spatial coordinates for a dim-dimensional space with
+// sides 2^ord in row-major order.
+func coord(n, ord, dim int) []int {
+	v := make([]int, dim)
+	for i := 0; i < dim; i++ {
+		v[i] = n % (1 << ord)
+		n /= (1 << ord)
+	}
+	return v
+}
+
+func noError(t testing.TB, err error) {
+	t.Helper()
+	if err != nil {
+		t.Fatal(err)
+	}
+}
+
+func errorIs(t testing.TB, err, target error) {
+	t.Helper()
+	if err == nil {
+		t.Fatal("Expected an error")
+	}
+	if !errors.Is(err, target) {
+		t.Fatalf("Expected %v to be %v", err, target)
+	}
+}