diff --git a/Project.toml b/Project.toml index e0d9275d..205cf836 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "SpecialFunctions" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "1.7.0" +version = "1.8.0" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/gamma.jl b/src/gamma.jl index ab8c4eed..c0cb034d 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -1,6 +1,6 @@ # This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license -using Base.MPFR: ROUNDING_MODE +using Base.MPFR: MPFRRoundingMode, ROUNDING_MODE export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, logabsbinomial @@ -643,12 +643,20 @@ See also [`logabsgamma`](@ref) for real `x`. """ loggamma(x::Number) = _loggamma(float(x)) -function _loggamma(x::Union{Real,BigFloat}) +function _loggamma(x::Real) (y, s) = logabsgamma(x) s < 0 && throw(DomainError(x, "`gamma(x)` must be non-negative")) return y end +function _loggamma(x::BigFloat) + isnan(x) && return x + y = BigFloat() + ccall((:mpfr_lngamma, :libmpfr), Cint, (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), y, x, ROUNDING_MODE[]) + isnan(y) && throw(DomainError(x, "`gamma(x)` must be non-negative")) + return y +end + # Compute the logΓ(z) function using a combination of the asymptotic series, # the Taylor series around z=1 and z=2, the reflection formula, and the shift formula. # Many details of these techniques are discussed in D. E. G. Hare, diff --git a/test/gamma.jl b/test/gamma.jl index 8f49290f..15cd4716 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -80,31 +80,33 @@ @test_throws MethodError logfactorial(1.0) end - # loggamma & logabsgamma test cases (from Wolfram Alpha) - @test loggamma(-300im) ≅ -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im - @test loggamma(3.099) ≅ loggamma(3.099+0im) ≅ 0.786413746900558058720665860178923603134125854451168869796 - @test loggamma(1.15) ≅ loggamma(1.15+0im) ≅ -0.06930620867104688224241731415650307100375642207340564554 - @test logabsgamma(0.89)[1] ≅ loggamma(0.89+0im) ≅ 0.074022173958081423702265889979810658434235008344573396963 - @test loggamma(0.91) ≅ loggamma(0.91+0im) ≅ 0.058922567623832379298241751183907077883592982094770449167 - @test loggamma(0.01) ≅ loggamma(0.01+0im) ≅ 4.599479878042021722513945411008748087261001413385289652419 - @test loggamma(-3.4-0.1im) ≅ -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im - @test loggamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im - @test loggamma(-13.4+0.0im) ≅ conj(loggamma(-13.4-0.0im)) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im - @test loggamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im - @test logabsgamma(1+exp2(-20))[1] ≅ loggamma(1+exp2(-20)+0im) ≅ -5.504750066148866790922434423491111098144565651836914e-7 - @test loggamma(1+exp2(-20)+exp2(-19)*im) ≅ -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im - @test loggamma(-300+2im) ≅ -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im - @test loggamma(300+2im) ≅ 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im - @test loggamma(1-6im) ≅ -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im - @test loggamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im - @test loggamma(1+6.5im) ≅ conj(loggamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im - @test loggamma(1+1im) ≅ conj(loggamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im - @test loggamma(-pi*1e7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im - @test loggamma(-pi*1e7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im - @test loggamma(-pi*1e14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im - @test loggamma(-pi*1e14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im - @test loggamma(2.05 + 0.03im) ≅ conj(loggamma(2.05 - 0.03im)) ≅ 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im - @test loggamma(2+exp2(-20)+exp2(-19)*im) ≅ 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im + # values taken from Wolfram Alpha + @testset "loggamma & logabsgamma test cases" begin + @test loggamma(-300im) ≅ -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im + @test loggamma(3.099) ≅ loggamma(3.099+0im) ≅ 0.786413746900558058720665860178923603134125854451168869796 + @test loggamma(1.15) ≅ loggamma(1.15+0im) ≅ -0.06930620867104688224241731415650307100375642207340564554 + @test logabsgamma(0.89)[1] ≅ loggamma(0.89+0im) ≅ 0.074022173958081423702265889979810658434235008344573396963 + @test loggamma(0.91) ≅ loggamma(0.91+0im) ≅ 0.058922567623832379298241751183907077883592982094770449167 + @test loggamma(0.01) ≅ loggamma(0.01+0im) ≅ 4.599479878042021722513945411008748087261001413385289652419 + @test loggamma(-3.4-0.1im) ≅ -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im + @test loggamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im + @test loggamma(-13.4+0.0im) ≅ conj(loggamma(-13.4-0.0im)) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im + @test loggamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im + @test logabsgamma(1+exp2(-20))[1] ≅ loggamma(1+exp2(-20)+0im) ≅ -5.504750066148866790922434423491111098144565651836914e-7 + @test loggamma(1+exp2(-20)+exp2(-19)*im) ≅ -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im + @test loggamma(-300+2im) ≅ -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im + @test loggamma(300+2im) ≅ 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im + @test loggamma(1-6im) ≅ -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im + @test loggamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im + @test loggamma(1+6.5im) ≅ conj(loggamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im + @test loggamma(1+1im) ≅ conj(loggamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im + @test loggamma(-pi*1e7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im + @test loggamma(-pi*1e7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im + @test loggamma(-pi*1e14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im + @test loggamma(-pi*1e14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im + @test loggamma(2.05 + 0.03im) ≅ conj(loggamma(2.05 - 0.03im)) ≅ 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im + @test loggamma(2+exp2(-20)+exp2(-19)*im) ≅ 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im + end @testset "loggamma for non-finite arguments" begin @test loggamma(Inf + 0im) === Inf + 0im @@ -118,6 +120,38 @@ @test loggamma(Inf + Inf*im) === loggamma(NaN + 0im) === loggamma(NaN*im) === NaN + NaN*im end + @testset "BigFloat" begin + # test cases (taken from WolframAlpha, computed to 78 digits ≈ 256 bits) + @test loggamma(big"3.099") ≈ big"0.78641374690055805872066586017892360313412585445116886979672329071050823224651" rtol=1e-75 + @test loggamma(big"1.15") ≈ big"-0.06930620867104688224241731415650307100375642207340564554412494594148673455871" rtol=1e-75 + @test logabsgamma(big"0.89")[1] ≈ big"0.0740221739580814237022658899798106584342350083445733969634566129726762260738245" rtol=1e-75 + @test loggamma(big"0.91") ≈ big"0.0589225676238323792982417511839070778835929820947704491677379048793029707373113" rtol=1e-75 + @test loggamma(big"0.01") ≈ big"4.59947987804202172251394541100874808726100141338528965241917138771477998920321" rtol=1e-75 + @test loggamma(1 + exp2(big"-20.0")) ≈ big"-5.50475006614886679092243442349111109814456565183691425527816079744208067935466e-7" rtol=1e-75 + + # consistency + @test loggamma(big(3124.0)) == log(gamma(big(3124.0))) + @test loggamma(big(3124.0)) ≈ loggamma(3124.0) + @test logabsgamma(big(3124.0)) == (loggamma(big(3124.0)), 1) + @test logabsgamma(big(3124.0))[1] ≈ logabsgamma(3124.0)[1] + + # promotions + @test loggamma(big(3124)) == loggamma(big(3124.0)) + @test loggamma(big(3//2)) == loggamma(big(1.5)) + @test logabsgamma(big(3124)) == logabsgamma(big(3124.0)) + @test logabsgamma(big(3//2)) == logabsgamma(big(1.5)) + + # negative values + @test loggamma(big(-3.0)) == big(Inf) + @test_throws DomainError loggamma(big(-2.76)) + + # non-finite values + @test isnan(loggamma(big(NaN))) + @test isnan(logabsgamma(big(NaN))[1]) + @test loggamma(big(Inf)) == big(Inf) + @test logabsgamma(big(Inf))[1] == big(Inf) + end + @testset "Other float types" begin struct NotAFloat <: AbstractFloat end