cyl_bessel_j returns incorrect sign when using --fast-math on GCC or -std=c++11 on Intel Compiler

[Dustin Lanham <lanhamd@…>]

The cyl_bessel_j function returns the incorrect sign for certain inputs when using GCC's --fast-math flag, or the Intel compiler's -std=c++11.

This bug has been shown to affect GCC 4.4.7 as delivered in RHEL 6.3 and the Intel compiler version 13.1.3 20130607 on the same OS.

The following minimal case program can replicate the problem:

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

int main(int argc, char** argv)
{
  double besselTerm = boost::math::cyl_bessel_j(1.5, 7.90559);
  std::cout << besselTerm << std::endl;
}

When this program is compiled under GCC with the --fast-math flag, or ICPC with the -std=c++11 flag, it gives a result of -0.0504874. Checking this answer against a reference MATLAB implementation, or GCC 4.1.2 shows that the correct answer should be 0.0504874. When compiled on GCC 4.4.7 *without* the --fast-math flag, or ICPC without the -std=c++11 flag, the correct answer is returned.

This behavior has been found for all inputs of P from 1.5 to 30.5 with the same Z input of 7.90559. Other values of Z can cause the problem as well, but have not been included here for brevity. Note that not all combinations of P and Z give invalid results, just some do, and the behavior is always that the sign is wrong (the absolute value is correct). Using the --fast-math flag on older versions of GCC produces the correct result.

[Dustin Lanham <lanhamd@…>]

Addendum: This behavior was noted in two separate versions of boost: 1.50.0 and 1.55.0.

[John Maddock]

To a degree, this is not our bug, a reduced test case is simply:

   long double ld = -5.3165867831218916301793863361917824e-2467L;
   std::cout << std::signbit(ld) << std::endl;

Which incorrectly outputs 0 when compiling with --fast-math.

I suspect there may be other std lib functions which have similar issues too, I'll investigate working around this one, but you may just get whacked by another...

[John Maddock]

Fix is here: ​https://github.com/boostorg/math/commit/9f8ffee4b7a3f82b1c582735d43522d7d0cde746

As an alternative, define BOOST_MATH_DISABLE_STD_FPCLASSIFY when building.