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

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

[vincent-srcware at vinc17 dot net]

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

Vincent Lefèvre <vincent-srcware at vinc17 dot net> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
            Summary|pow() produces inaccurate   |pow() produces inaccurate
                   |results for base ~ 1.0, and |results for base ~ 1.0, and
                   |large exponent              |large exponent on 32-bit
                   |                            |x86

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jsm28 at gcc dot gnu.org]

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

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|SUSPENDED                   |NEW
         AssignedTo|aj at suse dot de           |unassigned at sourceware
                   |                            |dot org

--- Comment #7 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-02-22 21:29:35 UTC ---
Still present.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[jsm28 at gcc dot gnu.org]

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

--- Comment #8 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-02-22 21:52:38 UTC ---
I should add: part of the ABI is that the libm functions are expected to be
called with the FPU in 64-bit mode, so the results with it in 53-bit mode are
not a bug; it's only the 64-bit mode case where we ought to do better.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[calixte.denizet@scilab-enterprises.com]

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

--- Comment #9 from Calixte <calixte.denizet@scilab-enterprises.com> 2012-02-25 16:06:18 UTC ---
Created attachment 6245
  --> http://sourceware.org/bugzilla/attachment.cgi?id=6245
C code with a better algorithm

Hi all,

When the exponent is an integer, a fast exponentiation is used and when it is
negative the base is inverted before (cf e_pow.S). So in this case the error
induced by the inversion is cumulated in squaring several times. In the present
case base=(2^53-1)*2^(-53), exponent=2^54 and 1/base=(2^52+1)*2^(-52),
classical approximation shows that (1/base)^(2^54))~exp(4) and
1/(base^(2^54))~exp(2). So an easy way to avoid to cumulated error induced by
inversion is to invert the result rather than invert the base.

The attached code (in C) shows that clearly. I'm not able to fix the assembly
code in e_pow.S, so I let experimented people do that.

Thanks

Calixte

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[vincent-srcware at vinc17 dot net]

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

