M.Hiroi's Home Page

Common Lisp Programming

お気楽 ISLisp プログラミング超入門

Copyright (C) 2021 Makoto Hiroi
All rights reserved.

簡単なプログラム

●複素数

今回は ISLisp で「複素数 (complex number)」の演算プログラムを作ってみましょう。数学では複素数 z を \(x + iy\) と表記します。x を実部、y を虚部、i を虚数単位といいます。虚数単位は 2 乗すると -1 になる数です。実部と虚部の 2 つの数値を格納するデータ構造を用意すれば、プログラミング言語で複素数を表すことができます。今回は ISLisp のクラスを使って複素数を表すことにしましょう。 

リスト : 複素数型

;;; クラス定義
(defclass <complex> ()
  ((re :accessor realpart :initform 0.0 :initarg re)
   (im :accessor imagpart :initform 0.0 :initarg im)))

;;; コンストラクタ
(defun complex (x y)
  (create (class <complex>) 're (convert x <float>) 'im (convert y <float>)))

;;; 表示
(defun cprint (z)
  (format (standard-output) "#C(~G ~G)~%" (realpart z) (imagpart z)))

クラス名は <complex> としました。実部をスロット re に、虚部をスロット im にセットします。格納する数値の型は実数 (float) とします。実部はメソッド realpart で、虚部はメソッド imagpart で取得します。名前は Common Lisp から拝借しました。

複素数は関数 complex で生成します。引数 x が実部で、y が虚部です。x, y をスロット re, im にセットするとき、convert を使って float に変換しています。関数 cprint は標準出力に複素数 #C(x y) を出力します。表記は Common Lisp と同じです。

複素数 \(z = x + iy\) の虚部の符号を反転した数 \(x - iy\) を複素共役といいます。複素数を極形式 \(z = |z|(\cos \theta + i \sin \theta)\) で表した場合、\(|z|\) を絶対値、\(\theta\) を偏角といいます。絶対値 \(|z|\) の定義は \(\sqrt{x^2 + y^2}\) です。偏角は数学関数 atan2(y, x) で求めることができます。これをプログラムすると次のようになります。 

リスト : 複素共役, 絶対値, 偏角

;;; 複素共役
(defun conj (z)
  (complex (realpart z) (- (imagpart z))))

;;; 絶対値
(defun cabs (z)
  (cond
   ((= (realpart z) 0.0) (abs (imagpart z)))
   ((= (imagpart z) 0.0) (abs (realpart z)))
   ((> (abs (imagpart z)) (abs (realpart z)))
    (let ((temp (quotient (realpart z) (imagpart z))))
      (* (abs (imagpart z)) (sqrt (+ 1 (* temp temp))))))
   (t
    (let ((temp (quotient (imagpart z) (realpart z))))
      (* (abs (realpart z)) (sqrt (+ 1 (* temp temp))))))))

;;; 偏角
(defun carg (z) (atan2 (imagpart z) (realpart z)))

複素共役を求める関数 conj と偏角を求める関数 carg は簡単ですね。(atan2 y x) は直交座標においてベクトル (x, y) と x 軸との角度を求める関数です。角度 \(\theta\) の範囲は \(-\pi \lt \theta \leq \pi\) になります。簡単な例を示しましょう。 

ISLisp>(atan2 0 1)
0.0
ISLisp>(atan2 1 1)
0.785398
ISLisp>(atan2 1 0)
1.5708
ISLisp>(atan2 0 -1)
3.14159
ISLisp>(atan2 -1 1)
-0.785398
ISLisp>(atan2 -1 0)
-1.5708
ISLisp>(atan2 -1 -1)
-2.35619
ISLisp>(atan2 -0.1 -1)
-3.04192
ISLisp>(atan2 -0.01 -1)
-3.13159
ISLisp>(atan2 -0.001 -1)
-3.14059

\(y \geq 0\) の場合、atan2 の返り値は \(0 \leq \theta \leq \pi\) になり、\(y \lt 0\) の場合は \(-\pi \lt \theta \lt 0\) になります。

絶対値を求める関数 cabs は、定義をそのままプログラムすると二乗の計算でオーバーフローすることがあります。たとえば、#C(1e300 1e300) の絶対値を求めてみましょう。このとき、1e300 の二乗でオーバーフローします。 

ISLisp>(* 1e300 1e300)
> Error at *
> Floating point overflow: #<SFUNCTION 00000D32: *>, (1e+300 1e+300)

絶対値を求める場合、定義をそのままプログラムすると二乗の計算でオーバーフローすることがあります。参考文献『Ｃ言語による最新アルゴリズム事典』によると、次のように場合分けすることで、\(x^2\) や \(y^2\) で生じ得る上位桁あふれを回避することができるそうです。

\( \sqrt{x^2 + y^2} = \begin{cases} |x| \sqrt{1 + \left(\dfrac{y}{x}\right)^2}\ \quad & \mathrm{if} \ |x| \geq |y| \\ |y| \sqrt{1 + \left(\dfrac{x}{y}\right)^2}\ \quad & \mathrm{other} \end{cases} \)

今回は参考文献 1 のプログラムを ISLisp に移植しました。それでは実際に試してみましょう。 

ISLisp>(cprint (complex 1 1))
#C(1.0 1.0)
NIL
ISLisp>(cprint (conj (complex 1 1)))
#C(1.0 -1.0)
NIL
ISLisp>(cprint (conj (conj (complex 1 1))))
#C(1.0 1.0)
NIL
ISLisp>(cabs (complex 1 1))
1.41421
ISLisp>(cabs (complex 1 -1))
1.41421
ISLisp>(cabs (complex 1e300 1e300))
1.41421e+300
ISLisp>(cabs (complex 1e301 1e300))
1.00499e+301
ISLisp>(cabs (complex 1e300 1e301))
1.00499e+301

ISLisp>(carg (complex 1 0))
0.0
ISLisp>(carg (complex 1 1))
0.785398
ISLisp>(carg (complex 0 1))
1.5708
ISLisp>(carg (complex -1 1))
2.35619
ISLisp>(carg (complex -1 0))
3.14159
ISLisp>(carg (complex 1 -1))
-0.785398
ISLisp>(carg (complex 0 -1))
-1.5708
ISLisp>(carg (complex -1 -1))
-2.35619

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

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

除算の場合、絶対値の計算と同様にオーバーフローの対策が必要になります。今回は参考文献 1 のプログラムを ISLisp に移植しました。プログラムは次のようになります。 

リスト : 複素数の四則演算

;;; 加算
(defun cadd (x y)
  (complex (+ (realpart x) (realpart y)) (+ (imagpart x) (imagpart y))))

;;; 減算
(defun csub (x y)
  (complex (- (realpart x) (realpart y)) (- (imagpart x) (imagpart y))))

;;; 乗算
(defun cmul (x y)
  (complex (- (* (realpart x) (realpart y)) (* (imagpart x) (imagpart y)))
           (+ (* (realpart x) (imagpart y)) (* (imagpart x) (realpart y)))))

;;; 除算
(defun cdiv (x y)
  (if (>= (abs (realpart y)) (abs (imagpart y)))
      (let* ((u (quotient (imagpart y) (realpart y)))
             (v (+ (realpart y) (* (imagpart y) u))))
        (complex (quotient (+ (realpart x) (* (imagpart x) u)) v)
                 (quotient (- (imagpart x) (* (realpart x) u)) v)))
    (let* ((u (quotient (realpart y) (imagpart y)))
           (v (+ (* (realpart y) u) (imagpart y))))
      (complex (quotient (+ (* (realpart x) u) (imagpart x)) v)
               (quotient (- (* (imagpart x) u) (realpart x)) v)))))

それでは実際に試してみましょう。 

ISLisp>(defglobal a (complex 1 2))
A
ISLisp>(defglobal b (complex 3 4))
B
ISLisp>(cprint a)
#C(1.0 2.0)
NIL
ISLisp>(cprint b)
#C(3.0 4.0)
NIL
ISLisp>(cprint (cadd a b))
#C(4.0 6.0)
NIL
ISLisp>(cprint (csub a b))
#C(-2.0 -2.0)
NIL
ISLisp>(cprint (cmul a b))
#C(-5.0 10.0)
NIL
ISLisp>(cprint (cdiv a b))
#C(0.44 0.08)
NIL
ISLisp>(cprint (cdiv (complex 1 0) (complex 1e300 1e300)))
#C(5e-301 -5e-301)
NIL
ISLisp>(cprint (cdiv (complex 1 0) (complex 1e301 1e300)))
#C(9.90099e-302 -9.90099e-303)
NIL
ISLisp>(cprint (cdiv (complex 1 0) (complex 1e300 1e301)))
#C(9.90099e-303 -9.90099e-302)
NIL

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

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

\( \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 (x + iy) &=& \log |z| e^{i\theta} \\ &=& \log |z| + \log_e e^{i\theta} \\ &=& \log |z| + i\theta, \quad (-\pi \leq \theta \leq \pi) \end{eqnarray} \end{array} \)

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

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

プログラムと簡単な実行例を示します。 

リスト : 複素関数

