[issue27761] Private _nth_root function loses accuracy

[Tim Peters]

Tim Peters added the comment:

Let's spell one of these out, to better understand why sticking to native precision is inadequate.  Here's one way to write the Newton step in "guess + relatively_small_correction" form:

    def plain(x, n):
        g = x**(1.0/n)
        return g - (g - x/g**(n-1))/n

Alas, it's pretty much useless.  That's because when g is a really good guess, `g` and `x/g**(n-1)` are, in native precision, almost the same.  So the subtraction cancels out almost all the bits, leaving only a bit or two of actual information.  For this kind of approach to be helpful in native precision, it generally requires a clever way of rewriting the correction computation that _doesn't_ suffer massive cancellation.

Example:

>>> x = 5.283415219603294e-14

and we want the square root.  Mark's `nroot` always gives the best possible result:

>>> nroot(x, 2)
2.298568080262861e-07

In this case, so does sqrt(), and also x**0.5 (on my box - pow() may not on yours):

>>> sqrt(x)
2.298568080262861e-07
>>> x**0.5
2.298568080262861e-07

You're going to think it needs "improvement", because the square of that is not x:

>>> sqrt(x)**2 < x
True

Let's see what happens in the `plain()` function above:

>>> g = x**0.5
>>> temp = x/g**(2-1)

>>> g.hex()
'0x1.ed9d1dd7ce57fp-23'
>>> temp.hex()
'0x1.ed9d1dd7ce580p-23'

They differ by just 1 in the very last bit.  There's nothing wrong with "the math", but _in native precision_ the subtraction leaves only 1 bit of information.

Then:

>>> correction = (g - temp)/2
>>> correction
-1.3234889800848443e-23

is indeed small compared to g, but it only has one "real" bit:

>>> correction.hex()
'-0x1.0000000000000p-76'

Finally,

>>> g - correction
2.2985680802628612e-07

is _worse_ than the guess we started with.  Not because of "the math", but because sticking to native precision leaves us with an extremely crude approximation to the truth.

This can be repaired in a straightforward way by computing the crucial subtraction with greater than native precision.  For example,

    import decimal
    c = decimal.DefaultContext.copy()
    c.prec = 25

    def blarg(x, n,
              D=decimal.Decimal,
              pow=c.power,
              sub=c.subtract,
              div=c.divide):
        g = x**(1.0/n)
        Dg = D(g)
        return g - float(sub(Dg, div(D(x), pow(Dg, n-1)))) / n

    del decimal, c

Then the difference is computed as

Decimal('-2.324439989147273024835272E-23')

and the correction (convert that to float and divide by 2) as:

-1.1622199945736365e-23

The magnitude of that is smaller than of the -1.3234889800848443e-23 we got in native precision, and adding the new correction to g makes no difference:  the correction is (correctly!) viewed as being too small (compared with g) to matter.

So, bottom line:  there is no known way of writing the Newton step that won't make things worse at times, so long as it sticks to native precision.  Newton iterations in native precision are wonderful when, e.g., you know you have about 20 good bits and want to get about 40 good bits quickly.  The last bit or two remain pure noise, unless you can write the correction computation in a way that retains "a bunch" of significant bits.  By far the easiest way to do the latter is simply to use more precision.

----------

_______________________________________
Python tracker <report at bugs.python.org>
<https://bugs.python.org/issue27761>
_______________________________________