Skip to content

Commit cb31d53

Browse files
TSGutTSGut
committed
add loggamma tests & resolve name conflict
1 parent 66a6c99 commit cb31d53

4 files changed

Lines changed: 405 additions & 124 deletions

File tree

‎src/core.jl‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
# Float64 version adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
2+
# Copied from Bessels.jl (MIT license https://github.com/JuliaMath/Bessels.jl/)
3+
"""
4+
gamma(x::T) where T <: Union{Float16, Float32, Float64}
5+
6+
Returns the gamma function, ``Γ(x)``, for real `x`.
7+
8+
```math
9+
\\Gamma(x) = \\int_x^\\infty e^{-t} t^{x-1} dt \\,
10+
```
11+
12+
External links: [DLMF](https://dlmf.nist.gov/5.2.1), [Wikipedia](https://en.wikipedia.org/wiki/Gamma_function)
13+
"""
14+
function gamma(_x::Float64)
15+
x = _x
16+
if x < 0 || !isfinite(x)
17+
isinteger(x) && throw(DomainError(_x, "NaN result for non-NaN input."))
18+
if !isfinite(x)
19+
x == -Inf && throw(DomainError(_x, "NaN result for non-NaN input."))
20+
return x
21+
end
22+
s = sinpi(x)
23+
x = -x # Use this rather than the traditional x = 1-x to avoid roundoff.
24+
s *= x
25+
end
26+
if x > 11.5
27+
w = inv(x)
28+
coefs = (1.0,
29+
8.333333333333331800504e-2, 3.472222222230075327854e-3, -2.681327161876304418288e-3, -2.294719747873185405699e-4,
30+
7.840334842744753003862e-4, 6.989332260623193171870e-5, -5.950237554056330156018e-4, -2.363848809501759061727e-5,
31+
7.147391378143610789273e-4
32+
)
33+
w = evalpoly(w, coefs)
34+
# avoid overflow
35+
v = x ^ muladd(0.5, x, -0.25)
36+
res = sqrt(2pi) * v * (v / exp(x)) * w
37+
38+
if _x < 0
39+
return π / (res * s)
40+
else
41+
return res
42+
end
43+
end
44+
P = (
45+
1.000000000000000000009e0, 8.378004301573126728826e-1, 3.629515436640239168939e-1, 1.113062816019361559013e-1,
46+
2.385363243461108252554e-2, 4.092666828394035500949e-3, 4.542931960608009155600e-4, 4.212760487471622013093e-5
47+
)
48+
Q = (
49+
9.999999999999999999908e-1, 4.150160950588455434583e-1, -2.243510905670329164562e-1, -4.633887671244534213831e-2,
50+
2.773706565840072979165e-2, -7.955933682494738320586e-4, -1.237799246653152231188e-3, 2.346584059160635244282e-4,
51+
-1.397148517476170440917e-5
52+
)
53+
54+
z = 1.0
55+
while x >= 3.0
56+
x -= 1.0
57+
z *= x
58+
end
59+
while x < 2.0
60+
z /= x
61+
x += 1.0
62+
end
63+
64+
x -= 2.0
65+
p = evalpoly(x, P)
66+
q = evalpoly(x, Q)
67+
return _x < 0 ? π * q / (s * z * p) : z * p / q
68+
end
69+
70+
function gamma(_x::Float32)
71+
x = Float64(_x)
72+
if x < 0 || !isfinite(x)
73+
isinteger(x) && throw(DomainError(_x, "NaN result for non-NaN input."))
74+
if !isfinite(x)
75+
x == -Inf32 && throw(DomainError(_x, "NaN result for non-NaN input."))
76+
return _x
77+
end
78+
s = sinpi(x)
79+
x = 1 - x
80+
end
81+
if x < 5
82+
z = 1.0
83+
while x > 1
84+
x -= 1
85+
z *= x
86+
end
87+
num = evalpoly(x, (1.0, 0.41702538904450015, 0.24081703455575904, 0.04071509011391178, 0.015839573267537377))
88+
den = x*evalpoly(x, (1.0, 0.9942411061082665, -0.17434932941689474, -0.13577921102050783, 0.03028452206514555))
89+
res = z * num / den
90+
else
91+
x -= 1
92+
w = evalpoly(inv(x), (2.506628299028453, 0.20888413086840676, 0.008736513049552962, -0.007022997182153692, 0.0006787969600290756))
93+
res = @fastmath sqrt(x) * exp(muladd(log(x), x, -x)) * w
94+
end
95+
return Float32(_x < 0 ? π / (s * res) : res)
96+
end
97+
98+
function gamma(_x::Float16)
99+
x = Float32(_x)
100+
if x < 0 || !isfinite(x)
101+
isinteger(x) && throw(DomainError(_x, "NaN result for non-NaN input."))
102+
if !isfinite(x)
103+
x == -Inf && throw(DomainError(_x, "NaN result for non-NaN input."))
104+
return _x
105+
end
106+
s = sinpi(x)
107+
x = 1 - x
108+
end
109+
x > 14 && return Float16(ifelse(_x > 0, Inf32, 0f0))
110+
z = 1f0
111+
while x > 1
112+
x -= 1
113+
z *= x
114+
end
115+
num = evalpoly(x, (1.0f0, 0.4170254f0, 0.24081704f0, 0.04071509f0, 0.015839573f0))
116+
den = x*evalpoly(x, (1.0f0, 0.9942411f0, -0.17434932f0, -0.13577922f0, 0.030284522f0))
117+
return Float16(_x < 0 ? Float32(π)*den / (s*z*num) : z * num / den)
118+
end
119+
120+
function gamma(n::Integer)
121+
n < 0 && throw(DomainError(n, "`n` must not be negative."))
122+
n == 0 && return Inf
123+
n > 20 && return gamma(Float64(n))
124+
@inbounds return Float64(factorial(n-1))
125+
end
126+
127+
gamma_near_1(x) = evalpoly(x-one(x), (1.0, -0.5772156649015329, 0.9890559953279725, -0.23263776388631713))

