M.Hiroi's Home Page

お気楽 Common Lisp プログラミング入門

入門編

Copyright (C) 2021 Makoto Hiroi
All rights reserved.

複素数

近年、多くのプログラミング言語で「複素数 (complex number)」がサポートされるようになりました。たとえば、Ｃ言語では 1999 年に発行された規格 C99 で複素数型が導入されました。Go 言語や Python でも複素数型をサポートしていますし、複素数用の標準ライブラリを用意している言語 (C++, Ruby, Haskell など) も多くあります。

他の言語では、FORTRAN や Common Lisp が昔から複素数型をサポートしています。Common Lisp の場合、基本的な数学関数でも複素数を適用できるのであれば、引数に複素数を渡して計算することができます。今回は Common Lisp の複素数についてまとめてみました。

●Common Lisp の数

Common Lisp の数は大きく分けると、整数 (integer)、分数 (ratio)、実数 (float)、複素数 (complex) の 4 種類があります。Common Lisp の整数はメモリの許す限り任意の精度で扱うことができます。分数は分子と分母を整数で保持しているので、整数と同様に任意の精度で扱うことができます。

実数は浮動小数点数 (floating point number) として表現されます。浮動小数点数には IEEE 754 という標準仕様があり、近代的なプログラミング言語のほとんどは、IEEE 754 に準拠した浮動小数点数をサポートしています。浮動小数点数はすべての小数を正確に表現することはできません。このため、実数は近似的な値になります。

IEEE 754 には通常の数値以外にも、負のゼロ (-0.0)、正負の無限大 (∞, -∞)、NaN (Not a Number, 非数) といった値が定義されています。これらの値は Common Lisp の仕様 (ANSI Common Lisp) に規定されていませんが、処理系によってはサポートしているものがあります。

●無限大

一般に、無限大は値のオーバーフロー、ゼロ除算 (数値 / 0.0)、数学関数の計算結果 (たとえば log(0.0)) などで発生します。ANSI Common Lisp に無限大は規定されていないので、浮動小数点数のゼロ除算や log(0.0) はエラー (例外) を送出するのが普通です。ただし、処理系によっては無限大を表す値を用意しているものがあります。

SBCL の場合、無限大は以下に示す定数に格納されています。

sb-ext:double-float-positive-infinity ; \(+\infty\)
sb-ext:double-float-negative-infinity ; \(-\infty\)

無限大は他の数値と比較したり演算することもできますが、結果が非数 (NaN) になることもあります。ただし、ANSI Common Lisp は NaN を規定していないので、動作は処理系に依存します。 

* (defconstant inf sb-ext:double-float-positive-infinity)

INF
* (defconstant -inf sb-ext:double-float-negative-infinity)

-INF
* (= inf inf)

T
* (= inf -inf)

NIL
* (/= inf -inf)

T
* (< -inf inf)

T
* (< inf -inf)

NIL
* (< 0 inf)

T
* (< -inf 0)

T
* (+ inf 10)

#.SB-EXT:DOUBLE-FLOAT-POSITIVE-INFINITY
* (- inf 10)

#.SB-EXT:DOUBLE-FLOAT-POSITIVE-INFINITY
* (* inf 10)

#.SB-EXT:DOUBLE-FLOAT-POSITIVE-INFINITY
* (/ inf 10)

#.SB-EXT:DOUBLE-FLOAT-POSITIVE-INFINITY
* (* inf inf)

#.SB-EXT:DOUBLE-FLOAT-POSITIVE-INFINITY
* (+ inf inf)

#.SB-EXT:DOUBLE-FLOAT-POSITIVE-INFINITY

* (- inf inf)
=> 処理系が落ちる

* (+ inf -inf)
=> 処理系が落ちる

* (* inf 0.0)
=> 処理系が落ちる

一般に、inf - inf, inf + -inf, inf * 0.0 は NaN になるのですが、SBCL (version 1.4.5) では処理系が終了してしまいます。SBCL では NaN をサポートしていないのかもしれません。

●負のゼロ

負のゼロ (-0.0) は、計算結果が負の極めて小さな値でアンダーフローになったとき発生します。また、正の値を負の無限大で除算する、負の値を正の無限大で除算する、負の値と 0.0 を乗算しても -0.0 が得られます。 

* 0.0

