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
Javascript: https://github.com/stdlib-js/math-base-special-exp2
Numpy: https://numpy.org/doc/stable/reference/generated/numpy.exp2.html
C++: https://en.cppreference.com/w/cpp/numeric/math/exp2 (Not authoritative)
Ruby: https://www.rubydoc.info/gems/ruby-mpfi/MPFI%2FMath.exp2

Similar Issues:
https://bugs.python.org/issue44357
https://bugs.python.org/issue31980

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.)

> a libm exp2 is supposedly more accurate than pow(2.0, x), though I don’t really see how this would be the case

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 https://github.com/python/cpython/pull/29829.

I'd just like to add that the whole Python team is amazing. Thank you for doing what you do!

New changeset 6266e4af873a27c9d352115f2f7a1ad0885fc031 by Gideon in branch 'main':
bpo-45917: Add math.exp2() method - return 2 raised to the power of x (GH-29829)
https://github.com/python/cpython/commit/6266e4af873a27c9d352115f2f7a1ad0885fc031

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
0x1.de4f9662d84f8p+9   956.62177691995657369 -1.00976, -0.00976
-0x1.60f9152be0a09p+4  -22.06081120624188330 1.02472, 0.02472
-0x1.687b056d7a81ap+8 -360.48055156937482479 1.48743, 0.48743
0x1.8e97e9d405622p+9   797.18682337057930454 1.05224, 0.05224
-0x1.2d1e3a03eda7fp+9 -602.23614548782632028 -1.21876, -0.21876
0x1.3af55e79cd45dp+8   314.95847283612766887 -1.10044, -0.10044
0x1.0e7fba610cde6p+9   540.99787533882476964 -1.39782, -0.39782
0x1.9c7d0258e460dp+9   824.97663413192060489 1.19690, 0.19690
0x1.3de5064eb1598p+9   635.78925498637818237 1.75376, -0.24624
-0x1.d5189d23da3d0p+9 -938.19229553371587826 1.07734, 0.07734
0x1.967d0857aa500p+9   812.97681709114112891 1.23630, 0.23630
-0x1.30ee89e018914p+6  -76.23294782781550794 -1.10275, -0.10275
-0x1.e35eb8936dddbp+9 -966.74000780930089149 -1.02686, -0.02686
-0x1.28d40d7693088p+6  -74.20708260795993283 1.00352, 0.00352
-0x1.e965d067d1084p+7 -244.69885563303625986 1.21136, 0.21136
-0x1.b1fbeec1c1ba3p+7 -216.99205594529948371 -1.05536, -0.05536
-0x1.543d715a5824cp+9 -680.48002175620922571 1.24955, 0.24955
0x1.526829d46c034p+9   676.81377654336984051 -1.17826, -0.17826
-0x1.bdaf1d7850c74p+6 -111.42101085656196346 1.08670, 0.08670
-0x1.48218d1605dd0p+9 -656.26211810385029821 1.06103, 0.06103
-0x1.16298bcdb9103p+9 -556.32457896744051595 -1.23732, -0.23732
0x1.39ff24b1a7573p+8   313.99665365539038930 -1.20931, -0.20931
0x1.9cdf1d0101646p+8   412.87153631481157845 -1.23481, -0.23481
differ=38,452
really_bad=7,306
worst exp2 ulp error -1.91748

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]

> on Windows, exp2(x) is way worse then pow(2, x)

Darn.

> I expect we should just live with it.

Agreed.