nx_geo_eqarea.nx source

1// nx_geo_eqarea.nx -- GEO-005 (other half): EQUAL-AREA cells, the fix for the geohash's CELL-SIZE 2// weakness. A lat/lon grid (geohash) partitions latitude by equal Δlat, but a cell's GROUND AREA = 3// Δlon * (sin(lat_hi) - sin(lat_lo)) shrinks by cos(lat) toward the poles -- so a "uniform" geohash 4// grid actually has cells ~7x smaller near the poles than at the equator (density bias). The 5// CYLINDRICAL-EQUAL-AREA fix: project lat -> y = sin(lat) and partition by equal Δy. Because the 6// spherical area element is exactly dlon * d(sin lat), equal Δlon*Δsin rectangles have EQUAL ground 7// area at every latitude. Composes the Hilbert ordering (nx_geo_hilbert) -> equal-area cells in 8// Hilbert order = the S2/H3 cell-uniformity target. REUSES fx.nx CORDIC sin (no float, no reinvention). 9// EXCEED is MEASURED: lat/lon band areas vary >2x (via real fx_sin) while equal-area bands are uniform. 10import "fx.nx" 11import "nx_itoa_lib.nx" // shared MSB-first emitter (zero-alloc) 12import "nx_syscalls.nx" 13 14func ea_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 } 15// MIGRATED to the shared emitter (debt 1785563586). The old body mmapped a scratch buffer 16// per call and never freed it. At PAGE granularity that is 4096B leaked PER CALL -- the 17// defect that took 28.5GB of a 36GB host in nx_ts_lumadiff (2MB input, ~3.66M calls). 18// nxi_* is MSB-first, allocates NOTHING, and emits identical bytes including the sign. 19func ea_n(fd: i64, v: i64) -> i64 { nxi_fd(fd, v); return 0 } 20// full-degree angle -> Q16.16 radians (multiply before divide; full angles don't underflow Q16.16). 21func ea_deg_q16(deg: i64) -> i64 { return deg * FX_TWO_PI / 360 } 22func main(argc: i64, argv: *i64) -> i64 { 23 let N: i64 = 6 // bands over the 0..90 deg quadrant (sin 0..1) 24 let ll: *i64 = sys_mmap(N * 8) as *i64 // lat/lon (geohash) band ground-areas (= Δsin via fx_sin) 25 let ea: *i64 = sys_mmap(N * 8) as *i64 // equal-area band ground-areas (= equal Δsin) 26 var k: i64 = 0 27 while k < N { 28 // lat/lon: equal Δlat partition; ground area = Δlon * (sin(lat_hi) - sin(lat_lo)), Δlon=1. 29 let lat_lo: i64 = k * 90 / N 30 let lat_hi: i64 = (k + 1) * 90 / N 31 ll[k] = fx_sin(ea_deg_q16(lat_hi)) - fx_sin(ea_deg_q16(lat_lo)) 32 // equal-area: equal Δ(sin lat) partition; ground area = Δlon * Δsin = FX_ONE/N (uniform). 33 let sin_lo: i64 = k * FX_ONE / N 34 let sin_hi: i64 = (k + 1) * FX_ONE / N 35 ea[k] = sin_hi - sin_lo 36 k = k + 1 37 } 38 var ll_max: i64 = ll[0]; var ll_min: i64 = ll[0] 39 var ea_max: i64 = ea[0]; var ea_min: i64 = ea[0] 40 k = 1 41 while k < N { 42 if ll[k] > ll_max { ll_max = ll[k] } 43 if ll[k] < ll_min { ll_min = ll[k] } 44 if ea[k] > ea_max { ea_max = ea[k] } 45 if ea[k] < ea_min { ea_min = ea[k] } 46 k = k + 1 47 } 48 // verdict: lat/lon cells MUST vary >2x (the measured cos-shrinkage weakness; also proves fx_sin 49 // is really cos-weighted -- a broken/constant sin makes them uniform -> RED), AND equal-area cells 50 // MUST be uniform (max==min -- the fix). No-false-green on both sides. 51 let ll_spread: i64 = ll_max - ll_min 52 let ea_spread: i64 = ea_max - ea_min // 0 in real arithmetic; <=1 from FX_ONE/N integer rounding 53 var ok: i64 = 1 54 if ll_max <= ll_min { ok = 0 } 55 if ll_max < 2 * ll_min { ok = 0 } 56 if ea_spread > 1 { ok = 0 } // equal-area uniform within integer-rounding tolerance 57 if (ea_spread * 1000) >= ll_spread { ok = 0 } // MEASURED: equal-area >1000x more uniform than lat/lon 58 ea_w(1, "GEOEQAREA N=" as *u8); ea_n(1, N) 59 ea_w(1, " latlon_area_max=" as *u8); ea_n(1, ll_max); ea_w(1, " latlon_area_min=" as *u8); ea_n(1, ll_min) 60 ea_w(1, " eqarea_max=" as *u8); ea_n(1, ea_max); ea_w(1, " eqarea_min=" as *u8); ea_n(1, ea_min); ea_w(1, "\n" as *u8) 61 let lf: i64 = sys_openat_append("knowledge/status/geo_eqarea.log" as *u8, 420) 62 if lf >= 0 { 63 ea_w(lf, "GEOEQAREA GEO-005-eqcells latlon_area_max=" as *u8); ea_n(lf, ll_max) 64 ea_w(lf, " latlon_area_min=" as *u8); ea_n(lf, ll_min) 65 ea_w(lf, " latlon_spread=" as *u8); ea_n(lf, ll_spread) 66 ea_w(lf, " eqarea_spread=" as *u8); ea_n(lf, ea_spread) 67 ea_w(lf, " exceed=equal-ground-area-vs-geohash-pole-shrinkage" as *u8) 68 if ok == 1 { ea_w(lf, " verdict=GREEN\n" as *u8) } else { ea_w(lf, " verdict=RED\n" as *u8) } 69 sys_close(lf) 70 } 71 if ok == 1 { sys_exit(0) } 72 sys_exit(1) 73 return 1 74}