user0/rust
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add fmaf128

Browse source 
Resolve all remaining `f64`-specific items in the generic version of
`fma`, then expose `fmaf128`.
This commit is contained in:
Trevor Gross 2025-02-05 23:45:14 +00:00 • committed by Trevor Gross
parent bbdcc7ef89
commit 9223d60dfa
14 changed files with 238 additions and 68 deletions
Download patch file Download diff file Expand all files Collapse all files

7
library/compiler-builtins/libm/crates/libm-macros/src/shared.rs
View file

@ -106,6 +106,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
None,
&["fma"],
),
(
// `(f128, f128, f128) -> f128`
FloatTy::F128,
Signature { args: &[Ty::F128, Ty::F128, Ty::F128], returns: &[Ty::F128] },
None,
&["fmaf128"],
),
(
// `(f32) -> i32`
FloatTy::F32,

1
library/compiler-builtins/libm/crates/libm-test/benches/icount.rs
View file

@ -108,6 +108,7 @@ main!(
icount_bench_floorf16_group,
icount_bench_floorf_group,
icount_bench_fma_group,
icount_bench_fmaf128_group,
icount_bench_fmaf_group,
icount_bench_fmax_group,
icount_bench_fmaxf128_group,

1
library/compiler-builtins/libm/crates/libm-test/benches/random.rs
View file

@ -127,6 +127,7 @@ libm_macros::for_each_function! {
| fdimf16
| floorf128
| floorf16
| fmaf128
| fmaxf128
| fmaxf16
| fminf128

23
library/compiler-builtins/libm/crates/libm-test/src/gen/case_list.rs
View file

@ -6,6 +6,9 @@
//!
//! This is useful for adding regression tests or expected failures.
#[cfg(f128_enabled)]
use libm::hf128;
use crate::{CheckBasis, CheckCtx, GeneratorKind, MathOp, op};
pub struct TestCase<Op: MathOp> {
@ -250,7 +253,7 @@ fn fma_cases() -> Vec<TestCase<op::fma::Routine>> {
TestCase::append_pairs(
&mut v,
&[
// Previously failure with incorrect sign
// Previous failure with incorrect sign
((5e-324, -5e-324, 0.0), Some(-0.0)),
],
);
@ -261,6 +264,24 @@ fn fmaf_cases() -> Vec<TestCase<op::fmaf::Routine>> {
vec![]
}
#[cfg(f128_enabled)]
fn fmaf128_cases() -> Vec<TestCase<op::fmaf128::Routine>> {
let mut v = vec![];
TestCase::append_pairs(
&mut v,
&[(
// Tricky rounding case that previously failed in extensive tests
(
hf128!("-0x1.1966cc01966cc01966cc01966f06p-25"),
hf128!("-0x1.669933fe69933fe69933fe6997c9p-16358"),
hf128!("-0x0.000000000000000000000000048ap-16382"),
),
Some(hf128!("0x0.c5171470a3ff5e0f68d751491b18p-16382")),
)],
);
v
}
fn fmax_cases() -> Vec<TestCase<op::fmax::Routine>> {
vec![]
}

2
library/compiler-builtins/libm/crates/libm-test/src/mpfloat.rs
View file

@ -196,7 +196,7 @@ libm_macros::for_each_function! {
expm1 | expm1f => exp_m1,
fabs | fabsf => abs,
fdim | fdimf | fdimf16 | fdimf128 => positive_diff,
fma | fmaf => mul_add,
fma | fmaf | fmaf128 => mul_add,
fmax | fmaxf | fmaxf16 | fmaxf128 => max,
fmin | fminf | fminf16 | fminf128 => min,
lgamma | lgammaf => ln_gamma,

2
library/compiler-builtins/libm/crates/libm-test/src/precision.rs
View file

@ -560,3 +560,5 @@ impl MaybeOverride<(f128, i32)> for SpecialCase {}
impl MaybeOverride<(f32, f32, f32)> for SpecialCase {}
impl MaybeOverride<(f64, f64, f64)> for SpecialCase {}
#[cfg(f128_enabled)]
impl MaybeOverride<(f128, f128, f128)> for SpecialCase {}

1
library/compiler-builtins/libm/crates/libm-test/tests/compare_built_musl.rs
View file

@ -99,6 +99,7 @@ libm_macros::for_each_function! {
fdimf16,
floorf128,
floorf16,
fmaf128,
fmaxf128,
fmaxf16,
fminf128,

1
library/compiler-builtins/libm/crates/util/src/main.rs
View file

@ -96,6 +96,7 @@ fn do_eval(basis: &str, op: &str, inputs: &[&str]) {
| fdimf16
| floorf128
| floorf16
| fmaf128
| fmaxf128
| fmaxf16
| fminf128

7
library/compiler-builtins/libm/etc/function-definitions.json
View file

@ -356,6 +356,13 @@
],
"type": "f32"
},
"fmaf128": {
"sources": [
"src/math/fmaf128.rs",
"src/math/generic/fma.rs"
],
"type": "f128"
},
"fmax": {
"sources": [
"src/math/fmax.rs",

1
library/compiler-builtins/libm/etc/function-list.txt
View file

@ -53,6 +53,7 @@ floorf128
floorf16
fma
fmaf
fmaf128
fmax
fmaxf
fmaxf128

1
library/compiler-builtins/libm/src/libm_helper.rs
View file

@ -208,6 +208,7 @@ libm_helper! {
(fn fabs(x: f128) -> (f128); => fabsf128);
(fn fdim(x: f128, y: f128) -> (f128); => fdimf128);
(fn floor(x: f128) -> (f128); => floorf128);
(fn fmaf128(x: f128, y: f128, z: f128) -> (f128); => fmaf128);
(fn fmax(x: f128, y: f128) -> (f128); => fmaxf128);
(fn fmin(x: f128, y: f128) -> (f128); => fminf128);
(fn fmod(x: f128, y: f128) -> (f128); => fmodf128);

7
library/compiler-builtins/libm/src/math/fmaf128.rs Normal file
View file

@ -0,0 +1,7 @@
/// Fused multiply add (f128)
///
/// Computes `(x*y)+z`, rounded as one ternary operation (i.e. calculated with infinite precision).
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
pub fn fmaf128(x: f128, y: f128, z: f128) -> f128 {
return super::generic::fma(x, y, z);
}

250
library/compiler-builtins/libm/src/math/generic/fma.rs
View file

@ -1,10 +1,11 @@
/* SPDX-License-Identifier: MIT */
/* origin: musl src/math/fma.c. Ported to generic Rust algorithm in 2025, TG. */
use core::{f32, f64};
use super::super::support::{DInt, HInt, IntTy};
use super::super::{CastFrom, CastInto, Float, Int, MinInt};
const ZEROINFNAN: i32 = 0x7ff - 0x3ff - 52 - 1;
/// Fused multiply-add that works when there is not a larger float size available. Currently this
/// is still specialized only for `f64`. Computes `(x * y) + z`.
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
@ -18,79 +19,99 @@ where
{
let one = IntTy::<F>::ONE;
let zero = IntTy::<F>::ZERO;
let magic = F::from_parts(false, F::BITS - 1 + F::EXP_BIAS, zero);
/* normalize so top 10bits and last bit are 0 */
// Normalize such that the top of the mantissa is zero and we have a guard bit.
let nx = Norm::from_float(x);
let ny = Norm::from_float(y);
let nz = Norm::from_float(z);
if nx.e >= ZEROINFNAN || ny.e >= ZEROINFNAN {
if nx.is_zero_nan_inf() || ny.is_zero_nan_inf() {
// Value will overflow, defer to non-fused operations.
return x * y + z;
}
if nz.e >= ZEROINFNAN {
if nz.e > ZEROINFNAN {
/* z==0 */
if nz.is_zero_nan_inf() {
if nz.is_zero() {
// Empty add component means we only need to multiply.
return x * y;
}
// `z` is NaN or infinity, which sets the result.
return z;
}
/* mul: r = x*y */
// multiply: r = x * y
let zhi: F::Int;
let zlo: F::Int;
let (mut rlo, mut rhi) = nx.m.widen_mul(ny.m).lo_hi();
/* either top 20 or 21 bits of rhi and last 2 bits of rlo are 0 */
/* align exponents */
// Exponent result of multiplication
let mut e: i32 = nx.e + ny.e;
// Needed shift to align `z` to the multiplication result
let mut d: i32 = nz.e - e;
let sbits = F::BITS as i32;
/* shift bits z<<=kz, r>>=kr, so kz+kr == d, set e = e+kr (== ez-kz) */
// Scale `z`. Shift `z <<= kz`, `r >>= kr`, so `kz+kr == d`, set `e = e+kr` (== ez-kz)
if d > 0 {
// The magnitude of `z` is larger than `x * y`
if d < sbits {
// Maximum shift of one `F::BITS` means shifted `z` will fit into `2 * F::BITS`. Shift
// it into `(zhi, zlo)`. No exponent adjustment necessary.
zlo = nz.m << d;
zhi = nz.m >> (sbits - d);
} else {
// Shift larger than `sbits`, `z` only needs the top half `zhi`. Place it there (acts
// as a shift by `sbits`).
zlo = zero;
zhi = nz.m;
e = nz.e - sbits;
d -= sbits;
// `z`'s exponent is large enough that it now needs to be taken into account.
e = nz.e - sbits;
if d == 0 {
// Exactly `sbits`, nothing to do
} else if d < sbits {
rlo = (rhi << (sbits - d))
| (rlo >> d)
| IntTy::<F>::from((rlo << (sbits - d)) != zero);
// Remaining shift fits within `sbits`. Leave `z` in place, shift `x * y`
rlo = (rhi << (sbits - d)) | (rlo >> d);
// Set the sticky bit
rlo |= IntTy::<F>::from((rlo << (sbits - d)) != zero);
rhi = rhi >> d;
} else {
// `z`'s magnitude is enough that `x * y` is irrelevant. It was nonzero, so set
// the sticky bit.
rlo = one;
rhi = zero;
}
}
} else {
// `z`'s magnitude once shifted fits entirely within `zlo`
zhi = zero;
d = -d;
if d == 0 {
// No shift needed
zlo = nz.m;
} else if d < sbits {
zlo = (nz.m >> d) | IntTy::<F>::from((nz.m << (sbits - d)) != zero);
// Shift s.t. `nz.m` fits into `zlo`
let sticky = IntTy::<F>::from((nz.m << (sbits - d)) != zero);
zlo = (nz.m >> d) | sticky;
} else {
// Would be entirely shifted out, only set the sticky bit
zlo = one;
}
}
/* add */
/* addition */
let mut neg = nx.neg ^ ny.neg;
let samesign: bool = !neg ^ nz.neg;
let mut nonzero: i32 = 1;
let mut rhi_nonzero = true;
if samesign {
/* r += z */
// r += z
rlo = rlo.wrapping_add(zlo);
rhi += zhi + IntTy::<F>::from(rlo < zlo);
} else {
/* r -= z */
// r -= z
let (res, borrow) = rlo.overflowing_sub(zlo);
rlo = res;
rhi = rhi.wrapping_sub(zhi.wrapping_add(IntTy::<F>::from(borrow)));
@ -99,129 +120,226 @@ where
rhi = rhi.signed().wrapping_neg().unsigned() - IntTy::<F>::from(rlo != zero);
neg = !neg;
}
nonzero = (rhi != zero) as i32;
rhi_nonzero = rhi != zero;
}
/* set rhi to top 63bit of the result (last bit is sticky) */
if nonzero != 0 {
/* Construct result */
// Shift result into `rhi`, left-aligned. Last bit is sticky
if rhi_nonzero {
// `d` > 0, need to shift both `rhi` and `rlo` into result
e += sbits;
d = rhi.leading_zeros() as i32 - 1;
/* note: d > 0 */
rhi = (rhi << d) | (rlo >> (sbits - d)) | IntTy::<F>::from((rlo << d) != zero);
rhi = (rhi << d) | (rlo >> (sbits - d));
// Update sticky
rhi |= IntTy::<F>::from((rlo << d) != zero);
} else if rlo != zero {
// `rhi` is zero, `rlo` is the entire result and needs to be shifted
d = rlo.leading_zeros() as i32 - 1;
if d < 0 {
// Shift and set sticky
rhi = (rlo >> 1) | (rlo & one);
} else {
rhi = rlo << d;
}
} else {
/* exact +-0 */
// exact +/- 0.0
return x * y + z;
}
e -= d;
/* convert to double */
let mut i: F::SignedInt = rhi.signed(); /* i is in [1<<62,(1<<63)-1] */
// Use int->float conversion to populate the significand.
// i is in [1 << (BITS - 2), (1 << (BITS - 1)) - 1]
let mut i: F::SignedInt = rhi.signed();
if neg {
i = -i;
}
let mut r: F = F::cast_from_lossy(i); /* |r| is in [0x1p62,0x1p63] */
// `|r|` is in `[0x1p62,0x1p63]` for `f64`
let mut r: F = F::cast_from_lossy(i);
if e < -(F::EXP_BIAS as i32 - 1) - (sbits - 2) {
/* result is subnormal before rounding */
if e == -(F::EXP_BIAS as i32 - 1) - (sbits - 1) {
let mut c: F = magic;
/* Account for subnormal and rounding */
// Unbiased exponent for the maximum value of `r`
let max_pow = F::BITS - 1 + F::EXP_BIAS;
if e < -(max_pow as i32 - 2) {
// Result is subnormal before rounding
if e == -(max_pow as i32 - 1) {
let mut c = F::from_parts(false, max_pow, zero);
if neg {
c = -c;
}
if r == c {
/* min normal after rounding, underflow depends
* on arch behaviour which can be imitated by
* a double to float conversion */
return r.raise_underflow();
// Min normal after rounding,
return r.raise_underflow_ret_self();
}
/* one bit is lost when scaled, add another top bit to
* only round once at conversion if it is inexact */
if (rhi << F::SIG_BITS) != zero {
let iu: F::Int = (rhi >> 1) | (rhi & one) | (one << 62);
if (rhi << (F::SIG_BITS + 1)) != zero {
// Account for truncated bits. One bit will be lost in the `scalbn` call, add
// another top bit to avoid double rounding if inexact.
let iu: F::Int = (rhi >> 1) | (rhi & one) | (one << (F::BITS - 2));
i = iu.signed();
if neg {
i = -i;
}
r = F::cast_from_lossy(i);
r = F::cast_from(2i8) * r - c; /* remove top bit */
/* raise underflow portably, such that it
* cannot be optimized away */
r += r.raise_underflow2();
r = F::cast_from_lossy(i);
// Remove the top bit
r = F::cast_from(2i8) * r - c;
r += r.raise_underflow_ret_zero();
}
} else {
/* only round once when scaled */
d = 10;
i = (((rhi >> d) | IntTy::<F>::from(rhi << (F::BITS as i32 - d) != zero)) << d)
.signed();
// Only round once when scaled
d = F::EXP_BITS as i32 - 1;
let sticky = IntTy::<F>::from(rhi << (F::BITS as i32 - d) != zero);
i = (((rhi >> d) | sticky) << d).signed();
if neg {
i = -i;
}
r = F::cast_from(i);
r = F::cast_from_lossy(i);
}
}
// Use our exponent to scale the final value.
super::scalbn(r, e)
}
/// Representation of `F` that has handled subnormals.
#[derive(Clone, Copy, Debug)]
struct Norm<F: Float> {
/// Normalized significand with one guard bit.
/// Normalized significand with one guard bit, unsigned.
m: F::Int,
/// Unbiased exponent, normalized.
/// Exponent of the mantissa such that `m * 2^e = x`. Accounts for the shift in the mantissa
/// and the guard bit; that is, 1.0 will normalize as `m = 1 << 53` and `e = -53`.
e: i32,
neg: bool,
}
impl<F: Float> Norm<F> {
/// Unbias the exponent and account for the mantissa's precision, including the guard bit.
const EXP_UNBIAS: u32 = F::EXP_BIAS + F::SIG_BITS + 1;
/// Values greater than this had a saturated exponent (infinity or NaN), OR were zero and we
/// adjusted the exponent such that it exceeds this threashold.
const ZERO_INF_NAN: u32 = F::EXP_SAT - Self::EXP_UNBIAS;
fn from_float(x: F) -> Self {
let mut ix = x.to_bits();
let mut e = x.exp() as i32;
let neg = x.is_sign_negative();
if e == 0 {
// Normalize subnormals by multiplication
let magic = F::from_parts(false, F::BITS - 1 + F::EXP_BIAS, F::Int::ZERO);
let scaled = x * magic;
let scale_i = F::BITS - 1;
let scale_f = F::from_parts(false, scale_i + F::EXP_BIAS, F::Int::ZERO);
let scaled = x * scale_f;
ix = scaled.to_bits();
e = scaled.exp() as i32;
e = if e != 0 { e - (F::BITS as i32 - 1) } else { 0x800 };
e = if e == 0 {
// If the exponent is still zero, the input was zero. Artifically set this value
// such that the final `e` will exceed `ZERO_INF_NAN`.
1 << F::EXP_BITS
} else {
// Otherwise, account for the scaling we just did.
e - scale_i as i32
};
}
e -= F::EXP_BIAS as i32 + 52 + 1;
e -= Self::EXP_UNBIAS as i32;
// Absolute value, set the implicit bit, and shift to create a guard bit
ix &= F::SIG_MASK;
ix |= F::IMPLICIT_BIT;
ix <<= 1; // add a guard bit
ix <<= 1;
Self { m: ix, e, neg }
}
/// True if the value was zero, infinity, or NaN.
fn is_zero_nan_inf(self) -> bool {
self.e >= Self::ZERO_INF_NAN as i32
}
/// The only value we have
fn is_zero(self) -> bool {
// The only exponent that strictly exceeds this value is our sentinel value for zero.
self.e > Self::ZERO_INF_NAN as i32
}
}
/// Type-specific helpers that are not needed outside of fma.
pub trait FmaHelper {
fn raise_underflow(self) -> Self;
fn raise_underflow2(self) -> Self;
fn raise_underflow_ret_self(self) -> Self;
fn raise_underflow_ret_zero(self) -> Self;
}
impl FmaHelper for f64 {
fn raise_underflow(self) -> Self {
let x0_ffffff8p_63 = f64::from_bits(0x3bfffffff0000000); // 0x0.ffffff8p-63
let fltmin: f32 = (x0_ffffff8p_63 * f32::MIN_POSITIVE as f64 * self) as f32;
fn raise_underflow_ret_self(self) -> Self {
/* min normal after rounding, underflow depends
* on arch behaviour which can be imitated by
* a double to float conversion */
let fltmin: f32 = (hf64!("0x0.ffffff8p-63") * f32::MIN_POSITIVE as f64 * self) as f32;
f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * fltmin as f64
}
fn raise_underflow2(self) -> Self {
fn raise_underflow_ret_zero(self) -> Self {
/* raise underflow portably, such that it
* cannot be optimized away */
let tiny: f64 = f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * self;
(tiny * tiny) * (self - self)
}
}
#[cfg(f128_enabled)]
impl FmaHelper for f128 {
fn raise_underflow_ret_self(self) -> Self {
self
}
fn raise_underflow_ret_zero(self) -> Self {
f128::ZERO
}
}
#[cfg(test)]
mod tests {
use super::*;
fn spec_test<F>()
where
F: Float + FmaHelper,
F: CastFrom<F::SignedInt>,
F: CastFrom<i8>,
F::Int: HInt,
u32: CastInto<F::Int>,
{
let x = F::from_bits(F::Int::ONE);
let y = F::from_bits(F::Int::ONE);
let z = F::ZERO;
// 754-2020 says "When the exact result of (a × b) + c is non-zero yet the result of
// fusedMultiplyAdd is zero because of rounding, the zero result takes the sign of the
// exact result"
assert_biteq!(fma(x, y, z), F::ZERO);
assert_biteq!(fma(x, -y, z), F::NEG_ZERO);
assert_biteq!(fma(-x, y, z), F::NEG_ZERO);
assert_biteq!(fma(-x, -y, z), F::ZERO);
}
#[test]
fn spec_test_f64() {
spec_test::<f64>();
}
#[test]
#[cfg(f128_enabled)]
fn spec_test_f128() {
spec_test::<f128>();
}
}

2
library/compiler-builtins/libm/src/math/mod.rs
View file

@ -385,6 +385,7 @@ cfg_if! {
mod fabsf128;
mod fdimf128;
mod floorf128;
mod fmaf128;
mod fmaxf128;
mod fminf128;
mod fmodf128;
@ -402,6 +403,7 @@ cfg_if! {
pub use self::fabsf128::fabsf128;
pub use self::fdimf128::fdimf128;
pub use self::floorf128::floorf128;
pub use self::fmaf128::fmaf128;
pub use self::fmaxf128::fmaxf128;
pub use self::fminf128::fminf128;
pub use self::fmodf128::fmodf128;