【Rustで作る電卓】sin関数実装 Rustで挑む三角関数の実装(1) マクローリン展開と範囲還元 #18

今回から三角関数の実装に進みたいと思います。いろいろ触ってみていると一筋縄ではいかないことが多く大変そうです。今回はとりあえず超シンプルにマクローリン展開で実装し、問題点を随時修正していくことにします。

1. sin関数(正弦波関数)の計算

sin関数の近似を行う手法として、一番初めに思いつくものとして、マクローリン展開(テイラー展開でx=0)にしたものがある。sin関数をマクローリン展開すると以下のような式になります。
\[
\sin(x)
= x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!} + \cdots
= \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{(2n+1)!}
\]
理想的に桁落ちなどがない完全な計算機なら無限回までやればいいですが、無限回計算するは現実的ではないのでf64精度を出せる16~17項目まで計算させてみます。

\[
\sin(x)
\simeq \sum_{n=0}^{17} (-1)^n \frac{x^{2n+1}}{(2n+1)!}
\]

# sin 関数を超シンプルに実装してみた例
pub fn sin(x: f64) -> f64 {
    let mut sum = 0.0;
    let mut term = x; // 第1項は x^1 / 1!
    let x2 = x * x;

    // 11項でxの21乗までになるのでちょっと多めに17項目まで加算する
    for i in 1..=17 {
        sum += term;
        
        // 次の項を計算するための漸化式:
        // 次の項 = 前の項 * (-x^2) / ((2n) * (2n+1))
        let n = i as f64;
        term *= -x2 / ((2.0 * n) * (2.0 * n + 1.0));
    }

    sum
}

1.1 実験：シンプルなマクローリン展開での実装

シンプルなマクローリン展開での実装と、rustで標準搭載されているf64:sin()を動かしてみて結果(result)と計算にかかった時間(1万回　duration)を比較してみました。

sin(0.5 * π)の結果

sin(0.5 [rad])の結果

sin(3 [rad])の結果

sin(50 [rad])の結果

結果まとめ
計算結果を見ると、1.5[rad]程度の入力であれば、シンプルなマクローリン展開での実装でも標準関数と同じ精度(f64の限界精度)を出せることが確認できた。しかし、値が3を入れた時点で、小数15位以降の計算結果に誤りが生じました。50など大きな数を入れると実用的に問題になるどころか、とてつもなく大きな数が出るような計算結果となり実用性は0です。

1.2 精度が破壊されれる原因

マクローリン展開の実装で、大きな入力値（例えば \(x=50\)）を入れたときに結果が暴走した原因は、主に2つの数値計算上の問題にあります。

① 巨大な中間項による「桁落ち (Cancellation)」

マクローリン展開の各項 \(\frac{x^{2n+1}}{(2n+1)!}\) を見てみましょう。
\(x=50\) のとき、分子の \(x^{2n+1}\) はとてつもない大きさになります。例えば \(n=25\) あたりでは、分子は \(50^{51} \approx 4.5 \times 10^{86}\) という、f64 が扱える最大値（約 \(1.8 \times 10^{308}\)）には収まるものの、非常に巨大な数値です。

この巨大な数値を計算したあと、次の項を引く……という動作を繰り返すと、「非常に大きな数同士の引き算」が発生します。浮動小数点数は有効数字（約15〜17桁）が決まっているため、大きな数同士を引くと、本来残したかったはずの微小な桁の情報が消し飛んでしまいます。これが「桁落ち」です。

例：X=4 n=10で 80桁精度とf64精度で計算したときの差

high(計算精度80桁) : 8.6082705155892152154811431214216716548470129261336799149411740473722988192930382E-8
f64 : 1.9572941063391263e-20
abs err : 8.6082705155872579213748039951621518591087832995654166697297820800627704301968089E-8

