Merge pull request #522 from TSGut/master

Implement loggamma(z::Complex{BigFloat})

Author: JeffreySarnoff

src/gamma.jl

@@ -759,6 +759,100 @@ function loggamma_asymptotic(z::Complex{Float64})
                          6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02)
 end
 
+# Return scaled Sterling coefficients (B_{2k}/(2k*(2k-1)) * zr^(1-2k) for k=0,1,...,n)
+function _scaled_sterling_coeffs(n::Integer, zr::Complex{BigFloat})
+    mmax = 2n
+    A = Vector{Rational{BigInt}}(undef, mmax+1)
+    E = Vector{Complex{BigFloat}}(undef, n+1)
+
+    # Compute Bernoulli numbers using Akiyama-Tanigawa algorithm, and store the scaled Stirling coefficients in E
+    @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
+            # store term B_{2k}/(2k*(2k-1)) * zr^(1-2k)
+            E[k + 1] = A[1] / (2k * (2k - 1)) * zr^(1-2k)
+        end
+    end
+
+    return E
+end
+
+function _loggamma_complex_bigfloat(z::Complex{BigFloat})
+    # We use branch correction (offset by multiples of 2πi) 
+    # using Float64 logic instead of complicated manual high precision branch tracking
+    bigpi = big(pi)
+
+    # Reflection formula
+    rez = real(z)
+    if rez < 0.5
+        val = log(bigpi) - log(sinpi(z)) - _loggamma_complex_bigfloat(1 - z)
+        return _loggamma_branchcorrect(val, z)
+    end
+
+    # Upward recurrence: shift z to the Stirling region
+    p = precision(BigFloat)
+    r = max(0, Int(ceil(p - abs(z))))
+    zr = z + r
+
+    # Stirling
+    N = max(10, p÷15)
+    B = _scaled_sterling_coeffs(N, zr)
+    lg = sum(B[2:end]) + (zr - big"0.5")*log(zr) - zr + log(sqrt(2*bigpi))
+
+    # Undo the upward shift via recurrence in single product form
+    prodarg = prod([(z + (i-1)) for i in 1:r])
+    s = log(prodarg)
+
+    # Apply branch correction
+    return _loggamma_branchcorrect(lg - s, z)
+end
+
+function loggamma(z::Complex{BigFloat})
+    imz = imag(z)
+    rez = real(z)
+    if iszero(imz)
+        # if z is purely real, then we can just call the real loggamma
+        return Complex(loggamma(rez))
+    end
+    # Guardrail precision by 16 bits for complex BigFloat inputs,
+    # compute using the internal implementation, then round the result back.
+    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
+
+# branch correct loggamma by offsetting by multiples of 2πi to match the Float64 version
+function _loggamma_branchcorrect(val_big::Complex{BigFloat}, z::Complex{BigFloat})
+    zf = _loggamma_oracle64_point(z)
+    val_f = loggamma(zf)
+    imz = imag(val_big)
+    k = round(Int, (Float64(imz) - imag(val_f)) / (2*pi))
+    return Complex{BigFloat}(real(val_big), imz - 2*big(pi)*k)
+end
+
+# map a BigFloat complex point to the appropriate Float64 complex point, so that we can use the Float64 loggamma to determine the correct branch cut for the BigFloat loggamma
+function _loggamma_oracle64_point(z::Complex{BigFloat})
+    rez = real(z)
+    xr = Float64(rez)
+    xi = Float64(imag(z))
+    n = round(Int, xr)
+    if n ≤ 0 && isapprox(xr, Float64(n); atol=eps(Float64)) && abs(xi) ≤ 2eps(Float64)
+        xr = rez > n ? nextfloat(xr) : prevfloat(xr)
+    end
+    return Complex{Float64}(xr, xi)
+end
+
+
 @doc raw"""
     beta(x, y)
 

test/gamma.jl

@@ -153,6 +153,76 @@
         @test logabsgamma(big(Inf))[1] == big(Inf)
     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
+
     @testset "Other float types" begin
         struct NotAFloat <: AbstractFloat
         end
@@ -165,11 +235,6 @@
         @test_throws MethodError gamma(NotAFloat())
         @test_throws MethodError logabsgamma(NotAFloat())
         @test_throws MethodError loggamma(NotAFloat())
-
-        # https://github.com/JuliaMath/SpecialFunctions.jl/issues/339
-        # https://github.com/JuliaMath/SpecialFunctions.jl/issues/233
-        @test_throws MethodError gamma(complex(big(1.0)))
-        @test_throws MethodError loggamma(complex(big(1.0)))
     end
 end
 

test/runtests.jl

@@ -11,7 +11,7 @@ using Aqua: Aqua
 using ExplicitImports: ExplicitImports
 using JET: JET
 
-using SpecialFunctions: AmosException, f64
+using SpecialFunctions: AmosException, f64, _loggamma_oracle64_point
 
 # useful test functions for relative error, which differ from isapprox
 # relerr separately looks at the real and imaginary parts if one of the arguments is complex