U(a,b,z) returns results different from Mathematica/MATLAB when z is negative #63

Activity

1. MikaelSlevinsky commented on May 3, 2023

   

   MikaelSlevinsky

   on May 3, 2023

   Collaborator

   More actions

   Thanks for filing the issue! Evaluating U near the branch cut (-∞,0] requires a different approach, such as Kummer's transformation #61.
2. AshtonSBradley

   added a commit that references this issue on Apr 12, 2026

   resolve JuliaMath#61 , adding tests from JuliaMath#63

   691da1b
3. AshtonSBradley commented on Apr 14, 2026

   

   AshtonSBradley

   on Apr 14, 2026

   More actions

   @MikaelSlevinsky I hope #92 resolves this in the way you had in mind.
4. hm271 commented on Aug 4, 2026

   

   hm271

   on Aug 4, 2026

   More actions

   I can confirm this is still present in v0.3.30 (Julia 1.12.6), and I have some
   additional characterization that may help narrow it down.

   1. The error is independent of working precision.

   Evaluating at increasing BigFloat precision gives bit-for-bit identical
   relative errors, so this is not a rounding/conditioning issue but an
   algorithmic one:

   precision	rel. err, U(0.5, 0.5, -2.5)	rel. err, U(-0.5, -0.5, -2.5)
   128	1.9788244846232889e-1	6.551202723008086e-2
   256	1.9788244846232889e-1	6.551202723008086e-2
   512	1.9788244846232889e-1	6.551202723008086e-2

   Reference values are from Arb via Nemo.hypergeometric_u, which carries
   rigorous error bounds.

   2. The failure region is a narrow wedge around arg z = π, not just the axis.

   With z = -2.5 + iδ at 256-bit working precision (a correct result here is
   ~1e-73):

   δ	arg(z)/π	rel. err
   0 … 1e-6	1.000000	1.98e-1
   1e-3	0.999873	2.15e-5
   1e-2	0.998727	2.51e-43
   1e-1	0.987274	1.21e-73 (correct)

   3. Everywhere else in the z-plane is fine at full precision.

   Scanning z = r*exp(im*θ) over θ ∈ [0, 2π) in steps of π/12, for r = 2.5, 10,
   40, and for three parameter sets — (a,b) = (0.5, 0.5), (1.5, 2.0), and a
   complex pair (5.18+0.85im, 6.0+0.1im) — every point except θ = π agrees with
   Arb to ~1e-74…1e-77 at 256 bits. So the algorithm is sound in general; the
   problem is specific to the negative real direction.

   4. The package already detects the failure, but returns a value anyway.

   Every bad point emits

   ┌ Warning: Rational approximation to ₂F₀ reached the maximum type of (1048576, 1048576).
   └ @ HypergeometricFunctions .../src/weniger.jl:221
   

   and the warning fires on exactly the bad region (plus a slightly wider margin
   where accuracy is degraded but still usable). Since U returns a
   plausible-looking number regardless, this is a silent-wrong-answer path for any
   caller not watching stderr — which is how I ran into it, in a black-hole
   perturbation theory code where U is called inside a series sum.

   Suggestion

   I understand #61 tracks the missing convergent path via Kummer's
   transformation, which is the real fix. Ahead of that, it would already help a
   lot if U either threw or returned NaN when the Weniger rational
   approximation hits its maximum type, rather than returning the unconverged
   value — that turns a silent numerical error into a visible one.

   I'm not familiar enough with the internals to propose a patch, but I'm happy to
   share the full scan script or test any additional parameter ranges that would
   be useful.