BUG? DOC? Why is exp(p*log(x)) faster than power(x, p)?

Describe the issue:

I assume this is a documentation issue rather than a real bug.

power() is very slow.
I suppose it must be slow in order to achieve accuracy?
Maybe the documentation can explain how much accuracy is lost by using the faster approach exp(p*log(x)), so that users can decide between the two approaches.

Reproduce the code example:

import numpy as np

def power(x: np.ndarray, p: float):
    temp = np.empty_like(x)
    np.log(x, out=temp)
    np.multiply(p, temp, out=temp)
    np.exp(temp, out=temp)
    return temp

x = np.logspace(-4, 4, 50_000_000)

%%timeit
power(x, 2.1)
"477 ms ± 34.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)"

%%timeit
np.power(x, 2.1)
"2.71 s ± 32.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)"

np.allclose(power(x, 2.1), np.power(x, 2.1))
"True"

Error message:

No response

NumPy/Python version information:

1.22.4 3.9.7 (default, Sep 16 2021, 13:09:58)
[GCC 7.5.0]

Context for the issue:

No response

Start with the reproduction code comparing the custom exp(log()) approach with np.power(x, 2.1), then investigate the accuracy difference across the shown input range. Done means documenting why np.power is slower and quantifying the accuracy tradeoff so users can choose between the approaches.