-1

I found the following answer here on Stackoverflow:

https://stackoverflow.com/a/356187/1829329

But it only works for integers as n in nth root:

import gmpy2 as gmpy

result = gmpy.root((1/0.213), 31.5).real
print('result:', result)

results in:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-eb4628226deb> in <module>()
      8 
----> 9 result = gmpy.root((1/0.213), 31.5).real
     10 
     11 print('result:', result)

TypeError: root() requires 'mpfr','int' arguments

What is a good and precise way to calculate such a root? (This is the python code representation of some formular, which I need to use to calculate in a lecture.)

EDIT#1

Here is my solution based on Spektre's answer and information from the people over here at http://math.stackexchange.com.

import numpy as np

def naive_root(nth, a, datatype=np.float128):
    """This function can only calculate the nth root, if the operand a is positive."""
    logarithm = np.log2(a, dtype=datatype)
    exponent = np.multiply(np.divide(1, nth, dtype=datatype), logarithm, dtype=datatype)
    result = np.exp2(exponent, dtype=datatype)
    return result

def nth_root(nth, a, datatype=np.float128):
    if a == 0:
        print('operand is zero')
        return 0
    elif a > 0:
        print('a > 0')
        return naive_root(nth, a, datatype=datatype)
    elif a < 0:
        if a % 2 == 1:
            print('a is odd')
            return -naive_root(nth, np.abs(a))
        else:
            print('a is even')
            return naive_root(nth, np.abs(a))
Community
  • 1
  • 1
Zelphir Kaltstahl
  • 5,722
  • 10
  • 57
  • 86
  • When a is less than 0, and you are taking a non-integer root, the result is a complex number. Your `nth_root` function does not do that. I think you've made it more complex than needed. – casevh Dec 12 '15 at 05:12

3 Answers3

1

see Power by squaring for negative exponents

anyway as I do not code in python or gmpy here some definitions first:

  • pow(x,y) means x powered by y
  • root(x,y) means x-th root of y

As these are inverse functions we can rewrite:

  • pow(root(x,y),x)=y

equation

You can use this to check for correctness. As the functions are inverse you can write also this:

  • pow(x,1/y)=root(y,x)
  • root(1/x,y)=pow(y,x)

So if you got fractional (rational) root or power you can compute it as integer counterpart with inverse function.

Also if you got for example something like root(2/3,5) then you need to separate to integer operands first:

root(2/3,5)=pow(root(2,5),3)
 ~11.18034 = ~2.236068   ^3
 ~11.18034 = ~11.18034

For irational roots and powers you can not obtain precise result. Instead you round the root or power to nearest possible representation you can to minimize the error or use pow(x,y) = exp2(y*log2(x)) approach. If you use any floating point or fixed point decimal numbers then you can forget about precise results and go for pow(x,y) = exp2(y*log2(x)) from the start ...

[Notes]

I assumed only positive operand ... if you got negative number powered or rooted then you need to handle the sign for integer roots and powers (odd/even). For irational roots and powers have the sign no meaning (or at least we do not understand any yet).

Community
  • 1
  • 1
Spektre
  • 49,595
  • 11
  • 110
  • 380
  • Thanks for the mathematical background knowledge. Is it a bad idea to naively (simply typing what that formular says under "Derivation from Newton's method") implement the Newton Algorithm from https://en.wikipedia.org/wiki/Nth_root_algorithm using high precision floating numbers? If I understand correctly, it can also handle negative roots, so I'd not have to handle the sign, if I ever want to have a negative root. – Zelphir Kaltstahl Dec 11 '15 at 13:44
  • @Zelphir that is not a good idea (from my point of view) as the root can be directly computed by single power. In Newtons algorithm you use iteartion through multiple power calls further increasing the error of approximation in comparison to direct computation: `root(x,y)=pow(y,1/x)=exp2(x*log2(1/y)` ,,. to handle the negative `y` just compute with `|y|` and set result to negative if `1/y` or `y` is odd (`y mod 2==1` or `1/y mod 2 ==1`) the functions `exp2,log2` are directly implemented on most FPU and also directly computable ... – Spektre Dec 11 '15 at 13:54
1

If you are willing to use Python 3.x, the native pow() will do exactly what you want by just using root(x,y) = pow(x,1/y). It will automatically return a complex result if that is appropriate.

Python 3.4.3 (default, Sep 27 2015, 20:37:11)
[GCC 5.2.1 20150922] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> pow(1/0.213, 1/31.5)
1.0503191465568489
>>> pow(1/0.213, -1/31.5)
0.952091565004975
>>> pow(-1/0.213, -1/31.5)
(0.9473604081457588-0.09479770688958634j)
>>> pow(-1/0.213, 1/31.5)
(1.045099874779588+0.10457801566102139j)
>>>

Returning a complex result instead of raising a ValueError is one of changes in Python 3. If you want the same behavior with Python 2, you can use gmpy2 and enable complex results.

>>> import gmpy2
>>> gmpy2.version()
'2.0.5'
>>> gmpy2.get_context().allow_complex=True
>>> pow(1/gmpy2.mpfr("0.213"), 1/gmpy2.mpfr("31.5"))
mpfr('1.0503191465568489')
>>> pow(-1/gmpy2.mpfr("0.213"), 1/gmpy2.mpfr("31.5"))
mpc('1.0450998747795881+0.1045780156610214j')
>>> pow(-1/gmpy2.mpfr("0.213"), -1/gmpy2.mpfr("31.5"))
mpc('0.94736040814575884-0.094797706889586358j')
>>> pow(1/gmpy2.mpfr("0.213"), -1/gmpy2.mpfr("31.5"))
mpfr('0.95209156500497505')
>>> 
casevh
  • 11,093
  • 1
  • 24
  • 35
  • Huh, I expected a bigger loss of precision when using `pow`, maybe I should've tried that way first. I already updated my question to contain the function I wrote for calculating the nth root. I am actually always using Python 3. Does Python internally use any high speed high precision libraries for calculation? I used NumPy now, with float128 precision, simply because I can : ) – Zelphir Kaltstahl Dec 12 '15 at 12:57
  • I am aware that there is gmpy2 on PyPi for Python 3 as well - stumbled upon some posts on SO you wrote during my search. I am not sure how NumPy works internally and how similar that is to gmpy2. I tried some gmpy2 as well and couldn't see any difference, except gmpy being a bit more complicated to read (imo). – Zelphir Kaltstahl Dec 12 '15 at 12:59
  • 1
    `numpy` is optimized to work with arrays/vectors of types that recognized by the underlying hardware platform. the `128` in `float128` indicates the amount of memory used by the data type when stored in an array. It usually corresponds to a smaller hardware type. See http://stackoverflow.com/questions/9062562/what-is-the-internal-precision-of-numpy-float128 .`gmpy2` does arbitrary-precision calculations and does not rely on the platform's native floating-point types. – casevh Dec 12 '15 at 15:45
  • Nice to know, thank you! Next time, when I want to have better precision, I'll have a closer look on gmpy2. – Zelphir Kaltstahl Dec 13 '15 at 01:42
0

Here is something I use that seems to work with any number just fine:

root = number**(1/nthroot)
print(root)

It works with any number data type.

Akinni
  • 67
  • 1
  • 14