From 889429a6c45cb107e3f3c2338a7b61ce2b9bc2a9 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 29 Sep 2021 01:43:41 +0200 Subject: [PATCH 1/7] Use `mpfr_lngamma` --- src/gamma.jl | 10 ++++++++-- test/gamma.jl | 11 +++++++++++ 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index ab8c4eed..f8d9f48f 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,18 @@ 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) + z = BigFloat() + ccall((:mpfr_lngamma, :libmpfr), Cint, (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), z, x, ROUNDING_MODE[]) + return z +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..3fb48ba9 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -80,6 +80,17 @@ @test_throws MethodError logfactorial(1.0) end + @testset "BigFloat" begin + loggamma_mpfr(x::BigFloat) = invoke(SpecialFunctions._loggamma, Tuple{BigFloat}, x) + logabsgamma_mpfr(x::BigFloat) = invoke(SpecialFunctions._logabsgamma, Tuple{BigFloat}, x) + @test loggamma(big(3124)) == loggamma_mpfr(big(3124.0)) + @test loggamma(big(3124)) ≈ loggamma(3124) + @test loggamma(big(1e6)) == loggamma_mpfr(big(1e6)) + @test logabsgamma(big(1e7)) == logabsgamma_mpfr(big(1e7)) + @test logabsgamma(big(1e8)) == (loggamma(big(1e8)), 1) + @test logabsgamma(big(1e8))[1] ≈ logabsgamma(1e8)[1] + 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 From b76891f6ef48c9cbef47e4dd7c2d8c00ba32a2a4 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 29 Sep 2021 01:44:14 +0200 Subject: [PATCH 2/7] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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" From c3edd0a57753bb4141d04f1be1dd7bd9bad62b3f Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 29 Sep 2021 01:52:24 +0200 Subject: [PATCH 3/7] Rename variable --- src/gamma.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index f8d9f48f..b16f1e6d 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -650,9 +650,9 @@ function _loggamma(x::Real) end function _loggamma(x::BigFloat) - z = BigFloat() - ccall((:mpfr_lngamma, :libmpfr), Cint, (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), z, x, ROUNDING_MODE[]) - return z + y = BigFloat() + ccall((:mpfr_lngamma, :libmpfr), Cint, (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), y, x, ROUNDING_MODE[]) + return y end # Compute the logΓ(z) function using a combination of the asymptotic series, From 6470467de49f0c8c54cb8da18431656eb60b8129 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 29 Sep 2021 09:39:37 +0200 Subject: [PATCH 4/7] Throw error if `gamma(x) < 0` --- src/gamma.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/gamma.jl b/src/gamma.jl index b16f1e6d..c0cb034d 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -650,8 +650,10 @@ function _loggamma(x::Real) 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 From a4ccd8ca54d16b0476b91f3946de8c2352e6327d Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 29 Sep 2021 09:40:10 +0200 Subject: [PATCH 5/7] Add more tests --- test/gamma.jl | 93 ++++++++++++++++++++++++++++++++------------------- 1 file changed, 58 insertions(+), 35 deletions(-) diff --git a/test/gamma.jl b/test/gamma.jl index 3fb48ba9..878e5558 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -80,43 +80,34 @@ @test_throws MethodError logfactorial(1.0) end - @testset "BigFloat" begin - loggamma_mpfr(x::BigFloat) = invoke(SpecialFunctions._loggamma, Tuple{BigFloat}, x) - logabsgamma_mpfr(x::BigFloat) = invoke(SpecialFunctions._logabsgamma, Tuple{BigFloat}, x) - @test loggamma(big(3124)) == loggamma_mpfr(big(3124.0)) - @test loggamma(big(3124)) ≈ loggamma(3124) - @test loggamma(big(1e6)) == loggamma_mpfr(big(1e6)) - @test logabsgamma(big(1e7)) == logabsgamma_mpfr(big(1e7)) - @test logabsgamma(big(1e8)) == (loggamma(big(1e8)), 1) - @test logabsgamma(big(1e8))[1] ≈ logabsgamma(1e8)[1] + # 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 - # 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 - @testset "loggamma for non-finite arguments" begin @test loggamma(Inf + 0im) === Inf + 0im @test loggamma(Inf - 0.0im) === Inf - 0.0im @@ -129,6 +120,38 @@ @test loggamma(Inf + Inf*im) === loggamma(NaN + 0im) === loggamma(NaN*im) === NaN + NaN*im end + @testset "BigFloat" begin + # test cases (values taken from WolframAlpha) + @test loggamma(big(3.099)) ≈ 0.786413746900558058720665860178923603134125854451168869796 + @test loggamma(big(1.15)) ≈ -0.06930620867104688224241731415650307100375642207340564554 + @test logabsgamma(big(0.89))[1] ≈ 0.074022173958081423702265889979810658434235008344573396963 + @test loggamma(big(0.91)) ≈ 0.058922567623832379298241751183907077883592982094770449167 + @test loggamma(big(0.01)) ≈ 4.599479878042021722513945411008748087261001413385289652419 + @test loggamma(1 + exp2(big(-20.0))) ≈ -5.504750066148866790922434423491111098144565651836914e-7 + + # 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 From aeb8e9720454cb05afd5dd635e4d5c85c1b3143d Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 29 Sep 2021 09:41:31 +0200 Subject: [PATCH 6/7] Fix whitespace --- test/gamma.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/gamma.jl b/test/gamma.jl index 878e5558..edee7d35 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -80,7 +80,7 @@ @test_throws MethodError logfactorial(1.0) end - # values taken from Wolfram Alpha + # 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 From 35c1cdd3c09834b55a76ba8077d12dd7ddd1f77b Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 29 Sep 2021 23:29:20 +0200 Subject: [PATCH 7/7] Improve tests --- test/gamma.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/gamma.jl b/test/gamma.jl index edee7d35..15cd4716 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -121,13 +121,13 @@ end @testset "BigFloat" begin - # test cases (values taken from WolframAlpha) - @test loggamma(big(3.099)) ≈ 0.786413746900558058720665860178923603134125854451168869796 - @test loggamma(big(1.15)) ≈ -0.06930620867104688224241731415650307100375642207340564554 - @test logabsgamma(big(0.89))[1] ≈ 0.074022173958081423702265889979810658434235008344573396963 - @test loggamma(big(0.91)) ≈ 0.058922567623832379298241751183907077883592982094770449167 - @test loggamma(big(0.01)) ≈ 4.599479878042021722513945411008748087261001413385289652419 - @test loggamma(1 + exp2(big(-20.0))) ≈ -5.504750066148866790922434423491111098144565651836914e-7 + # 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)))