Skip to content

Commit

Permalink
Merge pull request #510 from tgross35/replace-fenv
Browse files Browse the repository at this point in the history 
Migrate away from nonfunctional `fenv` stubs
  • Loading branch information
@tgross35
tgross35 authored Feb 10, 2025
2 parents 6fab367 + e8fbb05 commit 5151be0
Show file tree
Hide file tree
Showing 12 changed files with 470 additions and 178 deletions.
1 change: 1 addition & 0 deletions crates/libm-test/src/f8_impl.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ impl Float for f8 {
const INFINITY: Self = Self(0b0_1111_000);
const NEG_INFINITY: Self = Self(0b1_1111_000);
const NAN: Self = Self(0b0_1111_100);
const MIN_POSITIVE_NORMAL: Self = Self(1 << Self::SIG_BITS);
// FIXME: incorrect values
const EPSILON: Self = Self::ZERO;
const PI: Self = Self::ZERO;
Expand Down
25 changes: 13 additions & 12 deletions src/math/cbrt.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,15 @@
*/

use super::Float;
use super::fenv::Rounding;
use super::support::cold_path;
use super::support::{FpResult, Round, cold_path};

/// Compute the cube root of the argument.
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
pub fn cbrt(x: f64) -> f64 {
cbrt_round(x, Round::Nearest).val
}

pub fn cbrt_round(x: f64, round: Round) -> FpResult<f64> {
const ESCALE: [f64; 3] = [
1.0,
hf64!("0x1.428a2f98d728bp+0"), /* 2^(1/3) */
Expand All @@ -33,8 +36,6 @@ pub fn cbrt(x: f64) -> f64 {

let off = [hf64!("0x1p-53"), 0.0, 0.0, 0.0];

let rm = Rounding::get();

/* rm=0 for rounding to nearest, and other values for directed roundings */
let hx: u64 = x.to_bits();
let mut mant: u64 = hx & f64::SIG_MASK;
Expand All @@ -51,7 +52,7 @@ pub fn cbrt(x: f64) -> f64 {
to that for x a signaling NaN, it correctly triggers
the invalid exception. */
if e == f64::EXP_SAT || ix == 0 {
return x + x;
return FpResult::ok(x + x);
}

let nz = ix.leading_zeros() - 11; /* subnormal */
Expand Down Expand Up @@ -124,8 +125,8 @@ pub fn cbrt(x: f64) -> f64 {
* from ulp(1);
* for rounding to nearest, ady0 is tiny when dy is near from 1/2 ulp(1),
* or from 3/2 ulp(1). */
let mut ady0: f64 = (ady - off[rm as usize]).abs();
let mut ady1: f64 = (ady - (hf64!("0x1p-52") + off[rm as usize])).abs();
let mut ady0: f64 = (ady - off[round as usize]).abs();
let mut ady1: f64 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs();

if ady0 < hf64!("0x1p-75") || ady1 < hf64!("0x1p-75") {
cold_path();
Expand All @@ -140,8 +141,8 @@ pub fn cbrt(x: f64) -> f64 {
dy = (y1 - y) - dy;
y1 = y;
ady = dy.abs();
ady0 = (ady - off[rm as usize]).abs();
ady1 = (ady - (hf64!("0x1p-52") + off[rm as usize])).abs();
ady0 = (ady - off[round as usize]).abs();
ady1 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs();

if ady0 < hf64!("0x1p-98") || ady1 < hf64!("0x1p-98") {
cold_path();
Expand All @@ -157,7 +158,7 @@ pub fn cbrt(x: f64) -> f64 {
y1 = hf64!("0x1.de87aa837820fp+0").copysign(zz);
}

if rm != Rounding::Nearest {
if round != Round::Nearest {
let wlist = [
(hf64!("0x1.3a9ccd7f022dbp+0"), hf64!("0x1.1236160ba9b93p+0")), // ~ 0x1.1236160ba9b930000000000001e7e8fap+0
(hf64!("0x1.7845d2faac6fep+0"), hf64!("0x1.23115e657e49cp+0")), // ~ 0x1.23115e657e49c0000000000001d7a799p+0
Expand All @@ -170,7 +171,7 @@ pub fn cbrt(x: f64) -> f64 {

for (a, b) in wlist {
if azz == a {
let tmp = if rm as u64 + sign == 2 { hf64!("0x1p-52") } else { 0.0 };
let tmp = if round as u64 + sign == 2 { hf64!("0x1p-52") } else { 0.0 };
y1 = (b + tmp).copysign(zz);
}
}
Expand All @@ -194,7 +195,7 @@ pub fn cbrt(x: f64) -> f64 {
}
}

f64::from_bits(cvt3)
FpResult::ok(f64::from_bits(cvt3))
}

fn fmaf64(x: f64, y: f64, z: f64) -> f64 {
Expand Down
49 changes: 0 additions & 49 deletions src/math/fenv.rs
View file

This file was deleted.

91 changes: 76 additions & 15 deletions src/math/generic/ceil.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,30 +7,35 @@
//! performance seems to be better (based on icount) and it does not seem to experience rounding
//! errors on i386.
use super::super::support::{FpResult, Status};
use super::super::{Float, Int, IntTy, MinInt};

pub fn ceil<F: Float>(x: F) -> F {
ceil_status(x).val
}

pub fn ceil_status<F: Float>(x: F) -> FpResult<F> {
let zero = IntTy::<F>::ZERO;

let mut ix = x.to_bits();
let e = x.exp_unbiased();

// If the represented value has no fractional part, no truncation is needed.
if e >= F::SIG_BITS as i32 {
return x;
return FpResult::ok(x);
}

if e >= 0 {
let status;
let res = if e >= 0 {
// |x| >= 1.0

let m = F::SIG_MASK >> e.unsigned();
if (ix & m) == zero {
// Portion to be masked is already zero; no adjustment needed.
return x;
return FpResult::ok(x);
}

// Otherwise, raise an inexact exception.
force_eval!(x + F::MAX);
status = Status::INEXACT;

if x.is_sign_positive() {
ix += m;
Expand All @@ -40,7 +45,11 @@ pub fn ceil<F: Float>(x: F) -> F {
F::from_bits(ix)
} else {
// |x| < 1.0, raise an inexact exception since truncation will happen (unless x == 0).
force_eval!(x + F::MAX);
if ix & F::SIG_MASK == F::Int::ZERO {
status = Status::OK;
} else {
status = Status::INEXACT;
}

if x.is_sign_negative() {
// -1.0 < x <= -0.0; rounding up goes toward -0.0.
Expand All @@ -52,18 +61,30 @@ pub fn ceil<F: Float>(x: F) -> F {
// +0.0 remains unchanged
x
}
}
};

FpResult::new(res, status)
}

#[cfg(test)]
mod tests {
use super::*;
use crate::support::Hexf;

/// Test against https://en.cppreference.com/w/cpp/numeric/math/ceil
fn spec_test<F: Float>() {
// Not Asserted: that the current rounding mode has no effect.
for f in [F::ZERO, F::NEG_ZERO, F::INFINITY, F::NEG_INFINITY].iter().copied() {
assert_biteq!(ceil(f), f);
fn spec_test<F: Float>(cases: &[(F, F, Status)]) {
let roundtrip = [F::ZERO, F::ONE, F::NEG_ONE, F::NEG_ZERO, F::INFINITY, F::NEG_INFINITY];

for x in roundtrip {
let FpResult { val, status } = ceil_status(x);
assert_biteq!(val, x, "{}", Hexf(x));
assert_eq!(status, Status::OK, "{}", Hexf(x));
}

for &(x, res, res_stat) in cases {
let FpResult { val, status } = ceil_status(x);
assert_biteq!(val, res, "{}", Hexf(x));
assert_eq!(status, res_stat, "{}", Hexf(x));
}
}

Expand All @@ -72,7 +93,17 @@ mod tests {
#[test]
#[cfg(f16_enabled)]
fn spec_tests_f16() {
spec_test::<f16>();
let cases = [
(0.1, 1.0, Status::INEXACT),
(-0.1, -0.0, Status::INEXACT),
(0.9, 1.0, Status::INEXACT),
(-0.9, -0.0, Status::INEXACT),
(1.1, 2.0, Status::INEXACT),
(-1.1, -1.0, Status::INEXACT),
(1.9, 2.0, Status::INEXACT),
(-1.9, -1.0, Status::INEXACT),
];
spec_test::<f16>(&cases);
}

#[test]
Expand All @@ -83,7 +114,17 @@ mod tests {

#[test]
fn spec_tests_f32() {
spec_test::<f32>();
let cases = [
(0.1, 1.0, Status::INEXACT),
(-0.1, -0.0, Status::INEXACT),
(0.9, 1.0, Status::INEXACT),
(-0.9, -0.0, Status::INEXACT),
(1.1, 2.0, Status::INEXACT),
(-1.1, -1.0, Status::INEXACT),
(1.9, 2.0, Status::INEXACT),
(-1.9, -1.0, Status::INEXACT),
];
spec_test::<f32>(&cases);
}

#[test]
Expand All @@ -94,12 +135,32 @@ mod tests {

#[test]
fn spec_tests_f64() {
spec_test::<f64>();
let cases = [
(0.1, 1.0, Status::INEXACT),
(-0.1, -0.0, Status::INEXACT),
(0.9, 1.0, Status::INEXACT),
(-0.9, -0.0, Status::INEXACT),
(1.1, 2.0, Status::INEXACT),
(-1.1, -1.0, Status::INEXACT),
(1.9, 2.0, Status::INEXACT),
(-1.9, -1.0, Status::INEXACT),
];
spec_test::<f64>(&cases);
}

#[test]
#[cfg(f128_enabled)]
fn spec_tests_f128() {
spec_test::<f128>();
let cases = [
(0.1, 1.0, Status::INEXACT),
(-0.1, -0.0, Status::INEXACT),
(0.9, 1.0, Status::INEXACT),
(-0.9, -0.0, Status::INEXACT),
(1.1, 2.0, Status::INEXACT),
(-1.1, -1.0, Status::INEXACT),
(1.9, 2.0, Status::INEXACT),
(-1.9, -1.0, Status::INEXACT),
];
spec_test::<f128>(&cases);
}
}
Loading

0 comments on commit 5151be0

Please sign in to comment.