fix bug in rem_pio2_large

This commit is contained in:
Andrey Zgarbul 2018-07-14 08:13:35 +03:00
parent f40b8096aa
commit 45d9d8caa6
5 changed files with 51 additions and 51 deletions

View file

@ -1,4 +1,4 @@
use super::{cosdf, sindf, rem_pio2f};
use super::{k_cosf, k_sinf, rem_pio2f};
use core::f64::consts::FRAC_PI_2;
@ -14,52 +14,52 @@ pub fn cosf(x: f32) -> f32 {
let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
let mut ix = x.to_bits();
let sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
let mut ix = x.to_bits();
let sign = (ix >> 31) != 0;
ix &= 0x7fffffff;
if ix <= 0x3f490fda { /* |x| ~<= pi/4 */
if ix < 0x39800000 { /* |x| < 2**-12 */
/* raise inexact if x != 0 */
force_eval!(x + x1p120);
return 1.;
}
return cosdf(x64);
}
if ix <= 0x407b53d1 { /* |x| ~<= 5*pi/4 */
if ix > 0x4016cbe3 { /* |x| ~> 3*pi/4 */
return -cosdf(if sign { x64+C2_PIO2 } else { x64-C2_PIO2 });
if ix <= 0x3f490fda { /* |x| ~<= pi/4 */
if ix < 0x39800000 { /* |x| < 2**-12 */
/* raise inexact if x != 0 */
force_eval!(x + x1p120);
return 1.;
}
return k_cosf(x64);
}
if ix <= 0x407b53d1 { /* |x| ~<= 5*pi/4 */
if ix > 0x4016cbe3 { /* |x| ~> 3*pi/4 */
return -k_cosf(if sign { x64+C2_PIO2 } else { x64-C2_PIO2 });
} else {
if sign {
return sindf(x64 + C1_PIO2);
} else {
return sindf(C1_PIO2 - x64);
if sign {
return k_sinf(x64 + C1_PIO2);
} else {
return k_sinf(C1_PIO2 - x64);
}
}
}
if ix <= 0x40e231d5 { /* |x| ~<= 9*pi/4 */
if ix > 0x40afeddf { /* |x| ~> 7*pi/4 */
return cosdf(if sign { x64+C4_PIO2 } else { x64-C4_PIO2 });
} else {
if sign {
return sindf(-x64 - C3_PIO2);
} else {
return sindf(x64 - C3_PIO2);
}
}
if ix <= 0x40e231d5 { /* |x| ~<= 9*pi/4 */
if ix > 0x40afeddf { /* |x| ~> 7*pi/4 */
return k_cosf(if sign { x64+C4_PIO2 } else { x64-C4_PIO2 });
} else {
if sign {
return k_sinf(-x64 - C3_PIO2);
} else {
return k_sinf(x64 - C3_PIO2);
}
}
}
/* cos(Inf or NaN) is NaN */
if ix >= 0x7f800000 {
return x-x;
}
}
/* general argument reduction needed */
let (n, y) = rem_pio2f(x);
/* cos(Inf or NaN) is NaN */
if ix >= 0x7f800000 {
return x-x;
}
/* general argument reduction needed */
let (n, y) = rem_pio2f(x);
match n&3 {
0 => { cosdf( y) },
1 => { sindf(-y) },
2 => { -cosdf( y) },
_ => { sindf( y) },
0 => { k_cosf( y) },
1 => { k_sinf(-y) },
2 => { -k_cosf( y) },
_ => { k_sinf( y) },
}
}

View file

@ -5,7 +5,7 @@ const C2 : f64 = -0.00138867637746099294692; /* -0x16c087e80f1e27.0p-62 */
const C3 : f64 = 0.0000243904487962774090654; /* 0x199342e0ee5069.0p-68 */
#[inline]
pub(crate) fn cosdf(x : f64) -> f32 {
pub(crate) fn k_cosf(x : f64) -> f32 {
let z = x*x;
let w = z*z;
let r = C2+z*C3;

View file

@ -5,7 +5,7 @@ const S3 : f64 = -0.000198393348360966317347; /* -0x1a00f9e2cae774.0p-65 */
const S4 : f64 = 0.0000027183114939898219064; /* 0x16cd878c3b46a7.0p-71 */
#[inline]
pub(crate) fn sindf(x : f64) -> f32 {
pub(crate) fn k_sinf(x : f64) -> f32 {
let z = x*x;
let w = z*z;
let r = S3 + z*S4;

View file

@ -36,14 +36,14 @@ pub use self::{
truncf::truncf,
};
mod sindf;
mod cosdf;
mod k_cosf;
mod k_sinf;
mod rem_pio2f;
mod rem_pio2_large;
use self::{
sindf::sindf,
cosdf::cosdf,
k_cosf::k_cosf,
k_sinf::k_sinf,
rem_pio2f::rem_pio2f,
rem_pio2_large::rem_pio2_large,
};

View file

@ -1,5 +1,5 @@
use super::scalbn;
use super::floor;
use math::scalbn;
use math::floor;
// initial value for jk
const INIT_JK : [usize; 4] = [3,4,4,6];
@ -263,7 +263,7 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize)
let mut fw : f64;
let mut n : i32;
let mut ih : i32;
let mut z = 0f64;
let mut z : f64;
let mut f : [f64;20] = [0.;20];
let mut fq : [f64;20] = [0.;20];
let mut q : [f64;20] = [0.;20];
@ -308,7 +308,7 @@ pub(crate) fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize)
'recompute: loop {
/* distill q[] into iq[] reversingly */
let mut i = 0i32;
let mut z = q[jz];
z = q[jz];
for j in (1..=jz).rev() {
fw = (x1p_24*z) as i32 as f64;
iq[i as usize] = (z - x1p24*fw) as i32;