-
-
Notifications
You must be signed in to change notification settings - Fork 31.7k
Add math.exp2() function: 2^x #90075
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Comments
Dear Python Support Team, I was looking through Python’s list of supported methods in the math module, and I noticed that C99’s exp2 method was not implemented. This method raises 2 to the power of the supplied argument. I understand that it’s pretty trivial to so this in Python using 2**x or math.pow(x, 2), but I think there are a few reasons why we might want to incorporate it: Uniformity: This method exists most other programming languages and libraries, including numpy. Consistency: Every math method from C99 except exp2 is in python’s math or cmath module (math.cbrt will be added as of python 3.11). Triviality: this method is a part of C99 and is also supported by Visual Studio, so it’s very easy to implement. Accuracy(?): a libm exp2 is supposedly more accurate than pow(2.0, x), though I don’t really see how this would be the case (See https://bugs.python.org/issue31980) That said, this method is a little redundant, so I completely understand if this request is rejected Non-exhaustive list of other languages / libraries that use this method: Rust: https://docs.rs/libm/0.1.1/libm/fn.exp2.html Similar Issues: |
Sounds good to me, provided that all the common platforms that we care about have a reasonable quality implementation. This should be a straightforward wrapping of the C99 function, and with sufficient tests the buildbots should tell us if there are any issues on common platforms. @Gideon: are you're interested in working on a pull request? I'd be happy to review. (Ideally I'd like to have exp10 too, but that's not in C99 so platform support is likely to be spotty. If anyone's interested in pursuing that, we should make it a separate issue.)
pow is a difficult function to implement at high accuracy, and there are a good number of low quality pow implementations around in system math libraries. It's much easier to come up with a high accuracy implementation of a single-argument function - there are well known techniques for generating approximating polynomials that simply don't extend well to functions of two arguments. sqrt is similar: pow(x, 0.5) is very often not correctly rounded even on systems where sqrt(x) _is_. (Though that one's a bit of a cheat, since common processors have dedicated instructions for a correctly-rounded sqrt.) |
See also previous discussion towards the end of https://bugs.python.org/issue3366. FWIW, I don't think there's value in adding exp2 to the cmath module too: we'd have to write our own implementation, and it's just not a function that appears often in the complex world. |
On the subject of accuracy, there doesn't seem to be much in it on my mac laptop, and it looks as though pow(2.0, x) is giving correctly rounded results as often as (if not more often than) exp2(x). Here's the log of a terminal session, after recompiling Python to add exp2. It shows the ulps error (tested against a high-precision Decimal computation, which we're treating as representing the "exact" result) for both exp2(x) and pow(2.0, x) when the two results differ, for a selection of randomly chosen x in the range(-1000.0, 1000.0). Columns in the output are: x (in hex), x (in decimal), ulps error in exp2(x), ulps error in pow(2.0, x) >>> from decimal import getcontext, Decimal
>>> from math import exp2, pow, ulp
>>> import random
>>> getcontext().prec = 200
>>> def exp2_error_ulps(x):
... libm = exp2(x)
... exactish = 2**Decimal(x)
... return float(Decimal(libm) - exactish) / ulp(libm)
...
>>> def pow2_error_ulps(x):
... libm = pow(2.0, x)
... exactish = 2**Decimal(x)
... return float(Decimal(libm) - exactish) / ulp(libm)
...
>>> for n in range(10000):
... x = random.uniform(-1000.0, 999.0) + random.random()
... if exp2(x) != pow(2.0, x):
... print(f"{x.hex():21} {x:22.17f} {exp2_error_ulps(x): .5f}, {pow2_error_ulps(x): .5f}")
...
0x1.e28f2ad3da122p+5 60.31990590581177969 0.50669, -0.49331
-0x1.929e790e1d293p+9 -805.23806930946227567 0.50082, -0.49918
-0x1.49803564f5b8ap+8 -329.50081473349621319 0.49736, -0.50264
-0x1.534cf08081f4bp+8 -339.30054476902722627 -0.50180, 0.49820
-0x1.b430821fb4ad2p+8 -436.18948553238908517 -0.49883, 0.50117
0x1.2c87a8431bd8fp+8 300.52991122655743084 -0.50376, 0.49624
0x1.3e476f9a09c8cp+7 159.13952332848964488 0.50062, -0.49938
0x1.cb8b9c61e7e89p+9 919.09070991347937252 0.49743, -0.50257
0x1.ab86ed0e6c7f6p+9 855.05410938546879152 0.49742, -0.50258
0x1.97bc9af3cbf85p+9 815.47347876986952997 -0.50076, 0.49924
-0x1.b5434441ba11bp+8 -437.26276026528074681 -0.50062, 0.49938
-0x1.0ead35218910ep+9 -541.35318392937347198 0.50192, -0.49808
-0x1.dbae0b861b89cp+9 -951.35972668022759535 0.50601, -0.49399
0x1.522f005d2dcc4p+6 84.54589982597377684 -0.50704, 0.49296
0x1.398ff48d53ee1p+9 627.12465063665524667 -0.50102, 0.49898
-0x1.381307fbd89f5p+5 -39.00929257159069863 -0.50526, 0.49474
0x1.9dc4c85f7c53ap+9 827.53736489840161994 -0.50444, 0.49556
0x1.b357f6012d3c2p+9 870.68719496449216422 -0.50403, 0.49597
-0x1.a6446703677bbp+9 -844.53439371636284250 0.50072, -0.49928
0x1.e3dd54b28998bp+7 241.93228681497234334 0.49897, -0.50103
0x1.b4f77f18a233ep+8 436.96678308448815642 0.49593, -0.50407
-0x1.578c4ce7a7c1bp+3 -10.73587651486564276 -0.50505, 0.49495
0x1.25a9540e1ee65p+5 36.70767985374258302 0.49867, -0.50133
-0x1.6e220f7db7668p+8 -366.13304887511776542 -0.49904, 0.50096
-0x1.94214ed3e5264p+9 -808.26021813095985635 0.50420, -0.49580
0x1.9dcc3d281da18p+5 51.72472602215219695 -0.50423, 0.49577
-0x1.3ba66909e6a40p+7 -157.82502013149678532 -0.50077, 0.49923
-0x1.9eac2c52a1b47p+9 -829.34510262389892432 -0.50540, 0.49460 |
Sounds good. I've already made the necessary code changes on my own build, so I'll just finish writing the tests + documentation and submit a PR. |
I've submitted a PR at #29829. I'd just like to add that the whole Python team is amazing. Thank you for doing what you do! |
All done. Many thanks, Gideon! |
Bad news: on Windows, exp2(x) is way worse then pow(2, x). Here I changed the loop of Mark's little driver like so: differ = really_bad = 0
worst = 0.0
for n in range(100_000):
x = random.uniform(-1000.0, 999.0) + random.random()
if exp2(x) != pow(2.0, x):
differ += 1
exp2err = exp2_error_ulps(x)
pow2err = pow2_error_ulps(x)
assert abs(pow2err) < 0.52
if abs(exp2err) >= 1.0:
if abs(exp2err) > abs(worst):
worst = exp2err
really_bad += 1
if really_bad < 25:
print(f"{x.hex():21} {x:22.17f} {exp2err:.5f}, {pow2err:.5f}")
print(f"{differ=:,}")
print(f"{really_bad=:,}")
print(f"worst exp2 ulp error {worst:.5f}") Then output from one run: 0x1.0946680d45f28p+9 530.55005041041749791 -1.04399, -0.04399 So they differed in more than a third of the cases; in about a fifth of the differing cases, the exp2 error was at least 1 ulp, and nearly 2 ulp at worst; while in all the differing cases the pow(2, x) error was under 0.52 ulp. |
Across millions of tries, same thing: Windows exp2 is off by at least 1 ulp over a third of the time, and by over 2 ulp about 3 times per million. Still haven't seen pow(2, x) off by as much as 0.52 ulp. From its behavior, it appears Windows implements exp2(x) like so: i = floor(x)
x -= i # now 0 <= x < 1
return ldexp(exp2(x), i) So it's apparently using some sub-state-of-the-art approximation to 2**x over the domain [0, 1]. But a consequence is that it gets it exactly right whenever x is an integer, so it's unlikely anyone will notice it's sloppy ;-) I expect we should just live with it. |
[Tim]
Darn.
Agreed. |
Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.
Show more details
GitHub fields:
bugs.python.org fields:
The text was updated successfully, but these errors were encountered: