Skip to content

Commit 59010a3

Browse files
committed
fix sqrt and cbrt for subnormal numbers
1 parent ef689cc commit 59010a3

File tree

4 files changed

+8
-5
lines changed

4 files changed

+8
-5
lines changed

.DS_Store

-6 KB
Binary file not shown.

.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,3 +2,4 @@
22
*.jl.*.cov
33
*.jl.mem
44
deps/deps.jl
5+
.DS_Store

src/.DS_Store

-6 KB
Binary file not shown.

src/math/ops/op_dd_dd.jl

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -57,14 +57,16 @@ end
5757

5858

5959
function sqrt_dd_dd(x::Tuple{T,T}) where {T<:IEEEFloat}
60-
(isnan(HI(x)) | iszero(HI(x))) && return x
61-
signbit(HI(x)) && throw(DomainError("sqrt(x) expects x >= 0"))
60+
hi, lo = x
61+
(isnan(hi) | iszero(hi)) && return x
62+
signbit(hi) && throw(DomainError("sqrt(x) expects x >= 0"))
63+
issubnormal(hi) && return DoubleFloat{T}(sqrt(hi), zero(T))
6264

6365
half = T(0.5)
6466
dhalf = (half, zero(T))
6567

66-
r = inv(sqrt(HI(x)))
67-
h = (HI(x) * half, LO(x) * half)
68+
r = inv(sqrt(hi))
69+
h = (hi * half, lo * half)
6870

6971
r2 = mul_fpfp_dd(r, r)
7072
hr2 = mul_dddd_dd(h, r2)
@@ -98,7 +100,7 @@ invcuberootsquared(A) is found iteratively using Newton's method with a final ap
98100
=#
99101

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

0 commit comments

Comments
 (0)