Bug 2445

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.