;;; 指数関数
(defun cexp (z)
  (let ((a (* (exp (realpart z)))))
    (complex (* (cos (imagpart z)) a) (* (sin (imagpart z)) a))))

;;; 対数関数
(defun clog (z) (complex (log (cabs z)) (carg z)))

;;; べき乗
(defun cpow (x y) (cexp (cmul y (clog x))))

ISLisp>*pi*
3.14159
ISLisp>(cprint (cexp (complex 0 0)))
#C(1.0 0.0)
NIL
ISLisp>(cprint (cexp (complex 0 (quotient *pi* 4))))
#C(0.707107 0.707107)
NIL
ISLisp>(cprint (cexp (complex 0 (quotient *pi* 2))))
#C(6.12303e-017 1.0)
NIL
ISLisp>(cprint (cexp (complex 0 *pi*)))
#C(-1.0 1.22461e-016)
NIL
ISLisp>(cprint (cexp (complex 1 1)))
#C(1.46869 2.28736)
NIL

ISLisp>(cprint (clog (complex 1 1)))
#C(0.346574 0.785398)
NIL
ISLisp>(cprint (clog (complex 1 0)))
#C(0.0 0.0)
NIL
ISLisp>(cprint (clog (complex 0 1)))
#C(0.0 1.5708)
NIL
ISLisp>(cprint (clog (complex 1 -1)))
#C(0.346574 -0.785398)
NIL
ISLisp>(cprint (clog (complex 1e300 1e300)))
#C(691.122 0.785398)
NIL

ISLisp>(cprint (cpow (complex 1 1) (complex 0 0)))
#C(1.0 0.0)
NIL
ISLisp>(cprint (cpow (complex 1 1) (complex 1 0)))
#C(1.0 1.0)
NIL
ISLisp>(cprint (cpow (complex 1 1) (complex 2 0)))
#C(1.22461e-016 2.0)
NIL
ISLisp>(cprint (cpow (complex 1 1) (complex 3 0)))
#C(-2.0 2.0)
NIL
ISLisp>(cprint (cpow (complex 1 2) (complex 3 4)))
#C(0.12901 0.0339241)
NIL
ISLisp>(cprint (cpow (complex 1 1) (complex 1 1)))
#C(0.273957 0.583701)
NIL

●複素数の三角関数

複素数の三角関数の定義は、オイラーの公式から導かれる式の \(\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}\)

プログラムは次のようになります。 

リスト : 複素数の三角関数

(defun csin (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))

(defun ccos (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))

(defun ctan (z)
  (let* ((x (realpart z))
         (y (imagpart z))
         (d (+ (cos (* 2 x)) (cosh (* 2 y)))))
    (complex (quotient (sin (* 2 x)) d) (quotient (sinh (* 2 y)) d))))

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

ISLisp>(cprint (csin (complex 0 1)))
#C(0.0 1.1752)
NIL
ISLisp>(cabs (csin (complex 0 1)))
1.1752
ISLisp>(cprint (csin (complex 1 1)))
#C(1.29846 0.634964)
NIL
ISLisp>(cabs (csin (complex 1 1)))
1.4454
ISLisp>(cprint (csin (complex 0 0)))
#C(0.0 0.0)
NIL

ISLisp>(cprint (ccos (complex 0 1)))
#C(1.54308 0.0)
NIL
ISLisp>(cabs (ccos (complex 0 1)))
1.54308
ISLisp>(cprint (ccos (complex 1 1)))
#C(0.83373 -0.988898)
NIL
ISLisp>(cabs (ccos (complex 1 1)))
1.29345
ISLisp>(cprint (ccos (complex 0 0)))
#C(1.0 0.0)
NIL

ISLisp>(cprint (ctan (complex 0 1)))
#C(0.0 0.761594)
NIL
ISLisp>(cabs (ctan (complex 0 1)))
0.761594
ISLisp>(cprint (ctan (complex 1 1)))
#C(0.271753 1.08392)
NIL
ISLisp>(cabs (ctan (complex 1 1)))
1.11747
ISLisp>(cprint (ctan (complex 0 0)))
#C(0.0 0.0)
NIL
-- 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}\)

プログラムと簡単な実行例を示します。 

リスト : 複素数の双曲線関数

(defun csinh (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (sinh x) (cos y)) (* (cosh x) (sin y)))))

(defun ccosh (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (cosh x) (cos y)) (* (sinh x) (sin y)))))

(defun ctanh (z)
  (let* ((x (realpart z))
          (y (imagpart z))
         (d (+ (cosh (* 2 x)) (cos (* 2 y)))))
    (complex (quotient (sinh (* 2 x)) d) (quotient (sin (* 2 y)) d))))

