Base: make TwicePrecision more accurate

Update the arithmetic code to use algorithms published since
twiceprecision.jl got written.

Handle complex numbers.

Prevent some harmful promotions.

Some of the introduced extended precision arithmetic code will also be
used in `base/div.jl` to tackle issue #49450 with PR #49561.

Results:

The example in #33677 still gives the same result.

I wasn't able to evaluate #23497 because of how much Julia changed in
the meantime.

Proof that #26553 is fixed:
```julia-repl
julia> const TwicePrecision = Base.TwicePrecision
Base.TwicePrecision

julia> a = TwicePrecision{Float64}(0.42857142857142855, 2.3790493384824782e-17)
Base.TwicePrecision{Float64}(0.42857142857142855, 2.3790493384824782e-17)

julia> b = TwicePrecision{Float64}(3.4, 8.881784197001253e-17)
Base.TwicePrecision{Float64}(3.4, 8.881784197001253e-17)

julia> a_minus_b_inexact = Tuple(a - b)
(-2.9714285714285715, 1.0150610510858574e-16)

julia> BigFloat(a) == sum(map(BigFloat, Tuple(a)))
true

julia> BigFloat(b) == sum(map(BigFloat, Tuple(b)))
true

julia> a_minus_b = BigFloat(a) - BigFloat(b)
-2.971428571428571428571428571428577679589762353999797347402693647078382455095635

julia> Float64(a_minus_b) == first(a_minus_b_inexact)
true

julia> Float64(a_minus_b - Float64(a_minus_b)) == a_minus_b_inexact[begin + 1]
true
```

Fixes #26553

Updates #49589

Author: nsajko

base/complex.jl

@@ -1125,3 +1125,190 @@ function complex(A::AbstractArray{T}) where T
     end
     convert(AbstractArray{typeof(complex(zero(T)))}, A)
 end
