Why is `muladd` not consistent even though it has the effect?

General Usage
1

This inconsistency also seems to infect functions that call muladd like evalpoly. This means that currently PhiloxRNG.jl’s randn_f32 which calls evalpoly is also impacted and returns different random numbers even if the seed is exactly the same and the function is running on the exact same computer. Is this a bug in the compiler, or should I just replace evalpoly and such with versions that use fma directly?

Here is an example:

using PhiloxRNG
@inline g(k) = randn_f32(UInt64(k), UInt64(0), UInt64(0))

@noinline full(k)       = g(k)
@noinline first_only(k) = g(k)[1]

k = 93

full(k)[1] - first_only(k)
# -1.1920929f-7

Base.infer_effects(randn_f32)
# (+c,+e,+n,+t,+s,+m,+u,+o,+r)

julia> versioninfo()
Julia Version 1.13.1
Commit 96ca370cf0e (2026-09-25 19:34 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 7 9800X3D 8-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver5)
  GC: Built with stock GC
Threads: 4 default, 1 interactive, 4 GC (on 16 virtual cores)
2

The muladd chooses either the fused or split implementation depending on which one is faster.

If you can tolerate some precision issues, I think avoiding explicit fma would allow your code to be faster.

3

Interestingly in most cases I have tested explicit fma is slightly faster.

4

As far as I know fma is supposed to be faster, unless it doesn’t exist on the processor and has to be emulated.

5

Yes, that is why it is so strange that the Julia compiler sometimes decides to not convert it to fma depending on the context of the surrounding code, not just based on the processor.

6

If this is true, it seems more likely to be an LLVM bug than a Julia thing. Can you construct a minimal example where this occurs?

Something that’s come up in the past for me is that common-subexpression elimination can sometimes mistakenly eliminate fma instructions, as in e.g. this 2016 gcc bug. I tried something simple along these lines, but it seemed to work in Julia 1.12.6 and 1.13.3 on Apple silicon: for f(x,y,z) = muladd(x,y,z) * muladd(x,y,-z), doing @code_native f(1.0,2.0,3.0) correctly emits fmadd and fnmsub instructions (rather than trying to hoist the x*y).

7

I agree it’s surprising at first that randn_f32 gets a +c but isn’t actually :consistent in your usage here… but I’m actually not sure it’s a real compiler problem because it’s inlined! And given that, I’m not surprised to see that there’d exist some functions that optimize it differently (particularly if they throw away some bits of the computation!).

8

Here’s a more minimal example, stemming from a favorite old post of mine:

julia> discriminant(a,b,c) = muladd(b, b, -4*a*c)
discriminant (generic function with 1 method)

julia> g1(x,y,z) = discriminant(x,y,z), y^2.0
g1 (generic function with 1 method)

julia> g2(x,y,z) = discriminant(x,y,z), y^2
g2 (generic function with 1 method)

julia> g1(.2, .1, .0125)
(-8.326672684688674e-19, 0.010000000000000002)

julia> g2(.2, .1, .0125)
(8.326672684688674e-19, 0.010000000000000002)

julia> Base.infer_effects(discriminant, (Float64, Float64, Float64))
(+c,+e,+n,+t,+s,+m,+u,+o,+r)

This then has potential real-world impacts in a quadratic-equation root-solver:

julia> @inline function roots(a,b,c)
           if discriminant(a,b,c) < 0
               "no real roots"
           else
               roots = [(-b ± sqrt(discriminant(a,b,c)))/(2*a) for (±) in (+, -)]
               "root at $(unique(roots))"
           end
       end
roots (generic function with 1 method)

julia> roots(.2, .1, .0125)
"no real roots"

julia> f1(a,b,c) = roots(a,b,c), b^2.0
f1 (generic function with 1 method)

julia> f1(.2, .1, .0125)
("no real roots", 0.010000000000000002)

julia> f2(a,b,c) = roots(a,b,c), b*b
f2 (generic function with 1 method)

julia> f2(.2, .1, .0125)
ERROR: DomainError with -8.326672684688674e-19:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
[...]


julia> versioninfo()
Julia Version 1.13.3
Commit 67ac47b9cd8 (2026-10-25 15:28 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin27.0.0)
  CPU: 8 × Apple M1 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, apple-m1)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)

discriminant is +c but is giving different values, even inside the same method body!

9

It does look like marking muladd as !c is a long-standing todo item… I just asked Claude to run some benchmarks to see if uncommenting it from the list there is now safe to do.

10

I’ve updated PhiloxRNG to use explicit fma instead of muladd or regular adds and multiplies. It turns out that in addition to fusing muladd, CUDA.jl will sometimes fuse regular adds and multiplies, and that was a big part of why the random numbers were occasionally different on different platforms.