float_backend 设计

本页推导 float_backend 为 Complex[Double] 实现的公式，解释每个公式如何避免上溢、下溢和相消，给出分支切割线与主值，并解释为何将标量能力拆分为三个 trait。

设计目标

为 Complex[Double] 提供 $\mathbb C$ 的初等函数：在明确记录的分支切割线上取主值，在整个 Double 范围内使用数值稳定的公式，并在要紧之处处理 IEEE 特殊值，同时让这一切都不进入泛型核心。

数学背景

极坐标形式

每个 $z = x + iy \ne 0$ 都有极坐标形式 $z = r e^{i\theta}$，其中 $r = |z| = \sqrt{x^2 + y^2}$。角度 $\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}$$

设计决策

三个能力 trait

问题。 这些算法对实数标量的需求超出了 luna-generic 和 arithmetic trait，但并非每种标量类型都具备所有额外能力。

备选方案。 一个庞大的“浮点实数” trait；不用 trait（硬编码 Double）；一组分层的 trait。

选择。 三个 trait，各自负责一个关注点：

• FloatingAnalyticScalar 命名代数与解析能力（带有 Num、Compare、常量以及实数初等函数及其反函数的 Field）。它不对 IEEE 754 做任何假设，因此精确类型或区间类型也可能满足它。
• FloatingSpecialValues 命名 IEEE 754 这一层：NaN、带符号的无穷大以及零的符号。对于没有这些值的类型，它们没有意义，因此不属于解析 trait。
• FloatingBackendScalar 组合了二者，并加入下面稳定算法所用的基本操作：hypot（不会上溢的模）、log1p（接近一处的对数）、trunc 与 to_int（在 pow 中检测整数指数）以及 from_double（常量）。

这样代码就可以只要求能表达其需求的最小 trait，这正是 Luna Flow 相对于单一“实数” trait 所偏好的做法。三者都为 Float 和 Double 实现。公开函数目前仍针对 Complex[Double] 编写，并直接调用 Double 的基本操作；还没有对这些 trait 泛型的函数。

自由函数

MoonBit 不允许一个包为另一个包中定义的类型添加方法或 trait 实例，因此解析函数不能是 Complex[T] 的方法。它们是自由函数，如 @fb.log(z)，泛型核心也因此不含浮点语义。

不会上溢的模

直接计算 $\sqrt{x^2 + y^2}$ 时，当 $|x| > 1.34 \times 10^{154}$ 会上溢，尽管 $|z|$ 本身可以表示；对极小的输入则会下溢。取 $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 使用第二种形式，因此对每个有限非零的 $z$，$\ln|z|$ 都是有限的；abs_sqr 使用第三种形式，只有当 $|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|$）可得根为 $\frac{|y|}{2w} \pm wi$，符号与 $y$ 相同。实部永不为负，符合主分支的要求。在割线上，$y = -0$ 与 $y = +0$ 同样处理，因此两侧都映射到 $+\sqrt{|x|}\,i$。

Smith 除法

教科书式商要除以 $c^2 + d^2$，存在核心设计中所述的上溢问题。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$。只有 $r^2 \le 1$ 被平方，因此分母仅在商本身上溢时才会上溢。div 只计算一次 $1/c$ 再与之相乘。若 $w$ 有无穷分量而 $z$ 没有 NaN 分量，div 遵循 C99 附录 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|$ 很大时，$\sinh 2y$ 和 $\cosh 2y$ 都会上溢，而 $\tan z \to \pm i$。将分子分母同乘以 $2d$（$d = e^{-2|y|}$），并利用 $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),$$

其他象限可由各分量的奇性得出。这两个公式在两个区域会损失精度，算法对它们分别处理：

• 当 $B$ 接近 $1$（$B > 0.6417$）时，$\arcsin B$ 是病态的。此时实部按 $\arctan(x/\sqrt D)$ 计算，其中量 $D$ 由 $r + x + 1$ 和 $s \pm (1 - x)$ 构成，不涉及减法。
• 当 $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$。

反正切与反双曲正切

记 $q = (1 + iz)/(1 - iz)$，其中 $z = x + iy$：

$$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} .$$

对数内的比值等于 $(1 + u)/(1 - u)$，其中 $u = 2y/(1 + |z|^2)$；当 $|u| < 0.1$ 时，本包将其按 $\operatorname{log1p}(u) - \operatorname{log1p}(-u)$ 求值，以保持实轴附近的精度。当 $|z|$ 很大时，实部在 $x/m, y/m$ 上求值，其中 $m = \max(|x|, |y|)$。由 $\operatorname{atanh} z = -i\operatorname{atan}(iz)$，atanh 是交换 $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、它们的双曲对应函数以及反倒数类函数都是与核心的 Complex::inv 的复合。它们继承了其未缩放的公式以及在模为零时中止的行为，而不是在极点处产生无穷大。acot 对 $z = 0$ 做特殊处理，返回 $\pi/2$。

正确性与不变量

测试套件检查的恒等式

在第一象限的采样点上检查 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 语义带入泛型类型。
• 单一的浮点实数 trait。 这会把 IEEE 特殊值强加给每个解析标量。
• 处处使用教科书公式。 更简单，但如上文推导所示，它们在 $|z| \gtrsim 10^{154}$ 时上溢，并在分支切割线附近发生相消。
• 完整的 C99 附录 G 特殊值表。 只处理 API 中列出的情形；其余的无穷大和 NaN 通过普通算术传播。

边界

• 仅为 Complex[Double] 提供解析函数；虽然这些 trait 已为 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$ 取自该文。 ↩