10889 extended domain coverage of betaIncomplete and betaIncompleteCompl #10909
+223
−24
Conversation
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
Bring std.mathspecial.betaIncomplete's handling of NaN parameter values in line with how D's operators like operator!"+" handle them. I.e., the NaN with the largest payload is returned.
Support was added to betaIncomplete and betaIncompleteCompl for the cases when a or b is +0 or infinity.
When x > 1/b and x is larger than the mean, the complement is computed instead. x -> 1-x. If prior to reversal, x < real.epsilon, 1-x = 1, so x = 1 after reversal. Since algorthm requires x < 1, the computation of complement fails. To resolve this failure, when the reversal causes x -> 1, use x -> nextDown(1-x) instead. This fix resolved a problematic case for betaIncompleteInv that is tracked by a unittest. The unittest was updated.
Combine the Stirling approximations of the individual log gamma terms when one of them is large and the difference between the arguments is significant. This reduces the impact of subtracting the two larger terms.
In betaDistPowerSeries, when s*a is close to 1, extreme loss of precision can occur when combining terms in log space. This fix resolves this issue.
The Beta CDF has a hard upper limit of 1. Due to round-off error, the result of betaDistPowerSeries can be very slightly larger than 1 when it should be 1. In this case, a value of 1 is returned. To prevent clamping from hiding a more serious bug, if the result would be greater than 1 + 2 epsilon, NaN is returned to indicated the computation failed.
Contributor
|
Thanks for your pull request and interest in making D better, @tedgin! We are looking forward to reviewing it, and you should be hearing from a maintainer soon.
Please see CONTRIBUTING.md for more information. If you have addressed all reviews or aren't sure how to proceed, don't hesitate to ping us with a simple comment. Bugzilla referencesYour PR doesn't reference any Bugzilla issue. If your PR contains non-trivial changes, please reference a Bugzilla issue or create a manual changelog. Testing this PR locallyIf you don't have a local development environment setup, you can use Digger to test this PR: dub run digger -- build "master + phobos#10909" |
thewilsonator
approved these changes
Dec 6, 2025
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
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.
This PR addresses #10889. It makes the following changes to
std.mathspecial.betaIncompleteandstd.mathspecial.betaIncompleteCompl.operator!"+"handle them. I.E., when multiple NaN arguments are provided, the one with the largest payload is returned.a,b,x) parameter space whereaandbwere varied logarithmicly from 10^-38 to 10^38, andxwas varied from 0 to 1, concentrating on values near 0 or near 1. This revealed several places in parameter space where betaIncomplete returned NaN or a value outside of its range [0, 1]. Its algorithm is modified to handle these places.a=bandx=0.5.gamma(a)*gamma(b) / gamma(a+b), are replaced with calls tobeta(a, b)to all betaIncomplete to benefit from beta`s handling of special cases.logGamma(a) + logGamma(b) - logGamma(a+b), were replaced with a newly created functionlogBetathat improves the precision of the computation by reducing the impact of subtraction when a and b have an order of magnitude size difference and b is large.betaDistPowerSeriesis modified to handle the case where naive subtraction results in an extreme loss of precision.betaDistPowerSeriesreturning1+real.epsiloninstead of1. Values that are slightly larger than1are clamped to1. If a value that is more than slightly larger than 1 would be returned, this indicatesbetaDistPowerSeriesfailed, so NaN is returned.