具有多个根的编程语言
Programming language with multiple roots
2^(-1/3)的答案是三个根:
0.79370、-0.39685-0.68736i 和 0.39685+0.68736i(大约)
在 Wolfram Alpha 查看正确答案。
我知道几种支持复数的语言,但它们都只是return三个结果中的第一个:
Python:
>>> complex(2,0)**(-1/3)
(0.7937005259840998-0j)
八度:
>> (2+0i)^(-1/3)
ans = 0.79370
朱莉娅:
julia> complex(2,0)^(-1/3)
0.7937005259840998 + 0.0im
我正在寻找的是:
>> 2^(-1/3)
[0.79370+0i, -0.39685-0.68736i, 0.39685+0.68736i]
是否有一种编程语言(带有 REPL)可以正确地 return 所有三个根,而不必求助于任何特殊的模块或库,并且还有可用的开源实现?
下面是求解 b1/n 的所有根的方法通过多项式 xn - b 的根Matlab 的 roots
or Octave's roots
:
b = 2;
n = -3; % for b^(1/n)
c = [1 zeros(1,abs(n)-1) -b];
r = roots(c).^sign(n);
哪个returns
r =
-0.396850262992050 - 0.687364818499301i
-0.396850262992050 + 0.687364818499301i
0.793700525984100 + 0.000000000000000i
或者,使用 roots of unity(不确定这在数值上有多稳健):
b = 2;
n = -3;
n0 = abs(n);
r0 = b^(1/n0);
w = exp(2*pi*1i/n0);
r = (r0*w.^(0:n0-1).').^sign(n)
或者使用 Matlab 的符号数学工具箱:
b = 2;
n = -3;
c = [1 zeros(1,abs(n)-1) -b];
r = solve(poly2sym(c)).^sign(n)
哪个returns:
r =
2^(2/3)/2
2^(2/3)/(2*((3^(1/2)*1i)/2 - 1/2))
-2^(2/3)/(2*((3^(1/2)*1i)/2 + 1/2))
在某些情况下,您可能还会发现 nthroot
helpful (Octave documentation)。
正如许多评论所解释的那样,想要一种通用语言默认给出复杂根函数每个分支的结果可能是一项艰巨的任务。但是 Julia 非常自然地允许 specializing/overloading 运算符(因为即使是开箱即用的实现也通常是用 Julia 编写的)。具体来说:
using Roots,Polynomials # Might need to Pkg.add("Roots") first
import Base: ^
^{T<:AbstractFloat}(b::T, r::Rational{Int64}) =
roots(poly([0])^r.den - b^abs(r.num)).^sign(r.num)
现在当尝试将浮点数提高到有理数时:
julia> 2.0^(-1//3)
3-element Array{Complex{Float64},1}:
-0.39685-0.687365im
-0.39685+0.687365im
0.793701-0.0im
请注意,将 ^
的定义专门化为有理指数可以解决评论中提到的舍入问题。
2^(-1/3)的答案是三个根:
0.79370、-0.39685-0.68736i 和 0.39685+0.68736i(大约)
在 Wolfram Alpha 查看正确答案。
我知道几种支持复数的语言,但它们都只是return三个结果中的第一个:
Python:
>>> complex(2,0)**(-1/3)
(0.7937005259840998-0j)
八度:
>> (2+0i)^(-1/3)
ans = 0.79370
朱莉娅:
julia> complex(2,0)^(-1/3)
0.7937005259840998 + 0.0im
我正在寻找的是:
>> 2^(-1/3)
[0.79370+0i, -0.39685-0.68736i, 0.39685+0.68736i]
是否有一种编程语言(带有 REPL)可以正确地 return 所有三个根,而不必求助于任何特殊的模块或库,并且还有可用的开源实现?
下面是求解 b1/n 的所有根的方法通过多项式 xn - b 的根Matlab 的 roots
or Octave's roots
:
b = 2;
n = -3; % for b^(1/n)
c = [1 zeros(1,abs(n)-1) -b];
r = roots(c).^sign(n);
哪个returns
r =
-0.396850262992050 - 0.687364818499301i
-0.396850262992050 + 0.687364818499301i
0.793700525984100 + 0.000000000000000i
或者,使用 roots of unity(不确定这在数值上有多稳健):
b = 2;
n = -3;
n0 = abs(n);
r0 = b^(1/n0);
w = exp(2*pi*1i/n0);
r = (r0*w.^(0:n0-1).').^sign(n)
或者使用 Matlab 的符号数学工具箱:
b = 2;
n = -3;
c = [1 zeros(1,abs(n)-1) -b];
r = solve(poly2sym(c)).^sign(n)
哪个returns:
r =
2^(2/3)/2
2^(2/3)/(2*((3^(1/2)*1i)/2 - 1/2))
-2^(2/3)/(2*((3^(1/2)*1i)/2 + 1/2))
在某些情况下,您可能还会发现 nthroot
helpful (Octave documentation)。
正如许多评论所解释的那样,想要一种通用语言默认给出复杂根函数每个分支的结果可能是一项艰巨的任务。但是 Julia 非常自然地允许 specializing/overloading 运算符(因为即使是开箱即用的实现也通常是用 Julia 编写的)。具体来说:
using Roots,Polynomials # Might need to Pkg.add("Roots") first
import Base: ^
^{T<:AbstractFloat}(b::T, r::Rational{Int64}) =
roots(poly([0])^r.den - b^abs(r.num)).^sign(r.num)
现在当尝试将浮点数提高到有理数时:
julia> 2.0^(-1//3)
3-element Array{Complex{Float64},1}:
-0.39685-0.687365im
-0.39685+0.687365im
0.793701-0.0im
请注意,将 ^
的定义专门化为有理指数可以解决评论中提到的舍入问题。