Add Wright's generalized Bessel function

[nalimilan]

Translation of SciPy's Cython implementation (scipy/scipy#11313). The only changes are to use idiomatic and efficient Julia, and to call Julia equivalents of SciPy functions. All constants are the same except for rgamma_zero and exp_inf which are slightly adjusted to what openlibm returns.

(This function is useful to implement the Tweedie distribution in Distributions.jl.)

CI failures appear to be unrelated (#513).

[nalimilan]

@devmotion Are you the key dev in this repo too? :-)

[JeffreySarnoff]

Test Summary: | Pass Broken Total Time
wrightbessel | 3842 15 3857 1.0s

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.51%. Comparing base (f534b8c) to head (1020511).
⚠️ Report is 27 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #512      +/-   ##
==========================================
+ Coverage   94.17%   94.51%   +0.34%     
==========================================
  Files          14       15       +1     
  Lines        2969     3174     +205     
==========================================
+ Hits         2796     3000     +204     
- Misses        173      174       +1     

Flag	Coverage Δ
unittests	94.51% <100.00%> (+0.34%)	⬆️

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

☔ View full report in Codecov by Harness.
📢 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.

[devmotion]

Thanks for the PR. I checked the branch against a BigFloat evaluation of the defining series and compared it with current upstream (scipy/xsf, include/xsf/wright_bessel.h). The reference values in test/data/wrightbessel.txt are fine, I verified all 3516 rows and the worst relative error is 1.8e-15. Some comments below, the first five I consider blocking.

1. Two InexactErrors, both reachable. floor(Int, ...) throws where the C++ cast is silently UB:

wrightbessel(1e-20, 175.0, 0.5)   # InexactError: Int64(3.472411158866379e20), _wb_series:62
wrightbessel(5.0,   1.0,  1e300)  # InexactError: Int64(2.615320972023645e49), _wb_large_a:80

The first is reached for b in (170, 178.47), which passes none of the three _wb_small_a conditions and falls through to _wb_series. In _wb_series a clamp suffices:

kf = (rgamma_zero - b) / a
kf < nstop && (nstop = floor(Int, kf))

In _wb_large_a I would drop the Int conversion entirely since loggamma accepts Float64:

k_max   = floor((a^(-a) * x)^(1 / (1 + a)))
n_start = max(0.0, k_max - (n ÷ 2))
for i in 0:(n - 1)
    k = n_start + i
    res += exp(k*lnx - loggamma(k+1) - loggamma(a*k+b) - max_exponent)
end

This is bit-identical on (5,1,1e4), (10,1,1e12), (5,100,1e5), (100,0,1e20), does not allocate, and returns 1.5691925832141792e50 for x = 1e300.

2. Negative return values. Φ is positive but _wb_integral can produce a negative res1 + res2, and then exp(exp_term) overflows:

wrightbessel(0.1, 50.0, 3000.0)   # -Inf

On a grid with a in [0.1, 4], b in [10, 150], x in [300, 1e5] (936 points) I get 63 negative results and 173 NaNs, most of them outside the documented 0.5 <= a <= 1.8 && b >= 100 && x >= 1e5 hole. Error in log(Φ) for a = 0.3:

b \ x	1e3	3e3	1e4	3e4
20	4e-13	2e-6	3e-10	-2e-12
50	-3e-14	15	380	3e-6
100	-2e-7	NaN	450	2300

I assume this comes from upstream (your res > 0 ? ... : NaN is equivalent to log of a negative number in C), but it would be new behaviour here. At the very least the non-log branch should return NaN as well so that the two branches agree, and the NaN domains have to be documented. Currently the docstring does not mention them at all, so users just see a NaN. Extending the hard-domain check to cover the region above would be preferable to returning a wrong finite value.

3. Why Float64 only? These all throw a MethodError:

wrightbessel(1, 2, 3)
wrightbessel(1.0, 2.0, 3)
wrightbessel(1f0, 2f0, 3f0)

Elsewhere we promote in the public method and dispatch to an internal one, see zeta in src/gamma.jl:233, logabsbeta in src/gamma.jl:925, and the Float32/Float16 methods in src/erf.jl:40-41:

wrightbessel(a::Real, b::Real, x::Real)    = _wrightbessel(map(float, promote(a, b, x))..., false)
logwrightbessel(a::Real, b::Real, x::Real) = _wrightbessel(map(float, promote(a, b, x))..., true)

_wrightbessel(a::T, b::T, x::T, logret::Bool) where {T<:Union{Float16,Float32}} =
    T(_wrightbessel(Float64(a), Float64(b), Float64(x), logret))

map(float, promote(a, b, x)) rather than promote(float(a), float(b), float(x)): float(1) === 1.0, so the latter widens Float32 arguments to Float64. Dispatching to _wrightbessel and not back to wrightbessel also matters, map(float, promote(...)) is a fixed point for BigFloat and a self-recursive fallback would not terminate. The docstring signature needs updating too.

4. NaN arguments should return NaN, not throw. This is what upstream does and what we do elsewhere (gamma(NaN), besselj(0, NaN)). Throwing also makes the function inconvenient for the Distributions.jl use case that motivates the PR.

5. DomainError instead of ArgumentError for negative arguments, cf. src/bessel.jl:438, src/gamma.jl:103, src/ellip.jl:51. We use ArgumentError only for things like unsupported orders. The message says "arguments must be positive" but the check is < 0.

Accuracy

rgamma(y) = exp(-loggamma(y)) is less accurate than the cephes::rgamma used upstream, since the absolute error of the exponent becomes the relative error of the result:

