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

<p><strong>EDIT:</strong> So, you were giving me this solution: <code>n ** (1.0 / exp)</code> which I knew when I asked this question, but it just doesn't work for, for example, <code>exp = 3</code>. You can't express <code>1/3</code> in terms of rational numbers, so <code>125 ** (1/3)</code> gives incorrect result <code>4.999999...</code>. 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 <code>exp</code>. If there isn't such function or algorithm, I will use this (<code>n ** (1/exp)</code>).</p>

<p>In Squeak Smalltalk, there is a <code>nthRoot:</code> message that answers the exact <code>Integer</code> 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 <code>n**(1/exp)</code>; the method rounds to nearest float by appropriate care of residual.</p>

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

<pre><code>Integer&gt;&gt;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 &gt;= guess ] whileFalse:
            [ guess := nextGuess ].
    ( guess raisedTo: aPositiveInteger) &gt; self  ifTrue:
            [ guess := guess - 1 ].
    ^guess
</code></pre>

<p>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:</p>

<pre><code>Integer&gt;&gt;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 &gt; ((guess + 1 raisedTo: aPositiveInteger) + (guess raisedTo: aPositiveInteger))
        ifTrue: [guess + 1]
        ifFalse: [guess]
</code></pre>

<p>Then exactness is tested in nthRoot:</p>

<pre><code>Integer&gt;&gt;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 &lt;= 0 ifTrue: [ ^(Fraction numerator: guess * 4 - excess sign denominator: 4) asFloat].

    scaled := self &lt;&lt; (nBits * aPositiveInteger).
    guess := scaled nthRootRounded: aPositiveInteger.
    excess := (guess raisedTo: aPositiveInteger) - scaled.
    ^(Fraction numerator: guess * 4 - excess sign denominator: 1 &lt;&lt; (nBits + 2)) asFloat
</code></pre>

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

<p>It works for large integers like:</p>

<ul>
<li><code>(10 raisedTo: 600) nthRoot: 300</code> -> <code>100</code> "exact"</li>
<li><code>(10 raisedTo: 600) + 1 nthRoot: 300</code> -> <code>100.0</code> "inexact"</li>
</ul>

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

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

<p>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).</p>

<ul>
<li><a href="http://smallissimo.blogspot.fr/2011/09/clarifying-and-optimizing.html" rel="nofollow">http://smallissimo.blogspot.fr/2011/09/clarifying-and-optimizing.html</a></li>
<li><a href="http://smallissimo.blogspot.fr/2011/09/reviewing-fraction-asfloat.html" rel="nofollow">http://smallissimo.blogspot.fr/2011/09/reviewing-fraction-asfloat.html</a></li>
</ul>

python algorithm floating-point precision nth-root