結果は歴然で、x=4 (πよりちょっと大きめ）程度でも大きな差が生じる。

② 収束の遅さ

マクローリン展開は \(x=0\) 付近での近似精度は最強ですが、\(x\) が 0 から離れるほど、正確な値を出すためにより多くの項数が必要になります。\(x=50\) ほどの距離になると、たとえ17項計算しても誤差を埋めきることはできず、級数が収束する前に数値が破綻してしまいます。

1.3 解決策：範囲還元 (Range Reduction)

「大きな \(x\) で精度が出ないなら、マクローリン展開が得意な \(x \approx 0\) の範囲まで \(x\) を小さくしてしまえばいい」というのが一般的な解決手法です。
sin(正弦波関数)は「\(2\pi\) ごとに同じ値を繰り返す」という周期性があるためこれを利用します。

どんなに大きな \(x\) が来ても、まずは \(x \pmod{2\pi}\) を計算して \([0, 2\pi)\) の範囲に押し込めます。
さらに、\(\sin\) の対称性を利用すると、計算範囲をさらに小さく限定できます。

• \(\sin(x + \pi) = -\sin(x)\)
• \(\sin(x) = \sin(\pi – x)\)

これらを組み合わせることで、最終的に計算機が解くべき範囲を \([-\pi/4, \pi/4]\) （約 \(-0.78\) 〜 \(0.78\)）まで絞り込むことができます。この範囲であれば、マクローリン展開は数項で f64 の限界精度を叩き出せます。

コラム:なぜ \(\pi/4\) を目指すのか

ここまでの解説で、どんな巨大な数でも \(\pm \pi/2\) （約 \(\pm 1.57\)）の範囲に収められることが分かりました。しかし、近代的な実装ではここからさらにもう一歩踏み込み、\([-\pi/4, \pi/4]\) （約 \(\pm 0.78\)）を計算の基本単位とします。
なぜこれほどまでに「小さくすること」に執着するのでしょうか？そのステップを順を追って見ていきましょう。

ステップ1：周期性による「一周」への押し込み

まず、\(\sin(x) = \sin(x + 2n\pi)\) を利用します。

• 操作: \(x \pmod{2\pi}\)
• 結果: どんな巨大な数も \(0 \le x < 2\pi\) の範囲に収まります。

ステップ2：対称性による「半円」から「左右」への追い込み

次に、\(\sin(x)\) の山と谷が鏡合わせになっている性質（\(\sin(x) = \sin(\pi – x)\) など）を利用します。

• 操作: グラフの後半分（\(\pi\) 以降）を負の範囲へ、右半分（\(\pi/2\) 以降）を折り返します。
• 結果: \([-\pi/2, \pi/2]\) の範囲に収まります。

ステップ3：余角の公式による「\(\pi/4\)」への最終還元

ここからどうやって\(\pi/4\)に押し込むかというと、余角の公式（\(\sin(x) = \cos(\pi/2 – x)\)）という関係を使います。
もし、入力された \(x\) が \(\pi/4\)（約0.78）を超えて \(\pi/2\) に近い場合、そのまま計算する代わりに、「\(\pi/2\) からどれだけ離れているか」を計算の対象にするのです。(単純に言えば45度を越えたらy軸から見た時の角度の方が小さいことを生かす方法です)

1. もし \(x \in [0, \pi/4]\) なら: そのまま \(\sin(x)\) を計算する。
2. もし \(x \in [\pi/4, \pi/2]\) なら: \(t = \pi/2 – x\) と置く。すると \(t\) は必ず \([0, \pi/4]\) に収まる。
   • このとき、\(\sin(x) = \sin(\pi/2 – t) = \mathbf{\cos(t)}\) を計算すれば良い。

なぜここまでやるのか？（精度とコストの損益分岐点）

マクローリン展開（多項式近似）において、「近似する範囲を半分にする」ことは、計算コストを劇的に下げ、精度を安定させることに直結します。

計算範囲	必要項数（f64精度）	最大値のべき乗 (xn)
\([0, \pi/2]\)	約11〜13項	\(1.57^{25} \approx 8.2 \times 10^{4}\)
\([0, \pi/4]\)	約7〜8項	\(0.78^{15} \approx 0.025\)

• 項数が減る: ループ回数や乗算回数が減り、純粋に速くなります。
• 桁落ちを防ぐ: \(x\) が \(1\) を超えない範囲（\(0.78..\) 以下）であれば、べき乗を重ねても値が巨大化しません。つまり、前述した「巨大な数同士の引き算」が発生せず、精度が極めて安定します。

結論 \(\pi/4\)を目指す理由

「\(\sin\) だけを頑張って広い範囲で解く」のではなく、「\(\sin\) と \(\cos\) をセットで用意し、どちらも一番得意な \([0, \pi/4]\) という極小範囲だけで戦わせる」。これが、一般的な計算機で\(\pi/4\)にまで範囲還元する理由です。

1.4 実装：範囲還元付き sin 関数

周期性と対称性を利用した、改良版の sin 関数を実装してみました。(とりあえず今回は余角の関係を使わずsinだけで実装しています)

use std::f64::consts::PI;

/// 改良版：範囲還元（Range Reduction）を取り入れた sin 関数
pub fn sin(mut x: f64) -> f64 {
    let tau = 2.0 * PI;

    // 1. 周期性を利用して [0, 2π) に収める
    x %= tau;
    if x < 0.0 { x += tau; }

    // 2. 対称性を利用して [-π/2, π/2] に収める(シンプル版余角を使わない)
    let mut sign = 1.0;
    if x > 1.5 * PI {
        x -= tau;        // [1.5π, 2π) -> [-0.5π, 0)
    } else if x > 0.5 * PI {
        x = PI - x;      // [0.5π, 1.5π] -> [-0.5π, 0.5π]
    }

    // 3. 小さくなった x でマクローリン展開
    let mut sum = 0.0;
    let mut term = x;
    let x2 = x * x;

    for i in 1..=10 {
        sum += term;
        let n = i as f64;
        term *= -x2 / ((2.0 * n) * (2.0 * n + 1.0));
    }

    sum
}

1.5 範囲還元後実験

\(\pi/2\)までの範囲還元を行った実装とstdとの比較結果を示します。

sin(0.5 * π)の結果

sin(0.5 [rad])の結果

sin(3 [rad])の結果

sin(50 [rad])の結果

結果まとめ
計算結果を見ると、3[rad]を入れると小数第16位以降の計算結果に誤りが生じました。これは、項数を12に減らした結果ではなく、f64を途中の計算に用いたためによる丸め誤差である。また、範囲還元を行ったことで前回は見当はずれな結果となった50を入れても誤差を小数第15位程度に抑えることができている。

1.6 考察と次回の課題

範囲還元を導入したことで、入力値が \(50\) という大きな数であっても、計算結果が「爆発」することなく \(\sin\) 関数らしい値を返せるようになりました。しかし、厳密に f64::sin() と比較すると、まだ無視できない微小な誤差が残っています。

残された課題1：浮動小数点数による「\(\pi\)」の限界

今回の実装では、範囲還元に x % (2.0 * PI) を使用しました。しかし、この \(2\pi\)（tau）自体が f64 という有限の精度で表現された近似値であるという点が問題になります。

入力 \(x\) が大きくなればなるほど、この「わずかな \(\pi\) の精度のズレ」が蓄積し、剰余（あまり）の結果が本来の値から大きく乖離していきます。数千、数万という大きな入力を扱うには、単なる mod 演算ではない、より高度な引数還元アルゴリズム（Payne-Hanek Reductionなど）の検討が必要になります。

残された課題2：計算コストと「\(\pi/4\)」への還元

前述のコラムで触れた通り、現在は \(\pm \pi/2\) までしか還元できていません。さらに精度を高め、項数を減らして高速化するためには、\(\pm \pi/4\) への還元が不可欠です。

そのためには、\(\sin\) だけでなく \(\cos\) のマクローリン展開も実装し、状況に応じて「\(\sin\) として計算するか、\(\cos\) として計算するか」を切り替えるロジックを組む必要があります。

まとめ

今回は、Rustを使って \(\sin\) 関数の基礎的な実装に挑戦しました。

• マクローリン展開の限界: \(x=0\) から離れると、桁落ちと収束速度の問題で精度が崩壊する。
• 範囲還元の威力: 周期性を利用して入力を小さくすることで、巨大な入力でも計算が可能になる。
• 実用化への壁: 単純な f64 の \(\pi\) を使った還元では、大きな数で精度が低下し始める。

当初は \(\sin\) を単体で完成させる予定でしたが、精度と効率を追求すると、\(\cos\) との同時実装が最適な実装パスであることが分かってきました。

次回は、\(\cos\) の実装を加え、「余角の公式」を用いた \(\pm \pi/4\) への完全な範囲還元、そして大きな入力値でも精度を落とさないための工夫について掘り下げていきたいと思います。

次回 cosの実装と範囲還元改良 :関連記事は、2026年5月27日に公開予定 (あと16時間)