M.Hiroi's Home Page

F# Junk Scripts

Copyright (C) 2022 Makoto Hiroi
All rights reserved.

複素関数

今回は引数に複素数を受け取る数学関数 (複素関数) を取り上げます。複素関数は .NET のライブラリ System.Numerics.Complex に用意されていますが、私達でもプログラムすることができます。ここでは、拙作のページ F# 入門: 「複素数」で作成した構造体 Complex に複素関数を追加することにします。なお、複素数の四則演算は float と Complex でも可能なように多重定義しています。詳細はプログラムリストをお読みくださいませ。

●複素数の指数関数と対数関数

\( \begin{array}{l} e^{i\theta} = \cos \theta + i \sin \theta \quad (Euler's formula) \\ \begin{eqnarray} e^{x+iy} &=& e^x e^{iy} \\ &=& e^x (\cos y + i \sin y) \end{eqnarray} \end{array} \)

複素数の対数関数は複素数 z を絶対値 \(|z|\) と偏角 \(\theta\) を使って導くことができます。

\( \begin{array}{l} x + iy = |z| e^{i\theta} \\ \begin{eqnarray} \log_e (x + iy) &=& log_e |z| e^{i\theta} \\ &=& log_e |z| + log_e e^{i\theta} \\ &=& log_e |z| + i\theta, \quad (-\pi \leq \theta \leq \pi) \end{eqnarray} \end{array} \)

複素数 x, y のべき乗 \(x^y\) は次式で求めることができます。

\(x^y = e^{y\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) { ... }

> Complex.Exp(Complex(0.0, pi/2.0));;
val it: Complex = complex(6.12323e-17, 1) { ... }

> Complex.Exp(Complex(0.0, pi));;
val it: Complex = complex(-1, 1.22465e-16) { ... }

> Complex.Exp(Complex(1.0, 1.0));;
val it: Complex = complex(1.46869, 2.28736) { ... }

> Complex.Log(Complex(1.0, 1.0));;
val it: Complex = complex(0.346574, 0.785398) { ... }

> Complex.Log(Complex(1.0, 0.0));;
val it: Complex = complex(0, 0) { ... }

> Complex.Log(Complex(0.0, 1.0));;
val it: Complex = complex(0, 1.5708) { ... }

> Complex.Log(Complex(1.0, -1.0));;
val it: Complex = complex(0.346574, -0.785398) { ... }

> Complex.Log(Complex(1e300, 1e300));;
val it: Complex = complex(691.122, 0.785398) { ... }

> let a = Complex(1, 1);;
val a: Complex = complex(1, 1)

> Complex.Pow(a, Complex(0, 0));;
val it: Complex = complex(1, 0) { ... }

> Complex.Pow(a, Complex(1, 0));;
val it: Complex = complex(1, 1) { ... }

> Complex.Pow(a, Complex(2, 0));;
val it: Complex = complex(1.22465e-16, 2) { ... }

> Complex.Pow(a, Complex(3, 0));;
val it: Complex = complex(-2, 2) { ... }

> Complex.Pow(a, a);;
val it: Complex = complex(0.273957, 0.583701) { ... }

> Complex.Pow(Complex(1, 2), Complex(3, 4));;
val it: Complex = complex(0.12901, 0.0339241) { ... }

関数 \(\log z \ (z = x + iy)\) は負の実軸 (\(-\infty \lt x \lt 0\)) において、x + 0.0j と x - 0.0j では値が異なります。 

> Complex.Log(Complex(-1.0, 0.0));;
val it: Complex = complex(0, 3.14159) { ... }

> Complex.Log(Complex(-1.0, -0.0));;
val it: Complex = complex(0, -3.14159) { ... }

> Complex.Log(Complex(-1e300, 0.0));;
val it: Complex = complex(690.776, 3.14159) { ... }

> Complex.Log(Complex(-1e300, -0.0));;
val it: Complex = complex(690.776, -3.14159) { ... }

このように、関数 \(\log z\) は負の実軸上で 2 つの値を持ちます。数学では値を一つ返す関数を「一価関数」、複数の値を返す関数を「多価関数」といいます。ここで、定義域を制限することで多価関数を一価関数にみなすことを考えます。関数 \(\log z\) の場合、負の実軸を定義域から取り除けば、\(\log z\) を一価関数とみなすことができるわけです。

参考 URL 『逆双曲線関数と逆三角関数の branch cut | 雑記帳』によると、この取り除いた領域を branch cut と呼ぶそうです。プログラミングでは branch cut を定義域から取り除くのではなく、その領域では不連続な関数とするそうです。モジュール cmath のドキュメントには、\(\log\) は「分枝切断を一つもち、0 から負の実数軸に沿って \(-\infty\) へと延びており」と記述されています。Python のマニュアルに合わせて、本稿でも branch cut を「分枝切断」と記述することにします。

プログラミング言語の場合、0.0 と -0.0 を区別する処理系であれば、Python のように 2 つの値を区別することができます。0.0 と -0.0 を区別しない処理系では、偏角 \(\theta\) の範囲を \(-\pi \lt \theta \leq \pi\) に制限することで、\(\log z\) の返り値を (\(-\pi\) を取り除いて) 一つにすることができます。

●複素数の三角関数

複素数の三角関数の定義は、オイラーの公式から導かれる式の \(\theta\) を複素数 \(z\) に変えたものになります。

\(\sin -\theta = -\sin \theta, \ \cos -\theta = \cos \theta \) より
\( \begin{eqnarray} e^{i(-\theta)} &=& \cos -\theta + i \sin -\theta \\ &=& \cos \theta - i \sin \theta \end{eqnarray} \)

\( \begin{array}{l} e^{i\theta} + e^{-i\theta} = 2 \cos \theta \\ \cos \theta = \dfrac{e^{i\theta} + e^{-i\theta}}{2} \\ e^{i\theta} - e^{-i\theta} = 2i \sin \theta \\ \sin \theta = \dfrac{e^{i\theta} - e^{-i\theta}}{2i} \end{array} \)

\(\theta\) を複素数 z に置き換えた式が三角関数の定義になる

\( \sin z = \dfrac{e^{iz} - e^{-iz}}{2i}, \quad \cos z = \dfrac{e^{iz} + e^{-iz}}{2} \)

\(\sin z, \cos z\) に純虚数 \(ix\) を与えると双曲線関数 (\(\sinh x, \cosh x\)) になります。

双曲線関数の定義
\(\begin{eqnarray} \sinh x = \dfrac{e^x - e^{-x}}{2} \\ \cosh x = \dfrac{e^x + e^{-x}}{2} \end{eqnarray}\)
\(\begin{eqnarray} \sin ix &=& \dfrac{e^{iix} - e^{-iix}}{2i} \\ &=& \dfrac{e^{-x} - e^x}{2i} \times \dfrac{-i}{-i} \\ &=& i \dfrac{e^x - e^{-x}}{2} \\ &=& i \sinh x \\ \cos ix &=& \dfrac{e^{iix} + e^{-iix}}{2} \\ &=& \dfrac{e^{-x} + e^x}{2} \\ &=& \cosh x \end{eqnarray}\)

これに三角関数の加法定理 [*1] を使うと次の式が導かれます。

\(\begin{eqnarray} \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) &=& \dfrac{\sin 2x + \sin 2iy}{\cos 2x + \cos 2iy} \\ &=& \dfrac{\sin 2x + i \sinh 2y}{\cos 2x + \cosh 2y} \end{eqnarray}\)

これをそのままプログラムすると、次のようになります。 

リスト : 三角関数

  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 = ()
-- note --------
[*1] 三角関数の公式は引数が複素数でも成り立ちます。ただし、\(|\sin x| \leq 1, |\cos x| \leq 1\) という関係式は、x が実数だと成立しますが複素数では成立しません。

●複素数の双曲線関数

複素数の双曲線関数の定義は、実数の定義で引数 x を複素数 z に変えたものになります。

双曲線関数の定義 (z は複素数)
\(\begin{eqnarray} \sinh z = \dfrac{e^{z} - e^{-z}}{2} \\ \cosh z = \dfrac{e^{z} + e^{-z}}{2} \end{eqnarray}\)

sinh z, cosh z に純虚数 ix を与えると三角関数 (sin x, cos x) になります。

\(\begin{eqnarray} \sinh ix &=& \dfrac{e^{ix} - e^{-ix}}{2} \\ &=& \dfrac{e^{ix} - e^{-ix}}{2} \times \dfrac{i}{i} \\ &=& i \dfrac{e^{ix} - e^{-ix}}{2i} \\ &=& i \sin x \\ \cosh ix &=& \dfrac{e^{ix} + e^{ix}}{2} \\ &=& \cos x \end{eqnarray}\)

これに双曲線関数の加法定理を使うと、次の式が導かれます。

双曲線関数の加法定理
\(\begin{eqnarray} \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) &=& \dfrac{\sinh(x + y)}{\cosh(x + y)} \\ &=& \dfrac{\sinh 2x + \sinh 2y}{\cosh 2x + \cosh 2y} \end{eqnarray}\)
\(\begin{eqnarray} \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) &=& \dfrac{\sinh(x + iy)}{\cosh(x + iy)} \\ &=& \dfrac{\sinh 2x + i \sin 2y}{\cosh 2x + \cos 2y} \end{eqnarray}\)

