特殊関数: 楕円積分・楕円関数・テータ関数

概要

楕円積分・楕円関数・テータ関数は、互いに逆関数や恒等式で密に結びついた一群をなす。 ざっくりした関係は次のとおり:

• 楕円積分 $F(\varphi, k)$ は、振幅 $\varphi$ を引数に取り「弧長」を返す積分。
• Jacobi 楕円関数 $\operatorname{sn}, \operatorname{cn}, \operatorname{dn}$ は、その逆関数。 $u = F(\varphi, k)$ に対し $\operatorname{sn}(u, k) = \sin\varphi$ 等が成り立つ。
• テータ関数 $\vartheta_1 \ldots \vartheta_4$ は、急速収束する $q$ 級数で楕円関数を組み立てる「素材」。 Jacobi 関数も Weierstrass 関数もテータ関数の比で書ける。
• Weierstrass の $\wp$ 関数は、周期格子 $2\omega_1 \mathbb{Z} + 2\omega_2 \mathbb{Z}$ 上の二重周期関数で、 格子点に 2 位の極を持つ。微分方程式 $(\wp')^2 = 4\wp^3 - g_2\wp - g_3$ を満たす。

評価アルゴリズムは大きく 2 系統に分かれる。楕円積分は Carlson 対称形 + 倍化定理、 または AGM (算術幾何平均) の反復で計算する。テータ関数・Weierstrass 関数は $q$ 級数 (必要に応じてモジュラー還元) で計算する。いずれも反復・級数が二次収束ないし超指数収束し、少数のステップで高精度に達する。

関連 API: ellipticK, carlsonRF, jacobiSn ほか (jacobiTheta1 …)。

Carlson 対称形

Carlson は、引数について対称な 4 つの基本積分で楕円積分を統一的に表した:

$$R_F(x,y,z) = \frac{1}{2}\int_0^\infty \frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}$$

$$R_D(x,y,z) = \frac{3}{2}\int_0^\infty \frac{dt}{(t+z)^{3/2}\sqrt{(t+x)(t+y)}}$$

$$R_J(x,y,z,p) = \frac{3}{2}\int_0^\infty \frac{dt}{(t+p)\sqrt{(t+x)(t+y)(t+z)}}, \qquad R_C(x,y) = \frac{1}{2}\int_0^\infty \frac{dt}{(t+y)\sqrt{t+x}}$$

$R_C(x,y) = R_F(x,y,y)$ は退化形で、$R_J$ の内部で繰り返し使われる。

倍化定理による反復評価

Carlson 形の威力は倍化定理 (duplication theorem) にある。 各反復で $\lambda = \sqrt{x}\sqrt{y} + \sqrt{y}\sqrt{z} + \sqrt{z}\sqrt{x}$ を作り、引数を

$$x \leftarrow \frac{x+\lambda}{4}, \quad y \leftarrow \frac{y+\lambda}{4}, \quad z \leftarrow \frac{z+\lambda}{4}$$

と更新する。これを繰り返すと $x, y, z$ は互いに近づき、共通の平均 $A = (x+y+z)/3$ のまわりに集まる。 相対偏差 $\delta_x = (A-x)/A$ 等がすべて小さくなったら反復を止め、$A$ のまわりの補正多項式を 1 回足す:

$$R_F \approx \frac{1}{\sqrt{A}}\left(1 - \frac{E_2}{10} + \frac{E_3}{14} + \frac{E_2^2}{24} - \frac{3\,E_2 E_3}{44}\right)$$

ここで $E_2, E_3$ は偏差 $\delta_x, \delta_y, \delta_z$ の対称式。 偏差が小さいので多項式の項数はわずかで済み、全体として反復回数に対して指数的に精度が上がる。

$R_D, R_J$ では各反復で総和項 $\sigma$ も蓄積する点が異なる ($R_J$ の総和項は $R_C$ を 1 回呼ぶ)。sangi の carlsonRF / carlsonRD / carlsonRJ / carlsonRC はいずれもこの反復 (収束判定は機械イプシロン、上限 100 反復) を実装する。

完全・不完全楕円積分を Carlson 形で表す

Legendre 標準形は Carlson 形の組み合わせで書ける。これにより、端点近傍での桁落ちを避けつつ統一的に評価できる:

$$K(k) = R_F(0,\, 1-k^2,\, 1), \qquad E(k) = R_F(0,\, 1-k^2,\, 1) - \frac{k^2}{3} R_D(0,\, 1-k^2,\, 1)$$

