[Bug math/706] New: pow() produces inaccurate results for base ~ 1.0, and large exponent

["todd dot allen at attglobal dot net" <sourceware-bugzilla at sources dot redhat dot com>]

The libm version of pow() is producing inaccurate results when the base is very
close to 1.0 and the exponent has a very large magnitude.  Once I get this bug
reported, I'll attach an example source that demonstrates the problem.  The
constants used in this example come from the paranoia test program.

I'm running this on Fedora Core 3 with the glibc-2.3.4-2.fc3 rpm.  The output of
the test program is:

pow( 0.9999999999999999, -18014398509482000.0)
   ideal          ~  7.3890560989307099 (diff = 5.9508e-14)
   paranoia       =  7.3890560989306504
   64-bit log/exp =  7.3890560989306637 (diff = 1.33227e-14)
   64-bit pow     =  7.3890560954893578 (diff = -3.44129e-09)
   53-bit log/exp =  7.3890560989306637 (diff = 1.33227e-14)
   53-bit pow     = 54.5981484040811011 (diff = 47.2091)

The "ideal" value was computed using the arbitrary-precision calculator calc
(from ftp://ftp.uu.net/pub/calc).  The value expected by paranoia is reported
next, and it's pretty close (deviation of about 6e-14).  The values produced by
log,*,exp, are similarly close.  The values produced by pow are the problematic
ones.

Fir the 64-bit pow case, the deviation from the paranoia value is around 3e-9,
or at least 10,000 times too much deviation.

I've also included code that sets the x86 rounding mode to 53-bit precision. 
The log,*,exp code still produces a pretty good value.  The pow function call's
error gets really out of hand, producing a result which doesn't even have the
right order of magnitude.  I don't know if libm is supposed to be able to
support the x87 53-bit precision rounding mode, but it would be nice if the
values were a bit closer than this.  :)

Regardless, the error in the default 64-bit rounding mode is what I'm reporting,
at least primarily.

I see that the ieee754_pow routine is using an iterative approach to compute the
value, because the exponent is an integer.  Note that the exponent is an integer
only because it's so large that all of the mantissa bits are over in the
integral part of the value.  I suspect that this approach just accumulates too
much multiplication error.  However, I also suspect that it was coded that way
for a good reason, presumably because there are cases where it produces less
error than the log,*,exp approach.  So, perhaps there should just be some
additional condition that precludes that code from being used for these values.
 Maybe some a cutoff if the base is very close to 1.0.  I'm not sure.

I ran across one other discussion of a bug very similar to this one by some
python maintainers at the following URL's:

https://sourceforge.net/tracker/?func=detail&atid=105470&aid=705231&group_id=5470
http://mail.python.org/pipermail/python-bugs-list/2003-May/017987.html

They mentioned that they intended to report the issue, but perhaps they didn't
do so.

-- 
           Summary: pow() produces inaccurate results for base ~ 1.0, and
                    large exponent
           Product: glibc
           Version: 2.3.4
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: aj at suse dot de
        ReportedBy: todd dot allen at attglobal dot net
                CC: glibc-bugs at sources dot redhat dot com


http://sources.redhat.com/bugzilla/show_bug.cgi?id=706

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