ISLisp>(cprint (csinh (complex 0 1)))
#C(0.0 0.841471)
NIL
ISLisp>(cprint (csinh (complex 1 1)))
#C(0.634964 1.29846)
NIL
ISLisp>(cprint (csinh (complex 1 0)))
#C(1.1752 0.0)
NIL
ISLisp>(cprint (csinh (complex 0 0)))
#C(0.0 0.0)
NIL

ISLisp>(cprint (ccosh (complex 0 1)))
#C(0.540302 0.0)
NIL
ISLisp>(cprint (ccosh (complex 1 1)))
#C(0.83373 0.988898)
NIL
ISLisp>(cprint (ccosh (complex 1 0)))
#C(1.54308 0.0)
NIL
ISLisp>(cprint (ccosh (complex 0 0)))
#C(1.0 0.0)
NIL

ISLisp>(cprint (ctanh (complex 0 1)))
#C(0.0 1.55741)
NIL
ISLisp>(cprint (ctanh (complex 1 1)))
#C(1.08392 0.271753)
NIL
ISLisp>(cprint (ctanh (complex 1 0)))
#C(0.761594 0.0)
NIL
ISLisp>(cprint (ctanh (complex 0 0)))
#C(0.0 0.0)
NIL

●複素数の平方根

複素数 z の平方根は次の式で求めることができます。

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

プログラムと実行例を示します。 

リスト : 複素数の平方根

(defun csqrt (z)
  (let* ((x (realpart z))
         (a (cabs z))
         (b (sqrt (quotient (- a x) 2))))
    (complex (sqrt (quotient (+ a x) 2))
             (if (>= (imagpart z) 0) b (- b)))))

ISLisp>(cprint (csqrt (complex 1 0)))
#C(1.0 0.0)
NIL
ISLisp>(cprint (csqrt (complex 2 0)))
#C(1.41421 0.0)
NIL
ISLisp>(cprint (csqrt (complex 3 0)))
#C(1.73205 0.0)
NIL

ISLisp>(cprint (csqrt (complex -1 0)))
#C(0.0 1.0)
NIL
ISLisp>(cprint (csqrt (complex -2 0)))
#C(0.0 1.41421)
NIL
ISLisp>(cprint (csqrt (complex -3 0)))
#C(0.0 1.73205)
NIL

ISLisp>(defglobal a (csqrt (complex 1 1)))
A
ISLisp>(cprint a)
#C(1.09868 0.45509)
NIL
ISLisp>(cprint (cmul a a))
#C(1.0 1.0)
NIL
ISLisp>(defglobal a (csqrt (complex 1 -1)))
A
ISLisp>(cprint a)
#C(1.09868 -0.45509)
NIL
ISLisp>(cprint (cmul a a))
#C(1.0 -1.0)
NIL
ISLisp>(defglobal a (csqrt (complex -1 1)))
A
ISLisp>(cprint a)
#C(0.45509 1.09868)
NIL
ISLisp>(cprint (cmul a a))
#C(-1.0 1.0)
NIL
ISLisp>(defglobal a (csqrt (complex -1 -1)))
A
ISLisp>(cprint a)
#C(0.45509 -1.09868)
NIL
ISLisp>(cprint (cmul a a))
#C(-1.0 -1.0)
NIL

●逆三角関数

三角関数の逆関数を「逆三角関数 (inverse trigonometric function)」といいます。Lisp / Scheme で用いられる逆三角関数を以下に示します。

ここでは引数 x を実数とします。asin x は引数 x が与えられたとき sin w = x となる角度 w を求めます。同様に acos x は cos w = x となる角度 w を、atan x は tan x = w となる角度 w を求めます。Web 上で数式を表示する JavaScript ライブラリ 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}\)

ISLisp の仕様で定義されている逆三角関数には atan と atan2 があります。atan2 は他のプログラミング言語、たとえばＣ言語の数学関数 atan2 と同じです。 

(atan2 y x) => 角度 (ラジアン)

引数 x, y は実数です。atan2 は直交座標系においてベクトル (x, y) と x 軸との角度を求める関数です。複素平面で考えると、複素数 \(x + iy\) の偏角 \(\theta\) を求めることと同じです。OKI-ISLisp の場合、返り値 (角度 \(\theta\)) の範囲は \(-\pi \lt \theta \leq \pi\) になります。

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

