float_backend 設計

このページでは float_backend が Complex[Double] 向けに実装する式を導出し、それぞれがオーバーフロー、アンダーフロー、桁落ちをどう回避するかを説明し、分岐切断と主値を示し、スカラーの能力を 3 つのトレイトに分けた理由を説明します。

設計目標

Complex[Double] に対して $\mathbb C$ の初等関数を提供します。文書化された分岐切断上での主値、Double の全範囲で数値的に安定な式、そして重要な箇所での IEEE 特殊値の処理を備え、これらすべてをジェネリックな core から切り離します。

数学的背景

極形式

すべての $z = x + iy \ne 0$ は、$r = |z| = \sqrt{x^2 + y^2}$ として極形式 $z = r e^{i\theta}$ を持ちます。角度 $\theta$ は $2\pi$ を法として定まります。主偏角 $\operatorname{Arg} z$ は $(-\pi, \pi]$ に入る代表元で、負の実軸以外では $\operatorname{atan2}(y, x)$ に等しくなります。

主値の対数とその分岐切断

$e^w = z$ の解は $w = \ln|z| + i(\operatorname{Arg} z + 2k\pi)$ です。主値の対数は $k = 0$ をとります:

$$\operatorname{Log} z = \ln|z| + i\operatorname{Arg} z .$$

$\operatorname{Arg}$ は負の実軸をまたぐと $2\pi$ だけ跳ぶため、$\operatorname{Log}$ は $\mathbb C \setminus (-\infty, 0]$ 上で解析的であり、$(-\infty, 0]$ がその分岐切断です。他の多価関数はすべて $\operatorname{Log}$ を通じて定義され、分岐切断をそこから受け継ぎます。¹1 W. Kahan, “Branch cuts for complex elementary functions, or much ado about nothing’s sign bit”, in The State of the Art in Numerical Analysis, Clarendon Press, 1987. 符号付きゼロによって分岐切断のどちら側かを選べることについても論じています。

他の関数の主値

$$\begin{aligned} \sqrt z &= e^{\frac12 \operatorname{Log} z}, & z^w &= e^{w \operatorname{Log} z}, \\ \operatorname{asin} z &= -i\operatorname{Log}\big(iz + \sqrt{1 - z^2}\big), & \operatorname{acos} z &= \tfrac{\pi}{2} - \operatorname{asin} z, \\ \operatorname{atan} z &= \tfrac{i}{2}\big(\operatorname{Log}(1 - iz) - \operatorname{Log}(1 + iz)\big), & \operatorname{asinh} z &= -i\operatorname{asin}(iz), \\ \operatorname{acosh} z &= 2\operatorname{Log}\Big(\sqrt{\tfrac{z + 1}{2}} + \sqrt{\tfrac{z - 1}{2}}\Big), & \operatorname{atanh} z &= \tfrac12\big(\operatorname{Log}(1 + z) - \operatorname{Log}(1 - z)\big). \end{aligned}$$

関数	分岐切断	主値の範囲
log	$(-\infty, 0]$	$\operatorname{Im} \in (-\pi, \pi]$
sqrt	$(-\infty, 0)$	$\operatorname{Re} \ge 0$
pow	$z$ について $(-\infty, 0]$($w$ が整数でない場合)	$\operatorname{Log}$ から
asin	$(-\infty, -1) \cup (1, \infty)$	$\operatorname{Re} \in [-\pi/2, \pi/2]$
acos	$(-\infty, -1) \cup (1, \infty)$	$\operatorname{Re} \in [0, \pi]$
atan	$i(-\infty, -1) \cup i(1, \infty)$	$\operatorname{Re} \in [-\pi/2, \pi/2]$
asinh	$i(-\infty, -1) \cup i(1, \infty)$	$\operatorname{Im} \in [-\pi/2, \pi/2]$
acosh	$(-\infty, 1)$	$\operatorname{Re} \ge 0$, $\operatorname{Im} \in [-\pi, \pi]$
atanh	$(-\infty, -1) \cup (1, \infty)$	$\operatorname{Im} \in [-\pi/2, \pi/2]$

逆数系の関数は合成です: $\sec z = 1/\cos z$、$\operatorname{asec} z = \operatorname{acos}(1/z)$ など。

実部と虚部

三角関数と双曲線関数は、$\cos(iy) = \cosh y$ と $\sin(iy) = i\sinh y$ を用いた加法定理によって分解されます:

