libm/math/expm1f.rs
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.c */
/*
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
const O_THRESHOLD: f32 = 8.8721679688e+01; /* 0x42b17180 */
const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
const INV_LN2: f32 = 1.4426950216e+00; /* 0x3fb8aa3b */
/*
* Domain [-0.34568, 0.34568], range ~[-6.694e-10, 6.696e-10]:
* |6 / x * (1 + 2 * (1 / (exp(x) - 1) - 1 / x)) - q(x)| < 2**-30.04
* Scaled coefficients: Qn_here = 2**n * Qn_for_q (see s_expm1.c):
*/
const Q1: f32 = -3.3333212137e-2; /* -0x888868.0p-28 */
const Q2: f32 = 1.5807170421e-3; /* 0xcf3010.0p-33 */
/// Exponential, base *e*, of x-1 (f32)
///
/// Calculates the exponential of `x` and subtract 1, that is, *e* raised
/// to the power `x` minus 1 (where *e* is the base of the natural
/// system of logarithms, approximately 2.71828).
/// The result is accurate even for small values of `x`,
/// where using `exp(x)-1` would lose many significant digits.
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
pub fn expm1f(mut x: f32) -> f32 {
let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
let mut hx = x.to_bits();
let sign = (hx >> 31) != 0;
hx &= 0x7fffffff;
/* filter out huge and non-finite argument */
if hx >= 0x4195b844 {
/* if |x|>=27*ln2 */
if hx > 0x7f800000 {
/* NaN */
return x;
}
if sign {
return -1.;
}
if x > O_THRESHOLD {
x *= x1p127;
return x;
}
}
let k: i32;
let hi: f32;
let lo: f32;
let mut c = 0f32;
/* argument reduction */
if hx > 0x3eb17218 {
/* if |x| > 0.5 ln2 */
if hx < 0x3F851592 {
/* and |x| < 1.5 ln2 */
if !sign {
hi = x - LN2_HI;
lo = LN2_LO;
k = 1;
} else {
hi = x + LN2_HI;
lo = -LN2_LO;
k = -1;
}
} else {
k = (INV_LN2 * x + (if sign { -0.5 } else { 0.5 })) as i32;
let t = k as f32;
hi = x - t * LN2_HI; /* t*ln2_hi is exact here */
lo = t * LN2_LO;
}
x = hi - lo;
c = (hi - x) - lo;
} else if hx < 0x33000000 {
/* when |x|<2**-25, return x */
if hx < 0x00800000 {
force_eval!(x * x);
}
return x;
} else {
k = 0;
}
/* x is now in primary range */
let hfx = 0.5 * x;
let hxs = x * hfx;
let r1 = 1. + hxs * (Q1 + hxs * Q2);
let t = 3. - r1 * hfx;
let mut e = hxs * ((r1 - t) / (6. - x * t));
if k == 0 {
/* c is 0 */
return x - (x * e - hxs);
}
e = x * (e - c) - c;
e -= hxs;
/* exp(x) ~ 2^k (x_reduced - e + 1) */
if k == -1 {
return 0.5 * (x - e) - 0.5;
}
if k == 1 {
if x < -0.25 {
return -2. * (e - (x + 0.5));
}
return 1. + 2. * (x - e);
}
let twopk = f32::from_bits(((0x7f + k) << 23) as u32); /* 2^k */
if (k < 0) || (k > 56) {
/* suffice to return exp(x)-1 */
let mut y = x - e + 1.;
if k == 128 {
y = y * 2. * x1p127;
} else {
y = y * twopk;
}
return y - 1.;
}
let uf = f32::from_bits(((0x7f - k) << 23) as u32); /* 2^-k */
if k < 23 { (x - e + (1. - uf)) * twopk } else { (x - (e + uf) + 1.) * twopk }
}