+
+const ComplexFloat = Complex{<:AbstractFloat}
+
+## Multi-word floating-point for `Complex`
+
+complex_twiceprec(real::Tw, imag::Tw) where {Tw<:TwicePrecision{<:AbstractFloat}} =
+    TwicePrecision(
+        complex(real.hi, imag.hi),
+        complex(real.lo, imag.lo),
+    )
+
+real(c::TwicePrecision{<:ComplexFloat}) =
+    TwicePrecision(real(c.hi), real(c.lo))
+
+imag(c::TwicePrecision{<:ComplexFloat}) =
+    TwicePrecision(imag(c.hi), imag(c.lo))
+
+# +
+
+function +(x::TwicePrecision{<:AbstractFloat}, y::ComplexFloat)
+    y_re = real(y)
+    y_im = imag(y)
+    sum = x + y_re
+    TwicePrecision(complex(sum.hi, y_im), complex(sum.lo))
+end
+
+function +(
+    x::TwicePrecision{<:ComplexFloat},
+    y::Union{AbstractFloat,TwicePrecision{<:AbstractFloat}},
+)
+    x_re = real(x)
+    x_im = imag(x)
+    sum = x_re + y
+    complex_twiceprec(sum, x_im)
+end
+
++(x::TwicePrecision{T}, y::T) where {T<:ComplexFloat} =
+    TwicePrecision(mw_operate(+, Tuple(x), (y,))...)
+
++(x::TwicePrecision{<:AbstractFloat}, y::TwicePrecision{<:ComplexFloat}) =
+    y + x
+
++(x::Tw, y::Tw) where {Tw<:TwicePrecision{<:ComplexFloat}} =
+    TwicePrecision(mw_operate(+, Tuple(x), Tuple(y))...)
+
+# -
+
+-(x::TwicePrecision{<:AbstractFloat}, y::ComplexFloat) =
+    x + -y
+
+-(x::TwicePrecision{<:ComplexFloat}, y::Union{AbstractFloat,ComplexFloat}) =
+    x + -y
+
+-(
+    x::TwicePrecision{<:ComplexFloat},
+    y::TwicePrecision{<:Union{AbstractFloat,ComplexFloat}},
+) =
+    x + -y
+
+-(x::TwicePrecision{<:ComplexFloat}, y::TwicePrecision{<:AbstractFloat}) =
+    x + -y
+
+# *
+
+function *(x::TwicePrecision{<:AbstractFloat}, y::ComplexFloat)
+    y_re = real(y)
+    y_im = imag(y)
+    res_re = x * y_re
+    res_im = x * y_im
+    complex_twiceprec(res_re, res_im)
+end
+
+function *(
+    x::TwicePrecision{<:ComplexFloat},
+    y::Union{AbstractFloat,TwicePrecision{<:AbstractFloat}},
+)
+    x_re = real(x)
+    x_im = imag(x)
+    res_re = x_re * y
+    res_im = x_im * y
+    complex_twiceprec(res_re, res_im)
+end
+
+function *(x::Tw, y::Union{T,Tw}) where {T<:ComplexFloat,Tw<:TwicePrecision{T}}
+    # TODO: evaluate whether it makes sense to write custom code, with
+    # separate methods for `*(x::Tw, y::T)` and `*(x::Tw, y::Tw)`.
+
+    # TODO: try preventing some overflows and underflows?
+
+    # TODO: for `*(x::Tw, y::T)`: evaluate whether it makes sense to
+    # use Algorithm 3 from the paper "Accurate Complex Multiplication
+    # in Floating-Point Arithmetic", https://hal.science/hal-02001080v2
+    # The algorithm must first be adapted to return double-word results
+    # as explained in the same paper.
+
+    x_re = real(x)
+    x_im = imag(x)
+    y_re = real(y)
+    y_im = imag(y)
+    res_re = x_re*y_re - x_im*y_im
+    res_im = x_im*y_re + x_re*y_im
+    complex_twiceprec(res_re, res_im)
+end
+
+*(x::TwicePrecision{<:AbstractFloat}, y::TwicePrecision{<:ComplexFloat}) =
+    y * x
+
+# abs2
+
+function abs2(x::TwicePrecision{<:ComplexFloat})
+    # TODO: evaluate whether it makes sense to write custom code
+
+    # TODO: try preventing some overflows and underflows?
+
+    re = real(x)
+    im = imag(x)
+    re*re + im*im
+end
+
+# /
+
+function div_complex_twiceprec_impl(
+    x::Union{AbstractFloat,TwicePrecision{<:AbstractFloat}},
+    y,
+)
+    # TODO: use the inverse of abs2 instead for better performance?
+    a = abs2(y)
+
+    y_re = real(y)
+    y_im = imag(y)
+    res_re = x*y_re/a
+    res_im = -x*y_im/a
+    complex_twiceprec(res_re, res_im)
+end
+
+function div_complex_twiceprec_impl(
+    x::Union{ComplexFloat,TwicePrecision{<:ComplexFloat}},
+    y,
+)
+    # TODO: use the inverse of abs2 instead for better performance?
+    a = abs2(y)
+
+    x_re = real(x)
+    x_im = imag(x)
+    y_re = real(y)
+    y_im = imag(y)
+    res_re = (x_re*y_re + x_im*y_im) / a
+    res_im = (x_im*y_re - x_re*y_im) / a
+    complex_twiceprec(res_re, res_im)
+end
+
+/(
+    x::TwicePrecision{<:AbstractFloat},
+    y::Union{ComplexFloat,TwicePrecision{<:ComplexFloat}},
+) =
+    # TODO: evaluate whether it makes sense to write custom code
+    # TODO: try preventing some overflows and underflows?
+    div_complex_twiceprec_impl(x, y)
+
+/(
+    x::AbstractFloat,
+    y::TwicePrecision{<:ComplexFloat},
+) =
+    # TODO: evaluate whether it makes sense to write custom code
+    # TODO: try preventing some overflows and underflows?
+    div_complex_twiceprec_impl(x, y)
+
+function /(
+    x::TwicePrecision{<:ComplexFloat},
+    y::Union{AbstractFloat,TwicePrecision{<:AbstractFloat}},
+)
+    x_re = real(x)
+    x_im = imag(x)
+    res_re = x_re / y
+    res_im = x_im / y
+    complex_twiceprec(res_re, res_im)
+end
+
+/(x::Tw, y::Union{T,Tw}) where {T<:ComplexFloat,Tw<:TwicePrecision{T}} =
+    # TODO: evaluate whether it makes sense to write custom code
+    # TODO: try preventing some overflows and underflows?
+    div_complex_twiceprec_impl(x, y)
+
+/(x::ComplexFloat, y::TwicePrecision{<:ComplexFloat}) =
+    # TODO: evaluate whether it makes sense to write custom code
+    # TODO: try preventing some overflows and underflows?
+    div_complex_twiceprec_impl(x, y)

base/mpfr.jl

@@ -986,6 +986,11 @@ isfinite(x::BigFloat) = !isinf(x) && !isnan(x)
 iszero(x::BigFloat) = x == Clong(0)
 isone(x::BigFloat) = x == Clong(1)
 
+# The branchy version should be faster for `BigFloat` because MPFR
+# arithmetic is expensive, and the branchless version requires more
+# arithmetic.
+Base.add12(x::T, y::T) where {T<:BigFloat} = Base.add12_branchful(x, y)
+
 @eval typemax(::Type{BigFloat}) = $(BigFloat(Inf))
 @eval typemin(::Type{BigFloat}) = $(BigFloat(-Inf))