Merge pull request #167 from t-bltg/rounding

fix subnormal `sqrt` and `cbrt`

Author: JeffreySarnoff

.DS_Store

@@ -2,3 +2,4 @@
 *.jl.*.cov
 *.jl.mem
 deps/deps.jl
+.DS_Store

.gitignore

@@ -57,14 +57,15 @@ 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"))
+    hi, lo = x
+    (isnan(hi) | iszero(hi)) && return x
+    signbit(hi) && throw(DomainError("sqrt(x) expects x >= 0"))
 
     half = T(0.5)
     dhalf = (half, zero(T))
 
-    r = inv(sqrt(HI(x)))
-    h = (HI(x) * half, LO(x) * half)
+    r = inv(sqrt(hi))
+    h = (hi * half, lo * half)
 
     r2 = mul_fpfp_dd(r, r)
     hr2 = mul_dddd_dd(h, r2)
@@ -79,6 +80,7 @@ function sqrt_dd_dd(x::Tuple{T,T}) where {T<:IEEEFloat}
     r = add_dddd_dd(r, radj)
 
     r = mul_dddd_dd(r, x)
+    isnan(HI(r)) && return DoubleFloat{T}(sqrt(hi), zero(T))  # subnormal numbers
 
     return r
 end
@@ -98,7 +100,7 @@ invcuberootsquared(A) is found iteratively using Newton's method with a final ap
 =#
 
 function cbrt_dd_dd(a::Tuple{T,T}) where {T<:IEEEFloat}
-    hi, lo = HILO(a)
+    isnan(HI(a)) && return a
     a2 = mul_dddd_dd(a,a)
     one1 = one(T)
     onethird = (0.3333333333333333, 1.850371707708594e-17)
@@ -149,5 +151,7 @@ function cbrt_dd_dd(a::Tuple{T,T}) where {T<:IEEEFloat}
     amax3 = mul_dddd_dd(amax3, onethird)
 
     ax = add_dddd_dd(ax, amax3)
+
+    isnan(HI(ax)) && return DoubleFloat{T}(cbrt(HI(a)), zero(T))  # subnormal numbers
     return ax
 end

src/.DS_Store

@@ -46,3 +46,18 @@
     @test typeof(rand(Double32)) == Double32
     @test typeof(rand(Double16)) == Double16
 end
+
+@testset "Subnormal numbers" begin
+    for (F, D) in (
+        (Float64, Double64),
+        (Float32, Double32),
+        (Float16, Double16),
+    )
+        # NOTE: using `prevfloat(floatmin(F))` doesn't expose `sqrt(D(sn))` returning a `NaN`
+        sn = D(floatmin(F) / 10)
+
+        @test issubnormal(sn)
+        @test sqrt(F(sn)) == F(sqrt(D(sn)))
+        @test cbrt(F(sn)) == F(cbrt(D(sn)))
+    end
+end