spatial/curve: new package to constuct 2-, 3- and 4-D Hilbert curves
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)
+ }
+}
