Fix DivisionByZeroError in the hypergeometric family for alternating series - #492
Open
gaoflow wants to merge 1 commit into
Open
Conversation
The shared generalizedHypergeometric() convergence loop tested abs(product / sum) > tol, dividing by the running partial sum. For alternating series (e.g. z < 0) that sum passes through exactly 0 on an early partial sum, so the test raised DivisionByZeroError on finite, well-defined values. This reached the public confluentHypergeometric (1F1) and hypergeometric (2F1) entry points, e.g. 1F1(1,1,-1), which is the identity e^z and should return e^-1 = 0.3679. Test convergence by comparing the term magnitude to tol * |sum| instead: multiplying can never divide by zero, and for a non-zero sum it is algebraically identical to the old test, so already-correct results are unchanged. When the partial sum is momentarily 0 the loop simply keeps iterating. The existing MAX_ITERATIONS guard is retained. The sibling lowerIncompleteGammaSeries() in the same file already guards its identical divide-by-accumulator (abs(sum) > 0 && ...); this back-ports that safety to the hypergeometric loop. Adds regression tests across 1F1 and 2F1 for inputs whose partial sum crosses zero, plus the 1F1(a;a;z) = e^z axiom over the sign of z. Expected values from mpmath (dps 45 and 50 agree).
Author
|
The red jobs aren't from this change: all 53 failures are |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The shared
generalizedHypergeometric()convergence loop testswhich divides by the running partial sum
$sum. For any alternating series the partial sum passes through exactly0on an early iteration, so this raisesDivisionByZeroErroron finite, well-defined inputs. It reaches the publicconfluentHypergeometric(₁F₁) andhypergeometric(₂F₁) entry points:₁F₁(a;a;z) ≡ eᶻis a textbook identity, so the first case unambiguously should returne⁻¹.A grid sweep of ₁F₁ (
a ∈ {-0.5, 0.5, 1, 1.5, 2, 3},b ∈ {0.5, 1, 1.5, 2, 3},z ∈ -3..5) and ₂F₁ (|z| < 1) found 15 crashing inputs, all the same root cause, all on finite values.Fix
Compare the term magnitude to
$tol * \abs($sum)(multiply) instead of\abs($product / $sum)(divide). Multiplying can never divide by zero, and for a non-zero$sumit is algebraically identical to the old test, so every already-correct result is unchanged. When the partial sum is momentarily zero the loop simply keeps iterating. TheMAX_ITERATIONSguard is retained.The sibling
lowerIncompleteGammaSeries()in the same file already guards the identical divide-by-accumulator (\abs($sum) > 0 && \abs($term / $sum) < $tol); this brings the hypergeometric loop in line with it.No regression
I compared the fixed output against the unpatched output across the full grid: all 500+ currently-working inputs are bit-for-bit identical (as expected, since the convergence test is only algebraically rearranged), and the 15 formerly-crashing inputs now match
mpmathhyp1f1/hyp2f1(dps 30 and 35 in agreement) to the series' ~1e-8 tolerance.Tests
Added data-provider regression tests across ₁F₁ and ₂F₁ for the zero-crossing inputs, plus the
₁F₁(a;a;z) = eᶻaxiom over the sign ofz. Expected values are frommpmath(dps 45 and 50 agree). Full suite is green (phpunit, phpcs, phpstan).