Julia の小ネタ--045 3乗根

a の 3 乗根（立方根）は，3乗すると a になる数だ。Julia には cbrt という関数がある。また，定義により a ^ (1/3) でもある。

どちらも同じか？

正の数の場合は同じだ。

julia> cbrt(8)

2.0

julia> 8^(1/3)

2.0

負の数の場合は違う。負の数を 1/3 乗することはできない。

julia> cbrt(-8)

-2.0

julia> (-8)^(1/3)

ERROR: 

DomainError with -8.0:

Exponentiation yielding a complex result requires a complex argument.

Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

指示に従って以下のようにする。

julia> (-8 + 0im)^(1/3)

1.0 + 1.7320508075688772im

得られる数を 3 乗すると（誤差の範囲内で）ちゃんと -8 になる。

julia> a = (-8.0 + 0.0im)^(1/3)

1.0 + 1.7320508075688772im

julia> a^3

-7.999999999999998 + 8.881784197001252e-16im

-------

他の言語ではどうなるか?

Octave

cbrt は Julia と同じ結果を返す。

>> cbrt(-8)
ans = -2

負の数を 1/3 乗しようとしても，Julia のようには文句を言わず，答えを返してくれる。

>> a = (-8)^(1/3)
a =  1.000000000000000 + 1.732050807568877i
>> a^3
ans = -7.999999999999995e+00 + 9.797174393178820e-16i

Python

numpy.cbrt も Julia と同じ結果を返す。

>>> import numpy as np
>>> np.cbrt(-8)
-2.0

負の数の 1/3 乗のときは，Octave と同じふるまいをする。

>>> a = (-8)**(1/3)
>>> a
(1.0000000000000002+1.7320508075688772j)
>>> a**3
(-8+3.1086244689504383e-15j)

R

R には cbrt 関数はない(SparkR というライブラリにあるようだが，R 4.2.1 には対応していないようだ）。

負の数の 1/3 乗のときは，Julia のように親切なアドバイスをくれるわけでもなく，ぶっきらぼうに NaN を返す。

> options(digits=16)
> (-8)^(1/3)
[1] NaN

虚数にすれば，答えを出してくれる。

> a = (-8 + 0i)^(1/3)
> a
[1] 1+1.732050807568877i
> a^3
[1] -7.999999999999996e+00+2e-15i