Add forced rounding to storage format for x87 to rem_pio2.rs as well.

This commit is contained in:
Peter Michael Green 2021-12-21 23:53:06 +00:00
parent 874209b56c
commit e95ea2b11d
2 changed files with 79 additions and 1 deletions

View file

@ -50,7 +50,11 @@ pub(crate) fn rem_pio2(x: f64) -> (i32, f64, f64) {
fn medium(x: f64, ix: u32) -> (i32, f64, f64) {
/* rint(x/(pi/2)), Assume round-to-nearest. */
let f_n = x as f64 * INV_PIO2 + TO_INT - TO_INT;
let tmp = x as f64 * INV_PIO2 + TO_INT;
// force rounding of tmp to it's storage format on x87 to avoid
// excess precision issues.
#[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]force_eval!(tmp);
let f_n = tmp - TO_INT;
let n = f_n as i32;
let mut r = x - f_n * PIO2_1;
let mut w = f_n * PIO2_1T; /* 1st round, good to 85 bits */

View file

@ -57,3 +57,77 @@ pub fn sincos(x: f64) -> (f64, f64) {
_ => (0.0, 1.0),
}
}
// These tests are based on those from sincosf.rs
#[cfg(test)]
mod tests {
use super::sincos;
const TOLERANCE: f64 = 1e-6;
#[test]
fn with_pi() {
let (s, c) = sincos(core::f64::consts::PI);
assert!(
(s - 0.0).abs() < TOLERANCE,
"|{} - {}| = {} >= {}",
s,
0.0,
(s - 0.0).abs(),
TOLERANCE
);
assert!(
(c + 1.0).abs() < TOLERANCE,
"|{} + {}| = {} >= {}",
c,
1.0,
(s + 1.0).abs(),
TOLERANCE
);
}
#[test]
fn rotational_symmetry() {
use core::f64::consts::PI;
const N: usize = 24;
for n in 0..N {
let theta = 2. * PI * (n as f64) / (N as f64);
let (s, c) = sincos(theta);
let (s_plus, c_plus) = sincos(theta + 2. * PI);
let (s_minus, c_minus) = sincos(theta - 2. * PI);
assert!(
(s - s_plus).abs() < TOLERANCE,
"|{} - {}| = {} >= {}",
s,
s_plus,
(s - s_plus).abs(),
TOLERANCE
);
assert!(
(s - s_minus).abs() < TOLERANCE,
"|{} - {}| = {} >= {}",
s,
s_minus,
(s - s_minus).abs(),
TOLERANCE
);
assert!(
(c - c_plus).abs() < TOLERANCE,
"|{} - {}| = {} >= {}",
c,
c_plus,
(c - c_plus).abs(),
TOLERANCE
);
assert!(
(c - c_minus).abs() < TOLERANCE,
"|{} - {}| = {} >= {}",
c,
c_minus,
(c - c_minus).abs(),
TOLERANCE
);
}
}
}