b	exp(-loggamma(b))	inv(gamma(b))
50	9.2e-15	1.4e-16
100	5.6e-14	6.1e-17
171	1.2e-13	2.0e-19

This affects wrightbessel(a, b, 0.0) and every term of _wb_series. Something like

_rgamma(y::Real) = y <= 171.6 ? inv(gamma(y)) : exp(-loggamma(y))

is more accurate and cheaper for the common case. The exp(-loggamma) branch is still needed for subnormal results, where inv(gamma) returns zero.

The claim "accuracy is generally higher than 1e-11" does not hold outside the cases you list, e.g. (1.0, 100, 1e4) gives 1.0e-3, (0.5, 80, 3e3) gives 1.3e-4, (1.0, 50, 1e4) gives 3.5e-9 and (0.001, 50, 8.0) gives 5.3e-12. Can you make the statement domain-dependent, or restrict the domains?

sin(u + π*b) in _Kmod loses about b*eps of phase, sin(u)*cospi(b) + cos(u)*sinpi(b) is exact in b (3.8e-14 vs 4.5e-18 at b = 170). I measured this end to end and it only improves points that are already accurate, so it is cosmetic and does not address point 2.

Type stability

_wb_asymptotic allocates 1152 bytes per call, _wb_small_a 96 bytes for order >= 3. Two causes:

• ntuple(i -> a^(i-1), 15) and 17 in lines 226-228. ntuple is only unrolled up to n = 10, above that it builds a Vector and splats it, and infers as Tuple{Vararg{Float64}}. Use Val(15) and Val(17).
• (A1, ..., A6)[1:(order+1)] in lines 180 and 207. Indexing a tuple with a non-constant range loses the length. Initialise the unused coefficients to 0.0 instead of NaN and always evaluate the full polynomial, they are multiplied by a^k so the result is unchanged. I checked this is bit-identical for orders 2 to 5.

res += (-1)^(k+1) * C[k] * Zp in lines 551-554 can be evalpoly(inv(Z), (C1, -C2, C3, ...)), which agrees to 1 ulp and is unrolled.

Style

• rgamma, rgamma_zero, exp_inf, x_laguerre, w_laguerre, x_legendre, w_legendre are added to the SpecialFunctions namespace without a prefix. x_legendre in particular is too generic. Can you prefix them or put the file in a submodule? rgamma is also annotated ::Real although it is only valid for non-negative arguments.
• In the else branch of _wb_small_a, order >= 3 holds and order = min(order, 5) has already run, so if order >= 2 is always true and the nesting is an artefact of the translation. Line 186 is indented by one space too many.
• eps in line 723 shadows Base.eps.
• log(π) in line 769 and sqrt(2.0 * π * ...) in line 232 can use logπ and sqrt2π, we already depend on IrrationalConstants.
• eachindex(x_laguerre, w_laguerre, x_legendre) omits w_legendre.
• nstart of _wb_series is 0 at both call sites.

Docs

• logwrightbessel is added to functions_list.md but not to functions_overview.md. The overview does list loggamma, logbeta and logabsgamma.
• For a >= rgamma_zero all terms with k >= 1 underflow, so Φ = 1/Γ(b) exactly. Returning rgamma(b) and -loggamma(b) would be cheap and better than NaN.
• functions_overview.md adds a blank line after the new row, and line 787 has trailing whitespace.
• The PR description says rgamma_zero and exp_inf were adjusted, but both are bit-identical to upstream (178.47241115886637 and 0x40862e42fefa39f0).

Tests

• @test wrightbessel(a, b, 0.0) == SpecialFunctions.rgamma(b) in line 27 is trivially true, rgamma(b) is the implementation for x == 0. Comparing with inv(gamma(b)) would also show the accuracy issue above.
• I would not use @test_broken in line 103. These are documented accuracy limits, not known bugs, and @test_broken fails as Unexpectedly Pass if the accuracy improves, which is a realistic flakiness risk across our CI matrix. Just test the reduced tolerance.
• Can the reference values be computed on the fly with BigFloat instead of shipping test/data/wrightbessel.txt? The grid needs at most 500 series terms, so the plain defining series converges everywhere in it. This generator needs nothing beyond Base and our own loggamma:

  function wrightbessel_ref(a, b, x; prec = 128)
      setprecision(BigFloat, prec) do
          A, B, X = BigFloat(a), BigFloat(b), BigFloat(x)
          iszero(X) && return exp(-loggamma(B))
          s = zero(BigFloat); k = 0; small = 0
          while true
              t = exp(k*log(X) - loggamma(BigFloat(k)+1) - loggamma(A*k+B))
              s += t
              if !iszero(s) && t < eps(BigFloat)*s
                  small += 1
                  small > 5 && break
              end
              k += 1
          end
          return s
      end
  end

  It reproduces every row of the file to 1.8e-15 and takes 7.5 s for all 3516 rows at prec=128 (5.7 s at 96 bits), which seems acceptable for the test suite. That removes 344 KB from the tarball, makes the reference self-documenting, drops the DelimitedFiles dependency, and lets the grid be extended without regenerating a file. wrightbessel_nemo could then go away as well.
• The coverage gaps correspond to the problems above: no x > 1e20, no b in (170, 178.47) with small a, no a = 0.3 and no x = 1e4 in the grid, no Int or Float32 arguments, no @inferred and no allocation tests. A positivity test would have caught the 63 negative results:

  @test all(wrightbessel(a, b, x) >= 0 for a in ..., b in ..., x in ...)

A ChainRulesCore rule would be nice in a follow-up, ∂/∂x Φ(a, b, x) = Φ(a, b-1, x) is already used as a test identity here.