What is the numerical stability of std::pow() compared to iterated multiplication?

Question

What sort of stability issues arise or are resolved by using std::pow()?

• Will it be more stable (or faster, or at all different) in general to implement a simple function to perform log(n) iterated multiplies if the exponent is known to be an integer?

• How does std::sqrt(x) compare, stability-wise, to something of the form std::pow(x, k/2)? Would it make sense to choose the method preferred for the above to raise to an integer power, then multiply in a square root, or should I assume that std::pow() is fast and accurate to machine precision for this? If k = 1, is there a difference from std::sqrt()?

• How would std::pow(x, k/2) or the method above compare, stability-wise, to an integer exponentiation of std::sqrt(x)?

And as a bonus, what are the speed differences likely to be?

Festus

You're asking a question of which parts are amenable to experimentation and analysis. What are the results of your experiments and analysis ? – Pimiento

Pascal Cuoq's answer is excellent, but I would add that lots of traditional and popular implementations of pow are outrageously crappy. – Beaverette

@tmyklebu, which ones? – Festus

@Festus I almost mentioned this but it just complicated the message. There are a number of confused questions on this site caused by pow(10.0, 2.0) evaluating to a double other than 100.0 on Windows (example: https://mcmap.net/q/493614/-why-is-my-integer-math-with-std-pow-giving-the-wrong-answer ). Since the mathematical result is exactly representable as the double 100.0, it means that that pow function commits an error of more than 1ULP. But even if you use a botched pow function that is sometimes wrong by, say, 1.5 ULP, it remains more accurate than four multiplications. – Meyers

@trbabb: I haven't catalogued them, but I've seen a fairly recent Linux machine's pow be garbage. – Beaverette

Answer 1

Will it be more stable (or faster, or at all different) in general to implement a simple function to perform log(n) iterated multiplies if the exponent is known to be an integer?

The result of exponentiation by squaring for integer exponents is in general less accurate than pow, but both are stable in the sense that close inputs produce close results. You can expect exponentiation by squaring to introduce 0.5 ULP of relative error by multiplication (for instance, 1 ULP of error for computing x³ as x * x * x).

When the second argument n is statically known to be 2, then by all means implement xⁿ as x * x. In that case it is faster and more accurate than any possible alternative.

How does std::sqrt(x) compare, stability-wise, to something of the form std::pow(x, k/2)

First, the accuracy of sqrt cannot be beat for an IEEE 754 implementation, because sqrt is one of the basic operations that this standard mandates to be as accurate as possible.

But you are not asking about sqrt, you are asking (I think) about <computation of xⁿ> * sqrt(x) as opposed to pow(x, n + 0.5). Again, in general, for a quality implementation of pow, you can expect pow(x, n + 0.5) to be more accurate than the alternatives. Although sqrt(x) would be computed to 0.5 ULP, the multiplication introduces its own approximation of up to 0.5 ULP, and all in all, it is better to obtain the result you are interested in in a single call to a well-implemented function. A quality implementation of pow will give you 1 ULP of accuracy for its result, and the best implementations will “guarantee” 0.5 ULP.

And as a bonus, what are the speed differences likely to be?

If you know in advance that the exponent is going to be a small integer or multiple of 0.5, then you have information that the implementer of pow did not have, so you can beat them by at least the cost of the test to determine that the second argument is a small integer. Plus, the implementer of a quality implementation is aiming for a more accurate result than simple exponentiation by squaring provides. On the other hand, the implementer of pow can use extremely sophisticated techniques to minimize the average execution time despite the better accuracy: see for instance CRlibm's implementation. I put the verb “guarantee” above inside quotes when talking about the best implementations of pow because pow is one function for which CRlibm's 0.5 ULP accuracy guarantee is only “with astronomical probability”.

Meyers

link to CRlibm seems broken or site down? – Berkley