Fix errors reported by Clippy in libm
This commit is contained in:
Trevor Gross
parent
dbf8a4ebe5
commit
c2402e4d3c
29 changed files with 307 additions and 306 deletions
|
|
@ -5,12 +5,15 @@
|
|||
#![allow(clippy::assign_op_pattern)]
|
||||
#![allow(clippy::deprecated_cfg_attr)]
|
||||
#![allow(clippy::eq_op)]
|
||||
#![allow(clippy::excessive_precision)]
|
||||
#![allow(clippy::float_cmp)]
|
||||
#![allow(clippy::int_plus_one)]
|
||||
#![allow(clippy::many_single_char_names)]
|
||||
#![allow(clippy::mixed_case_hex_literals)]
|
||||
#![allow(clippy::needless_late_init)]
|
||||
#![allow(clippy::needless_return)]
|
||||
#![allow(clippy::unreadable_literal)]
|
||||
#![allow(clippy::zero_divided_by_zero)]
|
||||
|
||||
mod libm_helper;
|
||||
mod math;
|
||||
|
|
|
|||
|
|
@ -90,7 +90,7 @@ pub fn asin(mut x: f64) -> f64 {
|
|||
/* |x| < 0.5 */
|
||||
if ix < 0x3fe00000 {
|
||||
/* if 0x1p-1022 <= |x| < 0x1p-26, avoid raising underflow */
|
||||
if ix < 0x3e500000 && ix >= 0x00100000 {
|
||||
if (0x00100000..0x3e500000).contains(&ix) {
|
||||
return x;
|
||||
} else {
|
||||
return x + x * comp_r(x * x);
|
||||
|
|
|
|||
|
|
@ -54,7 +54,7 @@ pub fn asinf(mut x: f32) -> f32 {
|
|||
if ix < 0x3f000000 {
|
||||
/* |x| < 0.5 */
|
||||
/* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */
|
||||
if (ix < 0x39800000) && (ix >= 0x00800000) {
|
||||
if (0x00800000..0x39800000).contains(&ix) {
|
||||
return x;
|
||||
}
|
||||
return x + x * r(x * x);
|
||||
|
|
|
|||
|
|
@ -42,9 +42,9 @@ pub fn atan2f(y: f32, x: f32) -> f32 {
|
|||
/* when y = 0 */
|
||||
if iy == 0 {
|
||||
return match m {
|
||||
0 | 1 => y, /* atan(+-0,+anything)=+-0 */
|
||||
2 => PI, /* atan(+0,-anything) = pi */
|
||||
3 | _ => -PI, /* atan(-0,-anything) =-pi */
|
||||
0 | 1 => y, /* atan(+-0,+anything)=+-0 */
|
||||
2 => PI, /* atan(+0,-anything) = pi */
|
||||
_ => -PI, /* atan(-0,-anything) =-pi */
|
||||
};
|
||||
}
|
||||
/* when x = 0 */
|
||||
|
|
@ -55,17 +55,17 @@ pub fn atan2f(y: f32, x: f32) -> f32 {
|
|||
if ix == 0x7f800000 {
|
||||
return if iy == 0x7f800000 {
|
||||
match m {
|
||||
0 => PI / 4., /* atan(+INF,+INF) */
|
||||
1 => -PI / 4., /* atan(-INF,+INF) */
|
||||
2 => 3. * PI / 4., /* atan(+INF,-INF)*/
|
||||
3 | _ => -3. * PI / 4., /* atan(-INF,-INF)*/
|
||||
0 => PI / 4., /* atan(+INF,+INF) */
|
||||
1 => -PI / 4., /* atan(-INF,+INF) */
|
||||
2 => 3. * PI / 4., /* atan(+INF,-INF)*/
|
||||
_ => -3. * PI / 4., /* atan(-INF,-INF)*/
|
||||
}
|
||||
} else {
|
||||
match m {
|
||||
0 => 0., /* atan(+...,+INF) */
|
||||
1 => -0., /* atan(-...,+INF) */
|
||||
2 => PI, /* atan(+...,-INF) */
|
||||
3 | _ => -PI, /* atan(-...,-INF) */
|
||||
0 => 0., /* atan(+...,+INF) */
|
||||
1 => -0., /* atan(-...,+INF) */
|
||||
2 => PI, /* atan(+...,-INF) */
|
||||
_ => -PI, /* atan(-...,-INF) */
|
||||
}
|
||||
};
|
||||
}
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@ pub fn atanhf(mut x: f32) -> f32 {
|
|||
if u < 0x3f800000 - (32 << 23) {
|
||||
/* handle underflow */
|
||||
if u < (1 << 23) {
|
||||
force_eval!((x * x) as f32);
|
||||
force_eval!(x * x);
|
||||
}
|
||||
} else {
|
||||
/* |x| < 0.5, up to 1.7ulp error */
|
||||
|
|
|
|||
|
|
@ -95,7 +95,7 @@ pub fn exp2f(mut x: f32) -> f32 {
|
|||
/* NaN */
|
||||
return x;
|
||||
}
|
||||
if ui >= 0x43000000 && ui < 0x80000000 {
|
||||
if (0x43000000..0x80000000).contains(&ui) {
|
||||
/* x >= 128 */
|
||||
x *= x1p127;
|
||||
return x;
|
||||
|
|
@ -127,7 +127,7 @@ pub fn exp2f(mut x: f32) -> f32 {
|
|||
let z: f64 = (x - uf) as f64;
|
||||
/* Compute r = exp2(y) = exp2ft[i0] * p(z). */
|
||||
let r: f64 = f64::from_bits(i!(EXP2FT, i0 as usize));
|
||||
let t: f64 = r as f64 * z;
|
||||
let t: f64 = r * z;
|
||||
let r: f64 = r + t * (p1 as f64 + z * p2 as f64) + t * (z * z) * (p3 as f64 + z * p4 as f64);
|
||||
|
||||
/* Scale by 2**k */
|
||||
|
|
|
|||
|
|
@ -115,7 +115,7 @@ pub fn expm1(mut x: f64) -> f64 {
|
|||
}
|
||||
ui = ((0x3ff + k) as u64) << 52; /* 2^k */
|
||||
let twopk = f64::from_bits(ui);
|
||||
if k < 0 || k > 56 {
|
||||
if !(0..=56).contains(&k) {
|
||||
/* suffice to return exp(x)-1 */
|
||||
y = x - e + 1.0;
|
||||
if k == 1024 {
|
||||
|
|
|
|||
|
|
@ -115,7 +115,7 @@ pub fn expm1f(mut x: f32) -> f32 {
|
|||
return 1. + 2. * (x - e);
|
||||
}
|
||||
let twopk = f32::from_bits(((0x7f + k) << 23) as u32); /* 2^k */
|
||||
if (k < 0) || (k > 56) {
|
||||
if !(0..=56).contains(&k) {
|
||||
/* suffice to return exp(x)-1 */
|
||||
let mut y = x - e + 1.;
|
||||
if k == 128 {
|
||||
|
|
|
|||
|
|
@ -1,5 +1,3 @@
|
|||
use core::u64;
|
||||
|
||||
/// Absolute value (magnitude) (f64)
|
||||
/// Calculates the absolute value (magnitude) of the argument `x`,
|
||||
/// by direct manipulation of the bit representation of `x`.
|
||||
|
|
|
|||
|
|
@ -3,9 +3,9 @@ use core::f64;
|
|||
/// Positive difference (f64)
|
||||
///
|
||||
/// Determines the positive difference between arguments, returning:
|
||||
/// * x - y if x > y, or
|
||||
/// * +0 if x <= y, or
|
||||
/// * NAN if either argument is NAN.
|
||||
/// * x - y if x > y, or
|
||||
/// * +0 if x <= y, or
|
||||
/// * NAN if either argument is NAN.
|
||||
///
|
||||
/// A range error may occur.
|
||||
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
|
||||
|
|
|
|||
|
|
@ -3,9 +3,9 @@ use core::f32;
|
|||
/// Positive difference (f32)
|
||||
///
|
||||
/// Determines the positive difference between arguments, returning:
|
||||
/// * x - y if x > y, or
|
||||
/// * +0 if x <= y, or
|
||||
/// * NAN if either argument is NAN.
|
||||
/// * x - y if x > y, or
|
||||
/// * +0 if x <= y, or
|
||||
/// * NAN if either argument is NAN.
|
||||
///
|
||||
/// A range error may occur.
|
||||
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
|
||||
|
|
|
|||
|
|
@ -71,7 +71,7 @@ pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 {
|
|||
underflow may not be raised correctly, example:
|
||||
fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
|
||||
*/
|
||||
if e < 0x3ff - 126 && e >= 0x3ff - 149 && fetestexcept(FE_INEXACT) != 0 {
|
||||
if ((0x3ff - 149)..(0x3ff - 126)).contains(&e) && fetestexcept(FE_INEXACT) != 0 {
|
||||
feclearexcept(FE_INEXACT);
|
||||
// prevent `xy + vz` from being CSE'd with `xy + z` above
|
||||
let vz: f32 = unsafe { read_volatile(&z) };
|
||||
|
|
|
|||
|
|
@ -1,5 +1,3 @@
|
|||
use core::u64;
|
||||
|
||||
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
|
||||
pub fn fmod(x: f64, y: f64) -> f64 {
|
||||
let mut uxi = x.to_bits();
|
||||
|
|
@ -74,7 +72,7 @@ pub fn fmod(x: f64, y: f64) -> f64 {
|
|||
} else {
|
||||
uxi >>= -ex + 1;
|
||||
}
|
||||
uxi |= (sx as u64) << 63;
|
||||
uxi |= sx << 63;
|
||||
|
||||
f64::from_bits(uxi)
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
use core::{f32, u32};
|
||||
use core::f32;
|
||||
|
||||
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
|
||||
pub fn fmodf(x: f32, y: f32) -> f32 {
|
||||
|
|
|
|||
|
|
@ -21,7 +21,7 @@ pub fn ilogb(x: f64) -> i32 {
|
|||
e
|
||||
} else if e == 0x7ff {
|
||||
force_eval!(0.0 / 0.0);
|
||||
if (i << 12) != 0 { FP_ILOGBNAN } else { i32::max_value() }
|
||||
if (i << 12) != 0 { FP_ILOGBNAN } else { i32::MAX }
|
||||
} else {
|
||||
e - 0x3ff
|
||||
}
|
||||
|
|
|
|||
|
|
@ -21,7 +21,7 @@ pub fn ilogbf(x: f32) -> i32 {
|
|||
e
|
||||
} else if e == 0xff {
|
||||
force_eval!(0.0 / 0.0);
|
||||
if (i << 9) != 0 { FP_ILOGBNAN } else { i32::max_value() }
|
||||
if (i << 9) != 0 { FP_ILOGBNAN } else { i32::MAX }
|
||||
} else {
|
||||
e - 0x7f
|
||||
}
|
||||
|
|
|
|||
|
|
@ -104,7 +104,8 @@ pub fn jn(n: i32, mut x: f64) -> f64 {
|
|||
0 => -cos(x) + sin(x),
|
||||
1 => -cos(x) - sin(x),
|
||||
2 => cos(x) - sin(x),
|
||||
3 | _ => cos(x) + sin(x),
|
||||
// 3
|
||||
_ => cos(x) + sin(x),
|
||||
};
|
||||
b = INVSQRTPI * temp / sqrt(x);
|
||||
} else {
|
||||
|
|
@ -118,130 +119,128 @@ pub fn jn(n: i32, mut x: f64) -> f64 {
|
|||
a = temp;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if ix < 0x3e100000 {
|
||||
/* x < 2**-29 */
|
||||
/* x is tiny, return the first Taylor expansion of J(n,x)
|
||||
* J(n,x) = 1/n!*(x/2)^n - ...
|
||||
*/
|
||||
if nm1 > 32 {
|
||||
/* underflow */
|
||||
b = 0.0;
|
||||
} else {
|
||||
temp = x * 0.5;
|
||||
b = temp;
|
||||
a = 1.0;
|
||||
i = 2;
|
||||
while i <= nm1 + 1 {
|
||||
a *= i as f64; /* a = n! */
|
||||
b *= temp; /* b = (x/2)^n */
|
||||
i += 1;
|
||||
}
|
||||
b = b / a;
|
||||
}
|
||||
} else if ix < 0x3e100000 {
|
||||
/* x < 2**-29 */
|
||||
/* x is tiny, return the first Taylor expansion of J(n,x)
|
||||
* J(n,x) = 1/n!*(x/2)^n - ...
|
||||
*/
|
||||
if nm1 > 32 {
|
||||
/* underflow */
|
||||
b = 0.0;
|
||||
} else {
|
||||
/* use backward recurrence */
|
||||
/* x x^2 x^2
|
||||
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
|
||||
* 2n - 2(n+1) - 2(n+2)
|
||||
*
|
||||
* 1 1 1
|
||||
* (for large x) = ---- ------ ------ .....
|
||||
* 2n 2(n+1) 2(n+2)
|
||||
* -- - ------ - ------ -
|
||||
* x x x
|
||||
*
|
||||
* Let w = 2n/x and h=2/x, then the above quotient
|
||||
* is equal to the continued fraction:
|
||||
* 1
|
||||
* = -----------------------
|
||||
* 1
|
||||
* w - -----------------
|
||||
* 1
|
||||
* w+h - ---------
|
||||
* w+2h - ...
|
||||
*
|
||||
* To determine how many terms needed, let
|
||||
* Q(0) = w, Q(1) = w(w+h) - 1,
|
||||
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
|
||||
* When Q(k) > 1e4 good for single
|
||||
* When Q(k) > 1e9 good for double
|
||||
* When Q(k) > 1e17 good for quadruple
|
||||
*/
|
||||
/* determine k */
|
||||
let mut t: f64;
|
||||
let mut q0: f64;
|
||||
let mut q1: f64;
|
||||
let mut w: f64;
|
||||
let h: f64;
|
||||
let mut z: f64;
|
||||
let mut tmp: f64;
|
||||
let nf: f64;
|
||||
|
||||
let mut k: i32;
|
||||
|
||||
nf = (nm1 as f64) + 1.0;
|
||||
w = 2.0 * nf / x;
|
||||
h = 2.0 / x;
|
||||
z = w + h;
|
||||
q0 = w;
|
||||
q1 = w * z - 1.0;
|
||||
k = 1;
|
||||
while q1 < 1.0e9 {
|
||||
k += 1;
|
||||
z += h;
|
||||
tmp = z * q1 - q0;
|
||||
q0 = q1;
|
||||
q1 = tmp;
|
||||
temp = x * 0.5;
|
||||
b = temp;
|
||||
a = 1.0;
|
||||
i = 2;
|
||||
while i <= nm1 + 1 {
|
||||
a *= i as f64; /* a = n! */
|
||||
b *= temp; /* b = (x/2)^n */
|
||||
i += 1;
|
||||
}
|
||||
t = 0.0;
|
||||
i = k;
|
||||
while i >= 0 {
|
||||
t = 1.0 / (2.0 * ((i as f64) + nf) / x - t);
|
||||
b = b / a;
|
||||
}
|
||||
} else {
|
||||
/* use backward recurrence */
|
||||
/* x x^2 x^2
|
||||
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
|
||||
* 2n - 2(n+1) - 2(n+2)
|
||||
*
|
||||
* 1 1 1
|
||||
* (for large x) = ---- ------ ------ .....
|
||||
* 2n 2(n+1) 2(n+2)
|
||||
* -- - ------ - ------ -
|
||||
* x x x
|
||||
*
|
||||
* Let w = 2n/x and h=2/x, then the above quotient
|
||||
* is equal to the continued fraction:
|
||||
* 1
|
||||
* = -----------------------
|
||||
* 1
|
||||
* w - -----------------
|
||||
* 1
|
||||
* w+h - ---------
|
||||
* w+2h - ...
|
||||
*
|
||||
* To determine how many terms needed, let
|
||||
* Q(0) = w, Q(1) = w(w+h) - 1,
|
||||
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
|
||||
* When Q(k) > 1e4 good for single
|
||||
* When Q(k) > 1e9 good for double
|
||||
* When Q(k) > 1e17 good for quadruple
|
||||
*/
|
||||
/* determine k */
|
||||
let mut t: f64;
|
||||
let mut q0: f64;
|
||||
let mut q1: f64;
|
||||
let mut w: f64;
|
||||
let h: f64;
|
||||
let mut z: f64;
|
||||
let mut tmp: f64;
|
||||
let nf: f64;
|
||||
|
||||
let mut k: i32;
|
||||
|
||||
nf = (nm1 as f64) + 1.0;
|
||||
w = 2.0 * nf / x;
|
||||
h = 2.0 / x;
|
||||
z = w + h;
|
||||
q0 = w;
|
||||
q1 = w * z - 1.0;
|
||||
k = 1;
|
||||
while q1 < 1.0e9 {
|
||||
k += 1;
|
||||
z += h;
|
||||
tmp = z * q1 - q0;
|
||||
q0 = q1;
|
||||
q1 = tmp;
|
||||
}
|
||||
t = 0.0;
|
||||
i = k;
|
||||
while i >= 0 {
|
||||
t = 1.0 / (2.0 * ((i as f64) + nf) / x - t);
|
||||
i -= 1;
|
||||
}
|
||||
a = t;
|
||||
b = 1.0;
|
||||
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
|
||||
* Hence, if n*(log(2n/x)) > ...
|
||||
* single 8.8722839355e+01
|
||||
* double 7.09782712893383973096e+02
|
||||
* long double 1.1356523406294143949491931077970765006170e+04
|
||||
* then recurrent value may overflow and the result is
|
||||
* likely underflow to zero
|
||||
*/
|
||||
tmp = nf * log(fabs(w));
|
||||
if tmp < 7.09782712893383973096e+02 {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = b * (2.0 * (i as f64)) / x - a;
|
||||
a = temp;
|
||||
i -= 1;
|
||||
}
|
||||
a = t;
|
||||
b = 1.0;
|
||||
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
|
||||
* Hence, if n*(log(2n/x)) > ...
|
||||
* single 8.8722839355e+01
|
||||
* double 7.09782712893383973096e+02
|
||||
* long double 1.1356523406294143949491931077970765006170e+04
|
||||
* then recurrent value may overflow and the result is
|
||||
* likely underflow to zero
|
||||
*/
|
||||
tmp = nf * log(fabs(w));
|
||||
if tmp < 7.09782712893383973096e+02 {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = b * (2.0 * (i as f64)) / x - a;
|
||||
a = temp;
|
||||
i -= 1;
|
||||
}
|
||||
} else {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = b * (2.0 * (i as f64)) / x - a;
|
||||
a = temp;
|
||||
/* scale b to avoid spurious overflow */
|
||||
let x1p500 = f64::from_bits(0x5f30000000000000); // 0x1p500 == 2^500
|
||||
if b > x1p500 {
|
||||
a /= b;
|
||||
t /= b;
|
||||
b = 1.0;
|
||||
}
|
||||
i -= 1;
|
||||
} else {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = b * (2.0 * (i as f64)) / x - a;
|
||||
a = temp;
|
||||
/* scale b to avoid spurious overflow */
|
||||
let x1p500 = f64::from_bits(0x5f30000000000000); // 0x1p500 == 2^500
|
||||
if b > x1p500 {
|
||||
a /= b;
|
||||
t /= b;
|
||||
b = 1.0;
|
||||
}
|
||||
i -= 1;
|
||||
}
|
||||
z = j0(x);
|
||||
w = j1(x);
|
||||
if fabs(z) >= fabs(w) {
|
||||
b = t * z / b;
|
||||
} else {
|
||||
b = t * w / a;
|
||||
}
|
||||
}
|
||||
z = j0(x);
|
||||
w = j1(x);
|
||||
if fabs(z) >= fabs(w) {
|
||||
b = t * z / b;
|
||||
} else {
|
||||
b = t * w / a;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -315,7 +314,8 @@ pub fn yn(n: i32, x: f64) -> f64 {
|
|||
0 => -sin(x) - cos(x),
|
||||
1 => -sin(x) + cos(x),
|
||||
2 => sin(x) + cos(x),
|
||||
3 | _ => sin(x) - cos(x),
|
||||
// 3
|
||||
_ => sin(x) - cos(x),
|
||||
};
|
||||
b = INVSQRTPI * temp / sqrt(x);
|
||||
} else {
|
||||
|
|
|
|||
|
|
@ -64,128 +64,126 @@ pub fn jnf(n: i32, mut x: f32) -> f32 {
|
|||
b = b * (2.0 * (i as f32) / x) - a;
|
||||
a = temp;
|
||||
}
|
||||
} else if ix < 0x35800000 {
|
||||
/* x < 2**-20 */
|
||||
/* x is tiny, return the first Taylor expansion of J(n,x)
|
||||
* J(n,x) = 1/n!*(x/2)^n - ...
|
||||
*/
|
||||
if nm1 > 8 {
|
||||
/* underflow */
|
||||
nm1 = 8;
|
||||
}
|
||||
temp = 0.5 * x;
|
||||
b = temp;
|
||||
a = 1.0;
|
||||
i = 2;
|
||||
while i <= nm1 + 1 {
|
||||
a *= i as f32; /* a = n! */
|
||||
b *= temp; /* b = (x/2)^n */
|
||||
i += 1;
|
||||
}
|
||||
b = b / a;
|
||||
} else {
|
||||
if ix < 0x35800000 {
|
||||
/* x < 2**-20 */
|
||||
/* x is tiny, return the first Taylor expansion of J(n,x)
|
||||
* J(n,x) = 1/n!*(x/2)^n - ...
|
||||
*/
|
||||
if nm1 > 8 {
|
||||
/* underflow */
|
||||
nm1 = 8;
|
||||
}
|
||||
temp = 0.5 * x;
|
||||
b = temp;
|
||||
a = 1.0;
|
||||
i = 2;
|
||||
while i <= nm1 + 1 {
|
||||
a *= i as f32; /* a = n! */
|
||||
b *= temp; /* b = (x/2)^n */
|
||||
i += 1;
|
||||
}
|
||||
b = b / a;
|
||||
} else {
|
||||
/* use backward recurrence */
|
||||
/* x x^2 x^2
|
||||
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
|
||||
* 2n - 2(n+1) - 2(n+2)
|
||||
*
|
||||
* 1 1 1
|
||||
* (for large x) = ---- ------ ------ .....
|
||||
* 2n 2(n+1) 2(n+2)
|
||||
* -- - ------ - ------ -
|
||||
* x x x
|
||||
*
|
||||
* Let w = 2n/x and h=2/x, then the above quotient
|
||||
* is equal to the continued fraction:
|
||||
* 1
|
||||
* = -----------------------
|
||||
* 1
|
||||
* w - -----------------
|
||||
* 1
|
||||
* w+h - ---------
|
||||
* w+2h - ...
|
||||
*
|
||||
* To determine how many terms needed, let
|
||||
* Q(0) = w, Q(1) = w(w+h) - 1,
|
||||
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
|
||||
* When Q(k) > 1e4 good for single
|
||||
* When Q(k) > 1e9 good for double
|
||||
* When Q(k) > 1e17 good for quadruple
|
||||
*/
|
||||
/* determine k */
|
||||
let mut t: f32;
|
||||
let mut q0: f32;
|
||||
let mut q1: f32;
|
||||
let mut w: f32;
|
||||
let h: f32;
|
||||
let mut z: f32;
|
||||
let mut tmp: f32;
|
||||
let nf: f32;
|
||||
let mut k: i32;
|
||||
/* use backward recurrence */
|
||||
/* x x^2 x^2
|
||||
* J(n,x)/J(n-1,x) = ---- ------ ------ .....
|
||||
* 2n - 2(n+1) - 2(n+2)
|
||||
*
|
||||
* 1 1 1
|
||||
* (for large x) = ---- ------ ------ .....
|
||||
* 2n 2(n+1) 2(n+2)
|
||||
* -- - ------ - ------ -
|
||||
* x x x
|
||||
*
|
||||
* Let w = 2n/x and h=2/x, then the above quotient
|
||||
* is equal to the continued fraction:
|
||||
* 1
|
||||
* = -----------------------
|
||||
* 1
|
||||
* w - -----------------
|
||||
* 1
|
||||
* w+h - ---------
|
||||
* w+2h - ...
|
||||
*
|
||||
* To determine how many terms needed, let
|
||||
* Q(0) = w, Q(1) = w(w+h) - 1,
|
||||
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
|
||||
* When Q(k) > 1e4 good for single
|
||||
* When Q(k) > 1e9 good for double
|
||||
* When Q(k) > 1e17 good for quadruple
|
||||
*/
|
||||
/* determine k */
|
||||
let mut t: f32;
|
||||
let mut q0: f32;
|
||||
let mut q1: f32;
|
||||
let mut w: f32;
|
||||
let h: f32;
|
||||
let mut z: f32;
|
||||
let mut tmp: f32;
|
||||
let nf: f32;
|
||||
let mut k: i32;
|
||||
|
||||
nf = (nm1 as f32) + 1.0;
|
||||
w = 2.0 * (nf as f32) / x;
|
||||
h = 2.0 / x;
|
||||
z = w + h;
|
||||
q0 = w;
|
||||
q1 = w * z - 1.0;
|
||||
k = 1;
|
||||
while q1 < 1.0e4 {
|
||||
k += 1;
|
||||
z += h;
|
||||
tmp = z * q1 - q0;
|
||||
q0 = q1;
|
||||
q1 = tmp;
|
||||
}
|
||||
t = 0.0;
|
||||
i = k;
|
||||
while i >= 0 {
|
||||
t = 1.0 / (2.0 * ((i as f32) + nf) / x - t);
|
||||
nf = (nm1 as f32) + 1.0;
|
||||
w = 2.0 * nf / x;
|
||||
h = 2.0 / x;
|
||||
z = w + h;
|
||||
q0 = w;
|
||||
q1 = w * z - 1.0;
|
||||
k = 1;
|
||||
while q1 < 1.0e4 {
|
||||
k += 1;
|
||||
z += h;
|
||||
tmp = z * q1 - q0;
|
||||
q0 = q1;
|
||||
q1 = tmp;
|
||||
}
|
||||
t = 0.0;
|
||||
i = k;
|
||||
while i >= 0 {
|
||||
t = 1.0 / (2.0 * ((i as f32) + nf) / x - t);
|
||||
i -= 1;
|
||||
}
|
||||
a = t;
|
||||
b = 1.0;
|
||||
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
|
||||
* Hence, if n*(log(2n/x)) > ...
|
||||
* single 8.8722839355e+01
|
||||
* double 7.09782712893383973096e+02
|
||||
* long double 1.1356523406294143949491931077970765006170e+04
|
||||
* then recurrent value may overflow and the result is
|
||||
* likely underflow to zero
|
||||
*/
|
||||
tmp = nf * logf(fabsf(w));
|
||||
if tmp < 88.721679688 {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = 2.0 * (i as f32) * b / x - a;
|
||||
a = temp;
|
||||
i -= 1;
|
||||
}
|
||||
a = t;
|
||||
b = 1.0;
|
||||
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
|
||||
* Hence, if n*(log(2n/x)) > ...
|
||||
* single 8.8722839355e+01
|
||||
* double 7.09782712893383973096e+02
|
||||
* long double 1.1356523406294143949491931077970765006170e+04
|
||||
* then recurrent value may overflow and the result is
|
||||
* likely underflow to zero
|
||||
*/
|
||||
tmp = nf * logf(fabsf(w));
|
||||
if tmp < 88.721679688 {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = 2.0 * (i as f32) * b / x - a;
|
||||
a = temp;
|
||||
i -= 1;
|
||||
}
|
||||
} else {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = 2.0 * (i as f32) * b / x - a;
|
||||
a = temp;
|
||||
/* scale b to avoid spurious overflow */
|
||||
let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
|
||||
if b > x1p60 {
|
||||
a /= b;
|
||||
t /= b;
|
||||
b = 1.0;
|
||||
}
|
||||
i -= 1;
|
||||
} else {
|
||||
i = nm1;
|
||||
while i > 0 {
|
||||
temp = b;
|
||||
b = 2.0 * (i as f32) * b / x - a;
|
||||
a = temp;
|
||||
/* scale b to avoid spurious overflow */
|
||||
let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
|
||||
if b > x1p60 {
|
||||
a /= b;
|
||||
t /= b;
|
||||
b = 1.0;
|
||||
}
|
||||
i -= 1;
|
||||
}
|
||||
z = j0f(x);
|
||||
w = j1f(x);
|
||||
if fabsf(z) >= fabsf(w) {
|
||||
b = t * z / b;
|
||||
} else {
|
||||
b = t * w / a;
|
||||
}
|
||||
}
|
||||
z = j0f(x);
|
||||
w = j1f(x);
|
||||
if fabsf(z) >= fabsf(w) {
|
||||
b = t * z / b;
|
||||
} else {
|
||||
b = t * w / a;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -160,7 +160,8 @@ fn sin_pi(mut x: f64) -> f64 {
|
|||
1 => k_cos(x, 0.0),
|
||||
2 => k_sin(-x, 0.0, 0),
|
||||
3 => -k_cos(x, 0.0),
|
||||
0 | _ => k_sin(x, 0.0, 0),
|
||||
// 0
|
||||
_ => k_sin(x, 0.0, 0),
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -95,7 +95,8 @@ fn sin_pi(mut x: f32) -> f32 {
|
|||
1 => k_cosf(y),
|
||||
2 => k_sinf(-y),
|
||||
3 => -k_cosf(y),
|
||||
0 | _ => k_sinf(y),
|
||||
// 0
|
||||
_ => k_sinf(y),
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -10,8 +10,8 @@ pub fn nextafter(x: f64, y: f64) -> f64 {
|
|||
return y;
|
||||
}
|
||||
|
||||
let ax = ux_i & !1_u64 / 2;
|
||||
let ay = uy_i & !1_u64 / 2;
|
||||
let ax = ux_i & (!1_u64 / 2);
|
||||
let ay = uy_i & (!1_u64 / 2);
|
||||
if ax == 0 {
|
||||
if ay == 0 {
|
||||
return y;
|
||||
|
|
|
|||
|
|
@ -98,8 +98,8 @@ pub fn pow(x: f64, y: f64) -> f64 {
|
|||
let (hx, lx): (i32, u32) = ((x.to_bits() >> 32) as i32, x.to_bits() as u32);
|
||||
let (hy, ly): (i32, u32) = ((y.to_bits() >> 32) as i32, y.to_bits() as u32);
|
||||
|
||||
let mut ix: i32 = (hx & 0x7fffffff) as i32;
|
||||
let iy: i32 = (hy & 0x7fffffff) as i32;
|
||||
let mut ix: i32 = hx & 0x7fffffff_i32;
|
||||
let iy: i32 = hy & 0x7fffffff_i32;
|
||||
|
||||
/* x**0 = 1, even if x is NaN */
|
||||
if ((iy as u32) | ly) == 0 {
|
||||
|
|
@ -355,7 +355,7 @@ pub fn pow(x: f64, y: f64) -> f64 {
|
|||
}
|
||||
|
||||
/* compute 2**(p_h+p_l) */
|
||||
let i: i32 = j & (0x7fffffff as i32);
|
||||
let i: i32 = j & 0x7fffffff_i32;
|
||||
k = (i >> 20) - 0x3ff;
|
||||
let mut n: i32 = 0;
|
||||
|
||||
|
|
|
|||
|
|
@ -13,6 +13,8 @@
|
|||
* ====================================================
|
||||
*/
|
||||
|
||||
use core::cmp::Ordering;
|
||||
|
||||
use super::{fabsf, scalbnf, sqrtf};
|
||||
|
||||
const BP: [f32; 2] = [1.0, 1.5];
|
||||
|
|
@ -115,15 +117,13 @@ pub fn powf(x: f32, y: f32) -> f32 {
|
|||
/* special value of y */
|
||||
if iy == 0x7f800000 {
|
||||
/* y is +-inf */
|
||||
if ix == 0x3f800000 {
|
||||
match ix.cmp(&0x3f800000) {
|
||||
/* (-1)**+-inf is 1 */
|
||||
return 1.0;
|
||||
} else if ix > 0x3f800000 {
|
||||
Ordering::Equal => return 1.0,
|
||||
/* (|x|>1)**+-inf = inf,0 */
|
||||
return if hy >= 0 { y } else { 0.0 };
|
||||
} else {
|
||||
Ordering::Greater => return if hy >= 0 { y } else { 0.0 },
|
||||
/* (|x|<1)**+-inf = 0,inf */
|
||||
return if hy >= 0 { 0.0 } else { -y };
|
||||
Ordering::Less => return if hy >= 0 { 0.0 } else { -y },
|
||||
}
|
||||
}
|
||||
if iy == 0x3f800000 {
|
||||
|
|
|
|||
|
|
@ -50,7 +50,7 @@ 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 tmp = x as f64 * INV_PIO2 + TO_INT;
|
||||
let tmp = x * 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")))]
|
||||
|
|
|
|||
|
|
@ -425,8 +425,6 @@ pub(crate) fn rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) ->
|
|||
for i in (0..=jz).rev() {
|
||||
fw += i!(fq, i);
|
||||
}
|
||||
// TODO: drop excess precision here once double_t is used
|
||||
fw = fw as f64;
|
||||
i!(y, 0, =, if ih == 0 { fw } else { -fw });
|
||||
fw = i!(fq, 0) - fw;
|
||||
for i in 1..=jz {
|
||||
|
|
|
|||
|
|
@ -67,14 +67,12 @@ pub fn sincosf(x: f32) -> (f32, f32) {
|
|||
}
|
||||
}
|
||||
/* -sin(x+c) is not correct if x+c could be 0: -0 vs +0 */
|
||||
else {
|
||||
if sign {
|
||||
s = -k_sinf(x as f64 + S2PIO2);
|
||||
c = -k_cosf(x as f64 + S2PIO2);
|
||||
} else {
|
||||
s = -k_sinf(x as f64 - S2PIO2);
|
||||
c = -k_cosf(x as f64 - S2PIO2);
|
||||
}
|
||||
else if sign {
|
||||
s = -k_sinf(x as f64 + S2PIO2);
|
||||
c = -k_cosf(x as f64 + S2PIO2);
|
||||
} else {
|
||||
s = -k_sinf(x as f64 - S2PIO2);
|
||||
c = -k_cosf(x as f64 - S2PIO2);
|
||||
}
|
||||
|
||||
return (s, c);
|
||||
|
|
@ -91,14 +89,12 @@ pub fn sincosf(x: f32) -> (f32, f32) {
|
|||
s = -k_cosf(x as f64 - S3PIO2);
|
||||
c = k_sinf(x as f64 - S3PIO2);
|
||||
}
|
||||
} else if sign {
|
||||
s = k_sinf(x as f64 + S4PIO2);
|
||||
c = k_cosf(x as f64 + S4PIO2);
|
||||
} else {
|
||||
if sign {
|
||||
s = k_sinf(x as f64 + S4PIO2);
|
||||
c = k_cosf(x as f64 + S4PIO2);
|
||||
} else {
|
||||
s = k_sinf(x as f64 - S4PIO2);
|
||||
c = k_cosf(x as f64 - S4PIO2);
|
||||
}
|
||||
s = k_sinf(x as f64 - S4PIO2);
|
||||
c = k_cosf(x as f64 - S4PIO2);
|
||||
}
|
||||
|
||||
return (s, c);
|
||||
|
|
|
|||
|
|
@ -144,13 +144,15 @@ pub fn sqrt(x: f64) -> f64 {
|
|||
ix0 = (ix0 & 0x000fffff) | 0x00100000;
|
||||
if (m & 1) == 1 {
|
||||
/* odd m, double x to make it even */
|
||||
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
|
||||
ix0 *= 2;
|
||||
ix0 += ((ix1 & sign) >> 31).0 as i32;
|
||||
ix1 += ix1;
|
||||
}
|
||||
m >>= 1; /* m = [m/2] */
|
||||
|
||||
/* generate sqrt(x) bit by bit */
|
||||
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
|
||||
ix0 *= 2;
|
||||
ix0 += ((ix1 & sign) >> 31).0 as i32;
|
||||
ix1 += ix1;
|
||||
q = 0; /* [q,q1] = sqrt(x) */
|
||||
q1 = Wrapping(0);
|
||||
|
|
@ -165,7 +167,8 @@ pub fn sqrt(x: f64) -> f64 {
|
|||
ix0 -= t;
|
||||
q += r.0 as i32;
|
||||
}
|
||||
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
|
||||
ix0 *= 2;
|
||||
ix0 += ((ix1 & sign) >> 31).0 as i32;
|
||||
ix1 += ix1;
|
||||
r >>= 1;
|
||||
}
|
||||
|
|
@ -186,7 +189,8 @@ pub fn sqrt(x: f64) -> f64 {
|
|||
ix1 -= t1;
|
||||
q1 += r;
|
||||
}
|
||||
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
|
||||
ix0 *= 2;
|
||||
ix0 += ((ix1 & sign) >> 31).0 as i32;
|
||||
ix1 += ix1;
|
||||
r >>= 1;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -136,6 +136,9 @@ macro_rules! int_impl_common {
|
|||
}
|
||||
|
||||
fn ilog2(self) -> u32 {
|
||||
// On our older MSRV, this resolves to the trait method. Which won't actually work,
|
||||
// but this is only called behind other gates.
|
||||
#[allow(clippy::incompatible_msrv)]
|
||||
<Self>::ilog2(self)
|
||||
}
|
||||
};
|
||||
|
|
|
|||
|
|
@ -45,7 +45,8 @@ fn sinpi(mut x: f64) -> f64 {
|
|||
1 => k_cos(x, 0.0),
|
||||
2 => k_sin(-x, 0.0, 0),
|
||||
3 => -k_cos(x, 0.0),
|
||||
0 | _ => k_sin(x, 0.0, 0),
|
||||
// 0
|
||||
_ => k_sin(x, 0.0, 0),
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -143,7 +144,7 @@ pub fn tgamma(mut x: f64) -> f64 {
|
|||
/* special cases */
|
||||
if ix >= 0x7ff00000 {
|
||||
/* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
|
||||
return x + core::f64::INFINITY;
|
||||
return x + f64::INFINITY;
|
||||
}
|
||||
if ix < ((0x3ff - 54) << 20) {
|
||||
/* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue