NumPy负数情况下的第n个奇次方根

Question

我希望在Python中计算某些数字的第n个奇次方根。Numpy提供了一个立方根函数，使用该函数可以计算x^(1/3)。

x = np.linspace(-100,100,100)
np.cbrt(x)
>>> array([-4.64158883, -4.26859722, -3.81571414, -3.21829795, -2.23144317,
    2.23144317,  3.21829795,  3.81571414,  4.26859722,  4.64158883])

然而，如果我想以直接的方式计算其他k阶奇根的相同内容，我有些困难。我无法直接使用np.power，甚至无法计算立方根：

np.power(x,1./3)
>>> array([       nan,        nan,        nan,        nan,        nan,
   2.23144317, 3.21829795, 3.81571414, 4.26859722, 4.64158883])

(-100.)**(1./3)
>>> ValueError: negative number cannot be raised to a fractional power

我可以计算x的绝对值的第k个奇次方根，然后根据x中负数项的符号进行更改，但我想知道是否有更直接的方法。这是我的当前解决方案：

def kth_root(x,k):
    if k % 2 != 0:
        res = np.power(np.abs(x),1./k)
        return res*np.sign(x)
    else:
        return np.power(np.abs(x),1./k)

kth_root(x,3)
>>> array([-4.64158883, -4.26859722, -3.81571414, -3.21829795, -2.23144317,
    2.23144317,  3.21829795,  3.81571414,  4.26859722,  4.64158883])

- benbo

---

np.sign 比 abs(x)/x 更高效。顺便说一下，你需要这样做的原因是通常计算幂时使用对数。 - PM 2Ring

你正在使用Python 2吗？在Python 3中，print((-100) ** (1/3)) # (2.3207944168063896+4.019733843830847j)。 - DeepSpace

@PM2Ring 我可能漏掉了什么，但是 OP 的 (-100.)**(1./3) 也不是 Numpy 吗？ - DeepSpace

顺便提一下，记住 ** 运算符比 - 运算符优先级高。 - DeepSpace

@DeepSpace 是的，但问题标记了Numpy，并且标题以Numpy开头。所以我想他们想要使用Numpy。 :) - PM 2Ring

Answer 1

我将用我的当前解决方案回答我的问题。这并不是说不存在更简单或更快的方法。

def kth_root(x,k):
    if k % 2 != 0:
        res = np.power(np.abs(x),1./k)
        return res*np.sign(x)
    else:
        return np.power(np.abs(x),1./k)


x = np.linspace(-100,100,100)
kth_root(x,3)
>>> array([-4.64158883, -4.26859722, -3.81571414, -3.21829795, -2.23144317,
2.23144317,  3.21829795,  3.81571414,  4.26859722,  4.64158883])

- benbo

Answer 2

我遇到了类似的问题，并从你开始定义一个函数来计算正确实域下 x 的有理指数 n/d 次幂。

def r_pow(x, n, d):
    """
    Compute x to the power of n/d (not reduced to lowest
    expression) with the correct function real domains.
    
    ARGS:
        x (int,float,array): base
        n (int)            : exponent numerator
        d (int)            : exponent denominator
        
    RETURNS:
        x to the power of n/d
    """
    
    # list to array
    if type(x) == list:
        x = np.array(x)
    # check inputs
    if type(n) != int or type(d) != int:
        raise Exception("Exponent numerator and denominator must be integers")
    # if denominator is zero
    if not d:
        raise Exception("Exponent denominator cannot be 0")
        
    # raise x to power of n
    X = x**n
    # even denominator
    if not d % 2:
        # domain is for X>=0 only
        if type(x) == np.ndarray:
            X[X<0] = np.nan
        elif X < 0:
            X = np.nan
        res = np.power(X, 1./d)
        return res
    # odd denominator
    else:
        # domain is all R
        res = np.power(np.abs(X), 1./d)
        res *= np.sign(X)
        return res

因此，例如，我们可以有 x^(2/3)，其中定义域为所有实数，函数符号始终为正（因为 x^2 始终为正）。

x = np.linspace(-10, 10, 1000)
plt.plot(x, r_pow(x, 2, 3));

或者 x^(1/3)，其中定义域为所有实数，但对于 x<0，函数符号为负。

x = np.linspace(-10, 10, 1000)
plt.plot(x, r_pow(x, 1, 3));

或 x^(1/6)，其中实数域为 [0，+∞)

x = np.linspace(-10, 10, 1000)
plt.plot(x, r_pow(x, 1, 6));

- Max Pierini