0.0
* -0.0

-0.0
* 0d0

0.0d0
* -0d0

-0.0d0
* -1d-323

-9.881312916824931d-324
* -1d-324

-0.0d0
* (/ -1d-323 2)

-4.9406564584124654d-324
* (/ -1d-323 4)

-0.0d0
* (/ 1d0 -inf)

-0.0d0
* (/ -1d0 inf)

-0.0d0

IEEE 754 では、演算子 (Common Lisp では =) による 0.0 と -0.0 の比較は等しいと判定されます。なお SBCL の場合、(eql 0.0 -0.0) は偽を返します。 

* (= 0.0 -0.0)

T
* (/= 0.0 -0.0)

NIL
* (= 0d0 -0d0)

T
* (/= 0d0 -0d0)

NIL
* (eql 0.0 -0.0)

NIL
* (eql -0.0 -0.0)

T
* (eql 0d0 -0d0)

NIL
* (eql -0d0 -0d0)

T

なお、-0.0 は数学関数 (関数 atan など) や複素数の演算処理などで使われます。 

* (sqrt 0d0)

0.0d0
* (sqrt -0d0)

-0.0d0
* (atan 0d0 -1d0)

3.141592653589793d0
* (atan -0d0 -1d0)

-3.141592653589793d0

●Common Lisp の複素数

数学では複素数 z を \(x + yi\) と表記します。x を実部、y を虚部、i を虚数単位といいます。虚数単位は 2 乗すると -1 になる数です。実部と虚部の 2 つの数値を格納するデータ構造を用意すれば、プログラミング言語でも複素数を表すことができます。

Common Lisp では複素数 \(z = x + yi\) を ＃C(x y) と表記します。複素数のデータ型は complex で、型述語は complexp です。実部と虚部の指定には、整数、分数、実数を使うことができます。実部と虚部が異なる種類 (型) の数ならば、同じ種類の数になるよう変換されます。計算結果で実部と虚部が整数または分数になった場合、虚部が 0 だと実部のデータ型 (整数または分数) に変換されます。

複素数は関数 complex で生成することもできます。 

complex a &optional b

a が実部で、b が虚部です。a が整数または分数のとき、虚部 b を省略すると a のデータ型になります。複素数 z の実部は関数 realpart で、虚部は関数 imagpart で取得することができます。複素数の虚部の符号を反転することを「複素共役」といいます。複素共役は関数 conjugate で求めることができます。

簡単な例を示しましょう。 

* #c(1 2)

#C(1 2)
* (complex 1/2 1/3)

#C(1/2 1/3)
* (complex 1 0)

1
* (complex 1/2 0)

1/2
* (complex 1.0 0)

