Base: TwicePrecision: improve constructors

[nsajko]

EDIT: when someone gets around to reviewing this, it'd probably be best to skip all the comments, seeing as none are relevant any more after I reduced the PR, leaving the tricky concerns for future PRs

Correctness is improved for constructing TwicePrecision{T} from TwicePrecision{S}. Previously a call like
TwicePrecision{Float64}(::TwicePrecision{Float32}) would leave hi and lo unnormalized; while a call like
TwicePrecision{Float32}(::TwicePrecision{Float64}) would additionally introduce unnecessary truncation.

Accuracy is improved for constructing TwicePrecision from types like BigFloat or Rational. Previously a call like
TwicePrecision{Float64}(::BigFloat) would unconditionally leave lo as zero.

Example code that tests the improvements (all Boolean values should evaluate to true):

const TwicePrecision = Base.TwicePrecision
const canonicalize2 = Base.canonicalize2

function twiceprecision_roundtrip_is_not_lossy(
    ::Type{S},
    x::T,
) where {S<:Number, T<:Union{Number,TwicePrecision}}
    tw = TwicePrecision{S}(x)
    x == T(tw)
end

function twiceprecision_is_normalized(tw::Tw) where {Tw<:TwicePrecision}
    (hi, lo) = (tw.hi, tw.lo)
    normalized = Tw(canonicalize2(hi, lo)...)
    (abs(lo) ≤ abs(hi)) & (tw == normalized)
end

rand_twiceprecision(::Type{T}) where {T<:Number} = TwicePrecision{T}(rand(widen(T)))

rand_twiceprecision_is_ok(::Type{T}) where {T<:Number} = !iszero(rand_twiceprecision(T).lo)

setprecision(BigFloat, 70)  # The mantissa needs to be just a bit larger than for `Float64`

rand_twiceprecision_is_ok(Float32)
rand_twiceprecision_is_ok(Float32)
rand_twiceprecision_is_ok(Float32)

rand_twiceprecision_is_ok(Float64)
rand_twiceprecision_is_ok(Float64)
rand_twiceprecision_is_ok(Float64)

twiceprecision_roundtrip_is_not_lossy(Float64, rand(BigFloat))
twiceprecision_roundtrip_is_not_lossy(Float64, rand(BigFloat))
twiceprecision_roundtrip_is_not_lossy(Float64, rand(BigFloat))

twiceprecision_roundtrip_is_not_lossy(Float64, rand_twiceprecision(Float32))
twiceprecision_roundtrip_is_not_lossy(Float64, rand_twiceprecision(Float32))
twiceprecision_roundtrip_is_not_lossy(Float64, rand_twiceprecision(Float32))

twiceprecision_is_normalized(TwicePrecision{Float32}(rand_twiceprecision(Float64)))
twiceprecision_is_normalized(TwicePrecision{Float32}(rand_twiceprecision(Float64)))
twiceprecision_is_normalized(TwicePrecision{Float32}(rand_twiceprecision(Float64)))

twiceprecision_is_normalized(TwicePrecision{Float64}(rand_twiceprecision(Float32)))
twiceprecision_is_normalized(TwicePrecision{Float64}(rand_twiceprecision(Float32)))
twiceprecision_is_normalized(TwicePrecision{Float64}(rand_twiceprecision(Float32)))

Updates #49589

[nsajko]

FTR I don't run the tests locally because make testall (or even nice -n 20 make testall) locks up my laptop and it doesn't seem possible to run the ranges tests in isolation.

[KristofferC]

make test-ranges might do it

[Seelengrab]

I think the idea of this PR is good (leaving fields uninitialized is generally bad), but if I'm not mistaken, existing calls with a single argument were already zero initialized with this, right?

TwicePrecision{T}(x::T) where {T} = TwicePrecision{T}(x, zero(T))

Adding more specialized constructors for creating a TwicePrecision{Float32} from a Float64 without loss in accuracy is a good thing, though I'm not sure the dance with to_twiceprecision and twiceprecision_to_number could not be done with fewer dispatches to regular constructors.

[nsajko]

I think the idea of this PR is good (leaving fields uninitialized is generally bad), but if I'm not mistaken, existing calls with a single argument were already zero initialized with this, right?

