On Sun, Aug 14, 2016 at 10:52:42AM +0000, Mark Dickinson wrote:
> Same deal here: those aren't the actual errors; they're approximations 
> to the errors, since the computations of the epsilons depends on (a) 
> the usual floating-point rounding, and more significantly (b) the 
> accuracy of the `y**n` computation. It's entirely possible that the 
> value giving the smaller epsilon is actually the worse of the two 
> approximations.

Let's call the two calculated roots y and w, where y is generated by 
pow() and w is the "improved" version, and the actual true 
mathematical root is r (which might fall between floats).

You're saying that it might turn out that abs(y**n - x) < abs(w**n - x), 
but abs(w - r) < abs(y - r), so I return the wrong calculated root.

I can see how that could happen, but not what to do about it. I'm 
tempted to say "that platform's pow() is weird, and the best I can do 
is return whatever it returns". That way you're no worse than if I just 
used pow() and didn't try to improve the result at all.

I think this would be a problem if I wanted to make nth_root() a public 
function with a guarantee that it will be better than pow(), but at 
worst here it just slows the geometric mean down a bit compared to a 
naive pow(product, 1/n).

What else can I do? Since I'm only dealing with integer powers, should I 
try using my own ipow(y, n) for testing?

def ipow(y, n):
    x = 1.0
    while n > 0:
        if n % 2 == 1:
            x *= y
        n //= 2
        y *= y
    return x