powm1(0.0, -2.0) doesn't throw an exception. Should it?

[WarrenWeckesser]

Here's a test program that I compiled with the develop branch:

#include <iostream>
#include <boost/math/special_functions/powm1.hpp>

using namespace std;
using boost::math::powm1;


int main()
{
    try {
        double p = powm1(0.0, -2.0);
        cout << "p = " << p << endl;
    } catch (const exception& e) {
        cerr << "*** Caught: " << e.what() << endl;
        cerr << "*** Type: " << typeid(e).name() << endl;
    }
}

I expected powm1(0.0, -2.0) to throw an exception, but it does not; the program outputs p = inf. Should it throw an exception?

[mborland]

I think the behavior is correct. If you compare to the error handling of std::pow returning infinity is expected.

[WarrenWeckesser]

OK, I'm still figuring out the conventions followed by the functions in Boost/Math/Special Functions. Based on the policy docs, I thought the default behavior would be to throw a std::domain_error.

I see that log1p(-1.0) and atanh(±1.0) also do not throw, so I guess this is not unusual.

On the other hand, erf_inv(±1.0) and cyl_neumann(0, 0.0) each throw an overflow error. (Shouldn't those exceptions be domain errors, as the closest thing to a pole error?) So in general, is the rule "It depends, you just have to try it"?

[NAThompson]

@WarrenWeckesser : Lol, I like that rule. Personally I find it preferable to use all the semantic space that IEEE floats gives us and like the inf, but I think that you could definitely argue for a domain_error.

[jzmaddock]

It's a bug, but I'm not sure in which direction: wolframalpha reports 0^-2 as complex-infinity which would be a domain error, but if we treat it as 1/0^2 then it should be a straight overflow error.