Float 三角関数アルゴリズム

概要

calx の Float は以下の 10 種の三角関数・双曲線関数を提供する:

基本戦略は全関数に共通する三段構成である:

1. 引数縮約 (argument reduction): 入力を Taylor 級数の収束が速い小さな範囲に縮小する
2. Taylor 級数: Paterson-Stockmeyer 評価で級数を計算する
3. 倍角復元 (doubling): 倍角公式を反復適用して元のスケールに戻す

入力精度に応じてディスパッチを行い、低精度では高速なファストパスを使用する:

逆三角関数は atan に帰着させる:

• $\operatorname{asin}(x) = \operatorname{atan}\!\bigl(x / \sqrt{1 - x^2}\bigr)$
• $\operatorname{acos}(x) = \pi/2 - \operatorname{asin}(x)$

sin・cos

引数縮約

入力 $x$ を $[0, \pi/2]$ の範囲に縮約する:

1. $2\pi$ による縮約: $x \bmod 2\pi$ で $[0, 2\pi)$ に収める
2. 象限判定: $[0, 2\pi)$ を 4 象限に分割し、符号と関数の対応を決定する
3. $[0, \pi/2]$ への縮小: 対称性を利用して $[0, \pi/2]$ に射影する

さらに象限に応じて sin Taylor と cos Taylor のどちらを使うか選択する。 第 1・第 3 象限では sin の Taylor 級数、 第 2・第 4 象限では cos の Taylor 級数に帰着できる。

Taylor 級数

$u = x^2$ として:

$$\sin(x) = x \cdot \sum_{k=0}^{\infty} \frac{(-1)^k \, u^k}{(2k+1)!} = x \left(1 - \frac{u}{3!} + \frac{u^2}{5!} - \cdots \right)$$

$$\cos(x) = \sum_{k=0}^{\infty} \frac{(-1)^k \, u^k}{(2k)!} = 1 - \frac{u}{2!} + \frac{u^2}{4!} - \cdots$$

Paterson-Stockmeyer 評価

$n$ 次多項式を通常の Horner 法で評価すると $O(n)$ 回の乗算が必要だが、 Paterson-Stockmeyer 法では $O(\sqrt{n})$ 回に削減できる。

多項式 $P(u) = c_0 + c_1 u + c_2 u^2 + \cdots + c_n u^n$ を ブロックサイズ $l$ で分割する。まず $u, u^2, \ldots, u^{l-1}$ のべき乗テーブルを構築し、 $u^l$ を「大ステップ」の基底として Horner 的に評価する:

$$P(u) = \bigl((\cdots(B_{s} \cdot u^l + B_{s-1}) \cdot u^l + B_{s-2}) \cdots\bigr) \cdot u^l + B_0$$

各ブロック $B_i = c_{il} + c_{il+1} u + \cdots + c_{il+l-1} u^{l-1}$ は 事前計算済みのべき乗テーブルとの内積で $O(l)$ 回の乗算で評価できる。 全体で $O(l + n/l) = O(\sqrt{n})$ 回 ($l = \sqrt{n}$ のとき最小) の乗算となる。

cos 経由の sin 計算 (高精度パス)

sin を直接 Taylor 展開する代わりに、以下の手順で計算する:

1. cosDoubling で $\cos(\theta)$ を計算する (Taylor 1 本)
2. $\sin(\theta) = \sqrt{1 - \cos^2(\theta)}$ で sin を求める

Taylor 級数の評価が 1 本分で済むため、sin と cos を独立に展開するより高速である。

倍角復元

引数縮約で $K$ 回半減した角度 $\theta / 2^K$ に対して Taylor 級数を評価した後、 倍角公式を $K$ 回適用して元のスケールに復元する。

倍角公式

$$\sin(2\theta) = 2 \sin(\theta) \cos(\theta)$$

$$\cos(2\theta) = 2\cos^2(\theta) - 1$$

cosDoubling

