libm/math/
expm1.rs

1/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.c */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13const O_THRESHOLD: f64 = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
14const LN2_HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
15const LN2_LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
16const INVLN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
17/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
18const Q1: f64 = -3.33333333333331316428e-02; /* BFA11111 111110F4 */
19const Q2: f64 = 1.58730158725481460165e-03; /* 3F5A01A0 19FE5585 */
20const Q3: f64 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
21const Q4: f64 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
22const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
23
24/// Exponential, base *e*, of x-1 (f64)
25///
26/// Calculates the exponential of `x` and subtract 1, that is, *e* raised
27/// to the power `x` minus 1 (where *e* is the base of the natural
28/// system of logarithms, approximately 2.71828).
29/// The result is accurate even for small values of `x`,
30/// where using `exp(x)-1` would lose many significant digits.
31#[cfg_attr(assert_no_panic, no_panic::no_panic)]
32pub fn expm1(mut x: f64) -> f64 {
33    let hi: f64;
34    let lo: f64;
35    let k: i32;
36    let c: f64;
37    let mut t: f64;
38    let mut y: f64;
39
40    let mut ui = x.to_bits();
41    let hx = ((ui >> 32) & 0x7fffffff) as u32;
42    let sign = (ui >> 63) as i32;
43
44    /* filter out huge and non-finite argument */
45    if hx >= 0x4043687A {
46        /* if |x|>=56*ln2 */
47        if x.is_nan() {
48            return x;
49        }
50        if sign != 0 {
51            return -1.0;
52        }
53        if x > O_THRESHOLD {
54            x *= f64::from_bits(0x7fe0000000000000);
55            return x;
56        }
57    }
58
59    /* argument reduction */
60    if hx > 0x3fd62e42 {
61        /* if  |x| > 0.5 ln2 */
62        if hx < 0x3FF0A2B2 {
63            /* and |x| < 1.5 ln2 */
64            if sign == 0 {
65                hi = x - LN2_HI;
66                lo = LN2_LO;
67                k = 1;
68            } else {
69                hi = x + LN2_HI;
70                lo = -LN2_LO;
71                k = -1;
72            }
73        } else {
74            k = (INVLN2 * x + if sign != 0 { -0.5 } else { 0.5 }) as i32;
75            t = k as f64;
76            hi = x - t * LN2_HI; /* t*ln2_hi is exact here */
77            lo = t * LN2_LO;
78        }
79        x = hi - lo;
80        c = (hi - x) - lo;
81    } else if hx < 0x3c900000 {
82        /* |x| < 2**-54, return x */
83        if hx < 0x00100000 {
84            unsafe { ::core::ptr::read_volatile(&x) };force_eval!(x);
85        }
86        return x;
87    } else {
88        c = 0.0;
89        k = 0;
90    }
91
92    /* x is now in primary range */
93    let hfx = 0.5 * x;
94    let hxs = x * hfx;
95    let r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
96    t = 3.0 - r1 * hfx;
97    let mut e = hxs * ((r1 - t) / (6.0 - x * t));
98    if k == 0 {
99        /* c is 0 */
100        return x - (x * e - hxs);
101    }
102    e = x * (e - c) - c;
103    e -= hxs;
104    /* exp(x) ~ 2^k (x_reduced - e + 1) */
105    if k == -1 {
106        return 0.5 * (x - e) - 0.5;
107    }
108    if k == 1 {
109        if x < -0.25 {
110            return -2.0 * (e - (x + 0.5));
111        }
112        return 1.0 + 2.0 * (x - e);
113    }
114    ui = ((0x3ff + k) as u64) << 52; /* 2^k */
115    let twopk = f64::from_bits(ui);
116    if !(0..=56).contains(&k) {
117        /* suffice to return exp(x)-1 */
118        y = x - e + 1.0;
119        if k == 1024 {
120            y = y * 2.0 * f64::from_bits(0x7fe0000000000000);
121        } else {
122            y = y * twopk;
123        }
124        return y - 1.0;
125    }
126    ui = ((0x3ff - k) as u64) << 52; /* 2^-k */
127    let uf = f64::from_bits(ui);
128    if k < 20 {
129        y = (x - e + (1.0 - uf)) * twopk;
130    } else {
131        y = (x - (e + uf) + 1.0) * twopk;
132    }
133    y
134}
135
136#[cfg(test)]
137mod tests {
138    #[test]
139    fn sanity_check() {
140        assert_eq!(super::expm1(1.1), 2.0041660239464334);
141    }
142}
libm/math/expm1.rs

libm/math/
expm1.rs