これをそのままプログラムすると、次のようになります。 

リスト : 双曲線関数

  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| = \sqrt{x^2 + y^2}, \ \theta = \arg z \ (-\pi \leq \theta \leq \pi)\) とすると

\( \sqrt{x + iy} = \begin{cases} \sqrt{|z| e^{i\theta}} = \sqrt{|z|} e^{i\theta/2} & (1) \\ \sqrt{|z| e^{i\theta + 2\pi}} = \sqrt{|z|} e^{i\theta/2 + \pi} & (2) \end{cases} \)

式 (1) を平方根の主値といいます。角度は \(2\pi\) を足すと同じ角度になるので、式 (2) がもう一つの解になります。三角関数の半角の公式を使うと、式 (1) から次の式が導かれます。

三角関数の半角の公式
\(\begin{eqnarray} \sin^2{\left(\frac{\theta}{2}\right)} = \dfrac{1 - \cos \theta}{2} \\ \cos^2{\left(\frac{\theta}{2}\right)} = \dfrac{1 + \cos \theta}{2} \end{eqnarray}\)
\(y \geq 0\) の場合
\(\begin{eqnarray} \sqrt{|z|} e^{i\theta/2} &=& \sqrt{|z|} \left(\cos{\frac{\theta}{2}} + i \sin{\frac{\theta}{2}}\right) \\ &=& \sqrt{|z|} \left(\sqrt{\frac{1 + \cos \theta}{2}}) + i \sqrt{\frac{1 - \cos \theta}{2}}\right) \\ &=& \dfrac{\sqrt{|z| + |z|\cos \theta}}{2} + i \sqrt{\frac{|z| - |z| \cos \theta}{2}} \\ &=& \sqrt{\frac{|z| + x}{2}} + i \sqrt{\frac{|z| - x}{2}}, \quad (|z|\cos \theta = x) \end{eqnarray}\)

