【Rustで作る電卓】ランベルト（Lambert）の連分数展開でtan(正接)を実装 三角関数の実装(5) #22

今回は正接関数(tan)をフルスクラッチ実装していこうと思います。これまでsinやcosを実装してきており一見、これまでの延長線上でいけそうに見えますが、実は tan には特有の難しさがあり面白いところでもあります。

1.マクローリン展開は不向き？

前回まで扱ってきたsinとcosはマクローリン展開(テイラー展開 x=0)を用いて実装してきましたが、tanはマクローリン展開での実装が困難な特性を持っています。

これまでの sin や cos は、項を増やせば増やすほど精度が上がるマクローリン展開が非常に有効でした。しかし、tan のマクローリン展開を見てみましょう。導出はちょこっと面倒(煩雑)ですが、普通に\( \tan x= \sin x /cos x \)を微分していけば求まります。手元で計算してみると項の大変さが身をもって体験できます。

\[\tan x = \sum_{n=1}^{\infty} \frac{|B_{2n}| 2^{2n} (2^{2n} – 1)}{(2n)!} x^{2n-1} = x + \frac{1}{3}x^3 + \frac{2}{15}x^5 + \frac{17}{315}x^7 + \dots\]

ここで登場する \(B_{2n}\) はベルヌーイ数と呼ばれるもので、項が進むにつれて分子・分母が爆発的に巨大化し、規則性も複雑です。
(ベルヌーイ数はここでは詳しく紹介しませんが興味があれば、面白いので調べてみてください)

さらに致命的なのが収束半径です。tan は \(\pi/2\) で無限大に発散するため、マクローリン展開では \(\pi/2\) を超えるどころか、近づくことすら困難です。無理に計算しようとすると、とんでもない数の項数が必要になり、計算効率も精度も「絶望的」になります。

tanをマクローリン展開で近似するのは現実的でない理由

• 項数が増えるほど複雑　∵ベルヌーイ数が出てくる
• 収束半径が大きすぎる

2. Lambertの公式

マクローリン展開ではないtanの近似手法として挙げられるのが、18世紀の数学者ランベルト（Lambert）が導き出した連分数展開です。
\[\tan x = \cfrac{x}{1 – \cfrac{x^2}{3 – \cfrac{x^2}{5 – \cfrac{x^2}{7 – \dots}}}}\]

分母の中に分数が入り続けるこの形。一見複雑に見えますが、実はマクローリン展開よりも圧倒的に収束が早く、特に tan のような特異点（無限に飛ぶ点）を持つ関数と非常に相性が良いという特徴があります。

なぜ？この変な連分数がtanになるのという話：πが無理数であることの証明：ランベルトの功績と連分数展開の歴史【ランベルトの連分数展開 #1】

3. Rustでの反復アルゴリズムの実装

この連分数をコンピュータで計算する場合、「一番下の深い部分」から上に向かって計算していきます。f64精度であれば,12~13項目までやれば十分であるのでそのように実装してみます。

実装コード

/// ランベルトの連分数展開を用いた tan のカーネル関数
/// n: 反復回数 (項数)
fn kernel_tan(x: f64, n: usize) -> f64 {
    if x == 0.0 { return 0.0; }
    
    let x2 = x * x;
    let mut res = 0.0;

    // 下から上へ向かって反復計算
    // (2k + 1) の部分を n から 1 まで回す
    for i in (1..=n).rev() {
        let odd = (2 * i + 1) as f64;
        res = x2 / (odd - res);
    }

    x / (1.0 - res)
}

ループを逆順（rev()）に回すことで、深い階層から順に分母を確定させていく手法です。わずか10回程度の反復で、f64 の精度限界に近い値が得られます。
sinやcosでも触れましたが、\(|x|\)が非常に大きいときの精度を上げるため、payne_hanek法で範囲還元をしています。
範囲還元とは？：【Rustで作る電卓】sin関数実装 Rustで挑む三角関数の実装(1) マクローリン展開と範囲還元 #18

pub fn tan(x: f64) -> f64 {
    let is_negative = x < 0.0;
    let (r, q) = reduce_payne_hanek(x.abs());//範囲還元

    let val = kernel_tan(r, 13);

    let mut res = if q % 2 == 0 {
        val // q=0, 2 (0, pi 付近)
    } else {
        -1.0 / val // q=1, 3 (pi/2, 3pi/2 付近)
    };

    if is_negative { -res } else { res }
}

tanの実装例

4. 精度評価と結果

tan(3.14)の計算結果

実行結果

巨大な数に対しても、小数第14〜15位までの精度を維持できていることが確認できました。マクローリン展開ではこの精度を出そうとしたら50項とか計算しなければならないため、ここまで単純に見える連分数でここまでの精度が出せているのはとても興味深いです。
※123456.7は大きめの数、且つ、\(\pi/2\)で割るとほぼ0である(=極に近い)数です。

5. 速度比較と改良

今回の手法が十分な精度を出せるということはわかりましたが、速度面ではf64::tanと比べると約5倍強程度遅い(5ms vs 43ms)ということがわかりました。

この時間の内訳を見てみると　範囲還元のpayne_hanek法が一回当たり平均17ns、100万で17ms。残りがtanの計算とその他のヘッドなどで26msなっていることがわかりました。

payne_hanek法は比較的重い処理なのでこれに時間がかかるのは仕方ないため、tan周りを少し改良してみることにしました。

for文のunroll

簡単に実装できるのが、for文の展開です。プログラミングをしているとべた書きせず構文を使ってスマートに書くことが良いとされます。しかし、低レイヤーでの実装だと、明示的に何回処理するのか、for文のイテレーションをわざわざ保持させないような方式の方が早かったりします。

fn kernel_tan_unrolled(x: f64) -> f64 {
    if x == 0.0 { return 0.0; }
    
    let x2 = x * x;

    // 下から上へ、2k+1 の値を 27, 25, ..., 3 と直接計算
    // ループによるカウンタの更新や分岐を排除
    let mut res = x2 / 27.0;
    res = x2 / (25.0 - res);
    res = x2 / (23.0 - res);
    res = x2 / (21.0 - res);
    res = x2 / (19.0 - res);
    res = x2 / (17.0 - res);
    res = x2 / (15.0 - res);
    res = x2 / (13.0 - res);
    res = x2 / (11.0 - res);
    res = x2 / (9.0 - res);
    res = x2 / (7.0 - res);
    res = x2 / (5.0 - res);
    res = x2 / (3.0 - res);

    x / (1.0 - res)
}

結果：3ms程度早くなった

for文をunroll(展開)した結果( 3ms早くなった）

さいごに

今回は、ランベルトの連分数展開を用いてtanを実装してみました。連分数でtanを近似することで、極でも安定してtanを求められるのが実感でき、とても面白い体験ができました。

for文のunrollという単純な高速化をやってみましたが、次回はもう少しちゃんとlogなどの実装で触れたホーナー法 + FMAでマシンフレンドリーな実装にすることで高速化し、f64::tanにどの程度肉薄できるか挑戦してみようと思います。

次回：関連記事は、2026年6月1日に公開予定 (あと1日)