Skip to content
This repository was archived by the owner on Apr 28, 2025. It is now read-only.

Commit 75e1b59

Browse files
AmanieuAmanieu
authored
Merge pull request #249 from plugwash/master
2 parents 74a208a + 7862895 commit 75e1b59

File tree

11 files changed

+189
-16
lines changed

11 files changed

+189
-16
lines changed

‎libm/src/math/ceil.rs

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
#![allow(unreachable_code)]
12
use core::f64;
23

34
const TOINT: f64 = 1. / f64::EPSILON;
@@ -15,6 +16,24 @@ pub fn ceil(x: f64) -> f64 {
1516
return unsafe { ::core::intrinsics::ceilf64(x) }
1617
}
1718
}
19+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
20+
{
21+
//use an alternative implementation on x86, because the
22+
//main implementation fails with the x87 FPU used by
23+
//debian i386, probablly due to excess precision issues.
24+
//basic implementation taken from https://github.com/rust-lang/libm/issues/219
25+
use super::fabs;
26+
if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() {
27+
let truncated = x as i64 as f64;
28+
if truncated < x {
29+
return truncated + 1.0;
30+
} else {
31+
return truncated;
32+
}
33+
} else {
34+
return x;
35+
}
36+
}
1837
let u: u64 = x.to_bits();
1938
let e: i64 = (u >> 52 & 0x7ff) as i64;
2039
let y: f64;

‎libm/src/math/floor.rs

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
#![allow(unreachable_code)]
12
use core::f64;
23

34
const TOINT: f64 = 1. / f64::EPSILON;
@@ -15,6 +16,24 @@ pub fn floor(x: f64) -> f64 {
1516
return unsafe { ::core::intrinsics::floorf64(x) }
1617
}
1718
}
19+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
20+
{
21+
//use an alternative implementation on x86, because the
22+
//main implementation fails with the x87 FPU used by
23+
//debian i386, probablly due to excess precision issues.
24+
//basic implementation taken from https://github.com/rust-lang/libm/issues/219
25+
use super::fabs;
26+
if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() {
27+
let truncated = x as i64 as f64;
28+
if truncated > x {
29+
return truncated - 1.0;
30+
} else {
31+
return truncated;
32+
}
33+
} else {
34+
return x;
35+
}
36+
}
1837
let ui = x.to_bits();
1938
let e = ((ui >> 52) & 0x7ff) as i32;
2039

‎libm/src/math/fma.rs

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -218,7 +218,11 @@ mod tests {
218218
-0.00000000000000022204460492503126,
219219
);
220220

221-
assert_eq!(fma(-0.992, -0.992, -0.992), -0.007936000000000007,);
221+
let result = fma(-0.992, -0.992, -0.992);
222+
//force rounding to storage format on x87 to prevent superious errors.
223+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
224+
let result = force_eval!(result);
225+
assert_eq!(result, -0.007936000000000007,);
222226
}
223227

224228
#[test]

‎libm/src/math/j1f.rs

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -369,6 +369,12 @@ mod tests {
369369
}
370370
#[test]
371371
fn test_y1f_2002() {
372-
assert_eq!(y1f(2.0000002_f32), -0.10703229_f32);
372+
//allow slightly different result on x87
373+
let res = y1f(2.0000002_f32);
374+
if cfg!(all(target_arch = "x86", not(target_feature = "sse2"))) && (res == -0.10703231_f32)
375+
{
376+
return;
377+
}
378+
assert_eq!(res, -0.10703229_f32);
373379
}
374380
}

‎libm/src/math/mod.rs

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,6 @@
11
macro_rules! force_eval {
22
($e:expr) => {
3-
unsafe {
4-
::core::ptr::read_volatile(&$e);
5-
}
3+
unsafe { ::core::ptr::read_volatile(&$e) }
64
};
75
}
86

