Add fmaf16 and fmaf128

Author: tgross35

crates/libm-macros/src/shared.rs

@@ -92,6 +92,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         None,
         &["copysignf128"],
     ),
+    (
+        // `(f16, f16, f16) -> f16`
+        FloatTy::F16,
+        Signature { args: &[Ty::F16, Ty::F16, Ty::F16], returns: &[Ty::F16] },
+        None,
+        &["fmaf16"],
+    ),
     (
         // `(f32, f32, f32) -> f32`
         FloatTy::F32,
@@ -106,6 +113,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,

crates/libm-test/src/mpfloat.rs

@@ -144,7 +144,7 @@ libm_macros::for_each_function! {
         expm1 | expm1f => exp_m1,
         fabs | fabsf => abs,
         fdim | fdimf => positive_diff,
-        fma | fmaf => mul_add,
+        fma | fmaf | fmaf16 | fmaf128 => mul_add,
         fmax | fmaxf => max,
         fmin | fminf => min,
         lgamma | lgammaf => ln_gamma,

crates/libm-test/src/precision.rs

@@ -401,7 +401,12 @@ fn bessel_prec_dropoff<F: Float>(
     None
 }
 
+#[cfg(f16_enabled)]
+impl MaybeOverride<(f16, f16, f16)> 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 {}
+
 impl MaybeOverride<(f32, i32)> for SpecialCase {}
 impl MaybeOverride<(f64, i32)> for SpecialCase {}

crates/libm-test/tests/compare_built_musl.rs

@@ -47,7 +47,7 @@ where
 libm_macros::for_each_function! {
     callback: musl_rand_tests,
     // Musl does not support `f16` and `f128` on all platforms.
-    skip: [copysignf16, copysignf128, fabsf16, fabsf128],
+    skip: [copysignf16, copysignf128, fabsf16, fabsf128, fmaf16, fmaf128],
     attributes: [
         #[cfg_attr(x86_no_sse, ignore)] // FIXME(correctness): wrong result on i586
         [exp10, exp10f, exp2, exp2f, rint]

crates/libm-test/tests/multiprecision.rs

@@ -136,6 +136,8 @@ libm_macros::for_each_function! {
         fdimf,
         fma,
         fmaf,
+        fmaf16,
+        fmaf128,
         fmax,
         fmaxf,
         fmin,

etc/function-definitions.json

@@ -318,16 +318,32 @@
     "fma": {
         "sources": [
             "src/libm_helper.rs",
-            "src/math/fma.rs"
+            "src/math/fma.rs",
+            "src/math/generic/fma.rs"
         ],
         "type": "f64"
     },
     "fmaf": {
         "sources": [
-            "src/math/fmaf.rs"
+            "src/math/fmaf.rs",
+            "src/math/generic/fma.rs"
         ],
         "type": "f32"
     },
+    "fmaf128": {
+        "sources": [
+            "src/math/fmaf128.rs",
+            "src/math/generic/fma.rs"
+        ],
+        "type": "f128"
+    },
+    "fmaf16": {
+        "sources": [
+            "src/math/fmaf16.rs",
+            "src/math/generic/fma.rs"
+        ],
+        "type": "f16"
+    },
     "fmax": {
         "sources": [
             "src/libm_helper.rs",

etc/function-list.txt

@@ -47,6 +47,8 @@ floor
 floorf
 fma
 fmaf
+fmaf128
+fmaf16
 fmax
 fmaxf
 fmin

src/math/fmaf.rs

@@ -47,6 +47,10 @@ use super::fenv::{
 /// according to the rounding mode characterized by the value of FLT_ROUNDS.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 {
+    if true {
+        return super::generic::fma_big::<f32, f64>(x, y, z);
+    }
+
     let xy: f64;
     let mut result: f64;
     let mut ui: u64;

src/math/fmaf128.rs

@@ -0,0 +1,4 @@
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn fmaf128(x: f128, y: f128, z: f128) -> f128 {
+    super::generic::fma(x, y, z)
+}

src/math/fmaf16.rs

@@ -0,0 +1,4 @@
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn fmaf16(x: f16, y: f16, z: f16) -> f16 {
+    super::generic::fma_big::<f16, f32>(x, y, z)
+}

src/math/generic/fma.rs

@@ -0,0 +1,58 @@
+#![allow(unused)]
+
+use super::super::fenv::{
+    FE_INEXACT, FE_TONEAREST, FE_UNDERFLOW, feclearexcept, fegetround, feraiseexcept, fetestexcept,
+};
+use super::super::{CastFrom, CastInto, Float, IntTy, MinInt};
+
+/// Fused multiply add.
+pub fn fma<F: Float>(x: F, y: F, z: F) -> F {
+    todo!()
+}
+
+pub fn fma_big<F: Float, B: Float>(x: F, y: F, z: F) -> F
+where
+    F: CastInto<B>,
+    B: CastInto<F>,
+    B::Int: CastInto<i32>,
+    i32: CastFrom<i32>,
+    B: CastFrom<F>,
+{
+    let one = IntTy::<B>::ONE;
+
+    let xy: B;
+    let mut result: B;
+    let mut ui: B::Int;
+    let e: i32;
+
+    xy = x.cast() * y.cast();
+    result = xy + z.cast();
+    ui = result.to_bits();
+    e = i32::cast_from(ui >> F::SIG_BITS) & F::EXP_MAX as i32;
+    let zb = B::cast_from(z);
+
+    let prec_diff = B::SIG_BITS - F::SIG_BITS;
+    let excess_prec = ui & ((one << prec_diff) - one);
+    let x = one << (prec_diff - 1);
+
+    // Common case: the larger precision is fine
+    if excess_prec != x
+        || e == i32::cast_from(F::EXP_MAX)
+        || (result - xy == zb && result - zb == xy)
+        || fegetround() != FE_TONEAREST
+    {
+        // TODO: feclearexcept
+
+        return result.cast();
+    }
+
+    let neg = ui & B::SIGN_MASK > IntTy::<B>::ZERO;
+    let err = if neg == (zb > xy) { xy - result + zb } else { zb - result + xy };
+    if neg == (err < B::ZERO) {
+        ui += one;
+    } else {
+        ui -= one;
+    }
+
+    B::from_bits(ui).cast()
+}

src/math/generic/mod.rs

@@ -1,5 +1,7 @@
 mod copysign;
 mod fabs;
+mod fma;
 
 pub use copysign::copysign;
 pub use fabs::fabs;
+pub use fma::{fma, fma_big};

src/math/mod.rs

@@ -114,7 +114,7 @@ use self::rem_pio2::rem_pio2;
 use self::rem_pio2_large::rem_pio2_large;
 use self::rem_pio2f::rem_pio2f;
 #[allow(unused_imports)]
-use self::support::{CastFrom, CastInto, DInt, Float, HInt, Int, MinInt};
+use self::support::{CastFrom, CastInto, DInt, Float, HInt, Int, IntTy, MinInt};
 
 // Public modules
 mod acos;
@@ -336,19 +336,23 @@ cfg_if! {
     if #[cfg(f16_enabled)] {
         mod copysignf16;
         mod fabsf16;
+        mod fmaf16;
 
         pub use self::copysignf16::copysignf16;
         pub use self::fabsf16::fabsf16;
+        pub use self::fmaf16::fmaf16;
     }
 }
 
 cfg_if! {
     if #[cfg(f128_enabled)] {
         mod copysignf128;
         mod fabsf128;
+        mod fmaf128;
 
         pub use self::copysignf128::copysignf128;
         pub use self::fabsf128::fabsf128;
+        pub use self::fmaf128::fmaf128;
     }
 }
 

src/math/support/float_traits.rs

@@ -147,7 +147,6 @@ pub trait Float:
 }
 
 /// Access the associated `Int` type from a float (helper to avoid ambiguous associated types).
-#[allow(dead_code)]
 pub type IntTy<F> = <F as Float>::Int;
 
 macro_rules! float_impl {

src/math/support/int_traits.rs

@@ -382,3 +382,22 @@ cast_into!(u64);
 cast_into!(i64);
 cast_into!(u128);
 cast_into!(i128);
+
+cast_into!(f32; f64);
+cast_into!(f64; f32);
+
+cfg_if! {
+    if #[cfg(f16_enabled)] {
+        cast_into!(f16; f32, f64);
+        cast_into!(f32; f16);
+        cast_into!(f64; f16);
+    }
+}
+
+cfg_if! {
+    if #[cfg(f128_enabled)] {
+        cast_into!(f128; f32, f64);
+        cast_into!(f32; f128);
+        cast_into!(f64; f128);
+    }
+}