diff --git a/benchmark/gamma_benchmark.jl b/benchmark/gamma_benchmark.jl new file mode 100644 index 0000000..e91daf5 --- /dev/null +++ b/benchmark/gamma_benchmark.jl @@ -0,0 +1,74 @@ +using BenchmarkTools +using Bessels +using Bessels.GammaFunctions: loggamma, logabsgamma +using ArbNumerics +import SpecialFunctions + +setworkingprecision(ArbFloat, 500) +setextrabits(128) + + + +# compute relative errors: abs(1 - Bessels.f(x)/ArbNumerics.f(ArbFloat(x))) +# 1000 random points in (0, 100), report mean and max + +println("=== loggamma Float64 accuracy ===") +x_vals = rand(1000) .* 100 +errs = [abs(1 - loggamma(x) / Float64(ArbNumerics.lgamma(ArbFloat(x)))) for x in x_vals] +println(" mean relative error: ", sum(errs)/length(errs)) +println(" max relative error: ", maximum(errs)) + +println("\n=== loggamma ComplexF64 accuracy ===") +z_vals = [rand()*100 + rand()*100*im for _ in 1:1000] +errs_c = Float64[] +for z in z_vals + ref = Complex{Float64}(ArbNumerics.lgamma(ArbFloat(real(z)) + ArbFloat(imag(z))*im)) + push!(errs_c, abs(1 - loggamma(z) / ref)) +end +println(" mean relative error: ", sum(errs_c)/length(errs_c)) +println(" max relative error: ", maximum(errs_c)) + +println("\n=== logabsgamma Float64 accuracy ===") +errs_la = Float64[] +for x in x_vals + ref = Float64(ArbNumerics.lgamma(ArbFloat(x))) + y, _ = logabsgamma(x) + push!(errs_la, abs(1 - y / ref)) +end +println(" mean relative error: ", sum(errs_la)/length(errs_la)) +println(" max relative error: ", maximum(errs_la)) + + + + + +# performance benchmarks + +suite = BenchmarkGroup() + +suite["loggamma"] = BenchmarkGroup() +suite["loggamma"]["Bessels_Float64"] = @benchmarkable loggamma(x) setup=(x = rand()*100) +suite["loggamma"]["SpecialFunctions_Float64"] = @benchmarkable SpecialFunctions.loggamma(x) setup=(x = rand()*100) +suite["loggamma"]["Bessels_ComplexF64"] = @benchmarkable loggamma(z) setup=(z = rand(ComplexF64)*100) +suite["loggamma"]["SpecialFunctions_ComplexF64"] = @benchmarkable SpecialFunctions.loggamma(z) setup=(z = rand(ComplexF64)*100) + +suite["logabsgamma"] = BenchmarkGroup() +suite["logabsgamma"]["Bessels_Float64"] = @benchmarkable logabsgamma(x) setup=(x = rand()*100) +suite["logabsgamma"]["SpecialFunctions_Float64"] = @benchmarkable SpecialFunctions.logabsgamma(x) setup=(x = rand()*100) + +results = run(suite, verbose=true) + +println("\n=== loggamma Float64 timing ===") +println("Bessels: ", results["loggamma"]["Bessels_Float64"]) +println("SpecialFunctions: ", results["loggamma"]["SpecialFunctions_Float64"]) +println("Speedup (Bessels): ", round(median(results["loggamma"]["SpecialFunctions_Float64"]).time / median(results["loggamma"]["Bessels_Float64"]).time, digits=2), "x") + +println("\n=== loggamma ComplexF64 timing ===") +println("Bessels: ", results["loggamma"]["Bessels_ComplexF64"]) +println("SpecialFunctions: ", results["loggamma"]["SpecialFunctions_ComplexF64"]) +println("Speedup (Bessels): ", round(median(results["loggamma"]["SpecialFunctions_ComplexF64"]).time / median(results["loggamma"]["Bessels_ComplexF64"]).time, digits=2), "x") + +println("\n=== logabsgamma Float64 timing ===") +println("Bessels: ", results["logabsgamma"]["Bessels_Float64"]) +println("SpecialFunctions: ", results["logabsgamma"]["SpecialFunctions_Float64"]) +println("Speedup (Bessels): ", round(median(results["logabsgamma"]["SpecialFunctions_Float64"]).time / median(results["logabsgamma"]["Bessels_Float64"]).time, digits=2), "x") diff --git a/src/GammaFunctions/GammaFunctions.jl b/src/GammaFunctions/GammaFunctions.jl index 39a0399..f15e2d2 100644 --- a/src/GammaFunctions/GammaFunctions.jl +++ b/src/GammaFunctions/GammaFunctions.jl @@ -3,7 +3,9 @@ module GammaFunctions using ..Math: SQ2PI export gamma +export loggamma, logabsgamma, logfactorial include("gamma.jl") +include("loggamma.jl") end diff --git a/src/GammaFunctions/gamma.jl b/src/GammaFunctions/gamma.jl index 9a9df75..0f00490 100644 --- a/src/GammaFunctions/gamma.jl +++ b/src/GammaFunctions/gamma.jl @@ -114,4 +114,7 @@ function gamma(n::Integer) @inbounds return Float64(factorial(n-1)) end -gamma_near_1(x) = evalpoly(x-one(x), (1.0, -0.5772156649015329, 0.9890559953279725, -0.23263776388631713)) \ No newline at end of file +gamma_near_1(x) = evalpoly(x-one(x), (1.0, -0.5772156649015329, 0.9890559953279725, -0.23263776388631713)) + +gamma(x::BigFloat) = exp(loggamma(x)) +gamma(x::Complex) = exp(loggamma(x)) \ No newline at end of file diff --git a/src/GammaFunctions/loggamma.jl b/src/GammaFunctions/loggamma.jl new file mode 100644 index 0000000..a9619cf --- /dev/null +++ b/src/GammaFunctions/loggamma.jl @@ -0,0 +1,372 @@ +# Pure Julia loggamma and logabsgamma implementations +# Partly adapted from SpecialFunctions.jl (MIT license) +# using Stirling asymptotic series, Taylor series at z=1 and z=2, +# reflection formula, and shift recurrence. +# See: D. E. G. Hare, "Computing the principal branch of log-Gamma," +# J. Algorithms 25, pp. 221-236 (1997) + +""" + loggamma(x::Real) + +Returns the log of the absolute value of ``\\Gamma(x)`` for real `x`. +Throws a `DomainError` if ``\\Gamma(x)`` is negative. + +For complex arguments, `exp(loggamma(x))` matches `gamma(x)` up to floating-point error +but may differ from `log(gamma(x))` by an integer multiple of ``2\\pi i``. + +External links: [DLMF](https://dlmf.nist.gov/5.4), [Wikipedia](https://en.wikipedia.org/wiki/Gamma_function#The_log-gamma_function) +""" +loggamma(x::Float64) = _loggamma(x) +loggamma(x::Union{Float16, Float32}) = typeof(x)(_loggamma(Float64(x))) +loggamma(x::Rational) = loggamma(float(x)) +loggamma(x::Integer) = loggamma(float(x)) +loggamma(z::Complex{Float64}) = _loggamma(z) +loggamma(z::Complex{Float32}) = Complex{Float32}(_loggamma(Complex{Float64}(z))) +loggamma(z::Complex{Float16}) = Complex{Float16}(_loggamma(Complex{Float64}(z))) +loggamma(z::Complex{<:Integer}) = _loggamma(Complex{Float64}(z)) +loggamma(z::Complex{<:Rational}) = loggamma(float(z)) +function loggamma(x::BigFloat) + if isnan(x) + return x + elseif isinf(x) + return x > 0 ? x : BigFloat(NaN) + elseif x <= 0 + x == 0 && return BigFloat(Inf) + isinteger(x) && return BigFloat(Inf) # negative integer pole + y, sgn = _logabsgamma(x) + sgn < 0 && throw(DomainError(x, "`gamma(x)` must be non-negative")) + return y + end + return real(_loggamma_complex_bigfloat(Complex{BigFloat}(x, zero(BigFloat)))) +end +loggamma(z::Complex{BigFloat}) = _loggamma(z) + +""" + logfactorial(x) + +Compute the logarithmic factorial of a nonnegative integer `x` via loggamma. +""" +logfactorial(x::Integer) = x < 0 ? throw(DomainError(x, "`x` must be non-negative.")) : loggamma(float(x + oneunit(x))) + +""" + logabsgamma(x::Real) + +Returns a tuple `(log(abs(Γ(x))), sign(Γ(x)))` for real `x`. +""" +logabsgamma(x::Real) = _logabsgamma(float(x)) + +const HALF_LOG2PI_F64 = 9.1893853320467274178032927e-01 +const LOGPI_F64 = 1.1447298858494002 +const TWO_PI_F64 = 6.2831853071795864769252842 + +# Stirling asymptotic series for log(Γ(x)), valid for x > 0 sufficiently large +# coefficients are bernoulli[2k] / (2k*(2k-1)) for k = 1,...,8 +function _loggamma_stirling(x::Float64) + t = inv(x) + w = t * t + return muladd(x - 0.5, log(x), -x + HALF_LOG2PI_F64 + # log(2π)/2 + t * @evalpoly(w, + 8.333333333333333333333368e-02, -2.777777777777777777777778e-03, + 7.936507936507936507936508e-04, -5.952380952380952380952381e-04, + 8.417508417508417508417510e-04, -1.917526917526917526917527e-03, + 6.410256410256410256410257e-03, -2.955065359477124183006535e-02 + ) + ) +end + +# Asymptotic series for log(Γ(z)) for complex z with sufficiently large real(z) or |imag(z)| +function _loggamma_asymptotic(z::Complex{Float64}) + zinv = inv(z) + t = zinv * zinv + return (z - 0.5) * log(z) - z + HALF_LOG2PI_F64 + # log(2π)/2 + zinv * @evalpoly(t, + 8.333333333333333333333368e-02, -2.777777777777777777777778e-03, + 7.936507936507936507936508e-04, -5.952380952380952380952381e-04, + 8.417508417508417508417510e-04, -1.917526917526917526917527e-03, + 6.410256410256410256410257e-03, -2.955065359477124183006535e-02 + ) +end + +function _logabsgamma(x::Float64) + if isnan(x) + return x, 1 + elseif x > 0 + return _loggamma_unsafe_pos(x), 1 + elseif x == 0 + return Inf, Int(sign(1/x)) # ±0 → correct sign + else + # reflection formula: Γ(x) = π / (sin(πx) * Γ(1-x)) + s = sinpi(x) + s == 0 && return Inf, 1 + sgn = signbit(s) ? -1 : 1 + return LOGPI_F64 - log(abs(s)) - _loggamma(1.0 - x), sgn + end +end + +# Version of logabsgamma for Float64 without input checks, used for loggamma reduction to avoid double checks +function _logabsgamma_unsafe_sub0(x::Float64) + # Since this is only used from loggamma, we can assume x < 0 is a clean input + s = sinpi(x) + s == 0 && return Inf, 1 + sgn = signbit(s) ? -1 : 1 + return LOGPI_F64 - log(abs(s)) - _loggamma(1.0 - x), sgn +end + +function _logabsgamma(x::Float32) + y, s = _logabsgamma(Float64(x)) + return Float32(y), s +end + +function _logabsgamma(x::Float16) + y, s = _logabsgamma(Float64(x)) + return Float16(y), s +end + +function _logabsgamma(x::BigFloat) + if isnan(x) + return x, 1 + elseif isinf(x) + return x > 0 ? (x, 1) : (BigFloat(NaN), 1) + elseif x > 0 + return real(_loggamma_complex_bigfloat(Complex{BigFloat}(x, zero(BigFloat)))), 1 + elseif iszero(x) + return BigFloat(Inf), Int(sign(1 / x)) + end + + s = sinpi(x) + s == 0 && return BigFloat(Inf), 1 + return real(_loggamma_complex_bigfloat(Complex{BigFloat}(x, zero(BigFloat)))), (signbit(s) ? -1 : 1) +end + +# Lanczos-type rational approximation for loggamma on (2, 3) +# Used as the core for reduction-based approach +const _LOGGAMMA_P = ( + -2.44167345903529816830968e-01, 6.73523010531981020863696e-02, + -2.05808084277845478790009e-02, 7.38555102867398526627303e-03, + -2.89051033074153369901384e-03, 1.19275391170326097711398e-03, + -5.09669524743042422335582e-04, 2.23154759903498081132513e-04, + -9.94575127818085337147321e-05, 4.49262367382046739858373e-05, + -2.05077312586603517590604e-05 +) + +# loggamma for real Float64 +function _loggamma(x::Float64) + if isnan(x) + return x + elseif isinf(x) + return x > 0 ? Inf : NaN + elseif x ≤ 0 + x == 0 && return Inf + isinteger(x) && return Inf # negative integer pole + # reflection: log|Γ(x)| = log(π) - log|sin(πx)| - log(Γ(1-x)) + # but loggamma for real requires Γ(x)>0 + y, sgn = _logabsgamma_unsafe_sub0(x) + sgn < 0 && throw(DomainError(x, "`gamma(x)` must be non-negative")) + return y + elseif x < 7 + # shift x into asymptotic region [7,∞) + n = 7 - floor(Int, x) + z = x + prod = one(x) + for i in 0:n-1 + prod *= z + i + end + return _loggamma_stirling(z + n) - log(prod) + else + return _loggamma_stirling(x) + end +end + +# loggamma for positive x without checks, used for logabsgamma reduction to avoid double checks +function _loggamma_unsafe_pos(x::Float64) + if x < 7 + # shift x into asymptotic region [7,∞) + n = 7 - floor(Int, x) + z = x + prod = one(x) + for i in 0:n-1 + prod *= z + i + end + return _loggamma_stirling(z + n) - log(prod) + else + return _loggamma_stirling(x) + end +end + +# Complex loggamma for Float64 +# Combines the asymptotic series, Taylor series at z=1 and z=2, +# the reflection formula, and the shift recurrence. +function _loggamma(z::Complex{Float64}) + x, y = reim(z) + yabs = abs(y) + + if !isfinite(x) || !isfinite(y) + if isinf(x) && isfinite(y) + return Complex(x, x > 0 ? (y == 0 ? y : copysign(Inf, y)) : copysign(Inf, -y)) + elseif isfinite(x) && isinf(y) + return Complex(-Inf, y) + else + return Complex(NaN, NaN) + end + elseif x > 7 || yabs > 7 + return _loggamma_asymptotic(z) + elseif x < 0.1 + if x == 0 && y == 0 + return Complex(Inf, signbit(x) ? copysign(Float64(π), -y) : -y) + end + # reflection formula with correct branch cut. + return Complex(LOGPI_F64, copysign(TWO_PI_F64, y) * floor(0.5 * x + 0.25)) - + log(sinpi(z)) - _loggamma(1 - z) + elseif abs(x - 1) + yabs < 0.1 + # Taylor series at z=1 + # coefficients: [-γ; [(-1)^k * ζ(k)/k for k in 2:15]] + w = Complex(x - 1, y) + return w * @evalpoly(w, + -5.7721566490153286060651188e-01, 8.2246703342411321823620794e-01, + -4.0068563438653142846657956e-01, 2.705808084277845478790009e-01, + -2.0738555102867398526627303e-01, 1.6955717699740818995241986e-01, + -1.4404989676884611811997107e-01, 1.2550966952474304242233559e-01, + -1.1133426586956469049087244e-01, 1.000994575127818085337147e-01, + -9.0954017145829042232609344e-02, 8.3353840546109004024886499e-02, + -7.6932516411352191472827157e-02, 7.1432946295361336059232779e-02, + -6.6668705882420468032903454e-02 + ) + elseif abs(x - 2) + yabs < 0.1 + # Taylor series at z=2 + # coefficients: [1-γ; [(-1)^k * (ζ(k)-1)/k for k in 2:12]] + w = Complex(x - 2, y) + return w * @evalpoly(w, + 4.2278433509846713939348812e-01, 3.2246703342411321823620794e-01, + -6.7352301053198095133246196e-02, 2.0580808427784547879000897e-02, + -7.3855510286739852662729527e-03, 2.8905103307415232857531201e-03, + -1.1927539117032609771139825e-03, 5.0966952474304242233558822e-04, + -2.2315475845357937976132853e-04, 9.9457512781808533714662972e-05, + -4.4926236738133141700224489e-05, 2.0507212775670691553131246e-05 + ) + else + # shift using recurrence: loggamma(z) = loggamma(z+n) - log(∏(z+k)) + shiftprod = Complex(x, yabs) + x += 1 + sb = false + signflips = 0 + while x ≤ 7 + shiftprod *= Complex(x, yabs) + sb′ = signbit(imag(shiftprod)) + signflips += sb′ & (sb′ != sb) + sb = sb′ + x += 1 + end + shift = log(shiftprod) + if signbit(y) + shift = Complex(real(shift), signflips * -TWO_PI_F64 - imag(shift)) + else + shift = Complex(real(shift), imag(shift) + signflips * TWO_PI_F64) + end + return _loggamma_asymptotic(Complex(x, y)) - shift + end +end + +# Complex BigFloat loggamma +# Adapted from SpecialFunctions.jl (MIT license) +# Uses Stirling series with Bernoulli numbers computed via Akiyama-Tanigawa, +# reflection formula, upward recurrence, and branch correction via Float64 oracle. + +# Scaled Stirling coefficients B_{2k}/(2k*(2k-1)) * zr^(1-2k) for k=0,...,n +# Bernoulli numbers computed inline via the Akiyama-Tanigawa algorithm +function _scaled_stirling_coeffs(n::Integer, zr::Complex{BigFloat}) + mmax = 2n + A = Vector{Rational{BigInt}}(undef, mmax + 1) + E = Vector{Complex{BigFloat}}(undef, n + 1) + @inbounds for m = 0:mmax + A[m+1] = 1 // (m + 1) + for j = m:-1:1 + A[j] = j * (A[j] - A[j+1]) + end + if iseven(m) + k = m ÷ 2 + E[k+1] = A[1] / (2k * (2k - 1)) * zr^(1 - 2k) + end + end + return E +end + +function _loggamma_complex_bigfloat(z::Complex{BigFloat}) + bigpi = big(π) + x = real(z) + y = imag(z) + + if !isfinite(x) || !isfinite(y) + inf = BigFloat(Inf) + nan = BigFloat(NaN) + if isinf(x) && isfinite(y) + yim = x > 0 ? (iszero(y) ? y : copysign(inf, y)) : copysign(inf, -y) + return Complex{BigFloat}(x, yim) + elseif isfinite(x) && isinf(y) + return Complex{BigFloat}(-inf, y) + else + return Complex{BigFloat}(nan, nan) + end + end + + # reflection formula + if x < 0.5 + val = log(bigpi) - log(sinpi(z)) - _loggamma_complex_bigfloat(1 - z) + return _loggamma_branchcorrect(val, z) + end + + # upward recurrence: shift z into the Stirling region + p = precision(BigFloat) + r = max(0, Int(ceil(p - abs(z)))) + zr = z + r + + # Stirling series + N = max(10, p ÷ 15) + B = _scaled_stirling_coeffs(N, zr) + lg = sum(B[2:end]) + (zr - big"0.5") * log(zr) - zr + log(sqrt(2 * bigpi)) + + # undo upward shift via log of product + if r > 0 + prodarg = prod(z + (i - 1) for i in 1:r) + lg -= log(prodarg) + end + + return _loggamma_branchcorrect(lg, z) +end + +# Branch correction: offset by multiples of 2πi to match the Float64 branch +function _loggamma_branchcorrect(val::Complex{BigFloat}, z::Complex{BigFloat}) + zf = _loggamma_oracle64_point(z) + val_f = _loggamma(zf) + imv = imag(val) + k = round(Int, (Float64(imv) - imag(val_f)) / (2π)) + return Complex{BigFloat}(real(val), imv - 2 * big(π) * k) +end + +# Map a BigFloat complex point to Float64 for branch-cut determination +function _loggamma_oracle64_point(z::Complex{BigFloat}) + xr = Float64(real(z)) + xi = Float64(imag(z)) + n = round(Int, xr) + if n ≤ 0 && isapprox(xr, Float64(n); atol=eps(Float64)) && abs(xi) ≤ 2eps(Float64) + xr = real(z) > n ? nextfloat(xr) : prevfloat(xr) + end + return Complex{Float64}(xr, xi) +end + +# Complex{BigFloat} entry point with guard precision +function _loggamma(z::Complex{BigFloat}) + imz = imag(z) + rez = real(z) + if iszero(imz) + return Complex(loggamma(rez)) + end + p0 = precision(BigFloat) + guard = 16 + setprecision(p0 + guard) do + zhi = Complex{BigFloat}(rez, imz) + rhi = _loggamma_complex_bigfloat(zhi) + setprecision(p0) do + return Complex{BigFloat}(real(rhi), imag(rhi)) + end + end +end diff --git a/test/gamma_test.jl b/test/gamma_test.jl index d7ff6bf..4358a5e 100644 --- a/test/gamma_test.jl +++ b/test/gamma_test.jl @@ -14,3 +14,282 @@ end x = [0, 1, 2, 3, 8, 15, 20, 30] @test SpecialFunctions.gamma.(x) ≈ Bessels.gamma.(x) + +# loggamma tests +using Bessels.GammaFunctions: gamma, gamma_near_1, loggamma, logabsgamma, logfactorial, _loggamma_oracle64_point + +# exact at expansion point (constant term = Γ(1) = 1) +@test gamma_near_1(1.0) === 1.0 +# nearby values: degree-3 Taylor gives O((x-1)^4) error, +# ≲ 1e-4 for |x-1| ≤ 0.1 +for x in (0.90, 0.95, 1.0, 1.05, 1.10) + @test gamma_near_1(x) ≈ gamma(x) atol=1e-3 +end + +@testset "logfactorial" begin + @test logfactorial(0) ≈ 0.0 atol=5eps(Float64) + for n in (1, 2, 3, 10, 50) + @test logfactorial(n) ≈ loggamma(Float64(n + 1)) + @test logfactorial(n) ≈ SpecialFunctions.loggamma(Float64(n + 1)) atol=10eps(Float64) + end + @test logfactorial(Int32(7)) ≈ loggamma(8.0) + @test_throws DomainError logfactorial(-1) +end + +# real loggamma for Float64, Float32, Float16 against SpecialFunctions.jl +# Note: Stirling-based approach has higher relative error near loggamma zeros (x ≈ 1, 2) +for (T, max, rtol) in ((Float16, 13, 1.0), (Float32, 43, 1.0), (Float64, 170, 7)) + v = rand(T, 5000) * max + for x in v + ref = T(SpecialFunctions.loggamma(widen(x))) + @test isapprox(ref, loggamma(x), atol=1e-10, rtol=rtol*eps(T)) + end + @test isnan(loggamma(T(NaN))) + @test loggamma(T(Inf)) == T(Inf) +end + +@test gamma(0.29384) ≈ exp(loggamma(0.29384)) +@test gamma(0.29384+0.12938im) ≈ exp(loggamma(0.29384+0.12938im)) + +# logabsgamma +for (T, rtol) in ((Float16, 1.0), (Float32, 1.0), (Float64, 7)) + for x in T[0.5, 1.0, 1.5, 2.0, 5.0, 10.0, 50.0, -0.5, -1.5, -2.5, -3.5, -10.5] + y1, s1 = logabsgamma(x) + y2, s2 = SpecialFunctions.logabsgamma(Float64(x)) + @test isapprox(y1, T(y2), atol=T(1e-3), rtol=rtol*eps(T)) + @test s1 == s2 + end +end + +# logabsgamma edge cases and SpecialFunctions behavior consistency +@test logabsgamma(0.0) == (Inf, 1) +@test logabsgamma(-0.0) == (Inf, -1) +@test logabsgamma(-1.0) == (Inf, 1) +@test logabsgamma(-2.0) == (Inf, 1) +@test isnan(logabsgamma(NaN)[1]) +@test logabsgamma(NaN)[2] == 1 +# real loggamma should throw for negative gamma +@test_throws DomainError loggamma(-0.5) +@test loggamma(-1.5) == logabsgamma(-1.5)[1] + +# complex loggamma for Float64 +for z in [1.0+1.0im, 2.0+0.5im, 0.5+3.0im, 5.0+2.0im, 0.1+0.1im, + -1.5+0.5im, -0.5+2.0im, 3.0+0.01im, 10.0+5.0im, 0.5+0.01im] + @test isapprox(loggamma(z), SpecialFunctions.loggamma(z), rtol=7*eps(Float64)) +end + +# complex loggamma for Float32 and Float16 +for z in [1.0f0+1.0f0im, 5.0f0+2.0f0im, 0.5f0+3.0f0im] + @test isapprox(loggamma(Complex{Float32}(z)), Complex{Float32}(SpecialFunctions.loggamma(Complex{Float64}(z))), rtol=eps(Float32)) + @test isapprox(loggamma(Complex{Float16}(z)), Complex{Float16}(SpecialFunctions.loggamma(Complex{Float64}(z))), rtol=eps(Float16)) +end + +# complex loggamma edge cases and SpecialFunctions consistency +@test loggamma(Complex(Inf, 0.0)) == Complex(Inf, 0.0) +@test all(isnan, reim(loggamma(Complex(NaN, NaN)))) +@test loggamma(Complex(0.0, 0.0)) === Complex(Inf, -0.0) +@test loggamma(Complex(-0.0, 0.0)) == Complex(Inf, -Float64(π)) + +# BigFloat loggamma (real via complex with zero imaginary part) +for x in [big"0.5", big"1.5", big"5.0", big"10.0", big"50.0"] + @test isapprox(loggamma(x), SpecialFunctions.loggamma(x), rtol=1e-30) +end +@test isapprox(loggamma(big"1.0"), big"0.0", atol=1e-60) +@test isapprox(loggamma(big"2.0"), big"0.0", atol=1e-60) + +# Map Complex{Int64} to Complex{Float64} for loggamma tests +@test loggamma(Complex{Int64}(-300)) ≈ loggamma(Complex{Float64}(-300)) + + # 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 "Complex{BigFloat}" begin + for p in (128, 256, 512) + setprecision(p) do + # loggamma test cases (taken from WolframAlpha) + @test loggamma(Complex{BigFloat}(big"1.4", big"3.7")) ≈ big"-3.70940253309968418984314134854809331185736805646502631204023135069388352598997710462597918070702727437506783302173663095854186723146887577228573371509941123007270190112210546396844695302684314798287935259522555883298038366536077208173987529577169748995380430158372715939689026101797721568006330211802429958418317776826022950333285605816823534143341838578285140019591030846382087557833758489229193473180731236704742771447789943609553100247685639021852545170948778374393632642723407449801941928273507786" + big"2.45680905027686511843270815755752124912922749263425191557239423638341066818045733889313969102751920720477160888088575007996966892611180451906681270423482044538930858698226109815063261456588626497471947149577237077662427442208195627430175449469651634280548198970860137426524308382796138277607215437522597736608829802065700440231980717557305458696328442743516723459412895697201557510117612009902085541415052308628636399246616849495938682489375024940706179414665940019738206000502042295770027318071467942"*im + @test loggamma(Complex{BigFloat}(big"-3.4", big"-0.1")) ≈ big"-1.17333533220647794810490885589189574408477150036591434542448532049403984398309615180598419193372394147554168405641323544591892928508095342822832004280366916047823075044462434956352385622175279286874520918371463520339183833036366580824529184504285582306128824924251346911575648934153502523354232665207345984362970323798726840265637880699099675010674509364494927004863383425318768708942281798027350435165386600514393768950581065682200893415680265730090421639141142873001336343605907188519518351947777089750716521425258995341843362124378354032424181927250141434187420346538318521452730170494364862200575165169281864878658796523132562054308721448489115739548477753359680696950183723903855692026149555142451921517583039596054377943432700206626447104309600454668378148843293998091304506665118850678678831194150927534302652888726369436145442541195405764173555455442990408710747398811929878918982539083715025044202675895474356318799592540065071048146945557430938208873062999453334736614448777192367740242361544911" + big"12.3314655012478268428755861044159800943162689746718192812954820421814925428937569884922482899922182108501683937660262759142172101775968893536424040892149708658863325903803507911939566037154489951252238647647402247597088855701451028588020604254079765523864277968914769300229513694667769438632836597735443816632344312013588575240574588529169102277436974026980920640009348503482380251876347214375511385088076700003044257500058872364330837264152707906065680970651314444536833041753145865573422902672476837825912427417183184497877656780666843913107467161665891469009396090946299606182430537112098946313562060208498550350090478607894628838869289194986670248766434102637364025088229859997978069923936222838027747044361676772501869744770389248222868789957322667506469525282107480146222663211761905584280946561123905319927661045637168771975340291467464330952607239039478342492953711672563497381915457648833579449139976570359273070513506747282021008512940692328868075470257049078667016119521218353123785031967384424" * im + @test loggamma(Complex{BigFloat}(big"1.0", big"6.5")) ≈ big"-8.35533650251135956898874979636340693034277901250481133076775315558468225236054299379804391077446370723114013332576474380468988528404564621934767485982173074710822143101705650155694961710942206758861407956691226902034077206603599116440305209875779596048383199145404442671523667694954031176386914820878412589001648171807491290458567807444100704632567345208428201008150473023561239072738489117252232075537269146942018683897902370439498725767025305446910337402590297901357999640802408423888317102332950502" + big"6.43928161597598339481129290184076602632280364914798257448626986803948974762094722413114781643500457429277396907129309037120052638434079511888202779346673590756680534496688216048661275875197218642757656165236731041690196648409171135474570376929222374881526540613424000396812702483748804403461430606647372072917822435446045663422222644485591648971382453371140638905009688610340367386497741974035080606633217160834318470835793011218439321046372790514195519666101799949551367645485542165382267523023837261"*im + # sanity check against Complex{Float64} values + @test Complex{Float64}(loggamma(Complex{BigFloat}(big"1.4", big"3.7"))) ≈ loggamma(1.4+3.7im) + @test Complex{Float64}(loggamma(Complex{BigFloat}(big"-3.4", big"-0.1"))) ≈ loggamma(-3.4-0.1im) + @test Complex{Float64}(loggamma(Complex{BigFloat}(big"1.0", big"6.5"))) ≈ loggamma(1.0+6.5im) + @test Complex{Float64}(gamma(Complex{BigFloat}(big"1.4", big"3.7"))) ≈ gamma(1.4+3.7im) + @test Complex{Float64}(gamma(Complex{BigFloat}(big"-3.4", big"-0.1"))) ≈ gamma(-3.4-0.1im) + @test Complex{Float64}(gamma(Complex{BigFloat}(big"1.0", big"6.5"))) ≈ gamma(1.0+6.5im) + # consistency with exp(loggamma) + @test gamma(Complex{BigFloat}(big"1.4", big"3.7")) ≈ exp(loggamma(Complex{BigFloat}(big"1.4", big"3.7"))) + @test gamma(Complex{BigFloat}(big"-3.4", big"-0.1")) ≈ exp(loggamma(Complex{BigFloat}(big"-3.4", big"-0.1"))) + @test gamma(Complex{BigFloat}(big"1.0", big"6.5")) ≈ exp(loggamma(Complex{BigFloat}(big"1.0", big"6.5"))) + # zero-imaginary inputs should match the real BigFloat results + @test loggamma(Complex{BigFloat}(big"3.099", big"0.0")) ≈ loggamma(big"3.099") + @test loggamma(Complex{BigFloat}(big"1.15", big"0.0")) ≈ loggamma(big"1.15") + @test gamma(Complex{BigFloat}(big"3.099", big"0.0")) ≈ gamma(big"3.099") + # branch mapping + ε = BigFloat(2)^(-1400) + xs = BigFloat.(["-0.1", "-0.7", "-1.3", "-2.6", "-3.4", "-5.2"]) + for x in xs + z_minus = Complex{BigFloat}(x, -ε) + z_plus = Complex{BigFloat}(x, ε) + + Lm = loggamma(z_minus) + Lp = loggamma(z_plus) + + zf_minus = _loggamma_oracle64_point(z_minus) + zf_plus = _loggamma_oracle64_point(z_plus) + Lf_minus = loggamma(zf_minus) + Lf_plus = loggamma(zf_plus) + + @test isapprox(Float64(real(Lm)), real(Lf_minus); rtol=0, atol=1e-14) + @test isapprox(Float64(imag(Lm)), imag(Lf_minus); rtol=0, atol=1e-14) + @test isapprox(Float64(real(Lp)), real(Lf_plus); rtol=0, atol=1e-14) + @test isapprox(Float64(imag(Lp)), imag(Lf_plus); rtol=0, atol=1e-14) + + kdiff = round(Int, (imag(Lm) - imag(Lp)) / (2*big(pi))) + @test kdiff != 0 + @test isapprox(Lm - Lp, (2*big(pi))*BigFloat(kdiff)*im; rtol=0, atol=big"1e-70") + end + # Additional _loggamma_oracle64_point nudge tests at Float64 pole + n = -3 + x64 = Float64(n) + δp64 = nextfloat(x64) - x64 + δm64 = x64 - prevfloat(x64) + ε64 = eps(Float64) + + # right of pole: rez > n but Float64(rez) == n → nextfloat + rezR = BigFloat(x64) + BigFloat(δp64)/4 + zR = Complex{BigFloat}(rezR, BigFloat(ε64)/4) + zfR = _loggamma_oracle64_point(zR) + @test real(zfR) == nextfloat(x64) + @test abs(imag(zfR)) ≤ 2ε64 + + # left of pole: rez < n but Float64(rez) == n → prevfloat + rezL = BigFloat(x64) - BigFloat(δm64)/4 + zL = Complex{BigFloat}(rezL, -BigFloat(ε64)/4) + zfL = _loggamma_oracle64_point(zL) + @test real(zfL) == prevfloat(x64) + @test abs(imag(zfL)) ≤ 2ε64 + end + end + end + + # 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 + @test loggamma(Inf - 0.0im) === Inf - 0.0im + @test loggamma(Inf + 1im) === Inf + Inf*im + @test loggamma(Inf - 1im) === Inf - Inf*im + @test loggamma(-Inf + 0.0im) === -Inf - Inf*im + @test loggamma(-Inf - 0.0im) === -Inf + Inf*im + @test loggamma(Inf*im) === -Inf + Inf*im + @test loggamma(-Inf*im) === -Inf - Inf*im + @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-67 + @test loggamma(big"1.15") ≈ big"-0.06930620867104688224241731415650307100375642207340564554412494594148673455871" rtol=1e-67 + @test logabsgamma(big"0.89")[1] ≈ big"0.0740221739580814237022658899798106584342350083445733969634566129726762260738245" rtol=1e-67 + @test loggamma(big"0.91") ≈ big"0.0589225676238323792982417511839070778835929820947704491677379048793029707373113" rtol=1e-67 + @test loggamma(big"0.01") ≈ big"4.59947987804202172251394541100874808726100141338528965241917138771477998920321" rtol=1e-67 + @test loggamma(1 + exp2(big"-20.0")) ≈ big"-5.50475006614886679092243442349111109814456565183691425527816079744208067935466e-7" rtol=1e-67 + + # 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)) + @test loggamma(complex(3//2, 1//3)) ≈ loggamma(complex(1.5, 1 / 3)) + + # negative values + @test loggamma(big(-3.0)) == big(Inf) + @test loggamma(big(-1.5)) == logabsgamma(big(-1.5))[1] + @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 isnan(loggamma(big(-Inf))) + @test logabsgamma(big(Inf))[1] == big(Inf) + @test isnan(logabsgamma(big(-Inf))[1]) + + # BigFloat signed-zero edge cases for sign of Γ(x) + @test logabsgamma(BigFloat(0.0)) == (big(Inf), 1) + @test logabsgamma(BigFloat(-0.0)) == (big(Inf), -1) + + # Complex{BigFloat} non-finite branches + @test loggamma(Complex{BigFloat}(big(Inf), BigFloat(0.0))) == Complex{BigFloat}(big(Inf), BigFloat(0.0)) + @test loggamma(Complex{BigFloat}(big(Inf), big(1.0))) == Complex{BigFloat}(big(Inf), big(Inf)) + @test loggamma(Complex{BigFloat}(big(Inf), big(-1.0))) == Complex{BigFloat}(big(Inf), big(-Inf)) + @test loggamma(Complex{BigFloat}(big(-Inf), big(1.0))) == Complex{BigFloat}(big(-Inf), big(-Inf)) + @test loggamma(Complex{BigFloat}(big(-Inf), big(-1.0))) == Complex{BigFloat}(big(-Inf), big(Inf)) + @test loggamma(Complex{BigFloat}(big(1.0), big(Inf))) == Complex{BigFloat}(big(-Inf), big(Inf)) + @test loggamma(Complex{BigFloat}(big(1.0), big(-Inf))) == Complex{BigFloat}(big(-Inf), big(-Inf)) + z_nanbf = loggamma(Complex{BigFloat}(big(Inf), big(Inf))) + @test isnan(real(z_nanbf)) && isnan(imag(z_nanbf)) + end + +@test gamma(big"0.29384") ≈ exp(loggamma(big"0.29384")) +@test gamma(big"0.29384"+big"0.12938"*im) ≈ exp(loggamma(big"0.29384"+big"0.12938"*im)) +