‎libm/src/math/pow.rs

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -484,6 +484,10 @@ mod tests {
484484
let exp = expected(*val);
485485
let res = computed(*val);
486486

487+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
488+
let exp = force_eval!(exp);
489+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
490+
let res = force_eval!(res);
487491
assert!(
488492
if exp.is_nan() {
489493
res.is_nan()

‎libm/src/math/rem_pio2.rs

Lines changed: 18 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,12 @@ pub(crate) fn rem_pio2(x: f64) -> (i32, f64, f64) {
5050

5151
fn medium(x: f64, ix: u32) -> (i32, f64, f64) {
5252
/* rint(x/(pi/2)), Assume round-to-nearest. */
53-
let f_n = x as f64 * INV_PIO2 + TO_INT - TO_INT;
53+
let tmp = x as f64 * INV_PIO2 + TO_INT;
54+
// force rounding of tmp to it's storage format on x87 to avoid
55+
// excess precision issues.
56+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
57+
let tmp = force_eval!(tmp);
58+
let f_n = tmp - TO_INT;
5459
let n = f_n as i32;
5560
let mut r = x - f_n * PIO2_1;
5661
let mut w = f_n * PIO2_1T; /* 1st round, good to 85 bits */
@@ -190,20 +195,28 @@ mod tests {
190195

191196
#[test]
192197
fn test_near_pi() {
198+
let arg = 3.141592025756836;
199+
let arg = force_eval!(arg);
193200
assert_eq!(
194-
rem_pio2(3.141592025756836),
201+
rem_pio2(arg),
195202
(2, -6.278329573009626e-7, -2.1125998133974653e-23)
196203
);
204+
let arg = 3.141592033207416;
205+
let arg = force_eval!(arg);
197206
assert_eq!(
198-
rem_pio2(3.141592033207416),
207+
rem_pio2(arg),
199208
(2, -6.20382377148128e-7, -2.1125998133974653e-23)
200209
);
210+
let arg = 3.141592144966125;
211+
let arg = force_eval!(arg);
201212
assert_eq!(
202-
rem_pio2(3.141592144966125),
213+
rem_pio2(arg),
203214
(2, -5.086236681942706e-7, -2.1125998133974653e-23)
204215
);
216+
let arg = 3.141592979431152;
217+
let arg = force_eval!(arg);
205218
assert_eq!(
206-
rem_pio2(3.141592979431152),
219+
rem_pio2(arg),
207220
(2, 3.2584135866119817e-7, -2.1125998133974653e-23)
208221
);
209222
}

‎libm/src/math/rem_pio2f.rs

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,12 @@ pub(crate) fn rem_pio2f(x: f32) -> (i32, f64) {
4343
if ix < 0x4dc90fdb {
4444
/* |x| ~< 2^28*(pi/2), medium size */
4545
/* Use a specialized rint() to get fn. Assume round-to-nearest. */
46-
let f_n = x64 * INV_PIO2 + TOINT - TOINT;
46+
let tmp = x64 * INV_PIO2 + TOINT;
47+
// force rounding of tmp to it's storage format on x87 to avoid
48+
// excess precision issues.
49+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
50+
let tmp = force_eval!(tmp);
51+
let f_n = tmp - TOINT;
4752
return (f_n as i32, x64 - f_n * PIO2_1 - f_n * PIO2_1T);
4853
}
4954
if ix >= 0x7f800000 {

‎libm/src/math/sin.rs

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,5 +81,8 @@ pub fn sin(x: f64) -> f64 {
8181
fn test_near_pi() {
8282
let x = f64::from_bits(0x400921fb000FD5DD); // 3.141592026217707
8383
let sx = f64::from_bits(0x3ea50d15ced1a4a2); // 6.273720864039205e-7
84-
assert_eq!(sin(x), sx);
84+
let result = sin(x);
85+
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
86+
let result = force_eval!(result);
87+
assert_eq!(result, sx);
8588
}

‎libm/src/math/sincos.rs

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -57,3 +57,77 @@ pub fn sincos(x: f64) -> (f64, f64) {
5757
_ => (0.0, 1.0),
5858
}
5959
}
60+
61+
// These tests are based on those from sincosf.rs
62+
#[cfg(test)]
63+
mod tests {
64+
use super::sincos;
65+
66+
const TOLERANCE: f64 = 1e-6;
67+
68+
#[test]
69+
fn with_pi() {
70+
let (s, c) = sincos(core::f64::consts::PI);
71+
assert!(
72+
(s - 0.0).abs() < TOLERANCE,
73+
"|{} - {}| = {} >= {}",
74+
s,
75+
0.0,
76+
(s - 0.0).abs(),
77+
TOLERANCE
78+
);
79+
assert!(
80+
(c + 1.0).abs() < TOLERANCE,
81+
"|{} + {}| = {} >= {}",
82+
c,
83+
1.0,
84+
(s + 1.0).abs(),
85+
TOLERANCE
86+
);
87+
}
88+
89+
#[test]
90+
fn rotational_symmetry() {
91+
use core::f64::consts::PI;
92+
const N: usize = 24;
93+
for n in 0..N {
94+
let theta = 2. * PI * (n as f64) / (N as f64);
95+
let (s, c) = sincos(theta);
96+
let (s_plus, c_plus) = sincos(theta + 2. * PI);
97+
let (s_minus, c_minus) = sincos(theta - 2. * PI);
98+
99+
assert!(
100+
(s - s_plus).abs() < TOLERANCE,
101+
"|{} - {}| = {} >= {}",
102+
s,
103+
s_plus,
104+
(s - s_plus).abs(),
105+
TOLERANCE
106+
);
107+
assert!(
108+
(s - s_minus).abs() < TOLERANCE,
109+
"|{} - {}| = {} >= {}",
110+
s,
111+
s_minus,
112+
(s - s_minus).abs(),
113+
TOLERANCE
114+
);
115+
assert!(
116+
(c - c_plus).abs() < TOLERANCE,
117+
"|{} - {}| = {} >= {}",
118+
c,
119+
c_plus,
120+
(c - c_plus).abs(),
121+
TOLERANCE
122+
);
123+
assert!(
124+
(c - c_minus).abs() < TOLERANCE,
125+
"|{} - {}| = {} >= {}",
126+
c,
127+
c_minus,
128+
(c - c_minus).abs(),
129+
TOLERANCE
130+
);
131+
}
132+
}
133+
}

0 commit comments

Comments
 (0)