Strassen の 3 実数乗算とは何ですか？

複素数の積 (a+bi)(c+di) は素朴には 4 回の実数乗算と 2 回の加減算で計算するが、3 回の乗算と 5 回の加減算で同じ結果が得られる。これは Karatsuba 法を deg=1 多項式 a+bi に適用した特殊形で、sangi は対称な形 Re = ac - bd, Im = (a+b)(c+d) - ac - bd を用いる。実数乗算のコストが加減算より著しく大きい多倍長型 (Int, Float) では正味の高速化に繋がる。

Smith の除算法はなぜオーバーフローを防げるのですか？

素朴な (a+bi)/(c+di) = ((ac+bd)+(bc-ad)i) / (c²+d²) では中間量 c²+d² が表現範囲を超えてオーバーフローしうる。Smith の方法は |c| と |d| の大小で場合分けし、小さい方を大きい方で正規化することで、中間値の絶対値を入力の絶対値以下に抑える。

なぜ abs(z) は std::hypot を呼ぶのですか？

sqrt(a*a + b*b) は a*a または b*b の段階でオーバーフロー・アンダーフローしうる。std::hypot は a と b を内部で適切に正規化して |z| = sqrt(a²+b²) を計算するため、入力が表現可能であれば結果も表現可能となる。多倍長型では内部 abs()/normSq() を組み合わせた等価な処理を行う。

極形式 polar(r, θ) の内部実装は？

polar(r, θ) は r·cos(θ) + i·r·sin(θ) を計算して Complex を構築する。多倍長型では cos/sin のキャッシュ済み定数を使わず引数 θ ごとに計算するが、頻出する θ=2πk/N の形 (FFT 回転因子) には expI(ω) という専用関数が用意されている。

Complex 四則演算と内部表現

概要

Complex<T> は実部 re と虚部 im を直接保持する単純な 2 成分構造体として実装されている。テンプレート引数 T は double, float, Int, Float, Rational 等の ComplexScalar コンセプトを満たす任意の型を受け取れる。

四則演算の実装は T の特性に応じて切り替えられる。実数乗算が加減算より著しく重い多倍長型では Strassen の 3 実数乗算を用い、除算における中間値のオーバーフロー対策として浮動小数点型では Smith 法を使用する。

API の使い方については Complex API リファレンス — 算術演算子を参照。

内部表現とメモリレイアウト

クラス本体は以下の 2 つの public メンバ変数のみを持つ。

template <ComplexScalar T>
struct Complex {
    T re;
    T im;
    // ...
};

仮想関数や追加メタデータは持たないため、$\operatorname{sizeof}(\texttt{Complex}\langle T\rangle) = 2 \operatorname{sizeof}(T)$ が成立する。以下の static_assert によりレイアウトの正しさがコンパイル時に検証される。

sizeof(Complex<T>) == 2 * sizeof(T)
offsetof(Complex<T>, re) == 0
offsetof(Complex<T>, im) == sizeof(T)

このレイアウト保証により、Complex<double> の配列を std::complex<double> や C99 _Complex double の配列として reinterpret_cast で再解釈できる。 FFTW などの外部 C ライブラリへゼロコピーで渡すための要件である。

関連記事: 複素数の定義と基本性質 / 複素平面とガウス平面

加減算

加減算は成分ごとに行われる単純な演算である。

$$(a+bi) \pm (c+di) = (a \pm c) + (b \pm d) i$$

$T$ ごとに 2 回の加算 (または減算) を発行するだけで完了する。 T = double ではコンパイラが SSE / AVX の packed double 加算 1 命令に畳み込む可能性がある。多倍長型 (Int, Float) では T 側の加減算がそのまま呼び出され、実部・虚部それぞれの桁数に応じて線形時間 ($O(\max(\text{words}(a), \text{words}(c)))$) で完了する。

計算量: 実数加減算 2 回。

スカラー混合演算 $z \pm s$ は実部にのみスカラーを作用させる: $$(a+bi) \pm s = (a \pm s) + b i$$ 虚部の T オブジェクトに対して代入・コピーが発生しないよう、実装では実部側のみ更新する z.re += s の経路を取る。

乗算 — Strassen 3 実数乗算

素朴な 4 実数乗算

定義どおりの実装は 4 回の実数乗算と 2 回の加減算からなる。

$$(a+bi)(c+di) = (ac - bd) + (ad + bc) i$$

$T$ が double / float の場合、sangi はこの素朴な式を採用する。 FP 乗算は FP 加減算と同程度のスループットであり、加減算を 2 回増やしてまで乗算を 1 回減らす利点が少ないためである。

