nx_logadd_lib.nx source

1// nx_logadd_lib.nx -- FIXED-POINT LOG-DOMAIN ARITHMETIC, no floats anywhere: the primitive a linear-chain CRF's 2// forward-backward pass needs (log-sum-exp) and the probability it reads back (exp of a non-positive log-domain value). 3// Scale LA_S = 1024 units per nat (Q10). la_lse(a, b) = ln(e^a + e^b) in Q10; la_exp_neg(x), x <= 0 -> a Q10 probability. 4// The two tables are built ONCE from first principles rather than typed: e^(-1/S) in Q30 from its Taylor series (the 5// third term is already below one unit at this scale, so the constant is exact to the last bit), the powers e^(-d/S) by 6// recurrence, and softplus(d) = S*ln(1 + e^(-d/S)) through the estate's ONE integer logarithm (nx_intlog ilog2_1024) -- 7// composed, never re-implemented. The horizon LA past which softplus rounds to zero is DERIVED at init by scanning the 8// table, never declared. Minus infinity is LA_NEG_INF, the same sentinel the aspect model's Viterbi already uses, so a 9// masked transition reads identically in both decoders. license_tier: ORIGINAL No hw writes (Rule 26). LIB. 10import "nx_syscalls.nx" 11import "nx_intlog.nx" 12 13const LA_S: i64 = 1024 // Q10: units per nat 14const LA_Q30: i64 = 1073741824 // 2^30, the exp table's unit 15const LA_Q30_SHIFT: i64 = 30 16const LA_Q30_HALF: i64 = 536870912 // rounding half for a Q30 product 17const LA_Q20_SHIFT: i64 = 20 18const LA_Q20_HALF: i64 = 524288 19const LA_LN2_Q20: i64 = 726817 // ln 2 in Q20 (0.6931471806 * 2^20 = 726,817.0), turns a Q10 log2 into a Q10 ln 20const LA_LOG2_Q30_Q10: i64 = 30720 // ilog2_1024(2^30) = 30 * 1024: subtracted so ln(1+p) reads from a Q30 (1+p) 21const LA_NEG_INF: i64 = 0 - 1000000000 // log-domain minus infinity (== SP_NEG in nx_absa_seq) 22const LA_NEG_INF_GUARD: i64 = 0 - 500000000 // anything at or below this is treated as minus infinity (an -inf plus a finite) 23const LA_TABLE_CAP: i64 = 16384 // table entries; the derived horizon is ~S*ln(2S), a little under half of this 24const LA_I64_BYTES: i64 = 8 25const LA_TAYLOR_TERMS: i64 = 6 // terms of e^(-1/S) at Q30: the 3rd is already < 1 unit, six is exact by any measure 26 27static la_expt: *i64 // e^(-d/S) in Q30, d = 0..LA_TABLE_CAP-1 28static la_spt: *i64 // S * ln(1 + e^(-d/S)) in Q10, d = 0..LA_TABLE_CAP-1 29static la_horizon_v: i64 // first d where la_spt[d] == 0 (derived at init) 30static la_ready: i64 31static la_clamps: i64 // la_exp_neg calls whose argument was above zero (rounding drift in a caller), announced not hidden 32 33// e^(-1/S) in Q30 by its alternating Taylor series, integer arithmetic only 34func la_exp_unit_q30() -> i64 { 35 var sum: i64 = LA_Q30 36 var term: i64 = LA_Q30 37 var k: i64 = 1 38 while k <= LA_TAYLOR_TERMS { 39 term = term / (LA_S * k) // term_k = term_(k-1) / (S*k) 40 if (k % 2) == 1 { sum = sum - term } else { sum = sum + term } 41 k = k + 1 42 } 43 return sum 44} 45// build the tables; idempotent; returns 1 46func la_init() -> i64 { 47 if la_ready == 1 { return 1 } 48 la_expt = sys_mmap(LA_TABLE_CAP * LA_I64_BYTES) as *i64 49 la_spt = sys_mmap(LA_TABLE_CAP * LA_I64_BYTES) as *i64 50 let c: i64 = la_exp_unit_q30() 51 la_expt[0] = LA_Q30 52 var d: i64 = 1 53 while d < LA_TABLE_CAP { 54 la_expt[d] = (la_expt[d - 1] * c + LA_Q30_HALF) >> LA_Q30_SHIFT 55 d = d + 1 56 } 57 // softplus(d) = S * ln(1 + e^(-d/S)) = S * ln2 * (log2(Q30 + expt[d]) - 30): the ONE integer log, then ln2 in Q20 58 d = 0 59 while d < LA_TABLE_CAP { 60 let l2: i64 = ilog2_1024(LA_Q30 + la_expt[d]) - LA_LOG2_Q30_Q10 // Q10 log2(1+p), 0..1024 61 var v: i64 = (l2 * LA_LN2_Q20 + LA_Q20_HALF) >> LA_Q20_SHIFT // Q10 ln(1+p) (S already inside the Q10) 62 if v < 0 { v = 0 } 63 la_spt[d] = v 64 d = d + 1 65 } 66 la_horizon_v = LA_TABLE_CAP 67 d = 0 68 var found: i64 = 0 69 while (d < LA_TABLE_CAP) & (found == 0) { 70 if la_spt[d] == 0 { la_horizon_v = d; found = 1 } 71 d = d + 1 72 } 73 la_clamps = 0 74 la_ready = 1 75 return 1 76} 77func la_horizon() -> i64 { return la_horizon_v } 78func la_clamp_count() -> i64 { return la_clamps } 79func la_is_neg_inf(x: i64) -> i64 { if x <= LA_NEG_INF_GUARD { return 1 } return 0 } 80// ln(e^a + e^b) in Q10 81func la_lse(a: i64, b: i64) -> i64 { 82 if la_is_neg_inf(a) == 1 { return b } 83 if la_is_neg_inf(b) == 1 { return a } 84 var m: i64 = a 85 var d: i64 = a - b 86 if b > a { m = b; d = b - a } 87 if d >= la_horizon_v { return m } 88 return m + la_spt[d] 89} 90// e^x for x <= 0, as a Q10 probability (0..LA_S); x > 0 (a caller's rounding drift) clamps to LA_S and is counted 91func la_exp_neg(x: i64) -> i64 { 92 if x > 0 { la_clamps = la_clamps + 1; return LA_S } 93 let d: i64 = 0 - x 94 if d >= LA_TABLE_CAP { return 0 } 95 return (la_expt[d] + LA_Q20_HALF) >> LA_Q20_SHIFT 96} 97// e^x for x <= 0 in Q30 (the finer reading, for callers that accumulate many small probabilities) 98func la_exp_neg_q30(x: i64) -> i64 { 99 if x > 0 { la_clamps = la_clamps + 1; return LA_Q30 } 100 let d: i64 = 0 - x 101 if d >= LA_TABLE_CAP { return 0 } 102 return la_expt[d] 103}