$$\Pi(n,k) = R_F(0,\, 1-k^2,\, 1) + \frac{n}{3} R_J(0,\, 1-k^2,\, 1,\, 1-n)$$

不完全形は $\cos^2\varphi$ などを引数に入れる:

$$F(\varphi, k) = \sin\varphi \cdot R_F(\cos^2\varphi,\ 1 - k^2\sin^2\varphi,\ 1)$$

$$E(\varphi, k) = \sin\varphi \cdot R_F(\cdots) - \frac{k^2}{3}\sin^3\varphi \cdot R_D(\cdots)$$

sangi はこれらを ellipticK, ellipticE, ellipticPi (引数 1 個で完全、2〜3 個で不完全) として提供する。$k^2 \geq 1$ では $K(\pm 1) = \infty$ を返すなど、端点の特異性も扱う。 Carlson 形は複素引数にも自然に拡張でき、Complex<double> / Complex<Float> でも同じ反復が動く (収束判定を $|\delta| < \varepsilon$ のスカラー比較で行う)。

AGM (算術幾何平均)

$K(k)$ にはもう一つの古典的な高速評価がある。 算術幾何平均 $M(a, b)$ は、初期値 $a_0 = a,\ b_0 = b$ から

$$a_{n+1} = \frac{a_n + b_n}{2}, \qquad b_{n+1} = \sqrt{a_n b_n}$$

を反復した共通の極限である。算術平均と幾何平均が各ステップで近づくため、 $a_n - b_n$ は二次収束 (毎回ほぼ平方) でゼロに向かう。第一種完全楕円積分は AGM で次のように書ける:

