Merge branch 'master' into ellipk

Author: inkydragon

.github/workflows/ci.yml

@@ -17,18 +17,13 @@ jobs:
           - 'min'
           - 'lts'
           - '1'
+          - 'pre'
         os:
           - ubuntu-latest
           - macOS-latest
           - windows-latest
-        exclude:
-          - os: macOS-latest # Apple Silicon
-            version: 'min'
-        include:
-          - os: macOS-13 # Intel
-            version: 'min'
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
@@ -42,13 +37,13 @@ jobs:
         with:
           token: ${{ secrets.CODECOV_TOKEN }} # required
           fail_ci_if_error: true
-          file: ./lcov.info
+          files: ./lcov.info
           flags: unittests
   docs:
     name: Documentation
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: "1"

Project.toml

@@ -1,6 +1,6 @@
 name = "SpecialFunctions"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.5.1"
+version = "2.6.1"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -16,19 +16,27 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
 [compat]
+Aqua = "0.8.14"
 ChainRulesCore = "0.9.44, 0.10, 1"
 ChainRulesTestUtils = "0.6.8, 0.7, 1"
+ExplicitImports = "1.13.2"
 IrrationalConstants = "0.1, 0.2"
+JET = "0.9, 0.10"
 LogExpFunctions = "0.3.2"
 OpenLibm_jll = "0.7, 0.8"
 OpenSpecFun_jll = "0.5"
-julia = "1.5"
+Random = "<0.0.1, 1"
+Test = "<0.0.1, 1"
+julia = "1.10"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ChainRulesCore", "ChainRulesTestUtils", "Random", "Test"]
+test = ["Aqua", "ChainRulesCore", "ChainRulesTestUtils", "ExplicitImports", "JET", "Random", "Test"]

ext/SpecialFunctionsChainRulesCoreExt.jl

@@ -1,8 +1,9 @@
 module SpecialFunctionsChainRulesCoreExt
 
-using SpecialFunctions, ChainRulesCore
+using SpecialFunctions
+using ChainRulesCore: ChainRulesCore
 
-import SpecialFunctions: sqrtπ, invπ
+using SpecialFunctions: sqrtπ, invπ
 
 const BESSEL_ORDER_INFO = """
 derivatives of Bessel functions with respect to the order are not implemented currently:

src/SpecialFunctions.jl

@@ -15,8 +15,8 @@ using IrrationalConstants:
 
 import LogExpFunctions
 
-using OpenLibm_jll
-using OpenSpecFun_jll
+using OpenLibm_jll: libopenlibm
+using OpenSpecFun_jll: libopenspecfun
 
 export
     airyai,
@@ -103,9 +103,6 @@ include("gamma.jl")
 include("gamma_inc.jl")
 include("betanc.jl")
 include("beta_inc.jl")
-if !isdefined(Base, :get_extension)
-    include("../ext/SpecialFunctionsChainRulesCoreExt.jl")
-end
 include("deprecated.jl")
 
 for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, :logerfc, :logerfcx,
@@ -119,13 +116,9 @@ for f in (:beta, :lbeta)
 end
 polygamma(m::Integer, x::Missing) = missing
 
- # In future just use `fastabs` from Base.Math
- # https://github.com/JuliaLang/julia/blob/93fb785831dcfcc442f82fab8746f0244c5274ae/base/special/trig.jl#L1057
-if isdefined(Base.Math, :fastabs)
-    import Base.Math: fastabs
-else
-    fastabs(x::Number) = abs(x)
-    fastabs(x::Complex) = abs(real(x)) + abs(imag(x))
-end
+# `fastabs` is identical to `Base.Math.fastabs` which is not used here since it is not public
+# https://github.com/JuliaLang/julia/blob/93fb785831dcfcc442f82fab8746f0244c5274ae/base/special/trig.jl#L1057
+fastabs(x::Number) = abs(x)
+fastabs(x::Complex) = abs(real(x)) + abs(imag(x))
 
 end # module

src/ellip.jl

@@ -40,12 +40,10 @@ function _ellipk(m::Float64)
     (m === -Inf) && return 0.0
 
     flag_is_m_neg = false
-    if m < 0.0
+    if flag_is_m_neg
         x               = m / (m-1)         #dealing with negative args
-        flag_is_m_neg   = true
-    elseif m >= 0.0
+    else
         x               = m
-        flag_is_m_neg   = false
     end
 
     if x == 0.0
@@ -57,7 +55,7 @@ function _ellipk(m::Float64)
     elseif x > 1.0
         throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]"))
 
-    elseif 0.0 <= x < 0.1   #Table 2 from paper
+    elseif 0.0 < x < 0.1   #Table 2 from paper
         t = x-0.05
         t = @horner(t,
             1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 ,
@@ -149,7 +147,8 @@ function _ellipk(m::Float64)
             1170222242422.439893 , 8777948323668.937971 , 66101242752484.95041 ,
             499488053713388.7989 , 37859743397240299.20)
 
-    elseif x >= 0.9
+    else
+        @assert 0.9 <= x < 1
         td  = 1-x
         td1 = td-0.05
         qd  = @horner(td,
@@ -214,12 +213,10 @@ function _ellipe(m::Float64)
     (m === -Inf) && return Inf
 
     flag_is_m_neg = false
-    if m < 0.0
+    if flag_is_m_neg
         x               = m / (m-1)         #dealing with negative args
-        flag_is_m_neg   = true
-    elseif m >= 0.0
+    else
         x               = m
-        flag_is_m_neg   = false
     end
 
     if x == 0.0
@@ -230,7 +227,7 @@ function _ellipe(m::Float64)
     elseif x > 1.0
         throw(DomainError(m,"`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]"))
 
-    elseif 0.0 <= x < 0.1   #Table 2 from paper
+    elseif 0.0 < x < 0.1   #Table 2 from paper
         t = x-0.05
         t = @horner(t  ,
             +1.550973351780472328 , -0.400301020103198524 , -0.078498619442941939 ,
@@ -317,7 +314,8 @@ function _ellipe(m::Float64)
             -16120097.81581656797 , -109209938.5203089915 , -749380758.1942496220 ,
             -5198725846.725541393 , -36409256888.12139973)
 
-    elseif x >= 0.9
+    else
+        @assert 0.9 <= x < 1.0
         td  = 1-x
         td1 = td-0.05
 

src/gamma.jl

@@ -26,7 +26,7 @@ function _digamma(z::ComplexOrReal{Float64})
     # argument," Computer Phys. Commun.  vol. 4, pp. 221–226 (1972).
     x = real(z)
     if x <= 0 # reflection formula
-        ψ = -π * _cotpi(z)
+        ψ = -π / tanpi(z)
         z = 1 - z
         x = real(z)
     else
@@ -55,13 +55,6 @@ function _digamma(x::BigFloat)
     return z
 end
 
-"""
-    _cotpi(x) = cot(π * x)
-
-Accurate for integer arguments
-"""
-_cotpi(x) = @static VERSION >= v"1.10.0-DEV.525" ? inv(tanpi(x)) : cospi(x) / sinpi(x)
-
 """
     trigamma(x)
 
