隠れた高速asin()の存在

75日前原文(16bpp.net)

概要

PSRayTracing の最適化過程を通じて 三角関数近似 の検証
Taylor展開 や Padé Approximant による高速化の試行
AIツール（LLM） から得たアルゴリズムが最良の結果
実際のレンダリング速度 の比較と効果測定
最終的な推奨手法 とその評価

PSRayTracingにおける三角関数近似の最適化記録

PSRayTracing プロジェクトの最適化過程の記録
目標： asin(x) の高速かつ十分な精度の近似手法の探索
Taylor級数 による近似導入、精度と速度のバランスを検証
xが -0.8未満または0.8超 の場合、Taylor近似の誤差が大きいため std::asin() にフォールバック
4次のTaylor展開が最もパフォーマンス良好（+5％高速化）
- C++実装例：
  - xの範囲で分岐し、範囲外は標準関数を利用
  - 係数を定数で保持し、計算コスト低減

Padé Approximantによるさらなる高速化の試み

Padé Approximant は関数近似の数学的手法
Pythonで [3/4] や [5/4] Padé近似を実装し、Taylor展開と比較
- 誤差が減少し、特に中心付近で精度が向上
Half Angle Transform （半角変換）を併用し、端の誤差を補正
- |x| > 0.85で半角変換を用い、中心付近でPadé近似を利用
[1/2] Padé を端の補正用に追加し、さらに計算を簡略化
- C++でも同様のロジックで実装

実際のレンダリング速度比較

ベンチマーク環境 ：M4 Mac Mini, macOS Tahoe, GCC15（-O3最適化）
std::asin() ：約111秒
Taylor級数近似 ：約105秒（約5％高速化）
Padé Approximant + 半角補正 ：約105秒（Taylorと同等）
- 精度は向上したが、速度面で大きな差は見られず

AI（LLM）による最適解の発見

Gemini（LLM） に「C++で使える高速なasin(x)近似は？」と質問
Nvidia Cg Toolkit 由来の Minimax多項式近似 を返答
- Horner法による多項式評価＋平方根＋π/2オフセット
- PythonおよびC++で実装が容易
実際のレンダリング速度 ：約101秒（最速）
精度：標準asin関数とほぼ同一、違いは視認できないレベル
他環境（Intel i7, Ubuntu, GCC/clang） でも約1.5倍高速化を確認

結論と推奨手法

AIツールから得たMinimax多項式近似 が最も高速かつ高精度
Taylor級数やPadé Approximant は理論的には有効だが、実装・運用面ではAI提案手法が優秀
今後の最適化指針 ：AIツールの活用と、既存の数学的近似手法の併用検討

実装例 （C++）：

constexpr double HalfPi = 3.1415926535897932385 / 2.0;
double fast_asin_cg(const double x) {
    constexpr double a0 = 1.5707288;
    constexpr double a1 = -0.2121144;
    constexpr double a2 = 0.0742610;
    constexpr double a3 = -0.0187293;
    const double abs_x = fabs(x);
    double p = a3 * abs_x + a2;
    p = p * abs_x + a1;
    p = p * abs_x + a0;
    const double x_diff = sqrt(1.0 - abs_x);
    const double result = HalfPi - (x_diff * p);
    return copysign(result, x);
}

教訓：最適化前に十分なリサーチを行い、既存の知見やツールを活用することの重要性

Hackerたちの意見

速いアプローチってSIMDじゃない？[編集: それともGPU？] 1.05倍から1.90倍のスピードアップは素晴らしいけど、16倍のスピードアップはもっといいよね！おそらく、これらは独立した改善点かもしれないけど、優先順位をつけるならまずSIMDを選ぶかな。Intelのインストリンシックガイドでasinを探してみたら、AVX-512のインストリンシック _mm512_asin_ps があったけど、「シーケンス」って書いてあって、単一命令じゃないみたい。実際に使われているシーケンスはどこかのヘッダーファイルにあると思うけど、どこを見ればいいかわからないから、fast_asin_cgのSIMD版と比べることもできないな。 https://www.intel.com/content/www/us/en/docs/intrinsics-guid...

