invalid sqrt and cbrt result on subnormal numbers

[t-bltg]

I'm hitting NaNs in my codes, induced by calling sqrt on subnormal numbers, with inconsistencies w.r.t Float64:

julia> x = 2.967414e-317;
julia> issubnormal(x)
true
julia> sqrt(Float64(x))
5.447397634045551e-159  # plausible
julia> sqrt(Double64(x))
NaN

Maybe related to #151, but subnormal numbers are more annoying than overflows.

[JeffreySarnoff]

Thank you. This is an issue.
With the current implementation, 5.562684646268013e-309 is the smallest value that will sqrt
values from 5.0e-324 through 5.56268464626801e-309 do not sqrt properly
A reasonable, though inexact fix is to trap values 0 < x < 5.562684646268013e-309
and return Double64(sqrt(x))
To get the correct result requires more elaborate subnormal handling.
What are your thoughts?

[JeffreySarnoff]

!! cbrt has a similar and more widespread problem !!
also

x=prevfloat(0.0,1_125_899_906_842_626); (Double64(x))^(0.1)
ERROR: DomainError with -5.562684646268013e-309:

[JeffreySarnoff]

Worth noting
DoubleFloats are not designed to work well with subnormals because the usual primitives for adding, multiplying two floats that return a high part and a low part are not designed for subnormal inputs. To make them handle subnormals would slow down everything quite a bit.

[t-bltg]

Thanks for the comments.

What are your thoughts?

It's hard to propose something without knowing the internal more.
I'll have a closer look, but I really think this should be fixed because subnormal numbers are pretty common and usually correctly absorbed within multiplication or additions (contrarily to NaNs).

[JeffreySarnoff]

Would you consider sqrt(Double64(subnormal)) == Double64(sqrt(subnormal)) sufficient?

[t-bltg]

I'm not sure I follow (do you imply truncating ?).

The following fixes my issue for sqrt (haven't thought about cbrt):

diff --git a/src/math/ops/op_dd_dd.jl b/src/math/ops/op_dd_dd.jl
index ea3ea7a..1ad0818 100644
--- a/src/math/ops/op_dd_dd.jl
+++ b/src/math/ops/op_dd_dd.jl
@@ -59,6 +59,7 @@ end
 function sqrt_dd_dd(x::Tuple{T,T}) where {T<:IEEEFloat}
     (isnan(HI(x)) | iszero(HI(x))) && return x
     signbit(HI(x)) && throw(DomainError("sqrt(x) expects x >= 0"))
+    issubnormal(HI(x)) && return x

     half = T(0.5)
     dhalf = (half, zero(T))

[JeffreySarnoff]

I was suggesting something like this, so there is some square root used

+    issubnormal(HI(x)) && return DoubleFloat{T}(sqrt(HI(x)), zero(T))

[t-bltg]

That would enhance the current (consistent with Float64), and fix my issue yes.
However, I do not know what the performance impact of calling issubnormal might be.

[t-bltg]

This also seems to fix cbrt: do you want me to make a PR ?
EDIT: see #167

 function cbrt_dd_dd(a::Tuple{T,T}) where {T<:IEEEFloat}
-    hi, lo = HILO(a)
+    issubnormal(HI(a)) && return DoubleFloat{T}(cbrt(HI(a)), zero(T))
     a2 = mul_dddd_dd(a,a)
     one1 = one(T)
     onethird = (0.3333333333333333, 1.850371707708594e-17)

[JeffreySarnoff]

thanks for the PR

[nschaeff]

For sqrt, I think the alternative algorithm proposed in issue #187 handles subnormals without issue (and is faster).

That said, my opinion is that subnormals should not exist (flush to zero is sometimes implemented by hardware). If you hit the subnormal range, you are in hazardous terrain anyway and I would not mind DoubleFloats.jl giving no guarantee in subnormal range rather than being slowed down by many checks for special cases...