[Bug libc/2445] New: Bug in powl or expl?

[Bug libc/2445] New: Bug in powl or expl?

[john at johnmaddock dot co dot uk]

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.

-- 
           Summary: Bug in powl or expl?
           Product: glibc
           Version: unspecified
            Status: NEW
          Severity: normal
          Priority: P2
         Component: libc
        AssignedTo: drepper at redhat dot com
        ReportedBy: john at johnmaddock dot co dot uk
                CC: glibc-bugs at sources dot redhat dot com
 GCC build triplet: Linux 2.6.5-7.252-smp x86_64
  GCC host triplet: Linux 2.6.5-7.252-smp x86_64
GCC target triplet: Linux 2.6.5-7.252-smp x86_64


http://sourceware.org/bugzilla/show_bug.cgi?id=2445

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug libc/2445] Bug in powl or expl?

[drepper at redhat dot com]

------- Additional Comments From drepper at redhat dot com  2006-04-25 20:33 -------
If you have accuracy problems with libm, send a patch.  Until then the bugs is
suspended.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |SUSPENDED


http://sourceware.org/bugzilla/show_bug.cgi?id=2445

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/2445] Bug in powl or expl?

[pasky at suse dot cz]

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
          Component|libc                        |math


http://sourceware.org/bugzilla/show_bug.cgi?id=2445

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.