Skip to content

Commit c0442c2

Browse files
authored
Use mpfr_lngamma for loggamma(::BigFloat) (#351)
* Use `mpfr_lngamma` * Bump version * Rename variable * Throw error if `gamma(x) < 0` * Add more tests * Fix whitespace * Improve tests
1 parent 0c41812 commit c0442c2

File tree

3 files changed

+70
-28
lines changed

3 files changed

+70
-28
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "SpecialFunctions"
22
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
3-
version = "1.7.0"
3+
version = "1.8.0"
44

55
[deps]
66
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/gamma.jl

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license
22

3-
using Base.MPFR: ROUNDING_MODE
3+
using Base.MPFR: MPFRRoundingMode, ROUNDING_MODE
44

55
export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, logabsbinomial
66

@@ -643,12 +643,20 @@ See also [`logabsgamma`](@ref) for real `x`.
643643
"""
644644
loggamma(x::Number) = _loggamma(float(x))
645645

646-
function _loggamma(x::Union{Real,BigFloat})
646+
function _loggamma(x::Real)
647647
(y, s) = logabsgamma(x)
648648
s < 0 && throw(DomainError(x, "`gamma(x)` must be non-negative"))
649649
return y
650650
end
651651

652+
function _loggamma(x::BigFloat)
653+
isnan(x) && return x
654+
y = BigFloat()
655+
ccall((:mpfr_lngamma, :libmpfr), Cint, (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), y, x, ROUNDING_MODE[])
656+
isnan(y) && throw(DomainError(x, "`gamma(x)` must be non-negative"))
657+
return y
658+
end
659+
652660
# Compute the logΓ(z) function using a combination of the asymptotic series,
653661
# the Taylor series around z=1 and z=2, the reflection formula, and the shift formula.
654662
# Many details of these techniques are discussed in D. E. G. Hare,

test/gamma.jl

Lines changed: 59 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -80,31 +80,33 @@
8080
@test_throws MethodError logfactorial(1.0)
8181
end
8282

83-
# loggamma & logabsgamma test cases (from Wolfram Alpha)
84-
@test loggamma(-300im) -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im
85-
@test loggamma(3.099) loggamma(3.099+0im) 0.786413746900558058720665860178923603134125854451168869796
86-
@test loggamma(1.15) loggamma(1.15+0im) -0.06930620867104688224241731415650307100375642207340564554
87-
@test logabsgamma(0.89)[1] loggamma(0.89+0im) 0.074022173958081423702265889979810658434235008344573396963
88-
@test loggamma(0.91) loggamma(0.91+0im) 0.058922567623832379298241751183907077883592982094770449167
89-
@test loggamma(0.01) loggamma(0.01+0im) 4.599479878042021722513945411008748087261001413385289652419
90-
@test loggamma(-3.4-0.1im) -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im
91-
@test loggamma(-13.4-0.1im) -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im
92-
@test loggamma(-13.4+0.0im) conj(loggamma(-13.4-0.0im)) -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im
93-
@test loggamma(-13.4+8im) -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im
94-
@test logabsgamma(1+exp2(-20))[1] loggamma(1+exp2(-20)+0im) -5.504750066148866790922434423491111098144565651836914e-7
95-
@test loggamma(1+exp2(-20)+exp2(-19)*im) -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im
96-
@test loggamma(-300+2im) -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im
97-
@test loggamma(300+2im) 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im
98-
@test loggamma(1-6im) -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im
99-
@test loggamma(1-8im) -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im
100-
@test loggamma(1+6.5im) conj(loggamma(1-6.5im)) -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im
101-
@test loggamma(1+1im) conj(loggamma(1-1im)) -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im
102-
@test loggamma(-pi*1e7 + 6im) -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im
103-
@test loggamma(-pi*1e7 + 8im) -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im
104-
@test loggamma(-pi*1e14 + 6im) -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im
105-
@test loggamma(-pi*1e14 + 8im) -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im
106-
@test loggamma(2.05 + 0.03im) conj(loggamma(2.05 - 0.03im)) 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im
107-
@test loggamma(2+exp2(-20)+exp2(-19)*im) 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im
83+
# values taken from Wolfram Alpha
84+
@testset "loggamma & logabsgamma test cases" begin
85+
@test loggamma(-300im) -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im
86+
@test loggamma(3.099) loggamma(3.099+0im) 0.786413746900558058720665860178923603134125854451168869796
87+
@test loggamma(1.15) loggamma(1.15+0im) -0.06930620867104688224241731415650307100375642207340564554
88+
@test logabsgamma(0.89)[1] loggamma(0.89+0im) 0.074022173958081423702265889979810658434235008344573396963
89+
@test loggamma(0.91) loggamma(0.91+0im) 0.058922567623832379298241751183907077883592982094770449167
90+
@test loggamma(0.01) loggamma(0.01+0im) 4.599479878042021722513945411008748087261001413385289652419
91+
@test loggamma(-3.4-0.1im) -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im
92+
@test loggamma(-13.4-0.1im) -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im
93+
@test loggamma(-13.4+0.0im) conj(loggamma(-13.4-0.0im)) -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im
94+
@test loggamma(-13.4+8im) -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im
95+
@test logabsgamma(1+exp2(-20))[1] loggamma(1+exp2(-20)+0im) -5.504750066148866790922434423491111098144565651836914e-7
96+
@test loggamma(1+exp2(-20)+exp2(-19)*im) -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im
97+
@test loggamma(-300+2im) -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im
98+
@test loggamma(300+2im) 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im
99+
@test loggamma(1-6im) -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im
100+
@test loggamma(1-8im) -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im
101+
@test loggamma(1+6.5im) conj(loggamma(1-6.5im)) -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im
102+
@test loggamma(1+1im) conj(loggamma(1-1im)) -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im
103+
@test loggamma(-pi*1e7 + 6im) -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im
104+
@test loggamma(-pi*1e7 + 8im) -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im
105+
@test loggamma(-pi*1e14 + 6im) -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im
106+
@test loggamma(-pi*1e14 + 8im) -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im
107+
@test loggamma(2.05 + 0.03im) conj(loggamma(2.05 - 0.03im)) 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im
108+
@test loggamma(2+exp2(-20)+exp2(-19)*im) 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im
109+
end
108110

109111
@testset "loggamma for non-finite arguments" begin
110112
@test loggamma(Inf + 0im) === Inf + 0im
@@ -118,6 +120,38 @@
118120
@test loggamma(Inf + Inf*im) === loggamma(NaN + 0im) === loggamma(NaN*im) === NaN + NaN*im
119121
end
120122

123+
@testset "BigFloat" begin
124+
# test cases (taken from WolframAlpha, computed to 78 digits ≈ 256 bits)
125+
@test loggamma(big"3.099") big"0.78641374690055805872066586017892360313412585445116886979672329071050823224651" rtol=1e-75
126+
@test loggamma(big"1.15") big"-0.06930620867104688224241731415650307100375642207340564554412494594148673455871" rtol=1e-75
127+
@test logabsgamma(big"0.89")[1] big"0.0740221739580814237022658899798106584342350083445733969634566129726762260738245" rtol=1e-75
128+
@test loggamma(big"0.91") big"0.0589225676238323792982417511839070778835929820947704491677379048793029707373113" rtol=1e-75
129+
@test loggamma(big"0.01") big"4.59947987804202172251394541100874808726100141338528965241917138771477998920321" rtol=1e-75
130+
@test loggamma(1 + exp2(big"-20.0")) big"-5.50475006614886679092243442349111109814456565183691425527816079744208067935466e-7" rtol=1e-75
131+
132+
# consistency
133+
@test loggamma(big(3124.0)) == log(gamma(big(3124.0)))
134+
@test loggamma(big(3124.0)) loggamma(3124.0)
135+
@test logabsgamma(big(3124.0)) == (loggamma(big(3124.0)), 1)
136+
@test logabsgamma(big(3124.0))[1] logabsgamma(3124.0)[1]
137+
138+
# promotions
139+
@test loggamma(big(3124)) == loggamma(big(3124.0))
140+
@test loggamma(big(3//2)) == loggamma(big(1.5))
141+
@test logabsgamma(big(3124)) == logabsgamma(big(3124.0))
142+
@test logabsgamma(big(3//2)) == logabsgamma(big(1.5))
143+
144+
# negative values
145+
@test loggamma(big(-3.0)) == big(Inf)
146+
@test_throws DomainError loggamma(big(-2.76))
147+
148+
# non-finite values
149+
@test isnan(loggamma(big(NaN)))
150+
@test isnan(logabsgamma(big(NaN))[1])
151+
@test loggamma(big(Inf)) == big(Inf)
152+
@test logabsgamma(big(Inf))[1] == big(Inf)
153+
end
154+
121155
@testset "Other float types" begin
122156
struct NotAFloat <: AbstractFloat
123157
end

0 commit comments

Comments
 (0)