Why doesn't this Taylor series for exponentiation work?

Question

Don't know if this is more mathematics or programming. Anyway I want to implement an efficient way to approximate exponentiation involving non-integral exponents, faster than C++'s pow from cmath.

I just found this series online:

$a^x = 1 + x \ln a + \frac{(x \ln a) ^ 2} {2} + \frac{(x \ln a) ^ 3} {6} + \frac{(x \ln a) ^ 4} {24}...$

So I implemented it in Python:

import numba as nb
import numpy as np

f = 1
factorials = np.array([(f := f * i) for i in range(1, 17)], dtype=np.int64)

@nb.njit(cache=True, fastmath=True)
def taylor(base: float, exp: float, lim: int) -> float:
    c = exp * log(base)
    f = 1
    s = 1
    for i in range(lim):
        f *= c
        s += f / factorials[i]
    
    return s

And it seems to work well for exponents between 0 and 1, but doesn't work at all for exponents greater than 1:

In [397]: taylor(256, 1/3, 16)
Out[397]: 6.34960420776535

In [398]: 256**(1/3)
Out[398]: 6.3496042078727974

In [399]: taylor(256, 3, 16)
Out[399]: 8441719.024421709

In [400]: 256**3
Out[400]: 16777216

In [401]: taylor(256, 2, 16)
Out[401]: 61649.26195536944

In [402]: 256**2
Out[402]: 65536

In [403]: taylor(256, 1, 16)
Out[403]: 255.9821624637957

In [404]: taylor(256, 0.5, 16)
Out[404]: 15.999999887823558

I implemented the same logic as the series shown above, so the series must be the cause.

Was I lied to? Is the series wrong?

Answer 1

The given Taylor series is correct. However, when you truncate it at $n$ terms, the remaining error will be on a similar scale to the next term in the series - for example, if you calculate the value using three terms, i.e.

$$a^x \approx 1 + x \ln a + \frac{(x \ln a)^2}{2}$$

then the error will be on the order of $\frac{(x \ln a)^3}{6}$. This means that as $x$ gets larger, the error is going to also increase. And it also grows in relative terms - for any $a > 1$, the exponential function $a^x$ always eventually grows faster than any polynomial in $x$ as $x$ goes to infinity. If you're calculating 20 terms rather than 3 then it will take longer to reach a significant divergence, but there will always be a point where you lose control of the error.

One method to avoid this problem is roughly:

1. Shift and/or scale $x$ to put it in a particular range, e.g. find an integer $n$ such that $x_0 = \frac{x}{2^n} \in [1, 2]$.

2. Calculate $a^{x_0}$ to some desired accurately with a suitable method (as mentioned in the comments, there are often better options than Taylor approximations).

3. Calculate $(a^{x_0})^{2^n}$ via repeated squaring.

Even in this case, the squaring step is going to cause the error from step 2 to potentially accumulate. This is why it's common to calculate exponents via other methods, such as via a log-transformation.

Answer 2

The series is not wrong, you can obtain it from the exponential series. In fact, recall that $a^x=e^{x\log a}$, substituting in the exponential series you have: $$a^x=e^{x\log a}=\sum_{k=0}^{+\infty}\frac{(x\log a)^k}{k!}=1+x\log a+\frac{(x\log a)^2}{2}+\dots $$