Optimizing sums of products (dot products)

General Usage Performance
, ,
18

OK, this was just too juicy to resist. Here’s a barely-contrived example (Julia v1.10.5, Apple M1):

julia> discriminant(a,b,c) = @fastmath b^2 - 4*a*c
discriminant (generic function with 1 method)

julia> 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"

That works just fine. Sure, you might miss an exact root here or there, but those are the edge of numerical stability anyhow. What you probably wouldn’t expect is how easy it is for the above answer to change (and error!) depending on what you happen to do outside that function!

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

Is that scary enough? That said, clang and gcc do this sort of thing by default and you don’t see many folks complaining of issues like this.