\(y \lt 0\) の場合、虚部の符号が \(-\) になる
\( \sqrt{|z|} e^{i\theta/2} = \sqrt{\dfrac{|z| + x}{2}} - i \sqrt{\dfrac{|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) { ... }

> Complex.Sqrt(one * Complex(2, 0));;
val it: Complex = complex(1.41421, 0) { ... }

> Complex.Sqrt(one * Complex(3, 0));;
val it: Complex = complex(1.73205, 0) { ... }

> Complex.Sqrt(one * Complex(-1, 0));;
val it: Complex = complex(0, 1) { ... }

> Complex.Sqrt(one * Complex(-2, 0));;
val it: Complex = complex(0, 1.41421) { ... }

> Complex.Sqrt(one * Complex(-3, 0));;
val it: Complex = complex(0, 1.73205) { ... }

> 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 x\) は \(\log x\) と同じ分枝切断を持っています。x を負の整数とすると sqrt(x) の解は \(i \sqrt x\) になりますが、もうひとつ \(-i \sqrt x\) という解があります。 

> Complex.Sqrt(Complex(-2.0, 0.0));;
val it: Complex = complex(0, 1.41421) { ... }

> Complex.Sqrt(Complex(-2.0, -0.0));;
val it: Complex = complex(0, -1.41421) { ... }

> Complex.Sqrt(Complex(-1e300, 0.0));;
val it: Complex = complex(0, 1e+150) { ... }

> Complex.Sqrt(Complex(-1e300, -0.0));;
val it: Complex = complex(0, -1e+150) { ... }

●逆三角関数

三角関数の逆関数を「逆三角関数 (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 は無数に存在します。つまり、逆三角関数の返り値は無数にあることになりますが、通常は一つの値を返すように範囲を制限します。これを「主値」といいます。

逆三角関数の主値を以下に示します。

\(\begin{array}{lcc} \arcsin x = w & -1 \leq x \lt 1 & -\pi/2 \leq w \leq \pi/2 \\ \arccos x = w & -1 \leq x \lt 1 & 0 \leq w \leq \pi \\ \arctan x = w & x \in \mathbb{R} & -\pi/2 \leq w \leq \pi/2 \end{array}\)

本ページでは、数式の表示に JavaScript のライブラリ MathJax を使っています。MathJax では、逆三角関数を arcsin, arccos, arctan と表示します。

簡単な実行例を示します。 

> 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 = ()

●複素数の逆三角関数

複素数の逆三角関数の定義は、複素数の三角関数の定義から導くことができます。\(\arcsin z\) の定義は次のようになります。

\(z, w\) は複素数とする

\(\begin{array}{l} \arcsin z = w \\ z = \sin w = \dfrac{e^{iw} - e^{-iw}}{2i} \\ z \times 2ie^{iw} = \dfrac{e^{iw} - e^{-iw}}{2i} \times 2ie^{iw} \\ 2iz(e^{iw}) = (e^{iw})^2 - 1 \\ (e^{iw})^2 - 2iz(e^{iw}) - 1 = 0 \end{array}\)

\(e^{iw}\) の二次方程式と考えて解くと

\(e^{iw} = iz \pm \sqrt{1 - z^2}\)

両辺の対数をとって \(-i\) を掛け算し、平方根の主値 \(+\) を選ぶ

\(w = \arcsin z = -i \log{\left(iz + \sqrt{1 - z^2}\right)} \)

\(\arccos z, \arctan z\) は定義だけを示します。

\(\begin{array}{l} \arccos z = \dfrac{\pi}{2} - \arcsin z \\ \arctan z = i \dfrac{\log{\left(1 - iz\right)} - \log{\left(1 + iz\right)}}{2} \end{array}\)

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) { ... }

> Complex.Asin(Complex(2.0, -0.0));;
val it: Complex = complex(1.5708, -1.31696) { ... }

> Complex.Asin(Complex(-4.0, 0.0));;
val it: Complex = complex(-1.5708, 2.06344) { ... }

> Complex.Asin(Complex(-4.0, -0.0));;
val it: Complex = complex(-1.5708, -2.06344) { ... }

> 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) { ... }

