From c2402e4d3c2708cba97e463f7a863d08af8d51a1 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Fri, 1 Nov 2024 06:05:42 -0500 Subject: [PATCH] Fix errors reported by Clippy in `libm` --- library/compiler-builtins/libm/src/lib.rs | 3 + .../compiler-builtins/libm/src/math/asin.rs | 2 +- .../compiler-builtins/libm/src/math/asinf.rs | 2 +- .../compiler-builtins/libm/src/math/atan2f.rs | 22 +- .../compiler-builtins/libm/src/math/atanhf.rs | 2 +- .../compiler-builtins/libm/src/math/exp2f.rs | 4 +- .../compiler-builtins/libm/src/math/expm1.rs | 2 +- .../compiler-builtins/libm/src/math/expm1f.rs | 2 +- .../compiler-builtins/libm/src/math/fabs.rs | 2 - .../compiler-builtins/libm/src/math/fdim.rs | 6 +- .../compiler-builtins/libm/src/math/fdimf.rs | 6 +- .../compiler-builtins/libm/src/math/fmaf.rs | 2 +- .../compiler-builtins/libm/src/math/fmod.rs | 4 +- .../compiler-builtins/libm/src/math/fmodf.rs | 2 +- .../compiler-builtins/libm/src/math/ilogb.rs | 2 +- .../compiler-builtins/libm/src/math/ilogbf.rs | 2 +- library/compiler-builtins/libm/src/math/jn.rs | 240 +++++++++--------- .../compiler-builtins/libm/src/math/jnf.rs | 230 +++++++++-------- .../libm/src/math/lgamma_r.rs | 3 +- .../libm/src/math/lgammaf_r.rs | 3 +- .../libm/src/math/nextafter.rs | 4 +- .../compiler-builtins/libm/src/math/pow.rs | 6 +- .../compiler-builtins/libm/src/math/powf.rs | 12 +- .../libm/src/math/rem_pio2.rs | 2 +- .../libm/src/math/rem_pio2_large.rs | 2 - .../libm/src/math/sincosf.rs | 26 +- .../compiler-builtins/libm/src/math/sqrt.rs | 12 +- .../libm/src/math/support/int_traits.rs | 3 + .../compiler-builtins/libm/src/math/tgamma.rs | 5 +- 29 files changed, 307 insertions(+), 306 deletions(-) diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 98ac55988d0a..511ab598dbad 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -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; diff --git a/library/compiler-builtins/libm/src/math/asin.rs b/library/compiler-builtins/libm/src/math/asin.rs index 12fe08fc7452..12d0cd35fa58 100644 --- a/library/compiler-builtins/libm/src/math/asin.rs +++ b/library/compiler-builtins/libm/src/math/asin.rs @@ -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); diff --git a/library/compiler-builtins/libm/src/math/asinf.rs b/library/compiler-builtins/libm/src/math/asinf.rs index 2c785abe2baf..0ea49c0767cd 100644 --- a/library/compiler-builtins/libm/src/math/asinf.rs +++ b/library/compiler-builtins/libm/src/math/asinf.rs @@ -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); diff --git a/library/compiler-builtins/libm/src/math/atan2f.rs b/library/compiler-builtins/libm/src/math/atan2f.rs index fa33f54f6f76..95b466fff4e4 100644 --- a/library/compiler-builtins/libm/src/math/atan2f.rs +++ b/library/compiler-builtins/libm/src/math/atan2f.rs @@ -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) */ } }; } diff --git a/library/compiler-builtins/libm/src/math/atanhf.rs b/library/compiler-builtins/libm/src/math/atanhf.rs index 3545411bbd5b..80ccec1f67fe 100644 --- a/library/compiler-builtins/libm/src/math/atanhf.rs +++ b/library/compiler-builtins/libm/src/math/atanhf.rs @@ -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 */ diff --git a/library/compiler-builtins/libm/src/math/exp2f.rs b/library/compiler-builtins/libm/src/math/exp2f.rs index f4867b80eedb..f452b6a20f80 100644 --- a/library/compiler-builtins/libm/src/math/exp2f.rs +++ b/library/compiler-builtins/libm/src/math/exp2f.rs @@ -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 */ diff --git a/library/compiler-builtins/libm/src/math/expm1.rs b/library/compiler-builtins/libm/src/math/expm1.rs index 42608509a402..f25153f32a34 100644 --- a/library/compiler-builtins/libm/src/math/expm1.rs +++ b/library/compiler-builtins/libm/src/math/expm1.rs @@ -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 { diff --git a/library/compiler-builtins/libm/src/math/expm1f.rs b/library/compiler-builtins/libm/src/math/expm1f.rs index a862fe2558c4..12c6f532b96a 100644 --- a/library/compiler-builtins/libm/src/math/expm1f.rs +++ b/library/compiler-builtins/libm/src/math/expm1f.rs @@ -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 { diff --git a/library/compiler-builtins/libm/src/math/fabs.rs b/library/compiler-builtins/libm/src/math/fabs.rs index 8d3ea2fd6479..d7980eb65f2a 100644 --- a/library/compiler-builtins/libm/src/math/fabs.rs +++ b/library/compiler-builtins/libm/src/math/fabs.rs @@ -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`. diff --git a/library/compiler-builtins/libm/src/math/fdim.rs b/library/compiler-builtins/libm/src/math/fdim.rs index 014930097a02..7c58cb5a9a27 100644 --- a/library/compiler-builtins/libm/src/math/fdim.rs +++ b/library/compiler-builtins/libm/src/math/fdim.rs @@ -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)] diff --git a/library/compiler-builtins/libm/src/math/fdimf.rs b/library/compiler-builtins/libm/src/math/fdimf.rs index ea0b592d7ab9..2abd49a64c9f 100644 --- a/library/compiler-builtins/libm/src/math/fdimf.rs +++ b/library/compiler-builtins/libm/src/math/fdimf.rs @@ -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)] diff --git a/library/compiler-builtins/libm/src/math/fmaf.rs b/library/compiler-builtins/libm/src/math/fmaf.rs index 10bdaeab33de..79371c836c8f 100644 --- a/library/compiler-builtins/libm/src/math/fmaf.rs +++ b/library/compiler-builtins/libm/src/math/fmaf.rs @@ -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) }; diff --git a/library/compiler-builtins/libm/src/math/fmod.rs b/library/compiler-builtins/libm/src/math/fmod.rs index d892ffd8b72a..df16162bcae8 100644 --- a/library/compiler-builtins/libm/src/math/fmod.rs +++ b/library/compiler-builtins/libm/src/math/fmod.rs @@ -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) } diff --git a/library/compiler-builtins/libm/src/math/fmodf.rs b/library/compiler-builtins/libm/src/math/fmodf.rs index 1d80013842e7..671af85800f1 100644 --- a/library/compiler-builtins/libm/src/math/fmodf.rs +++ b/library/compiler-builtins/libm/src/math/fmodf.rs @@ -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 { diff --git a/library/compiler-builtins/libm/src/math/ilogb.rs b/library/compiler-builtins/libm/src/math/ilogb.rs index 9d58d06081b7..ccc4914be2bf 100644 --- a/library/compiler-builtins/libm/src/math/ilogb.rs +++ b/library/compiler-builtins/libm/src/math/ilogb.rs @@ -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 } diff --git a/library/compiler-builtins/libm/src/math/ilogbf.rs b/library/compiler-builtins/libm/src/math/ilogbf.rs index 85deb43c83b7..3585d6d36f16 100644 --- a/library/compiler-builtins/libm/src/math/ilogbf.rs +++ b/library/compiler-builtins/libm/src/math/ilogbf.rs @@ -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 } diff --git a/library/compiler-builtins/libm/src/math/jn.rs b/library/compiler-builtins/libm/src/math/jn.rs index aff051f24061..7f98ddc055a6 100644 --- a/library/compiler-builtins/libm/src/math/jn.rs +++ b/library/compiler-builtins/libm/src/math/jn.rs @@ -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 { diff --git a/library/compiler-builtins/libm/src/math/jnf.rs b/library/compiler-builtins/libm/src/math/jnf.rs index e5afda44896d..754f8f33b5af 100644 --- a/library/compiler-builtins/libm/src/math/jnf.rs +++ b/library/compiler-builtins/libm/src/math/jnf.rs @@ -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; } } diff --git a/library/compiler-builtins/libm/src/math/lgamma_r.rs b/library/compiler-builtins/libm/src/math/lgamma_r.rs index b26177e6ebf7..6becaad2ce91 100644 --- a/library/compiler-builtins/libm/src/math/lgamma_r.rs +++ b/library/compiler-builtins/libm/src/math/lgamma_r.rs @@ -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), } } diff --git a/library/compiler-builtins/libm/src/math/lgammaf_r.rs b/library/compiler-builtins/libm/src/math/lgammaf_r.rs index 723c90daf1ed..10cecee541cc 100644 --- a/library/compiler-builtins/libm/src/math/lgammaf_r.rs +++ b/library/compiler-builtins/libm/src/math/lgammaf_r.rs @@ -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), } } diff --git a/library/compiler-builtins/libm/src/math/nextafter.rs b/library/compiler-builtins/libm/src/math/nextafter.rs index 05762619109a..422bd7496d5e 100644 --- a/library/compiler-builtins/libm/src/math/nextafter.rs +++ b/library/compiler-builtins/libm/src/math/nextafter.rs @@ -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; diff --git a/library/compiler-builtins/libm/src/math/pow.rs b/library/compiler-builtins/libm/src/math/pow.rs index 7ecad291d182..736465cd16f6 100644 --- a/library/compiler-builtins/libm/src/math/pow.rs +++ b/library/compiler-builtins/libm/src/math/pow.rs @@ -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; diff --git a/library/compiler-builtins/libm/src/math/powf.rs b/library/compiler-builtins/libm/src/math/powf.rs index 2d9d1e4bbfdd..839c6c23d436 100644 --- a/library/compiler-builtins/libm/src/math/powf.rs +++ b/library/compiler-builtins/libm/src/math/powf.rs @@ -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 { diff --git a/library/compiler-builtins/libm/src/math/rem_pio2.rs b/library/compiler-builtins/libm/src/math/rem_pio2.rs index 4dfb8c658837..917e90819a50 100644 --- a/library/compiler-builtins/libm/src/math/rem_pio2.rs +++ b/library/compiler-builtins/libm/src/math/rem_pio2.rs @@ -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")))] diff --git a/library/compiler-builtins/libm/src/math/rem_pio2_large.rs b/library/compiler-builtins/libm/src/math/rem_pio2_large.rs index 1dfbba3b1b22..ec8397f4b6fc 100644 --- a/library/compiler-builtins/libm/src/math/rem_pio2_large.rs +++ b/library/compiler-builtins/libm/src/math/rem_pio2_large.rs @@ -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 { diff --git a/library/compiler-builtins/libm/src/math/sincosf.rs b/library/compiler-builtins/libm/src/math/sincosf.rs index 423845e44c2a..f3360767683e 100644 --- a/library/compiler-builtins/libm/src/math/sincosf.rs +++ b/library/compiler-builtins/libm/src/math/sincosf.rs @@ -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); diff --git a/library/compiler-builtins/libm/src/math/sqrt.rs b/library/compiler-builtins/libm/src/math/sqrt.rs index a443b7e4c613..3eaf52cda8ea 100644 --- a/library/compiler-builtins/libm/src/math/sqrt.rs +++ b/library/compiler-builtins/libm/src/math/sqrt.rs @@ -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; } diff --git a/library/compiler-builtins/libm/src/math/support/int_traits.rs b/library/compiler-builtins/libm/src/math/support/int_traits.rs index c5feef8d7268..b08907aa5b87 100644 --- a/library/compiler-builtins/libm/src/math/support/int_traits.rs +++ b/library/compiler-builtins/libm/src/math/support/int_traits.rs @@ -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)] ::ilog2(self) } }; diff --git a/library/compiler-builtins/libm/src/math/tgamma.rs b/library/compiler-builtins/libm/src/math/tgamma.rs index 3f38c0b1d9fd..60451416ab32 100644 --- a/library/compiler-builtins/libm/src/math/tgamma.rs +++ b/library/compiler-builtins/libm/src/math/tgamma.rs @@ -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 */