cos のみを倍角復元する場合:

1. 引数を $K$ 回半減: $\theta' = \theta / 2^K$
2. cosTaylor で $\cos(\theta')$ を計算する
3. $K$ 回反復: $c \leftarrow 2c^2 - 1$

各ステップは二乗 1 回 + 定数演算のみで、乗算は mpn::square の対称性最適化が効く。

sinCosDoubling

$(\sin, \cos)$ の組を同時に倍角復元する場合:

1. 引数を $K$ 回半減: $\theta' = \theta / 2^K$
2. sinTaylor と cosTaylor で $(\sin(\theta'), \cos(\theta'))$ を計算する
3. $K$ 回反復: $s' \leftarrow 2sc$、$c' \leftarrow 2c^2 - 1$、$(s, c) \leftarrow (s', c')$

Q128/Q256 ファストパス

sincos_small ($\leq 19$ 桁)

Float オブジェクトの構築を完全にバイパスし、Q128 固定小数点で計算する。

1. 引数縮約: $z = x - k \cdot (\pi/2)$ で $|z| < \pi/4$ に縮小する
2. 極小値の処理: $|z| < 2^{-64}$ のとき $\sin(z) \approx z$、$\cos(z) \approx 1$ として即座に返す
3. versine による桁落ち回避: $\cos(z)$ を直接計算する代わりに、versine $v$ を計算する: $$v = 1 - \cos(z) = \frac{z^2}{2!} - \frac{z^4}{4!} + \frac{z^6}{6!} - \cdots$$ $|z| < \pi/4$ では $\cos(z) > 1/\sqrt{2} \approx 0.707$ なので、 $\cos(z) = 1 - v$ は常に 2 ワードで正確に表現できる
4. Q128 Taylor: sin と versine の各係数を Q128 (128 ビット固定小数点) で評価する

38 桁パス

$\leq 38$ 桁の精度では、軽量な Float で sinTaylor / cosTaylor を直接呼び出す。 Paterson-Stockmeyer のオーバーヘッドが不要な項数で収束するため、 単純な Horner 評価で十分である。

atan

引数縮約

atan の引数縮約は 2 段階で行う:

1. 逆数変換: $|x| > 1$ のとき $$\operatorname{atan}(x) = \frac{\pi}{2} - \operatorname{atan}\!\left(\frac{1}{x}\right)$$ により $|x| < 1$ に縮小する
2. 引数半減: 以下の恒等式を $K$ 回適用し、引数をさらに小さくする: $$\operatorname{atan}(x) = 2 \cdot \operatorname{atan}\!\left(\frac{x}{1 + \sqrt{1 + x^2}}\right)$$

Taylor 級数

$u = x^2$ として:

$$\operatorname{atan}(x) = x \cdot \sum_{k=0}^{\infty} \frac{(-1)^k \, u^k}{2k+1} = x \left(1 - \frac{u}{3} + \frac{u^2}{5} - \cdots \right)$$

sin/cos と同じ構造の奇数次級数であり、Paterson-Stockmeyer 評価をそのまま適用できる。

Q128 ファストパス ($\leq 19$ 桁)

std::atan(double) で $\sim$53 ビットの初期近似 $y_0$ を得た後、 Newton 補正を 1 回適用する:

$$y_1 = y_0 + \frac{x \cos(y_0) - \sin(y_0)}{\cos(y_0) + x \sin(y_0)}$$

$\sin(y_0)$ と $\cos(y_0)$ は Q128 Taylor + 倍角復元で同時に計算する。 Newton 1 回で精度が倍増し、19 桁に到達する。

Q256 ファストパス ($\leq 38$ 桁)

Q128 パスと同様の構造で、Q256 sincos + Newton 補正により 38 桁の精度を得る。

逆三角関数

逆三角関数は atan に帰着させる。

asin

$$\operatorname{asin}(x) = \operatorname{atan}\!\left(\frac{x}{\sqrt{1 - x^2}}\right)$$

特殊値: $\operatorname{asin}(\pm 1) = \pm\pi/2$。

$|x|$ が 1 に近いとき $\sqrt{1 - x^2}$ の計算で桁落ちが生じるため、 ガードビットを追加して精度を確保する。

acos

$$\operatorname{acos}(x) = \frac{\pi}{2} - \operatorname{asin}(x)$$

特殊値: $\operatorname{acos}(1) = 0$、$\operatorname{acos}(-1) = \pi$。

atan2

$\operatorname{atan2}(y, x)$ は 4 象限の逆正接関数である。 象限を判定した後、$\operatorname{atan}(y/x)$ に帰着させる:

Q128/Q256 ファストパスも同様に適用される。

双曲線関数

Taylor 級数

双曲線関数の Taylor 級数は三角関数の交互符号を除いたものに一致する:

$$\sinh(x) = x + \frac{x^3}{3!} + \frac{x^5}{5!} + \cdots = x \cdot \sum_{k=0}^{\infty} \frac{u^k}{(2k+1)!}, \quad u = x^2$$

$$\cosh(x) = 1 + \frac{x^2}{2!} + \frac{x^4}{4!} + \cdots = \sum_{k=0}^{\infty} \frac{u^k}{(2k)!}, \quad u = x^2$$

符号が交互にならないため、すべての項が正であり、桁落ちの問題は生じない。

sinhCoshDoubling

三角関数と同様に引数半減 + 倍角復元を行う。倍角公式:

$$\sinh(2t) = 2 \sinh(t) \cosh(t)$$

$$\cosh(2t) = 2\cosh^2(t) - 1$$

$(\sinh, \cosh)$ の組を $K$ 回同時に倍角復元する。 構造は sinCosDoubling と同一である。

精度に応じたディスパッチ

• 19 桁以下: Q128 固定小数点の exp で高速に計算できる。
• 20 桁以上: sinhCoshDoubling (Taylor + 倍角復元) を使用する。 $\sinh = (e^x - e^{-x})/2$ の exp ベース方式は $e^{-x}$ の計算に Float 除算 (Newton 逆数反復) が必要であり、sinhCoshDoubling より遅い。

tanh は $\tanh(x) = \sinh(x) / \cosh(x)$ として計算する。 sinhCoshDoubling で $(\sinh, \cosh)$ の組を同時に得られるため、 除算 1 回のコストのみが追加される。

設計判断

Taylor 級数 vs AGM ベース

三角関数の計算には AGM (Arithmetic-Geometric Mean) ベースの手法も知られているが、 calx では Taylor 級数 + 引数縮約/倍角復元を一貫して採用している。 その理由は以下の通りである:

• 実用精度範囲での性能: AGM は漸近的には $O(M(n) \log n)$ ($M(n)$ は $n$ ビット乗算) だが、 定数因子が大きい。実用的な精度範囲 (数百〜数万桁) では Taylor + Paterson-Stockmeyer の方が高速である
• コードの共有: sin, cos, atan, sinh, cosh はすべて同じ構造 (引数縮約 → Paterson-Stockmeyer Taylor → 倍角復元) を持つため、 内部のべき乗テーブル構築・多項式評価ルーチンを共有できる
• ファストパスとの一貫性: Q128/Q256 ファストパスも Taylor 級数を使用しており、 精度閾値を跨いだとき結果の連続性が自然に保証される

参考文献

• Brent, R. P. and Zimmermann, P. (2010). Modern Computer Arithmetic. Cambridge University Press. Chapter 4: Elementary and Special Functions.
• Muller, J.-M. et al. (2018). Handbook of Floating-Point Arithmetic. 2nd ed. Birkhäuser. Chapter 11: Evaluating Floating-Point Elementary Functions.
• Paterson, M. S. and Stockmeyer, L. J. (1973). "On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials". SIAM J. Comput., 2(1), 60–66.