Calculation of maximum floating point error

Question

pow(x,y) is computed as e^(y*log(x)) Generally math libraries compute log(x) in quad precision (which is time consuming) in order to avoid loss of precision when computing y*log(x) as precision errors will magnify in the computation of e^(y*log(x)). Now, in case I would want to compute pow(x,y) in the following steps.

double pow(double x,double y) {
   return exp(y*log(x)); // Without any quad precision multiplication
}

What would be the maximum ULP error of this function. I do know that IEEE-754 standard says that any floating point operation should have less than 0.5 ULP error i.e 0.5*2^(-52). So if my operation y*log(x) suffers from a 0.5 ULP error, how do I compute the largest possible ULP error for e^(y*log(x))

Agreed that the computation of pow(x,y) is fairly complicated. Algorithms generally compute log(x) in a higher precision and the multiplication between y and log(x) is not straightforward. Since the ULP error depends on y*log(x), the maximum error would be for the largest value of Y*log(x) for which e^(y*log(x)) is not infinity. Right? How do I compute the number of ULP for such a case? What are the maximum number of bits of the mantissa in the double precision format that would vary form the actual value in case of the largest value of y*log(x) ?

Updated the question. Thanks for all the help!

So a 10 bit difference would result in how much ULP error? I calculated it as,

    ULP = (actual - computed)/ 2^(e-(p-1))

where e is the exponent of the actual number, p=53 for double precision. I read that I ULP = 2^(e-(p-1)) Let's assume,

    Actual = 1.79282279439444787915898270592 *10^308
    Computed = 1.79282279439451553814547593293 * 10^308 
    error= actual - computed = 6.7659e+294

Now

1 ULP = 2^(e- (p-1)) 
e = 1023 (exponent of the actual number - bias)
1 ULP = 2^(1023 - (53-1)) = 2^971
ULP error = error/ 1 ULP = 339

Is this correct?

Exponential functions convert an absolute error into a percentage error, so the ULP error in pow depends on the value of y log x. — meowgoesthedog
The statements in the question about pow are not generally true. Implementations of pow are fairly complicated and typically involve extra precision (to avoid loss of accuracy, not loss of precision), but the extra precision is often obtained by careful work with multiple floating-point or integer objects, not by “quad precision.” Additionally, the computation is not as simple as e^(y*•log(*x)). Aside from partitioning the problem into various cases, the exponent and significand of x are often separated. — Eric Postpischil
Note that the errors in y*log(x) are not limited to ½ ULP, as it consists of two operations—taking a logarithm and performing a multiplication. — Eric Postpischil
The limit on the error in IEEE-754 elementary operations in round-to-nearest-ties-to-even mode is less than or equal to ½ ULP, not less than ½ ULP. The latter is impossible, as a mathematical result that is exactly midway between two representable values cannot be rounded with less than ½ ULP error. — Eric Postpischil
Agreed that the computation of pow(x,y) is fairly complicated. Algorithms generally compute log(x) in a higher precision and the multiplication between y and log(x) is not straightforward. — Joseph Arnold

Eric Postpischil Eric Postpischil · Accepted Answer · 2018-10-01T17:22:55

I do not have time for a fuller answer now, but here is something to start:

Suppose the error in each individual operation—exp, *, and log—is at most ½ ULP. This means the computed value of log(x) is exactly log(x)•(1+e₀) for some e₀ satisfying -𝜖/2 ≤ e₀ ≤ 𝜖/2, where 𝜖 is ULP of 1 (2⁻⁵² IEEE-754 basic 64-bit binary format). And the computed value of y*log(x) is exactly y•log(x)•(1+e₀)•(1+e₁) for some -𝜖/2 ≤ e₀ ≤ 𝜖/2 and some -𝜖/2 ≤ e₁ ≤ 𝜖/2. And the computed value of exp(y*log(x)) is exactly e^{y•log(x)•(1+e₀)•(1+e₁)}•(1+e₂) for some -𝜖/2 ≤ e₀, e₁, e₂ ≤ 𝜖/2. And the error is of course e^{y•log(x)•(1+e₀)•(1+e₁)}•(1+e₂) − e^y•log(x).

Then you can start analyzing that expression for maximum and minimum values over the possible values of e₀, e₁, and e₂. If 0 < y and 1 < x, the greatest error occurs when e₀ = e₁ = e₂ = 𝜖/2.

Note that common implementations of exp and log are not usually correctly rounded. Errors of several ULP may be common. So you should not assume a ½ ULP bound unless the routine is documented to be correctly rounded.

Calculation of maximum floating point error

2 Answers