└

レイトレーシングについてはあまり詳しくないけど、あのasin呼び出しをうまく調整するのは難しいと思う。入力と出力のメモリが整列して連続している必要があるからね。俺の無知な直感では、どのピクセルがどの分岐を取るかはあまり規則性がないように思うし、どのasin呼び出しが必要になるかもわからないけど、間違ってるかもしれない。

└

浮動小数点の作業はあんまりやらないけど、普通のサイン命令はないと思う。古いx87浮動小数点スタックのやつしかないし。どんな「シーケンス」になるのか気になったけど、僕のコンパイラはその内蔵関数には古すぎるんだよね。Godboltでもgccやclangには役立たなかったけど、iccがコールを生成することはわかったよ。https://godbolt.org/z/a3EsKK4aY

たぶん速くはならないと思うけど、書くのは楽しかったよ： float asin(float x) { float x2 = 1.0f-fabs(x); u32 i = bitcast(x2); i = 0x5f3759df - (i>>1); float inv = bitcast(i); return copysign(pi/2-pi/2*(x2*inv),x); } 悪名高い浮動小数点ビットハックのおかげ。

└

浮動小数点ビットハック禁断の魔法

└

ザルゴを持ち込んだね。今の10年は君のせいだ。

└

float asinine(float x) { FTFY :P

└

// なんだこれ

└

字幕があったらいいかも。

└

この方法の悪いところは、ネイティブのCPU命令より遅いこと。良いところは、少なくともxの値が1.0と-1.0の時に結果が非常に正確だってことだね。JK

一般的に、ミニマックス近似はあまり評価されていないツールだと思う。特に、最適な多項式近似を生成するためのシンプルなレメズアルゴリズムはね。[0] いくつかの修正を加えれば、区間内の絶対誤差または相対誤差を最適化するように適応できるし、さらには有理関数近似も考案できるよ。（ただ、残念ながら、多くのアルゴリズムのプレゼンテーションでは、サンプルポイント選択があまりにも単純すぎて、特に小さな振動を含む非自明な入力曲線ではうまくいかないことがある。）[0] https://en.wikipedia.org/wiki/Remez_algorithm

└

Remezを「シンプル」って呼ぶのはどうかな…相対的なもんだし、僕はRemezよりもChebyshev近似の方がシンプルだと思うよ。

└

数学の授業ではテイラー/マクローリン級数をたくさん教えるけど（三角関数は「CORDIC」って呼ばれることもあるし、これも古い方法だよね）、実際のFPUやライブラリではあんまり使われてないんだよね。カリキュラムを更新して、もっと良い方法を知ってもらうべきかも。

OPが最後に良いミニマックス解を得られたのは嬉しいけど、記事が重要なポイントの一つを明確にしていないように思う。指定された区間内の誤差波形は重要で、特徴的なミニマックスのようなウィグルが見えないなら、改善のチャンスを無駄にしているってこと。テイラー級数は一般的に良い選択ではなく、テイラー級数のパデ近似も同様に良くない。パデ近似を使うなら、元の関数のものにすべきだと思う。俺はチェビシェフ近似が好きだな：https://www.embeddedrelated.com/showarticle/152.php これはしばしば複雑なレメズアルゴリズムにかなり近いから。

└

チェビシェフ多項式 cos(n arcos(x)) は、すべての連続関数 f:[0,1]->R が多項式関数で一様に近似できることの証明の一つを提供するんだ。バーンスタイン多項式は短い証明を提供するけど、最良の数値法とは言えないかもね: https://en.wikipedia.org/wiki/Bernstein_polynomial#See_also

この一文: > この素晴らしいコードのスニペットは、死んだソフトウェアのドキュメントに埋もれていて、その元の数式は60年代の数学の教科書に書かれていた、っていうのがなんか印象的だった。僕はこの手の仕事に少し関わったことがあって（既存のコードを改善する方法はほとんどないって結論に至った）、新しい特定のハードウェアやドメイン制約がない限り、あるいは何かのためにちょっと手軽にやりたい場合、もしくは研究論文レベルの時間と労力を投資する覚悟がない限りね。AbramowitzとStegunを「60年代の数学の教科書」って呼ぶのはちょっと面白いよね。Knuthの「コンピュータプログラミングの芸術」とかと同じくらいその分野では重要なものだし、マイナーなテキストじゃないよ。あんまり使わないと内容を忘れちゃうかもしれないけど、「ああ、もちろんそこに載ってるよね…」ってなるはず。

└

Abramowitz/Stegunは2010年に更新されて、今はここにあるよ: https://dlmf.nist.gov/

└

大学の授業で読まなかった本たちだよ。著名なテキストがいくつかスキップされてたのにはちょっと驚いてる。

└

この場合、AIは人間よりも歴史をよく理解してたね。

└

そうだね、ちょっとマイナーなものが欲しいなら、CodyとWaiteの「Software Manual for the Elementary Functions」を引っ張り出してみて。そしたら、ASINの実装がミニマックスだって分かるよ。

└

古典的な方法を改善する一つの方法は、分析的理想係数を最も近い浮動小数点数に近似するのではなく、その理想係数を出発点として、もう少し良いものを探すことだよ。SLEEFベクトル化数学ライブラリはこれを実現していて、通常、理論が予測するよりも低い多項式の次数で、全ての浮動小数点範囲に対して精度保証を提供できるんだ。asinf関数は、全ての単精度浮動小数点数に対して1 ULPの精度を持っていて、記事のasin_cgに似てるけど、主な違いは多項式の出力ではなく入力に対して平方根を計算することだね。

正確に言うと、これは元々1955年のHastingsからで、プリンストンの「デジタルコンピュータのための近似法 by Cecil Hastings」の159-163ページに載ってるんだ。実際には異なる定数を使った近似のバージョンがいくつかある。元の研究は1950年代のコンピュータに対してパフォーマンスを重視して行われたんだ。その後、有名なAbramowitzとStegunが許可を得てそれを4.4.45の式に入れたんだよね。それからnvidiaのCGライブラリがその式を基にしたコードを書いたんだ、たぶんいくつかの最適化もされてるだろうけど。

└

これを追いかけたのは、ベクトル化可能な関数近似に特に興味があったから。特に、範囲正規化を扱うためにビットバンギングを利用するやつね。（誰か良い参考文献持ってる？）残念ながら、これはHastings 1955のものじゃないんだ。Hastingsはテイラー級数とチェビシェフ多項式の近似を提供してるけど、OPの解法はパデ近似で、Hastingsでは全くカバーされてない。

これまでの作業とその話を考慮して、LLMに聞いてみることにしたんだ。7年前のStack Overflowの回答から自力で答えを出したのはすごいね！[1] この記事が公開される前は、「fast asin」の検索結果でこれが一番上に来てたはずだよ。[1]: https://stackoverflow.com/a/26030435

└

それは見たけど、そのページの大半はacos()について話してるんじゃないの？

テイラー多項式をホーナー法で評価する方がずっと良くない？（C++が自動的にこれをやってくれるかもしれないけど、丸め誤差があるかもしれないから、たぶん無理だろうね。）

どんなグラフィックスアプリケーションでも三角関数は頻繁に使われる。本人からの反論、「三角法を避ける」: https://iquilezles.org/articles/noacos/

DSP数学では、チェビシェフ多項式近似を使うのが一般的だよ。必要な範囲内で信じられないほど正確な結果が得られるんだ。

ハクソク