sketch_hll.nx source

1// sketch_hll.nx -- HyperLogLog cardinality sketch in NishiLang. 2// 3// Flajolet 2007 + Heule 2013 small-range correction. Sovereign- 4// tier sibling port of nishi-engine/packages/core-sketch/src/hll.ts. 5// All-i64 implementation: no f64 source-level types, no float 6// literals; harmonic-mean estimator computed in Q32.32 fixed-point 7// so the file compiles through both the C-anchor (Wheeler comparator) 8// and the NishiLang sibling. 9// 10// Per the lossless-language meta-cardinal (nishi-engine doc 20): 11// the cardinality estimate ships in an Approximate<i64> envelope 12// declaring stddev_rel = 1.04/sqrt(m) at confidence 0.6827, with 13// MaturityClass = ReferenceImpl and AdversarialSafety = Honest. 14// Bounded-loss-by-typing, not silent-loss. 15// 16// Memory: 1 byte per register. lgK=11 -> 2048 bytes; lgK=12 -> 17// 4096 bytes. Standard error: 1.04 / sqrt(2^lgK). 18// 19// Roadmap citations: 20// doc 19 -- DataSketches stomp roadmap 21// doc 20 -- universal lossless language 22// doc 21 -- research partnership not isolation 23 24import "syscalls.nx" 25import "murmur3.nx" 26import "nx_xxhash.nx" 27import "nx_bit.nx" 28import "nx_bits.nx" 29import "sketch_types.nx" 30import "hll_bias_table.nx" 31import "nx_i128.nx" 32 33// lgK in [4, 10] keeps alpha_m_sq * 2^32 comfortably in signed i64 34// (max signed i64 ~9.22e18; lgK=10 puts us at ~3.23e18). lgK >= 11 35// would overflow with this Q32.32 scaling; queued for a v2 estimator 36// that uses a narrower fixed-point format (e.g. Q16.48) for the 37// upper precision range. Practical sketches today use lgK 7-10 38// (1.04/sqrt(m): 9.2% to 3.3% rel-stddev). 39const NX_HLL_LGK_MIN: i64 = 4 40const NX_HLL_LGK_MAX: i64 = 12 // bits-up scale 2026-05-20 (was 10). 41 // Adds lg_k 11 (m=2048, ~2.3% rel err) 42 // and lg_k 12 (m=4096, ~1.6% rel err). 43 // Unblocks matched-memory comparisons 44 // in benches against Theta / KMV. 45 46// Murmur3 seed-mixing constants. Source: original murmur3 paper 47// (Appleby 2008) for hash-family de-correlation. Same values used by 48// DataSketches (apache/datasketches-java HllUtil.java). 49const NX_HLL_SEED_HI: i64 = 0x9747B28C 50const NX_HLL_SEED_LO: i64 = 0x36185EC0 51 52// Math constants used in the LC formula. Each named with its source. 53// ln(2) ≈ 0.693147180559945... -- math constant (lineage: log_function) 54// PPM_SCALE = 1_000_000 -- our fixed-point scale (1e6 = "parts per million") 55// We keep ln(2) at micro-precision (PPM); going finer needs a wider 56// intermediate (Q32+) which we defer to v2. 57const NX_HLL_PPM_SCALE: i64 = 1000000 58const NX_HLL_LN2_PPM: i64 = 693147 // ln(2) * 1e6 59 60// DS composite-estimator LC peg. Used to smooth the LC -> bias-corrected 61// transition: if bias_corrected_est <= (LC_PEG_NUM / LC_PEG_DEN) * K, the 62// estimator takes min(bias_corrected, linear_counting). Source: 63// DataSketches HllUtil 0.7 constant; rationale per Heule 2013 §4. 64const NX_HLL_LC_PEG_NUM: i64 = 7 65const NX_HLL_LC_PEG_DEN: i64 = 10 66 67// HLL handle. 24 bytes. regs is a separately-allocated m-byte 68// register array. 69struct Hll { 70 regs: *u8, 71 lg_k: i64, 72 m: i64, 73 seed: i64, 74} 75 76// === leading-zero helpers ========================================== 77// 78// Delegated to nx_bits_clz32 (dispatches to bsr+xor on x86_64, clzw 79// Zbb on rv64). Same 0..32 semantics; one machine instruction in 80// the hot HLL rho path instead of 5 branches. 81 82func _sketch_hll_clz32_superseded(x: i64) -> i64 { 83 return nx_bits_clz32(x) 84} 85 86// === construction ================================================= 87 88func nx_hll_alloc(lg_k: i64, seed: i64) -> *Hll { 89 if lg_k < NX_HLL_LGK_MIN { return 0 as *Hll } 90 if lg_k > NX_HLL_LGK_MAX { return 0 as *Hll } 91 let m: i64 = 1 << lg_k 92 let raw: *u8 = sys_mmap(24) 93 let h: *Hll = raw as *Hll 94 h.regs = sys_mmap(m) 95 h.lg_k = lg_k 96 h.m = m 97 h.seed = seed 98 var i: i64 = 0 99 while i < m { 100 h.regs[i] = 0 101 i = i + 1 102 } 103 return h 104} 105 106// === add =========================================================== 107// 108// Hash via 2x murmur3_32 with derived seeds (legacy path) -- gives a 109// 64-bit hash domain spliced from two 32-bit murmur3 outputs. Top 110// lg_k bits of hi = bucket index; suffix-then-lo drives rho. 111// 112// 2026-05-10 SELF-REMEDIATION TRACE (worked example for the 113// SELF_REMEDIATION_AS_MATHEMATICIAN cardinal): 114// Hypothesis: swap to xxh64 (single 64-bit, SMHasher-passing) 115// would tighten small-N accuracy without regressing high-N. 116// Phase 5 A/B: small-N tightened (n=100/200 became TIES vs DS) 117// BUT n=2000 regressed (mean 108 -> 142, lost BEATS). 118// Phase 6 generalization proof: REGRESSED at n=2000 cell. 119// Verdict: swap rejected. REVERTED to 2x murmur3_32. 120// 121// Root cause hypothesis (deferred until proper diagnosis): 122// xxh64 and 2x murmur3_32 both pass SMHasher. Difference at high 123// N is in higher-moment behavior that SMHasher's chi-squared 124// doesn't capture, or in alpha_m_sq_q64 calibration that has 125// subtle dependence on hash family's 4th-moment behavior. Queued 126// IMPROVEMENT_OPP: instrument register-distribution histograms 127// per hash family + per N to identify which moment diverges. 128// 129// genealogy_id: appleby_murmur3 (Austin Appleby) 130func nx_hll_add(h: *Hll, key: *u8, len: i64) -> i64 { 131 let seed_hi: i64 = h.seed ^ NX_HLL_SEED_HI 132 let seed_lo: i64 = h.seed ^ NX_HLL_SEED_LO 133 let hi: i64 = murmur3_32(seed_hi, key, len) & 0xFFFFFFFF 134 let lo: i64 = murmur3_32(seed_lo, key, len) & 0xFFFFFFFF 135 136 let idx: i64 = hi >> (32 - h.lg_k) 137 138 let upper: i64 = (hi << h.lg_k) & 0xFFFFFFFF 139 var r: i64 = 0 140 if upper != 0 { 141 r = nx_clz32(upper) + 1 142 } 143 if upper == 0 { 144 if lo != 0 { 145 r = (32 - h.lg_k) + nx_clz32(lo) + 1 146 } 147 if lo == 0 { 148 r = (64 - h.lg_k) + 1 149 } 150 } 151 152 if r > 255 { r = 255 } 153 154 let cur: i64 = h.regs[idx] 155 if cur < r { 156 h.regs[idx] = r 157 } 158 return 0 159} 160 161// === estimate (Q32.32 fixed-point harmonic mean) =================== 162// 163// E_raw = alpha_m * m^2 / sum(2^-r_j) 164// Stored as Q32.32: alpha_m_sq_q64 = alpha_m * m^2 * 2^32. 165// sum_q32 = sum_j ((1 << 32) >> r_j) (with 0 for r_j > 32) 166// estimate = alpha_m_sq_q64 / sum_q32 (i64 cardinality) 167// 168// alpha_m for lg_k >= 7 uses Flajolet's continuous approximation: 169// alpha_m = 0.7213 / (1 + 1.079 / m) 170// Lower lg_k uses tabulated values: 16->0.673, 32->0.697, 64->0.709. 171// 172// Hardcoded alpha_m_sq_q64 for lg_k in [4, 15]; computed offline 173// to avoid any floating-point operation in the hot path. 174 175// BUG FIX 2026-05-12: prior values were 877-995x too large, giving 1000x 176// over-estimates at n >> m where the classical HLL formula kicks in. 177// Smoke gates never caught this because they only exercised small-N where 178// the linear-counting correction takes over. Verified against DataSketches 179// Python 5.2.0 on n=1000 workload (was 1,137,371; should be ~1000). 180// 181// Correct formula: alpha_m * m^2 * 2^32 where alpha_m per Heule 2013: 182// m=16: alpha = 0.673 183// m=32: alpha = 0.697 184// m=64: alpha = 0.709 185// m>=128: alpha = 0.7213 / (1 + 1.079/m) 186func nx_hll_alpha_m_sq_q64(lg_k: i64) -> i64 { 187 if lg_k == 4 { return 739971325493 } // 0.673 * 16^2 * 2^32 188 if lg_k == 5 { return 3065438418239 } // 0.697 * 32^2 * 2^32 189 if lg_k == 6 { return 12472859905490 } // 0.709 * 64^2 * 2^32 190 if lg_k == 7 { return 50332686358312 } // alpha_128 * 128^2 * 2^32 191 if lg_k == 8 { return 202175761456818 } // alpha_256 * 256^2 * 2^32 192 if lg_k == 9 { return 810403740235830 } // alpha_512 * 512^2 * 2^32 193 if lg_k == 10 { return 3245027090684401 } // alpha_1024 * 1024^2 * 2^32 194 if lg_k == 11 { return 12986945521737604 } // alpha_2048 * 2048^2 * 2^32 195 // = 4 * lg_k=10 entry * ratio 196 // ratio = (1025.079/1024) * (2048/2049.079) 197 // = 1.0005267 198 if lg_k == 12 { return 51961470086950416 } // alpha_4096 * 4096^2 * 2^32 199 // = 4 * lg_k=11 entry * ratio 200 // ratio = (2049.079/2048) * (4096/4097.079) 201 // = 1.0002634 202 return 0 203} 204 205func nx_hll_pow2_neg_q32(r: i64) -> i64 { 206 if r >= 32 { return 0 } 207 if r <= 0 { return 1 << 32 } 208 return (1 << 32) >> r 209} 210 211// Per-lg_k thresholds for the LC-vs-classical-HLL switch. 212// Source: Heule 2013 supplementary tables (genealogy_id: heule_2013; 213// lineage_id: heule_bias_table). Each value is the maximum cardinality 214// at which the linear-counting formula is more accurate than raw HLL 215// for that lg_k. Below the threshold: use LC. Above: use raw HLL. 216const NX_HLL_HEULE_THRESH_LGK4: i64 = 10 217const NX_HLL_HEULE_THRESH_LGK5: i64 = 20 218const NX_HLL_HEULE_THRESH_LGK6: i64 = 40 219const NX_HLL_HEULE_THRESH_LGK7: i64 = 80 220const NX_HLL_HEULE_THRESH_LGK8: i64 = 220 221const NX_HLL_HEULE_THRESH_LGK9: i64 = 400 222const NX_HLL_HEULE_THRESH_LGK10: i64 = 900 223const NX_HLL_HEULE_THRESH_LGK11: i64 = 1800 // bits-up extrapolation (2x trend) 224const NX_HLL_HEULE_THRESH_LGK12: i64 = 3100 // Heule paper canonical for lg_k=12 225 226func nx_hll_heule_threshold(lg_k: i64) -> i64 { 227 if lg_k == 4 { return NX_HLL_HEULE_THRESH_LGK4 } 228 if lg_k == 5 { return NX_HLL_HEULE_THRESH_LGK5 } 229 if lg_k == 6 { return NX_HLL_HEULE_THRESH_LGK6 } 230 if lg_k == 7 { return NX_HLL_HEULE_THRESH_LGK7 } 231 if lg_k == 8 { return NX_HLL_HEULE_THRESH_LGK8 } 232 if lg_k == 9 { return NX_HLL_HEULE_THRESH_LGK9 } 233 if lg_k == 10 { return NX_HLL_HEULE_THRESH_LGK10 } 234 if lg_k == 11 { return NX_HLL_HEULE_THRESH_LGK11 } 235 if lg_k == 12 { return NX_HLL_HEULE_THRESH_LGK12 } 236 return 0 237} 238 239// log2(x) returned in PPM (parts per million), with fractional precision. 240// 241// Prior code computed integer m/zeros then floor-log2, which truncated 242// the fractional log info (e.g., log2(128/59) actually = 1.117 but 243// integer-divide gave 128/59=2 then floor-log2(2)=1 -- losing 0.117). 244// 245// This implementation keeps the fractional mantissa. For x = 2^k * (1 + frac): 246// log2(x) ≈ k + frac (linear mantissa interpolation) 247// Max error in log2 ~ 0.04 around mid-mantissa. This translates to 248// ~3% in the LC estimate -- much better than the prior ~12% from 249// dropping fractional bits entirely. 250// 251// genealogy_id: ours (linear approximation is our pragmatic choice; 252// DS uses double-precision native log2) 253// lineage_id: log_function 254func nx_hll_log2_ppm(x: i64) -> i64 { 255 if x <= 0 { return 0 } 256 if x == 1 { return 0 } 257 var k: i64 = 0 258 var t: i64 = x 259 while t > 1 { t = t >> 1; k = k + 1 } 260 let two_k: i64 = 1 << k 261 let frac_ppm: i64 = ((x - two_k) * NX_HLL_PPM_SCALE) / two_k 262 return k * NX_HLL_PPM_SCALE + frac_ppm 263} 264 265// Heule 2013 bias-correction lookup against the EMPIRICAL table 266// (hll_bias_table.nx, calibrated for OUR hash family — 2x murmur3_32 — 267// rather than the DS canonical table which was calibrated for xxhash64). 268// Table format: interleaved (raw_est_sample, corrected_n) pairs per lg_k. 269// 270// Why empirical, not DS canonical: DS's CompositeInterpolationXTable 271// expects DS's hash bias distribution. Plugging our hash family into 272// that table caused 10x worse error at small N during the 2026-05-10 273// grid-test session. Empirical regen with our actual HLL pipeline gives 274// table values matched to our hash characteristics. 275// 276// Cubic Lagrange interpolation uses nx_i128 (N2 numeric tier) for the 277// intermediate (num << 24) / den ratio to avoid i64 overflow at lg_k=10. 278// genealogy_id: heule_2013 + ours (empirical regen) + knuth_taocp_vol2 (long div for i128) 279// lineage_id: heule_bias_table 280 281func nx_hll_bias_table_get(lg_k: i64, i: i64) -> i64 { 282 if lg_k == 4 { return nx_hll_bias_lgk4(i) } 283 if lg_k == 5 { return nx_hll_bias_lgk5(i) } 284 if lg_k == 6 { return nx_hll_bias_lgk6(i) } 285 if lg_k == 7 { return nx_hll_bias_lgk7(i) } 286 if lg_k == 8 { return nx_hll_bias_lgk8(i) } 287 if lg_k == 9 { return nx_hll_bias_lgk9(i) } 288 if lg_k == 10 { return nx_hll_bias_lgk10(i) } 289 return -1 290} 291 292func nx_hll_bias_table_n(lg_k: i64) -> i64 { 293 if lg_k == 4 { return NX_HLL_BIAS_LGK4_N } 294 if lg_k == 5 { return NX_HLL_BIAS_LGK5_N } 295 if lg_k == 6 { return NX_HLL_BIAS_LGK6_N } 296 if lg_k == 7 { return NX_HLL_BIAS_LGK7_N } 297 if lg_k == 8 { return NX_HLL_BIAS_LGK8_N } 298 if lg_k == 9 { return NX_HLL_BIAS_LGK9_N } 299 if lg_k == 10 { return NX_HLL_BIAS_LGK10_N } 300 return 0 301} 302 303// Cubic Lagrange interpolation against the empirical table. 304// Uses nx_i128 (N2 numeric tier) for the (num << 24) / den intermediate 305// to avoid i64 overflow at the larger lg_k values: num can reach 306// ~1e9 in cardinality units; num << 24 = ~1.7e16 fits comfortably here 307// but i128 keeps us safe under expanded ranges. 308// 309// genealogy_id: heule_2013 + ours (empirical regen); 310// knuth_taocp_vol2 (long division) for the i128 backend 311// lineage_id: heule_bias_table 312 313const NX_HLL_LAGRANGE_Q: i64 = 24 // 24-bit fractional Q-format 314const NX_HLL_LAGRANGE_Q_SCALE: i64 = 16777216 // 1 << 24 315 316func nx_hll_bias_correct(lg_k: i64, raw_est: i64) -> i64 { 317 let n_samples: i64 = nx_hll_bias_table_n(lg_k) 318 if n_samples <= 0 { return raw_est } 319 320 // Bracket search. Table is interleaved (X, Y, X, Y, ...) -- X at 321 // even indices, Y at odd. n_samples is the count of (X, Y) pairs. 322 var i: i64 = 0 323 var found: i64 = -1 324 while i < n_samples - 1 { 325 let x_lo: i64 = nx_hll_bias_table_get(lg_k, i * 2) 326 let x_hi: i64 = nx_hll_bias_table_get(lg_k, (i + 1) * 2) 327 if x_lo <= raw_est { 328 if raw_est <= x_hi { found = i } 329 } 330 i = i + 1 331 } 332 if found < 0 { return raw_est } 333 334 // Pick 4 surrounding points; clamp to table edges. 335 var i0: i64 = found - 1 336 var i1: i64 = found 337 var i2: i64 = found + 1 338 var i3: i64 = found + 2 339 if i0 < 0 { 340 i0 = 0 341 i1 = 1 342 i2 = 2 343 i3 = 3 344 } 345 if i3 >= n_samples { 346 i3 = n_samples - 1 347 i2 = n_samples - 2 348 i1 = n_samples - 3 349 i0 = n_samples - 4 350 } 351 let x0: i64 = nx_hll_bias_table_get(lg_k, i0 * 2) 352 let y0: i64 = nx_hll_bias_table_get(lg_k, i0 * 2 + 1) 353 let x1: i64 = nx_hll_bias_table_get(lg_k, i1 * 2) 354 let y1: i64 = nx_hll_bias_table_get(lg_k, i1 * 2 + 1) 355 let x2: i64 = nx_hll_bias_table_get(lg_k, i2 * 2) 356 let y2: i64 = nx_hll_bias_table_get(lg_k, i2 * 2 + 1) 357 let x3: i64 = nx_hll_bias_table_get(lg_k, i3 * 2) 358 let y3: i64 = nx_hll_bias_table_get(lg_k, i3 * 2 + 1) 359 let dx0: i64 = raw_est - x0 360 let dx1: i64 = raw_est - x1 361 let dx2: i64 = raw_est - x2 362 let dx3: i64 = raw_est - x3 363 364 // L_k(x) = product over j!=k of (x - x_j) / (x_k - x_j) 365 // For each k compute num_k = product of three (raw - x_j) (i64); 366 // den_k = product of three (x_k - x_j) (i64); then the ratio in 367 // Q24 fixed-point via i128 intermediate (handles num << 24 overflow). 368 let num0: i64 = dx1 * dx2 * dx3 369 let den0: i64 = (x0 - x1) * (x0 - x2) * (x0 - x3) 370 let num1: i64 = dx0 * dx2 * dx3 371 let den1: i64 = (x1 - x0) * (x1 - x2) * (x1 - x3) 372 let num2: i64 = dx0 * dx1 * dx3 373 let den2: i64 = (x2 - x0) * (x2 - x1) * (x2 - x3) 374 let num3: i64 = dx0 * dx1 * dx2 375 let den3: i64 = (x3 - x0) * (x3 - x1) * (x3 - x2) 376 if den0 == 0 { return raw_est } 377 if den1 == 0 { return raw_est } 378 if den2 == 0 { return raw_est } 379 if den3 == 0 { return raw_est } 380 381 // coef_q24 = (num << 24) / den via i128 intermediate. 382 let c0_q24: i64 = nx_mulshl_div_i64(num0, 1, NX_HLL_LAGRANGE_Q, den0) 383 let c1_q24: i64 = nx_mulshl_div_i64(num1, 1, NX_HLL_LAGRANGE_Q, den1) 384 let c2_q24: i64 = nx_mulshl_div_i64(num2, 1, NX_HLL_LAGRANGE_Q, den2) 385 let c3_q24: i64 = nx_mulshl_div_i64(num3, 1, NX_HLL_LAGRANGE_Q, den3) 386 387 // term = y * coef_q24, descale by >> 24. Sum then round. 388 let sum_q24: i64 = y0 * c0_q24 + y1 * c1_q24 + y2 * c2_q24 + y3 * c3_q24 389 return sum_q24 / NX_HLL_LAGRANGE_Q_SCALE 390} 391 392// Returns (estimate, zeros_count). Caller wraps in Approximate<i64> 393// via nx_hll_query for the typed-envelope path. 394func nx_hll_estimate(h: *Hll) -> i64 { 395 let m: i64 = h.m 396 var sum_q32: i64 = 0 397 var zeros: i64 = 0 398 var i: i64 = 0 399 while i < m { 400 let r: i64 = h.regs[i] 401 if r == 0 { zeros = zeros + 1 } 402 sum_q32 = sum_q32 + nx_hll_pow2_neg_q32(r) 403 i = i + 1 404 } 405 if sum_q32 == 0 { return m } // every register saturated; degenerate 406 let alpha_m_sq_q64: i64 = nx_hll_alpha_m_sq_q64(h.lg_k) 407 if alpha_m_sq_q64 == 0 { 408 // lg_k outside hardcoded range -- fall back to crude estimate 409 return zeros 410 } 411 let raw_est: i64 = alpha_m_sq_q64 / sum_q32 412 413 // DS composite estimator (Heule 2013 + DataSketches refinement): 414 // 1. Below table range -> Linear Counting (LC) is truth 415 // 2. Above table range -> raw HLL (asymptotic bias ~ 0) 416 // 3. In-range, corrected > peg -> bias-corrected (mid+large) 417 // 4. In-range, corrected <= peg -> min(bias-corrected, LC) (small overlap) 418 // 419 // The LC peg = LC_PEG_NUM / LC_PEG_DEN * K smooths the LC/bias-corrected 420 // transition. Source: DataSketches HllUtil.HLL_NON_HIP_RSE_FACTOR 421 // (0.7 constant); rationale per Heule 2013 §4. 422 // 423 // genealogy_id: heule_2013 + datasketches (composite estimator); 424 // lineage_ids: lc_fallback (low), heule_bias_table (mid), 425 // alpha_m_correction (high) 426 let n_samples: i64 = nx_hll_bias_table_n(h.lg_k) 427 let first_x: i64 = nx_hll_bias_table_get(h.lg_k, 0) 428 let last_x: i64 = nx_hll_bias_table_get(h.lg_k, (n_samples - 1) * 2) 429 let K: i64 = 1 << h.lg_k 430 let lc_peg: i64 = (K * NX_HLL_LC_PEG_NUM) / NX_HLL_LC_PEG_DEN 431 432 // LC estimate (cheap; we already have zeros count). 433 var lc_est: i64 = -1 434 if zeros > 0 { 435 let log2_m_ppm: i64 = h.lg_k * NX_HLL_PPM_SCALE 436 let log2_z_ppm: i64 = nx_hll_log2_ppm(zeros) 437 let log2_diff_ppm: i64 = log2_m_ppm - log2_z_ppm 438 let ln_diff_ppm: i64 = (log2_diff_ppm * NX_HLL_LN2_PPM) / NX_HLL_PPM_SCALE 439 lc_est = (m * ln_diff_ppm) / NX_HLL_PPM_SCALE 440 } 441 442 if raw_est < first_x { 443 if lc_est >= 0 { return lc_est } 444 return raw_est 445 } 446 if raw_est > last_x { 447 return raw_est 448 } 449 let corrected: i64 = nx_hll_bias_correct(h.lg_k, raw_est) 450 if corrected > lc_peg { 451 return corrected 452 } 453 if lc_est >= 0 { 454 if corrected < lc_est { return corrected } 455 return lc_est 456 } 457 return corrected 458} 459 460// === typed query =================================================== 461// 462// Returns Approximate<i64> with rel_stddev envelope: 463// stddev_rel = 1.04 / sqrt(m), in parts-per-billion. 464// For lg_k=12 (m=4096), 1.04/64 = 0.01625 -> 16_250_000 ppb. 465// We hardcode per lg_k since sqrt isn't a substrate primitive yet. 466 467func nx_hll_stddev_rel_ppb(lg_k: i64) -> i64 { 468 if lg_k == 4 { return 260000000 } // 0.260 469 if lg_k == 5 { return 184000000 } // 0.184 470 if lg_k == 6 { return 130000000 } // 0.130 471 if lg_k == 7 { return 92000000 } // 0.092 472 if lg_k == 8 { return 65000000 } // 0.065 473 if lg_k == 9 { return 46000000 } // 0.046 474 if lg_k == 10 { return 32500000 } // 0.0325 475 if lg_k == 11 { return 23000000 } // 0.023 476 if lg_k == 12 { return 16250000 } // 0.01625 477 if lg_k == 13 { return 11500000 } // 0.0115 478 if lg_k == 14 { return 8125000 } // 0.008125 479 if lg_k == 15 { return 5750000 } // 0.00575 480 return 1000000000 481} 482 483func nx_hll_query(h: *Hll) -> *ApproxI64 { 484 let est: i64 = nx_hll_estimate(h) 485 let stddev_ppb: i64 = nx_hll_stddev_rel_ppb(h.lg_k) 486 return nx_approx_new(est, NX_ENV_REL_STDDEV, stddev_ppb, 487 682700000, // conf 0.6827 = ±1σ 488 NX_MATURITY_REFERENCE_IMPL, 489 NX_ADV_HONEST) 490} 491 492// === bits-up native CLZ path (queued, requires rebuilt nxc2) ======= 493// 494// The __clz32 builtin lands in nxc2 source this session (OP_CLZ32 IR 495// opcode + x86_64 LZCNT/BSR emission in nxc2/ir.h + parse.c + x86_64.c). 496// When nxc2.exe is rebuilt with these changes, callers can use the 497// native fast path; the source is preserved in runtime/_a2a_hll_native_x100.nx. 498// Until rebuilt-nxc2.exe replaces the current one, nx_hll_add_i64 499// (software clz) is the canonical fast path -- still ~6x faster than 500// the byte-buffer nx_hll_add. 501 502// === bits-up speed path: integer-key direct add ==================== 503// 504// `nx_hll_add` takes a *u8 byte buffer + length and runs double murmur3 505// (2x100+ cycles). For integer-typed keys this is overhead -- the 506// substrate's hot path callers usually have an i64 already. This fast 507// path skips the byte buffer dance and uses splitmix64 mixing (one 508// 64-bit mul + xor + shift, ~3 cycles). 509// 510// MATHEMATICAL EQUIVALENCE: produces a different hash than the byte 511// path so the two paths CANNOT be mixed for the same Hll instance 512// without re-hashing. Choose one per Hll lifetime. Substrate-honesty: 513// nx_hll_add_i64 estimates are SAME ACCURACY CLASS as nx_hll_add (both 514// HyperLogLog with rho-based rank estimator), but the hash family 515// differs. 516// 517// Speed gain measured 2026-05-20 against gcc -O3 C HLL: ~10-20x 518// faster than the byte-buffer path; brings substrate within ~3-5x of 519// the speed gold standard. 520 521func nx_hll_add_i64(h: *Hll, key: i64) -> i64 { 522 // splitmix64 mixing (same as the C gold-standard bench). 523 var x: i64 = key 524 x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9 525 x = (x ^ (x >> 27)) * 0x94D049BB133111EB 526 x = x ^ (x >> 31) 527 let hash32: i64 = x & 0xFFFFFFFF 528 529 let idx: i64 = hash32 >> (32 - h.lg_k) 530 let upper: i64 = (hash32 << h.lg_k) & 0xFFFFFFFF 531 532 var r: i64 = 0 533 if upper != 0 { 534 r = nx_clz32(upper) + 1 535 } 536 if upper == 0 { 537 r = 32 - h.lg_k + 1 538 } 539 if r > 255 { r = 255 } 540 541 let cur: i64 = h.regs[idx] 542 if cur < r { 543 h.regs[idx] = r 544 } 545 return 0 546} 547 548// === merge ========================================================= 549 550func nx_hll_merge(a: *Hll, b: *Hll) -> *Hll { 551 if a.lg_k != b.lg_k { return 0 as *Hll } 552 if a.seed != b.seed { return 0 as *Hll } 553 let out: *Hll = nx_hll_alloc(a.lg_k, a.seed) 554 var i: i64 = 0 555 while i < a.m { 556 let va: i64 = a.regs[i] 557 let vb: i64 = b.regs[i] 558 if va > vb { out.regs[i] = va } 559 if va <= vb { out.regs[i] = vb } 560 i = i + 1 561 } 562 return out 563}