Skip to content

Commit

Permalink
Base: twiceprecision: optimize mul12
Browse files Browse the repository at this point in the history
It is well known and obvious that the algorithm behind `Base.mul12`
(sometimes known as "Fast2Mult") doesn't require a "Fast2Sum" (known in
the Julia codebase as `Base.canonicalize2`) call at the end, so remove
it. This follows from the fact that IEEE-754 floating-point
multiplication is required to be well-rounded, so Fast2Sum can't change
the result.

Reference, for example, the beginning of
https://doi.org/10.1145/3121432 by Joldes, Muller, Popescu.

Furthermore, `Base.Math.two_mul` already exists, so use it as a kernel
function. This required adding a general fallback method for
`Base.Math.two_mul`, required, for example, for `BigFloat`.

Also removed the `iszero` check that is now obviously not necessary.
  • Loading branch information
nsajko committed May 4, 2023
1 parent 827d34a commit 195bae7
Show file tree
Hide file tree
Showing 2 changed files with 7 additions and 2 deletions.
5 changes: 5 additions & 0 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,11 @@ end

# non-type specific math functions

function two_mul(x::T, y::T) where {T<:Number}
xy = x*y
xy, fma(x, y, -xy)
end

@assume_effects :consistent @inline function two_mul(x::Float64, y::Float64)
if Core.Intrinsics.have_fma(Float64)
xy = x*y
Expand Down
4 changes: 2 additions & 2 deletions base/twiceprecision.jl
Original file line number Diff line number Diff line change
Expand Up @@ -112,8 +112,8 @@ julia> Float64(hi) + Float64(lo)
```
"""
function mul12(x::T, y::T) where {T<:AbstractFloat}
h = x * y
ifelse(iszero(h) | !isfinite(h), (h, h), canonicalize2(h, fma(x, y, -h)))
(h, l) = Base.Math.two_mul(x, y)
ifelse(!isfinite(h), (h, h), (h, l))
end
mul12(x::T, y::T) where {T} = (p = x * y; (p, zero(p)))
mul12(x, y) = mul12(promote(x, y)...)
Expand Down

0 comments on commit 195bae7

Please sign in to comment.