code wiki / _hdl_build / nx_geo_hilbert.nx
nx_geo_hilbert.nx source
↩ module page · 143 lines · 6078 B
1// nx_geo_hilbert.nx -- GEO-005: Hilbert space-filling-curve cell ordering, the S2/H3 superiority
2// target. The geohash (GEO-002) interleaves lon/lat bits = a Z-ORDER (Morton) curve, which has
3// BOUNDARY DISCONTINUITIES: consecutive indices can map to cells far apart (a jump at every quadrant
4// seam), so a range-scan over an index window can miss/overshoot spatially-near cells. The Hilbert
5// curve's DEFINING PROPERTY: consecutive indices are ALWAYS spatially adjacent (Manhattan distance
6// exactly 1) -- no discontinuities. PURE-INTEGER (division/modulo, NO bitwise -- the geohash's
7// no-float discipline + dodges the shift/or miscompile landmine), exact + deterministic. A grid of
8// n=2^order cells per axis maps bijectively to [0, n*n) Hilbert indices.
9// EXCEED is MEASURED, not asserted: hilbert_max_jump==1 vs zorder_max_jump>1 on the SAME grid.
10import "nx_syscalls.nx"
11
12func gh_w(fd: i64, s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(fd, s, n); return 0 }
13func gh_n(fd: i64, v: i64) -> i64 {
14 let bb: *u8 = sys_mmap(28); var m: i64 = v; if m < 0 { sys_write(fd, "-" as *u8, 1); m = 0 - m }
15 let t: *u8 = sys_mmap(28); var k: i64 = 0; if m == 0 { t[0] = 48 as u8; k = 1 }
16 while m > 0 { t[k] = (48 + (m % 10)) as u8; m = m / 10; k = k + 1 }
17 var i: i64 = 0; while i < k { bb[i] = t[k - 1 - i]; i = i + 1 } sys_write(fd, bb, k); return 0
18}
19func gh_abs(a: i64) -> i64 { if a < 0 { return 0 - a } return a }
20// rotate/reflect a quadrant (in-place on the xb/yb boxes). size = the rotation grid size.
21func hil_rot(size: i64, xb: *i64, yb: *i64, rx: i64, ry: i64) -> i64 {
22 if ry == 0 {
23 if rx == 1 { xb[0] = size - 1 - xb[0]; yb[0] = size - 1 - yb[0] }
24 let t: i64 = xb[0]; xb[0] = yb[0]; yb[0] = t
25 }
26 return 0
27}
28// (x,y) -> Hilbert distance d on an n x n grid. xor-of-low-bit done arithmetically.
29func hil_xy2d(n: i64, x: i64, y: i64) -> i64 {
30 let xb: *i64 = sys_mmap(8) as *i64
31 let yb: *i64 = sys_mmap(8) as *i64
32 xb[0] = x; yb[0] = y
33 var d: i64 = 0
34 var s: i64 = n / 2
35 while s > 0 {
36 var rx: i64 = 0
37 if (xb[0] / s) % 2 == 1 { rx = 1 }
38 var ry: i64 = 0
39 if (yb[0] / s) % 2 == 1 { ry = 1 }
40 var contrib: i64 = 3 * rx // (3*rx) xor ry, rx,ry in {0,1}
41 if ry == 1 { contrib = contrib + 1 - 2 * (contrib % 2) } // flip low bit: 0->1, 3->2
42 d = d + s * s * contrib
43 hil_rot(n, xb, yb, rx, ry)
44 s = s / 2
45 }
46 return d
47}
48// Hilbert distance d -> (x,y) into boxes.
49func hil_d2xy(n: i64, d: i64, xb: *i64, yb: *i64) -> i64 {
50 xb[0] = 0; yb[0] = 0
51 var t: i64 = d
52 var s: i64 = 1
53 while s < n {
54 let rx: i64 = (t / 2) % 2
55 let ry: i64 = ((t % 2) + rx) % 2 // (t xor rx) low bit = (t%2 + rx)%2
56 hil_rot(s, xb, yb, rx, ry)
57 xb[0] = xb[0] + s * rx
58 yb[0] = yb[0] + s * ry
59 t = t / 4
60 s = s * 2
61 }
62 return 0
63}
64// Z-order (Morton) index z -> (x,y), the geohash's ordering, for the exceed baseline.
65func mort_z2xy(z: i64, order: i64, xb: *i64, yb: *i64) -> i64 {
66 xb[0] = 0; yb[0] = 0
67 var i: i64 = 0
68 var pz: i64 = 1 // 2^(2i)
69 var px: i64 = 1 // 2^i
70 while i < order {
71 let zx: i64 = (z / pz) % 2
72 let zy: i64 = (z / (pz * 2)) % 2
73 xb[0] = xb[0] + zx * px
74 yb[0] = yb[0] + zy * px
75 pz = pz * 4
76 px = px * 2
77 i = i + 1
78 }
79 return 0
80}
81func main(argc: i64, argv: *i64) -> i64 {
82 let order: i64 = 3
83 let n: i64 = 8 // 2^order
84 let cells: i64 = n * n
85 let xb: *i64 = sys_mmap(8) as *i64
86 let yb: *i64 = sys_mmap(8) as *i64
87
88 // 1. BIJECTION: every Hilbert index round-trips xy<->d over the whole grid.
89 var rt_ok: i64 = 1
90 var d: i64 = 0
91 while d < cells {
92 hil_d2xy(n, d, xb, yb)
93 if hil_xy2d(n, xb[0], yb[0]) != d { rt_ok = 0 }
94 d = d + 1
95 }
96
97 // 2. HILBERT locality: max Manhattan jump between consecutive indices (defining property = 1).
98 var hmax: i64 = 0
99 d = 0
100 while d < cells - 1 {
101 hil_d2xy(n, d, xb, yb)
102 let x1: i64 = xb[0]; let y1: i64 = yb[0]
103 hil_d2xy(n, d + 1, xb, yb)
104 let jump: i64 = gh_abs(x1 - xb[0]) + gh_abs(y1 - yb[0])
105 if jump > hmax { hmax = jump }
106 d = d + 1
107 }
108 // 3. Z-ORDER locality: the geohash baseline -- has boundary discontinuities (jumps > 1).
109 var zmax: i64 = 0
110 var z: i64 = 0
111 while z < cells - 1 {
112 mort_z2xy(z, order, xb, yb)
113 let x1: i64 = xb[0]; let y1: i64 = yb[0]
114 mort_z2xy(z + 1, order, xb, yb)
115 let jump: i64 = gh_abs(x1 - xb[0]) + gh_abs(y1 - yb[0])
116 if jump > zmax { zmax = jump }
117 z = z + 1
118 }
119
120 // verdict: bijection + Hilbert PERFECT locality (hmax==1) + MEASURED exceed over Z-order (hmax<zmax).
121 var ok: i64 = 1
122 if rt_ok != 1 { ok = 0 }
123 if hmax != 1 { ok = 0 }
124 if zmax <= 1 { ok = 0 } // Z-order MUST show discontinuity (else no baseline)
125 if hmax >= zmax { ok = 0 } // exceed must be measured, not assumed
126
127 gh_w(1, "GEOHILBERT n=" as *u8); gh_n(1, n); gh_w(1, " cells=" as *u8); gh_n(1, cells)
128 gh_w(1, " bijection=" as *u8); gh_n(1, rt_ok)
129 gh_w(1, " hilbert_max_jump=" as *u8); gh_n(1, hmax)
130 gh_w(1, " zorder_max_jump=" as *u8); gh_n(1, zmax); gh_w(1, "\n" as *u8)
131 let lf: i64 = sys_openat_append("knowledge/status/geo_hilbert.log" as *u8, 420)
132 if lf >= 0 {
133 gh_w(lf, "GEOHILBERT GEO-005 bijection=" as *u8); gh_n(lf, rt_ok)
134 gh_w(lf, " hilbert_max_jump=" as *u8); gh_n(lf, hmax)
135 gh_w(lf, " zorder_max_jump=" as *u8); gh_n(lf, zmax)
136 gh_w(lf, " exceed=hilbert-no-boundary-discontinuity" as *u8)
137 if ok == 1 { gh_w(lf, " verdict=GREEN\n" as *u8) } else { gh_w(lf, " verdict=RED\n" as *u8) }
138 sys_close(lf)
139 }
140 if ok == 1 { sys_exit(0) }
141 sys_exit(1)
142 return 1
143}