$$\begin{aligned} \sin(x + iy) &= \sin x\cosh y + i\cos x\sinh y, & \cos(x + iy) &= \cos x\cosh y - i\sin x\sinh y, \\ \sinh(x + iy) &= \sinh x\cos y + i\cosh x\sin y, & \cosh(x + iy) &= \cosh x\cos y + i\sinh x\sin y . \end{aligned}$$

設計上の決定

3 つの能力トレイト

問題。 アルゴリズムは実スカラーに対して luna-generic や arithmetic のトレイト以上のものを必要としますが、すべてのスカラー型があらゆる追加能力を持つわけではありません。

選択肢。 1 つの大きな「浮動小数点実数」トレイト、トレイトを使わない(Double をハードコードする)、階層化されたトレイト群。

選択。 それぞれ別々の関心事を扱う 3 つのトレイト:

• FloatingAnalyticScalar は代数的・解析的な能力(Num、Compare、定数、実数の初等関数とその逆関数を備えた Field)に名前を付けます。IEEE 754 については何も仮定しないため、厳密な型や区間型でも満たすことができます。
• FloatingSpecialValues は IEEE 754 の層に名前を付けます: NaN、符号付き無限大、ゼロの符号です。これらを持たない型にとっては意味がないため、解析トレイトには含めていません。
• FloatingBackendScalar は両者を組み合わせ、以下の安定なアルゴリズムが使うプリミティブを追加します: hypot(オーバーフローしない絶対値)、log1p(1 付近の対数)、trunc と to_int(pow での整数の指数の検出)、from_double(定数)。

これにより、コードは自身の要件を表す最小のトレイトを要求できます。Luna Flow は単一の「実数」トレイトよりもこの方式を好みます。3 つとも Float と Double に対して実装されています。公開関数は依然として Complex[Double] 向けに書かれており、Double のプリミティブを直接呼び出します。これらのトレイトについてジェネリックな関数はまだありません。

自由関数

MoonBit では、パッケージが他のパッケージで定義された型にメソッドやトレイトインスタンスを追加できないため、解析関数を Complex[T] のメソッドにすることはできません。これらは @fb.log(z) のような自由関数であり、ジェネリックなコアは浮動小数点のセマンティクスから切り離されたままです。

オーバーフローしない絶対値

直接計算した $\sqrt{x^2 + y^2}$ は、$|z|$ が表現可能であっても $|x| > 1.34 \times 10^{154}$ でオーバーフローし、非常に小さい入力ではアンダーフローします。$m = \max(|x|, |y|)$、$t = \min(|x|, |y|)/m \in [0, 1]$ とすると:

$$|z| = m\sqrt{1 + t^2}, \qquad \ln|z| = \ln m + \tfrac12 \ln(1 + t^2), \qquad |z|^2 = m^2\big((x/m)^2 + (y/m)^2\big).$$

abs はこの種のスケーリングを使う hypot に委譲します。abs_log は 2 番目の形を使うため、ゼロでない有限の $z$ すべてに対して $\ln|z|$ は有限です。abs_sqr は 3 番目の形を使い、$|z|^2$ 自体がオーバーフローする場合にのみオーバーフローします。

安定な平方根

教科書どおりの式

$$\operatorname{Re}\sqrt z = \sqrt{\frac{|z| + x}{2}}, \qquad \operatorname{Im}\sqrt z = \operatorname{sign}(y)\sqrt{\frac{|z| - x}{2}}$$

は壊滅的な桁落ちを起こします: $x > 0$ かつ $|y| \ll x$ では $|z| - x$ がすべての桁を失い、$x < 0$ では $|z| + x$ が同様になります。このパッケージは桁落ちのない量だけを計算します

$$w = \sqrt{\frac{|x| + |z|}{2}} = \begin{cases} \sqrt{|x|}\,\sqrt{\tfrac12\big(1 + \sqrt{1 + t^2}\big)}, & |x| \ge |y|,\ t = |y|/|x|, \\ \sqrt{|y|}\,\sqrt{\tfrac12\big(t + \sqrt{1 + t^2}\big)}, & |x| < |y|,\ t = |x|/|y|, \end{cases}$$

ここでスケーリングした形はオーバーフローを避けます。もう一方の部分は $\operatorname{Re}\sqrt z \cdot \operatorname{Im}\sqrt z = y/2$ から復元します。$x \ge 0$ では根は $w + \frac{y}{2w}i$ です。実際、$w^2 = (x + |z|)/2$ とすると、

