Skip to content

Commit f2d5ea4

Browse files
committed
write logabsgamma in terms of lgamma, remove loggamma, and add tests
1 parent ce7a3c3 commit f2d5ea4

File tree

5 files changed

+35
-24
lines changed

5 files changed

+35
-24
lines changed

src/math/special/gamma_erf.jl

Lines changed: 2 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -2,24 +2,18 @@ erf(x::Double64) = Double64Float128(erf, x)
22
erfc(x::Double64) = Double64Float128(erfc, x)
33
gamma(x::Double64) = Double64Float128(gamma, x)
44
lgamma(x::Double64) = Double64Float128(lgamma, x)
5-
loggamma(x::Double64) = Double64Float128(lgamma, x)
65

76
erf(x::Double32) = Double32Float128(erf, x)
87
erfc(x::Double32) = Double32Float128(erfc, x)
98
gamma(x::Double32) = Double32Float128(gamma, x)
109
lgamma(x::Double32) = Double32Float128(lgamma, x)
11-
loggamma(x::Double32) = Double32Float128(lgamma, x)
1210

1311
erf(x::Double16) = Double16Float128(erf, x)
1412
erfc(x::Double16) = Double16Float128(erfc, x)
1513
gamma(x::Double16) = Double16Float128(gamma, x)
1614
lgamma(x::Double16) = Double16Float128(lgamma, x)
17-
loggamma(x::Double16) = Double16Float128(lgamma, x)
1815

1916
function logabsgamma(x::DoubleFloat{T}) where {T<:IEEEFloat}
20-
keep_precision = precision(BigFloat)
21-
setprecision(BigFloat, 128)
22-
result = logabsgamma(BigFloat(x))
23-
setprecision(BigFloat, keep_precision)
24-
return DoubleFloat{T}(result[1]), result[2]
17+
sign = x >= 0 ? 1 : 2*mod(ceil(Int64,x),2)-1
18+
return lgamma(x), sign
2519
end

src/math/special/specialfunctions.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
import .SpecialFunctions
2-
import .SpecialFunctions: erf, erfc, gamma, lgamma, loggamma, logabsgamma,
2+
import .SpecialFunctions: erf, erfc, gamma, lgamma, logabsgamma,
33
besselj0, besselj1, besselj, bessely0, bessely1, bessely,
44
ellipk
55

test/fma.jl

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,17 @@
11
@testset "fma" begin
22

3-
@test HILO(fma(Float64(pi), Float64(phi), Double64(gamma))) == (5.660419357216792, 4.1344240638440936e-16)
4-
@test HILO(fma(Float64(pi), Double64(phi), Float64(gamma))) == (5.660419357216792, 2.477303893634162e-16)
5-
@test HILO(fma(Double64(pi), Float64(phi), Float64(gamma))) == (5.660419357216793, -2.7164108363986695e-16)
6-
@test HILO(fma(Double64(pi), Double64(phi), Float64(gamma))) == (5.660419357216793, -4.422960158132908e-16)
7-
@test HILO(fma(Double64(pi), Float64(phi), Double64(gamma))) == (5.660419357216793, -2.765839987922976e-16)
8-
@test HILO(fma(Double64(pi), Double64(phi), Double64(gamma))) == (5.660419357216792, 4.4093948873440384e-16)
3+
@test HILO(fma(Float64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216792, 4.1344240638440936e-16)
4+
@test HILO(fma(Float64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216792, 2.477303893634162e-16)
5+
@test HILO(fma(Double64(pi), Float64(phi), Float64(eulergamma))) == (5.660419357216793, -2.7164108363986695e-16)
6+
@test HILO(fma(Double64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216793, -4.422960158132908e-16)
7+
@test HILO(fma(Double64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216793, -2.765839987922976e-16)
8+
@test HILO(fma(Double64(pi), Double64(phi), Double64(eulergamma))) == (5.660419357216792, 4.4093948873440384e-16)
99

10-
@test HILO(muladd(Float64(pi), Float64(phi), Double64(gamma))) == (5.660419357216792, 4.1344240638440936e-16)
11-
@test HILO(muladd(Float64(pi), Double64(phi), Float64(gamma))) == (5.660419357216792, 2.477303893634162e-16)
12-
@test HILO(muladd(Double64(pi), Float64(phi), Float64(gamma))) == (5.660419357216793, -2.7164108363986695e-16)
13-
@test HILO(muladd(Double64(pi), Double64(phi), Float64(gamma))) == (5.660419357216793, -4.422960158132908e-16)
14-
@test HILO(muladd(Double64(pi), Float64(phi), Double64(gamma))) == (5.660419357216793, -2.765839987922976e-16)
15-
@test HILO(muladd(Double64(pi), Double64(phi), Double64(gamma))) == (5.660419357216792, 4.4093948873440384e-16)
10+
@test HILO(muladd(Float64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216792, 4.1344240638440936e-16)
11+
@test HILO(muladd(Float64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216792, 2.477303893634162e-16)
12+
@test HILO(muladd(Double64(pi), Float64(phi), Float64(eulergamma))) == (5.660419357216793, -2.7164108363986695e-16)
13+
@test HILO(muladd(Double64(pi), Double64(phi), Float64(eulergamma))) == (5.660419357216793, -4.422960158132908e-16)
14+
@test HILO(muladd(Double64(pi), Float64(phi), Double64(eulergamma))) == (5.660419357216793, -2.765839987922976e-16)
15+
@test HILO(muladd(Double64(pi), Double64(phi), Double64(eulergamma))) == (5.660419357216792, 4.4093948873440384e-16)
1616

1717
end

test/runtests.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,12 +9,12 @@ using Printf
99
using Random
1010

1111
const phi = Base.MathConstants.golden
12-
const gamma = Base.MathConstants.eulergamma
12+
const eulergamma = Base.MathConstants.eulergamma
1313

1414
const phi_df64 = Double64(phi)
1515
const phi_df32 = Double32(phi)
16-
const gamma_df64 = Double64(gamma)
17-
const gamma_df32 = Double32(gamma)
16+
const gamma_df64 = Double64(eulergamma)
17+
const gamma_df32 = Double32(eulergamma)
1818

1919
const phi_df64hi = HI(phi_df64); const phi_df64lo = LO(phi_df64)
2020
const phi_df32hi = HI(phi_df32); const phi_df32lo = LO(phi_df32)

test/special_functions.jl

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1 +1,18 @@
1+
using SpecialFunctions
12

3+
@testset "special functions" begin
4+
# The intention here is not to check the accuracy of the libraries,
5+
# rather the sanity of library calls:
6+
piq = Double64(pi)
7+
halfq = Double64(0.5)
8+
@test gamma(halfq) sqrt(piq)
9+
for func in (erf, erfc, besselj0, besselj1, bessely0, bessely1, loggamma)
10+
@test func(halfq) func(0.5)
11+
end
12+
for func in (bessely, besselj)
13+
@test func(3,halfq) func(3,0.5)
14+
end
15+
@test gamma(Double64(2),3) gamma(2,3)
16+
@test all(logabsgamma(Double64(-0.5)) .≈ logabsgamma(-0.5))
17+
@test all(logabsgamma(Double64(-1.5)) .≈ logabsgamma(-1.5))
18+
end

0 commit comments

Comments
 (0)