Tim Peters added the comment:

Let me clarify something about the extended algorithm:  the starting guess is 
no longer the most significant source of error.  It's the `mul(D(x), pow(D2, 
e))` part.  `D(x)` is exact, but `pow(D2, e)` may not be exactly representable 
with 26 decimal digits, and the multiplication may also need to round.  Those 2 
operations may each suffer a half ulp error (but "ulp" in this case is 1 in the 
26th decimal digit).

The rough error analysis is straightforward to extend, but I used wxMaxima to 
do the Taylor expansions and to ensure I didn't drop any terms.  

So suppose that instead of exactly x*2**e we get x*2**e*(1+r1) for some small 
|r1| (on the order of 10**-25), and instead of the exact n'th root `t` the 
starting guess g = t*(1+r2) for some small |r2| (on the order of 2**-50).  Then 
the mathematical value of the corrected `g` is `t` times (dropping terms higher 
than quadratic in `r1` and `r2`):

1   + r1/n   - (n-1)/n * r1*r2   + (n-1)/2 * r2**2

The largest term (besides 1) by far is now `r1/n` (assuming a 
not-insanely-large `n`).  But that's still "way at the end" of the 
extended-precision decimal format, which _should be_ ;-) intuitively obvious.  
It's too small to have a direct effect on the last bit of the 53-bit result 
(but can certainly affect "close call" rounding cases).


Python tracker <rep...@bugs.python.org>
Python-bugs-list mailing list

Reply via email to