Karatsuba 対称形 (3 実数乗算)

$T$ が多倍長型 (Int, Float, Rational) の場合、実数乗算が加減算より著しく重い (典型的には数十〜数百倍) ため、乗算を 1 回減らすほうが正味で有利になる。この場合 sangi は Karatsuba 法を 1 次多項式 $a + b X$ ($X = i$, $X^2 = -1$) に適用した対称形の 3 実数乗算を用いる:

$$\begin{aligned} ac &= a \cdot c, \\ bd &= b \cdot d, \\ (a+bi)(c+di) &= (ac - bd) \;+\; \bigl((a+b)(c+d) - ac - bd\bigr)\,i. \end{aligned}$$

実数乗算 3 回 ($ac$, $bd$, $(a+b)(c+d)$) と加減算 5 回で同じ結果が得られる。 Strassen の 3 乗算式 $k_1 = (a+b)c$, $k_2 = a(d-c)$, $k_3 = b(c+d)$ と数学的に等価であり、 sangi は加減算が対称に並ぶこちらの形を採用する。

整合性の確認

展開すると確かに一致する:

$$\begin{aligned} \mathrm{Re} &= ac - bd, \\ \mathrm{Im} &= (a+b)(c+d) - ac - bd = ac + ad + bc + bd - ac - bd = ad + bc. \end{aligned}$$

計算量: 実数乗算 3 回 + 加減算 5 回。実数乗算 $M(n)$、加減算 $A(n)$ について、コストは $3 M(n) + 5 A(n)$。多倍長 $A(n) = O(n)$、$M(n)$ は Karatsuba 域以上では超線形であるため、 $4 M(n) + 2 A(n)$ より小さい。

除算 — Smith 法

素朴な式とオーバーフロー

定義どおりの実装は次の通り。

$$\frac{a+bi}{c+di} = \frac{(ac + bd) + (bc - ad) i}{c^2 + d^2}$$

$T = $ double でこの式を直接使うと、$c^2 + d^2$ の段階で $|c|, |d|$ が大きい場合にオーバーフロー、小さい場合にアンダーフローする。たとえば $c = d = 10^{200}$ では $c^2 + d^2 = 2 \times 10^{400}$ が double の表現範囲外になる。

Smith の方法

$T$ が浮動小数点型のとき、sangi は Smith (1962) の方法を採用する。 $|c|$ と $|d|$ の大きい方で正規化することで、すべての中間値の絶対値を入力の絶対値以下に抑える。

$|c| \geq |d|$ の場合:

$$r = d / c, \quad t = c + d \, r, \quad \mathrm{re} = (a + b r) / t, \quad \mathrm{im} = (b - a r) / t.$$

$|c| < |d|$ の場合:

$$r = c / d, \quad t = d + c \, r, \quad \mathrm{re} = (a r + b) / t, \quad \mathrm{im} = (b r - a) / t.$$

いずれの場合も $|r| \leq 1$ であり、$t$ も $\max(|c|, |d|)$ のスケールに収まる。条件分岐 1 回と 4 回の除算 + 3 回の乗算 + 3 回の加減算で済む。

多倍長型の場合

$T = $ Int / Float / Rational ではオーバーフローの概念が異なる。 Int は任意精度、Float は指数部が広いため、 $c^2 + d^2$ がそのまま安全に計算できる。この場合 sangi は素朴な式 $\mathrm{denom} = c^2 + d^2$ を採用する。実装上は denom = normSq(w) を呼び、結果として実部 $(a c + b d)/\mathrm{denom}$ と虚部 $(b c - a d)/\mathrm{denom}$ を構成する。

$c + di = 0$ の場合は NaN を返し、例外送出やプログラム停止は行わない。

Smith, R. L. (1962). "Algorithm 116: Complex division". Communications of the ACM, 5 (8), 435.

絶対値・ノルム — hypot 安全計算

ノルムの二乗 normSq

ノルムの二乗は単純に

$$|z|^2 = a^2 + b^2$$

で計算する。FFT 等で頻出する用途であり、平方根を取らない分だけ高速である。

絶対値 abs

絶対値 $|z| = \sqrt{a^2 + b^2}$ は素朴に sqrt(a*a + b*b) と書くと、中間量 $a^2$ または $b^2$ の段階でオーバーフロー・アンダーフローしうる。

$T$ が浮動小数点型のとき、sangi は std::hypot(a, b) を呼び出す。 std::hypot は IEEE 754 に基づいて $|z|$ を計算する間に中間値が範囲外にならないことを保証する。典型的な実装は $m = \max(|a|, |b|)$ で正規化し、 $m \sqrt{(a/m)^2 + (b/m)^2}$ を返す。

