sketch_geomean.nx source
↩ module page · 141 lines · 4746 B
1// sketch_geomean.nx -- streaming geometric mean via log-space accumulation.
2//
3// Geometric mean = (Π x_i)^(1/n) = 2^((Σ log_2 x_i) / n)
4//
5// Naive product overflows immediately (n=20, x=10 -> product ~ 10^20 >> i64).
6// We accumulate in log_2 domain instead: log space is additive and bounded.
7//
8// USE CASES (where geometric mean beats arithmetic mean):
9// - growth-rate averaging (CAGR / compound returns)
10// - ratios (speedup factors, multiplicative gain)
11// - normalized scores (e.g., averaging across scales)
12// - magnitude data (sizes, populations, intensities)
13//
14// INTEGER LOG-BASE-2 (no f64):
15// floor(log_2(x)) via bit_length(x) - 1
16// For finer fractional bits, mantissa interpolation could be added
17// (queued for v2).
18//
19// SUM_LOG_2 fixed-point: track in PPM scale. Each insert contributes
20// (bitlen - 1) * 1_000_000. Geometric mean recovered via 2^(sum/n_ppm).
21//
22// For 2^x where x is in PPM (fractional log_2): floor(x / 1_000_000) +
23// linear interpolation on the fractional part. Returns approximate
24// geometric mean.
25//
26// LOSSLESS-LANGUAGE DISCIPLINE: NX_ENV_REL_STDDEV with param_a
27// reflecting quantization (use 50_000_000 ppb = 5% as conservative
28// bound for floor-log estimator). ReferenceImpl tier.
29
30import "syscalls.nx"
31import "sketch_types.nx"
32
33struct Geomean {
34 sum_log_ppm: i64, // Σ log_2(x_i) * 1_000_000
35 n: i64,
36 n_positive: i64, // x_i must be > 0 for log; track invalid
37}
38
39// === construction =================================================
40
41func nx_gm_alloc() -> *Geomean {
42 let raw: *u8 = sys_mmap(32)
43 let g: *Geomean = raw as *Geomean
44 g.sum_log_ppm = 0
45 g.n = 0
46 g.n_positive = 0
47 return g
48}
49
50// === log_2 floor ==================================================
51
52func nx_gm_log2_floor(x: i64) -> i64 {
53 if x <= 1 { return 0 }
54 var n: i64 = 0
55 var t: i64 = x
56 while t > 1 {
57 t = t >> 1
58 n = n + 1
59 }
60 return n
61}
62
63// === add ==========================================================
64//
65// Returns 0 on success; -1 if value <= 0 (geometric mean undefined for
66// non-positive values).
67
68func nx_gm_add(g: *Geomean, value: i64) -> i64 {
69 if value <= 0 { return -1 }
70 g.sum_log_ppm = g.sum_log_ppm + nx_gm_log2_floor(value) * 1000000
71 g.n = g.n + 1
72 g.n_positive = g.n_positive + 1
73 return 0
74}
75
76// === 2^x reconstruction ===========================================
77//
78// Input: x_ppm (fractional log_2 in PPM, x = k + f where f in [0, 1)).
79// Output: 2^(k + f) = 2^k * 2^f.
80// 2^f for f in [0, 1) approximated via linear interpolation:
81// 2^f ≈ 1 + f (Bernoulli's inequality, conservative; underestimates by ~6%)
82// More accurate: 2^f ≈ 1 + 0.69 * f (using ln(2) factor)
83// We use the linear approximation for simplicity; v2 can use a 16-entry
84// lookup table.
85
86func nx_gm_pow2_floor(x_ppm: i64) -> i64 {
87 if x_ppm < 0 { return 0 }
88 let k: i64 = x_ppm / 1000000
89 let f_ppm: i64 = x_ppm - k * 1000000 // fractional, in PPM
90 // Cap k to prevent overflow (2^63 = 9.2e18).
91 if k > 62 { return 0x7FFFFFFFFFFFFFFF }
92 var base: i64 = 1
93 if k > 0 {
94 base = 1 << k
95 }
96 // 2^f ≈ (1_000_000 + f_ppm) / 1_000_000
97 return (base * (1000000 + f_ppm)) / 1000000
98}
99
100// === geometric mean ==============================================
101//
102// gmean = 2^(sum_log_ppm / n)
103
104func nx_gm_value(g: *Geomean) -> i64 {
105 if g.n_positive == 0 { return 0 }
106 let log_mean_ppm: i64 = g.sum_log_ppm / g.n_positive
107 return nx_gm_pow2_floor(log_mean_ppm)
108}
109
110// Equivalent log-domain query (useful for callers who want to combine).
111func nx_gm_log_mean_ppm(g: *Geomean) -> i64 {
112 if g.n_positive == 0 { return 0 }
113 return g.sum_log_ppm / g.n_positive
114}
115
116// === typed envelope ==============================================
117
118func nx_gm_query(g: *Geomean) -> *ApproxI64 {
119 let v: i64 = nx_gm_value(g)
120 return nx_approx_new(v, NX_ENV_REL_STDDEV, 50000000, // ~5% conservative
121 682700000,
122 NX_MATURITY_REFERENCE_IMPL,
123 NX_ADV_HONEST)
124}
125
126func nx_gm_count(g: *Geomean) -> i64 { return g.n }
127func nx_gm_count_positive(g: *Geomean) -> i64 { return g.n_positive }
128
129func nx_gm_memory_bytes(g: *Geomean) -> i64 {
130 return 32
131}
132
133// === merge =======================================================
134
135func nx_gm_merge(a: *Geomean, b: *Geomean) -> *Geomean {
136 let out: *Geomean = nx_gm_alloc()
137 out.sum_log_ppm = a.sum_log_ppm + b.sum_log_ppm
138 out.n = a.n + b.n
139 out.n_positive = a.n_positive + b.n_positive
140 return out
141}