Bug 2445

Summary: Bug in powl or expl?
Product: glibc Reporter: John Maddock <john>
Component: mathAssignee: Ulrich Drepper <drepper.fsp>
Status: RESOLVED WORKSFORME    
Severity: normal CC: glibc-bugs
Priority: P2 Flags: fweimer: security-
Version: unspecified   
Target Milestone: ---   
Host: Linux 2.6.5-7.252-smp x86_64 Target: Linux 2.6.5-7.252-smp x86_64
Build: Linux 2.6.5-7.252-smp x86_64 Last reconfirmed:

Description John Maddock 2006-03-11 19:13:02 UTC
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.
Comment 1 Ulrich Drepper 2006-04-25 20:33:21 UTC
If you have accuracy problems with libm, send a patch.  Until then the bugs is
suspended.
Comment 2 Joseph Myers 2012-02-23 01:46:37 UTC
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.