Overflow behaviour of Double64 does not match Float64

[KlausC]

Please re-open this issue:
The root cause seems to be in DoubleFloats:

julia> Double64(1e300)^2
NaN

julia> Float64(1e300)^2
Inf

Originally posted by @KlausC in #149 (comment)

[KlausC]

A fix could look like so: in DoubleFloats/src/math/errorfree.jl:120

"""
    two_hilo_sum(a, b)

*unchecked* requirement `|a| ≥ |b|`
Computes `s = fl(a+b)` and `e = err(a+b)`.
"""
@inline function two_hilo_sum(a::T, b::T) where {T<:FloatWithFMA}
    isfinite(b) || return a, a
    s = a + b
    e = b - (s - a)
    return s, e
end

[JeffreySarnoff]

This is a deep, heavily called function. I need to see if there is a way to resolve the issue that is less impactful.

[JeffreySarnoff]

Fortunately, this anomaly is not pervasive.

julia> exp(Double64(1e300))
Inf
julia> tanpi(Double64(1/2))
Inf

If the behavior is limited to one or just a few functions, it makes more sense to trap this case within ^ etc. I am checking other functions for the same behavior.

[JeffreySarnoff]

This issue is based in a corner case for multiplication or squaring, when the magnitudes get huge:

julia> a=Double64(sqrt(floatmax(Float64)))
1.3407807929942596e154

julia> a*a
1.7976931348623155e308

julia> a=Double64(1.000000000000001*sqrt(floatmax(Float64)))
1.340780792994261e154

julia> a*a
NaN

The next step is looking into this more deeply.

[JeffreySarnoff]

omitted

[KlausC]

I think, the issue is with all base operations (+, -, *, /), which all return NaN in cases when they should return Inf or -Inf.
``
julia> a = floatmax(Double64)
1.7976931348623157e308

julia> a + a
NaN

julia> -a - a
NaN

julia> a * a
NaN

julia> a / Double64(0.0)
Inf

julia> a / (1/a)
NaN

julia> 1/a
5.562684646268003e-309`

[JeffreySarnoff]

That is clear, and helpful.

[JeffreySarnoff]

I found the problem -- when the result is Inf the Double64 computation can result in HILO(result) == (Inf, -Inf) [or (-Inf, Inf) which prints as NaN when it should print as HI(result). I need to see if this is enough or if the result must be remade to the general form used in the package (Inf, NaN) or (-Inf, NaN) to work in carry through calculations.

[JeffreySarnoff]

THIS IS INCORRECT, see next

The results need to be remade -- or something else modified

julia> exp(Double64(Inf))
Inf

julia> exp(Double64((Inf,-Inf)) # !! incorrect initialization
NaN

julia> dump(ans)
Double64
  hi: Float64 NaN
  lo: Float64 NaN

[JeffreySarnoff]

Ignore last comment --

julia> b=Double64(Inf,-Inf)
Inf

julia> HILO(b)
(Inf, -Inf)

julia> exp(b)
Inf

[JeffreySarnoff]

still though, there is the problem you highlighted .. and unfortunately not a string display issue

a = floatmax(Float64)
b = floatmax(Double64)

julia> a = floatmax(Float64)
1.7976931348623157e308

julia> b = floatmax(Double64)
1.7976931348623157e308

julia> a+a
Inf

julia> b+b
NaN

julia> HILO(ans)
(NaN, NaN)

[KlausC]

I think we must first make the base operation correct. The exp example should be considered later.

julia> d = Double64(Inf, -Inf)
Inf

julia> d + d
NaN

julia> 

julia> d = Double64(Inf, Inf)
Inf

julia> d + d
NaN

[JeffreySarnoff]

It appears that I have to check for Inf after the initial arithmetic operation (establishing the HI part) and return that immediately in the event it is +/- Inf (or NaN .. although that will propagate). This changes the routines in op_dddd_dd.jl and probably those where one of the args is dd and the other Float64 (although I have not checked yet).

[KlausC]

I agree. Question is where to insert the checks to avoid runtime regression afa possible.
I my experience the intermediate terms are "contaminated" with NaN , when operations like Inf - Inf are induced.