ISLisp>(atan2 1.0 0.0)
1.5708
ISLisp>(atan2 0.0 1.0)
0.0
ISLisp>(atan2 1.0 1.0)
0.785398
ISLisp>(atan2 1.0 0.0)
1.5708
ISLisp>(atan2 1.0 -1.0)
2.35619
ISLisp>(atan2 0.0 -1.0)
3.14159
ISLisp>(atan2 -1.0 1.0)
-0.785398
ISLisp>(atan2 -1.0 0.0)
-1.5708
ISLisp>(atan2 -1.0 -1.0)
-2.35619
ISLisp>(atan2 -0.1 -1.0)
-2.35619
ISLisp>(atan2 -0.01 -1.0)
-3.04192
ISLisp>(atan2 -0.001 -1.0)
-3.13159
ISLisp>(atan2 -0.0001 -1.0)
-3.14059

\(\arcsin\) と \(\arccos\) は atan2 から求めることができます。

原点が中心で半径が 1 の円を考える
円周上の座標を (x, y) とし (x, y) - (0, 0) - (x, 0) の角度を \(\theta\) とする
三角関数の定義 \(\sin \theta = y, \ \cos \theta = x \ \) により

\( \begin{eqnarray} \arcsin y &=& \theta \\ &=& tan2(y, x) \\ &=& tan2\left(y, \sqrt{1 - y^2}\right) \\ \arccos x &=& \theta \\ &=& tan2(y, x) \\ &=& tan2\left(\sqrt{1 - x^2}, x\right) \end{eqnarray} \)

プログラムは次のようになります。 

リスト : 逆三角関数

(defun asin (x) (atan2 x (sqrt (- 1 (* x x)))))

(defun acos (x) (atan2 (sqrt (- 1 (* x x))) x))

逆三角関数には次の関係が成り立ちます。

\(\begin{array}{l} \arcsin x = \arctan{\dfrac{x}{\sqrt{1 - x^2}}} \\ \arccos x = \dfrac{\pi}{2} - \arcsin x \end{array}\)

上式を使って \(\arcsin, \ \arccos\) を定義することもできます。

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