$$K(k) = \frac{\pi}{2\, M(1,\ k')}, \qquad k' = \sqrt{1 - k^2}$$

二次収束ゆえ、目標精度 $B$ 桁に対し反復回数は $\mathcal{O}(\log B)$ にとどまる。 たとえば 1000 桁でも 10 数反復で足り、各反復は平方根 1 回なので、高精度計算では級数和より圧倒的に有利。

降下 Landen 変換との関係

AGM の各ステップは、母数 $k$ を小さくする降下 Landen 変換と等価である。 Landen 変換は

$$k_1 = \frac{1 - k'}{1 + k'} \quad (k' = \sqrt{1-k^2})$$

で新しい母数 $k_1 < k$ を作り、$K(k)$ を $K(k_1)$ に帰着させる。母数が反復ごとに $0$ に近づくと $K \to \pi/2$ なので、変換を遡って積分値を回復できる。 AGM 列 $(a_n, b_n)$ と Landen 列 $(k_n)$ は $b_n/a_n = k_n'$ で結びつき、本質的に同じ反復の二つの見方になっている。

sangi のテータ関数の高精度経路 (jacobiThetaZeroAGM) は、この AGM 構造を $z=0$ のテータ値に適用する。 $a_n = \vartheta_3(0, q^{2^n})^2,\ b_n = \vartheta_4(0, q^{2^n})^2$ が前進 AGM を満たす ($q^{2^n} \to 0$ なので $a_n, b_n \to 1$) ことを使い、十分小さい $q_N = q^{2^N}$ から級数で $a_N, b_N$ を出し、 逆 AGM $d = \sqrt{a_n^2 - b_n^2},\ a_{n-1} = a_n + d,\ b_{n-1} = a_n - d$ を $N$ 回遡って $\vartheta_3(0,q) = \sqrt{a_0}$, $\vartheta_4(0,q) = \sqrt{b_0}$, $\vartheta_2(0,q) = (a_0^2 - b_0^2)^{1/4}$ を得る。

Jacobi 楕円関数

Jacobi 楕円関数は楕円積分の逆関数として定義される。$u = F(\varphi, k)$ のとき、振幅 $\varphi = \operatorname{am}(u, k)$ とし

$$\operatorname{sn}(u, k) = \sin\varphi, \qquad \operatorname{cn}(u, k) = \cos\varphi, \qquad \operatorname{dn}(u, k) = \sqrt{1 - k^2 \operatorname{sn}^2(u, k)}$$

AGM / 降下 Landen による振幅の評価

sangi は振幅 $\operatorname{am}(u, k)$ を AGM 列と逆変換で計算する。手順は次のとおり:

• 前進 AGM: $a_0 = 1,\ b_0 = k',\ c_0 = k$ から $a_{n+1} = (a_n + b_n)/2,\ b_{n+1} = \sqrt{a_n b_n},\ c_{n+1} = (a_n - b_n)/2$ を反復し、 $c_n$ が十分小さくなる $N$ まで列 $\{a_n\}, \{c_n\}$ を保持する。
• 初期位相: $\varphi_N = 2^N a_N u$ をとる。
• 逆変換 (降下): $\varphi_{n-1} = \tfrac{1}{2}\bigl(\varphi_n + \arcsin(\tfrac{c_n}{a_n}\sin\varphi_n)\bigr)$ を $n = N$ から $1$ まで遡る。終端の $\varphi_0$ が $\operatorname{am}(u, k)$ になる。

あとは $\operatorname{sn} = \sin(\operatorname{am})$, $\operatorname{cn} = \cos(\operatorname{am})$, $\operatorname{dn} = \sqrt{1 - k^2 \operatorname{sn}^2}$ で得る。 境界の母数は閉じた形で扱う: $k = 0$ では $\operatorname{am}(u, 0) = u$、 $k = 1$ では $\operatorname{am}(u, 1) = \operatorname{gd}(u) = 2\arctan(\tanh(u/2))$ (Gudermann 関数)。

逆変換で $\frac{c_n}{a_n}\sin\varphi_n$ が数値誤差で $|\cdot| > 1$ になり得るため、 sangi の実装は $\arcsin$ の引数を $[-1, 1]$ にクランプして安全に評価する。 提供関数は jacobiSn, jacobiCn, jacobiDn (いずれも母数 $k$ を引数に取る)。

テータ関数

Jacobi のテータ関数は、ノーム $q$ (ただし $|q| < 1$) のべき級数で定義される。 指数が $n$ の二乗で増えるため、級数は超指数的に速く収束する:

$$\vartheta_1(z, q) = 2\sum_{n=0}^\infty (-1)^n q^{(n+1/2)^2} \sin((2n+1)z)$$

$$\vartheta_2(z, q) = 2\sum_{n=0}^\infty q^{(n+1/2)^2} \cos((2n+1)z)$$

$$\vartheta_3(z, q) = 1 + 2\sum_{n=1}^\infty q^{n^2} \cos(2nz), \qquad \vartheta_4(z, q) = 1 + 2\sum_{n=1}^\infty (-1)^n q^{n^2} \cos(2nz)$$

ノーム $q$ と半周期比 $\tau$

ノームは半周期比 $\tau$ から $q = e^{i\pi\tau}$ で定まる。$\tau$ が上半平面 ($\operatorname{Im}\tau > 0$) にあれば $|q| = e^{-\pi \operatorname{Im}\tau} < 1$ が保証される。sangi の qFromTau はこの変換を提供する。

級数の項は $|q|^{n^2}$ で減衰するので、$|q| \leq 0.5$ 程度なら数項で機械精度に達する。 実装は $q^{(n+1/2)^2}$ や $q^{n^2}$ を直接べき乗で求め、項が和に対して相対的に十分小さくなったら打ち切る。 複素ノーム $q$ では $q^k = \exp(k \log q)$ (半整数べき) と整数べき乗を使い分け、整数べきは安定な二進反復で計算する。

$\tau$ の還元 (モジュラー変換)

$|q|$ が $1$ に近い ($= \operatorname{Im}\tau$ が小さい) と級数の収束が遅くなる。 この場合は $\mathrm{SL}(2,\mathbb{Z})$ の生成元

• $T:\ \tau \to \tau + 1$ — $\vartheta_1, \vartheta_2$ に位相 $e^{i\pi/4}$ がつき、$\vartheta_3 \leftrightarrow \vartheta_4$ が入れ替わる
• $S:\ \tau \to -1/\tau$ — $(-i\tau)^{-1/2}$ と $\exp(i\pi z^2/\tau)$ の前因子がつき、添字が置換される

を繰り返し、$\tau$ を基本領域 ($|\tau| \geq 1$ かつ $|\operatorname{Re}\tau| \leq 1/2$、 すなわち $\operatorname{Im}\tau \geq \sqrt{3}/2 \approx 0.866$) に還元してから級数を評価する。 これにより常に十分小さい $|q|$ で計算でき、収束が保証される。 sangi の $\tau$ 形 API (jacobiTheta1Tau … jacobiTheta4Tau) はこの還元を内部で行う。

高精度の $z = 0$ 値については、前節の AGM 加速 (jacobiThetaZeroAGM) を使うと、 級数の項数を増やすより少ない反復で $\vartheta_2, \vartheta_3, \vartheta_4$ を同時に得られる。 提供関数は jacobiTheta1 … jacobiTheta4 (実引数・複素引数・複素ノーム・任意精度の各オーバーロード)。

Weierstrass の楕円関数

Weierstrass の $\wp$ 関数は、半周期 $\omega_1, \omega_2$ が張る格子上の二重周期関数で、 各格子点に 2 位の極を持つ偶関数である。微分方程式

$$(\wp')^2 = 4\wp^3 - g_2 \wp - g_3$$

を満たし、$g_2, g_3$ を不変量と呼ぶ。$\zeta_W$ と $\sigma_W$ は $\zeta_W'(z) = -\wp(z)$, $\sigma_W'(z)/\sigma_W(z) = \zeta_W(z)$ で結ばれる (これらは Riemann のゼータ・シグマとは別物)。

テータ関数経由の評価

格子和を直接取る代わりに、sangi はテータ関数による表現 (DLMF 23.6) で評価する。 まず $\tau = \omega_2/\omega_1$ が上半平面にあることを確認し、ノーム $q = e^{i\pi\tau}$ と $\alpha = \pi/(2\omega_1)$、 $v = \alpha z$ を用意する。$\wp$ は

$$\wp(z) = e_3 + \left(\alpha\, \vartheta_2(0,q)\, \vartheta_3(0,q)\, \frac{\vartheta_4(v,q)}{\vartheta_1(v,q)}\right)^2, \qquad e_3 = -\alpha^2\, \frac{\vartheta_2^4(0,q) + \vartheta_3^4(0,q)}{3}$$

で与えられる。格子点では $\vartheta_1(v) \approx 0$ となるので、そこを極として無限大を返す。 $\zeta_W$ と $\sigma_W$ も $\vartheta_1$ とその微分 $\vartheta_1', \vartheta_1'''$ の比で書ける:

$$\zeta_W(z) = \eta_1 \frac{z}{\omega_1} + \alpha\, \frac{\vartheta_1'(v,q)}{\vartheta_1(v,q)}, \qquad \eta_1 = -\frac{\pi^2}{12\,\omega_1}\, \frac{\vartheta_1'''(0,q)}{\vartheta_1'(0,q)}$$

$$\sigma_W(z) = \frac{2\omega_1}{\pi}\, \frac{\vartheta_1(v,q)}{\vartheta_1'(0,q)}\, \exp\!\left(\frac{\eta_1 z^2}{2\omega_1}\right)$$

テータ級数が速く収束するため、格子和を直接打ち切るより安定かつ高速に評価できる。 sangi の weierstrassP, weierstrassZeta, weierstrassSigma は半周期 $\omega_1, \omega_2$ を引数に取り、 Complex<double> で任意の半周期・任意の $z$ に対応する。 実引数版は $\omega_1$ 実・$\omega_2$ 純虚 (長方形格子) を扱い、極を避ければ実数値を返す。 一般のノーム $q$ (= $\operatorname{Re}\tau \neq 0$) に対応するため、内部では複素 $q$ のテータ級数を直接展開する。

比較表

大まかな指針: 楕円積分は引数対称で頑健な Carlson 形が既定。$K(k)$ を多数回・高桁で要るなら AGM。 Jacobi 関数は AGM + 降下 Landen の逆変換。テータ関数・Weierstrass 関数は $q$ 級数 (収束が遅ければ $\tau$ をモジュラー還元)。

参考文献

• Carlson, B. C. (1995). "Numerical computation of real or complex elliptic integrals". Numerical Algorithms, 10(1), 13–26.
• Olver, F. W. J. et al. (eds.). NIST Digital Library of Mathematical Functions (DLMF), §19 (楕円積分), §20 (テータ関数), §22 (Jacobi 楕円関数), §23 (Weierstrass 楕円関数). https://dlmf.nist.gov/
• Abramowitz, M. & Stegun, I. A. (1964). Handbook of Mathematical Functions. National Bureau of Standards. (§16 Jacobi 楕円関数, §17 楕円積分)
• Whittaker, E. T. & Watson, G. N. (1927). A Course of Modern Analysis. 4th ed. Cambridge University Press. (Ch. 20–22)
• Borwein, J. M. & Borwein, P. B. (1987). Pi and the AGM. Wiley.