boostorg/math

Unexpected exception thrown by `hypergeometric_1F1`

WarrenWeckesser opened this issue · 2 comments

Calling (for example) hypergeometric_1F1(13.0, 1.5, 61.0) throws the exception:

libc++abi: terminating due to uncaught exception of type std::domain_error: Error in function boost::math::tgamma<long double>(long double): Evaluation of tgamma at a negative integer -1.

This is on macos 13.6 with Apple clang version 14.0.3 (clang-1403.0.22.14.1).

I'm using the develop branch, commit 1722ef9.

My test code is hyp1f1.cpp.

The exception is not thrown when I change b to something that is not exactly 1.5. For example,

% ./hyp1f1 13 1.499999 61                                                                   
1F1 =   1.35509150477570084e+39
pFq =   1.35509150509447403e+39
% ./hyp1f1 13 1.500001 61
1F1 =   1.35508003684581251e+39
pFq =   1.35508003682636603e+39

Note that the result computed with hypergeometric_pFq is much more accurate in both cases. The expected values for these calculations can be computed with mpmath in Python:

In [28]: from mpmath import mp

In [29]: mp.dps = 250

In [30]: float(mp.hyp1f1(13.0, 1.499999, 61))
Out[30]: 1.3550915050944761e+39

In [31]: float(mp.hyp1f1(13.0, 1.500001, 61))
Out[31]: 1.3550800368263657e+39

The exception is thrown in the call of tgamma here:

result *= boost::math::tgamma(b_minus_a - 0.5f, pol) * pow(z / 4, -b_minus_a + T(0.5f));

The exception is thrown when b - a is a half-integer and the parameters are sufficiently large (I haven't tried to figure out the exact conditions). For example,

$ ./hyp1f1 14.25 2.75 54
terminate called after throwing an instance of 'boost::wrapexcept<std::domain_error>'
  what():  Error in function boost::math::tgamma<long double>(long double): Evaluation of tgamma at a negative integer -1.
Aborted (core dumped)
$ ./hyp1f1 14.125 2.625 100
terminate called after throwing an instance of 'boost::wrapexcept<std::domain_error>'
  what():  Error in function boost::math::tgamma<long double>(long double): Evaluation of tgamma at a negative integer -1.
Aborted (core dumped)

Confirmed, looking for the best fix now...