nx_f32_exp.nx source
↩ module page · 177 lines · 7048 B
1// nx_f32_exp.nx -- IEEE 754 binary32 exp(x) bits-up.
2//
3// L6 of the bits-up numeric tower. Composes L4 mul/add/sub/cvt to
4// produce a real-valued exp(x) without libm.
5//
6// Algorithm (standard, no code borrowed):
7//
8// 1. Special cases per IEEE 754:
9// exp(NaN) = NaN
10// exp(+0) = +1
11// exp(-0) = +1
12// exp(+inf) = +inf
13// exp(-inf) = +0
14// Overflow (x >> 88) -> +inf
15// Underflow (x << -87) -> +0
16//
17// 2. Range reduction: x = k * ln(2) + r, where k is an integer
18// and r in [-ln(2)/2, ln(2)/2] ~= [-0.347, 0.347].
19// k = round(x / ln(2))
20// r = x - k * ln(2)
21// Then exp(x) = 2^k * exp(r).
22//
23// 3. Polynomial approximation of exp(r) on reduced range:
24// Taylor through r^6 (Horner form, 6 multiplies + 6 adds).
25// Worst-case ULP error on r in [-0.347, 0.347]: ~4-8 ULPs in f32
26// (deemed acceptable for ML softmax; tighter Remez minimax is
27// a v2 lift).
28//
29// 4. Multiply by 2^k: add k to f32 biased exponent. Overflow ->
30// inf; underflow to subnormal range -> conservative 0 (v1).
31//
32// References absorbed clean-room (no code borrowed):
33// Hart 1968, "Computer Approximations" (range reduction
34// + polynomial design)
35// Cody+Waite 1980, "Software Manual for Elementary Functions"
36// Muller 2016, "Elementary Functions: Algorithms and Implementation"
37//
38// genealogy_id: range_reduction_canonical + taylor_horner_form
39// lineage_id: substrate_f32_exp_v1_taylor6
40
41import "nx_syscalls.nx"
42import "nx_tier.nx"
43import "nx_f32.nx"
44import "nx_f32_div.nx"
45import "nx_f32_cvt.nx"
46
47// f32 bit patterns of the constants we use (clean-room from spec).
48const NX_F32_ONE: i64 = 0x3F800000 // 1.0
49const NX_F32_LN2: i64 = 0x3F317218 // ln(2) = 0.6931472
50const NX_F32_INV_LN2: i64 = 0x3FB8AA3B // 1 / ln(2) = 1.4426950
51const NX_F32_INV_2: i64 = 0x3F000000 // 1/2 = 0.5
52const NX_F32_INV_3: i64 = 0x3EAAAAAB // 1/3 = 0.33333334
53const NX_F32_INV_4: i64 = 0x3E800000 // 1/4 = 0.25
54const NX_F32_INV_5: i64 = 0x3E4CCCCD // 1/5 = 0.2
55const NX_F32_INV_6: i64 = 0x3E2AAAAB // 1/6 = 0.16666667
56const NX_F32_OVF_THRESH: i64 = 0x42B0C0A5 // 88.0 approx
57const NX_F32_UNF_THRESH: i64 = 0xC2AEAC50 // -87.336 approx
58
59// ===== f32 -> i32 round-to-nearest-even (private helper) ==========
60//
61// Used for range reduction (k = round(x / ln2)).
62// Returns truncated i64 on overflow (no inf path needed for our
63// bounded usage 88 / 0.693 < 128).
64
65func _f32_to_i32_rne(value: i64) -> i64 {
66 let cls: nx_int = nx_f32_classify(value)
67 if cls == NX_F32_CLS_ZERO { return 0 }
68 if cls == NX_F32_CLS_NAN { return 0 }
69 if cls == NX_F32_CLS_INF { return 0 } // bounded usage
70
71 let sign: i64 = nx_f32_sign(value)
72 let exp_field: i64 = nx_f32_exp_field(value)
73 let mant: i64 = nx_f32_mant_field(value)
74 let sig: i64 = mant | NX_F32_IMPLICIT_1
75 let real_e: i64 = exp_field - NX_F32_EXP_BIAS
76
77 var result: i64 = 0
78 if real_e < 0 - 1 {
79 return 0 // |value| < 0.5 -> rounds to 0
80 }
81 if real_e >= 23 {
82 let shift_up: i64 = real_e - 23
83 if shift_up > 30 { return 0 }
84 result = sig << shift_up
85 } else {
86 let shift_down: i64 = 23 - real_e
87 let lost_mask: i64 = (1 << shift_down) - 1
88 let lost: i64 = sig & lost_mask
89 result = sig >> shift_down
90 let halfway: i64 = 1 << (shift_down - 1)
91 if lost > halfway { result = result + 1 }
92 if lost == halfway {
93 if (result & 1) == 1 { result = result + 1 }
94 }
95 }
96 if sign == 1 { return 0 - result }
97 return result
98}
99
100// ===== Multiply f32 value by 2^k (private helper) =================
101//
102// Equivalent to ldexp. Adjusts biased exponent. Overflow -> inf,
103// underflow -> conservative zero (v1).
104
105func _f32_ldexp(value: i64, k: i64) -> i64 {
106 if value == 0 { return 0 }
107 let cls: nx_int = nx_f32_classify(value)
108 if cls == NX_F32_CLS_NAN { return NX_F32_NAN_RAW }
109 if cls == NX_F32_CLS_INF { return value }
110 if cls == NX_F32_CLS_SUBNORMAL { return value } // v1 conservative
111
112 let exp_field: i64 = nx_f32_exp_field(value)
113 let new_exp: i64 = exp_field + k
114
115 if new_exp >= 255 {
116 let sign: i64 = nx_f32_sign(value)
117 return (sign << 31) | 0x7F800000
118 }
119 if new_exp <= 0 { return 0 }
120
121 let cleared: i64 = value & 0x807FFFFF // clear exp field
122 return cleared | (new_exp << 23)
123}
124
125// ===== exp(x) =====================================================
126
127func nx_f32_exp(x: i64) -> i64 {
128 let cls: nx_int = nx_f32_classify(x)
129
130 if cls == NX_F32_CLS_NAN { return NX_F32_NAN_RAW }
131 if cls == NX_F32_CLS_ZERO { return NX_F32_ONE }
132 if cls == NX_F32_CLS_INF {
133 if nx_f32_sign(x) == 1 { return 0 } // exp(-inf) = +0
134 return 0x7F800000 // exp(+inf) = +inf
135 }
136
137 // Overflow / underflow short-circuits (avoid polynomial garbage).
138 let sign_x: i64 = nx_f32_sign(x)
139 if sign_x == 0 {
140 // Positive: compare with overflow threshold (88).
141 let abs_x: i64 = x & 0x7FFFFFFF
142 if abs_x > NX_F32_OVF_THRESH { return 0x7F800000 }
143 } else {
144 // Negative: compare magnitude with underflow threshold (87.336).
145 let abs_x: i64 = x & 0x7FFFFFFF
146 let abs_uthresh: i64 = NX_F32_UNF_THRESH & 0x7FFFFFFF
147 if abs_x > abs_uthresh { return 0 }
148 }
149
150 // Range reduction: x = k*ln(2) + r. HARDWARE SSE f32 arithmetic
151 // (__f32_mul/__f32_add, 2026-07-10): IEEE-754 round-to-nearest, so
152 // BIT-IDENTICAL to the software nx_f32_mul/add here (r is small ~<0.35
153 // -> intermediates ~1.0, no subnormals/FTZ edge), but ~1 instr vs ~30.
154 // exp is the hottest transcendental in the LLM (silu ACT + softmax
155 // ATTN, ~100ms/token). nx_f32_exp_test gates the bit-exactness.
156 // (no __f32_sub intrinsic -> subtract via sign-flipped __f32_add.)
157 let y: i64 = __f32_mul(x, NX_F32_INV_LN2) // y = x / ln(2)
158 let k: i64 = _f32_to_i32_rne(y) // k = round(y)
159 let kf: i64 = nx_i32_to_f32(k)
160 let r: i64 = __f32_add(x, __f32_mul(kf, NX_F32_LN2) ^ 0x80000000)
161
162 // Polynomial: exp(r) ~= 1 + r(1 + r/2(1 + r/3(1 + r/4(1 + r/5(1 + r/6)))))
163 let r_inv6: i64 = __f32_mul(r, NX_F32_INV_6)
164 let h6: i64 = __f32_add(NX_F32_ONE, r_inv6)
165 let r_inv5: i64 = __f32_mul(r, NX_F32_INV_5)
166 let h5: i64 = __f32_add(NX_F32_ONE, __f32_mul(r_inv5, h6))
167 let r_inv4: i64 = __f32_mul(r, NX_F32_INV_4)
168 let h4: i64 = __f32_add(NX_F32_ONE, __f32_mul(r_inv4, h5))
169 let r_inv3: i64 = __f32_mul(r, NX_F32_INV_3)
170 let h3: i64 = __f32_add(NX_F32_ONE, __f32_mul(r_inv3, h4))
171 let r_inv2: i64 = __f32_mul(r, NX_F32_INV_2)
172 let h2: i64 = __f32_add(NX_F32_ONE, __f32_mul(r_inv2, h3))
173 let exp_r: i64 = __f32_add(NX_F32_ONE, __f32_mul(r, h2))
174
175 // Multiply by 2^k.
176 return _f32_ldexp(exp_r, k)
177}