ISLisp>(defun display (x) (format (standard-output) "~S~%" x))
DISPLAY
ISLisp>(mapc (lambda (x) (display (asin x))) '(-1.0 -0.5 0.0 0.5 1.0))
-1.5708
-0.523599
0.0
0.523599
1.5708
(-1.0 -0.5 0.0 0.5 1.0)
ISLisp>(mapc (lambda (x) (display (acos x))) '(-1.0 -0.5 0.0 0.5 1.0))
3.14159
2.0944
1.5708
1.0472
0.0
(-1.0 -0.5 0.0 0.5 1.0)
ISLisp>(mapc (lambda (x) (display (atan x))) '(-1e300 -1.0 0.0 1.0 1e300))
-1.5708
-0.785398
0.0
0.785398
1.5708
(-1e+300 -1.0 0.0 1.0 1e+300)

●複素数の逆三角関数

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

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

リスト : 複素数の逆三角関数

;;; 定数
(defconstant iunit  (complex 0 1))
(defconstant iunit- (complex 0 -1))
(defconstant cone   (complex 1 0))
(defconstant ctwo   (complex 2 0))
(defconstant cpi2   (complex (quotient *pi* 2) 0))

;;; 複素数の逆三角関数
(defun casin (z)
  (let ((z1 (csqrt (cadd cone z)))
        (z2 (csqrt (csub cone z))))
    (cmul iunit- (clog (cadd (cmul iunit z) (cmul z1 z2))))))

(defun cacos (z) (csub cpi2 (casin z)))

(defun catan (z)
  (let ((z1 (cmul iunit z)))
    (cdiv (cmul iunit
                (csub (clog (csub cone z1)) (clog (cadd cone z1))))
          ctwo)))

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

ISLisp>(cprint (casin ctwo))
#C(1.5708 -1.31696)
NIL
ISLisp>(cprint (csin (casin ctwo)))
#C(2.0 -1.06054e-016)
NIL
ISLisp>(cprint (casin iunit))
#C(0.0 0.881374)
NIL
ISLisp>(cprint (csin (casin iunit)))
#C(0.0 1.0)
NIL
ISLisp>(cprint (casin (complex 1 1)))
#C(0.666239 1.06128)
NIL
ISLisp>(cprint (csin (casin (complex 1 1))))
#C(1.0 1.0)
NIL

ISLisp>(cprint (cacos ctwo))
#C(0.0 1.31696)
NIL
ISLisp>(cprint (ccos (cacos ctwo)))
#C(2.0 0.0)
NIL
ISLisp>(cprint (cacos iunit))
#C(1.5708 -0.881374)
NIL
ISLisp>(cprint (ccos (cacos iunit)))
#C(8.65927e-017 1.0)
NIL
ISLisp>(cprint (cacos (complex 1 1)))
#C(0.904557 -1.06128)
NIL
ISLisp>(cprint (ccos (cacos (complex 1 1))))
#C(1.0 1.0)
NIL

ISLisp>(cprint (catan (complex 1 1)))
#C(1.01722 0.402359)
NIL
ISLisp>(cprint (ctan (catan (complex 1 1))))
#C(1.0 1.0)
NIL

●逆双曲線関数

双曲線関数の逆関数を「逆双曲線関数 (inverse hyperbolic function)」といいます。Lisp / Scheme で用いられる逆三角関数には asinh, acosh, atanh があります。ここでは逆双曲線関数を \(\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} \\ \end{array}\)
\(\begin{array}{l} 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}\)

関数 atanh は ISLisp の仕様に規定されているので、次の関係式から asinh と acosh を求めることができます。

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

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

プログラムと簡単な実行例を示します。 

リスト : 逆双曲線関数

(defun asinh (x) (atanh (quotient x (sqrt (+ 1 (* x x))))))

(defun acosh (x)
  (let ((a (sqrt (+ x 1.0)))
        (b (sqrt (- x 1.0))))
    (atanh (quotient (* a b) x))))

ISLisp>(mapcar (lambda (x) (asinh x))
 '(-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0))
(-1.44364 -1.19476 -0.881374 -0.481212 0.0 0.481212 0.881374 1.19476 1.44364)
ISLisp>(mapcar #'sinh (mapcar (lambda (x) (asinh x))
 '(-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0)))
(-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0)

ISLisp>(mapcar (lambda (x) (acosh x)) '(1.0 1.5 2.0 2.5 3.0 3.5 4.0))
(0.0 0.962424 1.31696 1.5668 1.76275 1.92485 2.06344)
ISLisp>(mapcar #'cosh (mapcar (lambda (x) (acosh x))
 '(1.0 1.5 2.0 2.5 3.0 3.5 4.0)))
(1.0 1.5 2.0 2.5 3.0 3.5 4.0)

ISLisp>(mapcar (lambda (x) (atanh x)) '(-0.75 -0.5 -0.25 0.0 0.25 0.5 0.75))
(-0.972955 -0.549306 -0.255413 0.0 0.255413 0.549306 0.972955)
ISLisp>(mapcar #'tanh (mapcar (lambda (x) (atanh x))
 '(-0.75 -0.5 -0.25 0.0 0.25 0.5 0.75)))
(-0.75 -0.5 -0.25 0.0 0.25 0.5 0.75)

●複素数の逆双曲線関数

複素数の逆双曲線関数の定義は、実数の定義で引数 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) \)

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

リスト : 複素数の逆双曲線関数

(defun casinh (z)
  (clog (cadd z (csqrt (cadd cone (cmul z z))))))

(defun cacosh (z)
  (let ((z1 (csqrt (cadd z cone)))
	(z2 (csqrt (csub z cone))))
    (clog (cadd z (cmul z1 z2)))))

(defun catanh (z)
  (cdiv (csub (clog (cadd cone z)) (clog (csub cone z))) ctwo))

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

ISLisp>(mapc (lambda (x) (cprint (casinh x)))
 (mapcar (lambda (x) (complex x 0.0))
 '(-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0)))
#C(-1.44364 0.0)
#C(-1.19476 0.0)
#C(-0.881374 0.0)
#C(-0.481212 0.0)
#C(0.0 0.0)
#C(0.481212 0.0)
#C(0.881374 0.0)
#C(1.19476 0.0)
#C(1.44364 0.0)
(... 略 ...)
ISLisp>(mapc (lambda (x) (cprint (csinh (casinh x))))
 (mapcar (lambda (x) (complex x 0.0))
 '(-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0)))
#C(-2.0 0.0)
#C(-1.5 0.0)
#C(-1.0 0.0)
#C(-0.5 0.0)
#C(0.0 0.0)
#C(0.5 0.0)
#C(1.0 0.0)
#C(1.5 0.0)
#C(2.0 0.0)
(... 略 ...)
ISLisp>(cprint (casinh (complex 0 1)))
#C(0.0 1.5708)
NIL
ISLisp>(cprint (csinh (casinh (complex 0 1))))
#C(0.0 1.0)
NIL
ISLisp>(cprint (casinh (complex 0 -1)))
#C(0.0 -1.5708)
NIL
ISLisp>(cprint (csinh (casinh (complex 0 -1))))
#C(0.0 -1.0)
NIL
ISLisp>(cprint (casinh (complex 1 1)))
#C(1.06128 0.666239)
NIL
ISLisp>(cprint (csinh (casinh (complex 1 1))))
#C(1.0 1.0)
NIL

ISLisp>(mapc (lambda (x) (cprint (cacosh x))) (mapcar (lambda (x) (complex x 0))
 '(1.0 1.5 2.0 2.5 3.0 3.5 4.0)))
#C(0.0 0.0)
#C(0.962424 0.0)
#C(1.31696 0.0)
#C(1.5668 0.0)
#C(1.76275 0.0)
#C(1.92485 0.0)
#C(2.06344 0.0)
(... 略 ...)
ISLisp>(mapc (lambda (x) (cprint (ccosh (cacosh x))))
 (mapcar (lambda (x) (complex x 0))
 '(1.0 1.5 2.0 2.5 3.0 3.5 4.0)))
#C(1.0 0.0)
#C(1.5 0.0)
#C(2.0 0.0)
#C(2.5 0.0)
#C(3.0 0.0)
#C(3.5 0.0)
#C(4.0 0.0)
(... 略 ...)
ISLisp>(cprint (cacosh (complex 0 1)))
#C(0.881374 1.5708)
NIL
ISLisp>(cprint (ccosh (cacosh (complex 0 1))))
#C(8.65927e-017 1.0)
NIL
ISLisp>(cprint (cacosh (complex 0 -1)))
#C(0.881374 -1.5708)
NIL
ISLisp>(cprint (ccosh (cacosh (complex 0 -1))))
#C(8.65927e-017 -1.0)
NIL
ISLisp>(cprint (cacosh (complex 1 1)))
#C(1.06128 0.904557)
NIL
ISLisp>(cprint (ccosh (cacosh (complex 1 1))))
#C(1.0 1.0)
NIL

ISLisp>(mapc (lambda (x) (cprint (catanh x))) (mapcar (lambda (x) (complex x 0))
 '(-0.75 -0.5 -0.25 0.0 0.25 0.5 0.75)))
#C(-0.972955 0.0)
#C(-0.549306 0.0)
#C(-0.255413 0.0)
#C(0.0 0.0)
#C(0.255413 0.0)
#C(0.549306 0.0)
#C(0.972955 0.0)
(... 略 ...)
ISLisp>(mapc (lambda (x) (cprint (ctanh (catanh x))))
 (mapcar (lambda (x) (complex x 0)) '(-0.75 -0.5 -0.25 0.0 0.25 0.5 0.75)))
#C(-0.75 0.0)
#C(-0.5 0.0)
#C(-0.25 0.0)
#C(0.0 0.0)
#C(0.25 0.0)
#C(0.5 0.0)
#C(0.75 0.0)
(... 略 ...)
ISLisp>(cprint (catanh (complex 0 1)))
#C(0.0 0.785398)
NIL
ISLisp>(cprint (ctanh (catanh (complex 0 1))))
#C(0.0 1.0)
NIL
ISLisp>(cprint (catanh (complex 0 -1)))
#C(0.0 -0.785398)
NIL
ISLisp>(cprint (ctanh (catanh (complex 0 -1))))
#C(0.0 -1.0)
NIL
ISLisp>(cprint (catanh (complex 1 1)))
#C(0.402359 1.01722)
NIL
ISLisp>(cprint (ctanh (catanh (complex 1 1))))
#C(1.0 1.0)
NIL

●参考文献, URL

  1. 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. 逆双曲線関数と逆三角関数の branch cut | 雑記帳

●プログラムリスト

;;;
;;; complex.l : ISLisp 用複素数
;;;
;;;             Copyright (C) 2021 Makoto Hiroi
;;;

;;; 複素数
(defclass <complex> ()
  ((re :accessor realpart :initform 0.0 :initarg re)
   (im :accessor imagpart :initform 0.0 :initarg im)))

;;; コンストラクタ
(defun complex (x y)
  (create (class <complex>) 're (convert x <float>) 'im (convert y <float>)))

;;; 表示
(defun cprint (z)
  (format (standard-output) "#C(~G ~G)~%" (realpart z) (imagpart z)))

;;; 複素共役
(defun conj (z)
  (complex (realpart z) (- (imagpart z))))

;;; 絶対値
(defun cabs (z)
  (cond
   ((= (realpart z) 0.0) (abs (imagpart z)))
   ((= (imagpart z) 0.0) (abs (realpart z)))
   ((> (abs (imagpart z)) (abs (realpart z)))
    (let ((temp (quotient (realpart z) (imagpart z))))
      (* (abs (imagpart z)) (sqrt (+ 1 (* temp temp))))))
   (t
    (let ((temp (quotient (imagpart z) (realpart z))))
      (* (abs (realpart z)) (sqrt (+ 1 (* temp temp))))))))

;;; 偏角
(defun carg (z) (atan2 (imagpart z) (realpart z)))

;;; 四則演算
(defun cadd (x y)
  (complex (+ (realpart x) (realpart y)) (+ (imagpart x) (imagpart y))))

(defun csub (x y)
  (complex (- (realpart x) (realpart y)) (- (imagpart x) (imagpart y))))

(defun cmul (x y)
  ;; (a + bi)(c + di) = ac + adi + bci - bd = (ac - bd) + (ad + bc)i
  (complex (- (* (realpart x) (realpart y)) (* (imagpart x) (imagpart y)))
           (+ (* (realpart x) (imagpart y)) (* (imagpart x) (realpart y)))))

(defun cdiv (x y)
  (if (>= (abs (realpart y)) (abs (imagpart y)))
      (let* ((u (quotient (imagpart y) (realpart y)))
             (v (+ (realpart y) (* (imagpart y) u))))
        (complex (quotient (+ (realpart x) (* (imagpart x) u)) v)
                 (quotient (- (imagpart x) (* (realpart x) u)) v)))
    (let* ((u (quotient (realpart y) (imagpart y)))
           (v (+ (* (realpart y) u) (imagpart y))))
      (complex (quotient (+ (* (realpart x) u) (imagpart x)) v)
               (quotient (- (* (imagpart x) u) (realpart x)) v)))))

;;; 指数関数
(defun cexp (z)
  (let ((a (* (exp (realpart z)))))
    (complex (* (cos (imagpart z)) a) (* (sin (imagpart z)) a))))

;;; 対数関数
(defun clog (z) (complex (log (cabs z)) (carg z)))

;;; べき乗
(defun cpow (x y) (cexp (cmul y (clog x))))

;;; 三角関数
(defun csin (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))

(defun ccos (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))

(defun ctan (z)
  (let* ((x (realpart z))
         (y (imagpart z))
         (d (+ (cos (* 2 x)) (cosh (* 2 y)))))
    (complex (quotient (sin (* 2 x)) d) (quotient (sinh (* 2 y)) d))))

;;;双曲線関数
(defun csinh (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (sinh x) (cos y)) (* (cosh x) (sin y)))))

(defun ccosh (z)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (* (cosh x) (cos y)) (* (sinh x) (sin y)))))

(defun ctanh (z)
  (let* ((x (realpart z))
          (y (imagpart z))
         (d (+ (cosh (* 2 x)) (cos (* 2 y)))))
    (complex (quotient (sinh (* 2 x)) d) (quotient (sin (* 2 y)) d))))

;;; 平方根
(defun csqrt (z)
  (let* ((x (realpart z))
         (a (cabs z))
         (b (sqrt (quotient (- a x) 2))))
    (complex (sqrt (quotient (+ a x) 2))
             (if (>= (imagpart z) 0) b (- b)))))

;;; 逆三角関数 (atan と atan2 は ISLisp で定義されている)
(defun asin (x) (atan2 x (sqrt (- 1 (* x x)))))
(defun acos (x) (atan2 (sqrt (- 1 (* x x))) x))

;;; 定数
(defconstant iunit  (complex 0 1))
(defconstant iunit- (complex 0 -1))
(defconstant cone   (complex 1 0))
(defconstant ctwo   (complex 2 0))
(defconstant cpi2   (complex (quotient *pi* 2) 0))

;;; 複素数の逆三角関数 (定義式をそのままプログラムする)
(defun casin (z)
  (let ((z1 (csqrt (cadd cone z)))
        (z2 (csqrt (csub cone z))))
    (cmul iunit- (clog (cadd (cmul iunit z) (cmul z1 z2))))))

(defun cacos (z) (csub cpi2 (casin z)))

(defun catan (z)
  (let ((z1 (cmul iunit z)))
    (cdiv (cmul iunit
                (csub (clog (csub cone z1)) (clog (cadd cone z1))))
          ctwo)))

;;; 逆双曲線関数 (atanh は ISLisp に定義されている)
(defun asinh (x) (atanh (quotient x (sqrt (+ 1 (* x x))))))

(defun acosh (x)
  (let ((a (sqrt (+ x 1.0)))
        (b (sqrt (- x 1.0))))
    (atanh (quotient (* a b) x))))

;;; 複素数の逆双曲線関数 (定義式をそのままプログラムする)
(defun casinh (z)
  (clog (cadd z (csqrt (cadd cone (cmul z z))))))

(defun cacosh (z)
  (let ((z1 (csqrt (cadd z cone)))
	(z2 (csqrt (csub z cone))))
    (clog (cadd z (cmul z1 z2)))))

(defun catanh (z)
  (cdiv (csub (clog (cadd cone z)) (clog (csub cone z))) ctwo))