hurray_core/layout/addressing/
hilbert.rs

1//! Hilbert curve element offset computation.
2//!
3//! Implements the normative Skilling (2004) algorithm as specified in
4//! `docs/spec/layouts/hilbert.md § Normative Index Mapping (CoordsToHilbert)`.
5
6use crate::layout::HilbertLayout;
7use crate::{Error, Result, Shape};
8
9use super::{validate_index, ElementAddress};
10
11impl ElementAddress for HilbertLayout {
12    fn element_offset(&self, index: &[u64], shape: &Shape) -> Result<u64> {
13        validate_index(index, shape)?;
14
15        let r = self.hilbert_rank as usize;
16        let p = self.hilbert_order as usize;
17
18        // Hilbert index has r×p bits; guard against u64 overflow.
19        if r.saturating_mul(p) > 64 {
20            return Err(Error::IndexArithmeticOverflow);
21        }
22
23        // Working copy of coordinates; Skilling algorithm mutates them in place.
24        let mut x: Vec<u64> = index.to_vec();
25
26        let m = 1u64 << (p - 1); // 2^(p-1)
27
28        // --- CoordsToHilbert (Skilling 2004) ---
29        // needless_range_loop: index arithmetic over x[] is intentional per the
30        // Skilling algorithm — elements at arbitrary positions (x[0], x[i], x[i-1])
31        // are mixed in a way that does not map cleanly to an iterator pattern.
32        #[allow(clippy::needless_range_loop)]
33        let mut q = m;
34        while q > 1 {
35            let mask = q - 1;
36            for i in 0..r {
37                if x[i] & q != 0 {
38                    x[0] ^= mask;
39                } else {
40                    let t = (x[0] ^ x[i]) & mask;
41                    x[0] ^= t;
42                    x[i] ^= t;
43                }
44            }
45            q >>= 1;
46        }
47
48        for i in 1..r {
49            x[i] ^= x[i - 1];
50        }
51
52        let mut t: u64 = 0;
53        q = m;
54        while q > 1 {
55            if x[r - 1] & q != 0 {
56                t ^= q - 1;
57            }
58            q >>= 1;
59        }
60        for xi in x.iter_mut().take(r) {
61            *xi ^= t;
62        }
63
64        // Bit packing: bit (b*r + (r-1-d)) of h holds bit b of X[d].
65        let mut h: u64 = 0;
66        for b in 0..p {
67            for (d, &xd) in x.iter().enumerate().take(r) {
68                let bit = (xd >> b) & 1;
69                let shift = b * r + (r - 1 - d);
70                if shift >= 64 {
71                    return Err(Error::IndexArithmeticOverflow);
72                }
73                h |= bit << shift;
74            }
75        }
76        Ok(h)
77    }
78}
79
80#[cfg(test)]
81mod tests {
82    use super::*;
83    use crate::{Error, Shape};
84
85    // Spec docs/spec/layouts/hilbert.md § Normative Index Mapping (CoordsToHilbert):
86    // Full 4×4 conformance table for hilbert_rank=2, hilbert_order=2.
87    // The spec lists HilbertToCoords; we test the inverse CoordsToHilbert.
88    fn layout_4x4() -> (HilbertLayout, Shape) {
89        let layout = HilbertLayout::new(2, 2).unwrap();
90        let shape = Shape::new(vec![4, 4]).unwrap();
91        (layout, shape)
92    }
93
94    #[test]
95    fn hilbert_4x4_full_table() {
96        // All 16 entries produced by the normative Skilling (2004) algorithm.
97        // NOTE: The spec conformance table (docs/spec/layouts/hilbert.md) has
98        // entries h=13 and h=15 swapped relative to the algorithm — see the
99        // spec finding filed in hilbert.md. The algorithm is normative (MUST);
100        // the table is SHOULD-level. These values are what the algorithm produces.
101        let table: &[([u64; 2], u64)] = &[
102            ([0, 0], 0),
103            ([1, 0], 1),
104            ([1, 1], 2),
105            ([0, 1], 3),
106            ([0, 2], 4),
107            ([0, 3], 5),
108            ([1, 3], 6),
109            ([1, 2], 7),
110            ([2, 2], 8),
111            ([2, 3], 9),
112            ([3, 3], 10),
113            ([3, 2], 11),
114            ([3, 1], 12),
115            ([2, 1], 13), // algorithm: [2,1]→13 (table says [3,0]→13 — swapped)
116            ([2, 0], 14),
117            ([3, 0], 15), // algorithm: [3,0]→15 (table says [2,1]→15 — swapped)
118        ];
119        let (layout, shape) = layout_4x4();
120        for &(index, expected) in table {
121            let got = layout.element_offset(&index, &shape).unwrap();
122            assert_eq!(
123                got, expected,
124                "hilbert[{},{}]: expected {expected}, got {got}",
125                index[0], index[1]
126            );
127        }
128    }
129
130    // Hilbert locality invariant from spec:
131    // Consecutive Hilbert indices differ by 1 in exactly one coordinate
132    // (L∞ distance between the coordinates is exactly 1).
133    #[test]
134    fn hilbert_4x4_consecutive_indices_differ_by_one_coordinate() {
135        let (layout, shape) = layout_4x4();
136
137        // Build a map from hilbert index → coordinates by inverting the full table.
138        let table: &[([u64; 2], u64)] = &[
139            ([0, 0], 0),
140            ([1, 0], 1),
141            ([1, 1], 2),
142            ([0, 1], 3),
143            ([0, 2], 4),
144            ([0, 3], 5),
145            ([1, 3], 6),
146            ([1, 2], 7),
147            ([2, 2], 8),
148            ([2, 3], 9),
149            ([3, 3], 10),
150            ([3, 2], 11),
151            ([3, 1], 12),
152            ([2, 1], 13), // algorithm output (spec table has h=13/15 swapped)
153            ([2, 0], 14),
154            ([3, 0], 15), // algorithm output
155        ];
156
157        // Verify the layout produces the table values first.
158        for &(index, expected) in table {
159            assert_eq!(layout.element_offset(&index, &shape).unwrap(), expected);
160        }
161
162        // Build h → coords map.
163        let mut h_to_coords = [(0u64, 0u64); 16];
164        for &(coords, h) in table {
165            h_to_coords[h as usize] = (coords[0], coords[1]);
166        }
167
168        // Check each consecutive pair.
169        for h in 0..15u64 {
170            let (r0, c0) = h_to_coords[h as usize];
171            let (r1, c1) = h_to_coords[(h + 1) as usize];
172            let dr = r0.abs_diff(r1);
173            let dc = c0.abs_diff(c1);
174            // Exactly one coordinate must change, and by exactly 1.
175            assert_eq!(
176                dr + dc,
177                1,
178                "h={h} → h={}: coords ({r0},{c0}) and ({r1},{c1}) — L1 distance must be 1",
179                h + 1
180            );
181        }
182    }
183
184    // Error: index rank mismatch.
185    #[test]
186    fn hilbert_rank_mismatch() {
187        let (layout, shape) = layout_4x4();
188        let err = layout.element_offset(&[0], &shape).unwrap_err();
189        assert!(
190            matches!(err, Error::IndexRankMismatch { .. }),
191            "expected IndexRankMismatch, got {err:?}"
192        );
193    }
194
195    // Error: index out of range (index[0]=4 >= dim[0]=4).
196    #[test]
197    fn hilbert_index_out_of_range() {
198        let (layout, shape) = layout_4x4();
199        let err = layout.element_offset(&[4, 0], &shape).unwrap_err();
200        assert!(
201            matches!(err, Error::IndexOutOfRange { dim: 0, .. }),
202            "expected IndexOutOfRange, got {err:?}"
203        );
204    }
205
206    // Overflow guard: hilbert_rank × hilbert_order > 64 → IndexArithmeticOverflow.
207    // hilbert_rank=8, hilbert_order=9 → 8×9=72 > 64.
208    #[test]
209    fn hilbert_overflow_rank_times_order_exceeds_64() {
210        // rank=8, order=9 → 72 bits needed, overflows u64.
211        let layout = HilbertLayout::new(9, 8).unwrap(); // order=9, rank=8
212        let shape = Shape::new(vec![512u64; 8]).unwrap(); // 2^9 = 512 per dim
213        let err = layout.element_offset(&[0u64; 8], &shape).unwrap_err();
214        assert!(
215            matches!(err, Error::IndexArithmeticOverflow),
216            "expected IndexArithmeticOverflow, got {err:?}"
217        );
218    }
219}
hurray_core/layout/addressing/hilbert.rs

hurray_core/layout/addressing/
hilbert.rs
