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

General Usage
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!