Skip to content

Commit 4e7504a

Browse files
bors[bot]Veykrilerikdesjardins
committed
100: Implement expm1 r=japaric a=Veykril ~~Closes 13~~, closes rust-lang#18 and ~~closes 14~~. ~~I wasn't sure where to put `__expo2(x: f64) -> f64` so I left it in `src/math/cosh.rs` for now.~~ Moved the function into it's own module. Edit: Didn't see that `exp` was already done in a pull request, I'll take it out once rust-lang#90 lands then. 103: implement fma r=japaric a=erikdesjardins closes rust-lang#19 Co-authored-by: Lukas Wirth <[email protected]> Co-authored-by: Erik <[email protected]>
3 parents 07ca284 + a27975e + 6d67a06 commit 4e7504a

File tree

5 files changed

+332
-6
lines changed

5 files changed

+332
-6
lines changed

‎src/lib.rs

-4
Original file line numberDiff line numberDiff line change
@@ -366,7 +366,6 @@ pub trait F64Ext: private::Sealed {
366366
#[cfg(todo)]
367367
fn signum(self) -> Self;
368368

369-
#[cfg(todo)]
370369
fn mul_add(self, a: Self, b: Self) -> Self;
371370

372371
#[cfg(todo)]
@@ -425,7 +424,6 @@ pub trait F64Ext: private::Sealed {
425424
(self.sin(), self.cos())
426425
}
427426

428-
#[cfg(todo)]
429427
fn exp_m1(self) -> Self;
430428

431429
fn ln_1p(self) -> Self;
@@ -485,7 +483,6 @@ impl F64Ext for f64 {
485483
fabs(self)
486484
}
487485

488-
#[cfg(todo)]
489486
#[inline]
490487
fn mul_add(self, a: Self, b: Self) -> Self {
491488
fma(self, a, b)
@@ -604,7 +601,6 @@ impl F64Ext for f64 {
604601
atan2(self, other)
605602
}
606603

607-
#[cfg(todo)]
608604
#[inline]
609605
fn exp_m1(self) -> Self {
610606
expm1(self)

‎src/math/expm1.rs

+125
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
use core::f64;
2+
3+
const O_THRESHOLD: f64 = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
4+
const LN2_HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
5+
const LN2_LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
6+
const INVLN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
7+
/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
8+
const Q1: f64 = -3.33333333333331316428e-02; /* BFA11111 111110F4 */
9+
const Q2: f64 = 1.58730158725481460165e-03; /* 3F5A01A0 19FE5585 */
10+
const Q3: f64 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
11+
const Q4: f64 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
12+
const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
13+
14+
#[allow(warnings)]
15+
pub fn expm1(mut x: f64) -> f64 {
16+
let hi: f64;
17+
let lo: f64;
18+
let k: i32;
19+
let c: f64;
20+
let mut t: f64;
21+
let mut y: f64;
22+
23+
let mut ui = x.to_bits();
24+
let hx = ((ui >> 32) & 0x7fffffff) as u32;
25+
let sign = (ui >> 63) as i32;
26+
27+
/* filter out huge and non-finite argument */
28+
if hx >= 0x4043687A {
29+
/* if |x|>=56*ln2 */
30+
if x.is_nan() {
31+
return x;
32+
}
33+
if sign != 0 {
34+
return -1.0;
35+
}
36+
if x > O_THRESHOLD {
37+
x *= f64::from_bits(0x7fe0000000000000);
38+
return x;
39+
}
40+
}
41+
42+
/* argument reduction */
43+
if hx > 0x3fd62e42 {
44+
/* if |x| > 0.5 ln2 */
45+
if hx < 0x3FF0A2B2 {
46+
/* and |x| < 1.5 ln2 */
47+
if sign == 0 {
48+
hi = x - LN2_HI;
49+
lo = LN2_LO;
50+
k = 1;
51+
} else {
52+
hi = x + LN2_HI;
53+
lo = -LN2_LO;
54+
k = -1;
55+
}
56+
} else {
57+
k = (INVLN2 * x + if sign != 0 { -0.5 } else { 0.5 }) as i32;
58+
t = k as f64;
59+
hi = x - t * LN2_HI; /* t*ln2_hi is exact here */
60+
lo = t * LN2_LO;
61+
}
62+
x = hi - lo;
63+
c = (hi - x) - lo;
64+
} else if hx < 0x3c900000 {
65+
/* |x| < 2**-54, return x */
66+
if hx < 0x00100000 {
67+
force_eval!(x);
68+
}
69+
return x;
70+
} else {
71+
c = 0.0;
72+
k = 0;
73+
}
74+
75+
/* x is now in primary range */
76+
let hfx = 0.5 * x;
77+
let hxs = x * hfx;
78+
let r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
79+
t = 3.0 - r1 * hfx;
80+
let mut e = hxs * ((r1 - t) / (6.0 - x * t));
81+
if k == 0 {
82+
/* c is 0 */
83+
return x - (x * e - hxs);
84+
}
85+
e = x * (e - c) - c;
86+
e -= hxs;
87+
/* exp(x) ~ 2^k (x_reduced - e + 1) */
88+
if k == -1 {
89+
return 0.5 * (x - e) - 0.5;
90+
}
91+
if k == 1 {
92+
if x < -0.25 {
93+
return -2.0 * (e - (x + 0.5));
94+
}
95+
return 1.0 + 2.0 * (x - e);
96+
}
97+
ui = ((0x3ff + k) as u64) << 52; /* 2^k */
98+
let twopk = f64::from_bits(ui);
99+
if k < 0 || k > 56 {
100+
/* suffice to return exp(x)-1 */
101+
y = x - e + 1.0;
102+
if k == 1024 {
103+
y = y * 2.0 * f64::from_bits(0x7fe0000000000000);
104+
} else {
105+
y = y * twopk;
106+
}
107+
return y - 1.0;
108+
}
109+
ui = ((0x3ff - k) as u64) << 52; /* 2^-k */
110+
let uf = f64::from_bits(ui);
111+
if k < 20 {
112+
y = (x - e + (1.0 - uf)) * twopk;
113+
} else {
114+
y = (x - (e + uf) + 1.0) * twopk;
115+
}
116+
y
117+
}
118+
119+
#[cfg(test)]
120+
mod tests {
121+
#[test]
122+
fn sanity_check() {
123+
assert_eq!(super::expm1(1.1), 2.0041660239464334);
124+
}
125+
}

‎src/math/fma.rs

+201
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,201 @@
1+
use core::{f32, f64};
2+
3+
use super::scalbn;
4+
5+
const ZEROINFNAN: i32 = 0x7ff - 0x3ff - 52 - 1;
6+
7+
struct Num {
8+
m: u64,
9+
e: i32,
10+
sign: i32,
11+
}
12+
13+
#[inline]
14+
fn normalize(x: f64) -> Num {
15+
let x1p63: f64 = f64::from_bits(0x43e0000000000000); // 0x1p63 === 2 ^ 63
16+
17+
let mut ix: u64 = x.to_bits();
18+
let mut e: i32 = (ix >> 52) as i32;
19+
let sign: i32 = e & 0x800;
20+
e &= 0x7ff;
21+
if e == 0 {
22+
ix = (x * x1p63).to_bits();
23+
e = (ix >> 52) as i32 & 0x7ff;
24+
e = if e != 0 { e - 63 } else { 0x800 };
25+
}
26+
ix &= (1 << 52) - 1;
27+
ix |= 1 << 52;
28+
ix <<= 1;
29+
e -= 0x3ff + 52 + 1;
30+
Num { m: ix, e, sign }
31+
}
32+
33+
#[inline]
34+
fn mul(x: u64, y: u64) -> (u64, u64) {
35+
let t1: u64;
36+
let t2: u64;
37+
let t3: u64;
38+
let xlo: u64 = x as u32 as u64;
39+
let xhi: u64 = x >> 32;
40+
let ylo: u64 = y as u32 as u64;
41+
let yhi: u64 = y >> 32;
42+
43+
t1 = xlo * ylo;
44+
t2 = xlo * yhi + xhi * ylo;
45+
t3 = xhi * yhi;
46+
let lo = t1 + (t2 << 32);
47+
let hi = t3 + (t2 >> 32) + (t1 > lo) as u64;
48+
(hi, lo)
49+
}
50+
51+
#[inline]
52+
pub fn fma(x: f64, y: f64, z: f64) -> f64 {
53+
let x1p63: f64 = f64::from_bits(0x43e0000000000000); // 0x1p63 === 2 ^ 63
54+
let x0_ffffff8p_63 = f64::from_bits(0x3bfffffff0000000); // 0x0.ffffff8p-63
55+
56+
/* normalize so top 10bits and last bit are 0 */
57+
let nx = normalize(x);
58+
let ny = normalize(y);
59+
let nz = normalize(z);
60+
61+
if nx.e >= ZEROINFNAN || ny.e >= ZEROINFNAN {
62+
return x * y + z;
63+
}
64+
if nz.e >= ZEROINFNAN {
65+
if nz.e > ZEROINFNAN {
66+
/* z==0 */
67+
return x * y + z;
68+
}
69+
return z;
70+
}
71+
72+
/* mul: r = x*y */
73+
let zhi: u64;
74+
let zlo: u64;
75+
let (mut rhi, mut rlo) = mul(nx.m, ny.m);
76+
/* either top 20 or 21 bits of rhi and last 2 bits of rlo are 0 */
77+
78+
/* align exponents */
79+
let mut e: i32 = nx.e + ny.e;
80+
let mut d: i32 = nz.e - e;
81+
/* shift bits z<<=kz, r>>=kr, so kz+kr == d, set e = e+kr (== ez-kz) */
82+
if d > 0 {
83+
if d < 64 {
84+
zlo = nz.m << d;
85+
zhi = nz.m >> 64 - d;
86+
} else {
87+
zlo = 0;
88+
zhi = nz.m;
89+
e = nz.e - 64;
90+
d -= 64;
91+
if d == 0 {
92+
} else if d < 64 {
93+
rlo = rhi << 64 - d | rlo >> d | ((rlo << 64 - d) != 0) as u64;
94+
rhi = rhi >> d;
95+
} else {
96+
rlo = 1;
97+
rhi = 0;
98+
}
99+
}
100+
} else {
101+
zhi = 0;
102+
d = -d;
103+
if d == 0 {
104+
zlo = nz.m;
105+
} else if d < 64 {
106+
zlo = nz.m >> d | ((nz.m << 64 - d) != 0) as u64;
107+
} else {
108+
zlo = 1;
109+
}
110+
}
111+
112+
/* add */
113+
let mut sign: i32 = nx.sign ^ ny.sign;
114+
let samesign: bool = (sign ^ nz.sign) == 0;
115+
let mut nonzero: i32 = 1;
116+
if samesign {
117+
/* r += z */
118+
rlo += zlo;
119+
rhi += zhi + (rlo < zlo) as u64;
120+
} else {
121+
/* r -= z */
122+
let t = rlo;
123+
rlo -= zlo;
124+
rhi = rhi - zhi - (t < rlo) as u64;
125+
if (rhi >> 63) != 0 {
126+
rlo = (-(rlo as i64)) as u64;
127+
rhi = (-(rhi as i64)) as u64 - (rlo != 0) as u64;
128+
sign = (sign == 0) as i32;
129+
}
130+
nonzero = (rhi != 0) as i32;
131+
}
132+
133+
/* set rhi to top 63bit of the result (last bit is sticky) */
134+
if nonzero != 0 {
135+
e += 64;
136+
d = rhi.leading_zeros() as i32 - 1;
137+
/* note: d > 0 */
138+
rhi = rhi << d | rlo >> 64 - d | ((rlo << d) != 0) as u64;
139+
} else if rlo != 0 {
140+
d = rlo.leading_zeros() as i32 - 1;
141+
if d < 0 {
142+
rhi = rlo >> 1 | (rlo & 1);
143+
} else {
144+
rhi = rlo << d;
145+
}
146+
} else {
147+
/* exact +-0 */
148+
return x * y + z;
149+
}
150+
e -= d;
151+
152+
/* convert to double */
153+
let mut i: i64 = rhi as i64; /* i is in [1<<62,(1<<63)-1] */
154+
if sign != 0 {
155+
i = -i;
156+
}
157+
let mut r: f64 = i as f64; /* |r| is in [0x1p62,0x1p63] */
158+
159+
if e < -1022 - 62 {
160+
/* result is subnormal before rounding */
161+
if e == -1022 - 63 {
162+
let mut c: f64 = x1p63;
163+
if sign != 0 {
164+
c = -c;
165+
}
166+
if r == c {
167+
/* min normal after rounding, underflow depends
168+
on arch behaviour which can be imitated by
169+
a double to float conversion */
170+
let fltmin: f32 = (x0_ffffff8p_63 * f32::MIN_POSITIVE as f64 * r) as f32;
171+
return f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * fltmin as f64;
172+
}
173+
/* one bit is lost when scaled, add another top bit to
174+
only round once at conversion if it is inexact */
175+
if (rhi << 53) != 0 {
176+
i = (rhi >> 1 | (rhi & 1) | 1 << 62) as i64;
177+
if sign != 0 {
178+
i = -i;
179+
}
180+
r = i as f64;
181+
r = 2. * r - c; /* remove top bit */
182+
183+
/* raise underflow portably, such that it
184+
cannot be optimized away */
185+
{
186+
let tiny: f64 = f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * r;
187+
r += (tiny * tiny) * (r - r);
188+
}
189+
}
190+
} else {
191+
/* only round once when scaled */
192+
d = 10;
193+
i = ((rhi >> d | ((rhi << 64 - d) != 0) as u64) << d) as i64;
194+
if sign != 0 {
195+
i = -i;
196+
}
197+
r = i as f64;
198+
}
199+
}
200+
scalbn(r, e)
201+
}

‎src/math/mod.rs

+4
Original file line numberDiff line numberDiff line change
@@ -16,12 +16,14 @@ mod exp;
1616
mod exp2;
1717
mod exp2f;
1818
mod expf;
19+
mod expm1;
1920
mod fabs;
2021
mod fabsf;
2122
mod fdim;
2223
mod fdimf;
2324
mod floor;
2425
mod floorf;
26+
mod fma;
2527
mod fmod;
2628
mod fmodf;
2729
mod hypot;
@@ -55,12 +57,14 @@ pub use self::exp::exp;
5557
pub use self::exp2::exp2;
5658
pub use self::exp2f::exp2f;
5759
pub use self::expf::expf;
60+
pub use self::expm1::expm1;
5861
pub use self::fabs::fabs;
5962
pub use self::fabsf::fabsf;
6063
pub use self::fdim::fdim;
6164
pub use self::fdimf::fdimf;
6265
pub use self::floor::floor;
6366
pub use self::floorf::floorf;
67+
pub use self::fma::fma;
6468
pub use self::fmod::fmod;
6569
pub use self::fmodf::fmodf;
6670
pub use self::hypot::hypot;

0 commit comments

Comments
 (0)