A quick precision test

[Claudio L.]

(06-05-2014 12:06 PM)Paul Dale Wrote:  There are other rounding cases. What about nnnnm|500000000|1 when m is odd and you are rounding to even? Or nnnnn|000000000|1 when you are rounding towards +infinity? Assuming you support rounding modes of course.

Yes, there's one case for each rounding mode, so it's still 1 in 10^9 chance no matter which mode you select.

(06-05-2014 12:06 PM)Paul Dale Wrote:  Agreed, the chance of being incorrect - if your internal result is correct - is low but it is still there. Correct rounding is a worthwhile, if somewhat quixotic, goal. One ULP is far more achievable.

Are you sure all of your guard digits are correct or even all the digits in your result? In general they (the guard digits) won't all be. There are many pitfalls for the unwary here. Even changing the order of operations can cause problems. You need to be amazingly careful to get all digits correct, if indeed this is even possible (which it isn't in general).

No, I cannot guarantee that all my guard digits are all correct, and that's one more reason not to obsess about rounding.
As long as the most significant of my guard digits is good, there's only 0.5 ULP worst-case error (if I round the last digit up or down incorrectly).

On the other hand, when I implemented transcendentals, I thought it this way:
* The user should have 2000 good digits, that was my target. Since the library works in groups of 9 digits, you actually get 2007 (multiple of nine).
* For me to be able to get guard digits, I calculated all the CORDIC constants with 2016 digits (9 digits more).
* In the CORDIC loop, I set the precision to 9 digits more (2025) for all intermediate operations. It may seem dumb, because the constants have only 2016, but it actually helps get me a few more good digits.

With this, I can get transcendentals correct to the last digit that will be displayed to the user on all the test cases I tried (except the ones where you lose precision, like your COS(near PI) case, but that comes from the argument, not the algorithm).

Claudio