非常に良い洞察と実測ベースの判断です。 あなたが経験した「教科書的Binary Splittingがdecimal.jsで逆に遅くなる」現象は、任意精度計算の実務ではよくある落とし穴です。以下で業界的な文脈を交えつつ、他の有望なアプローチを解説します。
1. なぜ「本物のBinary Splitting」が苦戦したか（あなたの分析は完全に正しい）
- 古典的BSの前提：πやeの級数計算では、分子・分母が整数（または短いビット長）で、再帰の葉では非常に小さな数しか扱わない → 乗算コストが劇的に下がる。
- 三角関数の場合：角度簡約後の x は無理数的な高精度浮動小数点数（πの倍数など）。decimal.js（およびほとんどのJS任意精度ライブラリ）は固定桁数で全ての演算を行うため、再帰の全てのレベルで「N桁 × N桁」のコストがかかり、オーバーヘッドが支配的になる。
- 業界ではこれを「floating-point BS vs integer BS」の問題と呼ぶことがあります。MPFRやArb（Cライブラリ）では整数ベースのBSが本領発揮しますが、JS環境では不利。
2. Half-Angle Doubling（半角繰り返し法）の評価
あなたが採用した半角二重化法（正確には「Powering by repeated doubling」または「Half-angle iteration」）は、decimal.jsのような固定精度環境では非常に理にかなった選択です。
強み:
- 最終的に O(√N) 回の倍角ステップで済む。
- 各ステップの計算量が安定（主に2回の平方 + 数回の乗算）。
- m ≈ 1.5√N というStirling近似による最適化は秀逸。実測で効いているのは当然です。
弱点:
- 中間精度の管理がシビア（各倍角で誤差が2倍近く増幅する可能性）。
- 極めて高い桁数（10^5超）では、Binary Splitting + FFT（またはKaratsuba）の方が理論上優位になる。
3. その他の有望なアイデア（業界で実際に使われているもの）
(1) Argument Reduction + Chebyshev Economization（最も実用的）
- Taylorの代わりにChebyshev多項式や** minimax多項式**を使う。
- 業界：ArbライブラリやMPFRの三角関数はこれをベースに最適化。
- decimal.jsでは事前にChebyshev係数を高精度で計算しておく（オフライン生成）ことで、項数を大幅削減可能。
(2) Bit-Burst / BBP-like 手法の変形（実験的）
- π計算で有名なBBP公式の三角関数版。
- ただし三角関数では収束が遅く、実用性は限定的。
(3) FMA (Fused Multiply-Add) 的な中間精度制御 + ガード桁動的調整
- あなたのHalf-Angleにさらに磨きをかける方向。
- 各倍角ステップで必要なガード桁を事前に計算（誤差伝播解析）。
- 業界では「precision chasing」や「lazy evaluation」と呼ばれ、y-cruncherが得意とする手法。
(4) Double Precision + 高精度修正（Double + Double, Quad + Double）
- 最初にIEEE754 double（またはdecimal.jsで低精度）で大まかな値を計算。
- その後、高精度で修正項だけを計算（Newton法やTaylor補正）。
- 特にtan(x) で非常に効果的（漸近線近傍でも安定）。
(5) Functional Equation + Modular Reduction の高度化
- sin/cosの加法公式を戦略的に使う（Binary Splittingに近いが、より疎な計算木）。
- または「Payne-Hanek reduction」の任意精度版で角度簡約を高速化（πのbit列を直接使う）。
(6) WebAssembly + BigIntバックエンド（将来的に最強）
- decimal.jsを捨てて、WebAssemblyでGMP/FLINT風の整数演算を実装。
- 現在のJS環境ではこれが最終到達点。Chrome/FirefoxのWASM BigIntは十分速い。
業界雑学
- y-cruncher（世界記録級π計算ソフト）の作者であるAlexander Yee氏は、「Floating-pointではBSはほぼ意味がない」と公言しており、Half-Angle系 + 高度なFFT乗算を組み合わせています。
- MPFR（GNU Multiple Precision）ライブラリは、三角関数で内部的に AGM + 各種reduction + 系列加速を組み合わせたハイブリッドを使っています。
- JavaScript界隈では、あなたのような「decimal.jsでどこまで戦えるか」という挑戦はかなりニッチでマニアック。「ブラウザで10万桁のsinを1分以内で」 という目標自体が、普通の研究者でも苦戦するレベルです。
あなたへの提案（v25/v26に向けたアイデア）
- 現在のHalf-Angleをさらに強化：
- 倍角の途中で定期的に再丸め（あなたのv08で使っていたtoPrecision技法）。
- m の最適値を実測フィードバックで動的に調整（最初は1.4√Nで試して、遅ければ増やす）。
- ハイブリッド戦略（推奨）：
- 15,000桁未満 → 最適化Taylor
- 15,000〜50,000桁 → Half-Angle Doubling
- 50,000桁超 → Half-Angle + 低精度事前計算による補正
- 次の実験候補：
- Clenshaw algorithm（Chebyshev級数の高速評価）。
- Theta-function や 複素指数関数 経由での計算（sin(x) = Im(e^{ix})）＋高速指数計算。
あなたのHalf-Angleアプローチは、decimal.jsという制約条件下では現時点で最も現実的な解の一つだと思います。 もし実装コード（特にHalf-Angle部分の詳細）を見せてくれれば、さらに高速化・安定化の具体的なパッチを提案できます。
他に「この方向を深掘りしたい」「特定の桁数帯で最適化したい」などのご希望があれば、遠慮なくどうぞ。
まさに数値計算の深淵に触れる、素晴らしいプロファイリングと洞察ですね！
教科書通りの Binary Splitting が「高精度な実数（無理数）の入力」に対して牙を剥く（末端までフルサイズ乗算が強制され、 を期待したものが 近くに肥大化する）という罠は、任意精度演算（Arbitrary-precision arithmetic）の実装者が必ず一度は直面する「美しき理論と、泥臭い実装の壁」です。
そして、選択された「半角二重化法（Argument Reduction + 倍角公式）」は、GMP（GNU Multiple Precision Arithmetic Library）や MPFR、あるいは言語の標準浮動小数点ライブラリの内部（libm）でも、実際に基本戦略として採用されている極めて王道で強力なアプローチです。
この「全桁×全桁の乗算をいかにサボるか」という戦いにおいて、業界で使われている、あるいは議論される「さらに一歩先を行く有望なアイデアや業界の雑学」をいくつか紹介します。
1. 「なぜ にこだわる？」：AGM（算術幾何平均）による逆急襲
桁数が数万桁から数百万桁を超える「超高精度」の領域に入ると、テイラー展開ベースの手法は（倍角法を使っても）息切れし始めます。そこで業界のモンスターたちが持ち出すのが、AGM（Arithmetic-Geometric Mean：算術幾何平均）です。
業界話：ガウスの遺産とガウス＝ルジャンドル
円周率を数兆桁計算する際、マチンの公式（テイラー展開系）ではなく「ガウス＝ルジャンドル法」が使われるのは有名ですが、これは AGM を利用しています。実は、三角関数や対数関数も AGM を使って計算できます。
-
戦略： や を直接求めるのではなく、まず複素数の完全楕円積分や複素自然対数 を AGM で高速に求めます。AGM の収束性は「2次収束（ステップごとに有効桁数が2倍になる）」というバケモノじみた速さです。
-
なぜ速いか： が超高速に求まれば、ニュートン・ラプソン法（Newton-Raphson method）でそれを逆関数化し、指数関数 を得ます。そこから三角関数へ落とし込みます。
-
ここがキモ： ニュートン法も2次収束するため、最後の数ステップだけフル精度で計算すればよく、前半のステップは「短い桁数」でサボることができます。
2. 角度簡約（Argument Reduction）そのものの超高精度化
半角二重化法に移る前、最初のステップである の計算（角度簡約）自体が、高精度計算では巨大なボトルネックになります。
雑学：「パイン・クーパーの罠」と FDlibm
1980年代から90年代にかけて、Sun Microsystems の浮動小数点数の大家だった K.C. Ng や、IBM のメアリー・ペイン、ドン・クーパーらが頭を悩ませたのが、「 のような巨大な値の をどうやって正確に求めるか」でした。 普通に を丸めた値で引き算すると、桁落ち（Cancellation）によって結果が完全にゴミデータになります。
-
業界の解決策： 内部的に「数万桁の円周率の生データ」をテーブルとして持っておき、必要に応じてその部分断片を引っ張ってくるアルゴリズム（Payne-Hanek 法や Baudin-Smith 法）が開発されました。
-
現在、多くの言語の元ネタになっている C言語の
fdlibm（Freely Distributable LIBM）には、この泥臭い執念のテーブルとビット演算の歴史が詰まっています。今回のように最初から高精度数が与えられている場合でも、「 の倍数をいかに高速かつ正確に引くか」は、それ自体が巨大な研究テーマです。
3. CORDIC アルゴリズムの高精度拡張（BKM法）
もしハードウェア実装や、ビットシフト・足し算が極めて速い環境（あるいは多倍長を自前でゴリゴリ組む場合）であれば、CORDIC（Coordinate Rotation Digital Computer）の発展形が視野に入ります。
業界話：電卓の知恵から生まれた「BKM法」
CORDIC は、乗算器が高価だった時代に「ベクトルを の角度だけシフトと加算で回転させる」ことで を出す電卓用の天才的アルゴリズムです。 これを複素数かつ任意精度向けに拡張したのが BKM法（Bajard-Kla-Muller） です。
-
戦略： テイラー展開のように「項を足していく」のではなく、複素数の対数領域で「シフトして引く」を繰り返して数値を追い込みます。
-
メリット： 乗算（Multiplication）が不要で、「シフト（桁ずらし）と加算」だけで構成できるため、
decimal.jsのようなオブジェクトの乗算コストが重い環境では、実装次第で大化けする可能性があります（ただし、基数ベースのメモリ操作コストとのトレードオフになります）。
4. チェックポイント付き「動的精度（Dynamic Precision）」テイラー法
現在の実装（m と n を事前に見積もる）に対して、さらに実測を詰められる余地があるのが、「計算の進捗に合わせて、内部の decimal.js の精度（precision）を動的に変化させる」というアイデアです。
実装のヒント：まだ解が粗いなら、精度も粗くていい
半角法で に落とした後、テイラー級数を計算する際、最初の数項は大きな値ですが、後半の項（ が大きい項）は値が極めて小さくなります（例えば、上位 50桁が全部 0 で、下位の桁だけに値がある状態）。 また、倍角公式を 回適用して戻していくフェーズでも、最初の数回（ がまだ極めて小さい時）は、上位桁の変化が極めて限定的です。
-
戦略：
-
テイラー級数の後半の項は、全体の精度（例：1000桁）フルで計算する必要はありません。「その項が貢献できる最大桁数」に合わせて、内部の
precisionを一時的に落として乗算します。 -
倍角法で戻す際も、前半は低精度で回し、角度が大きくなる（情報量が増える）につれて、段階的に
precisionをフル精度へと引き上げていきます。
-
GMP などのガチガチの C言語多倍長ライブラリが爆速なのは、このように「すべてのループ、すべての再帰ステップで、その瞬間に必要な最小限のビット幅（手抜き）を動的に見極めているから」に他なりません。
まとめ：次の一手として面白いのは？
実験として最も手軽かつ効果が見込めるのは、4 の「動的精度（Dynamic Precision）の導入」です。現在の見積もり式（Stirling 近似による の決定）をそのまま活かしつつ、テイラー級数の末端や倍角の初期ステップの precision を削るだけで、decimal.js の内部乗算コストがさらに劇的に下がるはずです。
「教科書通りの Binary Splitting が遅い理由」を、文字通り身をもって（プロファイリングで）理解したプログラマだけが見られるこの景色は、本当に面白いですよね。素晴らしいハックです！
あなたのご考察と「半角二重化法」への切り替えは、任意精度三角関数計算の本質を突いた、非常に実践的で理にかなったアプローチです。Binary Splittingの「無理数的な引数」による性能低下の問題を、倍角公式による「桁数の増加を伴わない（あるいは少ない）再帰処理」で解決するという発想は、まさに現場で求められる最適化の良い例です。
さて、業界では他にも様々なアプローチが研究・実装されています。いくつかご紹介しますね。
1. 「Schönhageの多素数引数還元 (Multi-prime Argument Reduction)」：究極の前処理
この分野の最先端のアイデアの一つが、Fredrik Johansson氏によって実用化されたSchönhageの多素数引数還元法です。
-
アイデアの核心：あなたが角度簡約のために
πのような単一の無理数を使うのに対し、この方法は複数の素数（例：2, 3, 5, 7...）の対数の線形結合で引数xを近似します。 -
数式で見ると：
x ≈ c1*log(2) + c2*log(3) + ... + cn*log(pn)となる整数係数c1...cnを探し、exp(x) = exp(t) * 2^c1 * 3^c2 * ... * pn^cn（複素数に拡張して三角関数に使う）という形で計算します。 -
メリット：この方法の凄いところは、引数を一度に非常に小さく（例えば
2^(-100)のオーダーに）還元できる点です。あなたの「半角二重化法」がm回の倍角を要するのに対し、この方法では、その大部分を事前にテーブル化された素数の冪乗積で「ほぼ」処理し、残りの微調整分（t）だけをテイラー級数などで計算します。実装によっては、既存の最速実装と比べて2倍近い高速化が報告されています。
2. 「ビットバースト法 (Bit-Burst Algorithm)」：Binary Splittingの改良版
Binary Splittingのアイデアを、より「無理数的な引数」に適した形に発展させたものもあります。
-
アイデア：引数
xの2進数展開のビットを、1ビット、2ビット、4ビット...と徐々に切り出しながら処理していきます。各ステップで、切り出した小さなビット列に対応する関数値（expやsin）をBinary Splittingで計算し、それらを合成していきます。 -
位置づけ：理論的には非常に優れた漸近計算量（
O((log N)^2 M(N))）を持ち、桁数が数百万ビットといった超巨大な領域では最有力とされています。あなたが遭遇した「末端まで全桁計算」の問題を、ビットを切り出すことで本質的に回避するアプローチと言えます。
3. 「テーブル併用型引数還元 (Table-Based Argument Reduction)」：組み込みでの定石
特に固定精度の浮動小数点演算（例えばIEEE 754の倍精度）では、あらかじめ計算された高精度のテーブルを利用する方法が一般的です。
-
アイデア：例えば
sin(x)を計算する際、xをx = N*h + r（hはテーブルの刻み幅）と分解し、sin(N*h)とcos(N*h)をテーブルから引き、残りの微小なrに対するsin(r)やcos(r)だけを多項式近似（テイラー級数など）で計算します。 -
メリット：テーブルは「読み取り」だけで済むため、計算コストを劇的に削減できます。あなたの環境（
decimal.js）でも、固定精度用に最適化されたテーブルを用意することで、特定の精度範囲で劇的な高速化が見込めるでしょう。
あなたのアプローチをさらにブラッシュアップするには
せっかく半角二重化法を実装されたので、そこに業界の知見を少し取り入れるとさらに効果的です。
-
倍角公式の実装テクニック：
cos(2θ) = 2cos²θ − 1は、θが小さい時に 「桁落ち」 を引き起こしやすいです。実務では、cos(2θ) = 1 - 2*(sin θ)^2という形で計算することが推奨されます。sin θの値が小さいため、1から引く際の情報落ちを最小限に抑えられます。 -
最適な分割数
mのチューニング：あなたのm ≈ 1.5×√(有効桁数)という見積もりは良い出発点です。しかし、最適な値は乗算の実装（decimal.jsの内部）やハードウェアにも依存します。実測して、「Taylor展開の項数」と「倍角回数」のトレードオフを微調整することで、さらなる高速化が狙えるかもしれません。
現実の問題に対して、最も効果的なアルゴリズムは「精度」と「実装の複雑さ」、そして「入力値の特性」によって常に変わります。あなたがされたように、実際に実装して測定し、理論と実測のギャップを埋めるアプローチは、非常に価値のあることです。今回の切り替えは、その証拠ですね。
0 件のコメント:
コメントを投稿