#C(1.0 0.0)
* (defvar a #c(1.234 5.678))

A
* a

#C(1.234 5.678)
* (realpart a)

1.234
* (imagpart a)

5.678
* (conjugate a)

#C(1.234 -5.678)
* (complexp a)

T
* (complexp 1.2345)

NIL

複素数は極形式 \(z = r (\cos \theta + i \sin \theta)\) で表すことができます。このとき、\(r\) を絶対値、\(\theta\) を偏角といいます。絶対値は関数 abs で求めることができます。偏角は複素平面において正の実軸とベクトル (x, y) との角度を表します。偏角を求めるには関数 phase を使います。phase の返り値 \(\theta\) は、負のゼロをサポートしている処理系では \(-\pi \leq \theta \leq \pi\) に、サポートしていない処理系では \(-\pi \lt \theta \leq \pi\) になります。

複素数 z はオイラーの公式 \(e^{i\theta} = \cos \theta + i \sin \theta\) を使って \(z = re^{i\theta}\) と書くこともできます。\(\theta\) に pi を代入するとオイラーの等式 \(e^{i \pi} + 1 = 0\) になります。ド・モアヴルの公式も導出できます。

\( z^n = (re^{i\theta})^n = r^ne^{in\theta} = r^n(\cos n\theta + i \sin n\theta) \)

関数 cis x は絶対値が 1 で偏角が x の複素数を生成します。cis は cos + i sin の略です。関数 signum z に複素数を与えると、絶対値が 1 で偏角が z と等しい複素数を返します。基本的な数学関数 (sqrt, exp, log, 三角関数など) は複素数でも使用できます。

簡単な例を示しましょう。 

* (defvar b (complex 1d0 1d0))

B
* b

#C(1.0d0 1.0d0)
* (abs b)

1.4142135623730951d0
* (phase b)

0.7853981633974483d0
* (/ pi 4)

0.7853981633974483d0
* (cis (/ pi 4))

#C(0.7071067811865476d0 0.7071067811865475d0)
* (signum b)

#C(0.7071067811865475d0 0.7071067811865475d0)


* (phase #c(1d0 0d0))

0.0d0
* (phase #c(1d0 1d0))

0.7853981633974483d0
* (phase #c(0d0 1d0))

1.5707963267948966d0
* (phase #c(-1d0 1d0))

2.356194490192345d0
* (phase #c(-1d0 0d0))

3.141592653589793d0
* (phase #c(1d0 -1d0))

-0.7853981633974483d0
* (phase #c(0d0 -1d0))

-1.5707963267948966d0
* (phase #c(-1d0 -1d0))

-2.356194490192345d0
* (phase #c(-1d0 -0d0))

-3.141592653589793d0

SBCL の場合、\(-1.0+0.0i\) の偏角 \(\theta\) は \(\pi\) になり、\(-1.0-0.0i\) の偏角は \(-\pi\) になります。ゼロと負のゼロを区別しない処理系では、偏角 \(\theta\) の範囲を \(-\pi \lt \theta \leq \pi\) に制限して、\(-1.0 + 0.0i \ (= -1.0 - 0.0i)\) の偏角を \(\pi\) とします。

●複素数の四則演算

複素数の四則演算は次のようになります。

\( \begin{array}{l} (a + bi) + (c + di) = (a + c) + (b + d)i \\ (a + bi) - (c + di) = (a - c) + (b - d)i \\ (a + bi) \times (c + di) = (ac - bd) + (bc + ad)i \\ (a + bi) \div (c + di) = \dfrac{ac + bd + (bc - ad)i}{c^2 + d^2} \end{array} \)

これらの演算は Common Lisp の関数 +, -, * , / で行うことができます。複素数の場合、大小の比較演算は使えませんが、等値の判定は関数 =, /= で行うことができます。eql を使うと 0.0 と -0.0 を区別することができます。簡単な例を示しましょう。 

* (defvar a #c(1.0 2.0))

A
* a

#C(1.0 2.0)
* (defvar b #c(3.0 4.0))

B
* b

#C(3.0 4.0)
* (+ a b)

#C(4.0 6.0)
* (- a b)

#C(-2.0 -2.0)
* (* a b)

#C(-5.0 10.0)
* (/ a b)

#C(0.44 0.08)
* (= #c(1.0 -0.0) #c(1.0 0.0))

T
* (eql #c(1.0 -0.0) #c(1.0 0.0))

NIL
* (eql #c(1.0 0.0) #c(1.0 0.0))

T

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

複素数を引数にとる指数関数はオイラー (Euler) の公式から導くことができます。

\( \begin{array}{l} e^{i\theta} = \cos \theta + i \sin \theta \\ e^{x+iy} = e^x e^{iy} = e^x (\cos y + i \sin y) \end{array} \)

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

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

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

\(x^y = e^{y\log x}\)

Common Lisp の関数 exp と log は複素数に対応しています。べき乗も関数 expt で求めることができます。簡単な例を示しましょう。 

* pi

3.141592653589793d0
* (exp #c(0d0 0d0))

#C(1.0d0 0.0d0)
* (exp (complex 0d0 (/ pi 4)))

#C(0.7071067811865476d0 0.7071067811865475d0)
* (exp (complex 0d0 (/ pi 2)))

#C(6.123233995736766d-17 1.0d0)
* (exp (complex 0d0 pi))

#C(-1.0d0 1.2246467991473532d-16)
* (exp (complex 0d0 (- pi)))

#C(-1.0d0 -1.2246467991473532d-16)
* (exp (complex 1d0 1d0))

#C(1.4686939399158851d0 2.2873552871788423d0)


* (log #c(1d0 1d0))

#C(0.34657359027997264d0 0.7853981633974483d0)
* (log #c(1d0 0d0))

#C(0.0d0 0.0d0)
* (log #c(0d0 1d0))

#C(0.0d0 1.5707963267948966d0)
* (log #c(1d0 -1d0))

#C(0.34657359027997264d0 -0.7853981633974483d0)
* (log #c(1d300 1d300))

#C(691.1221014884936d0 0.7853981633974483d0)


* (expt #c(1d0 1d0) 0d0)

#C(1.0d0 0.0d0)
* (expt #c(1d0 1d0) 1d0)

#C(1.0d0 0.9999999999999998d0)
* (expt #c(1d0 1d0) 2d0)

#C(1.2246467991473532d-16 2.0d0)
* (expt #c(1d0 1d0) 3d0)

#C(-1.9999999999999996d0 2.0d0)
* (expt #c(1d0 2d0) #c(3d0 4d0))

#C(0.12900959407446694d0 0.03392409290517001d0)
* (expt #c(1d0 1d0) #c(1d0 1d0))

#C(0.2739572538301211d0 0.5837007587586147d0)

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

* (log #c(-1d0 0d0))

#C(0.0d0 3.141592653589793d0)
* (log #c(-1d0 -0d0))

#C(0.0d0 -3.141592653589793d0)
* (log #c(-1d300 0d0))

#C(690.7755278982137d0 3.141592653589793d0)
* (log #c(-1d300 -0d0))

#C(690.7755278982137d0 -3.141592653589793d0)

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

参考 URL 8 『逆双曲線関数と逆三角関数の branch cut | 雑記帳』によると、この取り除いた領域を branch cut と呼ぶそうです。プログラミングでは branch cut を定義域から取り除くのではなく、その領域では不連続な関数とするそうです。参考文献 1 『COMMON LISP 第 2 版』では「分枝切断線」、Python のドキュメントでは「分枝切断」と記述されています。本稿では branch cut を「分枝切断」と記述することにします。

プログラミング言語の場合、0.0 と -0.0 を区別する処理系であれば、SBCL のように 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\)) になります。

双曲線関数の定義
\( \sinh x = \dfrac{e^x - e^{-x}}{2}, \quad \cosh x = \dfrac{e^x + e^{-x}}{2} \)
\(\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}\)

Common Lisp の三角関数 (sin, cos, tan) は複素数にも対応しています。簡単な実行例を示します。 

* (sin #c(0d0 1d0))

#C(0.0d0 1.1752011936438014d0)
* (sin #c(0d0 -1d0))

#C(0.0d0 -1.1752011936438014d0)
* (sin #c(1d0 -1d0))

#C(1.2984575814159773d0 -0.6349639147847361d0)
* (sin #c(1d0 1d0))

#C(1.2984575814159773d0 0.6349639147847361d0)
* (abs (sin #c(1d0 1d0)))

1.4453965766582495d0


* (cos #c(0d0 1d0))

#C(1.5430806348152437d0 -0.0d0)
* (cos #c(0d0 -1d0))

#C(1.5430806348152437d0 0.0d0)
* (cos #c(1d0 -1d0))

#C(0.8337300251311491d0 0.9888977057628651d0)
* (cos #c(1d0 1d0))

#C(0.8337300251311491d0 -0.9888977057628651d0)
* (abs (cos #c(1d0 1d0)))

1.2934544550420957d0


* (tan #c(0d0 1d0))

#C(0.0d0 0.7615941559557649d0)
* (tan #c(0d0 -1d0))

#C(0.0d0 -0.7615941559557649d0)
* (tan #c(1d0 -1d0))

#C(0.27175258531951174d0 -1.0839233273386946d0)
* (tan #c(1d0 1d0))

#C(0.27175258531951174d0 1.0839233273386946d0)
* (abs (tan #c(1d0 1d0)))

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

●複素数の双曲線関数

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

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

\(\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}\)

Common Lisp の双曲線関数 (\(\sinh, \cosh, \tanh\)) は複素数にも対応しています。簡単な使用例を示します。 

* (sinh #c(0d0 1d0))

#C(0.0d0 0.8414709848078965d0)
* (sinh #c(0d0 -1d0))

#C(0.0d0 -0.8414709848078965d0)
* (sinh #c(1d0 1d0))

#C(0.6349639147847361d0 1.2984575814159773d0)
* (sinh #c(1d0 0d0))

#C(1.1752011936438014d0 0.0d0)
* (sinh #c(0d0 0d0))

#C(0.0d0 0.0d0)


* (cosh #c(0d0 1d0))

#C(0.5403023058681398d0 0.0d0)
* (cosh #c(0d0 -1d0))

#C(0.5403023058681398d0 -0.0d0)
* (cosh #c(1d0 1d0))

#C(0.8337300251311491d0 0.9888977057628651d0)
* (cosh #c(1d0 0d0))

#C(1.5430806348152437d0 0.0d0)
* (cosh #c(0d0 0d0))

#C(1.0d0 0.0d0)


* (tanh #c(0d0 1d0))

#C(0.0d0 1.5574077246549023d0)
* (tanh #c(0d0 -1d0))

#C(0.0d0 -1.5574077246549023d0)
* (tanh #c(1d0 1d0))

#C(1.0839233273386946d0 0.27175258531951174d0)
* (tanh #c(1d0 0d0))

#C(0.7615941559557649d0 0.0d0)
* (tanh #c(0d0 0d0))

#C(0.0d0 0.0d0)

●複素数の平方根

複素数 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}} \)

Common Lisp の関数 sqrt は複素数にも対応しています。簡単な実行例を示します。 

* (sqrt #c(1d0 0d0))

#C(1.0d0 0.0d0)
* (sqrt #c(2d0 0d0))

#C(1.4142135623730951d0 0.0d0)
* (sqrt #c(3d0 0d0))

#C(1.7320508075688772d0 0.0d0)
* (sqrt #c(-1d0 0d0))

#C(0.0d0 1.0d0)
* (sqrt #c(-2d0 0d0))

#C(0.0d0 1.4142135623730951d0)
* (sqrt #c(-3d0 0d0))

#C(0.0d0 1.7320508075688772d0)


* (defvar x (sqrt #c(1d0 1d0)))

X
* x

#C(1.09868411346781d0 0.45508986056222733d0)
* (* x x)

#C(1.0000000000000002d0 1.0d0)
* (setq x (sqrt #c(1d0 -1d0)))

#C(1.09868411346781d0 -0.45508986056222733d0)
* (* x x)

#C(1.0000000000000002d0 -1.0d0)
* (setq x (sqrt #c(-1d0 1d0)))

#C(0.45508986056222733d0 1.09868411346781d0)
* (* x x)

#C(-1.0000000000000002d0 1.0d0)
* (setq x (sqrt #c(-1d0 -1d0)))

#C(0.45508986056222733d0 -1.09868411346781d0)
* (* x x)

#C(-1.0000000000000002d0 -1.0d0)

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

* (sqrt #c(-1d0 0d0))

#C(0.0d0 1.0d0)
* (sqrt #c(-1d0 -0d0))

#C(0.0d0 -1.0d0)
* (sqrt #c(-2d0 0d0))

#C(0.0d0 1.4142135623730951d0)
* (sqrt #c(-2d0 -0d0))

#C(0.0d0 -1.4142135623730951d0)
* (sqrt #c(-1d300 0d0))

#C(0.0d0 1.0d150)
* (sqrt #c(-1d300 -0d0))

#C(0.0d0 -1.0d150)

●逆三角関数

三角関数の逆関数を「逆三角関数 (inverse trigonometric function)」といいます。Common Lisp に用意されている逆三角関数を示します。

ここでは引数 x を実数とします。asin x は引数 x が与えられたとき sin w = x となる角度 w を求めます。同様に acos x は cos w = x となる角度 w を、atan x は tan x = w となる角度 w を求めます。数式を表示するライブラリ MathJax では、\(\arcsin x, \arccos x, \arctan x\) と表示されます。

三角関数には周期性があるので、上式を満たす角度 w は無数に存在します。つまり、逆三角関数の返り値は無数にあることになりますが、通常は一つの値を返すように範囲を制限します。これを「主値」といいます。逆三角関数の主値を以下に示します。

\(\begin{array}{lll} \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}\)

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

* (dolist (x '(-1d0 -0.5d0 0d0 0.5d0 1d0)) (print (asin x)))

-1.5707963267948966d0
-0.5235987755982989d0
0.0d0
0.5235987755982989d0
1.5707963267948966d0
NIL
* (dolist (x '(-1d0 -0.5d0 0d0 0.5d0 1d0)) (print (acos x)))

3.141592653589793d0
2.0943951023931957d0
1.5707963267948966d0
1.0471975511965979d0
0.0d0
NIL
* (dolist (x '(-1d300 -1d0 0d0 1d0 1d300)) (print (atan x)))

-1.5707963267948966d0
-0.7853981633974483d0
0.0d0
0.7853981633974483d0
1.5707963267948966d0
NIL

関数 atan は引数を 2 つ受け取ることができます。これは他のプログラミング言語、たとえばＣ言語の数学関数 atan2 と同じです。 

atan y &optional x => 角度 (ラジアン)

引数 x, y は実数です。atan は直交座標系においてベクトル (x, y) と x 軸との角度を求める関数です。複素平面で考えると、複素数 \(x + yi\) の偏角 \(\theta\) を求めることと同じです。SBCL は負のゼロをサポートしているので、返り値 (角度 \(\theta\)) の範囲は \(-\pi \leq \theta \leq \pi\) になります。

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

* (dolist (a '((0d0 1d0) (1d0 1d0) (1d0 0d0) (1d0 -1d0) (0d0 -1d0))) (print (apply #'atan a)))

0.0d0
0.7853981633974483d0
1.5707963267948966d0
2.356194490192345d0
3.141592653589793d0
NIL
* (dolist (a '((0d0 1d0) (-1d0 1d0) (-1d0 0d0) (-1d0 -1d0) (-0d0 -1d0))) (print (apply #'atan a)))

0.0d0
-0.7853981633974483d0
-1.5707963267948966d0
-2.356194490192345d0
-3.141592653589793d0
NIL

●複素数の逆三角関数

複素数の逆三角関数の定義は、複素数の三角関数の定義から導くことができます。\(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}\)

\(\arcsin z, \arccos z, \arctan z\) は次に示す分枝切断を持っています。

Common Lisp の逆三角関数 (asin, acos, atan) は複素数に対応しています。簡単な実行例を示します。 

* (dolist (x '(-1d0 -0.5d0 0d0 0.5d0 1d0)) (print (asin (complex x 0d0))))

#C(-1.5707963267948966d0 0.0d0)
#C(-0.5235987755982989d0 0.0d0)
#C(0.0d0 0.0d0)
#C(0.5235987755982989d0 0.0d0)
#C(1.5707963267948966d0 0.0d0)
NIL
* (asin #c(0d0 1d0))

#C(0.0d0 0.881373587019543d0)
* (sin (asin #c(0d0 1d0)))

#C(0.0d0 1.0d0)
* (asin #c(1d0 1d0))

#C(0.6662394324925153d0 1.0612750619050357d0)
* (sin (asin #c(1d0 1d0)))

#C(1.0000000000000002d0 1.0d0)

* (asin #c(2d0 0d0))

#C(1.5707963267948966d0 -1.3169578969248166d0)
* (asin #c(2d0 -0d0))

#C(1.5707963267948966d0 -1.3169578969248166d0)  ; 0.0i と -0.0i で同じ値が返る
* (asin #c(-4d0 0d0))

#C(-1.5707963267948966d0 2.0634370688955608d0)
* (asin #c(-4d0 -0d0))

#C(-1.5707963267948966d0 2.0634370688955608d0)  ; 0.0i と -0.0i で同じ値が返る


* (dolist (x '(-1d0 -0.5d0 0d0 0.5d0 1d0)) (print (acos (complex x 0d0))))

#C(3.141592653589793d0 0.0d0)
#C(2.0943951023931953d0 0.0d0)
#C(1.5707963267948966d0 0.0d0)
#C(1.0471975511965979d0 0.0d0)
#C(0.0d0 0.0d0)
NIL
* (acos #c(0d0 1d0))

#C(1.5707963267948966d0 -0.881373587019543d0)
* (cos (acos #c(0d0 1d0)))

#C(8.659560562354932d-17 1.0d0)
* (acos #c(1d0 1d0))

#C(0.9045568943023813d0 -1.0612750619050357d0)
* (cos (acos #c(1d0 1d0)))

#C(1.0000000000000002d0 1.0d0)


* (acos #c(2d0 0d0))

#C(0.0d0 1.3169578969248166d0)
* (acos #c(2d0 -0d0))

#C(0.0d0 1.3169578969248166d0)  ; 0.0i と -0.0i で同じ値が返る
* (acos #c(-4d0 0d0))

#C(3.141592653589793d0 -2.0634370688955608d0)
* (acos #c(-4d0 -0d0))

#C(3.141592653589793d0 -2.0634370688955608d0)  ; 0.0i と -0.0i で同じ値が返る


* (dolist (x '(-2d0 -1d0 -0.5d0 0.0 0.5d0 1d0 2d0)) (print (atan x)))

-1.1071487177940904d0
-0.7853981633974483d0
-0.4636476090008061d0
0.0
0.4636476090008061d0
0.7853981633974483d0
1.1071487177940904d0
NIL
* (atan #c(1d0 1d0))

#C(1.0172219678978514d0 0.40235947810852507d0)
* (tan (atan #c(1d0 1d0)))

#C(1.0d0 0.9999999999999999d0)
* (atan #c(1d0 -1d0))

#C(1.0172219678978514d0 -0.40235947810852507d0)
* (tan (atan #c(1d0 -1d0)))

#C(1.0d0 -0.9999999999999999d0)


* (atan #c(0d0 2d0))

#C(1.5707963267948966d0 0.5493061443340549d0)
* (atan #c(-0d0 2d0))

#C(-1.5707963267948966d0 0.5493061443340549d0)
* (atan #c(0d0 -4d0))

#C(1.5707963267948966d0 -0.25541281188299536d0)
* (atan #c(-0d0 -4d0))

#C(-1.5707963267948966d0 -0.25541281188299536d0)

SBCL (version 1.4.5) の場合、asin と acos の分枝切断の処理に問題があるようです。

●逆双曲線関数

双曲線関数の逆関数を「逆双曲線関数 (inverse hyperbolic function)」といいます。ANSI Common Lisp で用意されている逆双曲線関数には 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 = \cosh^{-1} 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}} \)

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

* (dolist (x '(-2d0 -1d0 0d0 1d0 2d0)) (print (asinh x)))

-1.4436354751788103d0
-0.881373587019543d0
0.0d0
0.881373587019543d0
1.4436354751788103d0
NIL
* (dolist (x '(-2d0 -1d0 0d0 1d0 2d0)) (print (sinh (asinh x))))

-1.9999999999999998d0
-1.0d0
0.0d0
1.0d0
1.9999999999999998d0
NIL

* (dolist (x '(1d0 1.5d0 2d0 2.5d0 3d0)) (print (acosh x)))

0.0d0
0.9624236501192069d0
1.3169578969248166d0
1.566799236972411d0
1.762747174039086d0
NIL
* (dolist (x '(1d0 1.5d0 2d0 2.5d0 3d0)) (print (cosh (acosh x))))

1.0d0
1.5d0
1.9999999999999998d0
2.5d0
3.0000000000000004d0
NIL

* (dolist (x '(-0.5d0 -0.25d0 0d0 0.25d0 0.5d0)) (print (atanh x)))

-0.5493061443340548d0
-0.25541281188299536d0
0.0d0
0.25541281188299536d0
0.5493061443340548d0
NIL
* (dolist (x '(-0.5d0 -0.25d0 0d0 0.25d0 0.5d0)) (print (tanh (atanh x))))

-0.49999999999999994d0
-0.25d0
0.0d0
0.25d0
0.49999999999999994d0
NIL

●複素数の逆双曲線関数

複素数の双曲線関数の定義は、実数の定義で引数 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) \quad \quad \quad (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 8 『逆双曲線関数と逆三角関数の 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) \)

\(\sinh^{-1} z, \cosh^{-1} z, \tanh^{-1} z\) は次に示す分枝切断を持っています。

Common Lisp の逆双曲線関数 (asinh, acosh, atanh) は複素数に対応しています。簡単な実行例を示します。 

* (dolist (x '(-2d0 -1d0 0d0 1d0 2d0)) (print (asinh (complex x 0d0))))

#C(-1.4436354751788103d0 0.0d0)
#C(-0.881373587019543d0 0.0d0)
#C(0.0d0 0.0d0)
#C(0.881373587019543d0 0.0d0)
#C(1.4436354751788103d0 0.0d0)
NIL
* (dolist (x '(-2d0 -1d0 0d0 1d0 2d0)) (print (sinh (asinh (complex x 0d0)))))

#C(-1.9999999999999998d0 0.0d0)
#C(-1.0d0 0.0d0)
#C(0.0d0 0.0d0)
#C(1.0d0 0.0d0)
#C(1.9999999999999998d0 0.0d0)
NIL
* (asinh #c(1d0 1d0))

#C(1.0612750619050357d0 0.6662394324925153d0)
* (sinh (asinh #c(1d0 1d0)))

#C(1.0d0 1.0000000000000002d0)


* (asinh #c(0d0 2d0))

#C(1.3169578969248166d0 1.5707963267948966d0)
* (asinh #c(-0d0 2d0))

#C(1.3169578969248166d0 1.5707963267948966d0)  ; 0.0i と -0.0i で同じ値が返る
* (asinh #c(0d0 -4d0))

#C(-2.0634370688955608d0 -1.5707963267948966d0)
* (asinh #c(-0d0 -4d0))

#C(-2.0634370688955608d0 -1.5707963267948966d0)  ; 0.0i と -0.0i で同じ値が返る


* (dolist (x '(1d0 1.5d0 2d0 2.5d0 3d0)) (print (acosh (complex x 0d0))))

#C(0.0d0 0.0d0)
#C(0.9624236501192069d0 0.0d0)
#C(1.3169578969248166d0 0.0d0)
#C(1.5667992369724109d0 0.0d0)
#C(1.762747174039086d0 0.0d0)
NIL
* (dolist (x '(1d0 1.5d0 2d0 2.5d0 3d0)) (print (cosh (acosh (complex x 0d0)))))

#C(1.0d0 0.0d0)
#C(1.5d0 0.0d0)
#C(1.9999999999999998d0 0.0d0)
#C(2.499999999999999d0 0.0d0)
#C(3.0000000000000004d0 0.0d0)
NIL
* (acosh #c(1d0 1d0))

#C(1.0612750619050357d0 0.9045568943023813d0)
* (cosh (acosh #c(1d0 1d0)))

#C(1.0000000000000002d0 1.0d0)

* (acosh #c(0d0 0d0))

#C(0.0d0 1.5707963267948966d0)


* (acosh #c(0d0 0d0))

#C(0.0d0 1.5707963267948966d0)
* (acosh #c(0d0 -0d0))

#C(0.0d0 -1.5707963267948966d0)
* (acosh #c(-4d0 0d0))

#C(2.0634370688955608d0 3.141592653589793d0)
* (acosh #c(-4d0 -0d0))

#C(-2.0634370688955608d0 -3.141592653589793d0)  ; 実部の符号が逆


* (dolist (x '(-0.5d0 -0.25d0 0d0 0.25d0 0.5d0)) (print (atanh (complex x 0d0))))

#C(-0.5493061443340549d0 0.0d0)
#C(-0.25541281188299536d0 0.0d0)
#C(0.0d0 0.0d0)
#C(0.25541281188299536d0 0.0d0)
#C(0.5493061443340549d0 0.0d0)
NIL
* (dolist (x '(-0.5d0 -0.25d0 0d0 0.25d0 0.5d0)) (print (tanh (atanh (complex x 0d0)))))

#C(-0.5000000000000001d0 0.0d0)
#C(-0.25d0 0.0d0)
#C(0.0d0 0.0d0)
#C(0.25d0 0.0d0)
#C(0.5000000000000001d0 0.0d0)
NIL
* (atanh #c(1d0 1d0))

#C(0.40235947810852507d0 1.0172219678978514d0)
* (tanh (atanh #c(1d0 1d0)))

#C(0.9999999999999999d0 1.0d0)


* (atanh #c(2d0 0d0))

#C(0.5493061443340549d0 1.5707963267948966d0)
* (atanh #c(2d0 -0d0))

#C(0.5493061443340549d0 -1.5707963267948966d0)
* (atanh #c(-4d0 0d0))

#C(-0.25541281188299536d0 1.5707963267948966d0)
* (atanh #c(-4d0 -0d0))

#C(-0.25541281188299536d0 -1.5707963267948966d0)

SBCL (version 1.4.5) の場合、asinh と acosh の分枝切断の処理に問題があるようです。

●参考文献, 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 | 雑記帳
初版 2021 年 11 月 6 日