# Why doesn’t Julia fma automatically?

**URL:** <https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982>\
**Category:** General Usage\
**Created:** [October 13, 2023, 3:31am UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982 "2023-10-13T03:31:53Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [October 16, 2023, 3:32pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/21 "2023-10-16T15:32:32Z")

</div>

That’s not how `fadd` and `fmul` work, but higher-level languages have more options other than treating \* and + exactly the same as those.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [October 16, 2023, 3:39pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/22 "2023-10-16T15:39:04Z")

</div>

Changing `a*b + c` to not give the exact answer spelled out in excruciating detail by IEEE 754 would absolutely be a massively breaking change. The ability to implement algorithms like accurate interval arithmetic, transcendental functions, high precision types like DoubleDouble, compensated arithmetic (Kahan summation), etc. all depend on following IEEE precisely. No one _implements_ these kinds of things in Python AI/ML DSLs, which is why they can get away with not following IEEE exactly. (Also, they generally only have vectorized array operations which are an entirely different situation, so this isn’t even a valid comparison in the first place. Accurate summation of arrays is very expensive and generally not done for that reason.) People do implement these algorithms in C, C++ and Julia, which is why non-IEEE “optimizations” are opt-in in all of those languages.

There isn’t really much point in continuing this discussion. It has been explained repeatedly why we cannot, should not, and will not do FMA automatically. You keep disregarding those explanations and insisting that we should anyway.

---

<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:** [October 16, 2023, 4:04pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/23 "2023-10-16T16:04:31Z")

</div>

> [@StefanKarpinski](#):
>
> No one _implements_ these kinds of things in Python AI/ML DSLs, which is why they can get away with not following IEEE exactly.

But they do implement in C? So why can C/C++ get away with it (pragmas?)? Would it not be possible to get FMA semantics by default, if in some cases it’s left out? Might need to be an exception for libraries…? To me, it just seems important to know what those cases should or could be.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [October 16, 2023, 4:10pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/24 "2023-10-16T16:10:42Z")

</div>

Where did this idea come from that C and C++ do FMA automatically? All major C/C++ compilers follow IEEE 754 unless you ask for `-funsafe-math` which is local to a compilation unit. Literally no language that can be used for implementing serious low level numerical functionality does FMA without explicitly opting into it.

---

<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:** [October 16, 2023, 4:14pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/25 "2023-10-16T16:14:33Z")

</div>

> Where did this idea come from that C and C++ do FMA automatically?

From Clang, I assume that means C and C++.

See my comment above, and where I’ve answered you on that before (I confirmed personally for Clang for RISC-V with Godbolt):

> [@X \* y + z does not automatically use FMA instruction](https://discourse.julialang.org/t/x-y-z-does-not-automatically-use-fma-instruction/100640/33):
>
> EDIT: That’s just not true in all cases. Clang does it by default depending on the arch, e.g. uses fmadd.d for RISC-V rv64gc clang 14.0.0 (changed from older versions; same change for 32-bit). It’s not only for -O1, also with -O0, or no options which I assume is the same. [It’s just that I overlooked when looking at boilerplate, the adds you see then are not for the floating point.] That’s not entirely true, depending on what you mean by you “asked for” it, and “by default”. I see I can get it…

> 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.

I believe Clang has an exeptions to its new default for x86… at least for now.

> [@StefanKarpinski](#):
>
> which is local to a compilation unit

Julia doesn’t have such, except well packages (precompiled by now), why I proposed a possible exception for them/libraries.

---

<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:** [October 16, 2023, 4:16pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/26 "2023-10-16T16:16:49Z")

</div>

> [@StefanKarpinski](#):
>
> Where did this idea come from that C and C++ do FMA automatically?

Actually, see the `-ffp-contract` option in GCC or Clang. Clang switched to enabling it by default relatively recently, while GCC always had it on-by-default, I suppose. I think that only works in a single expression, though, so presumably it doesn’t work across statements, which is why it wouldn’t cause too much harm.

---

<div class="post-metadata">

**Author:** ![adienes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adienes/32/37459_2.png) [@adienes](https://discourse.julialang.org/u/adienes)\
**Post date:** [October 16, 2023, 4:21pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/27 "2023-10-16T16:21:28Z")

</div>

> [@StefanKarpinski](#):
>
> Literally no language that can be used for implementing serious low level numerical functionality does FMA without explicitly opting into it.

I agree with you in spirit

however, I (think, could be wrong) that OCaml does automatic FMA on hardware that supports it

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 16, 2023, 4:31pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/28 "2023-10-16T16:31:51Z")

</div>

Yes, the [default floating point model](https://clang.llvm.org/docs/UsersManual.html#cmdoption-ffp-model) in clang also permits FMA [since 2020](https://reviews.llvm.org/D74436) as @Palli described. I’m also not fully convinced that it’s a great idea. The idea that the math you write in a programming language represents some abstract mathematical truth and “just does what you mean” is so very seductive but I really don’t think it’s meaningful in a way you can count on because there’s always some place that the rubber hits the road and you’re no longer abstract.

We’ve discussed this tension many times; this comment still represents my views:

> [@When "if (a\<b) x=1/(a-b)" divides by zero](https://discourse.julialang.org/t/when-if-a-b-x-1-a-b-divides-by-zero/7154/5):
>
> If a sufficiently advanced compiler is indistinguishable from an adversary, then giving the compiler access to `-ffast-math` is gifting that enemy nukes.

---

<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:** [October 16, 2023, 4:40pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/29 "2023-10-16T16:40:39Z")

</div>

> [@mbauman](#):
>
> I’m also not fully convinced that it’s a great idea. The idea that the math you write in a programming language represents some abstract mathematical truth and “just does what you mean” is so very seductive

I’m not saying I’m convinced either, just that there are bad results with doing _either_. I just know other people were convinced enough to change the default in Clang.

When I ask for a\*b + c, then I didn’t specifically ask for rounding twice or even once, _explicitly_. With integers (or rationals) that does NOT happen. Julia wants to be generic, and rounding needs to happen with floating point (unless you widen…). It seems very plausible that rounding only once is the new _better default_.

What Julia could to is define new types. Float64\_no\_FMA etc. and Float64 (changed to mean float allowing FMA, one rounding). Would that be an idea people like? The representation would be the same, so a “conversion” between a “noop”. I so like type systems and Julia to explore such options.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 16, 2023, 4:51pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/30 "2023-10-16T16:51:04Z")

</div>

> [@nsajko](#):
>
> > [@StefanKarpinski](#):
> >
> > Where did this idea come from that C and C++ do FMA automatically?
> 
> Actually, see the `-ffp-contract` option in GCC or Clang. Clang switched to enabling it by default relatively recently

Apparently, this ultimately stems [from C99](https://port70.net/%7Ensz/c/c99/n1256.html#7.12.2), which added an `FP_CONTRACT` pragma that you can use to specify whether automatic FMA is allowed:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/e/1ee67dd8d8b480b8ddc93d61ce02b1328f2b8482.png)  
and the standard explicitly says (in section 7.12.2) that **The default state (‘‘on’’ or ‘‘off’’) for the pragma is implementation-defined**, which gave compilers explicit permission to fuse by default (but technically only within expressions that appear explicitly in the source code, not across inlined or merged expressions, I think).

For a long time the default was “off” in gcc unless you used `-funsafe-math-optimizations` or `-ffast-math`, but apparently it was changed [around 2010](https://gcc.gnu.org/git/?p=gcc.git;a=commit;h=1694907238eb106bf7ac0e4eaedaa77bc7719b6d) to default to `-ffp-contract=fast` (allowing automatic FMA on supporting hardware); however, the `FP_CONTRACT` pragma to disable it [isn’t supported](https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845) and this was “fixed” by switching to [`-fp-contract=off` in ISO C mode](https://gcc.gnu.org/git/?p=gcc.git&a=commit;h=6dbe09581a1349498f6b4546f68cb9fd3205e120). `clang` documented a change in the default to `-ffp-contract=on` [in clang 14](https://releases.llvm.org/14.0.0/tools/clang/docs/ReleaseNotes.html), but apparently it didn’t [actually change](https://github.com/llvm/llvm-project/commit/48ad446a0fb2c9b98cb7047e4daf8a84c29cef8f) on many platforms until 2021, and apparently they [still don’t fully support](https://groups.google.com/g/llvm-dev/c/Ge_c3hp6AzE) the pragma to disable it. In MSVC the [default is `off`](https://learn.microsoft.com/en-us/cpp/preprocessor/fp-contract?view=msvc-170) as of Visual Studio 2022 version 17.0 (but was `on` by default in previous versions). Intel `icc` enables [`-fp-contract=fast` by default](https://www.intel.com/content/www/us/en/docs/dpcpp-cpp-compiler/developer-guide-reference/2023-2/floating-point-optimizations.html), but they also [support the pragma](https://www.intel.com/content/www/us/en/docs/cpp-compiler/developer-guide-reference/2021-8/intel-supported-pragma-reference.html). So, it’s a mess right now and you can’t reliably indicate in the source code itself (as opposed to a compiler flag) that FMA contraction is not desired, thanks to inconsistent support for `#pragma STDC FP_CONTRACT OFF`.

One thing to keep in mind is that, even though this is a potentially breaking change, it was a change made with little fanfare by multiple compilers (albeit explicitly allowed by the language standard). So, there is some precedent if one wanted to make a similar change in Julia 1.x (presumably adding something like a `@strict_math` macro for code to opt out).

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [October 16, 2023, 5:12pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/31 "2023-10-16T17:12:47Z")

</div>

> [@StefanKarpinski](#):
>
> The ability to implement algorithms like accurate interval arithmetic, transcendental functions, high precision types like DoubleDouble, compensated arithmetic (Kahan summation), etc. all depend on following IEEE precisely.

I’m well aware that these algorithms exist. I think that if we changed FMA behavior, it would make sense to increment the major version to make sure the code that relies on this got the memo. I just think it’s a bad idea to make a promise that very few (if any) other languages do. As noted above, C/C++ does not make this promise, making Julia in this respect _a lower-level language than C/C++_. I also think it’s very bad practice to rely on this kind of exact order-of-operations without explicitly specifying it in your code, either with a macro, by defining a separate `exact_plus` function for this, or by using parentheses.

I’d like to note that I’m not disregarding these explanations. I fully understand them, and I recognize there are some costs to this optimization in the form of making these high-precision addition algorithms slightly clunkier to express. I just don’t find these explanations convincing. They argue that [someone, somewhere,](https://xkcd.com/1172/) would theoretically pay a nonzero cost for this optimization, not that those costs aren’t worth it. Addition algorithms like Kahan summation are a tiny sliver of code in the Julia ecosystem that could easily be fixed by slapping an `@no_fma` macro at the top of the functions that implement them, instead of requiring the other 99.9% of Julia code to remember to either remember to include `@fma` or accept the accuracy and performance hit.

The way defaults are supposed to work is that they’re supposed to work for most people, with a way to opt-out for the few cases where it doesn’t work, not the other way around.

---

<div class="post-metadata">

**Author:** ![adienes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adienes/32/37459_2.png) [@adienes](https://discourse.julialang.org/u/adienes)\
**Post date:** [October 16, 2023, 5:29pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/32 "2023-10-16T17:29:12Z")

</div>

> [@ParadaCarleton](#):
>
> tiny sliver of code in the Julia ecosystem that could easily be fixed

this might land with more gravity if you are volunteering yourself to contribute such an easy fix on that tiny corner of the ecosystem

---

<div class="post-metadata">

**Author:** ![vchuravy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vchuravy/32/8_2.png) [@vchuravy](https://discourse.julialang.org/u/vchuravy)\
**Post date:** [October 16, 2023, 5:36pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/33 "2023-10-16T17:36:48Z")

</div>

One of the packages the I maintain (KernelAbstractions.jl) used to do fp contract by default (with the argument that it improved GPU performance and most GPU compilers do fp contract by default).

I received multiple real bug reports were this caused surprising and hard to understand behavior. The funniest was one were the same expression was evaluated to two different results within the same function.

I am a strong advocate for global behavior that is safe and consistent, with the local ability to opt-in to more aggressive behavior.

A real step forward would be to expose a much more fine-grained set of options for the `fastmath` macro so that you could choose to only get `:contract` and `:reassoc`

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 16, 2023, 5:47pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/34 "2023-10-16T17:47:14Z")

</div>

Here’s [a fun real world example](https://discourse.llvm.org/t/fp-contraction-fma-on-by-default/64975/4) from Clang’s world: using FMA means that `x1*y2 - x2*y1` is _sometimes_ nonzero when `x1==x2` and `y1==y2`.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2023, 5:55pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/35 "2023-10-16T17:55:19Z")

</div>

> **[GitHub - SciML/MuladdMacro.jl: This package contains a macro for converting...](https://github.com/SciML/MuladdMacro.jl)**
>
> This package contains a macro for converting expressions to use muladd calls and fused-multiply-add (FMA) operations for high-performance in the SciML scientific machine learning ecosystem - GitHub...

This macro makes it very easy to apply `@muladd` at a function level. That’s syntactic, understood, and safe. In many cases, that’s going to be a much simpler option to reason about.

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [October 16, 2023, 7:39pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/36 "2023-10-16T19:39:51Z")

</div>

> [@adienes](#):
>
> this might land with more gravity if you are volunteering yourself to contribute such an easy fix on that tiny corner of the ecosystem

Sure, I’d happily make a PR if the Julia devs announced a plan to optimize FMA automatically.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [October 16, 2023, 8:33pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/37 "2023-10-16T20:33:39Z")

</div>

And that one causes `a*conj(a)` to potentially have nonzero imaginary part if complex multiplication uses FMA. Another example is that `b^2 - 4 a c` will be nonnegative when b^2 \geq 4 a c if the obvious operations are used. However `fma(-4*a, c, b^2)` can still be negative, which can be inconvenient when using the quadratic formula. It’s not just problems with compensated summation or Double-Doubles.

It’s reassuring to see some resistance to the idea.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [October 16, 2023, 8:41pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/38 "2023-10-16T20:41:13Z")

</div>

Yeah, that’s exactly it. I mean I get it, this sure feels like it should be such an obvious win — better performance and more accurate calculations, what’s not to like!? I _want_ this to work.

But it’s _not_ uniformly more accurate. Getting “more accuracy” sporadically in unexpected places can lead to very baffling — _less accurate_ and _less algebraic_ — behaviors once you go just one step beyond the trivial `a*b + c`.

---

<div class="post-metadata">

**Author:** ![ParadaCarleton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paradacarleton/32/20005_2.png) [@ParadaCarleton](https://discourse.julialang.org/u/ParadaCarleton)\
**Post date:** [October 16, 2023, 11:58pm UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/39 "2023-10-16T23:58:06Z")

</div>

Thank you for your contribution 🙂 I think your post is the first one to bring up an actual case of these kinds of bugs causing problems in real life. I’m not sure if there might be something unique to the usual uses of KernelAbstractions.jl here, but if it caused problems there, that suggests it might not be a good choice for Julia as a whole.

Were the bugs all caused by people writing things like `a*x + b*y`, and not simple cases of `a*x+b`? If so, would it make sense to just stick to the basic/case easy win `a*x + b = fma(a, x, b)`, but just leave `a*x + b*y` alone?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 17, 2023, 12:13am UTC](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982/40 "2023-10-17T00:13:53Z")

</div>

the problem is that compilers don’t work one line at a time. You basically only have 2 options: allow re-association and contraction or prevent them.

[Previous page](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982.md?page=1)

[Next page](https://discourse.julialang.org/t/why-doesn-t-julia-fma-automatically/104982.md?page=3)