c/
math_exp.rs

1// Copyright (c) 2026 Meng Zhuo
2// SPDX-License-Identifier: BSD-3-Clause
3
4// Translated from musl: exp.c, exp2.c, expf.c, exp_data.c, exp2f_data.c
5
6// ============================================================
7// Data tables and constants for exp (double precision)
8// ============================================================
9
10const EXP_TABLE_BITS: u32 = 7;
11const EXP_N: u64 = 1 << EXP_TABLE_BITS;
12
13const EXP_INVLN2N: f64 = asdouble(0x40671547652B82FE);
14const EXP_NEGLN2HI_N: f64 = asdouble(0xBF762E42FEFA0000);
15const EXP_NEGLN2LO_N: f64 = asdouble(0xBD0CF79ABC9E3B3A);
16const EXP_SHIFT: f64 = asdouble(0x4338000000000000);
17
18const EXP_C2: f64 = asdouble(0x3FDFFFFFFFFFFDBD);
19const EXP_C3: f64 = asdouble(0x3FC555555555543C);
20const EXP_C4: f64 = asdouble(0x3FA55555CF172B91);
21const EXP_C5: f64 = asdouble(0x3F81111167A4D017);
22
23const EXP_1P1009: f64 = asdouble(0x7F00000000000000);
24const EXP_1P_NEG1022: f64 = asdouble(0x0010000000000000);
25
26// 2^(k/N) ~= H[k]*(1 + T[k]) for k in [0, N)
27// tab[2*k] = asuint64(T[k]), tab[2*k+1] = asuint64(H[k]) - (k << 52)/N
28const EXP_TAB: [u64; 256] = [
29    0x0, 0x3ff0000000000000,
30    0x3c9b3b4f1a88bf6e, 0x3feff63da9fb3335,
31    0xbc7160139cd8dc5d, 0x3fefec9a3e778061,
32    0xbc905e7a108766d1, 0x3fefe315e86e7f85,
33    0x3c8cd2523567f613, 0x3fefd9b0d3158574,
34    0xbc8bce8023f98efa, 0x3fefd06b29ddf6de,
35    0x3c60f74e61e6c861, 0x3fefc74518759bc8,
36    0x3c90a3e45b33d399, 0x3fefbe3ecac6f383,
37    0x3c979aa65d837b6d, 0x3fefb5586cf9890f,
38    0x3c8eb51a92fdeffc, 0x3fefac922b7247f7,
39    0x3c3ebe3d702f9cd1, 0x3fefa3ec32d3d1a2,
40    0xbc6a033489906e0b, 0x3fef9b66affed31b,
41    0xbc9556522a2fbd0e, 0x3fef9301d0125b51,
42    0xbc5080ef8c4eea55, 0x3fef8abdc06c31cc,
43    0xbc91c923b9d5f416, 0x3fef829aaea92de0,
44    0x3c80d3e3e95c55af, 0x3fef7a98c8a58e51,
45    0xbc801b15eaa59348, 0x3fef72b83c7d517b,
46    0xbc8f1ff055de323d, 0x3fef6af9388c8dea,
47    0x3c8b898c3f1353bf, 0x3fef635beb6fcb75,
48    0xbc96d99c7611eb26, 0x3fef5be084045cd4,
49    0x3c9aecf73e3a2f60, 0x3fef54873168b9aa,
50    0xbc8fe782cb86389d, 0x3fef4d5022fcd91d,
51    0x3c8a6f4144a6c38d, 0x3fef463b88628cd6,
52    0x3c807a05b0e4047d, 0x3fef3f49917ddc96,
53    0x3c968efde3a8a894, 0x3fef387a6e756238,
54    0x3c875e18f274487d, 0x3fef31ce4fb2a63f,
55    0x3c80472b981fe7f2, 0x3fef2b4565e27cdd,
56    0xbc96b87b3f71085e, 0x3fef24dfe1f56381,
57    0x3c82f7e16d09ab31, 0x3fef1e9df51fdee1,
58    0xbc3d219b1a6fbffa, 0x3fef187fd0dad990,
59    0x3c8b3782720c0ab4, 0x3fef1285a6e4030b,
60    0x3c6e149289cecb8f, 0x3fef0cafa93e2f56,
61    0x3c834d754db0abb6, 0x3fef06fe0a31b715,
62    0x3c864201e2ac744c, 0x3fef0170fc4cd831,
63    0x3c8fdd395dd3f84a, 0x3feefc08b26416ff,
64    0xbc86a3803b8e5b04, 0x3feef6c55f929ff1,
65    0xbc924aedcc4b5068, 0x3feef1a7373aa9cb,
66    0xbc9907f81b512d8e, 0x3feeecae6d05d866,
67    0xbc71d1e83e9436d2, 0x3feee7db34e59ff7,
68    0xbc991919b3ce1b15, 0x3feee32dc313a8e5,
69    0x3c859f48a72a4c6d, 0x3feedea64c123422,
70    0xbc9312607a28698a, 0x3feeda4504ac801c,
71    0xbc58a78f4817895b, 0x3feed60a21f72e2a,
72    0xbc7c2c9b67499a1b, 0x3feed1f5d950a897,
73    0x3c4363ed60c2ac11, 0x3feece086061892d,
74    0x3c9666093b0664ef, 0x3feeca41ed1d0057,
75    0x3c6ecce1daa10379, 0x3feec6a2b5c13cd0,
76    0x3c93ff8e3f0f1230, 0x3feec32af0d7d3de,
77    0x3c7690cebb7aafb0, 0x3feebfdad5362a27,
78    0x3c931dbdeb54e077, 0x3feebcb299fddd0d,
79    0xbc8f94340071a38e, 0x3feeb9b2769d2ca7,
80    0xbc87deccdc93a349, 0x3feeb6daa2cf6642,
81    0xbc78dec6bd0f385f, 0x3feeb42b569d4f82,
82    0xbc861246ec7b5cf6, 0x3feeb1a4ca5d920f,
83    0x3c93350518fdd78e, 0x3feeaf4736b527da,
84    0x3c7b98b72f8a9b05, 0x3feead12d497c7fd,
85    0x3c9063e1e21c5409, 0x3feeab07dd485429,
86    0x3c34c7855019c6ea, 0x3feea9268a5946b7,
87    0x3c9432e62b64c035, 0x3feea76f15ad2148,
88    0xbc8ce44a6199769f, 0x3feea5e1b976dc09,
89    0xbc8c33c53bef4da8, 0x3feea47eb03a5585,
90    0xbc845378892be9ae, 0x3feea34634ccc320,
91    0xbc93cedd78565858, 0x3feea23882552225,
92    0x3c5710aa807e1964, 0x3feea155d44ca973,
93    0xbc93b3efbf5e2228, 0x3feea09e667f3bcd,
94    0xbc6a12ad8734b982, 0x3feea012750bdabf,
95    0xbc6367efb86da9ee, 0x3fee9fb23c651a2f,
96    0xbc80dc3d54e08851, 0x3fee9f7df9519484,
97    0xbc781f647e5a3ecf, 0x3fee9f75e8ec5f74,
98    0xbc86ee4ac08b7db0, 0x3fee9f9a48a58174,
99    0xbc8619321e55e68a, 0x3fee9feb564267c9,
100    0x3c909ccb5e09d4d3, 0x3feea0694fde5d3f,
101    0xbc7b32dcb94da51d, 0x3feea11473eb0187,
102    0x3c94ecfd5467c06b, 0x3feea1ed0130c132,
103    0x3c65ebe1abd66c55, 0x3feea2f336cf4e62,
104    0xbc88a1c52fb3cf42, 0x3feea427543e1a12,
105    0xbc9369b6f13b3734, 0x3feea589994cce13,
106    0xbc805e843a19ff1e, 0x3feea71a4623c7ad,
107    0xbc94d450d872576e, 0x3feea8d99b4492ed,
108    0x3c90ad675b0e8a00, 0x3feeaac7d98a6699,
109    0x3c8db72fc1f0eab4, 0x3feeace5422aa0db,
110    0xbc65b6609cc5e7ff, 0x3feeaf3216b5448c,
111    0x3c7bf68359f35f44, 0x3feeb1ae99157736,
112    0xbc93091fa71e3d83, 0x3feeb45b0b91ffc6,
113    0xbc5da9b88b6c1e29, 0x3feeb737b0cdc5e5,
114    0xbc6c23f97c90b959, 0x3feeba44cbc8520f,
115    0xbc92434322f4f9aa, 0x3feebd829fde4e50,
116    0xbc85ca6cd7668e4b, 0x3feec0f170ca07ba,
117    0x3c71affc2b91ce27, 0x3feec49182a3f090,
118    0x3c6dd235e10a73bb, 0x3feec86319e32323,
119    0xbc87c50422622263, 0x3feecc667b5de565,
120    0x3c8b1c86e3e231d5, 0x3feed09bec4a2d33,
121    0xbc91bbd1d3bcbb15, 0x3feed503b23e255d,
122    0x3c90cc319cee31d2, 0x3feed99e1330b358,
123    0x3c8469846e735ab3, 0x3feede6b5579fdbf,
124    0xbc82dfcd978e9db4, 0x3feee36bbfd3f37a,
125    0x3c8c1a7792cb3387, 0x3feee89f995ad3ad,
126    0xbc907b8f4ad1d9fa, 0x3feeee07298db666,
127    0xbc55c3d956dcaeba, 0x3feef3a2b84f15fb,
128    0xbc90a40e3da6f640, 0x3feef9728de5593a,
129    0xbc68d6f438ad9334, 0x3feeff76f2fb5e47,
130    0xbc91eee26b588a35, 0x3fef05b030a1064a,
131    0x3c74ffd70a5fddcd, 0x3fef0c1e904bc1d2,
132    0xbc91bdfbfa9298ac, 0x3fef12c25bd71e09,
133    0x3c736eae30af0cb3, 0x3fef199bdd85529c,
134    0x3c8ee3325c9ffd94, 0x3fef20ab5fffd07a,
135    0x3c84e08fd10959ac, 0x3fef27f12e57d14b,
136    0x3c63cdaf384e1a67, 0x3fef2f6d9406e7b5,
137    0x3c676b2c6c921968, 0x3fef3720dcef9069,
138    0xbc808a1883ccb5d2, 0x3fef3f0b555dc3fa,
139    0xbc8fad5d3ffffa6f, 0x3fef472d4a07897c,
140    0xbc900dae3875a949, 0x3fef4f87080d89f2,
141    0x3c74a385a63d07a7, 0x3fef5818dcfba487,
142    0xbc82919e2040220f, 0x3fef60e316c98398,
143    0x3c8e5a50d5c192ac, 0x3fef69e603db3285,
144    0x3c843a59ac016b4b, 0x3fef7321f301b460,
145    0xbc82d52107b43e1f, 0x3fef7c97337b9b5f,
146    0xbc892ab93b470dc9, 0x3fef864614f5a129,
147    0x3c74b604603a88d3, 0x3fef902ee78b3ff6,
148    0x3c83c5ec519d7271, 0x3fef9a51fbc74c83,
149    0xbc8ff7128fd391f0, 0x3fefa4afa2a490da,
150    0xbc8dae98e223747d, 0x3fefaf482d8e67f1,
151    0x3c8ec3bc41aa2008, 0x3fefba1bee615a27,
152    0x3c842b94c3a9eb32, 0x3fefc52b376bba97,
153    0x3c8a64a931d185ee, 0x3fefd0765b6e4540,
154    0xbc8e37bae43be3ed, 0x3fefdbfdad9cbe14,
155    0x3c77893b4d91cd9d, 0x3fefe7c1819e90d8,
156    0x3c5305c14160cc89, 0x3feff3c22b8f71f1,
157];
158
159// ============================================================
160// Data tables and constants for expf (single precision)
161// ============================================================
162
163const EXP2F_TABLE_BITS: u32 = 5;
164const EXP2F_N: u64 = 1 << EXP2F_TABLE_BITS;
165
166const EXP2F_SHIFT_SCALED: f64 = asdouble(0x42E8000000000000);
167const EXP2F_POLY: [f64; 3] = [
168    asdouble(0x3FAC6AF84B912394),
169    asdouble(0x3FCEBFCE50FAC4F3),
170    asdouble(0x3FE62E42FF0C52D6),
171];
172const EXP2F_SHIFT: f64 = asdouble(0x4338000000000000);
173const EXP2F_INVLN2_SCALED: f64 = asdouble(0x40471547652B82FE);
174const EXP2F_POLY_SCALED: [f64; 3] = [
175    asdouble(0x3EBC6AF84B912394),
176    asdouble(0x3F2EBFCE50FAC4F3),
177    asdouble(0x3F962E42FF0C52D6),
178];
179
180// tab[i] = uint(2^(i/N)) - (i << 52-BITS)
181const EXP2F_TAB: [u64; 32] = [
182    0x3ff0000000000000, 0x3fefd9b0d3158574, 0x3fefb5586cf9890f, 0x3fef9301d0125b51,
183    0x3fef72b83c7d517b, 0x3fef54873168b9aa, 0x3fef387a6e756238, 0x3fef1e9df51fdee1,
184    0x3fef06fe0a31b715, 0x3feef1a7373aa9cb, 0x3feedea64c123422, 0x3feece086061892d,
185    0x3feebfdad5362a27, 0x3feeb42b569d4f82, 0x3feeab07dd485429, 0x3feea47eb03a5585,
186    0x3feea09e667f3bcd, 0x3fee9f75e8ec5f74, 0x3feea11473eb0187, 0x3feea589994cce13,
187    0x3feeace5422aa0db, 0x3feeb737b0cdc5e5, 0x3feec49182a3f090, 0x3feed503b23e255d,
188    0x3feee89f995ad3ad, 0x3feeff76f2fb5e47, 0x3fef199bdd85529c, 0x3fef3720dcef9069,
189    0x3fef5818dcfba487, 0x3fef7c97337b9b5f, 0x3fefa4afa2a490da, 0x3fefd0765b6e4540,
190];
191
192const EXP2F_OVFL: f32 = asfloat(0x42B17217);
193const EXP2F_UVFL: f32 = asfloat(0xC2CFF1B4);
194
195// ============================================================
196// exp (double precision)
197// ============================================================
198
199#[inline]
200fn exp_top12(x: f64) -> u32 {
201    (asuint64(x) >> 52) as u32
202}
203
204#[inline]
205fn exp_specialcase(tmp: f64, sbits: u64, ki: u64) -> f64 {
206    if (ki & 0x80000000) == 0 {
207        // k > 0, the exponent of scale might have overflowed by <= 460.
208        let sbits = sbits.wrapping_sub(1009u64 << 52);
209        let scale = asdouble(sbits);
210        let y = EXP_1P1009 * (scale + scale * tmp);
211        return eval_as_double(y);
212    }
213    // k < 0, need special care in the subnormal range.
214    let sbits = sbits + (1022u64 << 52);
215    let scale = asdouble(sbits);
216    let mut y = scale + scale * tmp;
217    if y < 1.0 {
218        let lo = scale - y + scale * tmp;
219        let hi = 1.0 + y;
220        let lo = 1.0 - hi + y + lo;
221        y = eval_as_double(hi + lo) - 1.0;
222        if y == 0.0 {
223            y = 0.0; // avoid -0.0 with downward rounding
224        }
225        // signal underflow exception explicitly
226        fp_force_eval(fp_barrier(EXP_1P_NEG1022) * EXP_1P_NEG1022);
227    }
228    y = EXP_1P_NEG1022 * y;
229    eval_as_double(y)
230}
231
232#[no_mangle]
233pub extern "C" fn exp(x: f64) -> f64 {
234    let mut abstop = exp_top12(x) & 0x7ff;
235
236    if predict_false(abstop.wrapping_sub(0x3C9) >= 0x408 - 0x3C9) {
237        if abstop.wrapping_sub(0x3C9) >= 0x80000000 {
238            // tiny x (including 0): avoid spurious underflow
239            return 1.0 + x;
240        }
241        if abstop >= 0x409 {
242            if x == f64::NEG_INFINITY {
243                return 0.0;
244            }
245            if abstop >= 0x7FF {
246                return 1.0 + x; // inf or nan
247            }
248            if asuint64(x) >> 63 != 0 {
249                return __math_uflow(0);
250            } else {
251                return __math_oflow(0);
252            }
253        }
254        // large but finite: force specialcase below
255        abstop = 0;
256    }
257
258    let z = EXP_INVLN2N * x;
259    let kd = eval_as_double(z + EXP_SHIFT);
260    let ki = asuint64(kd);
261    let kd = kd - EXP_SHIFT;
262    let r = x + kd * EXP_NEGLN2HI_N + kd * EXP_NEGLN2LO_N;
263    let idx = (2 * (ki % EXP_N)) as usize;
264    let top = ki << (52 - EXP_TABLE_BITS);
265    let tail = asdouble(EXP_TAB[idx]);
266    let sbits = EXP_TAB[idx + 1].wrapping_add(top);
267    let r2 = r * r;
268    let tmp = tail + r + r2 * (EXP_C2 + r * EXP_C3) + r2 * r2 * (EXP_C4 + r * EXP_C5);
269
270    if predict_false(abstop == 0) {
271        return exp_specialcase(tmp, sbits, ki);
272    }
273    let scale = asdouble(sbits);
274    eval_as_double(scale + scale * tmp)
275}
276
277// ============================================================
278// exp2 (double precision)
279// ============================================================
280
281const EXP2_SHIFT: f64 = asdouble(0x42C8000000000000);
282const EXP2_1P_NEG54: f64 = asdouble(0x3C90000000000000);
283const EXP2_C1: f64 = asdouble(0x3FE62E42FEFA39EF);
284const EXP2_C2: f64 = asdouble(0x3FCEBFBDFF82C424);
285const EXP2_C3: f64 = asdouble(0x3FAC6B08D70CF4B5);
286const EXP2_C4: f64 = asdouble(0x3F83B2ABD24650CC);
287const EXP2_C5: f64 = asdouble(0x3F55D7E09B4E3A84);
288
289#[inline]
290fn exp2_specialcase(tmp: f64, sbits: u64, ki: u64) -> f64 {
291    if (ki & 0x80000000) == 0 {
292        // k > 0, the exponent of scale might have overflowed by 1.
293        let sbits = sbits.wrapping_sub(1u64 << 52);
294        let scale = asdouble(sbits);
295        let y = 2.0 * (scale + scale * tmp);
296        return eval_as_double(y);
297    }
298    // k < 0, need special care in the subnormal range.
299    let sbits = sbits.wrapping_add(1022u64 << 52);
300    let scale = asdouble(sbits);
301    let mut y = scale + scale * tmp;
302    if y < 1.0 {
303        let lo = scale - y + scale * tmp;
304        let hi = 1.0 + y;
305        let lo = 1.0 - hi + y + lo;
306        y = eval_as_double(hi + lo) - 1.0;
307        // Avoid -0.0 with downward rounding.
308        if y == 0.0 {
309            y = 0.0;
310        }
311        // The underflow exception needs to be signaled explicitly.
312        fp_force_eval(fp_barrier(EXP_1P_NEG1022) * EXP_1P_NEG1022);
313    }
314    y = EXP_1P_NEG1022 * y;
315    eval_as_double(y)
316}
317
318pub(crate) fn exp2_impl(x: f64) -> f64 {
319    let top12 = |v: f64| (asuint64(v) >> 52) as u32;
320    let mut abstop = top12(x) & 0x7ff;
321    if predict_false(abstop.wrapping_sub(top12(EXP2_1P_NEG54)) >= top12(512.0) - top12(EXP2_1P_NEG54)) {
322        if abstop.wrapping_sub(top12(EXP2_1P_NEG54)) >= 0x80000000 {
323            // Avoid spurious underflow for tiny x. Note: 0 is common input.
324            return 1.0 + x;
325        }
326        if abstop >= top12(1024.0) {
327            if asuint64(x) == asuint64(f64::NEG_INFINITY) {
328                return 0.0;
329            }
330            if abstop >= top12(f64::INFINITY) {
331                return 1.0 + x;
332            }
333            if asuint64(x) >> 63 == 0 {
334                return __math_oflow(0);
335            } else if asuint64(x) >= asuint64(-1075.0) {
336                return __math_uflow(0);
337            }
338        }
339        if 2u64.wrapping_mul(asuint64(x)) > 2u64.wrapping_mul(asuint64(928.0)) {
340            // Large x is special cased below.
341            abstop = 0;
342        }
343    }
344
345    // exp2(x) = 2^(k/N) * 2^r, with 2^r in [2^(-1/2N),2^(1/2N)].
346    let kd = eval_as_double(x + EXP2_SHIFT);
347    let ki = asuint64(kd);
348    let kd = kd - EXP2_SHIFT;
349    let r = x - kd;
350    let idx = (2 * (ki % EXP_N)) as usize;
351    let top = ki << (52 - EXP_TABLE_BITS);
352    let tail = asdouble(EXP_TAB[idx]);
353    let sbits = EXP_TAB[idx + 1].wrapping_add(top);
354    let r2 = r * r;
355    let tmp = tail
356        + r * EXP2_C1
357        + r2 * (EXP2_C2 + r * EXP2_C3)
358        + r2 * r2 * (EXP2_C4 + r * EXP2_C5);
359    if predict_false(abstop == 0) {
360        return exp2_specialcase(tmp, sbits, ki);
361    }
362    let scale = asdouble(sbits);
363    eval_as_double(scale + scale * tmp)
364}
365
366// ============================================================
367// expf (single precision)
368// ============================================================
369
370#[inline]
371fn expf_top12(x: f32) -> u32 {
372    asuint(x) >> 20
373}
374
375#[no_mangle]
376pub extern "C" fn expf(x: f32) -> f32 {
377    let xd = x as f64;
378    let abstop = expf_top12(x) & 0x7ff;
379
380    if predict_false(abstop >= 0x42B) {
381        if x == f32::NEG_INFINITY {
382            return 0.0;
383        }
384        if abstop >= 0x7F8 {
385            return x + x; // inf or nan
386        }
387        if x > EXP2F_OVFL {
388            return __math_oflowf(0);
389        }
390        if x < EXP2F_UVFL {
391            return __math_uflowf(0);
392        }
393    }
394
395    let z = EXP2F_INVLN2_SCALED * xd;
396    let kd = eval_as_double(z + EXP2F_SHIFT);
397    let ki = asuint64(kd);
398    let kd = kd - EXP2F_SHIFT;
399    let r = z - kd;
400
401    let t = EXP2F_TAB[(ki % EXP2F_N) as usize].wrapping_add(ki << (52 - EXP2F_TABLE_BITS));
402    let s = asdouble(t);
403    let z = EXP2F_POLY_SCALED[0] * r + EXP2F_POLY_SCALED[1];
404    let r2 = r * r;
405    let mut y = EXP2F_POLY_SCALED[2] * r + 1.0;
406    y = z * r2 + y;
407    y = y * s;
408    y as f32
409}
c/math_exp.rs

c/
math_exp.rs