@@ -139,12 +132,12 @@ const cotderiv_Q = [cotderiv_q(m) for m in 1:100]
 function cotderiv(m::Integer, z)
     isinf(imag(z)) && return zero(z)
     if m <= 0
-        m == 0 && return π * _cotpi(z)
+        m == 0 && return π / tanpi(z)
         throw(DomainError(m, "`m` must be nonnegative."))
     end
     if m <= length(cotderiv_Q)
         q = cotderiv_Q[m]
-        x = _cotpi(z)
+        x = inv(tanpi(z))
         y = x*x
         s = q[1] + q[2] * y
         t = y
@@ -816,8 +809,7 @@ function logabsbeta(a::T, b::T) where T<:Real
     if a <= 0 && isinteger(a)
         if a + b <= 0 && isinteger(b)
             r = logbeta(1 - a - b, b)
-            # in julia ≥ 1.7, iseven doesn't require Int (julia#38976)
-            sgn = iseven(@static VERSION ≥ v"1.7" ? b : Int(b)) ? 1 : -1
+            sgn = iseven(b) ? 1 : -1
             return r, sgn
         else
             return -log(zero(a)), 1

src/gamma_inc.jl

@@ -1015,6 +1015,7 @@ end
 
 function __gamma_inc_inv(a::Float64, minpq::Float64, pcase::Bool)
     haseta = false
+    fp = NaN
 
     logp = pcase ? log(minpq) : log1p(-minpq)
     loggamma1pa = a <= 1.0 ? loggamma1p(a) : loggamma(a + 1.0)

src/logabsgamma/e_lgamma_r.jl

@@ -92,32 +92,37 @@ function _logabsgamma(x::Float64)
     lx = ux % UInt32
 
     #= purge off +-inf, NaN, +-0, tiny and negative arguments =#
-    signgam = 1
-    isneg = hx < Int32(0)
     ix = hx & 0x7fffffff