$$\Big(w + \frac{y}{2w}i\Big)^2 = w^2 - \frac{y^2}{4w^2} + yi = \frac{(x + |z|)^2 - y^2}{2(x + |z|)} + yi = \frac{2x^2 + 2x|z|}{2(x + |z|)} + yi = x + yi .$$

$x < 0$ では、$|x|$ を使った同じ計算により、根は $y$ の符号に応じて $\frac{|y|}{2w} \pm wi$ です。主枝が要求するとおり、実部が負になることはありません。分岐切断上では $y = -0$ は $y = +0$ と同様に扱われるため、両側とも $+\sqrt{|x|}\,i$ に写ります。

Smith の除算

教科書どおりの商は $c^2 + d^2$ で割るため、core の設計 で述べたオーバーフローの問題があります。Smith の方法は大きい方の成分で先に割ります。²2 R. L. Smith, “Algorithm 116: Complex division”, Communications of the ACM 5(8), 1962. $|c| \ge |d|$ のとき $r = d/c$ とすると $|r| \le 1$ です:

$$\frac{a + bi}{c + di} = \frac{(a + bi)(c - di)}{c^2 + d^2} = \frac{(a + bi)(1 - ri)}{c(1 + r^2)} = \frac{(a + br) + (b - ar)i}{c(1 + r^2)} ,$$

$|d| > |c|$ のときは $r = c/d$ として対称的に計算します。2 乗されるのは $r^2 \le 1$ だけなので、分母がオーバーフローするのは商がオーバーフローする場合に限ります。div は $1/c$ を一度だけ計算し、それを掛けます。$w$ に無限大の部分があり、$z$ に NaN の部分がない場合、div は C99 Annex G に従います: 有限の分子はゼロの結果を与え、無限大の分子は無限大の符号の商を与えます。

厳密な整数べき

$e^{n\operatorname{Log} z}$ による $z^n$ は角度 $n\theta$ と絶対値 $e^{n\ln|z|}$ を丸めるため、$(1 + i)^2$ でさえちょうど $2i$ にはなりません。$|n| \le 2^{31} - 1$ の実整数の指数に対しては、pow と pow_real は代わりに二進累乗法を使います: $O(\log |n|)$ 回の複素乗算で、小さなガウス整数に対しては厳密であり、$z^{-n} = (z^{-1})^n$ です。その他の指数では $e^{w\operatorname{Log} z}$ の極形式を使います: abs_log による $\ell = \ln|z|$ と $\theta = \arg z$ を用いて、

$$z^w = e^{(u + iv)(\ell + i\theta)} = e^{u\ell - v\theta}\big(\cos(u\theta + v\ell) + i\sin(u\theta + v\ell)\big), \qquad w = u + iv .$$

$z = 0$ では、$0^0 = 1$、実数 $w > 0$ に対して $0^w = 0$ となり、それ以外の指数はすべて NaN を与えます。

オーバーフローしない正接

$\cos^2 x + \sinh^2 y = \frac12(\cos 2x + \cosh 2y)$ を用いると $\tan(x + iy) = \dfrac{\sin 2x + i\sinh 2y}{\cos 2x + \cosh 2y}$ です。$|y|$ が大きいと、$\tan z \to \pm i$ であるにもかかわらず $\sinh 2y$ と $\cosh 2y$ がともにオーバーフローします。$d = e^{-2|y|}$ として分子と分母に $2d$ を掛け、$2d\cosh 2y = 1 + d^2$、$2d\sinh 2y = \operatorname{sign}(y)(1 - d^2)$ を使うと:

$$\tan(x + iy) = \frac{2d\sin 2x + i\operatorname{sign}(y)(1 - d^2)}{1 + d^2 + 2d\cos 2x} .$$

tan は $|y| \ge 1$ でこの形を使い、それ未満では精度のよい直接の形を使います。tanh は $x$ と $y$ の役割を入れ替えて同じ導出を適用します。

Hull、Fairgrieve、Tang による逆正弦アルゴリズム

$x, y \ge 0$ に対し、$r = |z + 1|$、$s = |z - 1|$、$A = \frac{r + s}{2} \ge 1$、$B = x/A \le 1$ とします。このとき³3 T. E. Hull, T. F. Fairgrieve and P. T. P. Tang, “Implementing the complex arcsine and arccosine functions using exception handling”, ACM Transactions on Mathematical Software 23(3), 1997. 切り替えの値 $1.5$ と $0.6417$ は同論文によるものです。

