diff options
63 files changed, 560 insertions, 189 deletions
diff --git a/.cargo_vcs_info.json b/.cargo_vcs_info.json index f2062d3..9ccba83 100644 --- a/.cargo_vcs_info.json +++ b/.cargo_vcs_info.json @@ -1,5 +1,6 @@ { "git": { - "sha1": "3d729b7a85daad3b4d54b22e02d00681762bc602" - } -} + "sha1": "4c8a973741c014b11ce7f1477693a3e5d4ef9609" + }, + "path_in_vcs": "" +}
\ No newline at end of file diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 80ce4eb..decd71f 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -12,7 +12,7 @@ jobs: - arm-unknown-linux-gnueabi - arm-unknown-linux-gnueabihf - armv7-unknown-linux-gnueabihf - - i686-unknown-linux-gnu + # - i686-unknown-linux-gnu - mips-unknown-linux-gnu - mips64-unknown-linux-gnuabi64 - mips64el-unknown-linux-gnuabi64 @@ -43,7 +43,7 @@ rust_library { host_supported: true, crate_name: "libm", cargo_env_compat: true, - cargo_pkg_version: "0.2.1", + cargo_pkg_version: "0.2.6", srcs: ["src/lib.rs"], edition: "2018", features: ["default"], @@ -51,6 +51,8 @@ rust_library { "//apex_available:platform", "com.android.resolv", ], + product_available: true, + vendor_available: true, min_sdk_version: "29", } @@ -59,7 +61,7 @@ rust_test { host_supported: true, crate_name: "libm", cargo_env_compat: true, - cargo_pkg_version: "0.2.1", + cargo_pkg_version: "0.2.6", srcs: ["src/lib.rs"], test_suites: ["general-tests"], auto_gen_config: true, diff --git a/CHANGELOG.md b/CHANGELOG.md index 28e2705..e8e9acf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,30 @@ This project adheres to [Semantic Versioning](http://semver.org/). ... +## [v0.2.1] - 2019-11-22 + +### Fixed +- sincosf + +## [v0.2.0] - 2019-10-18 + +### Added +- Benchmarks +- signum +- remainder +- remainderf +- nextafter +- nextafterf + +### Fixed +- Rounding to negative zero +- Overflows in rem_pio2 and remquo +- Overflows in fma +- sincosf + +### Removed +- F32Ext and F64Ext traits + ## [v0.1.4] - 2019-06-12 ### Fixed @@ -90,7 +114,9 @@ This project adheres to [Semantic Versioning](http://semver.org/). - Initial release -[Unreleased]: https://github.com/japaric/libm/compare/v0.1.4...HEAD +[Unreleased]: https://github.com/japaric/libm/compare/v0.2.1...HEAD +[v0.2.1]: https://github.com/japaric/libm/compare/0.2.0...v0.2.1 +[v0.2.0]: https://github.com/japaric/libm/compare/0.1.4...v0.2.0 [v0.1.4]: https://github.com/japaric/libm/compare/0.1.3...v0.1.4 [v0.1.3]: https://github.com/japaric/libm/compare/v0.1.2...0.1.3 [v0.1.2]: https://github.com/japaric/libm/compare/v0.1.1...v0.1.2 @@ -3,26 +3,34 @@ # When uploading crates to the registry Cargo will automatically # "normalize" Cargo.toml files for maximal compatibility # with all versions of Cargo and also rewrite `path` dependencies -# to registry (e.g., crates.io) dependencies +# to registry (e.g., crates.io) dependencies. # -# If you believe there's an error in this file please file an -# issue against the rust-lang/cargo repository. If you're -# editing this file be aware that the upstream Cargo.toml -# will likely look very different (and much more reasonable) +# If you are reading this file be aware that the original Cargo.toml +# will likely look very different (and much more reasonable). +# See Cargo.toml.orig for the original contents. [package] edition = "2018" name = "libm" -version = "0.2.1" +version = "0.2.6" authors = ["Jorge Aparicio <jorge@japaric.io>"] description = "libm in pure Rust" documentation = "https://docs.rs/libm" -keywords = ["libm", "math"] +readme = "README.md" +keywords = [ + "libm", + "math", +] categories = ["no-std"] license = "MIT OR Apache-2.0" repository = "https://github.com/rust-lang/libm" + +[profile.release] +lto = "fat" + [dev-dependencies.no-panic] version = "0.1.8" + [build-dependencies.rand] version = "0.6.5" optional = true diff --git a/Cargo.toml.orig b/Cargo.toml.orig index d9d6680..f942fde 100644 --- a/Cargo.toml.orig +++ b/Cargo.toml.orig @@ -6,8 +6,9 @@ documentation = "https://docs.rs/libm" keywords = ["libm", "math"] license = "MIT OR Apache-2.0" name = "libm" +readme = "README.md" repository = "https://github.com/rust-lang/libm" -version = "0.2.1" +version = "0.2.6" edition = "2018" [features] @@ -32,3 +33,7 @@ no-panic = "0.1.8" [build-dependencies] rand = { version = "0.6.5", optional = true } + +# This is needed for no-panic to correctly detect the lack of panics +[profile.release] +lto = "fat" @@ -1,3 +1,7 @@ +# This project was upgraded with external_updater. +# Usage: tools/external_updater/updater.sh update rust/crates/libm +# For more info, check https://cs.android.com/android/platform/superproject/+/master:tools/external_updater/README.md + name: "libm" description: "libm in pure Rust" third_party { @@ -7,13 +11,13 @@ third_party { } url { type: ARCHIVE - value: "https://static.crates.io/crates/libm/libm-0.2.1.crate" + value: "https://static.crates.io/crates/libm/libm-0.2.6.crate" } - version: "0.2.1" + version: "0.2.6" license_type: NOTICE last_upgrade_date { - year: 2020 - month: 11 - day: 26 + year: 2022 + month: 12 + day: 12 } } @@ -2,7 +2,7 @@ A port of [MUSL]'s libm to Rust. -[MUSL]: https://www.musl-libc.org/ +[MUSL]: https://musl.libc.org/ ## Goals diff --git a/TEST_MAPPING b/TEST_MAPPING index bfdf123..45401cb 100644 --- a/TEST_MAPPING +++ b/TEST_MAPPING @@ -3,21 +3,18 @@ "imports": [ { "path": "external/rust/crates/quiche" + }, + { + "path": "packages/modules/DnsResolver" } ], "presubmit": [ { - "name": "doh_unit_test" - }, - { "name": "libm_test_src_lib" } ], "presubmit-rust": [ { - "name": "doh_unit_test" - }, - { "name": "libm_test_src_lib" } ] @@ -18,6 +18,7 @@ fn main() { mod musl_reference_tests { use rand::seq::SliceRandom; use rand::Rng; + use std::env; use std::fs; use std::process::Command; @@ -26,7 +27,19 @@ mod musl_reference_tests { // These files are all internal functions or otherwise miscellaneous, not // defining a function we want to test. - const IGNORED_FILES: &[&str] = &["fenv.rs"]; + const IGNORED_FILES: &[&str] = &[ + "fenv.rs", + // These are giving slightly different results compared to musl + "lgamma.rs", + "lgammaf.rs", + "tgamma.rs", + "j0.rs", + "j0f.rs", + "jn.rs", + "jnf.rs", + "j1.rs", + "j1f.rs", + ]; struct Function { name: String, @@ -48,6 +61,12 @@ mod musl_reference_tests { } pub fn generate() { + // PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 + let target_arch = env::var("CARGO_CFG_TARGET_ARCH").unwrap(); + if target_arch == "powerpc64" { + return; + } + let files = fs::read_dir("src/math") .unwrap() .map(|f| f.unwrap().path()) diff --git a/ci/run-docker.sh b/ci/run-docker.sh index e7b80c7..c7ad60f 100755 --- a/ci/run-docker.sh +++ b/ci/run-docker.sh @@ -18,7 +18,7 @@ run() { --user $(id -u):$(id -g) \ -e CARGO_HOME=/cargo \ -e CARGO_TARGET_DIR=/target \ - -v $(dirname $(dirname `which cargo`)):/cargo \ + -v "${HOME}/.cargo":/cargo \ -v `pwd`/target:/target \ -v `pwd`:/checkout:ro \ -v `rustc --print sysroot`:/rust:ro \ @@ -5,6 +5,9 @@ TARGET=$1 CMD="cargo test --all --target $TARGET" +# Needed for no-panic to correct detect a lack of panics +export RUSTFLAGS="$RUSTFLAGS -Ccodegen-units=1" + # stable by default $CMD $CMD --release diff --git a/patches/test_build_fix.patch b/patches/test_build_fix.patch index 82c3756..06ef1e2 100644 --- a/patches/test_build_fix.patch +++ b/patches/test_build_fix.patch @@ -1,17 +1,8 @@ -From c18c704856ffa8b3349401300e66796d4cc8d4f5 Mon Sep 17 00:00:00 2001 -From: Jethro Beekman <jethro@fortanix.com> -Date: Thu, 24 Jun 2021 15:58:36 +0200 -Subject: [PATCH] Fix build failure with latest nightly - ---- - src/math/pow.rs | 4 ++-- - 1 file changed, 2 insertions(+), 2 deletions(-) - diff --git a/src/math/pow.rs b/src/math/pow.rs -index c7fd0df..f79680a 100644 +index 3020b1a..6a19ae6 100644 --- a/src/math/pow.rs +++ b/src/math/pow.rs -@@ -604,7 +604,7 @@ mod tests { +@@ -608,7 +608,7 @@ mod tests { // Factoring -1 out: // (negative anything ^ integer should be (-1 ^ integer) * (positive anything ^ integer)) @@ -20,7 +11,7 @@ index c7fd0df..f79680a 100644 .iter() .for_each(|int_set| { int_set.iter().for_each(|int| { -@@ -616,7 +616,7 @@ mod tests { +@@ -620,7 +620,7 @@ mod tests { // Negative base (imaginary results): // (-anything except 0 and Infinity ^ non-integer should be NAN) @@ -1,10 +1,7 @@ //! libm in pure Rust #![deny(warnings)] #![no_std] -#![cfg_attr( - all(target_arch = "wasm32", feature = "unstable"), - feature(core_intrinsics) -)] +#![cfg_attr(all(feature = "unstable"), feature(core_intrinsics))] #![allow(clippy::unreadable_literal)] #![allow(clippy::many_single_char_names)] #![allow(clippy::needless_return)] @@ -54,5 +51,7 @@ pub fn _eq(a: f64, b: f64) -> Result<(), u64> { } } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(all(test, feature = "musl-reference-tests"))] include!(concat!(env!("OUT_DIR"), "/musl-tests.rs")); diff --git a/src/math/acosh.rs b/src/math/acosh.rs index ac7a5f1..d1f5b9f 100644 --- a/src/math/acosh.rs +++ b/src/math/acosh.rs @@ -7,6 +7,7 @@ const LN2: f64 = 0.693147180559945309417232121458176568; /* 0x3fe62e42, 0xfefa3 /// Calculates the inverse hyperbolic cosine of `x`. /// Is defined as `log(x + sqrt(x*x-1))`. /// `x` must be a number greater than or equal to 1. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn acosh(x: f64) -> f64 { let u = x.to_bits(); let e = ((u >> 52) as usize) & 0x7ff; diff --git a/src/math/acoshf.rs b/src/math/acoshf.rs index 0879e1e..ad3455f 100644 --- a/src/math/acoshf.rs +++ b/src/math/acoshf.rs @@ -7,6 +7,7 @@ const LN2: f32 = 0.693147180559945309417232121458176568; /// Calculates the inverse hyperbolic cosine of `x`. /// Is defined as `log(x + sqrt(x*x-1))`. /// `x` must be a number greater than or equal to 1. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn acoshf(x: f32) -> f32 { let u = x.to_bits(); let a = u & 0x7fffffff; diff --git a/src/math/asinh.rs b/src/math/asinh.rs index 1429535..0abd80c 100644 --- a/src/math/asinh.rs +++ b/src/math/asinh.rs @@ -7,6 +7,7 @@ const LN2: f64 = 0.693147180559945309417232121458176568; /* 0x3fe62e42, 0xfefa3 /// /// Calculates the inverse hyperbolic sine of `x`. /// Is defined as `sgn(x)*log(|x|+sqrt(x*x+1))`. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn asinh(mut x: f64) -> f64 { let mut u = x.to_bits(); let e = ((u >> 52) as usize) & 0x7ff; diff --git a/src/math/asinhf.rs b/src/math/asinhf.rs index e22a291..09c7782 100644 --- a/src/math/asinhf.rs +++ b/src/math/asinhf.rs @@ -7,6 +7,7 @@ const LN2: f32 = 0.693147180559945309417232121458176568; /// /// Calculates the inverse hyperbolic sine of `x`. /// Is defined as `sgn(x)*log(|x|+sqrt(x*x+1))`. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn asinhf(mut x: f32) -> f32 { let u = x.to_bits(); let i = u & 0x7fffffff; diff --git a/src/math/atanf.rs b/src/math/atanf.rs index 73f3352..d042b3b 100644 --- a/src/math/atanf.rs +++ b/src/math/atanf.rs @@ -56,7 +56,7 @@ pub fn atanf(mut x: f32) -> f32 { if x.is_nan() { return x; } - z = ATAN_HI[3] + x1p_120; + z = i!(ATAN_HI, 3) + x1p_120; return if sign { -z } else { z }; } let id = if ix < 0x3ee00000 { @@ -97,13 +97,13 @@ pub fn atanf(mut x: f32) -> f32 { z = x * x; let w = z * z; /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ - let s1 = z * (A_T[0] + w * (A_T[2] + w * A_T[4])); - let s2 = w * (A_T[1] + w * A_T[3]); + let s1 = z * (i!(A_T, 0) + w * (i!(A_T, 2) + w * i!(A_T, 4))); + let s2 = w * (i!(A_T, 1) + w * i!(A_T, 3)); if id < 0 { return x - x * (s1 + s2); } let id = id as usize; - let z = ATAN_HI[id] - ((x * (s1 + s2) - ATAN_LO[id]) - x); + let z = i!(ATAN_HI, id) - ((x * (s1 + s2) - i!(ATAN_LO, id)) - x); if sign { -z } else { diff --git a/src/math/atanh.rs b/src/math/atanh.rs index 79a989c..b984c4a 100644 --- a/src/math/atanh.rs +++ b/src/math/atanh.rs @@ -5,6 +5,7 @@ use super::log1p; /// /// Calculates the inverse hyperbolic tangent of `x`. /// Is defined as `log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2`. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn atanh(x: f64) -> f64 { let u = x.to_bits(); let e = ((u >> 52) as usize) & 0x7ff; diff --git a/src/math/atanhf.rs b/src/math/atanhf.rs index 7b2f34d..a1aa314 100644 --- a/src/math/atanhf.rs +++ b/src/math/atanhf.rs @@ -5,6 +5,7 @@ use super::log1pf; /// /// Calculates the inverse hyperbolic tangent of `x`. /// Is defined as `log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2`. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn atanhf(mut x: f32) -> f32 { let mut u = x.to_bits(); let sign = (u >> 31) != 0; diff --git a/src/math/ceil.rs b/src/math/ceil.rs index eda28b9..22d8929 100644 --- a/src/math/ceil.rs +++ b/src/math/ceil.rs @@ -1,3 +1,4 @@ +#![allow(unreachable_code)] use core::f64; const TOINT: f64 = 1. / f64::EPSILON; @@ -15,6 +16,24 @@ pub fn ceil(x: f64) -> f64 { return unsafe { ::core::intrinsics::ceilf64(x) } } } + #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] + { + //use an alternative implementation on x86, because the + //main implementation fails with the x87 FPU used by + //debian i386, probablly due to excess precision issues. + //basic implementation taken from https://github.com/rust-lang/libm/issues/219 + use super::fabs; + if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() { + let truncated = x as i64 as f64; + if truncated < x { + return truncated + 1.0; + } else { + return truncated; + } + } else { + return x; + } + } let u: u64 = x.to_bits(); let e: i64 = (u >> 52 & 0x7ff) as i64; let y: f64; diff --git a/src/math/ceilf.rs b/src/math/ceilf.rs index f1edbd0..7bcc647 100644 --- a/src/math/ceilf.rs +++ b/src/math/ceilf.rs @@ -40,6 +40,8 @@ pub fn ceilf(x: f32) -> f32 { f32::from_bits(ui) } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { use super::*; diff --git a/src/math/copysign.rs b/src/math/copysign.rs index 1527fb6..1f4a35a 100644 --- a/src/math/copysign.rs +++ b/src/math/copysign.rs @@ -2,6 +2,7 @@ /// /// Constructs a number with the magnitude (absolute value) of its /// first argument, `x`, and the sign of its second argument, `y`. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn copysign(x: f64, y: f64) -> f64 { let mut ux = x.to_bits(); let uy = y.to_bits(); diff --git a/src/math/copysignf.rs b/src/math/copysignf.rs index 3514856..6c346e3 100644 --- a/src/math/copysignf.rs +++ b/src/math/copysignf.rs @@ -2,6 +2,7 @@ /// /// Constructs a number with the magnitude (absolute value) of its /// first argument, `x`, and the sign of its second argument, `y`. +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn copysignf(x: f32, y: f32) -> f32 { let mut ux = x.to_bits(); let uy = y.to_bits(); diff --git a/src/math/erf.rs b/src/math/erf.rs index a2c617d..5e21ba5 100644 --- a/src/math/erf.rs +++ b/src/math/erf.rs @@ -219,6 +219,7 @@ fn erfc2(ix: u32, mut x: f64) -> f64 { /// Calculates an approximation to the “error function”, which estimates /// the probability that an observation will fall within x standard /// deviations of the mean (assuming a normal distribution). +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn erf(x: f64) -> f64 { let r: f64; let s: f64; diff --git a/src/math/erff.rs b/src/math/erff.rs index 3840522..f74d4b6 100644 --- a/src/math/erff.rs +++ b/src/math/erff.rs @@ -130,6 +130,7 @@ fn erfc2(mut ix: u32, mut x: f32) -> f32 { /// Calculates an approximation to the “error function”, which estimates /// the probability that an observation will fall within x standard /// deviations of the mean (assuming a normal distribution). +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn erff(x: f32) -> f32 { let r: f32; let s: f32; diff --git a/src/math/exp.rs b/src/math/exp.rs index 5b163f9..d499427 100644 --- a/src/math/exp.rs +++ b/src/math/exp.rs @@ -124,7 +124,7 @@ pub fn exp(mut x: f64) -> f64 { /* if |x| > 0.5 ln2 */ if hx >= 0x3ff0a2b2 { /* if |x| >= 1.5 ln2 */ - k = (INVLN2 * x + HALF[sign as usize]) as i32; + k = (INVLN2 * x + i!(HALF, sign as usize)) as i32; } else { k = 1 - sign - sign; } diff --git a/src/math/exp10.rs b/src/math/exp10.rs index 9537f76..559930e 100644 --- a/src/math/exp10.rs +++ b/src/math/exp10.rs @@ -6,16 +6,17 @@ const P10: &[f64] = &[ 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, ]; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn exp10(x: f64) -> f64 { let (mut y, n) = modf(x); let u: u64 = n.to_bits(); /* fabs(n) < 16 without raising invalid on nan */ if (u >> 52 & 0x7ff) < 0x3ff + 4 { if y == 0.0 { - return P10[((n as isize) + 15) as usize]; + return i!(P10, ((n as isize) + 15) as usize); } y = exp2(LN10 * y); - return y * P10[((n as isize) + 15) as usize]; + return y * i!(P10, ((n as isize) + 15) as usize); } return pow(10.0, x); } diff --git a/src/math/exp10f.rs b/src/math/exp10f.rs index d45fff3..1279bc6 100644 --- a/src/math/exp10f.rs +++ b/src/math/exp10f.rs @@ -6,16 +6,17 @@ const P10: &[f32] = &[ 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, ]; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn exp10f(x: f32) -> f32 { let (mut y, n) = modff(x); let u = n.to_bits(); /* fabsf(n) < 8 without raising invalid on nan */ if (u >> 23 & 0xff) < 0x7f + 3 { if y == 0.0 { - return P10[((n as isize) + 7) as usize]; + return i!(P10, ((n as isize) + 7) as usize); } y = exp2f(LN10_F32 * y); - return y * P10[((n as isize) + 7) as usize]; + return y * i!(P10, ((n as isize) + 7) as usize); } return exp2(LN10_F64 * (x as f64)) as f32; } diff --git a/src/math/exp2.rs b/src/math/exp2.rs index 8ea434d..e0e385d 100644 --- a/src/math/exp2.rs +++ b/src/math/exp2.rs @@ -374,14 +374,14 @@ pub fn exp2(mut x: f64) -> f64 { let mut i0 = ui as u32; i0 = i0.wrapping_add(TBLSIZE as u32 / 2); let ku = i0 / TBLSIZE as u32 * TBLSIZE as u32; - let ki = ku as i32 / TBLSIZE as i32; + let ki = div!(ku as i32, TBLSIZE as i32); i0 %= TBLSIZE as u32; let uf = f64::from_bits(ui) - redux; let mut z = x - uf; /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */ - let t = f64::from_bits(TBL[2 * i0 as usize]); /* exp2t[i0] */ - z -= f64::from_bits(TBL[2 * i0 as usize + 1]); /* eps[i0] */ + let t = f64::from_bits(i!(TBL, 2 * i0 as usize)); /* exp2t[i0] */ + z -= f64::from_bits(i!(TBL, 2 * i0 as usize + 1)); /* eps[i0] */ let r = t + t * z * (p1 + z * (p2 + z * (p3 + z * (p4 + z * p5)))); scalbn(r, ki) diff --git a/src/math/exp2f.rs b/src/math/exp2f.rs index 8a890b8..f4867b8 100644 --- a/src/math/exp2f.rs +++ b/src/math/exp2f.rs @@ -126,7 +126,7 @@ pub fn exp2f(mut x: f32) -> f32 { uf -= redux; let z: f64 = (x - uf) as f64; /* Compute r = exp2(y) = exp2ft[i0] * p(z). */ - let r: f64 = f64::from_bits(EXP2FT[i0 as usize]); + let r: f64 = f64::from_bits(i!(EXP2FT, i0 as usize)); let t: f64 = r as f64 * z; let r: f64 = r + t * (p1 as f64 + z * p2 as f64) + t * (z * z) * (p3 as f64 + z * p4 as f64); diff --git a/src/math/expf.rs b/src/math/expf.rs index 47c1b2c..a53aa90 100644 --- a/src/math/expf.rs +++ b/src/math/expf.rs @@ -70,7 +70,7 @@ pub fn expf(mut x: f32) -> f32 { /* if |x| > 0.5 ln2 */ if hx > 0x3f851592 { /* if |x| > 1.5 ln2 */ - k = (INV_LN2 * x + HALF[sign as usize]) as i32; + k = (INV_LN2 * x + i!(HALF, sign as usize)) as i32; } else { k = 1 - sign - sign; } diff --git a/src/math/fabsf.rs b/src/math/fabsf.rs index 6655c4c..23f3646 100644 --- a/src/math/fabsf.rs +++ b/src/math/fabsf.rs @@ -14,6 +14,8 @@ pub fn fabsf(x: f32) -> f32 { f32::from_bits(x.to_bits() & 0x7fffffff) } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { use super::*; diff --git a/src/math/fenv.rs b/src/math/fenv.rs index 652e603..c91272e 100644 --- a/src/math/fenv.rs +++ b/src/math/fenv.rs @@ -5,7 +5,6 @@ pub(crate) const FE_UNDERFLOW: i32 = 0; pub(crate) const FE_INEXACT: i32 = 0; pub(crate) const FE_TONEAREST: i32 = 0; -pub(crate) const FE_TOWARDZERO: i32 = 0; #[inline] pub(crate) fn feclearexcept(_mask: i32) -> i32 { @@ -26,8 +25,3 @@ pub(crate) fn fetestexcept(_mask: i32) -> i32 { pub(crate) fn fegetround() -> i32 { FE_TONEAREST } - -#[inline] -pub(crate) fn fesetround(_r: i32) -> i32 { - 0 -} diff --git a/src/math/floor.rs b/src/math/floor.rs index b2b7605..d09f9a1 100644 --- a/src/math/floor.rs +++ b/src/math/floor.rs @@ -1,3 +1,4 @@ +#![allow(unreachable_code)] use core::f64; const TOINT: f64 = 1. / f64::EPSILON; @@ -15,6 +16,24 @@ pub fn floor(x: f64) -> f64 { return unsafe { ::core::intrinsics::floorf64(x) } } } + #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] + { + //use an alternative implementation on x86, because the + //main implementation fails with the x87 FPU used by + //debian i386, probablly due to excess precision issues. + //basic implementation taken from https://github.com/rust-lang/libm/issues/219 + use super::fabs; + if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() { + let truncated = x as i64 as f64; + if truncated > x { + return truncated - 1.0; + } else { + return truncated; + } + } else { + return x; + } + } let ui = x.to_bits(); let e = ((ui >> 52) & 0x7ff) as i32; diff --git a/src/math/floorf.rs b/src/math/floorf.rs index 287f086..dfdab91 100644 --- a/src/math/floorf.rs +++ b/src/math/floorf.rs @@ -40,6 +40,8 @@ pub fn floorf(x: f32) -> f32 { f32::from_bits(ui) } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { use super::*; diff --git a/src/math/fma.rs b/src/math/fma.rs index 3219dbc..f9a86dc 100644 --- a/src/math/fma.rs +++ b/src/math/fma.rs @@ -122,12 +122,12 @@ pub fn fma(x: f64, y: f64, z: f64) -> f64 { rhi += zhi + (rlo < zlo) as u64; } else { /* r -= z */ - let t = rlo; - rlo = rlo.wrapping_sub(zlo); - rhi = rhi.wrapping_sub(zhi.wrapping_sub((t < rlo) as u64)); + let (res, borrow) = rlo.overflowing_sub(zlo); + rlo = res; + rhi = rhi.wrapping_sub(zhi.wrapping_add(borrow as u64)); if (rhi >> 63) != 0 { - rlo = (-(rlo as i64)) as u64; - rhi = (-(rhi as i64)) as u64 - (rlo != 0) as u64; + rlo = (rlo as i64).wrapping_neg() as u64; + rhi = (rhi as i64).wrapping_neg() as u64 - (rlo != 0) as u64; sign = (sign == 0) as i32; } nonzero = (rhi != 0) as i32; @@ -218,6 +218,26 @@ mod tests { -0.00000000000000022204460492503126, ); - assert_eq!(fma(-0.992, -0.992, -0.992), -0.00793599999988632,); + let result = fma(-0.992, -0.992, -0.992); + //force rounding to storage format on x87 to prevent superious errors. + #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] + let result = force_eval!(result); + assert_eq!(result, -0.007936000000000007,); + } + + #[test] + fn fma_sbb() { + assert_eq!( + fma(-(1.0 - f64::EPSILON), f64::MIN, f64::MIN), + -3991680619069439e277 + ); + } + + #[test] + fn fma_underflow() { + assert_eq!( + fma(1.1102230246251565e-16, -9.812526705433188e-305, 1.0894e-320), + 0.0, + ); } } diff --git a/src/math/fmaf.rs b/src/math/fmaf.rs index 03d371c..2848f2a 100644 --- a/src/math/fmaf.rs +++ b/src/math/fmaf.rs @@ -29,8 +29,7 @@ use core::f32; use core::ptr::read_volatile; use super::fenv::{ - feclearexcept, fegetround, feraiseexcept, fesetround, fetestexcept, FE_INEXACT, FE_TONEAREST, - FE_TOWARDZERO, FE_UNDERFLOW, + feclearexcept, fegetround, feraiseexcept, fetestexcept, FE_INEXACT, FE_TONEAREST, FE_UNDERFLOW, }; /* @@ -91,16 +90,28 @@ pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 { * If result is inexact, and exactly halfway between two float values, * we need to adjust the low-order bit in the direction of the error. */ - fesetround(FE_TOWARDZERO); - // prevent `vxy + z` from being CSE'd with `xy + z` above - let vxy: f64 = unsafe { read_volatile(&xy) }; - let mut adjusted_result: f64 = vxy + z as f64; - fesetround(FE_TONEAREST); - if result == adjusted_result { - ui = adjusted_result.to_bits(); + let neg = ui >> 63 != 0; + let err = if neg == (z as f64 > xy) { + xy - result + z as f64 + } else { + z as f64 - result + xy + }; + if neg == (err < 0.0) { ui += 1; - adjusted_result = f64::from_bits(ui); + } else { + ui -= 1; + } + f64::from_bits(ui) as f32 +} + +#[cfg(test)] +mod tests { + #[test] + fn issue_263() { + let a = f32::from_bits(1266679807); + let b = f32::from_bits(1300234242); + let c = f32::from_bits(1115553792); + let expected = f32::from_bits(1501560833); + assert_eq!(super::fmaf(a, b, c), expected); } - z = adjusted_result as f32; - z } diff --git a/src/math/ilogb.rs b/src/math/ilogb.rs index 0a380b7..7d74dcf 100644 --- a/src/math/ilogb.rs +++ b/src/math/ilogb.rs @@ -1,6 +1,7 @@ const FP_ILOGBNAN: i32 = -1 - 0x7fffffff; const FP_ILOGB0: i32 = FP_ILOGBNAN; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn ilogb(x: f64) -> i32 { let mut i: u64 = x.to_bits(); let e = ((i >> 52) & 0x7ff) as i32; diff --git a/src/math/ilogbf.rs b/src/math/ilogbf.rs index b384fa4..0fa5874 100644 --- a/src/math/ilogbf.rs +++ b/src/math/ilogbf.rs @@ -1,6 +1,7 @@ const FP_ILOGBNAN: i32 = -1 - 0x7fffffff; const FP_ILOGB0: i32 = FP_ILOGBNAN; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn ilogbf(x: f32) -> i32 { let mut i = x.to_bits(); let e = ((i >> 23) & 0xff) as i32; diff --git a/src/math/j1f.rs b/src/math/j1f.rs index 5095894..c39f8ff 100644 --- a/src/math/j1f.rs +++ b/src/math/j1f.rs @@ -357,6 +357,8 @@ fn qonef(x: f32) -> f32 { return (0.375 + r / s) / x; } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { use super::{j1f, y1f}; @@ -367,6 +369,12 @@ mod tests { } #[test] fn test_y1f_2002() { - assert_eq!(y1f(2.0000002_f32), -0.10703229_f32); + //allow slightly different result on x87 + let res = y1f(2.0000002_f32); + if cfg!(all(target_arch = "x86", not(target_feature = "sse2"))) && (res == -0.10703231_f32) + { + return; + } + assert_eq!(res, -0.10703229_f32); } } diff --git a/src/math/lgamma.rs b/src/math/lgamma.rs index 5bc87e8..a08bc5b 100644 --- a/src/math/lgamma.rs +++ b/src/math/lgamma.rs @@ -1,5 +1,6 @@ use super::lgamma_r; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn lgamma(x: f64) -> f64 { lgamma_r(x).0 } diff --git a/src/math/lgamma_r.rs b/src/math/lgamma_r.rs index 9533e88..b26177e 100644 --- a/src/math/lgamma_r.rs +++ b/src/math/lgamma_r.rs @@ -152,7 +152,7 @@ fn sin_pi(mut x: f64) -> f64 { x = 2.0 * (x * 0.5 - floor(x * 0.5)); /* x mod 2.0 */ n = (x * 4.0) as i32; - n = (n + 1) / 2; + n = div!(n + 1, 2); x -= (n as f64) * 0.5; x *= PI; @@ -164,6 +164,7 @@ fn sin_pi(mut x: f64) -> f64 { } } +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn lgamma_r(mut x: f64) -> (f64, i32) { let u: u64 = x.to_bits(); let mut t: f64; diff --git a/src/math/lgammaf.rs b/src/math/lgammaf.rs index dfdc87f..a9c2da7 100644 --- a/src/math/lgammaf.rs +++ b/src/math/lgammaf.rs @@ -1,5 +1,6 @@ use super::lgammaf_r; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn lgammaf(x: f32) -> f32 { lgammaf_r(x).0 } diff --git a/src/math/lgammaf_r.rs b/src/math/lgammaf_r.rs index c5e559f..723c90d 100644 --- a/src/math/lgammaf_r.rs +++ b/src/math/lgammaf_r.rs @@ -88,7 +88,7 @@ fn sin_pi(mut x: f32) -> f32 { x = 2.0 * (x * 0.5 - floorf(x * 0.5)); /* x mod 2.0 */ n = (x * 4.0) as isize; - n = (n + 1) / 2; + n = div!(n + 1, 2); y = (x as f64) - (n as f64) * 0.5; y *= 3.14159265358979323846; match n { @@ -99,6 +99,7 @@ fn sin_pi(mut x: f32) -> f32 { } } +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn lgammaf_r(mut x: f32) -> (f32, i32) { let u = x.to_bits(); let mut t: f32; diff --git a/src/math/mod.rs b/src/math/mod.rs index c8d7bd8..05ebb70 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -1,8 +1,6 @@ macro_rules! force_eval { ($e:expr) => { - unsafe { - ::core::ptr::read_volatile(&$e); - } + unsafe { ::core::ptr::read_volatile(&$e) } }; } @@ -58,6 +56,24 @@ macro_rules! i { }; } +// Temporary macro to avoid panic codegen for division (in debug mode too). At +// the time of this writing this is only used in a few places, and once +// rust-lang/rust#72751 is fixed then this macro will no longer be necessary and +// the native `/` operator can be used and panics won't be codegen'd. +#[cfg(any(debug_assertions, not(feature = "unstable")))] +macro_rules! div { + ($a:expr, $b:expr) => { + $a / $b + }; +} + +#[cfg(all(not(debug_assertions), feature = "unstable"))] +macro_rules! div { + ($a:expr, $b:expr) => { + unsafe { core::intrinsics::unchecked_div($a, $b) } + }; +} + macro_rules! llvm_intrinsically_optimized { (#[cfg($($clause:tt)*)] $e:expr) => { #[cfg(all(feature = "unstable", $($clause)*))] @@ -154,6 +170,8 @@ mod remainder; mod remainderf; mod remquo; mod remquof; +mod rint; +mod rintf; mod round; mod roundf; mod scalbn; @@ -268,6 +286,8 @@ pub use self::remainder::remainder; pub use self::remainderf::remainderf; pub use self::remquo::remquo; pub use self::remquof::remquof; +pub use self::rint::rint; +pub use self::rintf::rintf; pub use self::round::round; pub use self::roundf::roundf; pub use self::scalbn::scalbn; diff --git a/src/math/pow.rs b/src/math/pow.rs index 8b5989c..6a19ae6 100644 --- a/src/math/pow.rs +++ b/src/math/pow.rs @@ -299,8 +299,8 @@ pub fn pow(x: f64, y: f64) -> f64 { ax = with_set_high_word(ax, ix as u32); /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ - let u: f64 = ax - BP[k as usize]; /* bp[0]=1.0, bp[1]=1.5 */ - let v: f64 = 1.0 / (ax + BP[k as usize]); + let u: f64 = ax - i!(BP, k as usize); /* bp[0]=1.0, bp[1]=1.5 */ + let v: f64 = 1.0 / (ax + i!(BP, k as usize)); let ss: f64 = u * v; let s_h = with_set_low_word(ss, 0); @@ -309,7 +309,7 @@ pub fn pow(x: f64, y: f64) -> f64 { 0.0, ((ix as u32 >> 1) | 0x20000000) + 0x00080000 + ((k as u32) << 18), ); - let t_l: f64 = ax - (t_h - BP[k as usize]); + let t_l: f64 = ax - (t_h - i!(BP, k as usize)); let s_l: f64 = v * ((u - s_h * t_h) - s_h * t_l); /* compute log(ax) */ @@ -328,12 +328,12 @@ pub fn pow(x: f64, y: f64) -> f64 { let p_h: f64 = with_set_low_word(u + v, 0); let p_l = v - (p_h - u); let z_h: f64 = CP_H * p_h; /* cp_h+cp_l = 2/(3*log2) */ - let z_l: f64 = CP_L * p_h + p_l * CP + DP_L[k as usize]; + let z_l: f64 = CP_L * p_h + p_l * CP + i!(DP_L, k as usize); /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ let t: f64 = n as f64; - t1 = with_set_low_word(((z_h + z_l) + DP_H[k as usize]) + t, 0); - t2 = z_l - (((t1 - t) - DP_H[k as usize]) - z_h); + t1 = with_set_low_word(((z_h + z_l) + i!(DP_H, k as usize)) + t, 0); + t2 = z_l - (((t1 - t) - i!(DP_H, k as usize)) - z_h); } /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ @@ -484,6 +484,10 @@ mod tests { let exp = expected(*val); let res = computed(*val); + #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] + let exp = force_eval!(exp); + #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] + let res = force_eval!(res); assert!( if exp.is_nan() { res.is_nan() diff --git a/src/math/powf.rs b/src/math/powf.rs index f3cf76f..68d2083 100644 --- a/src/math/powf.rs +++ b/src/math/powf.rs @@ -238,8 +238,8 @@ pub fn powf(x: f32, y: f32) -> f32 { ax = f32::from_bits(ix as u32); /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ - u = ax - BP[k as usize]; /* bp[0]=1.0, bp[1]=1.5 */ - v = 1.0 / (ax + BP[k as usize]); + u = ax - i!(BP, k as usize); /* bp[0]=1.0, bp[1]=1.5 */ + v = 1.0 / (ax + i!(BP, k as usize)); s = u * v; s_h = s; is = s_h.to_bits() as i32; @@ -247,7 +247,7 @@ pub fn powf(x: f32, y: f32) -> f32 { /* t_h=ax+bp[k] High */ is = (((ix as u32 >> 1) & 0xfffff000) | 0x20000000) as i32; t_h = f32::from_bits(is as u32 + 0x00400000 + ((k as u32) << 21)); - t_l = ax - (t_h - BP[k as usize]); + t_l = ax - (t_h - i!(BP, k as usize)); s_l = v * ((u - s_h * t_h) - s_h * t_l); /* compute log(ax) */ s2 = s * s; @@ -267,13 +267,13 @@ pub fn powf(x: f32, y: f32) -> f32 { p_h = f32::from_bits(is as u32 & 0xfffff000); p_l = v - (p_h - u); z_h = CP_H * p_h; /* cp_h+cp_l = 2/(3*log2) */ - z_l = CP_L * p_h + p_l * CP + DP_L[k as usize]; + z_l = CP_L * p_h + p_l * CP + i!(DP_L, k as usize); /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ t = n as f32; - t1 = ((z_h + z_l) + DP_H[k as usize]) + t; + t1 = ((z_h + z_l) + i!(DP_H, k as usize)) + t; is = t1.to_bits() as i32; t1 = f32::from_bits(is as u32 & 0xfffff000); - t2 = z_l - (((t1 - t) - DP_H[k as usize]) - z_h); + t2 = z_l - (((t1 - t) - i!(DP_H, k as usize)) - z_h); }; /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ diff --git a/src/math/rem_pio2.rs b/src/math/rem_pio2.rs index 6b7dbd3..644616f 100644 --- a/src/math/rem_pio2.rs +++ b/src/math/rem_pio2.rs @@ -50,7 +50,12 @@ 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")))] + let tmp = 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 */ @@ -167,21 +172,21 @@ pub(crate) fn rem_pio2(x: f64) -> (i32, f64, f64) { let mut z = f64::from_bits(ui); let mut tx = [0.0; 3]; for i in 0..2 { - tx[i] = z as i32 as f64; - z = (z - tx[i]) * x1p24; + i!(tx,i, =, z as i32 as f64); + z = (z - i!(tx, i)) * x1p24; } - tx[2] = z; + i!(tx,2, =, z); /* skip zero terms, first term is non-zero */ let mut i = 2; - while i != 0 && tx[i] == 0.0 { + while i != 0 && i!(tx, i) == 0.0 { i -= 1; } let mut ty = [0.0; 3]; let n = rem_pio2_large(&tx[..=i], &mut ty, ((ix as i32) >> 20) - (0x3ff + 23), 1); if sign != 0 { - return (-n, -ty[0], -ty[1]); + return (-n, -i!(ty, 0), -i!(ty, 1)); } - (n, ty[0], ty[1]) + (n, i!(ty, 0), i!(ty, 1)) } #[cfg(test)] @@ -190,20 +195,28 @@ mod tests { #[test] fn test_near_pi() { + let arg = 3.141592025756836; + let arg = force_eval!(arg); assert_eq!( - rem_pio2(3.141592025756836), + rem_pio2(arg), (2, -6.278329573009626e-7, -2.1125998133974653e-23) ); + let arg = 3.141592033207416; + let arg = force_eval!(arg); assert_eq!( - rem_pio2(3.141592033207416), + rem_pio2(arg), (2, -6.20382377148128e-7, -2.1125998133974653e-23) ); + let arg = 3.141592144966125; + let arg = force_eval!(arg); assert_eq!( - rem_pio2(3.141592144966125), + rem_pio2(arg), (2, -5.086236681942706e-7, -2.1125998133974653e-23) ); + let arg = 3.141592979431152; + let arg = force_eval!(arg); assert_eq!( - rem_pio2(3.141592979431152), + rem_pio2(arg), (2, 3.2584135866119817e-7, -2.1125998133974653e-23) ); } diff --git a/src/math/rem_pio2_large.rs b/src/math/rem_pio2_large.rs index 002ce2e..db97a39 100644 --- a/src/math/rem_pio2_large.rs +++ b/src/math/rem_pio2_large.rs @@ -27,7 +27,7 @@ const INIT_JK: [usize; 4] = [3, 4, 4, 6]; // // NB: This table must have at least (e0-3)/24 + jk terms. // For quad precision (e0 <= 16360, jk = 6), this is 686. -#[cfg(target_pointer_width = "32")] +#[cfg(any(target_pointer_width = "32", target_pointer_width = "16"))] const IPIO2: [i32; 66] = [ 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, @@ -242,12 +242,12 @@ pub(crate) fn rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) -> let mut iq: [i32; 20] = [0; 20]; /* initialize jk*/ - let jk = INIT_JK[prec]; + let jk = i!(INIT_JK, prec); let jp = jk; /* determine jx,jv,q0, note that 3>q0 */ let jx = nx - 1; - let mut jv = (e0 - 3) / 24; + let mut jv = div!(e0 - 3, 24); if jv < 0 { jv = 0; } diff --git a/src/math/rem_pio2f.rs b/src/math/rem_pio2f.rs index 5d392ba..775f5d7 100644 --- a/src/math/rem_pio2f.rs +++ b/src/math/rem_pio2f.rs @@ -43,7 +43,12 @@ pub(crate) fn rem_pio2f(x: f32) -> (i32, f64) { if ix < 0x4dc90fdb { /* |x| ~< 2^28*(pi/2), medium size */ /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - let f_n = x64 * INV_PIO2 + TOINT - TOINT; + let tmp = x64 * INV_PIO2 + TOINT; + // 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")))] + let tmp = force_eval!(tmp); + let f_n = tmp - TOINT; return (f_n as i32, x64 - f_n * PIO2_1 - f_n * PIO2_1T); } if ix >= 0x7f800000 { diff --git a/src/math/rint.rs b/src/math/rint.rs new file mode 100644 index 0000000..0c6025c --- /dev/null +++ b/src/math/rint.rs @@ -0,0 +1,48 @@ +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] +pub fn rint(x: f64) -> f64 { + let one_over_e = 1.0 / f64::EPSILON; + let as_u64: u64 = x.to_bits(); + let exponent: u64 = as_u64 >> 52 & 0x7ff; + let is_positive = (as_u64 >> 63) == 0; + if exponent >= 0x3ff + 52 { + x + } else { + let ans = if is_positive { + x + one_over_e - one_over_e + } else { + x - one_over_e + one_over_e + }; + + if ans == 0.0 { + if is_positive { + 0.0 + } else { + -0.0 + } + } else { + ans + } + } +} + +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] +#[cfg(test)] +mod tests { + use super::rint; + + #[test] + fn negative_zero() { + assert_eq!(rint(-0.0_f64).to_bits(), (-0.0_f64).to_bits()); + } + + #[test] + fn sanity_check() { + assert_eq!(rint(-1.0), -1.0); + assert_eq!(rint(2.8), 3.0); + assert_eq!(rint(-0.5), -0.0); + assert_eq!(rint(0.5), 0.0); + assert_eq!(rint(-1.5), -2.0); + assert_eq!(rint(1.5), 2.0); + } +} diff --git a/src/math/rintf.rs b/src/math/rintf.rs new file mode 100644 index 0000000..d427793 --- /dev/null +++ b/src/math/rintf.rs @@ -0,0 +1,48 @@ +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] +pub fn rintf(x: f32) -> f32 { + let one_over_e = 1.0 / f32::EPSILON; + let as_u32: u32 = x.to_bits(); + let exponent: u32 = as_u32 >> 23 & 0xff; + let is_positive = (as_u32 >> 31) == 0; + if exponent >= 0x7f + 23 { + x + } else { + let ans = if is_positive { + x + one_over_e - one_over_e + } else { + x - one_over_e + one_over_e + }; + + if ans == 0.0 { + if is_positive { + 0.0 + } else { + -0.0 + } + } else { + ans + } + } +} + +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] +#[cfg(test)] +mod tests { + use super::rintf; + + #[test] + fn negative_zero() { + assert_eq!(rintf(-0.0_f32).to_bits(), (-0.0_f32).to_bits()); + } + + #[test] + fn sanity_check() { + assert_eq!(rintf(-1.0), -1.0); + assert_eq!(rintf(2.8), 3.0); + assert_eq!(rintf(-0.5), -0.0); + assert_eq!(rintf(0.5), 0.0); + assert_eq!(rintf(-1.5), -2.0); + assert_eq!(rintf(1.5), 2.0); + } +} diff --git a/src/math/round.rs b/src/math/round.rs index bf72f5b..46fabc9 100644 --- a/src/math/round.rs +++ b/src/math/round.rs @@ -1,38 +1,10 @@ +use super::copysign; +use super::trunc; use core::f64; -const TOINT: f64 = 1.0 / f64::EPSILON; - #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] -pub fn round(mut x: f64) -> f64 { - let i = x.to_bits(); - let e: u64 = i >> 52 & 0x7ff; - let mut y: f64; - - if e >= 0x3ff + 52 { - return x; - } - if e < 0x3ff - 1 { - // raise inexact if x!=0 - force_eval!(x + TOINT); - return 0.0 * x; - } - if i >> 63 != 0 { - x = -x; - } - y = x + TOINT - TOINT - x; - if y > 0.5 { - y = y + x - 1.0; - } else if y <= -0.5 { - y = y + x + 1.0; - } else { - y = y + x; - } - - if i >> 63 != 0 { - -y - } else { - y - } +pub fn round(x: f64) -> f64 { + trunc(x + copysign(0.5 - 0.25 * f64::EPSILON, x)) } #[cfg(test)] @@ -43,4 +15,14 @@ mod tests { fn negative_zero() { assert_eq!(round(-0.0_f64).to_bits(), (-0.0_f64).to_bits()); } + + #[test] + fn sanity_check() { + assert_eq!(round(-1.0), -1.0); + assert_eq!(round(2.8), 3.0); + assert_eq!(round(-0.5), -1.0); + assert_eq!(round(0.5), 1.0); + assert_eq!(round(-1.5), -2.0); + assert_eq!(round(1.5), 2.0); + } } diff --git a/src/math/roundf.rs b/src/math/roundf.rs index 497e88d..becdb56 100644 --- a/src/math/roundf.rs +++ b/src/math/roundf.rs @@ -1,38 +1,14 @@ +use super::copysignf; +use super::truncf; use core::f32; -const TOINT: f32 = 1.0 / f32::EPSILON; - #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] -pub fn roundf(mut x: f32) -> f32 { - let i = x.to_bits(); - let e: u32 = i >> 23 & 0xff; - let mut y: f32; - - if e >= 0x7f + 23 { - return x; - } - if e < 0x7f - 1 { - force_eval!(x + TOINT); - return 0.0 * x; - } - if i >> 31 != 0 { - x = -x; - } - y = x + TOINT - TOINT - x; - if y > 0.5f32 { - y = y + x - 1.0; - } else if y <= -0.5f32 { - y = y + x + 1.0; - } else { - y = y + x; - } - if i >> 31 != 0 { - -y - } else { - y - } +pub fn roundf(x: f32) -> f32 { + truncf(x + copysignf(0.5 - 0.25 * f32::EPSILON, x)) } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { use super::roundf; @@ -41,4 +17,14 @@ mod tests { fn negative_zero() { assert_eq!(roundf(-0.0_f32).to_bits(), (-0.0_f32).to_bits()); } + + #[test] + fn sanity_check() { + assert_eq!(roundf(-1.0), -1.0); + assert_eq!(roundf(2.8), 3.0); + assert_eq!(roundf(-0.5), -1.0); + assert_eq!(roundf(0.5), 1.0); + assert_eq!(roundf(-1.5), -2.0); + assert_eq!(roundf(1.5), 2.0); + } } diff --git a/src/math/sin.rs b/src/math/sin.rs index 1329b41..a53843d 100644 --- a/src/math/sin.rs +++ b/src/math/sin.rs @@ -81,5 +81,8 @@ pub fn sin(x: f64) -> f64 { fn test_near_pi() { let x = f64::from_bits(0x400921fb000FD5DD); // 3.141592026217707 let sx = f64::from_bits(0x3ea50d15ced1a4a2); // 6.273720864039205e-7 - assert_eq!(sin(x), sx); + let result = sin(x); + #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] + let result = force_eval!(result); + assert_eq!(result, sx); } diff --git a/src/math/sincos.rs b/src/math/sincos.rs index d49f65c..ff5d87a 100644 --- a/src/math/sincos.rs +++ b/src/math/sincos.rs @@ -12,6 +12,7 @@ use super::{get_high_word, k_cos, k_sin, rem_pio2}; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn sincos(x: f64) -> (f64, f64) { let s: f64; let c: f64; @@ -57,3 +58,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 + ); + } + } +} diff --git a/src/math/sincosf.rs b/src/math/sincosf.rs index 2725caa..9a4c361 100644 --- a/src/math/sincosf.rs +++ b/src/math/sincosf.rs @@ -23,6 +23,7 @@ const S2PIO2: f32 = 2.0 * PI_2; /* 0x400921FB, 0x54442D18 */ const S3PIO2: f32 = 3.0 * PI_2; /* 0x4012D97C, 0x7F3321D2 */ const S4PIO2: f32 = 4.0 * PI_2; /* 0x401921FB, 0x54442D18 */ +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn sincosf(x: f32) -> (f32, f32) { let s: f32; let c: f32; @@ -122,6 +123,8 @@ pub fn sincosf(x: f32) -> (f32, f32) { } } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { use super::sincosf; @@ -145,10 +148,38 @@ mod tests { let (s_minus, c_minus) = sincosf(theta - 2. * PI); const TOLERANCE: f32 = 1e-6; - assert!((s - s_plus).abs() < TOLERANCE); - assert!((s - s_minus).abs() < TOLERANCE); - assert!((c - c_plus).abs() < TOLERANCE); - assert!((c - c_minus).abs() < TOLERANCE); + 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 + ); } } } diff --git a/src/math/sqrtf.rs b/src/math/sqrtf.rs index ee868c8..00b20e5 100644 --- a/src/math/sqrtf.rs +++ b/src/math/sqrtf.rs @@ -128,6 +128,8 @@ pub fn sqrtf(x: f32) -> f32 { } } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { use super::*; diff --git a/src/math/tgamma.rs b/src/math/tgamma.rs index f8ccf66..e64eff6 100644 --- a/src/math/tgamma.rs +++ b/src/math/tgamma.rs @@ -38,7 +38,7 @@ fn sinpi(mut x: f64) -> f64 { /* reduce x into [-.25,.25] */ n = (4.0 * x) as isize; - n = (n + 1) / 2; + n = div!(n + 1, 2); x -= (n as f64) * 0.5; x *= PI; @@ -118,18 +118,19 @@ fn s(x: f64) -> f64 { /* to avoid overflow handle large x differently */ if x < 8.0 { for i in (0..=N).rev() { - num = num * x + SNUM[i]; - den = den * x + SDEN[i]; + num = num * x + i!(SNUM, i); + den = den * x + i!(SDEN, i); } } else { for i in 0..=N { - num = num / x + SNUM[i]; - den = den / x + SDEN[i]; + num = num / x + i!(SNUM, i); + den = den / x + i!(SDEN, i); } } return num / den; } +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn tgamma(mut x: f64) -> f64 { let u: u64 = x.to_bits(); let absx: f64; @@ -157,7 +158,7 @@ pub fn tgamma(mut x: f64) -> f64 { return 0.0 / 0.0; } if x <= FACT.len() as f64 { - return FACT[(x as usize) - 1]; + return i!(FACT, (x as usize) - 1); } } diff --git a/src/math/tgammaf.rs b/src/math/tgammaf.rs index a8f161f..23e3814 100644 --- a/src/math/tgammaf.rs +++ b/src/math/tgammaf.rs @@ -1,5 +1,6 @@ use super::tgamma; +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn tgammaf(x: f32) -> f32 { tgamma(x as f64) as f32 } diff --git a/src/math/truncf.rs b/src/math/truncf.rs index a4c0016..20d5b73 100644 --- a/src/math/truncf.rs +++ b/src/math/truncf.rs @@ -31,6 +31,8 @@ pub fn truncf(x: f32) -> f32 { f32::from_bits(i) } +// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 +#[cfg(not(target_arch = "powerpc64"))] #[cfg(test)] mod tests { #[test] |