今回は引数に複素数を受け取る数学関数 (複素関数) を取り上げます。複素関数は .NET のライブラリ System.Numerics.Complex に用意されていますが、私達でもプログラムすることができます。ここでは、拙作のページ F# 入門: 複素数 で作成した構造体 Complex に複素関数を追加することにします。なお、複素数の四則演算は float と Complex でも可能なように多重定義しています。詳細は プログラムリスト をお読みくださいませ。
複素数を引数にとる指数関数はオイラー (Euler) の公式から導くことができます。
eiθ = cos θ + i sin θ -- オイラーの公式 ex+iy = ex * eiy = ex * (cos y + i sin y)
複素数の対数関数は複素数 z を絶対値 |z| と偏角 θ を使って導くことができます。
x + iy = |z| * eiθ loge (x + iy) = loge (|z| * eiθ) = loge |z| + loge eiθ = loge |z| + iθ, (-pi <= θ <= pi)
複素数 x, y のべき乗 xy は ey*log x で求めることができます。これをそのままプログラムすると次のようになります。
リスト : 指数関数、対数関数、べき乗 // 指数関数 Exp(z) static member Exp(z: Complex) = let a = exp(z.Real) Complex(a * cos(z.Imag), a * sin(z.Imag)) // 対数関数 Log(z) static member Log(z: Complex) = Complex(log(z.Magnitude()), z.Phase()) // べき乗 Pow(x, y) static member Pow(x: Complex, y: Complex) = Complex.Exp(y * Complex.Log(x))
簡単な例を示しましょう。
> let pi = System.Math.PI;; val pi: float = 3.141592654 > Complex.Exp(Complex(0.0, pi/4.0));; val it: Complex = complex(0.707107, 0.707107) {Imag = 0.7071067812; Real = 0.7071067812;} > Complex.Exp(Complex(0.0, pi/2.0));; val it: Complex = complex(6.12323e-17, 1) {Imag = 1.0; Real = 6.123233996e-17;} > Complex.Exp(Complex(0.0, pi));; val it: Complex = complex(-1, 1.22465e-16) {Imag = 1.224646799e-16; Real = -1.0;} > Complex.Exp(Complex(1.0, 1.0));; val it: Complex = complex(1.46869, 2.28736) {Imag = 2.287355287; Real = 1.46869394;} > Complex.Log(Complex(1.0, 1.0));; val it: Complex = complex(0.346574, 0.785398) {Imag = 0.7853981634; Real = 0.3465735903;} > Complex.Log(Complex(1.0, 0.0));; val it: Complex = complex(0, 0) {Imag = 0.0; Real = 0.0;} > Complex.Log(Complex(0.0, 1.0));; val it: Complex = complex(0, 1.5708) {Imag = 1.570796327; Real = 0.0;} > Complex.Log(Complex(1.0, -1.0));; val it: Complex = complex(0.346574, -0.785398) {Imag = -0.7853981634; Real = 0.3465735903;} > Complex.Log(Complex(1e300, 1e300));; val it: Complex = complex(691.122, 0.785398) {Imag = 0.7853981634; Real = 691.1221015;} > let a = Complex(1, 1);; val a: Complex = complex(1, 1) > Complex.Pow(a, Complex(0, 0));; val it: Complex = complex(1, 0) {Imag = 0.0; Real = 1.0;} > Complex.Pow(a, Complex(1, 0));; val it: Complex = complex(1, 1) {Imag = 1.0; Real = 1.0;} > Complex.Pow(a, Complex(2, 0));; val it: Complex = complex(1.22465e-16, 2) {Imag = 2.0; Real = 1.224646799e-16;} > Complex.Pow(a, Complex(3, 0));; val it: Complex = complex(-2, 2) {Imag = 2.0; Real = -2.0;} > Complex.Pow(a, a);; val it: Complex = complex(0.273957, 0.583701) {Imag = 0.5837007588; Real = 0.2739572538;} > Complex.Pow(Complex(1, 2), Complex(3, 4));; val it: Complex = complex(0.12901, 0.0339241) {Imag = 0.03392409291; Real = 0.1290095941;}
関数 log z (z = x + yi) は負の実軸 (-∞ < x < 0) において、x + 0.0i と x - 0.0i では値が異なります。
> Complex.Log(Complex(-1.0, 0.0));; val it: Complex = complex(0, 3.14159) {Imag = 3.141592654; Real = 0.0;} > Complex.Log(Complex(-1.0, -0.0));; val it: Complex = complex(0, -3.14159) {Imag = -3.141592654; Real = 0.0;} > Complex.Log(Complex(-1e300, 0.0));; val it: Complex = complex(690.776, 3.14159) {Imag = 3.141592654; Real = 690.7755279;} > Complex.Log(Complex(-1e300, -0.0));; val it: Complex = complex(690.776, -3.14159) {Imag = -3.141592654; Real = 690.7755279;}
このように、関数 log z は負の実軸上で 2 つの値を持ちます。数学では値を一つ返す関数を「一価関数」、複数の値を返す関数を「多価関数」といいます。ここで、定義域を制限することで多価関数を一価関数にみなすことを考えます。関数 log z の場合、負の実軸を定義域から取り除けば、log z を一価関数とみなすことができるわけです。
参考 URL 8 によると、この取り除いた領域を branch cut と呼ぶそうです。プログラミングでは branch cut を定義域から取り除くのではなく、その領域では不連続な関数とするそうです。参考文献 1 では「分枝切断線」、Python のドキュメントでは「分枝切断」と記述されています。本稿では branch cut を「分枝切断」と記述することにします。
プログラミング言語の場合、0.0 と -0.0 を区別する処理系であれば、.NET のように 2 つの値を区別することができます。0.0 と -0.0 を区別しない処理系では、偏角 θ の範囲を -pi < θ <= pi に制限することで、log z の返り値を (-pi を取り除いて) 一つにすることができます。
複素数の三角関数の定義は、オイラーの公式から導かれる式の θ を複素数 z に変えたものになります。
eiθ = cos θ + i sin θ -- (1) ei(-θ) = cos θ + i sin -θ (sin -θ = -sin θ, cos -θ = cos θ) = cos θ - i sin θ -- (2) (1) + (2) eiθ + e-iθ = 2 cos θ cos θ = (eiθ + e-iθ) / 2 (1) - (2) eiθ - e-iθ = 2i sin θ sin θ = (eiθ - e-iθ) / 2i θ を複素数 z に置き換えた式が三角関数の定義になる sin z = (eiz - e-iz) / 2i cos z = (eiz + e-iz) / 2
sin z, cos z に純虚数 ix を与えると双曲線関数 (sinh x, cosh x) になります。
双曲線関数の定義 sinh x = (ex - e-x) / 2 cosh x = (ex + e-x) / 2 sin ix = (eiix - e-iix) / 2i = (e-x - ex) / 2i (分子と分母に -i を掛ける) = i (ex - e-x) / 2 = i sinh x cos ix = (eiix + e-iix) / 2 = (e-x + ex) / 2 = cosh x
これに三角関数の加法定理 [*1] を使うと次の式が導かれます。
sin(x + iy) = sin x cos iy + cos x sin iy = sin x cosh y + i cos x sinh y cos(x + iy) = cos x cos iy - sin x sin iy = cos x cosh y - i sin x sinh y tan(x + iy) = (sin 2x + sin 2iy) / (cos 2x + cos 2iy) = (sin 2x + i sinh 2y) / (cos 2x + cosh 2y)
これをそのままプログラムすると、次のようになります。
リスト : 三角関数 static member Sin(z: Complex) = let x = z.Real let y = z.Imag Complex(sin(x) * cosh(y), cos(x) * sinh(y)) static member Cos(z: Complex) = let x = z.Real let y = z.Imag Complex(cos(x) * cosh(y), - sin(x) * sinh(y)) static member Tan(z: Complex) = let x = z.Real let y = z.Imag let d = cos(2.0 * x) + cosh(2.0 * y) Complex(sin(2.0 * x) / d, sinh(2.0 * y) / d)
簡単な実行例を示します。
> let a = [ Complex(0, 1); Complex(0, -1); Complex(1, 1); Complex(1, -1)];; val a: Complex list = [complex(0, 1); complex(0, -1); complex(1, 1); complex(1, -1)] > for x in a do printfn "%A" (Complex.Sin(x));; complex(0, 1.1752) complex(0, -1.1752) complex(1.29846, 0.634964) complex(1.29846, -0.634964) val it: unit = () > for x in a do printfn "%g" (Complex.Sin(x).Magnitude());; 1.1752 1.1752 1.4454 1.4454 val it: unit = () > for x in a do printfn "%A" (Complex.Cos(x));; complex(1.54308, -0) complex(1.54308, 0) complex(0.83373, -0.988898) complex(0.83373, 0.988898) val it: unit = () > for x in a do printfn "%g" (Complex.Cos(x).Magnitude());; 1.54308 1.54308 1.29345 1.29345 val it: unit = () > for x in a do printfn "%A" (Complex.Tan(x));; complex(0, 0.761594) complex(0, -0.761594) complex(0.271753, 1.08392) complex(0.271753, -1.08392) val it: unit = () > for x in a do printfn "%g" (Complex.Tan(x).Magnitude());; 0.761594 0.761594 1.11747 1.11747 val it: unit = ()
複素数の双曲線関数の定義は、実数の定義で引数 x を複素数 z に変えたものになります。
双曲線関数の定義 (z は複素数) sinh z = (ez - e-z) / 2 cosh z = (ez + e-z) / 2
sinh z, cosh z に純虚数 ix を与えると三角関数 (sin x, cos x) になります。
sinh ix = (eix - e-ix) / 2 = i (eix - e-ix) / 2i (分子と分母に i を掛ける) = i sin x cosh ix = (eix + eix) / 2 = cos x
これに双曲線関数の加法定理を使うと、次の式が導かれます。
双曲線関数の加法定理 sinh(x + y) = sinh x cosh y + cosh x sinh y cosh(x + y) = cosh x cosh y + sinh x sinh y tanh(x + y) = sinh(x + y) / cosh(x + y) = (sinh 2x + sinh 2y) / (cosh 2x + cosh 2y) sinh(x + iy) = sinh x cos y + i cosh x sin y cosh(x + iy) = cosh x cos y + i sinh x sin y tanh(x + iy) = sinh(x + iy) / cosh(x + iy) = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y)
これをそのままプログラムすると、次のようになります。
リスト : 双曲線関数 static member Sinh(z: Complex) = let x = z.Real let y = z.Imag Complex(sinh(x) * cos(y), cosh(x) * sin(y)) static member Cosh(z: Complex) = let x = z.Real let y = z.Imag Complex(cosh(x) * cos(y), sinh(x) * sin(y)) static member Tanh(z: Complex) = let x = z.Real let y = z.Imag let d = cosh(2.0 * x) + cos(2.0 * y) Complex(sinh(2.0 * x) / d, sin(2.0 * y) / d)
簡単な使用例を示します。
> let a = [ Complex(0, 1); Complex(0, -1); Complex(1, 1); Complex(1, -1)];; val a: Complex list = [complex(0, 1); complex(0, -1); complex(1, 1); complex(1, -1)] > for x in a do printfn "%A" (Complex.Sinh(x));; complex(0, 0.841471) complex(0, -0.841471) complex(0.634964, 1.29846) complex(0.634964, -1.29846) val it: unit = () > for x in a do printfn "%A" (Complex.Cosh(x));; complex(0.540302, 0) complex(0.540302, -0) complex(0.83373, 0.988898) complex(0.83373, -0.988898) val it: unit = () > for x in a do printfn "%A" (Complex.Tanh(x));; complex(0, 1.55741) complex(0, -1.55741) complex(1.08392, 0.271753) complex(1.08392, -0.271753) val it: unit = ()
複素数 z の平方根は次の式で求めることができます。
z = x + iy, |z| = √(x2 + y2), θ = z の偏角 (-pi <= θ <= pi) とすると √(x + iy) = √(|z| * eiθ) = √|z| * eiθ/2 -- (1) √(|z| * eiθ+ 2pi) = √|z| * eiθ/2 + pi -- (2)
式 (1) を平方根の主値といいます。角度は 2pi を足すと同じ角度になるので、式 (2) がもう一つの解になります。三角関数の半角の公式を使うと、式 (1) から次の式が導かれます。
三角関数の半角の公式 sin2(θ/2) = (1 - cos θ) / 2 cos2(θ/2) = (1 + cos θ) / 2 y >= 0 の場合 √|z| * eiθ/2 = √|z| * (cos(θ/2) + i sin(θ/2)) = √|z| * (√((1 + cos θ)/2) + i √((1 - cos θ)/2)) = √((|z| + |z|cos θ)/2 + i √((|z| - |z|cos θ)/2) |z|cos θ = x だから = √((|z| + x) / 2) + i √((|z| - x) / 2) y < 0 の場合、虚部の符号が - になる = √((|z| + x) / 2) - i √((|z| - x) / 2)
これをそのままプログラムすると、次のようになります。
リスト : 平方根 static member Sqrt(z: Complex z) = let x = z.Real; let a = z.Magnitude(); let b = sqrt((a - x) / 2.0); Complex(sqrt((a + x) / 2.0), if System.Double.IsNegative(z.Imag)) then -b else b)
虚部の符号を判定するとき IsNegative を使うと、-0.0 にも対応することができます。
簡単な実行例を示します。
> let one = Complex(1, 0);; val one: Complex = complex(1, 0) > Complex.Sqrt(one);; val it: Complex = complex(1, 0) {Imag = 0.0; Real = 1.0;} > Complex.Sqrt(one * Complex(2, 0));; val it: Complex = complex(1.41421, 0) {Imag = 0.0; Real = 1.414213562;} > Complex.Sqrt(one * Complex(3, 0));; val it: Complex = complex(1.73205, 0) {Imag = 0.0; Real = 1.732050808;} > Complex.Sqrt(one * Complex(-1, 0));; val it: Complex = complex(0, 1) {Imag = 1.0; Real = 0.0;} > Complex.Sqrt(one * Complex(-2, 0));; val it: Complex = complex(0, 1.41421) {Imag = 1.414213562; Real = 0.0;} > Complex.Sqrt(one * Complex(-3, 0));; val it: Complex = complex(0, 1.73205) {Imag = 1.732050808; Real = 0.0;} > let a = [ Complex(0, 1); Complex(0, -1); Complex(1, 1); Complex(1, -1)];; val a: Complex list = [complex(0, 1); complex(0, -1); complex(1, 1); complex(1, -1)] > for x in a do printfn "%A" (Complex.Sqrt(x));; complex(0.707107, 0.707107) complex(0.707107, -0.707107) complex(1.09868, 0.45509) complex(1.09868, -0.45509) val it: unit = () > for x in a do - let y = Complex.Sqrt(x) - printfn "%A" (y * y);; complex(0, 1) complex(0, -1) complex(1, 1) complex(1, -1) val it: unit = ()
sqrt は log と同じ分枝切断を持っています。x を負の整数とすると sqrt(x) の解は i √x になりますが、もうひとつ -i √x という解があります。
> Complex.Sqrt(Complex(-2.0, 0.0));; val it: Complex = complex(0, 1.41421) {Imag = 1.414213562; Real = 0.0;} > Complex.Sqrt(Complex(-2.0, -0.0));; val it: Complex = complex(0, -1.41421) {Imag = -1.414213562; Real = 0.0;} > Complex.Sqrt(Complex(-1e300, 0.0));; val it: Complex = complex(0, 1e+150) {Imag = 1e+150; Real = 0.0;} > Complex.Sqrt(Complex(-1e300, -0.0));; val it: Complex = complex(0, -1e+150) {Imag = -1e+150; Real = 0.0;}
三角関数の逆関数を「逆三角関数 (inverse trigonometric function)」といいます。以下に F# に用意されている逆三角関数を示します。
ここでは引数 x を実数とします。asin x は引数 x が与えられたとき sin w = x となる角度 w を求めます。同様に acos x は cos w = x となる角度 w を、atan x は tan x = w となる角度 w を求めます。三角関数には周期性があるので、上式を満たす角度 w は無数に存在します。つまり、逆三角関数の返り値は無数にあることになりますが、通常は一つの値を返すように範囲を制限します。これを「主値」といいます。
逆三角関数の主値を以下に示します。
簡単な実行例を示します。
> for x in [-1.0; -0.5; 0.0; 0.5; 1.0] do printfn "%g" (asin(x));; -1.5708 -0.523599 0 0.523599 1.5708 val it: unit = () > for x in [-1.0; -0.5; 0.0; 0.5; 1.0] do printfn "%g" (acos(x));; 3.14159 2.0944 1.5708 1.0472 0 val it: unit = () > for x in [-1e300; -1.0; 0.0; 1.0; 1e300] do printfn "%g" (atan(x));; -1.5708 -0.785398 0 0.785398 1.5708 val it: unit = ()
複素数の逆三角関数の定義は、複素数の三角関数の定義から導くことができます。asin z の定義は次のようになります。
asin z = w (z, w は複素数) z = sin w = (eiw - e-iw) / 2i 両辺に 2ieiw を掛けると 2iz(eiw) = (eiw)2 - 1 (eiw)2 - 2iz(eiw) - 1 = 0 eiw の二次方程式と考えて解くと eiw = iz ± √(1 - z2) 両辺の対数をとって -i を掛け算し、平方根の主値 (+) を選ぶ w = asin z = -i log(iz + √(1 - z2))
acos z と atan z は定義だけを示します。
acos z = pi/2 - asin z atan z = i(log(1 - iz) - log(1 + iz)) / 2
asin, acos, atan は次に示す分枝切断を持っています。
プログラムは次のようになります。
リスト : 逆三角関数 static member Asin(z: Complex) = let z1 = Complex.Sqrt(1.0 + z) let z2 = Complex.Sqrt(1.0 - z) let i = Complex(0, 1) -1.0 * i * Complex.Log(i * z + z1 * z2) static member Acos(z: Complex) = System.Math.PI / 2.0 - Complex.Asin(z) static member Atan(z: Complex) = let i = Complex(0, 1) let y = i * z let x = i * (Complex.Log(1.0 - y) - Complex.Log(1.0 + y)) / 2.0 if z.Real = 0.0 && System.Double.IsNegative(z.Real) && z.Imag > 1.0 then Complex(-x.Real, x.Imag) else x
Atan の場合、公式のままプログラムすると分枝切断 (虚軸上の 1 < ∞ ) で不具合 [*2] が生じます。そのため、z の実部が負のゼロかつ虚部が 1.0 よりも大きいときは、計算結果 x の実部の符号を反転しています。
簡単な実行例を示します。
> let a = [ Complex(0, 1); Complex(0, -1); Complex(1, 1); Complex(1, -1)];; val a: Complex list = [complex(0, 1); complex(0, -1); complex(1, 1); complex(1, -1)] > for x in a do Complex.Asin(x) |> printfn "%A";; complex(0, 0.881374) complex(0, -0.881374) complex(0.666239, 1.06128) complex(0.666239, -1.06128) val it: unit = () > for x in a do Complex.Asin(x) |> Complex.Sin |> printfn "%A";; complex(0, 1) complex(0, -1) complex(1, 1) complex(1, -1) val it: unit = () > Complex.Asin(Complex(2.0, 0.0));; val it: Complex = complex(1.5708, 1.31696) {Imag = 1.316957897; Real = 1.570796327;} > Complex.Asin(Complex(2.0, -0.0));; val it: Complex = complex(1.5708, -1.31696) {Imag = -1.316957897; Real = 1.570796327;} > Complex.Asin(Complex(-4.0, 0.0));; val it: Complex = complex(-1.5708, 2.06344) {Imag = 2.063437069; Real = -1.570796327;} > Complex.Asin(Complex(-4.0, -0.0));; val it: Complex = complex(-1.5708, -2.06344) {Imag = -2.063437069; Real = -1.570796327;} > for x in a do Complex.Acos(x) |> printfn "%A";; complex(1.5708, -0.881374) complex(1.5708, 0.881374) complex(0.904557, -1.06128) complex(0.904557, 1.06128) val it: unit = () > for x in a do Complex.Acos(x) |> Complex.Cos |> printfn "%A";; complex(8.65956e-17, 1) complex(8.65956e-17, -1) complex(1, 1) complex(1, -1) val it: unit = () > Complex.Acos(Complex(2.0, 0.0));; val it: Complex = complex(0, -1.31696) {Imag = -1.316957897; Real = 0.0;} > Complex.Acos(Complex(2.0, -0.0));; val it: Complex = complex(0, 1.31696) {Imag = 1.316957897; Real = 0.0;} > Complex.Acos(Complex(-4.0, 0.0));; val it: Complex = complex(3.14159, -2.06344) {Imag = -2.063437069; Real = 3.141592654;} > Complex.Acos(Complex(-4.0, -0.0));; val it: Complex = complex(3.14159, 2.06344) {Imag = 2.063437069; Real = 3.141592654;} > for x in a do Complex.Atan(x) |> printfn "%A";; complex(NaN, Infinity) complex(NaN, -Infinity) complex(1.01722, 0.402359) complex(1.01722, -0.402359) val it: unit = () > for x in a do Complex.Atan(x) |> Complex.Tan |> printfn "%A";; complex(NaN, NaN) complex(NaN, NaN) complex(1, 1) complex(1, -1) val it: unit = () > Complex.Atan(Complex(0.0, 2.0));; val it: Complex = complex(1.5708, 0.549306) {Imag = 0.5493061443; Real = 1.570796327;} > Complex.Atan(Complex(-0.0, 2.0));; val it: Complex = complex(-1.5708, 0.549306) {Imag = 0.5493061443; Real = -1.570796327;} > Complex.Atan(Complex(0.0, -4.0));; val it: Complex = complex(1.5708, -0.255413) {Imag = -0.2554128119; Real = 1.570796327;} > Complex.Atan(Complex(-0.0, -4.0));; val it: Complex = complex(-1.5708, -0.255413) {Imag = -0.2554128119; Real = -1.570796327;}
(0.0 + 1.0 i) * (0.0 + b i) = 0.0 * 0.0 - 1.0 * b + (0.0 * b + 1.0 * 0.0) = -b + (0.0 + 0.0) i = -b + 0.0i (0.0 + 1.0 i) * (-0.0 + b i) = 0.0 * -0.0 - 1.0 * b + (0.0 * b + 1.0 * -0.0) = -b + (0.0 + -0.0) i = -b + 0.0iIEEE 754 では、-0.0 + -0.0 は -0.0 になるが、0.0 + -0.0 や -0.0 + 0.0 は 0.0 になる。
> 0.0 + 0.0;; val it: float = 0.0 > -0.0 + 0.0;; val it: float = 0.0 > 0.0 + -0.0;; val it: float = 0.0 > -0.0 + -0.0;; val it: float = -0.0
双曲線関数の逆関数を「逆双曲線関数 (inverse hyperbolic function)」といいます。F# には逆双曲線関数が定義されていませんが、.NET のライブラリ System.Math には逆双曲線関数 Asinh, Acosh, Atanh があります。双曲線関数と逆双曲線関数の定義域と値域を示します。
x と y は実数です。x = cosh y の逆関数 y = acosh x を満たす y の値は 2 つありますが、ここでは y >= 0 を主値として選ぶことにします。
逆双曲線関数の定義は双曲線関数の定義から導くことができます。
asinh x = y y = sinh x = (ex - e-x) / 2 両辺に 2ex を掛けると 2yex = e2x - 1 e2x - 2yex - 1 = 0, (ex の 2 次方程式として解く) ex = (2y ± √(4y2 + 4))/2 = y ± √(y2 + 1) ex > 0 だから平方根の符号に + を選んで対数をとる x = log (y + √(y2 + 1))
acosh x = y y = cosh x = (ex + e-x) / 2 両辺に 2ex を掛けると 2yex = e2x + 1 e2x - 2yex + 1 = 0, (ex の 2 次方程式として解く) ex = (2y ± √(4y2 - 4))/2 = y ± √(y2 - 1) x = log (y ± √(y2 - 1)) 平方根の符号 + を主値としてとると x = log (y + √(y2 - 1)), (y >= 1)
atanh x = y y = tanh x = sinh x / cosh x = (ex - e-x) / (ex + e-x) = (e2x - 1) / (e2x + 1) y(e2x + 1) = (e2x - 1) (1 - y)e2x = 1 + y e2x = (1 + y) / (1 - y) x = (1 / 2) * log((1 + y) / (1 - y))
また、次の関係式を使って atanh から asinh と acosh を求めることができます。
asinh x = atanh(x / √(x2 + 1)) acosh x = atanh(√(1 - x2) / x)
簡単な実行例を示します。
> for x in [1.5; 1.0; 0.5; 0.0; -0.5; -1.0; -1.5] do System.Math.Asinh(x) |> printfn "%g";; 1.19476 0.881374 0.481212 0 -0.481212 -0.881374 -1.19476 val it: unit = () > for x in [1.5; 1.0; 0.5; 0.0; -0.5; -1.0; -1.5] do System.Math.Asinh(x) |> sinh |> printfn "%g";; 1.5 1 0.5 0 -0.5 -1 -1.5 val it: unit = () > for x in [1.0; 1.5; 2.0; 2.5; 3.0; 3.5; 4.0] do System.Math.Acosh(x) |> printfn "%g";; 0 0.962424 1.31696 1.5668 1.76275 1.92485 2.06344 val it: unit = () > for x in [1.0; 1.5; 2.0; 2.5; 3.0; 3.5; 4.0] do System.Math.Acosh(x) |> cosh |> printfn "%g";; 1 1.5 2 2.5 3 3.5 4 val it: unit = () > for x in [-0.75; -0.5; -0.25; 0.0; 0.25; 0.5; 0.75] do System.Math.Atanh(x) |> printfn "%g";; -0.972955 -0.549306 -0.255413 0 0.255413 0.549306 0.972955 val it: unit = () > for x in [-0.75; -0.5; -0.25; 0.0; 0.25; 0.5; 0.75] do System.Math.Atanh(x) |> tanh |> printfn "%g";; -0.75 -0.5 -0.25 0 0.25 0.5 0.75 val it: unit = ()
複素数の逆双曲線関数の定義は、実数の定義で引数 x を複素数 z に変えたものになります。
asinh z = log (z + √(z2 + 1)) acosh z = log (z + √(z2 - 1)) -- (1) = log (z + (√(z + 1)) * (√(z - 1))) -- (2) [*3] atanh z = (1 / 2) * log((1 + z) / (1 - z)) = (1 / 2) * (log(1 + z) - log(1 - z))
acosh(z) = 2 * log(sqrt((z + 1)/2) + sqrt((z - 1)/2))
asinh, acosh, atanh は次に示す分枝切断を持っています。
プログラムは次のようになります。
リスト : 逆双曲線関数 static member Asinh(z: Complex) = Complex.Log(z + Complex.Sqrt(z * z + 1.0)) static member Acosh(z: Complex) = let z1 = Complex.Sqrt(z + 1.0) let z2 = Complex.Sqrt(z - 1.0) Complex.Log(z + z1 * z2) static member Atanh(z: Complex) = (Complex.Log(1.0 + z) - Complex.Log(1.0 - z)) / 2.0
簡単な実行例を示しましょう。
> let a = [ Complex(0, 1); Complex(0, -1); Complex(1, 1); Complex(1, -1)];; val a: Complex list = [complex(0, 1); complex(0, -1); complex(1, 1); complex(1, -1)] > for x in a do Complex.Asinh(x) |> printfn "%A";; complex(0, 1.5708) complex(0, -1.5708) complex(1.06128, 0.666239) complex(1.06128, -0.666239) val it: unit = () > for x in a do Complex.Asinh(x) |> Complex.Sinh |> printfn "%A";; complex(0, 1) complex(0, -1) complex(1, 1) complex(1, -1) val it: unit = () > Complex.Asinh(Complex(0.0, 2.0));; val it: Complex = complex(1.31696, 1.5708) {Imag = 1.570796327; Real = 1.316957897;} > Complex.Asinh(Complex(-0.0, 2.0));; val it: Complex = complex(-1.31696, 1.5708) {Imag = 1.570796327; Real = -1.316957897;} > Complex.Asinh(Complex(0.0, -4.0));; val it: Complex = complex(2.06344, -1.5708) {Imag = -1.570796327; Real = 2.063437069;} > Complex.Asinh(Complex(-0.0, -4.0));; val it: Complex = complex(-2.06344, -1.5708) {Imag = -1.570796327; Real = -2.063437069;} > for x in a do Complex.Acosh(x) |> printfn "%A";; complex(0.881374, 1.5708) complex(0.881374, -1.5708) complex(1.06128, 0.904557) complex(1.06128, -0.904557) val it: unit = () > for x in a do Complex.Acosh(x) |> Complex.Cosh |> printfn "%A";; complex(8.65956e-17, 1) complex(8.65956e-17, -1) complex(1, 1) complex(1, -1) val it: unit = () > Complex.Acosh(Complex(0.0, 0.0));; val it: Complex = complex(0, 1.5708) {Imag = 1.570796327; Real = 0.0;} > Complex.Acosh(Complex(0.0, -0.0));; val it: Complex = complex(0, -1.5708) {Imag = -1.570796327; Real = 0.0;} > Complex.Acosh(Complex(-4.0, 0.0));; val it: Complex = complex(2.06344, 3.14159) {Imag = 3.141592654; Real = 2.063437069;} > Complex.Acosh(Complex(-4.0, -0.0));; val it: Complex = complex(2.06344, -3.14159) {Imag = -3.141592654; Real = 2.063437069;} > for x in a do Complex.Atanh(x) |> printfn "%A";; complex(0, 0.785398) complex(0, -0.785398) complex(0.402359, 1.01722) complex(0.402359, -1.01722) val it: unit = () > for x in a do Complex.Atanh(x) |> Complex.Tanh |> printfn "%A";; complex(0, 1) complex(0, -1) complex(1, 1) complex(1, -1) val it: unit = () > Complex.Atanh(Complex(2.0, 0.0));; val it: Complex = complex(0.549306, 1.5708) {Imag = 1.570796327; Real = 0.5493061443;} > Complex.Atanh(Complex(2.0, -0.0));; val it: Complex = complex(0.549306, -1.5708) {Imag = -1.570796327; Real = 0.5493061443;} > Complex.Atanh(Complex(-4.0, 0.0));; val it: Complex = complex(-0.255413, 1.5708) {Imag = 1.570796327; Real = -0.2554128119;} > Complex.Atanh(Complex(-4.0, -0.0));; val it: Complex = complex(-0.255413, -1.5708) {Imag = -1.570796327; Real = -0.2554128119;}
// // complex.fsx : 複素数と複素関数 // // Copyright (C) 2022 Makoto Hiroi // [<Struct; CustomEquality; NoComparison>] type Complex = val Real: float val Imag: float // コンストラクタ new(a: float, b: float) = { Real = a; Imag = b } new(a: int, b: int) = { Real = float a; Imag = float b } // メソッド override this.ToString () = sprintf "complex(%g, %g)" this.Real this.Imag // 等値演算子 override this.Equals (other: obj) = let that = other :?> Complex this.Real = that.Real && this.Real = that.Real // とても適当なので実際に使ってはいけない override this.GetHashCode () = (this.Real / this.Imag).GetHashCode() // 複素共役 member this.Conjugate() = Complex(this.Real, -this.Imag) // 偏角 member this.Phase() = atan2 this.Imag this.Real // 絶対値 member this.Magnitude() = if this.Real = 0.0 then abs this.Imag else if this.Imag = 0.0 then abs this.Real else if abs this.Imag > abs this.Real then let temp = this.Real / this.Imag (abs this.Imag) * sqrt(1.0 + temp * temp) else let temp = this.Imag / this.Real (abs this.Real) * sqrt(1.0 + temp * temp) // // 四則演算 // static member ( + ) (x: Complex, y: Complex) = Complex(x.Real + y.Real, x.Imag + y.Imag) static member ( + ) (x: float, y: Complex) = Complex(x + y.Real, y.Imag) static member ( + ) (x: Complex, y: float) = Complex(x.Real + y, x.Imag) static member ( - ) (x: Complex, y: Complex) = Complex(x.Real - y.Real, x.Imag - y.Imag) static member ( - ) (x: float, y: Complex) = Complex(x - y.Real, -y.Imag) static member ( - ) (x: Complex, y: float) = Complex(x.Real - y, x.Imag) static member ( * ) (x: Complex, y: Complex) = let a = x.Real let b = x.Imag let c = y.Real let d = y.Imag Complex(a * c - b * d, b * c + a * d) static member ( * ) (x: float, y: Complex) = Complex(x * y.Real, x * y.Imag) static member ( * ) (x: Complex, y: float) = Complex(x.Real * y, x.Imag * y) static member ( / ) (x: Complex, y: Complex) = let a = x.Real let b = x.Imag let c = y.Real let d = y.Imag if abs c >= abs d then let u = d / c let v = c + d * u Complex((a + b * u) / v, (b - a * u) / v) else let u = c / d let v = c * u + d Complex((a * u + b) / v, (b * u - a) / v) static member ( / ) (x: float, y: Complex) = let c = y.Real let d = y.Imag if abs c >= abs d then let u = d / c let v = c + d * u Complex(x / v, (-x * u) / v) else let u = c / d let v = c * u + d Complex((x * u) / v, (-x) / v) static member ( / ) (x: Complex, y: float) = Complex(x.Real / y, x.Imag / y) // // 複素関数 // // 指数関数 Exp(z) static member Exp(z: Complex) = let a = exp(z.Real) Complex(a * cos(z.Imag), a * sin(z.Imag)) // 対数関数 Log(z) static member Log(z: Complex) = Complex(log(z.Magnitude()), z.Phase()) // べき乗 Pow(x, y) static member Pow(x: Complex, y: Complex) = Complex.Exp(y * Complex.Log(x)) // 三角関数 static member Sin(z: Complex) = let x = z.Real let y = z.Imag Complex(sin(x) * cosh(y), cos(x) * sinh(y)) static member Cos(z: Complex) = let x = z.Real let y = z.Imag Complex(cos(x) * cosh(y), - sin(x) * sinh(y)) static member Tan(z: Complex) = let x = z.Real let y = z.Imag let d = cos(2.0 * x) + cosh(2.0 * y) Complex(sin(2.0 * x) / d, sinh(2.0 * y) / d) // 双曲線関数 static member Sinh(z: Complex) = let x = z.Real let y = z.Imag Complex(sinh(x) * cos(y), cosh(x) * sin(y)) static member Cosh(z: Complex) = let x = z.Real let y = z.Imag Complex(cosh(x) * cos(y), sinh(x) * sin(y)) static member Tanh(z: Complex) = let x = z.Real let y = z.Imag let d = cosh(2.0 * x) + cos(2.0 * y) Complex(sinh(2.0 * x) / d, sin(2.0 * y) / d) // 平方根 static member Sqrt(z: Complex) = let x = z.Real; let a = z.Magnitude(); let b = sqrt((a - x) / 2.0); Complex(sqrt((a + x) / 2.0), if System.Double.IsNegative(z.Imag) then -b else b) // 複素数の逆三角関数 static member Asin(z: Complex) = let z1 = Complex.Sqrt(1.0 + z) let z2 = Complex.Sqrt(1.0 - z) let i = Complex(0, 1) -1.0 * i * Complex.Log(i * z + z1 * z2) static member Acos(z: Complex) = System.Math.PI / 2.0 - Complex.Asin(z) static member Atan(z: Complex) = let i = Complex(0, 1) let y = i * z let x = i * (Complex.Log(1.0 - y) - Complex.Log(1.0 + y)) / 2.0 if z.Real = 0.0 && System.Double.IsNegative(z.Real) && z.Imag > 1.0 then Complex(-x.Real, x.Imag) else x // 複素数の逆双曲線関数 static member Asinh(z: Complex) = Complex.Log(z + Complex.Sqrt(z * z + 1.0)) static member Acosh(z: Complex) = let z1 = Complex.Sqrt(z + 1.0) let z2 = Complex.Sqrt(z - 1.0) Complex.Log(z + z1 * z2) static member Atanh(z: Complex) = (Complex.Log(1.0 + z) - Complex.Log(1.0 - z)) / 2.0