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.

[simonbyrne]

Thanks for the enthusiasm!

As noted above, TwicePrecision code is for floating point range computation:

julia/base/twiceprecision.jl

Lines 168 to 173 in c90aa77

	!!! warning
	`TwicePrecision` is an internal type used to increase the
	precision of floating-point ranges, and not intended for external use.
	If you encounter them in real code, the most likely explanation is
	that you are directly accessing the fields of a range. Use
	the function interface instead, `step(r)` rather than `r.step`

Many of the original design decisions were based on that (in particular, the original design didn't rely on FMAs, as the hardware operation wasn't available on all architectures, but that has since become less of a concern). It isn't really intended to be used as a generic number type (and so doesn't implement all the necessary operations one would want). For a fully-featured type, there is the DoubleFloats.jl package.

The primary downside to making it a full-featured type in Base is that it is extremely difficult to upgrade user code: if you make improvements or fixes to DoubleFloats.jl, then any user can Pkg.update() to get these. Changes in Base can take more than 6 months to wind their way into the release version of Julia.

[nsajko]

I never said I want TwicePrecision made public or anything. I just want to fix it as much as possible.

[KristofferC]

I just want to fix it as much as possible.

I think a good way of going about that is to point to areas where the float range code have bugs and how improving TwicePrecision fixes those bugs

[nsajko]

When I said

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.

I was referring to this requirement stated in splitprecs documentation:

julia/base/twiceprecision.jl

Line 16 in c90aa77

- `hi` can be exactly multiplied by the `hi` component of another call to `splitprec`.

The requirement seems to be satisfied by ensuring that the least significant bits of the most significant word (hi) are zero, basically reserving low bits of the mantissa for future operations, so I guess that it's also related to these tests of splitprec:

julia/test/ranges.jl

Line 136 in 9201414

@test endswith(bitstring(hi), repeat('0', Base.Math.significand_bits(T) ÷ 2))

The other requirements satisfied by splitprec seem natural for double-word numbers, while this, last, requirement actually causes loss of precision (because of the mandatory zero bits), making it really stand out.

Can someone who remembers please explain the rationale behind this requirement?

Interestingly, the only place where splitprec is used (apart from in tests), is for constructing TwicePrecision objects from integers:

julia/base/twiceprecision.jl

Line 212 in c90aa77

TwicePrecision{T}(canonicalize2(splitprec(T, i)...)...)

It would, however, seem much more natural to just use the more general constructor immediately above, the effect would be the same except that splitprec reserves precision for later to satisfy the requirement in question.

[simonbyrne]

I'd have to dig through the diffs, but the original motivation was for what is known as Veltkamp splitting, as a way to get higher-precision multiplication. Now that we have hardware fma ops on most platforms, I think most of the code has switched to using that (but I haven't looked for sure).

[JeffreySarnoff]

I'd have to dig through the diffs, but the original motivation was for what is known as Veltkamp splitting, as a way to get higher-precision multiplication. Now that we have hardware fma ops on most platforms, I think most of the code has switched to using that (but I haven't looked for sure).

Most code has switched to using fma, either directly or through llvm allowed muladd->fma or via MuladdMacro.jl afaik