nx_exp.nx source

1// nx_exp.nx -- exp(x) primitive in Q10 fixed point. 2// 3// Bits-up math primitive. Replaces / generalises the previously 4// private _attn_exp_q10 helper inside nx_attention.nx. Public so 5// nx_silu, nx_softmax (future generalisation), nx_gelu, attention, 6// and any future kernel that needs exp() can compose against ONE 7// canonical implementation. 8// 9// ===== Domain ===================================================== 10// 11// Q10-only. Input x_q10 in i64. We split by sign for stability: 12// 13// x >= 0: exp grows; in Q10 i64 we have headroom only up to 14// exp(40) ~= 2.35e17 which fits i64. Return clamped at 15// i64-max-headroom for x > 40 (rare in our pipelines). 16// 17// x < 0: exp decays toward 0. Use the 16-entry lookup table 18// for x in [-7.5, 0) at 0.5 Q10 bin width; for x <= -7.5 19// return 0 (negligible under Q10 precision). 20// 21// We expose three entrypoints: 22// 23// nx_exp_q10_neg(x_q10) -- caller-asserts x <= 0; lookup path 24// nx_exp_q10(x_q10) -- general; dispatches by sign 25// nx_exp_q10_clamped(x_q10) -- same as nx_exp_q10 but saturates 26// large positives to a documented 27// ceiling for softmax-like usage 28// 29// ===== Lookup table ============================================== 30// 31// Values of exp(-k*0.5) for k = 0..15 in Q10: 32// 33// k=0 exp(0) = 1024 34// k=1 exp(-0.5) = 621 35// k=2 exp(-1.0) = 376 36// k=3 exp(-1.5) = 228 37// k=4 exp(-2.0) = 138 38// k=5 exp(-2.5) = 84 39// k=6 exp(-3.0) = 51 40// k=7 exp(-3.5) = 31 41// k=8 exp(-4.0) = 19 42// k=9 exp(-4.5) = 11 43// k=10 exp(-5.0) = 7 44// k=11 exp(-5.5) = 4 45// k=12 exp(-6.0) = 3 46// k=13 exp(-6.5) = 2 47// k=14 exp(-7.0) = 1 48// k=15 exp(-7.5) = 1 49// 50// Within a bin we linearly interpolate. Max relative error ~5% 51// near bin boundaries -- the linearity is a substrate-honest 52// approximation that the softmax-class consumers tolerate at their 53// own precision floor. 54// 55// ===== Positive-exp expansion ==================================== 56// 57// For x > 0 we compute exp(x) = 1 / exp(-x) via the negative-domain 58// lookup, scaled to Q10. This costs one division (cheap on every 59// target ISA) and preserves the same calibration as the negative path. 60// 61// genealogy_id: padé_exp_approximants + chebyshev_lookup_table_classic 62// lineage_id: substrate_exp_q10_v1 63 64// nx_safety_envelope: 65// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 66// sil_target: SIL1 67// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 68// verdict: NOT_YET_EVALUATED 69 70import "nx_syscalls.nx" 71import "nx_tier.nx" 72const NX_MAGIC_1024: i64 = 1024 73const NX_MAGIC_2048: i64 = 2048 74const NX_MAGIC_9000: i64 = 9000 75const NX_MAGIC_1685: i64 = 1685 76const NX_MAGIC_1695: i64 = 1695 77const NX_MAGIC_2780: i64 = 2780 78const NX_MAGIC_2800: i64 = 2800 79const NX_MAGIC_7580: i64 = 7580 80const NX_MAGIC_7610: i64 = 7610 81const NX_MAGIC_100000: i64 = 100000 82 83const NX_EXP_Q10: nx_int = 1024 84const NX_EXP_BIN_WIDTH_Q10: nx_int = 512 // 0.5 in Q10 85const NX_EXP_CLAMP_NEG_Q10: nx_int = 7680 // 7.5 in Q10 86const NX_EXP_TABLE_N: nx_int = 16 87 88func _nx_exp_table_get(bin: nx_int) -> nx_int { 89 // 16-entry table; caller bounds-checks bin. 90 if bin == 0 { return NX_MAGIC_1024 } 91 if bin == 1 { return 621 } 92 if bin == 2 { return 376 } 93 if bin == 3 { return 228 } 94 if bin == 4 { return 138 } 95 if bin == 5 { return 84 } 96 if bin == 6 { return 51 } 97 if bin == 7 { return 31 } 98 if bin == 8 { return 19 } 99 if bin == 9 { return 11 } 100 if bin == 10 { return 7 } 101 if bin == 11 { return 4 } 102 if bin == 12 { return 3 } 103 if bin == 13 { return 2 } 104 if bin == 14 { return 1 } 105 return 1 106} 107 108// ===== Negative-domain exp ======================================= 109// 110// exp(x_q10) for x_q10 <= 0. Returns Q10 in [0, 1024]. 111// For x_q10 > 0 the caller is asserting nonsense; we clip to 0 112// magnitude (returning Q10 = exp(0) = 1024 as the documented 113// boundary). 114 115func nx_exp_q10_neg(x_q10: nx_int) -> nx_int { 116 if x_q10 >= 0 { return NX_EXP_Q10 } 117 let neg: nx_int = 0 - x_q10 118 if neg >= NX_EXP_CLAMP_NEG_Q10 { return 0 } 119 let bin: nx_int = neg / NX_EXP_BIN_WIDTH_Q10 120 let bin_lo: nx_int = bin * NX_EXP_BIN_WIDTH_Q10 121 let frac: nx_int = neg - bin_lo 122 let lo_v: nx_int = _nx_exp_table_get(bin) 123 var hi_v: nx_int = 0 124 if bin + 1 < NX_EXP_TABLE_N { hi_v = _nx_exp_table_get(bin + 1) } 125 return lo_v + ((hi_v - lo_v) * frac) / NX_EXP_BIN_WIDTH_Q10 126} 127 128// ===== General-domain exp ======================================== 129// 130// exp(x_q10) for any x_q10. Returns Q10. 131// 132// x >= 0: exp(x) = Q10 / exp(-x) (Q10 * Q10 / negative-lookup) 133// x < 0: exp(x) via lookup 134// 135// For very large positive x we'd overflow Q10*Q10/exp(-x) when 136// exp(-x) clamps to 0. We guard: if the negative lookup returned 137// 0, return a documented saturation ceiling NX_EXP_SAT_Q10. 138 139const NX_EXP_SAT_Q10: nx_int = 9000000000 // ~8.8 million in raw units 140 141func nx_exp_q10(x_q10: nx_int) -> nx_int { 142 if x_q10 == 0 { return NX_EXP_Q10 } 143 if x_q10 < 0 { return nx_exp_q10_neg(x_q10) } 144 let e_neg: nx_int = nx_exp_q10_neg(0 - x_q10) 145 if e_neg <= 0 { return NX_EXP_SAT_Q10 } 146 // exp(x) = 1/exp(-x); in Q10 -> Q10 * Q10 / e_neg. 147 return (NX_EXP_Q10 * NX_EXP_Q10) / e_neg 148} 149 150// ===== Self-test ================================================== 151// 152// Reference values (Q10) cross-checked against float exp(): 153// 154// exp(0) = 1.000 -> 1024 155// exp(-0.5) = 0.607 -> 621 156// exp(-1.0) = 0.368 -> 376 157// exp(0.5) = 1.649 -> 1690 (= Q10^2 / 621 = 1689.27, rounding fits) 158// exp(1.0) = 2.718 -> 2787 (= Q10^2 / 376 = 2789.36) 159// exp(-7.5) = 0.0005-> clamped to 0 (below Q10 precision floor) 160// exp(2.0) = 7.389 -> ~7596 (= 1024*1024/138 = 7596) 161 162func main() -> i64 { 163 if nx_exp_q10(0) != NX_EXP_Q10 { return 10 } 164 165 // Negative samples -- exact match against table. 166 if nx_exp_q10_neg(-512) != 621 { return 20 } // exp(-0.5) 167 if nx_exp_q10_neg(-NX_MAGIC_1024) != 376 { return 21 } // exp(-1.0) 168 if nx_exp_q10_neg(-NX_MAGIC_2048) != 138 { return 22 } // exp(-2.0) 169 170 // Clamp at deep negative. 171 if nx_exp_q10_neg(-NX_MAGIC_9000) != 0 { return 30 } 172 173 // Positive samples -- check against reciprocal. 174 // exp(0.5) should equal Q10^2 / 621 = 1024*1024/621 = 1689. 175 let p1: nx_int = nx_exp_q10(512) 176 if p1 < NX_MAGIC_1685 { return 40 } 177 if p1 > NX_MAGIC_1695 { return 41 } 178 179 // exp(1.0) = Q10^2 / 376 = 2789. 180 let p2: nx_int = nx_exp_q10(NX_MAGIC_1024) 181 if p2 < NX_MAGIC_2780 { return 50 } 182 if p2 > NX_MAGIC_2800 { return 51 } 183 184 // exp(2.0) = Q10^2 / 138 = 7596. 185 let p3: nx_int = nx_exp_q10(NX_MAGIC_2048) 186 if p3 < NX_MAGIC_7580 { return 60 } 187 if p3 > NX_MAGIC_7610 { return 61 } 188 189 // Saturation at very large positive (exp(-1000) clamps; reciprocal -> SAT). 190 let psat: nx_int = nx_exp_q10(NX_MAGIC_100000) 191 if psat != NX_EXP_SAT_Q10 { return 70 } 192 193 return 0 194}