-    ix ≥ 0x7ff00000 && return x * x, signgam
-    ix | lx == 0x00000000 && return Inf, signgam
+    ix ≥ 0x7ff00000 && return x * x, 1
+    ix | lx == 0x00000000 && return Inf, 1
+    isneg = hx < Int32(0)
     if ix < 0x3b900000 #= |x|<2**-70, return -log(|x|) =#
         if isneg
-            signgam = -1
-            return -log(-x), signgam
+            return -log(-x), -1
         else
-            return -log(x), signgam
+            return -log(x), 1
         end
     end
+
+    # remaining cases
     if isneg
         # ix ≥ 0x43300000 && return Inf, signgam #= |x|>=2**52, must be -integer =#
         t = sinpi(x)
-        iszero(t) && return Inf, signgam #= -integer =#
-        nadj = logπ - log(abs(t * x))
-        if t < 0.0; signgam = -1; end
-        x = -x
+        iszero(t) && return Inf, 1 #= -integer =#
+        r = logπ - log(abs(t * x)) - _logabsgamma_pos(-x, ix)
+        signgam = t < 0.0 ? -1 : 1
+    else
+        r = _logabsgamma_pos(x, ix)
+        signgam = 1
     end
+    return r, signgam
+end
+function _logabsgamma_pos(x::Float64, ix::UInt32)
     if ix < 0x40000000     #= x < 2.0 =#
         i = round(x, RoundToZero)
         f = x - i
         if f == 0.0 #= purge off 1; 2 handled by x < 8.0 branch =#
-            return 0.0, signgam
+            return 0.0
         elseif i == 1.0
             r = 0.0
             c = 1.0
@@ -169,8 +174,5 @@ function _logabsgamma(x::Float64)
     else #= 2^58 ≤ x ≤ Inf =#
         r = muladd(x, log(x), -x)
     end
-    if isneg
-        r = nadj - r
-    end
-    return r, signgam
+    return r
 end

src/logabsgamma/e_lgammaf_r.jl

@@ -20,33 +20,37 @@ function _logabsgamma(x::Float32)
     hx = reinterpret(Int32, x)
 
     #= purge off +-inf, NaN, +-0, tiny and negative arguments =#
-    signgam = 1
-    isneg = hx < Int32(0)
     ix = hx & 0x7fffffff
-    ix ≥ 0x7f800000 && return x * x, signgam
-    ix == 0x00000000 && return Inf32, signgam
+    ix ≥ 0x7f800000 && return x * x, 1
+    ix == 0x00000000 && return Inf32, 1
+    isneg = hx < Int32(0)
     if ix < 0x35000000 #= |x|<2**-21, return -log(|x|) =#
         if isneg
-            signgam = -1
-            return -log(-x), signgam
+            return -log(-x), -1
         else
-            return -log(x), signgam
+            return -log(x), 1
         end
     end
+
+    # remaining cases
     if isneg
         # ix ≥ 0x4b000000 && return Inf32, signgam #= |x|>=2**23, must be -integer =#
         t = sinpi(x)
-        t == 0.0f0 && return Inf32, signgam #= -integer =#
-        nadj = logπ - log(abs(t * x))
-        if t < 0.0f0; signgam = -1; end
-        x = -x
+        t == 0.0f0 && return Inf32, 1 #= -integer =#
+        r = logπ - log(abs(t * x)) - _logabsgamma_pos(-x, ix)
+        signgam = copysign(1, t)
+    else
+        r = _logabsgamma_pos(x, ix)
+        signgam = 1
     end
-
+    return r, signgam
+end
+function _logabsgamma_pos(x::Float32, ix::UInt32)
     if ix < 0x40000000 #= x < 2.0 =#
         i = round(x, RoundToZero)
         f = x - i
         if f == 0.0f0 #= purge off 1; 2 handled by x < 8.0 branch =#
-            return 0.0f0, signgam
+            return 0.0f0
         elseif i == 1.0f0
             r = 0.0f0
             c = 1.0f0
@@ -99,8 +103,5 @@ function _logabsgamma(x::Float32)
         #= 2^58 ≤ x ≤ Inf =#
         r = muladd(x, log(x), -x)
     end
-    if isneg
-        r = nadj - r
-    end
-    return r, signgam
+    return r
 end

test/logabsgamma.jl

@@ -156,3 +156,11 @@ x = 8.000001f0
 # (i.e. prevfloat(8.000001f0) == 8.0f0)
 # We still check appropriate behavior at 8.0f0
 @test ulp(8.0f0) < 0.4006594736129046
+
+@testset "JET" begin
+    # issue #502
+    JET.@test_call logabsgamma(1.0)
+    JET.@test_opt logabsgamma(1.0)
+    JET.@test_call logabsgamma(1f0)
+    JET.@test_opt logabsgamma(1f0)
+end