I have a strange problem with powl that's very specific to AMD64 SMP machines (Suse Linux 9 on HP 2x Opteron, the actual machine is one of the HP testdrive machines: td190.testdrive.hp.com, so you can even try it on the same machine if you need to). Basically I have an expression that should underflow to zero, but which yields a non-zero result *only* on the second or subsequent execution of the expression. The gcc version info is: maddock@td190:~/libs/math/test> g++ -v Reading specs from /usr/lib64/gcc-lib/x86_64-suse-linux/3.3.3/specs Configured with: ../configure --enable-threads=posix --prefix=/usr --with-local- prefix=/usr/local --infodir=/usr/share/info --mandir=/usr/share/man --enable-lan guages=c,c++,f77,objc,java,ada --disable-checking --libdir=/usr/lib64 --enable-l ibgcj --with-gxx-include-dir=/usr/include/g++ --with-slibdir=/lib64 --with-syste m-zlib --enable-shared --enable-__cxa_atexit x86_64-suse-linux Thread model: posix gcc version 3.3.3 (SuSE Linux) glibc version is 20040214, not the latest by any means, but unfortunately I'm unable to test a newer version on this machine. The code to reproduce the issue is: #include <math.h> #include <iostream> template <class T, class L> T regularised_gamma_prefix(T a, T z, const L&) { using namespace std; T agh = a + 4.6875L - 0.5L; //T prefix; //T d = ((z - a) - 4.6875L + 0.5L) / agh; return powl((z * expl((a-z)/a)) / agh, a); } template <class T> void test(T) { T a = 1858.277099609375; T z = 185827.703125; T r; for(int i = 0; i < 2; ++i) { r = regularised_gamma_prefix( a, z, int()); std::cout << r << " " << i << std::endl; } } int main() { test((long double)(0)); return 0; } which when run yields: 0 0 4.00731e-560 1 and on every other platform I've tried yields as expected: 0 0 0 1 My apologies if this has been fixed or reported already, but I haven't been able to find anything similar reported. Thanks in advance for any insights you may have, Regards, John Maddock.
If you have accuracy problems with libm, send a patch. Until then the bugs is suspended.
Given the description it seems unlikely this is a glibc bug (a compiler bug might be a possibility, or a kernel bug relating to floating point state), and the test works for me. I don't think it can meaningfully be investigated without information from someone with a system where it shows up about the exact point in execution at which the two computations diverge - where they cease to have the same instructions executed on the same arguments.