Automatically choose a correct fma implementation at run time

base/float.jl

@@ -198,8 +198,6 @@ widen(::Type{Float32}) = Float64
 /(x::Float32, y::Float32) = box(Float32,div_float(unbox(Float32,x),unbox(Float32,y)))
 /(x::Float64, y::Float64) = box(Float64,div_float(unbox(Float64,x),unbox(Float64,y)))
 
-fma(x::Float32, y::Float32, z::Float32) = box(Float32,fma_float(unbox(Float32,x),unbox(Float32,y),unbox(Float32,z)))
-fma(x::Float64, y::Float64, z::Float64) = box(Float64,fma_float(unbox(Float64,x),unbox(Float64,y),unbox(Float64,z)))
 muladd(x::Float32, y::Float32, z::Float32) = box(Float32,muladd_float(unbox(Float32,x),unbox(Float32,y),unbox(Float32,z)))
 muladd(x::Float64, y::Float64, z::Float64) = box(Float64,muladd_float(unbox(Float64,x),unbox(Float64,y),unbox(Float64,z)))
 
@@ -387,6 +385,31 @@ end
     eps() = eps(Float64)
 end
 
+# fused multiply-add
+fma_libm(x::Float32, y::Float32, z::Float32) =
+    ccall(("fmaf", libm_name), Float32, (Float32,Float32,Float32), x, y, z)
+fma_libm(x::Float64, y::Float64, z::Float64) =
+    ccall(("fma", libm_name), Float64, (Float64,Float64,Float64), x, y, z)
+fma_llvm(x::Float32, y::Float32, z::Float32) =
+    box(Float32,fma_float(unbox(Float32,x),unbox(Float32,y),unbox(Float32,z)))
+fma_llvm(x::Float64, y::Float64, z::Float64) =
+    box(Float64,fma_float(unbox(Float64,x),unbox(Float64,y),unbox(Float64,z)))
+# Disable LLVM's fma if it is incorrect, e.g. because LLVM falls back
+# onto a broken system libm; if so, use openlibm's fma instead
+# 1.0000305f0 = 1 + 1/2^15
+if fma_llvm(1.0000305f0, 1.0000305f0, -1.0f0) == 6.103609f-5
+    fma(x::Float32, y::Float32, z::Float32) = fma_llvm(x,y,z)
+else
+    fma(x::Float32, y::Float32, z::Float32) = fma_libm(x,y,z)
+end
+# 1.0000000009313226 = 1 + 1/2^30
+if (fma_llvm(1.0000000009313226, 1.0000000009313226, -1.0) ==
+    1.8626451500983188e-9)
+    fma(x::Float64, y::Float64, z::Float64) = fma_llvm(x,y,z)
+else
+    fma(x::Float64, y::Float64, z::Float64) = fma_libm(x,y,z)
+end
+
 ## byte order swaps for arbitrary-endianness serialization/deserialization ##
 bswap(x::Float32) = box(Float32,bswap_int(unbox(Float32,x)))
 bswap(x::Float64) = box(Float64,bswap_int(unbox(Float64,x)))