> Complex.Acos(Complex(2.0, -0.0));;
val it: Complex = complex(0, 1.31696) { ... }

> Complex.Acos(Complex(-4.0, 0.0));;
val it: Complex = complex(3.14159, -2.06344) { ... }

> Complex.Acos(Complex(-4.0, -0.0));;
val it: Complex = complex(3.14159, 2.06344) { ... }

> 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) { ... }

> Complex.Atan(Complex(-0.0, 2.0));;
val it: Complex = complex(-1.5708, 0.549306) { ... }

> Complex.Atan(Complex(0.0, -4.0));;
val it: Complex = complex(1.5708, -0.255413) { ... }

> Complex.Atan(Complex(-0.0, -4.0));;
val it: Complex = complex(-1.5708, -0.255413) { ... }
-- note --------
[*2] let y = i * z を計算するとき、z が 0.0 + b i と -0.0 + b i では同じ値になるため。 
(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.0i
IEEE 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 があります。MathJax では逆双曲線関数を \(\sinh^{-1} x, \cosh^{-1} x, \tanh^{-1} x \) で表します。双曲線関数と逆双曲線関数の定義域と値域を示します。

\(\begin{array}{lcc} y = \sinh x & -\infty \lt x \lt \infty & -\infty \lt y \lt \infty \\ y = \cosh x & -\infty \lt x \lt \infty & 1 \leq y \lt \infty \\ y = \tanh x & -\infty \lt x \lt \infty & -1 \lt y \lt 1 \\ y = \sinh^{-1} x & -\infty \lt x \lt \infty & -\infty \lt y \lt \infty \\ y = \cosh^{-1} x & 1 \leq x \lt \infty & 0 \leq y \lt \infty \\ y = \tanh^{-1} x & -1 \lt x \lt 1 & -\infty \lt y \lt \infty \end{array}\)

x と y は実数です。x = cosh y の逆関数 y = acosh x を満たす y の値は 2 つありますが、ここでは \(y \geq 0\) を主値として選ぶことにします。

逆双曲線関数の定義は双曲線関数の定義から導くことができます。

\(\begin{array}{l} \sinh^{-1} x = y \\ y = \sinh x = \dfrac{e^{x} - e^{-x}}{2} \\ y \times 2e^{x} = \dfrac{e^{x} - e^{-x}}{2} \times 2e^{x} \\ 2ye^{x} = e^{2x} - 1 \\ e^{2x} - 2ye^{x} - 1 = 0 \end{array}\)

\(e^{x}\) の 2 次方程式として解く

\(e^{x} = \dfrac{2y \pm \sqrt{4y^2 + 4}}{2} = y \pm \sqrt{y^2 + 1} \)

\(e^{x} \gt 0\) だから平方根の符号に \(+\) を選んで対数をとる

\(x = \log \left(y + \sqrt{y^2 + 1}\right)\)
\(\begin{array}{l} \cosh^{-1} x = y \\ y = \cosh x = \dfrac{e^{x} + e^{-x}}{2} \\ y \times 2e^{x} = \dfrac{e^{x} + e^{-x}}{2} \times 2e^{x} \\ 2ye^{x} = e^{2x} + 1 \\ e^{2x} - 2ye^{x} + 1 = 0 \end{array}\)

\(e^{x}\) の 2 次方程式として解く

\(e^{x} = \dfrac{2y \pm \sqrt{\left(4y^2 - 4\right)}}{2} = y \pm \sqrt{y^2 - 1} \)

\(e^{x} \gt 0\) だから平方根の符号に \(+\) を選んで対数をとる

\(x = \log \left(y + \sqrt{y^2 - 1}\right), \quad (y \geq 1) \)
\(\begin{array}{l} \tanh^{-1} x = y \\ \begin{eqnarray} y &=& \tanh x \\ &=& \dfrac{\sinh x}{\cosh x} \\ &=& \dfrac{e^{x} - e^{-x}}{e^{x} + e^{-x}} \\ &=& \dfrac{e^{2x} - 1}{e^{2x} + 1} \end{eqnarray} \\ y(e^{2x} + 1) = (e^{2x} - 1) \\ (1 - y)e^{2x} = 1 + y \\ e^{2x} = \dfrac{1 + y}{1 - y} \\ x = \dfrac{1}{2} \log{\dfrac{1 + y}{1 - y}} \end{array}\)

関数 \(\tanh^{-1} x\) を使うと、次の関係式から \(\sinh^{-1} x, \ \cosh^{-1} x\) を求めることができます。

\(\sinh^{-1} x = \tanh^{-1}{\dfrac{x}{\sqrt{x^2 + 1}}} \)

\(\cosh^{-1} x = \tanh^{-1}{\dfrac{\sqrt{1 - x^2}}{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 に変えたものになります。

\(\begin{array}{l} \sinh^{-1} z = \log \left(z + \sqrt{z^2 + 1}\right) \\ \begin{eqnarray} \cosh^{-1} z &=& \log \left(z + \sqrt{z^2 - 1}\right) &(1)& \\ &=& \log \left(z + \sqrt{z + 1} \sqrt{z - 1}\right) \quad &(2)& \end{eqnarray} \\ \tanh^{-1} z = \dfrac{1}{2} \log\left(\dfrac{1 + z}{1 - z}\right) = \dfrac{1}{2} \left( \log {\left( 1 + z \right)} - \log {\left(1 - z\right)} \right) \end{array}\)

参考 URL 『逆双曲線関数と逆三角関数の branch cut | 雑記帳』によると、\(\cosh^{-1} z\) を式 (1) でプログラムすると「分枝切断線」が複雑になるため、他の式 (たとえば (2) など) でプログラムする処理系が多いようです。ちなみに、ANSI Common Lisp では \(\cosh^{-1} z\) を次の式で定義しています。

\( \cosh^{-1} z = 2 \log \left( \sqrt{\dfrac{z + 1}{2}} + \sqrt{\dfrac{z - 1}{2}} \right) \)

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) { ... }

> Complex.Asinh(Complex(-0.0, 2.0));;
val it: Complex = complex(-1.31696, 1.5708) { ... }

> Complex.Asinh(Complex(0.0, -4.0));;
val it: Complex = complex(2.06344, -1.5708) { ... }

> Complex.Asinh(Complex(-0.0, -4.0));;
val it: Complex = complex(-2.06344, -1.5708) { ... }

> 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) { ... }

> Complex.Acosh(Complex(0.0, -0.0));;
val it: Complex = complex(0, -1.5708) { ... }

> Complex.Acosh(Complex(-4.0, 0.0));;
val it: Complex = complex(2.06344, 3.14159) { ... }

> Complex.Acosh(Complex(-4.0, -0.0));;
val it: Complex = complex(2.06344, -3.14159) { ... }

> 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) { ... }

> Complex.Atanh(Complex(2.0, -0.0));;
val it: Complex = complex(0.549306, -1.5708) { ... }

> Complex.Atanh(Complex(-4.0, 0.0));;
val it: Complex = complex(-0.255413, 1.5708) { ... }

> Complex.Atanh(Complex(-4.0, -0.0));;
val it: Complex = complex(-0.255413, -1.5708) { ... }

●参考文献, URL

  1. Guy L. Steele Jr., 『COMMON LISP 第 2 版』, 共立出版, 1991
  2. 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
  3. IEEE 754 -- Wikipedia
  4. IEEE 754における負のゼロ - Wikipedia
  5. NaN - Wikipedia
  6. 逆三角関数 - Wikipedia
  7. 逆双曲線関数 - Wikipedia
  8. 逆双曲線関数と逆三角関数の branch cut | 雑記帳

●プログラムリスト

//
// 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
初版 2022 年 5 月 7 日