$$\operatorname{asin} z = \arcsin B + i\ln\big(A + \sqrt{A^2 - 1}\big), \qquad \operatorname{acos} z = \arccos B - i\ln\big(A + \sqrt{A^2 - 1}\big),$$

であり、他の象限は各部分の奇関数性から得られます。どちらの式も 2 つの領域で精度を失うため、アルゴリズムはそれらを別扱いします:

• $B$ が $1$ に近い($B > 0.6417$)とき、$\arcsin B$ は悪条件です。実部は、$r + x + 1$ と $s \pm (1 - x)$ から減算なしで作った量 $D$ を用いて $\arctan(x/\sqrt D)$ として計算します。
• $A$ が $1$ に近い($A \le 1.5$)とき、$\ln(A + \sqrt{A^2 - 1})$ は $A^2 - 1$ の影響を受けます。アルゴリズムは $y^2/(r + x + 1)$ と $s \pm (1 - x)$ から桁落ちなしに $A - 1$ を計算し、$\operatorname{log1p}\big((A - 1) + \sqrt{(A - 1)(A + 1)}\big)$ を使います。

$|x|$ または $|y|$ が $10^{150}$ を超え、$r$ と $s$ がオーバーフローする場合、asin は $\sqrt{1 - z^2} \approx -iz$ から得られる漸近形 $\operatorname{asin} z \approx \operatorname{atan2}(x, y) + i(\ln 2 + \ln|z|)$ を使います。そこでは acos は $\pi/2 - \operatorname{asin} z$ を使います。

逆正接と逆双曲線正接

$z = x + iy$ として $q = (1 + iz)/(1 - iz)$ と書くと:

$$q = \frac{(1 - y) + ix}{(1 + y) - ix}, \qquad |q|^2 = \frac{x^2 + (1 - y)^2}{x^2 + (1 + y)^2}, \qquad \arg q = \operatorname{atan2}\big(2x,\ 1 - x^2 - y^2\big),$$

となり、$\operatorname{atan} z = \frac{1}{2i}\operatorname{Log} q$ から次が得られます

$$\operatorname{atan} z = \tfrac12\operatorname{atan2}(2x, 1 - x^2 - y^2) + \tfrac{i}{4}\ln\frac{x^2 + (1 + y)^2}{x^2 + (1 - y)^2} .$$

対数の中の比は、$u = 2y/(1 + |z|^2)$ として $(1 + u)/(1 - u)$ に等しくなります。$|u| < 0.1$ のとき、実軸付近での精度を保つため、このパッケージはこれを $\operatorname{log1p}(u) - \operatorname{log1p}(-u)$ として評価します。$|z|$ が大きいときは、$m = \max(|x|, |y|)$ として実部を $x/m, y/m$ 上で評価します。atanh は $\operatorname{atanh} z = -i\operatorname{atan}(iz)$ から、$x$ と $y$ を入れ替えた同じ計算であり、$|4x/((1 - x)^2 + y^2)| < 1/4$ のときに $\operatorname{log1p}$ に切り替えます。

Kahan の式による逆双曲線余弦

acosh は安定な平方根から組み立てた $2\operatorname{Log}\big(\sqrt{(z + 1)/2} + \sqrt{(z - 1)/2}\big)$ を使います。$\operatorname{Log}(z + \sqrt{z^2 - 1})$ とは異なり、各平方根の分岐切断はその引数が負の実数となる場所にあるため、追加の符号調整なしで正しい分岐切断 $(-\infty, 1)$ を持ちます。

厳密な実数の高速パス

sin と cos は $y = 0$ ちょうどのとき実関数の値を返します。asin、acos、asinh、acosh、atanh は、実数の入力を _real 関数へ、純虚数の入力を $y$ の実逆関数へ振り分けます。これにより実数直線上の結果は厳密に実数のままとなり、_real 関数は Double から始める呼び出し側にも役立ちます。

コアの逆数を通じた逆数系の関数

sec、csc、cot、それらの双曲線版、および逆数系の逆関数は、core の Complex::inv との合成です。そのスケーリングしない式と、絶対値がゼロのときの中断を受け継ぎ、極で無限大を返すことはありません。acot は $z = 0$ を特別扱いして $\pi/2$ を返します。

正しさと不変条件

テストスイートで確認している恒等式

第 1 象限のサンプル点における exp(log z) = z、sin(asin z) = z、cos(acos z) = z、tan(atan z) = z、sinh(asinh z) = z、tanh(atanh z) = z(許容誤差 $10^{-10}$ ～ $10^{-12}$)。sqrt(-3 + 4i) = 1 + 2i。非常に大きな asin の引数、atan の分岐切断、無限大の atanh と acosh の入力、無限大による除算、$0^0$、$0^{-1}$ に対する回帰テスト。

