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