[issue27761] Private _nth_root function loses accuracy

Tim Peters report at bugs.python.org
Wed Aug 17 01:17:28 EDT 2016


Tim Peters added the comment:

Noting that `floor_nroot` can be sped a lot by giving it a better starting guess.  In the context of `nroot`, the latter _could_ pass `int(x**(1/n))` as an excellent starting guess.  In the absence of any help, this version figures that out on its own; an optional `a=None` argument (to supply a starting guess, if desired) would make sense.

    def floor_nroot(x, n):
        """For positive integers x, n, return the floor of the nth root of x."""

        bl = x.bit_length()
        if bl <= 1000:
            a = int(float(x) ** (1.0/n))
        else:
            xhi = x >> (bl - 53)
            # x ~= xhi * 2**(bl-53)
            # x**(1/n) ~= xhi**(1/n) * 2**((bl-53)/n)
            # Let(bl-53)/n = w+f where w is the integer part.
            # 2**(w+f) is then 2**w * 2**f, where 2**w is a shift.
            a = xhi ** (1.0 / n)
            t = (bl - 53) / n
            w = int(t)
            a *= 2.0 ** (t - w)
            m, e = frexp(a)
            a = int(m * 2**53)
            e += w - 53
            if e >= 0:
                a <<= e
            else:
                a >>= -e
        # A guess of 1 can be horribly slow, since then the next
        # guess is approximately x/n.  So force the first guess to
        # be at least 2.  If that's too large, fine, it will be
        # cut down to 1 right away.
        a = max(a, 2)
        a = ((n-1)*a + x // a**(n-1)) // n
        while True:
            d = x // a**(n-1)
            if a <= d:
                return a
            a = ((n-1) * a + d) // n

I haven't yet found a case in the context of `nroot` where it doesn't get out on the first `if a <= d:` test.  Of course you can provoke as many iterations as you like by passing `x` with a lot more than 53 "significant" bits (the large `x` passed by `nroot` are mostly long strings of trailing 0 bits).

----------

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


More information about the Python-bugs-list mailing list