‎src/gamma.jl‎

Lines changed: 6 additions & 124 deletions
Original file line numberDiff line numberDiff line change
@@ -1,127 +1,9 @@
1-
# Float64 version adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
2-
# Copied from Bessels.jl (MIT license https://github.com/JuliaMath/Bessels.jl/)
3-
"""
4-
gamma(x::T) where T <: Union{Float16, Float32, Float64}
1+
module Gamma
52

6-
Returns the gamma function, ``Γ(x)``, for real `x`.
3+
include("core.jl")
4+
include("precompile.jl")
5+
include("loggamma.jl")
76

8-
```math
9-
\\Gamma(x) = \\int_x^\\infty e^{-t} t^{x-1} dt \\,
10-
```
7+
export gamma, loggamma, logabsgamma, logfactorial
118

12-
External links: [DLMF](https://dlmf.nist.gov/5.2.1), [Wikipedia](https://en.wikipedia.org/wiki/Gamma_function)
13-
"""
14-
function gamma(_x::Float64)
15-
x = _x
16-
if x < 0 || !isfinite(x)
17-
isinteger(x) && throw(DomainError(_x, "NaN result for non-NaN input."))
18-
if !isfinite(x)
19-
x == -Inf && throw(DomainError(_x, "NaN result for non-NaN input."))
20-
return x
21-
end
22-
s = sinpi(x)
23-
x = -x # Use this rather than the traditional x = 1-x to avoid roundoff.
24-
s *= x
25-
end
26-
if x > 11.5
27-
w = inv(x)
28-
coefs = (1.0,
29-
8.333333333333331800504e-2, 3.472222222230075327854e-3, -2.681327161876304418288e-3, -2.294719747873185405699e-4,
30-
7.840334842744753003862e-4, 6.989332260623193171870e-5, -5.950237554056330156018e-4, -2.363848809501759061727e-5,
31-
7.147391378143610789273e-4
32-
)
33-
w = evalpoly(w, coefs)
34-
# avoid overflow
35-
v = x ^ muladd(0.5, x, -0.25)
36-
res = sqrt(2pi) * v * (v / exp(x)) * w
37-
38-
if _x < 0
39-
return π / (res * s)
40-
else
41-
return res
42-
end
43-
end
44-
P = (
45-
1.000000000000000000009e0, 8.378004301573126728826e-1, 3.629515436640239168939e-1, 1.113062816019361559013e-1,
46-
2.385363243461108252554e-2, 4.092666828394035500949e-3, 4.542931960608009155600e-4, 4.212760487471622013093e-5
47-
)
48-
Q = (
49-
9.999999999999999999908e-1, 4.150160950588455434583e-1, -2.243510905670329164562e-1, -4.633887671244534213831e-2,
50-
2.773706565840072979165e-2, -7.955933682494738320586e-4, -1.237799246653152231188e-3, 2.346584059160635244282e-4,
51-
-1.397148517476170440917e-5
52-
)
53-
54-
z = 1.0
55-
while x >= 3.0
56-
x -= 1.0
57-
z *= x
58-
end
59-
while x < 2.0
60-
z /= x
61-
x += 1.0
62-
end
63-
64-
x -= 2.0
65-
p = evalpoly(x, P)
66-
q = evalpoly(x, Q)
67-
return _x < 0 ? π * q / (s * z * p) : z * p / q
68-
end
69-
70-
function gamma(_x::Float32)
71-
x = Float64(_x)
72-
if x < 0 || !isfinite(x)
73-
isinteger(x) && throw(DomainError(_x, "NaN result for non-NaN input."))
74-
if !isfinite(x)
75-
x == -Inf32 && throw(DomainError(_x, "NaN result for non-NaN input."))
76-
return _x
77-
end
78-
s = sinpi(x)
79-
x = 1 - x
80-
end
81-
if x < 5
82-
z = 1.0
83-
while x > 1
84-
x -= 1
85-
z *= x
86-
end
87-
num = evalpoly(x, (1.0, 0.41702538904450015, 0.24081703455575904, 0.04071509011391178, 0.015839573267537377))
88-
den = x*evalpoly(x, (1.0, 0.9942411061082665, -0.17434932941689474, -0.13577921102050783, 0.03028452206514555))
89-
res = z * num / den
90-
else
91-
x -= 1
92-
w = evalpoly(inv(x), (2.506628299028453, 0.20888413086840676, 0.008736513049552962, -0.007022997182153692, 0.0006787969600290756))
93-
res = @fastmath sqrt(x) * exp(muladd(log(x), x, -x)) * w
94-
end
95-
return Float32(_x < 0 ? π / (s * res) : res)
96-
end
97-
98-
function gamma(_x::Float16)
99-
x = Float32(_x)
100-
if x < 0 || !isfinite(x)
101-
isinteger(x) && throw(DomainError(_x, "NaN result for non-NaN input."))
102-
if !isfinite(x)
103-
x == -Inf && throw(DomainError(_x, "NaN result for non-NaN input."))
104-
return _x
105-
end
106-
s = sinpi(x)
107-
x = 1 - x
108-
end
109-
x > 14 && return Float16(ifelse(_x > 0, Inf32, 0f0))
110-
z = 1f0
111-
while x > 1
112-
x -= 1
113-
z *= x
114-
end
115-
num = evalpoly(x, (1.0f0, 0.4170254f0, 0.24081704f0, 0.04071509f0, 0.015839573f0))
116-
den = x*evalpoly(x, (1.0f0, 0.9942411f0, -0.17434932f0, -0.13577922f0, 0.030284522f0))
117-
return Float16(_x < 0 ? Float32(π)*den / (s*z*num) : z * num / den)
118-
end
119-
120-
function gamma(n::Integer)
121-
n < 0 && throw(DomainError(n, "`n` must not be negative."))
122-
n == 0 && return Inf
123-
n > 20 && return gamma(Float64(n))
124-
@inbounds return Float64(factorial(n-1))
125-
end
126-
127-
gamma_near_1(x) = evalpoly(x-one(x), (1.0, -0.5772156649015329, 0.9890559953279725, -0.23263776388631713))
9+
end

‎test/runtests.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,3 +21,5 @@ end
2121
x = [0, 1, 2, 3, 8, 15, 20, 30]
2222
@test SpecialFunctions.gamma.(x) ≈ gamma.(x)
2323
@inferred gamma(1)
24+
25+
include("test/test_loggamma.jl")

0 commit comments

Comments
 (0)