faster hypot for Float32 and Float16

[Moelf]

Before:

julia> @benchmark hypot(x,y) setup=(begin x,y = rand(Float32, 2) end) evals=10000
BenchmarkTools.Trial: 10000 samples with 10000 evaluations.
 Range (min … max):  3.051 ns … 8.185 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     5.700 ns             ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.616 ns ± 0.208 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                                          ▁▇▇  ▂▇█▅▁▂▂▂     ▂
  ▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁████▇█████████████ █
  3.05 ns     Histogram: log(frequency) by time     6.27 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark hypot(x,y) setup=(begin x,y = rand(Float16, 2) end) evals=10000
BenchmarkTools.Trial: 10000 samples with 10000 evaluations.
 Range (min … max):  3.133 ns … 10.518 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     6.265 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   6.229 ns ±  0.280 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                                          ▇▆▃▆▁▁█▄▇▂▄▂▁▁▁    ▂
  ▆▅▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄█████████████████▇ █
  3.13 ns      Histogram: log(frequency) by time     7.06 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

After:

julia> @benchmark _hypot(x,y) setup=(begin x,y = rand(Float32, 2) end) evals=10000
BenchmarkTools.Trial: 10000 samples with 10000 evaluations.
 Range (min … max):  2.161 ns … 5.114 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.188 ns             ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.197 ns ± 0.082 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▂ ▂ ▇█ ▂                                                  ▁
  █▁█▃██▁█▄▃█▃▁▄▃▁▁▁▁▅▃▁▃▁▃▃▄▄▁▁▁▃▁▃▁▄▄▄▃▅▃▄▄▄▄▄▅▃▄▄▆▅▄▄▄▁▅ █
  2.16 ns     Histogram: log(frequency) by time     2.47 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark _hypot(x,y) setup=(begin x,y = rand(Float16, 2) end) evals=10000
BenchmarkTools.Trial: 10000 samples with 10000 evaluations.
 Range (min … max):  2.547 ns … 3.631 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.620 ns             ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.630 ns ± 0.070 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▁ ▆ ▄ █ ▆   ▃                                          ▂
  ▇▇▃█▁█▆█▃█▄█▁▇▄██▅▄▅▄▄▃▄▅▅▅▅▄▄▅▅▅▅▆▃▆▅▆▇▇▇▆▆▆▃▅▃▅▅▅▅▄▅▅▄▆ █
  2.55 ns     Histogram: log(frequency) by time        3 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[Moelf]

close #36353

now, I'm not sure the old _hypot is needed for real number since 0, NaN, Inf was already doing the correct thing but I guess we will see from test results.

[oscardssmith]

I think if you write @fastmath sqrt(x*x+y*y), it will let x^2+y^2 turn into fma(x, x, y*y) which should be faster.

[brett-mahar]

@oscardssmith there was talk of removing fastmath from Julia, apparently its a bit shonky:

#36246

[Moelf]

I think FMA makes it much slower on machine without FMA? don't remember the conclusion

[oscardssmith]

@brett-mahar while we might eventually get rid of it, we won't until there is a good replacement. The reason we still have @fastmath is that you need it to get fast code in a bunch of circumstances.
@Moelf That's why I recommended using @fastmath rather than explicitly writing an fma in. This way, LLVM will be able to use fma on the CPUs where there is hardware support, and just do x*x+y*y on cpus that don't have fma in hardware.

[Moelf]

@fastmath doesn't exist at this stage

julia/base/Base.jl

Line 334 in 8812c5c

include("fastmath.jl")

[dkarrasch]

Seems like you need an explicit inf handling. Since the Float16 case didn't fail in the tests, maybe we should add it? Is Float16 tested at all?

[KristofferC]

This way, LLVM will be able to use fma on the CPUs where there is hardware support, and just do xx+yy on cpus that don't have fma in hardware.

Isn't that just muladd?

I don't think we should add a bunch of @fastmath in these types of functions because it is underspecified what transformations @fastmath actually allows. As it is right now, @fastmath is not used anywhere in Base.

[oscardssmith]

Oh yeah, muladd is probably what we want

[dkarrasch]

Shall we have tests for the Float16 case, including infs and nans? Otherwise this is good to go, if there are no accuracy concerns left.

[KristofferC]

Saw this on Wikipedia on FMA (https://en.wikipedia.org/wiki/Multiply–accumulate_operation#Fused_multiply–add):

Fused multiply–add can usually be relied on to give more accurate results. However, William Kahan has pointed out that it can give problems if used unthinkingly. If x^2 − y^2 is evaluated as ((x × x) − y × y) (following Kahan's suggested notation in which redundant parentheses direct the compiler to round the (x × x) term first) using fused multiply–add, then the result may be negative even when x = y due to the first multiplication discarding low significance bits. This could then lead to an error if, for instance, the square root of the result is then evaluated.

Here we do that except x^2 + y^2. I guess that's fine then?

[oscardssmith]

The reason this works is because we are doing a higher precision fma. the multiplication of 2 Float32 has 46 bits of precision or less, so a Float64 can store y*y exactly.

[ViralBShah]

Bump. What do we need to do to get this in?

[oscardssmith]

Not very much.