Precision Tail-off?

[Oscar Benjamin]

On Sat, 18 Feb 2023 at 11:19, Peter J. Holzer <hjp-python at hjp.at> wrote:
>
> On 2023-02-18 03:52:51 +0000, Oscar Benjamin wrote:
> > On Sat, 18 Feb 2023 at 01:47, Chris Angelico <rosuav at gmail.com> wrote:
> > > On Sat, 18 Feb 2023 at 12:41, Greg Ewing via Python-list
> > > > To avoid it you would need to use an algorithm that computes nth
> > > > roots directly rather than raising to the power 1/n.
> > > >
> > >
> > > It's somewhat curious that we don't really have that. We have many
> > > other inverse operations - addition and subtraction (not just "negate
> > > and add"), multiplication and division, log and exp - but we have
> > > exponentiation without an arbitrary-root operation. For square roots,
> > > that's not a problem, since we can precisely express the concept
> > > "raise to the 0.5th power", but for anything else, we have to raise to
> > > a fractional power that might be imprecise.
> >
> > Various libraries can do this. Both SymPy and NumPy have cbrt for cube roots:
>
> Yes, but that's a special case. Chris was talking about arbitrary
> (integer) roots. My calculator has a button labelled [x√y], but my
> processor doesn't have an equivalent operation.

All three of SymPy, mpmath and gmpy2 can do this as accurately as
desired for any integer root:

  >>> n = 123456789000000000000000000000000000000000000000000000000000000

  >>> sympy.root(n, 6)
  1000000000*13717421**(1/6)*3**(1/3)
  >>> sympy.root(n, 6).evalf(50)
  22314431635.562095902499928269233656421704825692573

  >>> mpmath.root(n, 6)
  mpf('22314431635.562096')
  >>> mpmath.mp.dps = 50
  >>> mpmath.root(n, 6)
  mpf('22314431635.562095902499928269233656421704825692572746')

  >>> gmpy2.root(n, 6)
  mpfr('22314431635.562096')
  >>> gmpy2.get_context().precision = 100
  >>> gmpy2.root(n, 6)
  mpfr('22314431635.56209590249992826924',100)

There are also specific integer only root routines like
sympy.integer_nthroot or gmpy2.iroot.

  >>> gmpy2.iroot(n, 6)
  (mpz(22314431635), False)
  >>> sympy.integer_nthroot(n, 6)
  (22314431635, False)

Other libraries like the stdlib math module and numpy define some
specific examples like cbrt or isqrt but not a full root or iroot.
What is lacking is a plain 64-bit floating point routine like:

  def root(x: float, n: int) -> float:
      return x ** (1/n) # except more accurate than this

It could be a good candidate for numpy and/or the math module. I just
noticed from the docs that the math module has a new in 3.11 cbrt
function that I didn't know about which suggests that a root function
might also be considered a reasonable addition in future. Similarly
isqrt was new in 3.8 and it is not a big leap from there to see
someone adding iroot.

--
Oscar