Repository navigation
Conversation
…uchy/logcdf` `ln( 0.5 + atan(z)/pi )` cancels left of the median and returns -Infinity from about `x = x0 - 1e16*gamma`; `logcdf( -1e30, 0, 1 )` is about -70.22. Use the equivalent `atan2( gamma, x0-x )/pi` there, and `log1p( -atan2( gamma, x-x0 )/pi )` right of the median, where the result is `ln` of a number close to 1. The fixtures were generated with Distributions.jl's Float64 logcdf, which uses the same formula and is up to 782 ULP off; they are now evaluated in 2048-bit BigFloat, plus a new far-left-tail fixture. All four fixtures are within 2 ULP; develop fails 1211 of the assertions.
|
Hello! Thank you for your contribution to stdlib. We noticed that the contributing guidelines acknowledgment is missing from your pull request. Here's what you need to do:
This acknowledgment confirms that you've read the guidelines, which include:
We can't review or accept contributions without this acknowledgment. Thank you for your understanding and cooperation. We look forward to reviewing your contribution! |
Coverage Report
The above coverage report was generated for the changes in this PR. |
Near the median the original `ln( 0.5 + atan(z)/pi )` doesn't cancel and is the more accurate form; the tail identities only pay off once the result is far from `ln(0.5)`. Gate them at `z < -1` and `z > 1`: on `large_gamma` this PR was further from exact than develop on 325 points, now 2.
|
Hello! Thank you for your contribution to stdlib. We noticed that the contributing guidelines acknowledgment is missing from your pull request. Here's what you need to do:
This acknowledgment confirms that you've read the guidelines, which include:
We can't review or accept contributions without this acknowledgment. Thank you for your understanding and cooperation. We look forward to reviewing your contribution! |
Description
This pull request:
stats/base/dists/cauchy/logcdfwithout cancellation (JS, factory and C).ln( 0.5 + atan(z)/pi )cancels left of the median: the sum loses everything onceatan(z)/piis within an ULP of-0.5, so the function returns-Infinityfrom aboutx = x0 - 1e16*gammaon, and is off by hundreds of ULP well before that.z = -1:0.5 + atan(z)/piequalsatan2( gamma, x0-x )/piexactly, with no subtraction.z = 1: the result islnof a number approaching1, so it useslog1p( -atan2( gamma, x-x0 )/pi ).atan2is close topi/2and picks up the rounding of1/pi): onlarge_gammathe first version was further from exact thandevelopon 325 points; gated, it's 2.x = x0still givesln(0.5);±Infinityare unchanged.The fixtures had to change. They were generated with Distributions.jl's
logcdf, which evaluates the same formula inFloat64. Againstmpmath, on the 1,046 points where this PR and the fixtures disagree, the fixtures (anddevelop) are up to 782 ULP off, and this PR is within 2.developonly passed because it reproduces the reference's own error.gennow evaluates the textbook formula in 2048-bitBigFloat, as in #15772 and #15774. All three existing fixtures were regenerated with that, plus a newleft_tail.jsonwithxfrom1e6to1e296scale units belowx0.With the new fixtures, every point is within 2 ULP (the existing tolerances of 6–12 ULP are unchanged; the new left-tail blocks use 2).
test.logcdf.js/test.factory.js/test.native.jspass 4,014 / 4,018 / 4,014. Ondevelop, 1,211 assertions intest.logcdf.jsfail.Related Issues
None.
Questions