主値からの既知の逸脱

実装は、主値では $\pi$ が必要な箇所で定数 $2\pi$(tau)を使っています:

関数と入力	返される値	主値
arg(z), $y = \pm 0$, $x < 0$	$2\pi$	$\pi$(C99: $y$ の符号により $\pm\pi$)
log(z)、同じ入力	$\ln\lvert x\rvert + 2\pi i$	$\ln\lvert x\rvert + \pi i$
pow、pow_real、負の実数の底、整数でない指数	角度 $2\pi$	角度 $\pi$
acos(z), $x < 0$, $y \ne 0$	$2\pi - \rho + \dots$	$\pi - \rho + \dots$
acos_real(x), $x < -1$	$2\pi - i\operatorname{acosh}\lvert x\rvert$	$\pi - i\operatorname{acosh}\lvert x\rvert$
asec_real(x), $-1 < x < 0$	$2\pi - i\operatorname{acosh}\lvert 1/x\rvert$	$\pi - i\operatorname{acosh}\lvert 1/x\rvert$
acosh_real(x), $x < -1$	$\operatorname{acosh}\lvert x\rvert + 2\pi i$	$\operatorname{acosh}\lvert x\rvert + \pi i$

$e^{2\pi i} = 1$ であるのに対し $e^{\pi i} = -1$ なので、これらの値は入力の対数や逆関数の値ですらありません: exp(log(-1)) は $1$ であり、$x < 0$ では cos(acos(z)) は $-z$ です。テストスイートは現在、負の実軸上での arg の値 $2\pi$ を検証しています。これらは実装で修正すべき欠陥であり、このマニュアルは現在の振る舞いを記録しています。

その他の精度に関する注意

• asin は $A \le 1.5$ の分岐で $\sqrt{A^2 - 1}$ を評価しますが、Hull–Fairgrieve–Tang アルゴリズム(およびここでの acos)は $\sqrt{(A - 1)(A + 1)}$ を使います。分岐点 $\pm1$ の近くではこれにより桁が失われます: $1 + 10^{-10}i$ では虚部の相対誤差が $4 \times 10^{-8}$ 程度になります。
• exp は $\cos y$ と $\sin y$ を掛ける前に $e^x$ でオーバーフローするため、exp(710 + 0i) の虚部は NaN になります($\infty \cdot 0$)。
• log1p の Float インスタンスは $\ln(1 + x)$ であり、$|x| \ll 1$ で相対精度が失われます。これを使う公開関数はまだありません。

採用しなかった代替案

• Complex のメソッド。 別パッケージからは不可能であり、関数をコアに移すとジェネリック型に IEEE のセマンティクスが持ち込まれます。
• 単一の浮動小数点実数トレイト。 あらゆる解析スカラーに IEEE の特殊値を強制することになります。
• どこでも教科書どおりの式。 より単純ですが、上で導出したとおり $|z| \gtrsim 10^{154}$ でオーバーフローし、分岐切断の近くで桁落ちを起こします。
• C99 Annex G の特殊値表の完全な実装。 API に記載したケースだけを処理し、残りの無限大と NaN は通常の算術を通じて伝播します。

範囲外

• 解析関数は Complex[Double] 向けのみです。トレイトは Float に対して実装されていますが、Complex[Float] 向けの関数はありません。
• arg を除き、符号付きゼロによる分岐切断の側の選択は行いません。
• 逆数系の関数は極で無限大を返すのではなく中断します。
• 誤差限界の保証はありません。精度は回帰テストで確認しています。
• 検査付きや文脈付き(Result)の版はありません。

Footnotes

1. W. Kahan, “Branch cuts for complex elementary functions, or much ado about nothing’s sign bit”, in The State of the Art in Numerical Analysis, Clarendon Press, 1987. 符号付きゼロによって分岐切断のどちら側かを選べることについても論じています。 ↩

2. R. L. Smith, “Algorithm 116: Complex division”, Communications of the ACM 5(8), 1962. ↩

3. T. E. Hull, T. F. Fairgrieve and P. T. P. Tang, “Implementing the complex arcsine and arccosine functions using exception handling”, ACM Transactions on Mathematical Software 23(3), 1997. 切り替えの値 $1.5$ と $0.6417$ は同論文によるものです。 ↩