Setting matmul to fast makes all computations run in low precision

[Joel-Dahne]

Toggling IntervalArithmetic.configure(matmul = :fast) seems to make it so that matrix multiplication for BigFloat only runs in Float64 precision.

julia> using IntervalArithmetic

julia> A = sin.(interval.(rand(BigFloat, 4, 4))) # Apply sin to get non-zero radius

4×4 Matrix{Interval{BigFloat}}:
 [0.089218, 0.0892181]₂₅₆_com  [0.578398, 0.578399]₂₅₆_com  …  [0.732301, 0.732302]₂₅₆_com
 [0.50122, 0.501222]₂₅₆_com    [0.374433, 0.374434]₂₅₆_com     [0.677579, 0.67758]₂₅₆_com
 [0.499401, 0.499402]₂₅₆_com   [0.519317, 0.519319]₂₅₆_com     [0.52119, 0.521191]₂₅₆_com
 [0.470956, 0.470957]₂₅₆_com   [0.722737, 0.722738]₂₅₆_com     [0.245764, 0.245765]₂₅₆_com

julia> IntervalArithmetic.configure(matmul = :slow)
(Float64, :set_based, :correct, :fast, :slow)

julia> B = A * A
4×4 Matrix{Interval{BigFloat}}:
 [1.04141, 1.04142]₂₅₆_com    [1.212, 1.21201]₂₅₆_com     …  [1.05327, 1.05328]₂₅₆_com
 [0.835392, 0.835393]₂₅₆_com  [1.21502, 1.21504]₂₅₆_com      [1.08355, 1.08356]₂₅₆_com
 [0.964768, 0.96477]₂₅₆_com   [1.29097, 1.29098]₂₅₆_com      [1.27822, 1.27823]₂₅₆_com
 [0.601367, 0.601368]₂₅₆_com  [0.805239, 0.80524]₂₅₆_com     [0.979897, 0.979898]₂₅₆_com

julia> radius.(B)
4×4 Matrix{BigFloat}:
 3.45447e-77  5.1817e-77   5.1817e-77   3.45447e-77
 2.59085e-77  3.45447e-77  5.1817e-77   3.45447e-77
 2.59085e-77  3.45447e-77  5.1817e-77   5.1817e-77
 1.72723e-77  2.59085e-77  2.59085e-77  3.45447e-77

julia> IntervalArithmetic.configure(matmul = :fast)
(Float64, :set_based, :correct, :fast, :fast)

julia> B = A * A
4×4 Matrix{Interval{BigFloat}}:
 [1.04141, 1.04142]₂₅₆_com    [1.212, 1.21201]₂₅₆_com     …  [1.05327, 1.05328]₂₅₆_com
 [0.835392, 0.835393]₂₅₆_com  [1.21502, 1.21504]₂₅₆_com      [1.08355, 1.08356]₂₅₆_com
 [0.964768, 0.96477]₂₅₆_com   [1.29097, 1.29098]₂₅₆_com      [1.27822, 1.27823]₂₅₆_com
 [0.601367, 0.601368]₂₅₆_com  [0.805239, 0.80524]₂₅₆_com     [0.979897, 0.979898]₂₅₆_com

julia> radius.(B)
4×4 Matrix{BigFloat}:
 2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16
 1.11022e-16  2.22045e-16  3.33067e-16  2.22045e-16
 1.66533e-16  2.22045e-16  3.33067e-16  3.33067e-16
 1.11022e-16  1.66533e-16  1.66533e-16  2.22045e-16

I don't know if this is intended (it is certainly much faster), but it should probably be documented in that case.

On a related note, is there a specific reason for the implementing a custom version of matrix multiplication for the :slow path rather than just using the default? A quick benchmark seems to indicate that the default version is slightly faster.

julia> using IntervalArithmetic, BenchmarkTools

julia> IntervalArithmetic.configure(matmul = :slow)

(Float64, :set_based, :correct, :fast, :slow)

julia> M = sin.(interval(rand(100, 100))); # Satisfy m + n + p > 256

julia> res = similar(M);

julia> @benchmark mul!($res, $M, $M)
BenchmarkTools.Trial: 78 samples with 1 evaluation per sample.
 Range (min … max):  64.304 ms … 64.798 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     64.379 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   64.407 ms ± 99.221 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

    ▄▆▃▁▃  █▁▁                                                 
  ▇▆█████▆▆███▄▇▄▁▆▄▄▄▇▄▁▄▁▆▁▁▄▁▁▁▄▁▇▁▄▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▁
  64.3 ms         Histogram: frequency by time        64.8 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark invoke($mul!, $(Tuple{Any,Any,Any}), $res, $M, $M) # Force use of default implementation
BenchmarkTools.Trial: 88 samples with 1 evaluation per sample.
 Range (min … max):  56.995 ms … 57.341 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     57.111 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   57.122 ms ± 83.804 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

       ▁▃▁▃ ▃▁  ▃  ▃   ▃  █▁           ▃     ▁                 
  ▇▇▄▄▁████▁██▇▄█▄▁█▇▄▇█▁▄██▄▄▁▁▇▇▇▁▁▄▄█▇▁▁▇▁█▄▄▄▁▄▄▁▁▄▁▄▁▁▁▄ ▁
  57 ms           Histogram: frequency by time        57.3 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

All of this is on the currently latest commit.