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

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

[todd dot allen at attglobal dot net]

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.

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

[todd dot allen at attglobal dot net]

------- Additional Comments From todd dot allen at attglobal dot net  2005-02-02 20:40 -------
Created an attachment (id=395)
 --> (http://sources.redhat.com/bugzilla/attachment.cgi?id=395&action=view)
loss_of_power.c


-- 


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.

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

[bugsy at ccur dot com]

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |bugsy at ccur 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.

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

[olecom at gmail dot com]

------- Additional Comments From olecom at gmail dot com  2005-05-21 19:27 -------
Bugs with pow() are mentioned to be fixed in fdlibm 5.3

<ftp://ftp.netlib.org/fdlibm/changes>

so, maybe if will be easy to fix them, glibc was based on fdlibm.


-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |olecom at gmail 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.

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

[aj at suse dot de]

------- Additional Comments From aj at suse dot de  2005-09-16 12:26 -------
The current glibc routines for pow are not based on fdlibm anymore.

-- 


http://sourceware.org/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.

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

[drepper at redhat dot com]

------- Additional Comments From drepper at redhat dot com  2005-09-25 23:53 -------
Provide a patch.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |WAITING


http://sourceware.org/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.

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

[drepper at redhat dot com]

------- Additional Comments From drepper at redhat dot com  2006-04-25 22:52 -------
Suspended until somebody provides a patch.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|WAITING                     |SUSPENDED


http://sourceware.org/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.

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

[vincent+libc at vinc17 dot org]

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |vincent+libc at vinc17 dot
                   |                            |org


http://sourceware.org/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.

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

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2009-09-24 08:07 -------
*** Bug 2579 has been marked as a duplicate of this bug. ***

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |hjl dot tools at gmail dot
                   |                            |com


http://sourceware.org/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.