--- Comment #10 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> 2012-02-27 15:40:08 UTC ---
(In reply to comment #9)
> Created attachment 6245 [details]
> C code with a better algorithm

Note concerning your code: it assumes that char is signed. You should use int
for small integers, not char.

> When the exponent is an integer, a fast exponentiation is used
[...]

This is not perfect (you still won't get a very accurate result), but certainly
better than the current code (I haven't looked at it), which is apparently
affected by a cancellation. Note that if you want to compute pow(x,y) where y
is not necessarily an integer, you can write y = yi + yf, where yi is an
integer and |yf| < 1 (e.g. |yf| <= 0.5 or 0 <= yf < 1). Then mathematically,
pow(x,y) = pow(pow(x,yi),yf), and you can use that to get a more accurate
result than the current algorithm (I haven't checked in detail...).

For more accurate algorithms concerning the integer powers:

@Article{KLLLM2009a,
        AUTHOR          = {P. Kornerup and Ch. Lauter and V. Lefèvre and N.
Louvet and J.-M. Muller},
        TITLE           = {Computing Correctly Rounded Integer Powers in
Floating-Point Arithmetic},
        JOURNAL         = {{ACM} Transactions on Mathematical Software},
        VOLUME          = {37},
        NUMBER          = {1},
        MONTH           = jan,
        YEAR            = {2010},
        URL             = {http://doi.acm.org/10.1145/1644001.1644005},
        DOI             = {10.1145/1644001.1644005}
}

and for correct rounding of the general pow() function, assuming you already
have an accurate algorithm:

@Article{LauLef2009a,
        AUTHOR          = {Ch. Lauter and V. Lefèvre},
        TITLE           = {An efficient rounding boundary test for
\texttt{pow(x,y)} in double precision},
        JOURNAL         = {{IEEE} Transactions on Computers},
        PUBLISHER       = {{IEEE} Computer Society Press, Los Alamitos, CA},
        VOLUME          = {58},
        NUMBER          = {2},
        PAGES           = {197--207},
        MONTH           = feb,
        YEAR            = {2009},
        KEYWORDS        = {floating-point arithmetic,correct rounding,power
function},
        URL             =
{http://www.vinc17.net/research/papers/ieeetc2009-powr.pdf},
        DOI             = {10.1109/TC.2008.202}
}

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[calixte.denizet@scilab-enterprises.com]

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

Calixte <calixte.denizet@scilab-enterprises.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |calixte.denizet@scilab-ente
                   |                            |rprises.com

--- Comment #11 from Calixte <calixte.denizet@scilab-enterprises.com> 2012-02-28 12:34:51 UTC ---
Clearly the pow alogrithm should be reimplemented with a better one.
What I proposed is just an easy workaround which can be easily and quickly
implemented (and it would improve the actual results).
About the char in the C code: it is just used for its sign (+/- 1) not for its
value.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[vincent-srcware at vinc17 dot net]

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

--- Comment #12 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> 2012-02-28 14:03:00 UTC ---
(In reply to comment #11)
> What I proposed is just an easy workaround which can be easily and quickly
> implemented (and it would improve the actual results).

I'm not sure that improving the code for integer exponents only would be a good
idea. A consequence would be that the f(x) = pow(a,x) implementation would be
completely non-monotonous (with important differences near the integer inputs),
with possible bad side effects on some codes.

> About the char in the C code: it is just used for its sign (+/- 1) not for its
> value.

But a sign is part of the value. On platforms where char is unsigned, a char
value cannot have a negative sign.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

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

[vincent-srcware at vinc17 dot net]

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

--- Comment #13 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> 2012-02-28 14:57:42 UTC ---
(In reply to comment #10)
> Note that if you want to compute pow(x,y) where y is not necessarily an
> integer, you can write y = yi + yf, where yi is an integer and |yf| < 1
> (e.g. |yf| <= 0.5 or 0 <= yf < 1). Then mathematically,
> pow(x,y) = pow(pow(x,yi),yf), and you can use that to get a more accurate
> result than the current algorithm (I haven't checked in detail...).

This is incorrect. I should have re-read my message. There are 2 possible
simple range reductions that wouldn't require much rewriting, in order to
reduce the problem to integer exponents (if one is able to fix that easily).

1. y = yi + yf, where yi is an integer and |yf| < 1. Then mathematically,
pow(x,y) = pow(x,yi) * pow(x,yf), where pow(x,yi) can be computed with an
iterated exponentiation. Note that yi and yf are computed exactly.

2. y = y0 * 2^e, where 0.5 <= |y0| < 1 for instance (which can be obtained with
frexp()). Then mathematically, pow(x,y) = pow(pow(x,y0),2^e), where the ^(2^e)
can be computed with repeated squarings.

Now, I fear that if the iterated exponentiation is not done in multiple
precision, this won't solve the problem for very large values of n. Indeed,
pow(x,n) evaluated with such an algorithm (whatever iteration is used) is
approximated with a factor of about[*] (1+epsilon)^(n-1) > 1 + (n-1).epsilon,
where 1+epsilon is the error bound on the multiplication in Higham's notation.
And this would be slow.

[*] Actually this is an upper bound, but in average, you would probably get the
same kind of problems with large values of n.

And Calixte, I've looked at your comment #9 and your code more closely, and I
agree that in the case of a negative exponent, the inversion should be done
after the exponentiation if such an algorithm is used.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/706] pow() produces inaccurate results for base ~ 1.0, and large integer exponent on 32-bit x86

[vincent-srcware at vinc17 dot net]

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

Vincent Lefèvre <vincent-srcware at vinc17 dot net> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
            Summary|pow() produces inaccurate   |pow() produces inaccurate
                   |results for base ~ 1.0, and |results for base ~ 1.0, and
                   |large exponent on 32-bit    |large integer exponent on
                   |x86                         |32-bit x86

--- Comment #14 from Vincent Lefèvre <vincent-srcware at vinc17 dot net> 2012-02-28 16:36:13 UTC ---
I've done some tests, and it seems that pow(x,y) for large y is inaccurate only
when y is an integer. So, if I understand correctly, an iterative algorithm is
used for integers y, and a log-exp algorithm (which seems accurate enough on
the few values I've tested) is used for other values. So, a possible easy fix
would be to always use the log-exp algorithm, possibly except for small
integers y.

Note: I've updated the bug summary following this observation (large exponent →
large integer exponent).

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/706] pow() produces inaccurate results for base ~ 1.0, and large integer exponent on 32-bit x86

[jsm28 at gcc dot gnu.org]

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

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |FIXED

--- Comment #15 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-04-09 09:47:09 UTC ---
Fixed by:

commit c483f6b4a4277bc209820efc1ae35d976af57b4e
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Mon Apr 9 09:42:05 2012 +0000

    Fix x86 pow inaccuracy for large integer exponents (bug 706).

(A similar issue for x86/x86_64 powl is tracked in bug 13881.)

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/706] pow() produces inaccurate results for base ~ 1.0, and large integer exponent on 32-bit x86

[jackie.rosen at hushmail dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=706

Jackie Rosen <jackie.rosen at hushmail dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jackie.rosen at hushmail dot com

--- Comment #16 from Jackie Rosen <jackie.rosen at hushmail dot com> ---
*** Bug 260998 has been marked as a duplicate of this bug. ***
Seen from the domain http://volichat.com
Page where seen: http://volichat.com/adult-chat-rooms
Marked for reference. Resolved as fixed @bugzilla.

-- 
You are receiving this mail because:
You are on the CC list for the bug.

[Bug math/706] pow() produces inaccurate results for base ~ 1.0, and large integer exponent on 32-bit x86

[schwab at sourceware dot org]

https://sourceware.org/bugzilla/show_bug.cgi?id=706

Andreas Schwab <schwab at sourceware dot org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|jackie.rosen at hushmail dot com   |

-- 
You are receiving this mail because:
You are on the CC list for the bug.