# X \* y + z does not automatically use FMA instruction

**URL:** <https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640>\
**Category:** Performance\
**Tags:** fast-math\
**Created:** [June 21, 2023, 7:07am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640 "2023-06-21T07:07:22Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![GLOBEX\_CORP](https://avatars.discourse-cdn.com/v4/letter/g/e9bcb4/32.png) [@GLOBEX\_CORP](https://discourse.julialang.org/u/GLOBEX_CORP)\
**Post date:** [June 21, 2023, 7:07am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/1 "2023-06-21T07:07:22Z")

</div>

I am new to Julia and was wondering why the code `x * y + z` does not use `fmadd` where all values are floats, but instead `fmul` followed by `fadd`. On the other hand, if I use the `@fastmath` decorator, the `@code_native` changes to use the `fmadd` instruction. My question is, if FMA is faster and more precise, why won’t it use `fmadd` to begin with on supported machines?

Machine: M1 apple silicon on 13.4 Ventura  
Julia: 1.9.1  
Optimization: -O3 and -O2

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [June 21, 2023, 8:33am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/2 "2023-06-21T08:33:30Z")

</div>

> [@GLOBEX\_CORP](#):
>
> if FMA is faster and more precise, why won’t it use `fmadd` to begin with on supported machines?

Because it gives a _different_ result. You can explicitly tell the compiler that you are fine with that with `@fastmath` or by using either the `muladd` or `fma` function. See also the recent discussion [How to enable vectorized fma instruction for multiply-add vectors?](https://discourse.julialang.org/t/how-to-enable-vectorized-fma-instruction-for-multiply-add-vectors/100219).

---

<div class="post-metadata">

**Author:** ![GLOBEX\_CORP](https://avatars.discourse-cdn.com/v4/letter/g/e9bcb4/32.png) [@GLOBEX\_CORP](https://discourse.julialang.org/u/GLOBEX_CORP)\
**Post date:** [June 21, 2023, 8:54am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/3 "2023-06-21T08:54:10Z")

</div>

@GunnarFarneback Thank you for your reply and the link.

I was under the impression that because rounding only occurs at the end as opposed to at each step, `fma` is more precise. So because not all machines have native `fma` instructions, it is better to default to behavior that will be consistent between all machines?

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 21, 2023, 9:04am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/4 "2023-06-21T09:04:08Z")

</div>

No, I think the reasoning is that fma does not comply with the floating point standard IIUC

---

<div class="post-metadata">

**Author:** ![GLOBEX\_CORP](https://avatars.discourse-cdn.com/v4/letter/g/e9bcb4/32.png) [@GLOBEX\_CORP](https://discourse.julialang.org/u/GLOBEX_CORP)\
**Post date:** [June 21, 2023, 9:25am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/5 "2023-06-21T09:25:12Z")

</div>

> [@lrnv](#):
>
> floating point standard IIUC

@lrnv I’m not familiar with this standard, do you have a link?

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 21, 2023, 9:40am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/6 "2023-06-21T09:40:16Z")

</div>

This part of Julia docs details the behavior : [Integers and Floating-Point Numbers · The Julia Language](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/)

On this page, it is said that Julia floats comply with the IEEE 754 standard, with a link to this Wikipedia page : [IEEE 754-2008 revision - Wikipedia](https://en.wikipedia.org/wiki/IEEE_754-2008)

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [June 21, 2023, 11:06am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/7 "2023-06-21T11:06:31Z")

</div>

> [@lrnv](#):
>
> fma does not comply with the floating point standard

This is some misunderstanding. IEEE 754 _specifies_ FMA.

> [@GLOBEX\_CORP](#):
>
> I was under the impression that because rounding only occurs at the end as opposed to at each step, `fma` is more precise.

Yeah, FMA is more accurate, but the added accuracy may not be what the user wants. C/C++ compilers have flags to control this behavior (`-ffp-contract` in GCC and Clang).

---

<div class="post-metadata">

**Author:** ![GLOBEX\_CORP](https://avatars.discourse-cdn.com/v4/letter/g/e9bcb4/32.png) [@GLOBEX\_CORP](https://discourse.julialang.org/u/GLOBEX_CORP)\
**Post date:** [June 21, 2023, 11:11am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/8 "2023-06-21T11:11:11Z")

</div>

> [@nsajko](#):
>
> This is some misunderstanding. IEEE 754 _specifies_ FMA.

That corroborates what I had found while digging around a bit.

> [@nsajko](#):
>
> Yeah, FMA is more accurate, but the added accuracy may not be what the user wants. C/C++ compilers have flags to control this behavior (`-ffp-contract` in GCC and Clang).

I understand, thanks.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 21, 2023, 11:34am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/9 "2023-06-21T11:34:17Z")

</div>

More often than not, `@fastmath` will make code more accurate, e.g. by using `fma` instructions and multiple separate accumulators in reductions.

Different than what you explicitly wrote is the problem, not accuracy.

On average, `@fastmath` also makes code much less accurate, because of a small number of catastrophic cases when code relies on being executed explicitly as written.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 21, 2023, 11:39am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/10 "2023-06-21T11:39:06Z")

</div>

> [@nsajko](#):
>
> C/C++ compilers have flags to control this behavior (`-ffp-contract` in GCC and Clang).

Neither GCC nor Clang need explicit flags to use FMA instructions:

> **[Compiler Explorer - C++](https://godbolt.org/z/da3PGM3nn)**
>
> double fmadd(double x, double y, double z){
> return x\*y + z;
> }

I am personally in favor of automatic FMA, but enough people oppose it that I don’t expect it to happen.

---

<div class="post-metadata">

**Author:** ![GLOBEX\_CORP](https://avatars.discourse-cdn.com/v4/letter/g/e9bcb4/32.png) [@GLOBEX\_CORP](https://discourse.julialang.org/u/GLOBEX_CORP)\
**Post date:** [June 21, 2023, 11:48am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/11 "2023-06-21T11:48:28Z")

</div>

> [@Elrod](#):
>
> Neither GCC nor Clang need explicit flags to use FMA instructions:

Is that not because you specified `-march=` ? -O3 on its own is not enough to use FMA.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [June 21, 2023, 11:50am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/12 "2023-06-21T11:50:19Z")

</div>

> [@Elrod](#):
>
> Neither GCC nor Clang need explicit flags to use FMA instructions

I think Clang used to require setting `-ffp-contract` for this, however they changed the default a year or so ago. [Currently](https://clang.llvm.org/docs/UsersManual.html#cmdoption-ffp-contract) `on` is the default.

---

<div class="post-metadata">

**Author:** ![gbaraldi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gbaraldi/32/22101_2.png) [@gbaraldi](https://discourse.julialang.org/u/gbaraldi)\
**Post date:** [June 21, 2023, 11:59am UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/13 "2023-06-21T11:59:43Z")

</div>

We may want to do that on julia if we add an easy way to back off from it.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 21, 2023, 3:12pm UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/14 "2023-06-21T15:12:02Z")

</div>

Haswell CPUs have FMA. Older ones don’t.  
Code that can run on older CPUs thus can’t use it.

Julia starts with the equivalent of `march=native` by default.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 21, 2023, 3:13pm UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/15 "2023-06-21T15:13:14Z")

</div>

We already have `--math-mode=ieee`.  
But I think we should add a local overwrite, too.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [June 21, 2023, 4:57pm UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/16 "2023-06-21T16:57:27Z")

</div>

FMA is probably not used for (outdated) historical reasons. What Julia does currently, and IEEE standard allows, i.e. no FMA, can result in catastrophic cancellation (loss of all accuracy) that FMA would avoid:

> **[1. Introduction — Floating Point and IEEE 754 12.3 documentation](https://docs.nvidia.com/cuda/floating-point/index.html)**
>
> White paper covering the most common issues related to NVIDIA GPUs.

> Let’s consider an example to illustrate how the FMA operation works using decimal arithmetic first for clarity. […] the correct mathematical result is [not zero, but close, and zero result would be very bad]
> 
> Rounding the multiply and add separately yields a result that is off by 0.00064. The corresponding FMA computation is wrong by only 0.00004, and its result is closest to the correct mathematical answer. […]
> 
> Figure 1 shows CUDA C++ code and output corresponding to inputs A and B and operations from the example above. The code is executed on two different hardware platforms: an x86-class CPU using SSE in single precision, and an NVIDIA GPU with compute capability 2.0. At the time this paper is written (Spring 2011) there are no commercially available x86 CPUs which offer hardware FMA. Because of this, the computed result in single precision in SSE would be 0.

> [@GunnarFarneback](#):
>
> > [@GLOBEX\_CORP](#):
> >
> > if FMA is faster and more precise, why won’t it use `fmadd` to begin with on supported machines?
> 
> Because it gives a _different_ result.

That seems like a good reason, but that different FMA result is (usually) more accurate, sometime much better, especially with Float32 (that you want to use for speed), and Float16 I guess. On this, I see why FMA is the default in C compilers now:

> The C standard permits intermediate floating-point results within an expression to be computed with more precision than their type would normally allow. This permits operation fusing, and Clang takes advantage of this by default

Since it’s faster (and more accurate) and we want to compete with those other compilers I support changing the default. On non-FMA hardware it’s going to be way slower though (if we would insist on same result by emulation), that seems ok, as such hardware is increasingly rare, already obsolete in my view.

> [@lrnv](#):
>
> On this page, it is said that Julia floats comply with the IEEE 754 standard, with a link to this Wikipedia page : [IEEE 754-2008 revision - Wikipedia](https://en.wikipedia.org/wiki/IEEE_754-2008)

I’m not sure FMA by default contradicts IEEE. IEEE demands by now FMA availability, I believe (it’s mentioned in that 2008 revision, Clause 5: Operations, as new, though might be optional still) but not its use. If I’m wrong on using it usually more accurate, i.e. sometimes less (or catastrophically so?), then I also support:

> [@Elrod](#):
>
> We already have `--math-mode=ieee`.  
> But I think we should add a local overwrite, too.

You had the global --math-mode=fast option, but it’s dangerous (in Julia at least, so my PR to make it a no-op was merged; it’s the default I think with -O3 in other languages, not Julia, they get away with it). I’m actually not sure why that is, was it since it enabled something more than FMA? I’m assume just the FMA by default is mostly ok, why I support it by default, and likely ok with that possible local opt-out. The standard library (and all the current package ecosystem) could be made to work ok with the new default, fixed with an opt-out where needed. There’s some issue on this already on JuliaLang.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [June 21, 2023, 5:40pm UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/17 "2023-06-21T17:40:03Z")

</div>

> [@Palli](#):
>
> FMA is probably not used for (outdated) historical reasons.

No, Gunnar already gave the right reason.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 21, 2023, 7:26pm UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/18 "2023-06-21T19:26:32Z")

</div>

> [@Palli](#):
>
> I’m actually not sure why that is, was it since it enabled something more than FM,A?

Consider things like [Kahan summation algorithm - Wikipedia](https://en.wikipedia.org/wiki/Kahan_summation_algorithm)  
Translating their algorithm from pseudo-code into Julia:

```julia
julia> using AccurateArithmetic, StableRNGs

julia> function KahanSum(input)
           sum = zero(eltype(input)) # Prepare the accumulator.
           c = zero(eltype(input)) # A running compensation for lost low-order bits.
       
           for i = eachindex(input)
               y = input[i] - c # c is zero the first time around.
               t = sum + y # Alas, sum is big, y small, so low-order digits of y are lost.
               c = (t - sum) - y # (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
               sum = t # Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
           end # Next time around, the lost low part will be added to y in a fresh attempt.
       
           return sum
       end
KahanSum (generic function with 1 method)

julia> function KahanSumFast(input)
           sum = zero(eltype(input)) # Prepare the accumulator.
           c = zero(eltype(input)) # A running compensation for lost low-order bits.
       
           @fastmath for i = eachindex(input)
               y = input[i] - c # c is zero the first time around.
               t = sum + y # Alas, sum is big, y small, so low-order digits of y are lost.
               c = (t - sum) - y # (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
               sum = t # Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
           end # Next time around, the lost low part will be added to y in a fresh attempt.
       
           return sum
       end
KahanSumFast (generic function with 1 method)

julia> rng = StableRNG(1); x = 5000 |> N->randn(rng, N) .* exp.(10 .* randn(rng, N)) |> x->[x;-x;1.0] |> x->x[sortperm(rand(rng,length(x)))];

julia> foldl(+, x)
4.554186716906959

julia> sum(x)
1.875

julia> sum(big, x)
1.0

julia> KahanSum(x)
1.0893429687695757

julia> KahanSumFast(x)
0.985827341906959

julia> sum_kbn(x)
1.0

julia> sum_oro(x)
1.0

julia> function sum_fast(input)
           sum = zero(eltype(input))
           @fastmath for i = eachindex(input)
               @inbounds sum += input[i] end return sum endsum_fast (generic function with 1 method)julia> sum_fast(x)0.985827341906959

```

If floating point arithmetic were associative, we could transform

```julia
y = input[i] - c
t = sum + y
c = (t - sum) - y
sum = t

```

into

```julia
y = input[i] - c
t = sum + y
c = (sum + y - sum) - y = 0
sum = t 

```

into

```julia
sum += input[i] 

```

That is, fast math essentially lets the compiler delete the code that tries to track and compensate for accumulating floating point error.

However, because fast math allows the compiler to change the result, it also applies a transform – multiple separate accumulators – that massively decreases the floating point error. Hence `KahanSumFast` – while equivalent to the naive `@fastmath` loop – is way more accurate (and faster) than `foldl(+, x)`, and not actually much worse than kahan summation for the example above.  
[AccurateArithmetic.jl](https://github.com/JuliaMath/AccurateArithmetic.jl) is much better than the naive KahanSummation, combining the compensation of error with the separate parallel accumulators (making it of course also much faster than the naive Kahan Summation).

```julia
julia> @btime KahanSum($x)
  45.747 μs (0 allocations: 0 bytes)
1.0893429687695757

julia> @btime KahanSumFast($x) # not actually a Kahan sum at all!
  511.698 ns (0 allocations: 0 bytes)
0.985827341906959

julia> @btime sum_kbn($x) # kahan sum with parallel accumulators
  1.722 μs (0 allocations: 0 bytes)
1.0

julia> @btime sum_oro($x) # Ogita–Rump–Oishi with parallel accumulators
  1.505 μs (0 allocations: 0 bytes)
1.0

```

Thanks to the parallel accumulators, `AccurateArithmetic` is only around 3x slower than the straight `@fastmath` sum, instead of 90x slower like the naive kahan summation described in the Wikipedia article.

Anyway, to answer your question:  
There are error-free transforms like in this article. `@fastmath` can undo them, like I showed. In this particular example, it isn’t bad, because `@fastmath` also allows applying another transform that dramatically improves both accuracy and speed, which isn’t allowed without it (`@simd` also allows this).  
But in other cases, such as implementing functions such as `exp`, we lose a lot of accuracy with `@fastmath`, because it only undoes the compensation!

In the case of languages like C, C++, Fortran, etc, a global `@fastmath` isn’t as bad, because functions like `exp` are replaced with only slightly less accurate implementations that have already been compiled with their error compensations intact (they just have less of it than the more accurate variants). That is, the “fast” exp you’re calling from C with fast math was itself not compiled under fast math.

Similarly, in Julia, `@fastmath exp(x)` will call a faster and less accurate `exp` implementation, but that implementation wasn’t itself compiled under fast math.

Implementing functions like these is a very delicate process. The authors of these libraries (like @Oscar_Smith ) are very deliberate about floating point rounding and error, and would thus need strict IEEE adherence for their implementations.

---

<div class="post-metadata">

**Author:** ![arch.d.robison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arch.d.robison/32/17699_2.png) [@arch.d.robison](https://discourse.julialang.org/u/arch.d.robison)\
**Post date:** [June 22, 2023, 3:23pm UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/19 "2023-06-22T15:23:01Z")

</div>

Using FMA automatically can cause surprises. Consider the following:

```julia
function f(a,b,c)
    @assert a*b ≥ c
    return sqrt(a*b-c)
end

a = 1.0 + 0.5^27
b = 1.0 - 0.5^27
c = 1.0
f(a,b,c)

```

It works as expected. Now suppose the compiler automatically uses FMA inside f, like this:

```julia
function f(a,b,c)
    @assert a*b ≥ c
    return sqrt(fma(a,b,-c))
end

```

Now `f(a,b,c)` throws and exception, because `fma(a,b,-c)` is negative, despite the assertion `a*b ≥ c`.

---

<div class="post-metadata">

**Author:** ![GLOBEX\_CORP](https://avatars.discourse-cdn.com/v4/letter/g/e9bcb4/32.png) [@GLOBEX\_CORP](https://discourse.julialang.org/u/GLOBEX_CORP)\
**Post date:** [June 22, 2023, 3:45pm UTC](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/20 "2023-06-22T15:45:46Z")

</div>

Thanks for the example.

I guess my question is why does replacing fma with muladd work?

```julia
function f(a,b,c)
    @assert a*b ≥ c
    return sqrt(muladd(a,b,-c))
end

```

What is the difference between the two functions?

Also using @fastmath works.

```julia
function f(a,b,c)
    @assert a*b ≥ c
    return @fastmath sqrt(a*b-c)
end

```

[Next page](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640.md?page=2)