Fix for scipy issue 18302

[mborland]

See: scipy/scipy#18302

Adds bounds checking so we don't return a pdf larger than beta, and catches an overflow exception when x is near std::numeric_limits::min.

[WarrenWeckesser]

Thanks for the fast response, @mborland!

It looks like the code in test/scipy_issue_18302.cpp uses the default policy--is that correct? I wrote a short C++ code to check the case x = 1e-310, and I didn't get the overflow error with the default policy. I had to use promote_double<false>.

[mborland]

Thanks for the fast response, @mborland!

It looks like the code in test/scipy_issue_18302.cpp uses the default policy--is that correct? I wrote a short C++ code to check the case x = 1e-310, and I didn't get the overflow error with the default policy. I had to use promote_double<false>.

I did get overflow with the default policy. Now that you bring it up I ran it on an ARM Mac which means long double = double. I think in this case it's worth effectively ignoring the overflow exception because we know what the result is supposed to be when that happens.

[WarrenWeckesser]

Ah, I see. I tested on a Linux machine, where long double is extended precision. I was just worried that the test wasn't actually testing the fix because it would pass even without the fix. But that is not the case on platforms where long double is double.

[mborland]

Ah, I see. I tested on a Linux machine, where long double is extended precision. I was just worried that the test wasn't actually testing the fix because it would pass even without the fix. But that is not the case on platforms where long double is double.

We recently added ARM and s390x to the CI system so we have 64, 80, and 128 bit long doubles are all tested now. Its helped find overflows and precision issues.

[WarrenWeckesser]

Nice. I think it might be slightly better to use promote_double<false> so the test in fact fails on all platforms without the fix, and therefore the test is verifying that the change actually fixes the problem on all platforms, but obviously that's your call.

[jzmaddock]

@mborland I think we should fix this updatream in ibeta_derivative here:

math/include/boost/math/special_functions/beta.hpp

Line 1455 in 109a814

T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function);

It's the 1/y that fails and we need to decide whether to return the derivative at zero (which may well be infinite) or at tools::min_value<T>() which may well not be. I'm edging towards the latter as the easier fix.

[mborland]

@jzmaddock This is green and I think hits your comments from both threads.

[mborland]

@jzmaddock Any other changes for this one, or good to merge?

[jzmaddock]

Looks good to me!