FMA implementation producing incorrect results

[ajtribick]

I've run into an issue with the FMA implementation when investigating creating a no_std version of my double-double arithmetic crate. It breaks the algorithm for multiplying two f64 values to produce a double-double, for which I am using algorithm 3 from Joldes et al. (2017): the high word of the result is the standard floating point multiplication, while the low word uses fused-multiply-add.

The output of the algorithm should be two non-overlapping words, I'm testing this with code similar to the following (uses rand 0.7):

#[cfg(test)]
mod tests {
    use rand::Rng;

    #[test]
    fn fma_test() {
        let mut rng = rand::thread_rng();
        let dist = rand::distributions::Uniform::<f64>::new_inclusive(-10000.0, 10000.0);

        for _ in 0..1000000 {
            let a = rng.sample(dist);
            let b = rng.sample(dist);
            let hi = a * b;
            let lo = libm::fma(a, b, -hi); // 1

            if !hi.is_normal() || !lo.is_normal() { continue; }
            assert!(
                libm::ilogb(hi) - libm::ilogb(lo) >= f64::MANTISSA_DIGITS as i32,
                "fma({:?}, {:?}, {:?})", a, b, -hi
            );
        }
    }
}

The above test fails with libm 0.2.1. If the line marked // 1 is replaced with

            let lo = a.mul_add(b, -hi);

then the test passes.

The equivalent code for f32 using libm::fmaf passes.

I'm not sure if this is the same as #200 but so far I've not seen segfaults.

[jethrogb]

This may be related to rust-lang/rust#85925 (comment)