Merge pull request #195 from araujoms/incgamma

add logabsgamma and unrelated fixes

Author: JeffreySarnoff

Project.toml

@@ -18,7 +18,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 Polynomials = "1, 2, 3, 4, 5, 6"
-Quadmath = "0.4, 0.5, 0.6, 0.7, 0.8"
+Quadmath = "0.5.10, 0.6, 0.7, 0.8"
 Requires = "1"
 SpecialFunctions = "1.0, 1.1, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5"
 GenericLinearAlgebra = "0.2.5, 0.3, 0.4, 0.5, 0.6"

src/math/ops/op_dd_dd.jl

@@ -57,7 +57,7 @@ end
 
 function sqrt_dd_dd(x::Tuple{T,T}) where {T<:IEEEFloat}
     iszero(HI(x)) && return x
-    signbit(HI(x)) && throw(DomainError("sqrt(x) expects x >= 0")) 
+    signbit(HI(x)) && throw(DomainError("sqrt(x) expects x >= 0"))
 
     ahi, alo = HILO(x)
     s = sqrt(ahi)

src/math/special/gamma_erf.jl

@@ -1,17 +1,16 @@
 erf(x::Double64) = Double64Float128(erf, x)
 erfc(x::Double64) = Double64Float128(erfc, x)
 gamma(x::Double64) = Double64Float128(gamma, x)
-lgamma(x::Double64) = Double64Float128(lgamma, x)
-loggamma(x::Double64) = Double64Float128(lgamma, x)
 
 erf(x::Double32) = Double32Float128(erf, x)
 erfc(x::Double32) = Double32Float128(erfc, x)
 gamma(x::Double32) = Double32Float128(gamma, x)
-lgamma(x::Double32) = Double32Float128(lgamma, x)
-loggamma(x::Double32) = Double32Float128(lgamma, x)
 
 erf(x::Double16) = Double16Float128(erf, x)
 erfc(x::Double16) = Double16Float128(erfc, x)
 gamma(x::Double16) = Double16Float128(gamma, x)
-lgamma(x::Double16) = Double16Float128(lgamma, x)
-loggamma(x::Double16) = Double16Float128(lgamma, x)
+
+function logabsgamma(x::DoubleFloat{T}) where {T<:IEEEFloat}
+    result = logabsgamma(Float128(x))
+    return DoubleFloat{T}(result[1]), result[2]
+end

src/math/special/specialfunctions.jl

@@ -1,6 +1,7 @@
 import .SpecialFunctions
-import .SpecialFunctions: erf, erfc, gamma, lgamma,
-    besselj0, besselj1, besselj, bessely0, bessely1, bessely
+import .SpecialFunctions: erf, erfc, gamma, logabsgamma,
+    besselj0, besselj1, besselj, bessely0, bessely1, bessely,
+    ellipk
 
 include("agm.jl")
 include("bessel.jl")

test/fma.jl

@@ -1,17 +1,17 @@
 @testset "fma" begin
 
-  @test HILO(fma(Float64(pi), Float64(phi), Double64(gamma))) == (5.660419357216792, 4.1344240638440936e-16)
-  @test HILO(fma(Float64(pi), Double64(phi), Float64(gamma))) == (5.660419357216792, 2.477303893634162e-16)
-  @test HILO(fma(Double64(pi), Float64(phi), Float64(gamma))) == (5.660419357216793, -2.7164108363986695e-16)
-  @test HILO(fma(Double64(pi), Double64(phi), Float64(gamma))) == (5.660419357216793, -4.422960158132908e-16)
-  @test HILO(fma(Double64(pi), Float64(phi), Double64(gamma))) == (5.660419357216793, -2.765839987922976e-16)
-  @test HILO(fma(Double64(pi), Double64(phi), Double64(gamma))) == (5.660419357216792, 4.4093948873440384e-16)
+  @test HILO(fma(Float64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216792, 4.1344240638440936e-16)
+  @test HILO(fma(Float64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216792, 2.477303893634162e-16)
+  @test HILO(fma(Double64(pi), Float64(phi), Float64(eulergamma))) == (5.660419357216793, -2.7164108363986695e-16)
+  @test HILO(fma(Double64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216793, -4.422960158132908e-16)
+  @test HILO(fma(Double64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216793, -2.765839987922976e-16)
+  @test HILO(fma(Double64(pi), Double64(phi), Double64(eulergamma))) == (5.660419357216792, 4.4093948873440384e-16)
 
-  @test HILO(muladd(Float64(pi), Float64(phi), Double64(gamma))) == (5.660419357216792, 4.1344240638440936e-16)
-  @test HILO(muladd(Float64(pi), Double64(phi), Float64(gamma))) == (5.660419357216792, 2.477303893634162e-16)
-  @test HILO(muladd(Double64(pi), Float64(phi), Float64(gamma))) == (5.660419357216793, -2.7164108363986695e-16)
-  @test HILO(muladd(Double64(pi), Double64(phi), Float64(gamma))) == (5.660419357216793, -4.422960158132908e-16)
-  @test HILO(muladd(Double64(pi), Float64(phi), Double64(gamma))) == (5.660419357216793, -2.765839987922976e-16)
-  @test HILO(muladd(Double64(pi), Double64(phi), Double64(gamma))) == (5.660419357216792, 4.4093948873440384e-16)
+  @test HILO(muladd(Float64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216792, 4.1344240638440936e-16)
+  @test HILO(muladd(Float64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216792, 2.477303893634162e-16)
+  @test HILO(muladd(Double64(pi), Float64(phi), Float64(eulergamma))) == (5.660419357216793, -2.7164108363986695e-16)
+  @test HILO(muladd(Double64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216793, -4.422960158132908e-16)
+  @test HILO(muladd(Double64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216793, -2.765839987922976e-16)
+  @test HILO(muladd(Double64(pi), Double64(phi), Double64(eulergamma))) == (5.660419357216792, 4.4093948873440384e-16)
 
 end

test/runtests.jl

@@ -9,12 +9,12 @@ using Printf
 using Random
 
 const phi = Base.MathConstants.golden
-const gamma = Base.MathConstants.eulergamma
+const eulergamma = Base.MathConstants.eulergamma
 
 const phi_df64   = Double64(phi)
 const phi_df32   = Double32(phi)
-const gamma_df64 = Double64(gamma)
-const gamma_df32 = Double32(gamma)
+const gamma_df64 = Double64(eulergamma)
+const gamma_df32 = Double32(eulergamma)
 
 const phi_df64hi = HI(phi_df64); const phi_df64lo = LO(phi_df64)
 const phi_df32hi = HI(phi_df32); const phi_df32lo = LO(phi_df32)

test/special_functions.jl

@@ -1 +1,18 @@
+using SpecialFunctions
 
+@testset "special functions" begin
+    # The intention here is not to check the accuracy of the libraries,
+    # rather the sanity of library calls:
+    piq = Double64(pi)
+    halfq = Double64(0.5)
+    @test gamma(halfq) ≈ sqrt(piq)
+    for func in (erf, erfc, besselj0, besselj1, bessely0, bessely1, loggamma)
+        @test func(halfq) ≈ func(0.5)
+    end
+    for func in (bessely, besselj)
+        @test func(3,halfq) ≈ func(3,0.5)
+    end
+    @test gamma(Double64(2),3) ≈ gamma(2,3)
+    @test all(logabsgamma(Double64(-0.5)) .≈ logabsgamma(-0.5))
+    @test all(logabsgamma(Double64(-1.5)) .≈ logabsgamma(-1.5))
+end