Implement loggamma(z::Complex{BigFloat})

[TSGut]

I need this for my MeijerG function package which should hopefully launch in the next few days.

Luckily it should also fix stuff other people (@MikaelSlevinsky, @hm271) have had issues with.

Fixes #440

[TSGut]

Just noticed there is a test that now has to be removed which wants this to throw an error.. 😄

Also this addresses #339 too (though that was about error type maybe?)

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.28%. Comparing base (4fead0f) to head (521bda1).
⚠️ Report is 14 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #522      +/-   ##
==========================================
+ Coverage   94.10%   94.28%   +0.17%     
==========================================
  Files          14       14              
  Lines        2971     3026      +55     
==========================================
+ Hits         2796     2853      +57     
+ Misses        175      173       -2     

Flag	Coverage Δ
unittests	94.28% <100.00%> (+0.17%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[JeffreySarnoff]

This implementation is many times faster and at least a digit more accurate:

"""
    lngamma(x::BigFloat, r::Base.MPFR.MPFRRoundingMode = Base.Rounding.rounding_raw(BigFloat))

Compute the logarithm of the Gamma function of `x` at the current `BigFloat` precision.
"""
function lngamma(x::BigFloat, r::Base.MPFR.MPFRRoundingMode = Base.Rounding.rounding_raw(BigFloat))
    # Allocate destination with the current task-local precision
    z = BigFloat()
    
    # mpfr_lngamma C signature: int mpfr_lngamma(mpfr_t rop, mpfr_t op, mpfr_rnd_t rnd)
    ccall((:mpfr_lngamma, Base.MPFR.libmpfr), Int32,
          (Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode),
          z, x, r)
          
    return z
end

"""
    lngamma!(z::BigFloat, x::BigFloat, r::Base.MPFR.MPFRRoundingMode = Base.Rounding.rounding_raw(BigFloat))

In-place version of `lngamma`. Computes the logarithm of the Gamma function of `x` 
and stores the result directly in `z`, avoiding memory allocations.
"""
function lngamma!(z::BigFloat, x::BigFloat, r::Base.MPFR.MPFRRoundingMode = Base.Rounding.rounding_raw(BigFloat))
    ccall((:mpfr_lngamma, Base.MPFR.libmpfr), Int32,
          (Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode),
          z, x, r)
          
    return z
end

[TSGut]

Perhaps I am missing something but that doesn't do anything about Complex{BigFloat}? BigFloat itself already works in SpecialFunctions.jl. I would be quite happy if there is indeed a faster version.

[JeffreySarnoff]

oops wrong snippet

[TSGut]

In fairness there are levers one can pull even on my implementation. Here are some performance patches which bring the methods 'closer':

julia> using BenchmarkTools

julia> @btime loggamma(big"1.4");
  16.000 μs (47 allocations: 2.33 KiB)

julia> @btime loggamma(big"1.4", big"3.7");
  536.375 μs (12956 allocations: 1.18 MiB)

If we really care we could precompute Bernoulli numbers at some ridiculous precision and use a lookup but that's very clunky design.

julia> @btime _bernoulli_upto(precision(BigFloat)÷15);
  51.208 μs (3017 allocations: 83.93 KiB)

But like I said, quite happy if there is a version out there already we can just call (I am aware of mpmath and Arb ofc). It's my understanding that MPFR doesn't have one.

[JeffreySarnoff]

Your work is good. I am trying to find the correct snippet, so we can both look at it.

[JeffreySarnoff]

Ok, the speeds are about the same.

The accuracy is greater in the alternative version (see the [gist} (https://gist.github.com/JeffreySarnoff/a25cee5e562f92182cf5193b54065654))

At the end is the test file I used. On many runs it showed ~14 more accurate bits (relative to a higher precision version of the same implementation).

[Admittedly, this is not definitive -- a fuller comparison against a shared higher resolution value constrains the additional bit accuracy within 8-16, at 1 standard deviation from the mean.]

Note: The downcasting to Float64 can give an incorrect result with edge cases. See lines 275-276, and the function at 302 for one way to handle this.

I have not looked to refactor the gist to best meld with SpecialFunctions.jl .. there may well be function duplication, and where there is, the original version is much better tested.

[JeffreySarnoff]

Please borrow whatever you find useful from the gist, as is or adapted. When tests pass , I will merge your PR.

[TSGut]

Depending on the input that fails tests because its branch tracking is wrong.

Unfortunately branch tracking this approximation can get complicated as one can see from looking at what Arb and mpmath do. Because of prevalence, LLMs are keen to use localised log branch tracking but that doesn't work here.

Luckily we can still make it robust near branch cuts without that, I added a patched oracle as well as explicit near branch testing so I believe this is robust now.

My method is also actually faster than that one after some slight tweaks:

using BenchmarkTools

julia> loggamma_timer = @btime time_loggamma($rands);
  2.693 s (80703433 allocations: 2.79 GiB)

julia> alt_log_gamma_timer = @btime time_log_gamma($rands);
  2.741 s (80703433 allocations: 2.79 GiB)

Though as you say since it uses the same mathematical idea they might as well be the same timings.

I did take some inspiration from this, though, in that you seem okay with guardrails for precision which is something I personally like but other people in the past have told me they dislike - it's a design choice really. I added a 16 bit buffer precision - in principle one could increase this but it really isn't needed in practice. The accuracy is good as is.

The remaining CI failure is because of JET shenanigans in the Julia prerelease and unrelated to this PR.

Happy for another code review round or merge if that's desired. :)

[JeffreySarnoff]

:)

[MikaelSlevinsky]

Cool!

[JeffreySarnoff]

OK Need Help -- It turns out that issue with JET is now preventing the entire repo from passing the tests -- so registering the revised SpecialFuncions.jl cannot occur. Do you know how to remedy this?

[giordano]

so registering the revised SpecialFuncions.jl cannot occur.

How do you mean?

[JeffreySarnoff]

It appears I took an intermediate red x too seriously. It did merge.
We need a second ok for the new release to be accepted. Would you?

JuliaRegistries/General#151501

[giordano]

It's all good already