Merge branch 'main' into dependabot/github_actions/actions/checkout-4

ViralBShah · web-flow · commit c4eef3741bd1 · 2024-04-03T00:00:11.000-04:00
diff --git a/.DS_Store b/.DS_Store
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -22,11 +22,11 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:
@@ -39,6 +39,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
@@ -15,7 +15,7 @@ jobs:
         run: which julia
         continue-on-error: true
       - name: Install Julia, but only if it is not already available in the PATH
-        uses: julia-actions/setup-julia@v1
+        uses: julia-actions/setup-julia@v2
         with:
           version: '1'
           arch: ${{ runner.arch }}
diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml
@@ -16,7 +16,7 @@ jobs:
     if: github.base_ref == github.event.repository.default_branch
     runs-on: ubuntu-latest
     steps:
-    - uses: julia-actions/setup-julia@v1
+    - uses: julia-actions/setup-julia@v2
       with:
         version: '1'
     - uses: actions/checkout@v4
diff --git a/.gitignore b/.gitignore
@@ -2,3 +2,4 @@
 *.jl.*.cov
 *.jl.mem
 deps/deps.jl
+.DS_Store
diff --git a/Project.toml b/Project.toml
@@ -4,7 +4,7 @@ keywords = ["math", "floating-point", "doubledouble", "extended-precision", "acc
 license = "MIT"
 desc = "extended precision floating point types using value pairs"
 repo = "https://github.com/JuliaMath/DoubleFloats.jl.git"
-version = "1.3.3"
+version = "1.3.8"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -3,7 +3,7 @@
 ### Math with 85+ accurate bits.
 #### Extended precision float and complex types
 
-- N.B. `Double64` is the most performant type <sup>[β](#involvement)</sup>
+- N.B. `Double64` is the most performant type [^β]
 
 ----
 
@@ -136,12 +136,10 @@ Usage questions can be posted on the [Julia Discourse forum](https://discourse.j
 
 ## Contributions
 
-Contributions are very welcome, as are feature requests and suggestions. Please open an [issue](https://github.com/JuliaMath/DoubleFloats.jl/issues) if you encounter any problems. The [contributing page](https://juliamath.github.io/DoubleFloats.jl/latest/man/contributing/) has a few guidelines that should be followed when opening pull requests.
+Contributions are very welcome, as are feature requests and suggestions. Please open an [issue](https://github.com/JuliaMath/DoubleFloats.jl/issues) if you encounter any problems. The [contributing page](https://juliamath.github.io/DoubleFloats.jl/latest/#Contributions) has a few guidelines that should be followed when opening pull requests.
 
 ----
 
-<a name="involvement">β</a>: If you want to get involved with moving `Double32` performance forward, great. I would provide guidance. Otherwise, for most purposes you are better off using `Float64` than `Double32` (`Float64` has more significant bits, wider exponent range, and is much faster).
+[^β]: If you want to get involved with moving `Double32` performance forward, great. I would provide guidance. Otherwise, for most purposes you are better off using `Float64` than `Double32` (`Float64` has more significant bits, wider exponent range, and is much faster).
 
 ----
-
-
diff --git a/docs/src/linearalgebra.md b/docs/src/linearalgebra.md
@@ -14,7 +14,7 @@ using DoubleFloats, LinearAlgebra
 n = 25
 vector = rand(Double64, n)
 matrix = reshape(rand(Double64,n*n),n,n)
-
+```
 
 ## Matrix Predicates
 
@@ -54,5 +54,3 @@ matrix = reshape(rand(Double64,n*n),n,n)
 - `asinh`, `acosh`, `atanh`, `acsch`, `asech`, `acoth`
 
 - matrixfunction(function, matrix)
-
-
diff --git a/src/.DS_Store b/src/.DS_Store
diff --git a/src/Double64.jl b/src/Double64.jl
diff --git a/src/math/elementary/explog.jl b/src/math/elementary/explog.jl
@@ -56,7 +56,6 @@ end
     return z
 end
 
-
 function mul_by_half(r::DoubleFloat{T}) where {T<:IEEEFloat}
     frhi, xphi = frexp(HI(r))
     frlo, xplo = frexp(LO(r))
@@ -120,20 +119,20 @@ function Base.:(^)(r::DoubleFloat{T}, n::Int) where {T<:IEEEFloat}
     return s
 end
 
-function Base.:(^)(r::DoubleFloat{T}, n::DoubleFloat{T}) where {T<:IEEEFloat}
-   if isinteger(n)
-      return r^Int(n)
-   else
-      return exp(n * log(r))
-   end
+function Base.:(^)(r::DoubleFloat{T}, n::DoubleFloat{T}) where {T <: IEEEFloat}
+    if isinteger(n)
+        return r^Int64(Float64(n)) # convert n to Float64 first
+    else
+       return exp(n * log(r))
+    end
 end
 
 function Base.:(^)(r::Int, n::DoubleFloat{T}) where {T<:IEEEFloat}
-   if isinteger(n)
-      return r^Int(n)
-   else
-      return exp(n * log(r))
-   end
+    if isinteger(n)
+        return r^Int64(Float64(n)) # convert n to Float64 first
+    else
+        return exp(n * log(r))
+    end
 end
 
 function calc_exp(a::DoubleFloat{T}) where {T<:IEEEFloat}
diff --git a/src/math/ops/op_dd_dd.jl b/src/math/ops/op_dd_dd.jl
@@ -55,32 +55,16 @@ end
     return zhi, zlo
 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"))
-
-    half = T(0.5)
-    dhalf = (half, zero(T))
-
-    r = inv(sqrt(HI(x)))
-    h = (HI(x) * half, LO(x) * half)
-
-    r2 = mul_fpfp_dd(r, r)
-    hr2 = mul_dddd_dd(h, r2)
-    radj = sub_fpdd_dd(half, hr2)
-    radj = mul_ddfp_dd(radj, r)
-    r = add_fpdd_dd(r, radj)
-
-    r2 = mul_dddd_dd(r, r)
-    hr2 = mul_dddd_dd(h, r2)
-    radj = sub_fpdd_dd(half, hr2)
-    radj = mul_dddd_dd(radj, r)
-    r = add_dddd_dd(r, radj)
-
-    r = mul_dddd_dd(r, x)
-
-    return r
+    iszero(HI(x)) && return x
+    signbit(HI(x)) && throw(DomainError("sqrt(x) expects x >= 0")) 
+
+    ahi, alo = HILO(x)
+    s = sqrt(ahi)
+    d = fma(-s, s, ahi)    # ahi=s*s+d,  same order of magnitude than alo, so we can add alo safely below:
+    d += alo               #  ahi+alo = s*s+d = s*s*(1+d/(s*s))   ==> sqrt(ahi+alo) = s*sqrt(1+d/s2) approx= s*(1+d/(2s2)) = s + d/(2*s)
+    d = d/(s+s)
+    return s,d
 end
 
 #=
@@ -98,7 +82,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 +133,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
diff --git a/src/type/convert.jl b/src/type/convert.jl
@@ -11,7 +11,11 @@ convert(::Type{Float32}, x::Double32) = HI(x)
 convert(::Type{Float16}, x::Double16) = HI(x)
 
 convert(::Type{DoubleFloat{T}}, x::I) where {T<:IEEEFloat, I<:Integer} = DoubleFloat{T}(T(x))
-convert(::Type{I}, x::DoubleFloat{T}) where {T<:IEEEFloat, I<:Integer} = I(round(BigInt, BigFloat(HI(x))+ BigFloat(LO(x))))
+
+function convert(::Type{I}, x::DoubleFloat{T}) where {T<:IEEEFloat, I<:Integer}
+    bi = trunc(BigInt, BigFloat(HI(x)) + BigFloat(LO(x)))
+    convert(Int64, bi)
+end
 
 function convert(::Type{BigFloat}, x::DoubleFloat{T}) where {T<:IEEEFloat}
     res = Base.BigFloat(HI(x)) + Base.BigFloat(LO(x))
diff --git a/test/concrete_accuracy.jl b/test/concrete_accuracy.jl
@@ -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
