ibeta_inv with subnormal a or b sets the FE_UNDERFLOW exception

[WarrenWeckesser]

This is another case where I'm not sure if the behavior is expected, undefined or a bug. If I call ibeta_inv(a, b, p) with a or b a subnormal number with the default Policy (e.g. ibeta_inv(12.5, 1e-320, 0.9), on return the floating point exception flag FE_UNDERFLOW is set. Does the default policy of ignoring underflow imply that the floating point flag FE_UNDERFLOW should not be set on return?

Demonstration code

#include <iostream>
#include <iomanip>
#include <cfenv>
#include <boost/math/special_functions/beta.hpp>

using namespace std;
using namespace boost::math;

void show_fp_exception_flags()
{
    if (std::fetestexcept(FE_DIVBYZERO)) {
        cout << " FE_DIVBYZERO";
    }
    // FE_INEXACT is common and not interesting.
    // if (std::fetestexcept(FE_INEXACT)) {
    //     cout << " FE_INEXACT";
    // }
    if (std::fetestexcept(FE_INVALID)) {
        cout << " FE_INVALID";
    }
    if (std::fetestexcept(FE_OVERFLOW)) {
        cout << " FE_OVERFLOW";
    }
    if (std::fetestexcept(FE_UNDERFLOW)) {
        cout << " FE_UNDERFLOW";
    }
    cout << endl;
}

int main(int argc, char *argv[])
{
    double a = 12.5;
    double b = 1e-320;
    double p = 0.9;

    std::feclearexcept(FE_ALL_EXCEPT);

    double x = ibeta_inv(a, b, p);

    show_fp_exception_flags();

    std::cout << std::scientific << std::setw(24)
              << std::setprecision(17) << x << std::endl;

    return 0;
}

[WarrenWeckesser]

On the other hand, if I change the policy to disable promoting double to long double, the demo code above will throw an exception when ibeta_inv(12.5, 1e-320, 0.9) is called:

$ g++ -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 -I /home/warren/repos/git/forks/boost/math/include ibeta_inv_subnormal.cpp -o ibeta_inv_subnormal
$ ./ibeta_inv_subnormal 
terminate called after throwing an instance of 'boost::wrapexcept<std::overflow_error>'
  what():  Error in function boost::math::beta<double>(double,double): numeric overflow
Aborted (core dumped)

That overflow is the result of beta(a, b) eventually computing 1/1e-320, which overflows to inf. (Perhaps the handling of subnormals in beta() should be a separate issue.)

[jzmaddock]

OK, I have a fix in the works for the second case, but the first is more problematic. The question is this: do we want to store and reset FPU flags on every special function call? We can certainly do that, at some small extra cost, but note that this potentially effects every special function, as almost any expression can underflow harmlessly and the function still give the correct result. Likewise there may be (mostly) harmless overflows in all sorts of intermediate calculations... in this case I say "mostly" because some floating point types may not support infinity and almost anything might happen in that case.

What do folks think? It would be relatively trivial to write a scoped object that saves and restores FPU exception flags, maybe with a sanity check for genuine underflow. It's really just a question of whether it's worth the effort. @WarrenWeckesser do you have a situation where these flags being set is more than just a minor nuisance? I must admit having the overflow flag incorrectly set just looks plain wrong.

[WarrenWeckesser]

@WarrenWeckesser do you have a situation where these flags being set is more than just a minor nuisance?

It depends, but in this case, it is fine just to know what to expect from boost/math. I ran into these behaviors while tracking down a seg. fault in the SciPy wrapper, and I wasn't sure if boost/math made any promises about the state of the flags after a call to a boost function.

That seg. fault turned out to be a bug in the SciPy code that handled the overflow exception generated by a call such as ibeta_inv(12.5, 1e-320, 0.9). (In case you're interested: it wasn't clear from the docs that the message argument of the user-defined handler for overflow can be nullptr; dereferencing message caused the seg. fault.)

[WarrenWeckesser]

@jzmaddock, after #883, with double promotion disabled, the call ibeta_inv(12.5, 1e-320, 0.9) no longer throws an overflow error, but the call ibeta_inv(1e-320, 12.5, 0.9) (a and b swapped) does throw.

[jzmaddock]

after #883, with double promotion disabled, the call ibeta_inv(12.5, 1e-320, 0.9) no longer throws an overflow error, but the call ibeta_inv(1e-320, 12.5, 0.9) (a and b swapped) does throw.

Addressed by #884 along with some more methodical testing.

[jzmaddock]

Can this issue be closed now, or are there still things to address?