Merge "Upgrade libm to 0.2.6"
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 @@
- 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
diff --git a/Android.bp b/Android.bp
index 8eacade..fde9dab 100644
--- a/Android.bp
+++ b/Android.bp
@@ -43,7 +43,7 @@
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"],
@@ -59,7 +59,7 @@
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 @@
...
+## [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 @@
- 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
diff --git a/Cargo.toml b/Cargo.toml
index e75fbba..5752ce6 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -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 @@
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 @@
[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"
diff --git a/METADATA b/METADATA
index 8b844f0..5735ca3 100644
--- a/METADATA
+++ b/METADATA
@@ -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 @@
}
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
}
}
diff --git a/README.md b/README.md
index 5d9e9bd..b864b5d 100644
--- a/README.md
+++ b/README.md
@@ -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/build.rs b/build.rs
index 9af6dec..80145a9 100644
--- a/build.rs
+++ b/build.rs
@@ -18,6 +18,7 @@
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 @@
// 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 @@
}
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 @@
--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 \
diff --git a/ci/run.sh b/ci/run.sh
index ed253ab..d0cd42a 100755
--- a/ci/run.sh
+++ b/ci/run.sh
@@ -5,6 +5,9 @@
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 @@
.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)
diff --git a/src/lib.rs b/src/lib.rs
index e228af9..29742b4 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -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 @@
}
}
+// 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 @@
/// 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 @@
/// 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 @@
///
/// 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 @@
///
/// 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 @@
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 @@
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 @@
///
/// 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 @@
///
/// 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 @@
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 @@
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 @@
/// 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 @@
/// 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 @@
/* 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 @@
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 @@
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 @@
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 @@
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 @@
/* 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 @@
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_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 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 @@
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 @@
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 @@
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 @@
-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::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 @@
* 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;
}
- z = adjusted_result as f32;
- z
+ 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);
+ }
}
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 @@
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 @@
}
#[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 @@
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 @@
}
}
+#[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 @@
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 @@
}
}
+#[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 @@
};
}
+// 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 remainderf;
mod remquo;
mod remquof;
+mod rint;
+mod rintf;
mod round;
mod roundf;
mod scalbn;
@@ -268,6 +286,8 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
/* 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 @@
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 @@
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 @@
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 @@
#[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 @@
//
// 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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
_ => (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 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 @@
}
}
+// 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 @@
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 @@
}
}
+// 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 @@
/* 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 @@
/* 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 @@
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 @@
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]
