base/twiceprecision.jl needs work

[nsajko]

Algorithms for arithmetic (+, -, *, /)

The algorithms for real (as opposed to complex) double-word addition, multiplication and division were seemingly (there are no relevant explanatory comments or references in the code) all created ad-hoc and need to be updated to the standard algos, something I'll try to do in #49569.

Regarding the algorithms for complex numbers, they seem to be missing entirely, causing Base.TwicePrecision{<:Complex} to fall back to the algorithms for real numbers. This is OK for addition because Complex is a Cartesian (orthogonal) representation of complex numbers, but multiplication and division could be quite broken. This could be fixed by adding methods for Base.TwicePrecision{<:Complex} to base/complex.jl. A paper exists for double-word complex multiplication, and I guess I could work something out for division.

Design of the Base.TwicePrecision type

The Base.TwicePrecision type has an awkward and unsafe design. It accepts elements of any type, not just floating-point-based types. I propose this design instead as a replacement:

In base/twiceprecision.jl:

abstract type AbstractMultiWordFloatingPoint{n,T<:Number} end

# some methods for `AbstractMultiWordFloatingPoint` ...

struct MultiWordFloatingPoint{n,T<:AbstractFloat} <:
    AbstractMultiWordFloatingPoint{n,T}

    words::NTuple{n,T}
end

# some methods for `MultiWordFloatingPoint` ...

const DoubleWordFloatingPoint =
    MultiWordFloatingPoint{2,T} where {T<:AbstractFloat}

# some methods for `DoubleWordFloatingPoint`, including arithmetic ...

In base/complex.jl:

struct MultiWordFloatingPointComplex{n,T<:Complex{<:AbstractFloat}} <:
    AbstractMultiWordFloatingPoint{n,T}

    words::NTuple{n,T}
end

# some methods for `MultiWordFloatingPointComplex` ...

const DoubleWordFloatingPointComplex =
    MultiWordFloatingPointComplex{2,T} where {T<:AbstractFloat}

# some methods for `DoubleWordFloatingPointComplex`, including arithmetic ...

Note that MultiWordFloatingPointComplex only makes sense because Complex is an orthogonal coordinate representation of the complex numbers. A polar representation, for example, would be nonsensical in this context because a multi-word representation represents an unevaluated sum.

The n type parameter would be restricted to two in practice for the time being, however I think that having it as a parameter right away gives a more elegant design.

Is supporting complex numbers even necessary?

The point of Base.TwicePrecision is supporting ranges, and while ranges with complex elements do exist, I think the step is always real? This seems like it would mean that double-word algorithms for complex numbers are an unnecessary evil.

EDIT: the answer seems to be yes, supporting complex numbers is necessary:

julia> (1.0:5.0) .* im
0.0 + 1.0im:0.0 + 1.0im:0.0 + 5.0im

julia> (1.0:5.0) .* im .* im
-1.0 + 0.0im:-1.0 + 0.0im:-5.0 + 0.0im

splitprec

This is more of a question than an issue, but I find this dubious so I thought it'd be good to raise it here with the other stuff.

Why does TwicePrecision{<:AbstractFloat}(::Integer) use that splitprec function:

https://github.com/JuliaLang/julia/blob/master/base/twiceprecision.jl#L213-L214

The fallback conversion method TwicePrecision{<:Any}(::Any) is sufficient for converting integers to double-word floating-point:

https://github.com/JuliaLang/julia/blob/master/base/twiceprecision.jl#L203-L207

So it seems that the intention is to give up accuracy on purpose while converting integers to TwicePrecision? It's likely that I'm missing something, but I find this quite strange and I can't imagine why would one desire such behavior. Also, it seems that something similar is happening when the truncbits function is used.

[timholy]

As the author of twiceprecision, I would love it if someone who knew this stuff fixed it. The only reason I worked on this at all was because I wanted to make ranges generic (#18777) and I couldn't do that while breaking a bunch of floating-point precision tests. So the current implementation is me thrashing my way towards not breaking any tests, but a more principled effort would be appreciated. I have no useful input to provide about the design, I'm sure you'll come up with something better.

[Seelengrab]

So it seems that the intention is to give up accuracy on purpose while converting integers to TwicePrecision?

That'd be contradicting its docstring:

In particular, while convert(Float64, i) can be lossy since Float64 has only 53 bits of precision, splitprec(Float64, i) is exact for any Int64/UInt64.

Also, regarding more words:

The Base.TwicePrecision type has an awkward and unsafe design. It accepts elements of any type, not just floating-point-based types. I propose this design instead as a replacement:

abstract type AbstractMultiWordFloatingPoint{n,T<:Number} end
...
struct MultiWordFloatingPointComplex{n,T<:Complex{<:AbstractFloat}} <:

Does Base have the ambition to implement that, or should this be in a package? TwicePrecision is only an internal detail for (mostly accurate) ranges after all, not necessarily a starting point for hypercomplex numbers.

[nsajko]

Does Base have the ambition to implement that, or should this be in a package? TwicePrecision is only an internal detail for (mostly accurate) ranges after all, not necessarily a starting point for hypercomplex numbers.

Not sure whether you're aware that TwicePrecision{<:Complex} is already a thing, although the current implementation is obviously lacking. It's used in the code snippets I added above in an edit, or in (1.0:5.0) .+ im for example.

[Seelengrab]

I'm specifically referring to the n-dimensionality of AbstractMultiWordFloatingPoint{n,T<:Number}. Complex currently is decidedly two dimensional after all - making this more generic over N dimensions to me seems more suited for a package than Base.

The non-Number element type is, I think, intentional, to support usecases like Unitful.jl, which aren't <: Number.

[timholy]

Adding to @Seelengrab's points, note that the limited functionality is largely intentional. Were this functionality to become public, one then must ask how you create ranges of TwicePrecision objects? You might need TwicePrecision{TwicePrecision{T}} which seems potentially difficult. If TwicePrecision remains hidden away then it can continue (soon, be better at!) generating accurate ranges without having to worry about recursion.

[Seelengrab]

Yes - TwicePrecision is a tool for getting more accurate ranges of floating point numbers, not a tool for actually implementing generic DoubleWord arithmetic (though of course that's what it emulates under the hood..). The existence of the PhysQuantities tests involving them makes me a bit nervous in terms of what we can actually change about TwicePrecision, but my gut says that improvements purely related to TwicePrecision{<:IEEEFloat} and TwicePrecision{<:Complex{<:IEEEFloat}} and interactions among these two ought to be in reach (and is what I think @nsajko is trying to achieve).

[nsajko]

Just disregard MultiWordFloatingPoint then, it's really not important.