Bug in boost::math::special_functions::sign showing up in special_functions::bessel_jy()

[peter zwart <PHZwart@…>]

Hi,

When computing scattering curves, I encountered a show-stopping bug that I tracked back to the boost::math::special_functions::sign function .

In math::special_functions::bessel_jy a small number is generated via

sqrt(tools::min_value<T>())

which subsequently ends up in this statement: sign(current)

I see this happening:

variable current: -1.24911e-2468 variable sign(current): 1

which causes problems for math::special_functions::sph_bessel(1,x) for x between 7.845 and 8 or so on my machines.

The fix that worked for me is the following

T sign_current; sign_current = std::fabs(current); if (current*2.0 < sign_current){

sign_current = -1.0;

} else {

sign_current = 1.0;

}

Ju = sign_current * sqrt(W / (q + gamma * (p - t))); Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));

For now, I'll run with this local fix but it seems that on some platforms the current sign function causes problems when computing spherical Bessel functions of any order above 0.

[anonymous]

Hi Peter,

I'm unable to reproduce here, so I need your help as the implementation of boost::math::sign is platform/compiler dependent, so I have a bunch of questions for you I'm afraid:

1) What compiler / platform / architecture is this? 2) Can check the return value of boost::math::sign(-1.24911e-2468L); and verify that this returns an incorrect result? 3) If (2) does return an incorrect result, can you please step though it's code and let me know which of the implementations of signbit_impl in sign.hpp is in use?

Many thanks, John Maddock.

[phzwart@…]

Hi,

1.

uname -a gives

Linux xxx.yyy.zzz.gov 2.6.33 #1 SMP Mon Mar 1 16:05:57 PST 2010 x86_64 x86_64 x86_64 GNU/Linux

localhost 8% g++ -v Using built-in specs. Target: x86_64-redhat-linux Configured with: ../configure --prefix=/usr --mandir=/usr/share/man --infodir=/usr/share/info --with-bugurl=​http://bugzilla.redhat.com/bugzilla --enable-bootstrap --enable-shared --enable-threads=posix --enable-checking=release --with-system-zlib --enable-cxa_atexit --disable-libunwind-exceptions --enable-languages=c,c++,objc,obj-c++,java,fortran,ada --enable-java-awt=gtk --disable-dssi --enable-plugin --with-java-home=/usr/lib/jvm/java-1.5.0-gcj-1.5.0.0/jre --enable-libgcj-multifile --enable-java-maintainer-mode --with-ecj-jar=/usr/share/java/eclipse-ecj.jar --disable-libjava-multilib --with-cpu=generic --build=x86_64-redhat-linux Thread model: posix gcc version 4.3.2 20081105 (Red Hat 4.3.2-7) (GCC)

2.

fault is there:

-1.24911e-2468 1

3.

this one is used:

#ifdef BOOST_MATH_USE_STD_FPCLASSIFY

template<class T> inline int signbit_impl(T x, native_tag const&) {

return (std::signbit)(x);

}

#endif

[anonymous]

Thanks for the info,

So am I correct in thinking that (std::signbit)(-1.24911e-2468L) returns 0? If so can you please file a bug report with the GNU libc guys?

The fix will be to define BOOST_MATH_DISABLE_STD_FPCLASSIFY for those platforms that are effected... if we can figure out which those are. In the mean time I'll update our regression tests with this case and see what fails...

Regards, John Maddock.

[anonymous]

it returns 1.

if it is a gnu libc issue, i'll file a bug report there as well.

P

[anonymous]

it returns 1.

Then I'm confused as to why it's failing, are you in a position to debug

boost::math::sign(-1.24911e-2468L);

And figure out why it isn't returning -1? The code is so simple it's hard to see how it could go wrong really and I still can't reproduce here...

Thanks, John.

[anonymous]

#include <iostream> #include <cmath> int main(void){

std::cout << "hello world " << std::endl;

std::cout << (std::signbit)(-1.24911e-2468L)<<std::endl;;

std::cout << (std::signbit)(-2.24911e-2468L)<<std::endl;;

std::cout << (std::signbit)(-3.24911)<<std::endl;

std::cout << (std::signbit)(3.24911)<<std::endl;

}

gives

1 1 1 0

Not much time for debugging this right at the moment, i'll try tonight.

[John Maddock]

Any news on this?

I've updated our regression tests with both the sph_bessel test case you give, and the sign case, and both seem to be passing on our test machines.

[anonymous]

No, sorry not time yet to dig further. I hope to get to it in the next couple of days.

Sorry,

P

[John Maddock]

Ping?

[John Maddock]

I need more info to deal with this, so I'm closing for now, please reopen if you can provide the needed info.