706
2005-02-02 20:39:13 +0000
pow() produces inaccurate results for base ~ 1.0, and large integer exponent on 32-bit x86
2019-04-10 08:47:58 +0000
1
1
1
Unclassified
glibc
math
2.3.4
All
All
RESOLVED
FIXED
P2
normal
---
1
todd.allen
unassigned
bugsy
calixte.denizet
glibc-bugs
hjl.tools
olecom
vincent-srcware
oldest_to_newest
2728
0
todd.allen
2005-02-02 20:39:13 +0000
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.
2729
1
395
todd.allen
2005-02-02 20:40:19 +0000
Created attachment 395
loss_of_power.c
3657
2
olecom
2005-05-21 19:27:06 +0000
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.
5108
3
aj
2005-09-16 12:26:52 +0000
The current glibc routines for pow are not based on fdlibm anymore.
5309
4
drepper.fsp
2005-09-25 23:53:53 +0000
Provide a patch.
9768
5
drepper.fsp
2006-04-25 22:52:17 +0000
Suspended until somebody provides a patch.
38479
6
vincent-srcware
2009-09-24 08:07:02 +0000
*** Bug 2579 has been marked as a duplicate of this bug. ***
53429
7
jsm28
2012-02-22 21:29:35 +0000
Still present.
53432
8
jsm28
2012-02-22 21:52:38 +0000
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.
53551
9
6245
calixte.denizet
2012-02-25 16:06:18 +0000
Created attachment 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
53579
10
vincent-srcware
2012-02-27 15:40:08 +0000
(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}
}
53614
11
calixte.denizet
2012-02-28 12:34:51 +0000
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.
53615
12
vincent-srcware
2012-02-28 14:03:00 +0000
(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.
53616
13
vincent-srcware
2012-02-28 14:57:42 +0000
(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.
53620
14
vincent-srcware
2012-02-28 16:36:13 +0000
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).
54497
15
jsm28
2012-04-09 09:47:09 +0000
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.)
66693
16
jackie.rosen
2014-02-16 18:27:22 +0000
*** 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.
395
2005-02-02 20:40:19 +0000
2005-02-02 20:40:19 +0000
loss_of_power.c
it.c
text/plain
559
todd.allen
LyoKKiogbG9zc19vZl9wb3dlcgoqKgoqKiBOT1RFOiBpZGVhbCB2YWx1ZSBwcm9kdWNlZCB3aXRo
OgoqKiAgICAgICBjYWxjICdlcHNpbG9uKDFlLTUwKTsgZXhwKGxuKDEgKyAoLTEuMTEwMjIzMDI0
NjI1MTZFLTE2KSkgKiAoLTEuODAxNDM5ODUwOTQ4MjBFKzE2KSknCiovCgojaW5jbHVkZSA8bWF0
aC5oPgoKc3RhdGljIHZvaWQgc2V0X2ZwdV81Myh2b2lkKQp7CiAgIF9fYXNtX18gX192b2xhdGls
ZV9fKCJwdXNobCAkMCIpOwogICBfX2FzbV9fIF9fdm9sYXRpbGVfXygiZnN0Y3cgMCglZXNwKSIp
OwogICBfX2FzbV9fIF9fdm9sYXRpbGVfXygiYW5kdyAgJDB4ZmNmZiwgMCglZXNwKSIpOyAvKiBT
ZXQgUEM9MCAgICAgICAgICAgICovCiAgIF9fYXNtX18gX192b2xhdGlsZV9fKCJvcncgICAkMHgw
MjAwLCAwKCVlc3ApIik7IC8qIFNldCBQQz0xMCAoNTMgYml0cykgKi8KICAgX19hc21fXyBfX3Zv
bGF0aWxlX18oImZsZGN3IDAoJWVzcCkiKTsKICAgX19hc21fXyBfX3ZvbGF0aWxlX18oInBvcGwg
LTQoJWVzcCkiKTsKfQoKc3RhdGljIHZvaWQgc2V0X2ZwdV82NCh2b2lkKQp7CiAgIF9fYXNtX18g
X192b2xhdGlsZV9fKCJwdXNobCAkMCIpOwogICBfX2FzbV9fIF9fdm9sYXRpbGVfXygiZnN0Y3cg
MCglZXNwKSIpOwogICBfX2FzbV9fIF9fdm9sYXRpbGVfXygib3J3ICAgJDB4MDMwMCwgMCglZXNw
KSIpOyAvKiBTZXQgUEM9MTEgKDY0IGJpdHMpICovCiAgIF9fYXNtX18gX192b2xhdGlsZV9fKCJm
bGRjdyAwKCVlc3ApIik7CiAgIF9fYXNtX18gX192b2xhdGlsZV9fKCJwb3BsIC00KCVlc3ApIik7
Cn0KCm1haW4oKQp7CiAgIGRvdWJsZSAgYmFzZSAgICAgPSAxLjAgLSAxLjExMDIyMzAyNDYyNTE2
RS0xNjsKICAgZG91YmxlICBleHBvbmVudCA9IC0xLjgwMTQzOTg1MDk0ODIwRSsxNjsKICAgZG91
YmxlICBpZGVhbCAgICA9IDcuMzg5MDU2MDk4OTMwNzEwMjIzMzk7CiAgIGRvdWJsZSAgcGFyYW5v
aWEgPSA3LjM4OTA1NjA5ODkzMDY1OwoKICAgcHJpbnRmKCJwb3coJTE5LjE2bGYsICUxLjFsZilc
biIsIGJhc2UsIGV4cG9uZW50KTsKCiAgIHByaW50ZigiICAgaWRlYWwgICAgICAgICAgfiAlMTku
MTZsZiAoZGlmZiA9ICVsZylcbiIsIAogICAgICAgICAgaWRlYWwsIGlkZWFsIC0gcGFyYW5vaWEp
OwogICBwcmludGYoIiAgIHBhcmFub2lhICAgICAgID0gJTE5LjE2bGZcbiIsIAogICAgICAgICAg
cGFyYW5vaWEpOwoKICAgewogICAgICBkb3VibGUgIHJlc3VsdCA9IGV4cDIobG9nMihiYXNlKSAq
IGV4cG9uZW50KTsKICAgICAgcHJpbnRmKCIgICA2NC1iaXQgbG9nL2V4cCA9ICUxOS4xNmxmIChk
aWZmID0gJWxnKVxuIiwgCiAgICAgICAgICAgICByZXN1bHQsIHJlc3VsdCAtIHBhcmFub2lhKTsK
ICAgfQogICB7CiAgICAgIGRvdWJsZSAgcmVzdWx0ID0gcG93KGJhc2UsIGV4cG9uZW50KTsKICAg
ICAgcHJpbnRmKCIgICA2NC1iaXQgcG93ICAgICA9ICUxOS4xNmxmIChkaWZmID0gJWxnKVxuIiwg
CiAgICAgICAgICAgICByZXN1bHQsIHJlc3VsdCAtIHBhcmFub2lhKTsKICAgfQoKICAgc2V0X2Zw
dV81MygpOwogICB7CiAgICAgIGRvdWJsZSAgcmVzdWx0ID0gZXhwMihsb2cyKGJhc2UpICogZXhw
b25lbnQpOwogICAgICBwcmludGYoIiAgIDUzLWJpdCBsb2cvZXhwID0gJTE5LjE2bGYgKGRpZmYg
PSAlbGcpXG4iLCAKICAgICAgICAgICAgIHJlc3VsdCwgcmVzdWx0IC0gcGFyYW5vaWEpOwogICB9
CiAgIHsKICAgICAgZG91YmxlICByZXN1bHQgPSBwb3coYmFzZSwgZXhwb25lbnQpOwogICAgICBw
cmludGYoIiAgIDUzLWJpdCBwb3cgICAgID0gJTE5LjE2bGYgKGRpZmYgPSAlbGcpXG4iLCAKICAg
ICAgICAgICAgIHJlc3VsdCwgcmVzdWx0IC0gcGFyYW5vaWEpOwogICB9CiAgIHNldF9mcHVfNjQo
KTsKfQo=
6245
2012-02-25 16:06:18 +0000
2012-02-25 16:06:18 +0000
C code with a better algorithm
plop.c
text/x-csrc
597
calixte.denizet
I2luY2x1ZGUgPG1hdGguaD4KI2luY2x1ZGUgPHN0ZGlvLmg+CgpzdGF0aWMgdm9pZCBzZXRfZnB1
XzUzKHZvaWQpCnsKICAgX19hc21fXyBfX3ZvbGF0aWxlX18oInB1c2hsICQwIik7CiAgIF9fYXNt
X18gX192b2xhdGlsZV9fKCJmc3RjdyAwKCVlc3ApIik7CiAgIF9fYXNtX18gX192b2xhdGlsZV9f
KCJhbmR3ICAkMHhmY2ZmLCAwKCVlc3ApIik7IC8qIFNldCBQQz0wICAgICAgICAgICAgKi8KICAg
X19hc21fXyBfX3ZvbGF0aWxlX18oIm9ydyAgICQweDAyMDAsIDAoJWVzcCkiKTsgLyogU2V0IFBD
PTEwICg1MyBiaXRzKSAqLwogICBfX2FzbV9fIF9fdm9sYXRpbGVfXygiZmxkY3cgMCglZXNwKSIp
OwogICBfX2FzbV9fIF9fdm9sYXRpbGVfXygicG9wbCAtNCglZXNwKSIpOwp9CgpzdGF0aWMgdm9p
ZCBzZXRfZnB1XzY0KHZvaWQpCnsKICAgX19hc21fXyBfX3ZvbGF0aWxlX18oInB1c2hsICQwIik7
CiAgIF9fYXNtX18gX192b2xhdGlsZV9fKCJmc3RjdyAwKCVlc3ApIik7CiAgIF9fYXNtX18gX192
b2xhdGlsZV9fKCJvcncgICAkMHgwMzAwLCAwKCVlc3ApIik7IC8qIFNldCBQQz0xMSAoNjQgYml0
cykgKi8KICAgX19hc21fXyBfX3ZvbGF0aWxlX18oImZsZGN3IDAoJWVzcCkiKTsKICAgX19hc21f
XyBfX3ZvbGF0aWxlX18oInBvcGwgLTQoJWVzcCkiKTsKfQoKZG91YmxlIG15cG93KGRvdWJsZSB4
LCB1bnNpZ25lZCBsb25nIGxvbmcgeSwgY2hhciB5c2lnbikKewogICAgZG91YmxlIHN0ID0gMTsK
CiAgICAvLyB5IGlzIHN1cHBvc2VkIHRvIGJlIGRpZmZlcmVudCBvZiB6ZXJvCiAgICAvLyBzaW5j
ZSB0aGlzIGNhc2UgaGFzIGJlZW4gdHJlYXRlZCBiZWZvcmUgCiAgICBkbwogICAgewoJaWYgKHkg
JiAxKQoJewoJICAgIHN0ICo9IHg7Cgl9Cgl4ICo9IHg7Cgl5ID4+PSAxOwogICAgfQogICAgd2hp
bGUgKHkpOwogICAgCiAgICBpZiAoeXNpZ24gPiAwKQogICAgewoJcmV0dXJuIHN0OwogICAgfQog
ICAgZWxzZQogICAgewoJcmV0dXJuIDEgLyBzdDsKICAgIH0KfQoKaW50IG1haW4oKSAKewogICAg
ZG91YmxlIGJhc2UgPSAxLjAgLSAxLjExMDIyMzAyNDYyNTE2RS0xNjsKICAgIGRvdWJsZSBleHBv
bmVudCA9IC0xLjgwMTQzOTg1MDk0ODIwRSsxNjsKICAgIGRvdWJsZSBpZGVhbCA9IDcuMzg5MDU2
MDk4OTMwNzEwMjIzMzk7CiAgICAKICAgIHNldF9mcHVfNTMoKTsKCiAgICBkb3VibGUgbmV3ID0g
bXlwb3coYmFzZSwgKHVuc2lnbmVkIGxvbmcgbG9uZykoLWV4cG9uZW50KSwgLTEpOwogICAgZG91
YmxlIG9sZCA9IG15cG93KDEgLyBiYXNlLCAodW5zaWduZWQgbG9uZyBsb25nKSgtZXhwb25lbnQp
LCAxKTsKCiAgICBwcmludGYoIm5ldyBhbGdvOiByZWxhdGl2ZSBlcnJvcj0lbGcgYW5kIGFic29s
dXRlIGVycm9yPSVsZ1xuIiwgZmFicyhuZXcgLSBpZGVhbCkgLyBpZGVhbCwgZmFicyhuZXcgLSBp
ZGVhbCkpOwogICAgcHJpbnRmKCJvbGQgYWxnbzogcmVsYXRpdmUgZXJyb3I9JWxnIGFuZCBhYnNv
bHV0ZSBlcnJvcj0lbGdcbiIsIGZhYnMob2xkIC0gaWRlYWwpIC8gb2xkLCBmYWJzKG9sZCAtIGlk
ZWFsKSk7CiAgICAKICAgIHNldF9mcHVfNjQoKTsKCiAgICByZXR1cm4gMDsKfQo=