Data.Complex 包给出了 -8 + 0i ((-8):+0) 的立方根的奇怪结果
Data.Complex package gives strange result for cube root of -8 + 0i ((-8):+0)
不得不说,数学不是我的强项。我希望 How to find the cube root of a negative integer such that it does not return NaN? by using the Data.Complex 包能得到不错的结果,但当我喜欢
*Main> ((-8):+0) ** (1/3)
1.0 :+ 1.732050807568877
我希望得到类似 (-2.0):+0
的结果,其中 -2 是实部,0 是虚部。然而,结果证明有一个重要的虚部。我已经检查了 (**)
RealFloat
类型的 Complex
实例,其中指出;
x ** y = case (x,y) of
(_ , (0:+0)) -> 1 :+ 0
((0:+0), (exp_re:+_)) -> case compare exp_re 0 of
GT -> 0 :+ 0
LT -> inf :+ 0
EQ -> nan :+ nan
((re:+im), (exp_re:+_))
| (isInfinite re || isInfinite im) -> case compare exp_re 0 of
GT -> inf :+ 0
LT -> 0 :+ 0
EQ -> nan :+ nan
| otherwise -> exp (log x * y)
where
inf = 1/0
nan = 0/0
所以我们必须正常查看 exp (log x * y)
部分,其中 exp
和 log
实例 Complex
看起来像;
exp (x:+y) = expx * cos y :+ expx * sin y
where expx = exp x
log z = log (magnitude z) :+ phase z
然后我转到 magnitude
,其定义如下;
magnitude :: (RealFloat a) => Complex a -> a
magnitude (x:+y) = scaleFloat k
(sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y)))
where k = max (exponent x) (exponent y)
mk = - k
sqr z = z * z
我被困在哪里。我只是简单地喜欢 sqrt (realPart z ^ 2 + imagPart z ^ 2)
.
我做错了什么..?
嗯,你得到的答案 是 -8 的立方根,你可以通过立方看到它:
> let x = ((-8):+0)**(1/3)
> x
1.0 :+ 1.732050807568877
> x^3
(-7.9999999999999964) :+ 2.220446049250313e-15
> x*x*x
(-7.9999999999999964) :+ 2.220446049250313e-15
>
这不是您要查找的立方根。
其实-8的立方根有3个。您可以通过计算一个的三个立方根来计算它们:
> let j = (-1/2):+sqrt(3)/2 -- a cube root of 1
> let units = [1,j,j*j] -- that can generate them all
并将它们乘以上面的值:
> map (*x) units -- the cube roots of -8
[1.0 :+ 1.732050807568877,
(-1.9999999999999996) :+ 1.1102230246251565e-16,
0.9999999999999997 :+ (-1.7320508075688767)]
> map (^3) $ map (*x) units -- prove they cube to -8
[(-7.9999999999999964) :+ 2.220446049250313e-15,
(-7.999999999999995) :+ 1.3322676295501873e-15,
(-7.999999999999992) :+ 8.881784197001252e-16]
>
正如您所确定的,((-8):+0)**(1/3)
产生的立方根是 -8 的特定立方根,计算方式为:
> exp(log((-8):+0)*(1/3))
1.0 :+ 1.732050807568877
>
可能值得注意的是 (1) exp
和 log
在 Data.Complex
中的定义,虽然看起来很奇怪,但这些函数的标准数学定义是复数; (2) 对于复数,x ** y
定义为 exp(log(x)*y)
在数学上完全合理,并确保 x ** y
具有您应该期望的所有属性,包括 属性 的立方((-8):+0)**(1/3)
应该给你 -8。
至于 magnitude
的定义(在 log
的定义中使用),您的更简单的定义也可以:
> magnitude (-8)
8.0
> sqrt (realPart (-8) ^ 2 + imagPart (-8) ^ 2)
8.0
>
然而,magnitude
在 Data.Complex
中的定义已被选择,以便在涉及的数字非常大的情况下提供正确答案:
> magnitude 1e300
1.0e300
> sqrt (realPart 1e300 ^ 2 + imagPart 1e300 ^ 2)
Infinity
>
不得不说,数学不是我的强项。我希望 How to find the cube root of a negative integer such that it does not return NaN? by using the Data.Complex 包能得到不错的结果,但当我喜欢
*Main> ((-8):+0) ** (1/3)
1.0 :+ 1.732050807568877
我希望得到类似 (-2.0):+0
的结果,其中 -2 是实部,0 是虚部。然而,结果证明有一个重要的虚部。我已经检查了 (**)
RealFloat
类型的 Complex
实例,其中指出;
x ** y = case (x,y) of
(_ , (0:+0)) -> 1 :+ 0
((0:+0), (exp_re:+_)) -> case compare exp_re 0 of
GT -> 0 :+ 0
LT -> inf :+ 0
EQ -> nan :+ nan
((re:+im), (exp_re:+_))
| (isInfinite re || isInfinite im) -> case compare exp_re 0 of
GT -> inf :+ 0
LT -> 0 :+ 0
EQ -> nan :+ nan
| otherwise -> exp (log x * y)
where
inf = 1/0
nan = 0/0
所以我们必须正常查看 exp (log x * y)
部分,其中 exp
和 log
实例 Complex
看起来像;
exp (x:+y) = expx * cos y :+ expx * sin y
where expx = exp x
log z = log (magnitude z) :+ phase z
然后我转到 magnitude
,其定义如下;
magnitude :: (RealFloat a) => Complex a -> a
magnitude (x:+y) = scaleFloat k
(sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y)))
where k = max (exponent x) (exponent y)
mk = - k
sqr z = z * z
我被困在哪里。我只是简单地喜欢 sqrt (realPart z ^ 2 + imagPart z ^ 2)
.
我做错了什么..?
嗯,你得到的答案 是 -8 的立方根,你可以通过立方看到它:
> let x = ((-8):+0)**(1/3)
> x
1.0 :+ 1.732050807568877
> x^3
(-7.9999999999999964) :+ 2.220446049250313e-15
> x*x*x
(-7.9999999999999964) :+ 2.220446049250313e-15
>
这不是您要查找的立方根。
其实-8的立方根有3个。您可以通过计算一个的三个立方根来计算它们:
> let j = (-1/2):+sqrt(3)/2 -- a cube root of 1
> let units = [1,j,j*j] -- that can generate them all
并将它们乘以上面的值:
> map (*x) units -- the cube roots of -8
[1.0 :+ 1.732050807568877,
(-1.9999999999999996) :+ 1.1102230246251565e-16,
0.9999999999999997 :+ (-1.7320508075688767)]
> map (^3) $ map (*x) units -- prove they cube to -8
[(-7.9999999999999964) :+ 2.220446049250313e-15,
(-7.999999999999995) :+ 1.3322676295501873e-15,
(-7.999999999999992) :+ 8.881784197001252e-16]
>
正如您所确定的,((-8):+0)**(1/3)
产生的立方根是 -8 的特定立方根,计算方式为:
> exp(log((-8):+0)*(1/3))
1.0 :+ 1.732050807568877
>
可能值得注意的是 (1) exp
和 log
在 Data.Complex
中的定义,虽然看起来很奇怪,但这些函数的标准数学定义是复数; (2) 对于复数,x ** y
定义为 exp(log(x)*y)
在数学上完全合理,并确保 x ** y
具有您应该期望的所有属性,包括 属性 的立方((-8):+0)**(1/3)
应该给你 -8。
至于 magnitude
的定义(在 log
的定义中使用),您的更简单的定义也可以:
> magnitude (-8)
8.0
> sqrt (realPart (-8) ^ 2 + imagPart (-8) ^ 2)
8.0
>
然而,magnitude
在 Data.Complex
中的定义已被选择,以便在涉及的数字非常大的情况下提供正确答案:
> magnitude 1e300
1.0e300
> sqrt (realPart 1e300 ^ 2 + imagPart 1e300 ^ 2)
Infinity
>