f256/math/
exp.rs

1// ---------------------------------------------------------------------------
2// Copyright:   (c) 2025 ff. Michael Amrhein (michael@adrhinum.de)
3// License:     This program is part of a larger application. For license
4//              details please read the file LICENSE.TXT provided together
5//              with the application.
6// ---------------------------------------------------------------------------
7// $Source$
8// $Revision$
9
10use core::num::FpCategory;
11
12use super::{
13    bkm::{bkm_e, bkm_l},
14    pow::approx_powf,
15    Float256, Float512,
16};
17use crate::math::log::approx_ln;
18use crate::{
19    abs_bits,
20    big_uint::{BigUInt, U256},
21    consts, f256, norm_signif_exp, SIGNIFICAND_BITS,
22};
23
24// f256::MAX.ln()
25// 181704.374500706303191870897247532238261583907221734753336211540408636177
26const LN_MAX: f256 = f256 {
27    bits: U256::new(
28        85075909386780507726367091749182177406,
29        126058098762831240896900834817458532725,
30    ),
31};
32
33// f256::MAX.log2() - 2⁻²³⁶
34// 262143.999999999999999999999999999999999999999999999999999999999999999999
35const LOG2_MAX: f256 = f256 {
36    bits: U256::new(
37        0x40010fffffffffffffffffffffffffff,
38        0xffffffffffffffffffffffffffffffff,
39    ),
40};
41
42#[allow(clippy::cast_sign_loss)]
43pub(crate) fn approx_exp(x: &Float512) -> Float512 {
44    // 2⁻²³⁶ <= |x| <= LN_MAX
45    let mut m = x.abs();
46    let mut e = 0_i32;
47    // exp(|x|) = exp(m⋅2ᵉ)
48    // |x| <= LN_MAX => |e| <= 18
49    // Assure that m < 1.5 as pre-condition for using fn bkm_e.
50    if m > Float512::THREE_HALF {
51        let ms = m.signif();
52        if ms < Float512::THREE_HALF.signif() {
53            e = m.exp();
54            m = Float512::from(&ms);
55        } else {
56            e = m.exp() + 1;
57            m = Float512::from(&ms).mul_pow2(-1);
58        }
59    }
60    let mut res = match e {
61        0 => bkm_e(&m),
62        1.. => {
63            // e >= 1 => exp(|x|) = exp(m)ⁿ with n =⋅2ᵉ
64            let n = 1_i32 << e as u32;
65            bkm_e(&m).powi(n)
66        }
67        -1 => {
68            // e = -1 => exp(|x|) = exp(m)ʸ with y = ½ = √(exp(m))
69            bkm_e(&m).sqrt()
70        }
71        _ => {
72            // e < -1
73            let mut a = Float512::TWO.mul_pow2(e);
74            let mut t = bkm_e(&m);
75            let w = a * approx_ln(&t);
76            t * approx_exp(&w)
77        }
78    };
79    // x < 0 => eˣ = 1/e⁻ˣ
80    if x.signum() == -1 {
81        res = res.recip();
82    }
83    res
84}
85
86impl f256 {
87    /// Returns e^(self), (the exponential function).
88    #[must_use]
89    pub fn exp(&self) -> Self {
90        // x = 0  or x is subnornal => eˣ = 1
91        // x = ∞ => eˣ = ∞
92        // x = -∞ => eˣ = 0
93        // x is nan => eˣ is nan
94        match self.classify() {
95            FpCategory::Zero | FpCategory::Subnormal => Self::ONE,
96            FpCategory::Infinite => {
97                [Self::INFINITY, Self::ZERO][self.sign() as usize]
98            }
99            FpCategory::Nan => Self::NAN,
100            _ => {
101                // self is finite and != 0
102                if self == &Self::ONE {
103                    // x = 1 => eˣ = e
104                    return consts::E;
105                }
106                let self_abs = self.abs();
107                if self_abs <= Self::EPSILON {
108                    // for very small x, eˣ ≅ 1+x+½x²
109                    let x = Float512::from(self);
110                    return Self::from(
111                        &(Float512::ONE + x + x.square().mul_pow2(-1)),
112                    );
113                }
114                if self_abs > LN_MAX {
115                    return [Self::INFINITY, Self::ZERO]
116                        [self.sign() as usize];
117                }
118                Self::from(&approx_exp(&Float512::from(self)))
119            }
120        }
121    }
122
123    /// Returns e^(self) - 1 in a way that is accurate even if the number is
124    /// close to zero.
125    #[must_use]
126    pub fn exp_m1(&self) -> Self {
127        // x = 0  or x is subnornal => eˣ-1 = 0
128        // x = ∞ => eˣ-1 = ∞
129        // x = -∞ => eˣ-1 = -1
130        // x is nan => eˣ-1 is nan
131        match self.classify() {
132            FpCategory::Zero | FpCategory::Subnormal => Self::ZERO,
133            FpCategory::Infinite => {
134                [Self::INFINITY, Self::NEG_ONE][self.sign() as usize]
135            }
136            FpCategory::Nan => Self::NAN,
137            _ => {
138                // self is finite and != 0
139                if self == &Self::ONE {
140                    // x = 1 => eˣ-1 = e-1
141                    return consts::E - Self::ONE;
142                }
143                let self_abs = self.abs();
144                if self_abs <= Self::EPSILON {
145                    // for very small x, eˣ-1 ≅ x+½x²
146                    let x = Float512::from(self);
147                    return Self::from(&(x + x.square().mul_pow2(-1)));
148                }
149                if self_abs > LN_MAX {
150                    return [Self::INFINITY, Self::NEG_ONE]
151                        [self.sign() as usize];
152                }
153                Self::from(
154                    &(approx_exp(&Float512::from(self)) - Float512::ONE),
155                )
156            }
157        }
158    }
159
160    /// Returns 2^(self).
161    #[must_use]
162    pub fn exp2(&self) -> Self {
163        const LOG2_MIN: f256 = f256::power_of_two(-510);
164        // x = 0  or x is subnornal => 2ˣ = 1
165        // x = ∞ => 2ˣ = ∞
166        // x = -∞ => 2ˣ = 0
167        // x is nan => 2ˣ is nan
168        match self.classify() {
169            FpCategory::Zero | FpCategory::Subnormal => Self::ONE,
170            FpCategory::Infinite => {
171                [Self::INFINITY, Self::ZERO][self.sign() as usize]
172            }
173            FpCategory::Nan => Self::NAN,
174            _ => {
175                // self is finite and != 0
176                if let Ok(e) = i32::try_from(self) {
177                    return Self::power_of_two(e);
178                }
179                let self_abs = self.abs();
180                if self_abs < LOG2_MIN {
181                    return Self::ONE;
182                }
183                if self.abs() > LOG2_MAX {
184                    return [Self::INFINITY, Self::ZERO]
185                        [self.sign() as usize];
186                }
187                // 2ˣ = eʷ with w = x⋅logₑ 2
188                Self::from(&approx_exp(
189                    &(Float512::from(self) * Float512::LN_2),
190                ))
191            }
192        }
193    }
194}
195
196#[cfg(test)]
197mod exp_tests {
198    use super::*;
199    use crate::big_uint::HiLo;
200    use crate::consts::E;
201    use core::ops::Neg;
202
203    #[test]
204    fn calc_ln_max() {
205        let ln_max = f256::MAX.ln();
206        assert_eq!(ln_max, LN_MAX);
207        assert!(ln_max.exp().diff_within_n_bits(&f256::MAX, 17));
208    }
209
210    #[test]
211    fn test_specials() {
212        assert!(f256::NAN.exp().is_nan());
213        assert_eq!(f256::INFINITY.exp(), f256::INFINITY);
214        assert_eq!(f256::NEG_INFINITY.exp(), f256::ZERO);
215        assert_eq!(f256::ZERO.exp(), f256::ONE);
216        assert_eq!(f256::NEG_ZERO.exp(), f256::ONE);
217    }
218
219    #[test]
220    fn test_subnormal() {
221        assert_eq!(f256::MIN_GT_ZERO.exp(), f256::ONE);
222        let mut f = f256::MIN_POSITIVE;
223        f = f - f.ulp();
224        assert!(f.is_subnormal());
225        assert_eq!(f.exp(), f256::ONE);
226    }
227
228    #[test]
229    fn test_near_zero() {
230        assert_eq!(f256::MIN_POSITIVE.exp(), f256::ONE);
231        let mut f = f256::EPSILON.div_pow2(2);
232        assert_eq!(f.exp(), f256::ONE);
233        f += f.ulp();
234        assert_eq!(f.exp(), f256::ONE);
235    }
236
237    #[test]
238    fn test_near_one() {
239        assert_eq!(f256::ONE.exp(), E);
240        let mut f = f256::ONE + f256::EPSILON;
241        assert_eq!(f.exp(), E + E.ulp());
242        let mut f = f256::ONE - f256::EPSILON / f256::TWO;
243        assert_eq!(f.exp(), E - E.ulp());
244    }
245
246    #[test]
247    fn test_near_epsilon() {
248        let f = f256::EPSILON;
249        assert_eq!(f.exp(), f256::ONE + f);
250        let g = f - f.ulp().div2();
251        assert_eq!(g.exp(), f256::ONE + f);
252        let h = f.div2() - f.ulp().div2();
253        assert_eq!(h.exp(), f256::ONE);
254    }
255
256    #[test]
257    fn test_overflow() {
258        let f = LN_MAX + LN_MAX.ulp();
259        assert_eq!(f.exp(), f256::INFINITY);
260        assert_eq!(f.neg().exp(), f256::ZERO);
261    }
262}
263
264#[cfg(test)]
265mod exp_m1_tests {
266    use super::*;
267    // use crate::big_uint::HiLo;
268    use crate::consts::E;
269    use core::ops::Neg;
270
271    #[test]
272    fn test_specials() {
273        assert!(f256::NAN.exp_m1().is_nan());
274        assert_eq!(f256::INFINITY.exp_m1(), f256::INFINITY);
275        assert_eq!(f256::NEG_INFINITY.exp_m1(), f256::NEG_ONE);
276        assert_eq!(f256::ZERO.exp_m1(), f256::ZERO);
277        assert_eq!(f256::NEG_ZERO.exp_m1(), f256::ZERO);
278    }
279
280    #[test]
281    fn test_subnormal() {
282        assert_eq!(f256::MIN_GT_ZERO.exp_m1(), f256::ZERO);
283        let mut f = f256::MIN_POSITIVE;
284        f = f - f.ulp();
285        assert!(f.is_subnormal());
286        assert_eq!(f.exp_m1(), f256::ZERO);
287    }
288
289    #[test]
290    fn test_near_zero() {
291        assert_eq!(f256::MIN_POSITIVE.exp_m1(), f256::MIN_POSITIVE);
292        let mut f = f256::EPSILON.div_pow2(2);
293        assert_eq!(f.exp_m1(), f);
294        f += f.ulp();
295        assert_eq!(f.exp_m1(), f);
296    }
297
298    #[test]
299    fn test_near_one() {
300        assert_eq!(f256::ONE.exp_m1(), E - f256::ONE);
301        let mut f = f256::ONE + f256::EPSILON;
302        assert_eq!(f.exp_m1(), E + E.ulp() - f256::ONE);
303        let mut f = f256::ONE - f256::EPSILON.div2();
304        assert_eq!(f.exp_m1(), E - E.ulp() - f256::ONE);
305    }
306
307    #[test]
308    fn test_near_epsilon() {
309        let f = f256::EPSILON;
310        assert_eq!(f.exp_m1(), f);
311        let g = f - f.ulp().div2();
312        assert_eq!(g.exp_m1(), f);
313        let h = f.div2() - f.ulp().div2();
314        assert_eq!(h.exp_m1(), h);
315    }
316
317    #[test]
318    fn test_overflow() {
319        let f = LN_MAX + LN_MAX.ulp();
320        assert_eq!(f.exp_m1(), f256::INFINITY);
321        assert_eq!(f.neg().exp_m1(), f256::NEG_ONE);
322    }
323}
324
325#[cfg(test)]
326mod exp2_tests {
327    use super::*;
328    use crate::big_uint::HiLo;
329    use crate::consts::E;
330    use core::ops::Neg;
331
332    #[test]
333    fn calc_log2_max() {
334        let mut log2_max = f256::MAX.log2();
335        log2_max -= log2_max.ulp().div2();
336        assert_eq!(log2_max, LOG2_MAX);
337        assert!(log2_max.exp2().diff_within_n_bits(&f256::MAX, 18));
338    }
339
340    #[test]
341    fn test_specials() {
342        assert!(f256::NAN.exp2().is_nan());
343        assert_eq!(f256::INFINITY.exp2(), f256::INFINITY);
344        assert_eq!(f256::NEG_INFINITY.exp2(), f256::ZERO);
345        assert_eq!(f256::ZERO.exp2(), f256::ONE);
346        assert_eq!(f256::NEG_ZERO.exp2(), f256::ONE);
347    }
348
349    #[test]
350    fn test_subnormal() {
351        assert_eq!(f256::MIN_GT_ZERO.exp2(), f256::ONE);
352        let mut f = f256::MIN_POSITIVE;
353        f = f - f.ulp();
354        assert!(f.is_subnormal());
355        assert_eq!(f.exp2(), f256::ONE);
356    }
357
358    #[test]
359    fn test_near_zero() {
360        assert_eq!(f256::MIN_POSITIVE.exp2(), f256::ONE);
361        let mut f = f256::EPSILON.div_pow2(2);
362        assert_eq!(f.exp2(), f256::ONE);
363        f += f.ulp();
364        assert_eq!(f.exp2(), f256::ONE);
365    }
366
367    #[test]
368    fn test_near_one() {
369        assert_eq!(f256::ONE.exp2(), f256::TWO);
370        let mut f = f256::ONE + f256::EPSILON;
371        assert_eq!(f.exp2(), f256::TWO + f256::TWO.ulp());
372        let mut f = f256::ONE - f256::EPSILON.div2();
373        assert_eq!(f.exp2(), f256::TWO - f256::TWO.ulp().div2());
374    }
375
376    #[test]
377    fn test_near_epsilon() {
378        let f = f256::EPSILON;
379        assert_eq!(f.exp2(), f256::ONE + f);
380        let g = f - f.ulp().div2();
381        assert_eq!(g.exp2(), f256::ONE + f);
382        let h = f.div2() - f.ulp().div2();
383        assert_eq!(h.exp2(), f256::ONE);
384    }
385
386    #[test]
387    fn test_min() {
388        let f = -(LOG2_MAX + LOG2_MAX.ulp());
389        assert_eq!(f.exp2(), f256::MIN_POSITIVE.div_pow2(2));
390    }
391
392    #[test]
393    fn test_overflow() {
394        let f = LOG2_MAX + LOG2_MAX.ulp();
395        assert_eq!(f.exp2(), f256::INFINITY);
396    }
397}
f256/math/exp.rs

f256/math/
exp.rs
