10

I'm looking for Python Nth root function/algorithm but before you post: NO INTEGER ROOT, HELL!
Where could I obtain at least a guide how to program Nth root function that produces precise float/Decimal?
Such function that doesn't return 1 nor 0 for root(125, 1756482845) (1st argument is the number, 2nd is the root depth (or something)).

EDIT: So, you were giving me this solution: n ** (1.0 / exp) which I knew when I asked this question, but it just doesn't work for, for example, exp = 3. You can't express 1/3 in terms of rational numbers, so 125 ** (1/3) gives incorrect result 4.999999.... I was asking for some "smart" algorithm, which gives correct result for such nice numbers and at least 4-decimal-points-accurate result for rational exp. If there isn't such function or algorithm, I will use this (n ** (1/exp)).

R.O.S.S
  • 605
  • 5
  • 18
  • 1
    I think either using the `**` operator or `math.pow` with some type casting will do and that's just as much as you can get from builtin functions in Python. I'm not aware of anything that produces more precise floating-point calculation in the built-in modules and functions of the language and would love to know that. If you need to compute more precise result, try looking for computation methods to produce the N-th root. A starting point may be Newton's method. The point is you always have to specify your required threshold for absolute/relative error. – 怀春춘 Sep 30 '16 at 15:08
  • All roots are either integers or irrational numbers. How precise you want the result to be? The decimal (or binary, for that matter) expansion of an irrational number cannot be precisely represented in finite memory. – kfx Sep 30 '16 at 16:43
  • For an interesting read on how challenging this problem is, see http://bugs.python.org/issue27761 – casevh Sep 30 '16 at 20:39
  • Note that your question is a generalized version of http://stackoverflow.com/questions/18063755/computing-a-correctly-rounded-an-almost-correctly-rounded-floating-point-cubic (you ask for any n, it's nontrivial when n is not a power of 2, and the linked question asks about 3). A solution that I did not bother to write up assumes that you have one extra bit of precision available and can compute `pow(candidate, n)` rounded up or down for final adjustments. – Pascal Cuoq Sep 30 '16 at 21:57
  • take a look at [Power by squaring for negative exponents](http://stackoverflow.com/a/30962495/2521214) especially bullet **#5** you can do this on floats/decimals too – Spektre Oct 01 '16 at 07:37

5 Answers5

12

I would try the gmpy2 library.

>>> import gmpy2
>>> gmpy2.root(125,3)
mpfr('5.0')
>>>

gmpy2 uses the MPFR library to perform correctly rounded floating point operations. The default precision is 53 bits but that can be increased.

>>> gmpy2.root(1234567890123456789**11, 11)
mpfr('1.2345678901234568e+18')  # Last digits are incorrect.
>>> gmpy2.get_context().precision=200
>>> gmpy2.root(1234567890123456789**11, 11)
mpfr('1234567890123456789.0',200)
>>>

Disclaimer: I maintain gmpy2.

casevh
  • 11,093
  • 1
  • 24
  • 35
  • Amazing! One of the best libraries out there (if not THE best). I played around with it and it really is precise. Wow. – Stephen Muga Sep 23 '20 at 17:13
3

You can do a binary search on the answer. If you want to find the X that is equal to the kth root of N, you can do a binary search on X testing for each step of the binary search whether X^k equals N +- some small constant to avoid precision issues.

Here is the code:

import math

N,K = map(float,raw_input().split()) # We want Kth root of N
lo = 0.0
hi = N
while 1:
    mid = (lo+hi)/2
    if math.fabs(mid**K-N) < 1e-9: # mid^K is really close to N, consider mid^K == N
        print mid
        break
    elif mid**K < N: lo = mid
    else: hi = mid

For (N,K) = (125,3) it prints 5.0, the correct answer. You can make it more precise by changing the 1e-9 constant, but there is a precision limit related to the float variables precision limit in Python

Lucas Sampaio
  • 306
  • 1
  • 7
  • 1
    `return True if math.fabs(a-b) < 1e-9 else False` is a rather verbose way to say `return math.fabs(a-b) < 1e-9`; that'll return `True` and `False` just fine. Also note, as of Python 3.5, there is a more powerful/flexible version of that test available, [`math.isclose`](https://docs.python.org/3/library/math.html#math.isclose). – ShadowRanger Oct 01 '16 at 02:22
  • Thank you for the tip! Fixing it – Lucas Sampaio Oct 01 '16 at 02:39
  • If the `125 ** (1/3)` solution to compute the cubic root is not acceptable to the OP because of the approximation in the computation of 1/3, there is no reason that this answer, based on a hardcoded epsilon, would satisfy them. You need a correctly rounded `**` in either round-upwards or round-downwards at one extra bit of precision to make this work, and not use any epsilon, or you are just making matters worse on average. – Pascal Cuoq Oct 02 '16 at 09:22
1

You mean something like that:

>>> 125**(1/9.0)
1.7099759466766968

Something else that might interest you is bigfloat module (haven't used personally just know it exists :) - actually had problem installing it in the past-maybe an OS X fault)

coder
  • 12,832
  • 5
  • 39
  • 53
  • But wait, 1/9.0 isn't even one ninth to begin with, so how precise is this? – harold Sep 30 '16 at 15:00
  • @harold, just said an idea alternative to `math.pow`, now about how precise this is, this has to do with what you are looking for... I don't know if the OP is ok with the result, but don't forget that in math precision is something really relevant. Anyway do you have anything else to propose? – coder Sep 30 '16 at 15:03
  • @R.O.S.S, check out `bigfloat` module - maybe that does what you need - have updated the answer with the link – coder Sep 30 '16 at 15:18
1

It's the function pow of math module.

import math
math.pow(4, 0.5)

will return the square root of 4, which is 2.0.

For root(125, 1756482845), what you need to do is

math.pow(125, 1.0 / 1756482845)
怀春춘
  • 197
  • 5
  • 1
    This isn't precise, see http://stackoverflow.com/questions/15978781/how-to-find-integer-nth-roots - `5 ** 3 = 125`, but `125 ** (1.0 / 3) = 4`. – R.O.S.S Sep 30 '16 at 14:57
  • 1
    I try math.pow(125, 1.0 / 1756482845) ** 125 which returns 125.00000486047219 on my machine. It depends on what is your accepted level of relative/absolute error. For absolute error of 1e-5 the code words, doesn't it? Maybe you should elaborate your question regarding how precise you want the calculation to be. – 怀春춘 Sep 30 '16 at 15:03
  • `math.pow` is pointless; the `**` operator accomplishes the same result, faster, and even if you wanted a function call, there is already a built-in (no import required) `pow` function. About the only thing `math.pow` does is force a `float` result even if the base and exponent are integer; losing precision unnecessarily isn't much of a selling point. – ShadowRanger Oct 01 '16 at 02:24
1

In Squeak Smalltalk, there is a nthRoot: message that answers the exact Integer result if ever the Integer receiver is exact nth power of some whole number. However, if the solution is an algebraic root, then the implementation does not fallback to a naive n**(1/exp); the method rounds to nearest float by appropriate care of residual.

Relevant code (MIT license) is reproduced here. The base algorithm is searching for truncated nth root of an Integer with some Newton-Raphson:

Integer>>nthRootTruncated: aPositiveInteger
    "Answer the integer part of the nth root of the receiver."
    | guess guessToTheNthMinusOne nextGuess |
    self = 0 ifTrue: [^0].
    self negative
        ifTrue:
            [aPositiveInteger even ifTrue: [ ArithmeticError signal: 'Negative numbers don''t have even roots.' ].
            ^(self negated nthRootTruncated: aPositiveInteger) negated].
    guess := 1 bitShift: self highBitOfMagnitude + aPositiveInteger - 1 // aPositiveInteger.
    [
        guessToTheNthMinusOne := guess raisedTo: aPositiveInteger - 1.
        nextGuess := (aPositiveInteger - 1 * guess * guessToTheNthMinusOne + self) // (guessToTheNthMinusOne * aPositiveInteger).
        nextGuess >= guess ] whileFalse:
            [ guess := nextGuess ].
    ( guess raisedTo: aPositiveInteger) > self  ifTrue:
            [ guess := guess - 1 ].
    ^guess

It's not particularly clever, because convergence can be very slow in case of huge exponent, but well, it works. Then, the same root rounded away from zero:

Integer>>nthRootRounded: aPositiveInteger
    "Answer the integer nearest the nth root of the receiver."
    | guess |
    self = 0 ifTrue: [^0].
    self negative
        ifTrue:
            [aPositiveInteger even ifTrue: [ ArithmeticError signal: 'Negative numbers don''t have even roots.' ].
            ^(self negated nthRootRounded: aPositiveInteger) negated].
    guess := self nthRootTruncated: aPositiveInteger.
    ^self * 2 > ((guess + 1 raisedTo: aPositiveInteger) + (guess raisedTo: aPositiveInteger))
        ifTrue: [guess + 1]
        ifFalse: [guess]

Then exactness is tested in nthRoot:

Integer>>nthRoot: aPositiveInteger
    "Answer the nth root of the receiver.
    Answer an Integer if root is exactly this Integer, else answer the Float nearest the exact root."

    | guess excess scaled nBits |
    guess := self nthRootRounded: aPositiveInteger.
    excess := (guess raisedTo: aPositiveInteger) - self.
    excess = 0 ifTrue: [ ^ guess ].

    nBits := Float precision - guess highBitOfMagnitude.
    nBits <= 0 ifTrue: [ ^(Fraction numerator: guess * 4 - excess sign denominator: 4) asFloat].

    scaled := self << (nBits * aPositiveInteger).
    guess := scaled nthRootRounded: aPositiveInteger.
    excess := (guess raisedTo: aPositiveInteger) - scaled.
    ^(Fraction numerator: guess * 4 - excess sign denominator: 1 << (nBits + 2)) asFloat

This could be applied to Fraction too, but the nearest float is a bit more complex, and Squeak implementation is currently naive.

It works for large integers like:

  • (10 raisedTo: 600) nthRoot: 300 -> 100 "exact"
  • (10 raisedTo: 600) + 1 nthRoot: 300 -> 100.0 "inexact"

If you don't have such expectations, the initial guess could use inexact naive n**(1/exp).

The code should be easy to port in Python and leaves a lot of place for optimization.

I didn't check what was available in Python, but maybe you'll need correctly rounded LargeInteger -> Float, and Fraction -> Float, like explained here (Smalltalk too, sorry about that, but the language does not really matter).

aka.nice
  • 9,100
  • 1
  • 28
  • 40