This PR is not about uninitialized fields, I don't think there are any issues with uninitialized fields with TwicePrecision.

Adding more specialized constructors for creating a TwicePrecision{Float32} from a Float64 without loss in accuracy is a good thing

Essentially I removed bad constructors and slightly modified some existing ones and then rearranged the code. I don't think I actually added any new logic (except for dispatch).

though I'm not sure the dance with to_twiceprecision and twiceprecision_to_number could not be done with fewer dispatches to regular constructors.

As stated in the comment above them, the motivation for doing it like that is to reuse code with complex.jl.

In case it's not clear why I modified complex.jl at all: my wish was to increase safety with making it less likely that a nonsensical TwicePrecision type will be constructed or used in any way. For example, TwicePrecision{<:Rational} or TwicePrecision{<:Integer} is something that shouldn't exist so my goal was to ensure that. Though, as you pointed out in the code review, this seems to break Unitful, so I'll have to think more about how to ensure safety.

[nsajko]

About safety: the thing with multi-word floating-point, including double-word floating-point that's used here, is that the allowed element types have to be whitelisted. My goal (with this PR and later #49569) is to have a double-word implementation for <:AbstractFloat in twiceprecision.jl and an implementation for <:Complex{<:AbstractFloat} in complex.jl, with all other combinations disallowed. This is important because otherwise we could be silently getting inaccurate numbers as a result of TwicePrecision operations. The situation with ranges of Unitful types being allowed was something I didn't expect, however, so I'll have to think about this some more.

[nsajko]

I think the only solution is to make packages like Unitful not use TwicePrecision unless they somehow opt-in. This would mean that ranges of Unitful types would see decreased accuracy until the Unitful package maintainers "opt-in".

Does this seem acceptable?

[Seelengrab]

This PR is not about uninitialized fields, I don't think there are any issues with uninitialized fields with TwicePrecision.

Apologies, I misread unnormalized! 😅 You are correct.

As stated in the comment above them, the motivation for doing it like that is to reuse code with complex.jl.

In case it's not clear why I modified complex.jl at all: my wish was to increase safety with making it less likely that a nonsensical TwicePrecision type will be constructed or used in any way.

I don't follow - can't complex.jl just simply call the constructors? Complex{<:AbstractFloat} <: Number is true after all. I don't see the appeal of introducing a layer of indirection instead of just defining the constructor, I guess 🤷

For example, I think having the specialization +(TwicePrecision{<:ComplexFloat}, TwicePrecision{<:ComplexFloat}) (as well as the complex addition with TwicePrecision{<:AbstractFloat}) defined in complex.jl would make perfect sense. Dispatch will select that method, without the need to change the implementation in twiceprecision.jl at all.

For example, TwicePrecision{<:Rational} or TwicePrecision{<:Integer} is something that shouldn't exist so my goal was to ensure that.

I see - the current changes would not have achieved that goal:

julia> Int <: Number
true

julia> Rational <: Number
true

About safety: the thing with multi-word floating-point, including double-word floating-point that's used here, is that the allowed element types have to be whitelisted. My goal (with this PR and later #49569) is to have a double-word implementation for <:AbstractFloat in twiceprecision.jl and an implementation for <:Complex{<:AbstractFloat} in complex.jl, with all other combinations disallowed.

That sounds more like a case for specializing the interactions between TwicePrecision objects, rather than trying to prevent construction of them entirely. Have you looked into how promote works? It's used extensively in mixed type arithmetic in julia - if single types are important for correctness, it's usually preferred to enforce those invariants locally (especially if changing the current behavior would break package code).

I think the only solution is to make packages like Unitful not use TwicePrecision unless they somehow opt-in.

I can't comment on what Unitful should or should not do - I think the current situation arises due to the extreme genericity of our ranges. Keeping that around/working by adding specializations for types where we know we can do better (like +(::TwicePrecision{<:ComplexFloat}, ::TwicePrecision{<:ComplexFloat})) seems preferrable to me.

[nsajko]

I reduced the PR, leaving all the tricky issues for future PRs and making it more easily reviewable.