$T$ が多倍長型のとき、sangi は同じ正規化戦略を多倍長 abs() / normSq() と sqrt() の組み合わせで実装する。 Float の指数部は数十億規模を扱えるため、典型的な数値計算で範囲外になることはほぼ無いが、分母 denom = c^2 + d^2 から平方根を取る経路を共通化することで一貫性を保つ。

偏角と共役

偏角 arg

偏角は atan2 で計算する。

$$\arg(z) = \operatorname{atan2}(b, a) \in (-\pi, \pi]$$

主値 (principal value) を採用しており、$z = 0$ の場合は atan2(0, 0) = 0 を返す。分枝切断は負の実軸 $a < 0, b = 0$ に置かれる。 $b \to 0^+$ では $\arg \to \pi$、$b \to 0^-$ では $\arg \to -\pi$ となる。

共役 conj

共役 $\bar{z} = a - bi$ は単に虚部の符号を反転するだけで $O(1)$。

$$|z|^2 = z \cdot \bar{z}$$

はノルムの二乗の代替計算となるが、sangi は対称性のある normSq を優先する (実部・虚部のキャンセルが生じない $a^2 + b^2$ のほうが数値的に安定するため)。

関連記事: 複素平面 — 絶対値と偏角 / 分枝切断

極形式と expI

polar(r, θ)

極形式から複素数を構築する関数。

$$\operatorname{polar}(r, \theta) = r \cos\theta + i \, r \sin\theta$$

$r$ と $\theta$ から実部 $r \cos\theta$、虚部 $r \sin\theta$ を計算して Complex<T> を直接構築する。多倍長型では T::sin / T::cos が呼び出され、指定精度に応じた級数展開・引数縮約 (range reduction) が内部で行われる。

expI(ω) — FFT 回転因子

FFT・NTT で頻出する形 $e^{i\omega} = \cos\omega + i \sin\omega$ を直接計算する関数。 $r = 1$ の polar と等価であるが、専用関数として用意することで r * cos、r * sin の乗算 (多倍長では非自明なコスト) を省略する。

$$\operatorname{expI}(\omega) = e^{i\omega} = \cos\omega + i \sin\omega$$

FFT で必要な回転因子 $W_N^k = e^{-2\pi i k/N}$ は次のように得られる:

for (int k = 0; k < N; ++k) {
    auto W_k = expI(-2.0 * M_PI * k / N);  // 単位円上の k-th 等分点
}

オイラー公式

$\omega = \pi$ では expI(M_PI) ≈ -1 + 0i となる。浮動小数点丸めにより虚部はわずかに残るが、多倍長型では T::pi(precision) のキャッシュ済み値と高精度三角関数により極めて高い精度で $-1$ に到達する。

関連記事: オイラーの公式 / 複素指数関数と三角関数

設計判断

浮動小数点と多倍長で乗算式を切り替える理由

浮動小数点演算では乗算と加減算のスループットがほぼ同じであり、現代の x86 / ARM では FMA (Fused Multiply-Add) によって乗算と加算が 1 命令に融合できる場合もある。この環境で Strassen を使うと、加減算 3 回増の代わりに乗算 1 回減という交換になり、正味でわずかに不利または同等となる。したがって T = double/float では素朴な 4 乗算式を採用している。

多倍長型では実数乗算のコストが加減算より $10^1$ – $10^3$ 倍大きく、乗算回数の削減が支配的に効く。この場合 Strassen の 3 乗算式が常に有利となる。

条件分岐コストの考慮

除算の Smith 法は $|c|$ と $|d|$ の比較で 1 回分岐するが、この分岐は分岐予測されやすく (片側に偏った入力分布が多い)、現代 CPU のパイプラインに与える影響は無視できる程度である。オーバーフローを起こす実装に比べれば、分岐コストは桁違いに小さい代償と言える。

NaN 安全性

$z / 0$ の場合は IEEE 754 NaN を返し、例外送出は行わない。 std::complex や Python cmath と同じ振る舞いを取る。多倍長型でも同様に NaN 伝播するため、計算チェーンの中でゼロ除算が発生してもプログラムが停止することはない。

参考文献

Smith, R. L. (1962). "Algorithm 116: Complex division". Communications of the ACM, 5 (8), 435.
Karatsuba, A. and Ofman, Yu. (1962). "Multiplication of Multidigit Numbers on Automata". Doklady Akademii Nauk SSSR, 145, 293–294.
Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997.