A quick precision test

[Claudio L.]

(06-04-2014 10:09 PM)Paul Dale Wrote:  

(06-04-2014 02:17 PM)Claudio L. Wrote:  The power function is guaranteed to give you an accurate result with 34 good digits, for any 2 parameters (in this case, your guess that it's less precise is incorrect).

I very very much doubt this. Perhaps results with one ULP error but even then I doubt it is that good. The power function is a very difficult one to get right. The GNU C library implementation for 64 bit doubles, in bad cases, falls back to multiple precision calculations carrying over 1600 guard bits. I expect someone rather clever spent a lot of time and effort determining that number as an upper bound to guarantee correctness for the 53 mantissa bits in the result.

I don't know that much, but the documentation of mpdecimal says the error is 1 ULP, so we use 9 extra digits and then we round, so the 'bad rounding' case would have to be in the form 999999995 in order for the 1 ULP error to give a badly rounded result. Notice that I said 'accurate digits' not 'correctly rounded'. The 1 ULP case may still happen, with a chance of 1 in 10^9.

(06-04-2014 10:09 PM)Paul Dale Wrote:  

Quote:I don't know how decNumber does it in particular, but in general powers with small integer exponents are handled by a sequence of multiplications, while powers of real numbers or large numbers are usually done through logarithms:

Code:

x^a = exp( a * ln(x) )

Which is usually very accurate, as ln() and exp() are both guaranteed to give you 34+GUARD good digits on any number you throw at them.

Two points here:

1. Expressing power as x^a = exp( a * ln(x) ) is prone to digit loss. Don't expect an accurate result unless you have a lot of extra digits. See the 1600 above
2. The decNumber versions I've looked at do not guarantee ln() and exp() are correctly rounded, the source just claims close but usually within one ULP.
   

A good rule of thumb to follow is to assume correct rounding for the four arithmetic operations and square root -- these are mandated by the IEEE floating point standards these days. Everything else is close.

- Pauli

I agree, but "close" as in the last digit may be off by 1 due to rounding and only in some corner cases (like 1 in 10^9 chance), on most numbers you should get the right value.

I'd love to test some corner cases to see how the algorithms behave (like that case that needs 1600 digits, do you have any links to good material on this subject?). I did some tests of my own but I don't know what the corner cases are. I'd like to have some literature where they list cases that I can use to stress the algorithms and see the results.

Claudio